SHALLOW CRUSTAL STRUCTURE BENEATH THE JUAN DE FUCA RIDGE F R O M 2 SEISMIC R E F R A C T I O N D TOMOGRAPHY By Donald John White B. Sc. (Physics) University of Toronto M . Sc. (Geophysics) University of British Columbia A THESIS S U B M I T T E D IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS DOCTOR FOR T H E DEGREE OF OF PHILOSOPHY in T H E FACULTY OF GRADUATE GEOPHYSICS STUDIES AND ASTRONOMY We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA March 1989 © Donald John White, 1989 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. The University of Britisi Vancouver, Canada Date DE-6 (2/88) Abstract The formation of oceanic lithosphere along ocean ridges, and the role that crustal magma chambers play in the accretionary process, continue to be fundamental issues in plate tectonics. To address these issues, a multi-receiver airgun/ocean bottom seismograph refraction line, designed to allow definition of lateral velocity and attenuation variations within the shallow crust, was shot across the Endeavour segment of the Juan de Fuca Ridge near 48° N , 129° W. A tomographic inversion procedure has been developed to invert the first arrival travel times and amplitudes from this profile for 2 attenuation structure. D velocity and The inversion method is suited to multi-source, multi-receiver refraction profiles where source/receiver spacings are denser than for conventional profiles. The travel time-velocity inversion scheme is based on an iterative solution of the linearized problem and allows for determination of continuous velocity variations as well as geometry of subhorizontal interfaces. The iterative procedure requires a good initial estimate of the velocity model. In each iteration, two-point ray tracing is performed to construct a linear system relating travel time residuals to velocity perturbations. A damped least-squares algorithm is used to solve this system for a velocity perturbation which is used to update the current velocity estimate. Once the final velocity structure of the model has been determined, amplitudes can be inverted directly for attenuation. Tests to ascertain resolution of the method reveal horizontal smearing of the solution due to ray geometry, drop-off in resolution with depth, and the effects of source-receiver geometry and velocity structure on resolution. Parameter weighting is important in removing streaking effects (caused by inhomogeneous ray coverage) from the solution. For the purposes of ray tracing, the model is parameterized in terms of constant gradient (velocity and attenuation) cells, which allow use of analytic expressions for kinematic and dynamic n ray properties, attenuation and inversion quantities. This parameterization causes scat- ter in the amplitudes calculated using zero-order asymptotic ray theory, a problem which is remedied by smoothing the velocity models before amplitude calculation. Application of this 2 D tomographic inversion scheme to first arrival travel times and amplitudes for the cross-ridge refraction line produced a 4-layer model for the shallow _1 crust. Layer 1 is 250 — 650 m thick, with v = 2.5 km/s and = 0.5 s . t ~800 m thick, v = 4.8 km/s and V v 2 z _1 2 — 1.0 s . the sequence of extrusives whereas layer 3 (~800 and layer 4 (i> =6.3 km/s, V i ; = 0 . 3 4 z 4 massive gabbro sequence, respectively. _ 1 s ) Layer 2 is Layer 1 and layer 2 likely represent m thick, i> =5.8 k m / s , V i> =0.5 3 z 3 - 1 s ) are associated with the dike complex and A n abrupt velocity transition between layer 1 and layer 2 may be a metamorphic front within the pillow basalts. A low velocity- high attenuation anomaly (velocities decreased by < 0.4 km/s and Q ~20-100), which is interpreted as a zone of increased fracture porosity and/or permeability associated with axial hydrothermal circulation, exists beneath the ridge in layer 2 and upper layer 3. Smaller low velocity-attenuative zones in layer 2, located 8 km to either side of the ridge may be loci of off-axis hydrothermal circulation. No evidence is found for the existence of a crustal magma chamber in the depth range of 1.5 — 3.0 km below the seafloor. Tests indicate that a 1 X 1 km zone of partial melt represents the minimum dimension of such a feature that would be clearly detected by this refraction experiment. These results suggest that Endeavour Ridge may be experiencing a period of diminished magma supply with the magma chamber reduced or eliminated by hydrothermal circulation. Asymmetry of the velocity anomalies observed in layer 3 and layer 4 suggest that crustal temperatures are elevated by 125 — 2 0 0 ° C beneath the ridge and to the east relative to temperatures west of the ridge, indicating that a deep crustal or upper mantle melting anomaly may exist east of the ridge. Table of Contents Abstract List of Tables List of Figures Acknowledgement 1 2 ii viii ix xii Introduction 1 1.1 Endeavour Ridge Segment, Juan de Fuca Ridge 2 1.2 The SEISRIDG-85 experiment 7 1.3 Magma chamber detection 8 1.4 Thesis outline D 2 Travel Time & Amplitude Tomography 2.1 Introduction 2.2 Problem Formulation 2.3 2.4 12 13 13 , 14 2.2.1 Velocity Formulation 14 2.2.2 Interface Formulation . 18 2.2.3 Attenuation Formulation 20 2.2.4 Model Parameterization 22 2.2.5 Basis functions for Method A 23 Linear Inversion . . . 25 2.3.1 Regularization 27 2.3.2 Global Norms 2.3.3 Solution Method 30 : 31 Iterative Inversion Scheme 32 2.4.1 Grid Cell Size 34 2.4.2 Damping Factor 34 iv 2.5 Resolution 35 2.6 Inversion Examples 47 2.7 3 2.6.1 Velocity Inversion 47 2.6.2 Effect of Damping 50 2.6.3 Streak Removal 50 2.6.4 Interface Inversion 55 2.6.5 Combined Velocity-Attenuation Inversion Summary 57 63 D y n a m i c R a y Tracing 65 3.1 Introduction 65 3.2 67 Asymptotic Ray Theory 3.2.1 Restrictions of ART 68 3.3 Ray Kinematics 69 3.3.1 Analytic Ray Tracing Formulae 3.3.2 Cell Boundary Intersections 3.3.3 Interface Intersections and Snell's Law 3.3.4 Two-point Ray Tracing 70 73 74 75 3.4 Ray Dynamics 75 3.4.1 Connection Formula for Geometrical Rays 3.4.2 Amplitudes for Reversed Rays 3.4.3 Head Waves 75 84 85 3.5 Attenuation 86 3.6 Inversion Formulae 88 3.6.1 Travel time-Velocity Jacobian Formulae 88 3.6.2 Amplitude-Attenuation Inversion Formulae 90 3.7 Summary 90 v 4 Tomography Applied to the Cross-Ridge Line 4.1 The Cross-Ridge Line 4.1.1 Projection onto a Line 4.1.2 The Data 91 91 92 94 4.2 Travel Time Inversion 102 4.2.1 Starting Model 103 4.2.2 Inversion for Shallow Velocity Structure 104 4.2.3 Inversion for Deep Velocity Structure 108 4.2.4 Travel Time Inversion Results 4.3 Amplitude Inversion 112 115 4.3.1 Amplitude Estimation for Inversion 4.3.2 Scaling of the Calculated Amplitudes 4.3.3 Amplitude Calculation for the Velocity Model 4.3.4 Application to the Data 123 4.3.5 Amplitude Inversion Results 133 5 Discussion of Inversion Results 117 118 120 135 5.1 Layered Nature of the Velocity Model 135 5.2 146 Magma Chamber, Velocity and Attenuation Anomalies 6 Summary 8z Conclusions 150 References 156 Appendices 166 A Data Reduction 166 A.l OBS Timing Corrections 166 A.2 OBS Repositioning &; Final Timing Correction A.2.1 OBS Repositioning Algorithm vi 167 170 A.3 Data Playback and Digitization 171 A.4 Amplitude Compensation 171 A.5 OBS-Playback Gain Factors 174 vii List of Tables 4.1 Seafloor velocities 103 A.l OBS timing corrections: drift and skew 167 A.2 OBS repositioning and static time shifts 168 A.3 Amplitude compensation parameters 174 A.4 OBS/playback system gain factors 176 A.5 Normalization factors from water wave matching 176 vin L i s t of Figures 1.1 Bathymetric and profile location map for Juan de Fuca Ridge 3 1.2 Endeavour Ridge bathymetry 4 1.3 Magma chamber test model 10 1.4 Travel time delays for magma chamber model 11 2.1 Refraction of rays at an interface 19 2.2 Model parameterization 24 2.3 Inversion scheme 2.4 Ray diagram for resolution example 2.5 Singular values for resolution example 38 2.6 Singular vectors for resolution example 39 2.7 Resolution diagonals for resolution example 41 2.8 Resolution kernels for different damping factors 42 2.9 Resolution kernels for a reduced velocity gradient 43 flowchart 33 . 37 2.10 Resolution for various source-receiver geometries 45 2.11 Singular values for various source-receiver geometries 46 2.12 Velocity inversion example 48 2.13 Demonstration of damping effects 51 2.14 Removal of streaking effects 53 2.15 Interface inversion example 56 2.16 Combined velocity-interface inversion example 58 2.17 Observed-calculated amplitudes for attenuation example 60 2.18 Observed-calculated amplitudes for attenuation example 61 2.19 Attenuation inversion example 62 ix 3.1 Ray geometry 71 3.2 X versus 9 curve for a triplication 3.3 Neighbouring rays 82 4.1 Projected positions for the cross-ridge line 93 4.2 Comparison of cross-ridge bathymetry 95 4.3 Seismic sections for OBS 1 and OBS 2 96 4.4 Seismic sections for OBS 3 and OBS 4 97 4.5 Seismic sections for OBS 5 and OBS 6 98 4.6 Seismic sections for OBS 7 and OBS 8 99 4.7 Low velocity arrival branch 101 4.8 Shallow velocity solution 106 4.9 Shallow inversion residual map 107 76 4.10 Shallow resolution test 109 4.11 Deeper inversion residual map Ill 4.12 Complete velocity solution 113 4.13 First arrival power spectrum 116 4.14 Observed vs. calculated water wave amplitudes . 119 4.15 Smoothed models for amplitude calculation 121 4.16 Example travel time and amplitude curves 122 4.17 Observed vs. calculated amplitudes 1 124 4.18 Observed vs. calculated amplitudes II 125 4.19 Initial attenuation solution 127 4.20 Normalized amplitude map 128 4.21 Revised velocity solution 130 4.22 Final attenuation solution 131 5.1 136 Final velocity and attenuation models x 5.2 Comparison of reflection-refraction interface 2 138 5.3 Velocity vs. porosity curves 140 5.4 Velocity vs. pressure curves 142 A.l Travel time comparison for repositioning of OBS 3 169 A.2 Output vs. input voltages for the OBS/playback system 173 A.3 Comparison of compensated and uncompensated amplitudes 175 xi Acknowledgement I would like to sincerely thank Ron Clowes who originally suggested this research project, and has provided guidance and support in his capacity as thesis supervisor. Ron also took an active part in the planning and execution of the SEISRIDG-85 experiment. I am grateful to the many others who were involved in the SEISRIDG-85 experiment. Captain Alan Chamberlain and the officers and crew of the CSS Parize.au are acknowledged for their cooperation and assistance. Chief Al Woods, Les Rourke and Dan Desjardins of the Fleet Diving Unit, CFB Esquimalt, safely handled and deployed explosives during the experiment. Ben Ciammaichella, John Bennest and Bob Meldrum were primarily responsible for preparing and operating the 10 ocean bottom seismographs. Bob Meldrum and Ben Ciammaichella also spent much time developing an automatic digitization system for the SEISRIDG-85 seismic data. Andy Boland, Sonya Dehler and Chris Pike assisted in acquiring the data, and Connie Cudrak helped with digitization of the data. Bill Hill, Gail Jewsbury and Bob Macdonald (at the Pacific Geoscience Centre) were involved during preparation for the cruise. The people who comprise the Department of Geophysics and Astronomy at UBC have provided a stimulating and enjoyable environment within which to work. In particular, discussions with Andy Boland, Andy Calvert, Connie Cudrak, Tony Enders and Doug Oldenburg have all contributed to this thesis. I would like to thank the members of my thesis supervisory committee, Dick Chase, Ron Clowes, Doug Oldenburg and Matt Yedlin, as well as external examiner G.M. Purdy, for critically reviewing this thesis. Finally, and most importantly I thank my wife Ruth for the patience and endurance that she has shown over the past four and a half years. Her support, and the smiling faces of our children, Jordan, Calvin and Malcolm, have been a continuing source of inspiration. xii Financial assistance for the data acquisition (SEISPJDG-85) was provided by DSS Contract No. 11SB-23227-5-0209 (1985) through the Pacific Geoscience Centre (GSC). Research support was provided by EMR Research Agreement No. 206 (1987) through the Geological Survey of Canada and Infrastructure Grant No. 0855 and Operating Grant A7707, both from the Natural Sciences and Engineering Research Council, Canada. I was financially supported during this period by the Ph.D. Trainee Program of Energy, Mines and Resources, Canada. Xlll Chapter 1 Introduction The accretionary process by which oceanic lithosphere is formed at spreading centres continues to receive much attention. Central to this topic is the role that crustal magma chambers play. Ophiolite studies (Greenbaum, 1972; Dewey & Kidd, 1977; Pallister & Hcpson, 1981) as well as petrologic analysis of basalt samples (Bryan & Moore, 1977) suggest that crustal magma chambers are steady-state features which extend several kilometres vertically and 6-30 km laterally beneath the spreading axis. Seismic evidence from various segments of the slow spreading Mid-Atlantic Ridge indicate that an axial magma chamber is absent or at least very small, whereas a magma chamber has been observed seismically along different regions of the fast spreading East Pacific Rise, typically at 1.52.5 km depth. But, even where magma chambers are observed, they are generally quite narrow (~3-4 km) with maximum widths of < 8 km (for a review, see Macdonald, 1986 or Orcutt, 1987). The absence of an axial magma chamber in regions of slow spreading, in contrast to the presence of such a feature at fast spreading centres, suggests that the transient nature of a magma chamber is spreading rate dependent. Intermediate rate spreading centres tend to have properties similar to faster spreading centres (Macdonald, 1983), and thus might be expected to have an associated crustal magma chamber. The majority of seismic data collected to date in studies of young oceanic crust consists of multichannel reflection data or conventional seismic refraction data. The former provide good resolution of sharp contrasts within the crust and the lateral variability of these contrasts, thereby mapping structural geometry, but the vertical resolution is generally poor due to a lack of well defined velocity information. Conventional seismic 1 refraction profiles are adequate for determining pseudo-l * velocity structure, but are 1 Chapter 1. Introduction 2 insufficient for determination of velocity structure where there is substantial lateral variD ation. To define two-dimensional (2 ) velocity structure using seismic refraction data, a modified acquisition procedure is required which utilizes a multi-source, multi-receiver acquisition geometry. Such data can be interpreted using trial-and-error forward modelling methods (e.g. Cerveny, Molotkov, & Psencik, 1977; Spence, WTiittall, & Clowes, 1 but seismic tomography provides a more efficient approach. Using a high frequency approximation, both seismic travel time and amplitude can be expressed as projections (i.e. line integrals) of velocity and attenuation, respectively, and thus the velocity and attenuation of the medium can be reconstructed using travel time and amplitude data. In an attempt to better define lateral structural variations in the shallow crust beneath an intermediate rate spreading axis, and in particular to determine whether a crustal magma chamber is present, a multi-source, multi-receiver refraction profile was recorded across the Endeavour segment of the Juan de Fuca Ridge near 48° N, 129° W (see Figure 1.1). A tomographic inversion scheme has been developed and applied to the D travel time and amplitude data for this profile, to obtain a 2 velocity and attenuation model for the shallow crust. 1.1 E n d e a v o u r R i d g e S e g m e n t , J u a n de F u c a R i d g e The Juan de Fuca Ridge forms the accretionary boundary between the Pacific and Juan de Fuca plates, extending north from the Blanco Fracture Zone to the Sovanco Fracture Zone (see Figure 1.1). It consists of several en echelon spreading segments and has an intermediate full spreading rate of 5.8 cm/yr (Riddihough, 1984)- The Endeavour Ridge segment is 90 km long. It overlaps adjacent spreading segments to the north and south at the Endeavour Offset and Cobb Offset, respectively. A bathymetric profile along the length of the Endeavour Ridge segment (Figure 1.2 (b)) reveals an arcuate shape with a minimum depth reached near 47° 58' N. In this region the spreading centre is characterized by a well-developed axial ridge which gradually disappears into broad axial Chapter 1. Introduction 3 Figure 1.1. Bathymetric and profile location map in the region of northern Juan de Fuca Ridge. Solid dots are the sites of the 8 OBSs; the solid line is the airgun profile shot into the OBSs. The dashed line is the location of a multichannel reflection profile. The stars identify two hydrothermal vent sites along the Juan de Fuca Ridge. Bathymetry map is based on Sea Beam and other data (Geological Survey of Canada, 1987); detailed contours north of 47° 56' and surrounding the ridge are at 20 m intervals, otherwise 50 m contours are plotted. The inset shows the regional plate tectonic regime with the area of the bathymetric map indicated by the rectangle. Chapter 1. Introduction q w CN E Axial ridge Ayv 2 /\V v - £ 4 5 3 6 7 ^ ^ ^ ^ CL Q /\1 / V.E. 25:1 CN 8 •OBS i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—r 0 10 20 30 Position (km) q CN s SEISRIDG-85 Line N £ / Endeavour Ridge CL Q South Endeavour Valley q V.E. North N Endeavour Valley 30:1 20 40 60 Position (km) Figure 1.2. Endeavour Ridge bathymetry, (a) Cross-strike bathymetric profile with OBS positions marked. The positions indicated are measured from the western end of the refraction line (Figure 1.1). (b) Along-strike bathymetric profile. The positions indicated are measured moving north along the ridge from 47° 41.5' N , 129° 17.0' W. Note that the vertical and horizontal scales are different in (a) and (b). Chapter 1. Introduction 5 grabens to the north and south. The morphology of Endeavour Ridge near 47° 58' N (Figure 1.2 (a)) is similar to that of the East Pacific Rise near 13° N which has a fast full spreading rate of 11-12 cm/yr (Francheteau & Ballard, 1983). Endeavour Ridge is built on a broad plateau suggesting an enhanced magma supply in this region for the last 1 m.y. (Delaney, Karsten & Hammond, 1986). Basalts from Endeavour Ridge are potassium-enriched, moderately fractionated MORBs (Karsten & Delaney, 1987 ), distinguishing them from the low potassium basalts south of the Cobb Offset and indicating that two distinct magma sources are associated with these two spreading segments (Delaney et al., 1986). The axial ridge is approximately 4 km wide at its base and rises 200-250 m above the surrounding seafioor (Figure 1.2 (a)). The summit depression is < 1 km wide and is 75-150 m deep. Analyses of submersible map- ping (MERGE Group, 1984; Tivey & Delaney, 1986; Tivey & Johnson, 1987), deep-tow photography, Sea MARC I and Sea MARC II side-scan sonar data (Kappel & Ryan, 1986; Karsten et ai, 1986) indicate that the flanks of the axial ridge consist of unfractured pillow flows with interpillow sedimentation, which predate the lightly sedimented pillows and lobate flows of the summit depression. Faulting and fracturing are limited to the summit depression and the seafioor at the base of the axial ridge. The lack of Assuring and faulting on the flanks of the axial ridge suggest that it is a volcanic construction, and that the summit depression was formed during a stage when crustal extension occurred with little associated volcanism (Kappel & Ryan, 1986; Tivey & Johnson, 1987). According to this model, the smaller ridges paralleling the axial ridge are the remnants of a former axial ridge that were rafted away intact. There is substantial evidence indicating that the region near 47° 58' N is the site of most recent magmatic injection along the Endeavour Ridge segment, and as such is the most likely location of a crustal magma chamber. Macdonald, Sempere & Fox (1984) postulate that a long wavelength axial high (see Figure 1.2 (b)) between overlapping spreading centres is generally found over the region where a magma diapir intrudes the Chapter 1. Introduction 6 crust. This high is distinct from that due to volcanic construction. Kappel & Ryan (1986) suggest that the plan-view bow-form shape which characterizes the axial volcanic edifice on the Endeavour Ridge segment is to be expected when magma enters the crust at a discrete centre and flows along-axis from that point. A hydrothermal plume extending from the seafioor to 300 m above the seafioor, with water temperatures elevated from the ambient temperatures by 0.10-0.13° C, has been observed for distances of 10 km along the ridge axis in this region (Baker & Massoth, 1985; Crane et al, 1985), and water temperatures of 400° C were recorded at two hydrothermal vent sites along the axial valley wall at 47° 58' N (Delaney, McDuff & Lupton, 1984). Franklin & Kappel (1986) require the presence of a heat source, sufficiently large to sustain a hydrothermal system at 400° C, to account for the formation of the large sulfide deposits that characterize the Endeavour Ridge, and they suggest that there has been a recent replenishment of the crustal magma chamber there. If this is the case, there has been no associated volcanism as Davis et al. (1984) argue that sediment cover on axial valley flows dates them at ~ 10,000 yr. This is in contrast to the Southern Juan de Fuca Ridge where very young volcanics (< 2000 yr) have been identified on the floor of the axial depression (Normark et al., 1982; Normark et al., 1983). Also, the high fluid lava to pillow flow ratio associated with a replenished magma chamber (Francheteau & Ballard, 1983) is not observed on the Endeavour Ridge. In summary, these data suggest that the along-axis high near 47° 58' N is the site of most recent crustal magmatic injection along the Endeavour Ridge segment. A heat source is currently present, but whether a magma chamber remains or has been replenished is unknown. The lack of observed recent volcanism weighs against large scale replenishment. There are several previous seismic studies which provide information about the crustal structure of the Juan de Fuca Ridge. McClain & Lewis (1980) describe results from an experiment located south of the Cobb Offset near 47° 00' N on the Northern Symmetrical Ridge, where the ridge morphology is similar to that of Endeavour Ridge. They find Chapter 1. Introduction 7 no evidence for the existence of a crustal magma chamber or lateral velocity variation beneath the ridge from single-channel, deep-tow reflection data across the ridge, or from reversed refraction lines along strike and normal to the ridge (shot spacing of ~1 km for distances < 30 km). Morton (1984), interpreting multichannel reflection data (24-fold, CDP points at 50 m intervals) collected across Southern Juan de Fuca Pudge, identifies a short reflector at ~1.0 s two-way travel time sub-bottom as a reflection from the roof of a magma chamber. A similar reflection at ~1.0 s sub-bottom is observed by Rohr, Milhere.it & Yorath (1988) on a multichannel reflection line (30-fold, CDP points at 12.5 m intervals) across Endeavour Ridge. As an alternative to the magma chamber interpretation, they suggest that the reflection may be caused by a region of high vertical temperature gradient. Rohr et al. (1988) also observed a shallow intermittent reflector at 0.30-0.55 s sub-bottom which they associate with either the pillow basalt-diabase dike contact or a metamorphic front within the pillow basalts. A deeper reflector is observed which dips from 1.4 s sub-bottom 3-5 km east of the axis, to Moho depths at 12 km distance. These observations indicate that there are significant lateral variations in crustal structure beneath Endeavour Ridge. 1.2 The SEISRIDG-85 experiment In order to understand further the variations in crustal structure, we planned a detailed seismic refraction survey, SEISRIDG-85, which was conducted in November of 1985 on the Endeavour Ridge segment of the Juan de Fuca Ridge near 48° N, 129° W (Figure 1.1). The program consisted of two phases, each phase designed to address one of the following objectives: (1) definition of the shallow crustal structure and imaging of a crustal magma chamber (if one exists), and (2) determination of the regional crustal structure to Moho depths. Objective (2) was accommodated in the initial phase of the experiment using 90 kg explosive charges recorded by an areal array of 10 ocean bottom seismographs Chapter 1. Introduction 8 (OBS), with maximum dimensions of 60 km. In the second phase of the experiment, designed to satisfy objective (1), the 10 OBSs formed a smaller array (maximum dimensions of 20 km) recording airgun shots. Cudrak (1988) has analyzed the reversed airgun lines to determine shallow structure, and analysis of the explosion lines (for deeper structure) is underway. In this analysis, we only consider a single refraction line (referred to as the crossridge line) which is centred on 47° 59' N , 129° 05' W, normal to Endeavour Ridge, and coincident with the multichannel reflection line of Rohr et al. (1988); see Figure 1.1. Eight OBSs were deployed over a distance of 20 km along the cross-ridge line, to record airgun shots fired at 200 m intervals for 30 km along the line. The closely spaced array of receivers was designed to allow definition of lateral structural variations beneath the ridge, and to provide multipath ray coverage at a depth of 2-3 km beneath the ridge where a crustal magma chamber might exist. 1.3 M a g m a chamber detection Detection of a magma chamber using seismic refraction techniques relies primarily on the velocity contrast of the low velocity magma chamber relative to the surrounding crustal material. The expected velocity contrast is large whether the chamber contains molten magma (2.0-2.5 km/s, Murase & McBirney, 1973) or cumulate mush (4.5-5.0 km/s), as crustal velocities in the associated depth range are typically > 6.0 km/s. The magma chamber will cause travel time delays for rays which encounter it, and possibly shadow zones. Replacing 6.0 km/s material with 4.5-2.5 km/s material results in delays of 55230 ms/km of travel path through the, low velocity zone, producing expected travel time delays of 200-900 ms for a 4 km wide magma chamber. The manifestation in seismic data of such a low velocity zone may be more subtle. Travel time residuals associated with the 4.5 km/s low velocity zone shown in Figure 1.3, Chapter 1. Introduction 9 for refracted arrivals corresponding to the source-receiver geometry used for the crossridge line, are shown in Figure 1.4 (a). The low velocity zone extends 1 km in depth, widening from 1 km at the top to 3 km at its base, and is imbedded in a layer with a velocity of 5.8 km/s at the interface above the magma chamber and a vertical velocity _1 gradient of 0.5 s . Thus the low velocity zone has velocities reduced by < 2.0 km/s relative to the background model. It is surrounded by a high gradient transition zone which results from specifying the model on a discrete grid. Maximum travel time residuals are < 180 ms, and the majority of the residuals are < 100 ms (Figure 1.4 (a)). Referring to the ray diagram (Figure 1.3) the first arrivals which encounter the low velocity zone are geometric arrivals which emulate diffractions; it is faster to travel around the magma chamber than through it. These geometric arrivals are analagous to diffractions, as they are bent around the corners of the magma chamber in the high velocity gradient transition region. If the boundary were replaced by a velocity discontinuity, these geometric arrivals would not be observed. As they correspond to diffractions, the amplitudes should be small. Burnett, Orcutt & Olson (1988) used a finite element algorithm for seismic wave propagation to show that very low energy arrivals are associated with energy diffracted around a 4 km wide magma chamber. In the example presented here, rays which travel through the magma chamber (not shown in Figure 1.3) are secondary arrivals and have greater delay times (similar to the large delays of 200-900 ms initially expected), but their travel path through the low velocity zone and subsequent arrival position at the surface is strongly dependent on the shape of the magma chamber. Thus, the effect of the magma chamber is to produce first arrival travel time delays of < 180 ms as well as regions of reduced or negligible amplitude. Secondary arrivals due to ray paths travelling directly through the magma chamber would in general not be identifiable, due to the reverberatory nature of real data, and thus would not contribute to the subsequent interpretation. Figure 1.4 (b) shows travel time residuals associated with a 1 x 1 km low velocity zone. 10 Chapter 1. Introduction Position (km) Figure 1.3. Magma chamber test model. A magma chamber, 1 km wide at the top and 3 km wide at the base, extending 1 km in depth is included in an otherwise layered model. The magma chamber has a velocity of 4.5 km/s compared to velocities of 6.0-6.5 km/s in adjacent regions of the model. A high velocity gradient transition zone surrounds the low velocity region as a result of specifying the model on a discrete grid. The first arrival rays for OBS 3 are shown. Ray paths through the water column have been excluded, here and in all subsequent ray diagrams. 11 Chapter 1. Introduction 8 7 6 5 ml iiiiiillMllJllllllllll — A llllllllllliiinr w 300 1 3 2 1 lli llllllllllllllllll nm .Jilllii..... 8 7 6 5 i 300 niiiiii A II 3 2 1 n 1111 it 111 • •• llllll... i—i—i—i—i—i—i—i—r 5 T—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—r—r 10 15 20 25 30 Position (km) Figure 1.4. Travel time delays for magma chamber model, (a) First arrival travel time residuals for model in Figure 1.3. Residuals are travel times for the layered model subtracted from those calculated for the magma chamber model. The positive residuals indicate delays associated with the magma chamber. The OBS numbers are shown on the left and OBS positions are indicated by triangles, (b) Travel time residuals for a l x l km magma chamber. Chapter 1. Introduction 12 They are limited to < 80 ms. As in the previous example, first arrivals which encounter the low velocity zone correspond to diffracted arrivals, although the diffractions in this case are not as severe. A wavelength of 0.5-0.6 km (for v—Q.0-7.0 km/s) is associated with a 12 Hz signal in the expected depth range of a magma chamber. As the dimensions of the velocity anomaly approach the wavelength of the seismic signal, the regions where travel time delays and reduced amplitudes are observed will be less distinct. Because of this, as well as the smaller travel time delays associated with a zone of this size, we estimate that a 1 x 1 km zone of partial melt represents the minimum dimensions of such a feature that would be clearly detected in this refraction experiment. 1.4 Thesis outline The remainder of this thesis is concerned with the development and application of an inversion method to determine the 2 velocity and attenuation structure for the cross-ridge D profile, as well as an interpretation of the results. In Chapter 2, the travel time-velocity and amplitude-attenuation problems are formulated in a manner suitable for inversion. The complete inversion scheme is then presented followed by a discussion of resolution and synthetic inversion examples. Chapter 3 describes explicit formulae required to construct the linear systems that are inverted in Chapter 2. These include analytic formulae for dynamic ray tracing and attenuation calculation, as well as inversion formulae. Chapter 4 describes application of the inversion procedure to data from the cross-ridge line, and the resulting velocity and attenuation models are presented. Finally, these results are discussed and compared with the coincident reflection line of Rohr et al. (1988) in Chapter 5. Chapter 2 D 2 2.1 Travel Time & Amplitude Tomography Introduction Inversion methods which can generally be referred to as tomographic reconstruction techniques have recently found various seismic applications. The most common consideration or application of these techniques has been to velocity inversion using seismic travel times obtained from borehole-borehole and/or borehole-surface measurements (e.g. Bois et ai, 1972; Wong & West, 1983; Menke, 1984a; Peterson, Paulsso McEvilly, 1985; Ivansson, 1985;« Radcliff, 1984; Lytle & Dines, 1980; Gustavsson 1986; Bregman, 1986). Amplitude data from borehole experiments have been used to invert for attenuation (Wong & West, 1983; Bregman, 1986). Other applications include 2 D velocity and reflector depth determination using multi-offset reflection data (e.g. Bishop et al, 1985; Fawcett & Clayton, 1984; Stork & Clayton, 1986), inv D sion of earthquake travel time delays for 3 velocity structure (e.g. Aki, Christoffers- son & Husebye, 1977; Humphreys, Clayton & Hager, 1984; Nakanishi, 1985), i D sion for 2 refractor velocity and/or refractor topography (e.g. Hearn & Clayton, 1986; D Kanasewich & Chiu, 1985; Pavlis, 1986) and determination of 2 velocity structure from surface—surface refraction data (e.g. Firbas, 1987; Chapman, 1987). We are interested in the determination of 2 D velocity and attenuation structure using surface-to-surface refraction data. An iterative tomographic inversion scheme has been developed to determine the D 2 velocity and attenuation structure of a medium from first arrival travel times and amplitudes. Specifically, the method is developed for surface-to-surface source-receiver 13 D Chapter 2. 2 Travel Time 8c Amplitude Tomography 14 geometries typically used in refraction surveys, but for cases where there is a higher density of sources and receivers than in conventional surveys. Because this method is based on turning rays, ray tracing constitutes an integral part of the procedure, unlike other methods (i.e. borehole-borehole or multi-offset reflection methods) where straight ray approximations are adequate. A linearized formulation of the travel time-velocity problem will be presented for a continuous velocity medium as well as for the case where an interface is present, followed by a formulation of the amplitude-attenuation problem. The model parameterization and the basis functions used for interpolation are then described, followed by a discussion of the linear inverse problem. The method for solving the linearized problem and its incorporation into an iterative inversion scheme are outlined, and finally, the resolution characteristics of this method are considered and the inversion method is demonstrated with various synthetic examples. All quantities required to construct the linearized travel time-velocity and amplitude-attenuation systems for inversion are determined by ray tracing. Details of the ray tracing algorithm are postponed until Chapter 3. 2.2 Problem Formulation 2.2.1 Velocity Formulation Our present objective is to determine the velocity structure of a medium using seismic travel times. Because the relationship between these two quantities is nonlinear, we seek a linearized formulation which may be used in an iterative fashion to achieve this goal. The travel time along a ray in a continuous medium whose velocity is ^o(x) is (2.1) where L represents the ray path. If the velocity of the medium is perturbed by fo(x) 0 Chapter 2. 2 Travel Time <fe Amplitude Tomography D 15 such that v(x) = v0(x) + 8v(x.), then the perturbed travel time is (2.2) where L is the new ray path. The travel time through the perturbed medium can be expressed as an expansion about i>o(x) as t{v(x)) = t(v (x)) + ^ 8v(x) + 0 dt dL dL dv 2 6v(x) + 0(8 v(x)). (2.3) v=v v=v v=v 0 0 0 Using Fermat's principle which states that the travel time for the ray path is stationary i.e. at_\ ah j v=v 2 0), and discarding terms of 0(S v(x)) we obtain 0 St = t{v(x)) - t(v (x)) 0 dt 8v(x). dv v= v (2.4) 0 Differentiating (2.1) and substituting the result in (2.4), St _ Sv{x) dv dl. JL„ vl r (2.5) This linear relationship between St and 8v(x) has been used by various authors (cf. Ju- lian & Anderson, 1968: Backus & Gilbert, 1969; Dziewonski & Gilbert, 1976). problem has been formulated here in terms of velocity instead of slowness, as analytic expressions used to evaluate (2.5) are simpler for a model parameterized using constant velocity gradient cells as compared to one utilizing constant slowness gradient cells. Because of the different source-receiver geometries encountered in seismic refraction experiments, we require a flexible formulation. For this reason, we choose to discretize the velocity model at this point, as opposed to seeking an analytic inversion formula for 8v(x) based on some specific recording geometry. If the velocity model is specified at the Chapter 2. 2 Travel Time & Amplitude Tomography nodes of a rectangular grid by N discrete values Vj, 16 j = 1,...,N then the velocity field can be obtained using an interpolation scheme of the form (2-6) «(*) = E where the /3j(x) form a set of basis functions over the region of model definition. The form of the /3j(x) depends on the method of interpolation used (for details of the interpolation method employed, see Section 2.2.5). The velocity perturbation becomes 6v{ )=£,l3 (x)8v X i (2.7) j which when substituted into (2.5) results in Alternatively, the differential St is by definition N FT = dt ( - ) E « r ^ . 2 9 • - OVA and thus by comparison with (2.8) we recognize that IT = - / ^ 7 % ( ' » 2 10 For a set of M rays a linear system of equations is formed from (2.8) St = ASv (2.11) where St is the vector of travel time residuals (t (, „ £ — t i i i), 0 geP ec matrix J ^ , and Sv is the velocity perturbation vector. ca cu ate( A is the Jacobian To form (2.11) we must be Chapter 2. 2 Travei Time & Amplitude Tomography 17 able to calculate the travel times and partial derivatives for the set of rays using (2.1) and (2.10), respectively. The ray tracing algorithm used to evaluate these quantities is described in Chapter 3. Due to the linearization, both of these expressions are evaluated along the unperturbed ray path using the unperturbed velocity field. Thus linearization has provided a system which linearly relates travel time perturbations to changes in the model, and allows direct calculation of this linear system. This formulation will be referred to as method A . An alternate approach to the problem is one which is more akin to Series Expansion Techniques (SET) commonly used in medical imaging (see Censor, 1983). If we introduce -1 the wave slowness u(x) = v (x), then (2.1) becomes t(u (x)) = / uo(x)dl. 0 (2.12a) If the medium is discretized into a grid of /V cells, each with an associated slowness value Uj, j = 1 , N , then (2.12a) can be approximated by N t = ^2ujdlj, (2.126) where dlj is the path length of the ray in the jth cell. This corresponds to an expansion of the slowness field as u(x) — Y,f=ifij{ ) jJ where the basis functions are defined by x { u 1 if x is within the jth cell, 0 otherwise. (2.13) In this case the slowness values are specified for each cell, in contrast to method A where velocities are associated with the grid nodes. For a set of M rays we form the linear system t = Au (2.14) D Chapter2. 2 Travel Time & Amplitude Tomography 18 where t = vector of observed travel times, A = dlij — matrix of partial path lengths and u = vector of slowness values. Again, to form (2.14) we need to be able to calculate the travel times t as well as the partial path length matrix A . The element A^ is the path length of the ith ray in the jth cell. This will be referred to as method B. 2.2.2 Interface Formulation As well as determining velocity perturbations within a continuous medium, we would also like to be able to determine the depth and geometry of interfaces in a model across which the velocity is discontinuous. When considering the change of travel time caused by perturbations of interface depth, the associated change in ray path must be considered. To facilitate this, the problem can be formulated in the following way. A ray travelling from S to D is refracted at depth z by an interface separating media with velocities v\ x and i> (see Figure 2.1). The unit interface normal 1 is at an angle 4> relative to the 2 n vertical. If the depth of the interface is increased by Az, then another ray travelling from S to D is refracted at point Y on the new interface and travels from X to D along the same path as the first ray. Because the travel times to W (indicated by the plane wavefront) are the same, as are the travel times along the coincident segment from X to D , the change in travel time due to the interface depth change is the travel time difference of the rays from W to X . This can be expressed as (2.15a) (2.156) The unit vectors ni and h are the tangents to ray 1 at the interface in mediums 1 and 2 2, respectively, and ei and e are unit vectors along the x and z axes. Obviously this 3 formulation is strictly valid for plane waves impinging on a plane interface separating constant velocity media, but locally it provides a useful approximation which can be Chapter 2. 2 D Travel Time & Amplitude Tomography 19 D 1 / t " ~» _ Az \ \ \ Y *""**"" --.._ * i \A/=wavefront r e >2 Figure 2.1. Ray geometry for refraction at an interface. The original interface (solid line) is at depth zj, and the perturbed interface (dashed line) is at depth z . Symbols are defined in the text. 2 D Chapter 2. 2 Travel Time &: Amplitude Tomography 20 used in an inhomogeneous medium. If the interface is discretized into N segments each with an associated depth Zj, then for a set of M rays a linear system can be formed St = B6z (2.16) relating the travel time residuals St to the interface depth changes Sz. B is the Jacobian matrix J - ^ which can be approximated by (2.15). 2.2.3 Attenuation Formulation In general, amplitude variations of a propagating seismic wave are caused by elastic and anelastic effects. Elastic effects such as geometrical spreading and energy partitioning at interfaces are controlled by the large-scale velocity structure of the medium, whereas anelastic attenuation is due to internal friction. In addition, anelastic attenuation is mimicked by scattering caused by small-scale (distances of < a wavelength) velocity fluc- tuations (Schoenberger & Levin, 1974; Richards & Menke, 1983). Anelastic attenuatio and small-scale scattering are described collectively by the effective attenuation a(x) which is defined as «> = ^ = - ; K ' (2 17) where <5(x) is the effective quality factor, U(u>,x) is the wave amplitude for angular frequency UJ and A[/(u>,x) is the change in amplitude per cycle due to the effective attenuation. It is assumed that cc(x) is constant over the frequency range of interest, and that the medium is non-dispersive (i.e. v(tv,x) — v(x)). The definition (2.17) can be used to obtain an explicit expression for U in terms of a. The wave amplitude at position I + 81 (I is the path length) along a ray is 2 U(l + 81) = U(l) + —81 + 0(8l ). (2.18) D Chapter 2. 2 Travel Time & Amplitude Tomography Assuming that the amplitude varies slowly along the ray (i.e. ^ 21 << 1), then SU = U(l + 8l)-U(l) = —81. (2.19) Using (2.19) in (2.17) and A = 2m>/w, to a(x) SU ~W SI ~2~v~(xj = (2.20) which when integrated yields, U(u>,x) = r/ ( x)exp B Wl J ^ d / J (2.21) L where C/#(x) is the elastic wave amplitude incorporating geometrical spreading and the effect of interfaces. Equation (2.21) is utilized in the spectral ratio or frequency ratio method to determine the integrated attenuation from the observed amplitude spectrum (e.g. Bath, 1974, P- 334; Trappe, 1988). For a narrow-band wave source with dominant frequency w , the attenuation factor in (2.21) can be approximated by using u for all c c frequencies about iv . Then, c U'(x) = I ^\dl JL v(x) ' (2.22) where the normalized amplitude U'(x) is U V )= l (M^l). ( . in UJ C \U(UJ,X) 2 23) J If f(x) is known (and hence L is known), (2.22) provides a linear relationship between the amplitude U'(x) and the attenuation a(x). Defining the attenuation by N discrete values ctj, j — l,...,N, then a(x) is obtained using an interpolation formula «M = £ft(x)«j. (*24) D Chapter 2. 2 Travel Time &c Amplitude Tomography 22 where the basis functions /3j(x) are the same as those used in (2.6). Analagous to the method A velocity formulation, we obtain the following linear system for a set of M rays: U' = A a (2.25) where, (2.26) The vector U ' contains the M observed amplitudes (see (2.23)) and a contains the N attenuation values to be determined. This formulation can be repeated using an approach analagous to method B applied to the velocities (see Section 2.2.1). In contrast to the travel time-velocity relationship, the amplitude-attenuation relation is linear. Once the velocity v(x) is determined for the structure, the attenuation a is obtained directly from the amplitudes U ' ; no iteration is required. To form (2.25), A must be calculated using (2.26), and the amplitudes U ' must be formed using the observed amplitudes U and the amplitudes U^; which account for elastic propagation in the medium v(x). The details of calculating UE and. A by ray tracing are provided in Sections 3.4 and 3.6.2, respectively. 2.2.4 Model Parameterization The forward problem formulation requires specification of the model in terms of a discrete set of parameters; velocities or slownesses, and attenuations. This discretization must be flexible enough to allow adequate representation of different types of models. In addition, we require a model consisting of constant velocity gradient cells allowing use of analytic ray tracing formulae (see Section 3.3). The most obvious model parameterization accomD D modating these requirements is a 2 (3 ) rectangular grid where the model parameters are specified at the nodes. For the purpose of ray tracing the rectangular cells of this grid can be divided into triangular (tetrahedral) subcells within which the velocity and D Chapter 2. 2 Travel Time &c Amplitude Tomography 23 attenuation gradients are constant, and are determined uniquely from the values at the vertices of the triangle (tetrahedron). This parameterization is well suited to method A, as the velocity perturbations of equation (2.11) are determined at the nodes of the grid. However, it is not immediately compatible with method B. In this latter formulation a constant slowness is assigned to each cell of the grid. As this slowness value represents an average value for the cell it is most appropriately assigned to the centre of the cell. Thus the slowness values in (2.14) are determined at the cell centres, and an intermediate step of interpolation/extrapolation is required to define the velocity values at the nodes of the grid for ray tracing. An average of nearest neighbours is used for interpolation. D Interfaces included in the 2 model are represented using cubic splines as opposed to a piecewise linear representation to eliminate erratic ray behaviour near the vertices of linear segments. 2.2.5 Basis functions for Method A D The 2 parameterization used for ray tracing is illustrated in Figure 2.2. Velocities are specified at the nodes of the grid, and the velocity within each triangular cell is given by v(yz) = v + V v(x - x ) + V v(z - z ), 0 where v , V i> and 0 x V i> z x 0 z (2.27) 0 are determined using the velocities at the vertices of the triangle. For example, within triangle 1, v —Vi, 0 V v=(v x 2 — Vi)/Ax and V v—(v — v )/Az, where z 5 2 Ax and Az are the horizontal and vertical cell dimensions, respectively. The model attenuation is specified in an analagous manner. Using this method of interpolation, the basis functions in (2.6) and (2.24) are piecewise continuous functions (generally discontinuous across cell boundaries) each of which is zero everywhere except in a maximum of 6 triangles of the model. For instance, the node where v is defined is a vertex 5 for 6 triangles. Thus, 3 (x.) is nonzero within these 6 triangles, and is zero everywhere 5 Figure 2.2. Velocity parameterization grid. The velocities associated with each node are labelled, as well as the triangular cell numbers. D Chapter 2. 2 Travel Time Si Amplitude Tomography 25 in triangles 2 and 7. The basis function 8 (x) is defined by 5 (z — z )lAZ in cell 1, (x — xo)/Ax in cell 3, 0 (z - z )/Az 0 B (x) = \ (x - x )/Ax 0 5 + ((x + Ax) 0 - x)lAx - f ( ( 2 + Az) - z)/Az 0 in cell 4, in cell 5, ((x + Ax) - x)/Ax + 1 in cell 6, ((20 + Az) - z)/Az + 1 in cell 8, 0 I0 Contributions to the Jacobian (2.28) in cells 2,7. in (2.10) or to the attenuation matrix in (2.26) will be limited to segments of the ray path which encounter triangles within which f3j(x) is nonzero. 2.3 Linear Inversion The result of ray tracing is a linear system of equations represented as Ax (2.29) This corresponds to (2.11), (2.14), (2.16) or (2.25). The matrix A is dimensioned M x N where M—number of data, and A^=number of parameters. Typically this system is sparse (i.e. most elements of A are zero) and inconsistent (i.e. no exact solution exists). In the case of (2.11), (2.14) or (2.25) it is most commonly underdetermined and overconstrained (i.e. M < N, and M > rank(A)) while in (2.16) it is normally overconstrained. To compensate for the reliability of the data measurements, (2.29) can be multiplied by the data covariance matrix Q (Wiggins, 1972). If the measurement uncertainties are statistically independent (we assume this condition), then Q = diag {^-}- The standard deviations <Tj referred to here are the estimated measurement uncertainties. In what Chapter 2. 2 Travel Time &: Amplitude Tomography 26 follows it is assumed that (2.29) has been multiplied by Q. A common approach to solving such a system is to utilize the generalized inverse of A (Lanczos, 1961; Wiggins, 1972; van der Sluis & van der Vorst, 1987). Any such matrix can be written using the singular value decomposition (SVD) of A as A = UAV r (2.30) where U and V are orthogonal matrices of dimension M x M and N x N respectively, and A is an M X N diagonal matrix with diagonal elements Ai > A > ... > A > 0, 2 r + r < min(M, N). The generalized inverse of A denoted by A is T A+ = V A + U r r , (2.31) from which we obtain the particular solution + x = A b. (2.32) The subscripts r in (2.31) indicate that only those columns of U and V corresponding to nonzero singular values Aj are included in U and V . Likewise, A+ contains only r r nonzero singular values. Since b generally lies in the M-dimensional data space, A + projects b onto the column space of A (i.e. minimizing | | A x — b||) and then maps this point into x , lying in the row space of A (a subspace of the ^-dimensional model space). Thus, the particular solution x has no component orthogonal to the row space of A , and the only nonzero contribution to the prediction error e = A x — b comes from that component of b orthogonal to the column space of A . A and A + provide a one-to-one mapping between the column space and the row space of A , and thus x can always be found. The particular solution x reduces to the smallest solution when the problem is Chapter 2. 2 Travel Time <fc Amplitude Tomography 27 purely underdetermined, and to the least-squares solution when the problem is purely overconstrained (Menke, 1984b). The approximate solution can be expressed in terms of the true model as x = V V r T r x. (2.33) Wiggins (1972) defines the resolution matrix R = V V r T r , since the resemblance of R to the identity matrix I determines how well x resolves the true model x. The jth row of R (or column since R = R ) indicates the resolution of the jth model parameter T as it is the least-squares solution for maximizing the jth parameter. Perfect resolution of the jth parameter is achieved if the jth row contains zeros except for a 1 in the jth column, and the resolution is compact if the nonzero elements correspond to parameters physically adjacent to the jth parameter. If the rows of R are characterized by compact resolution, Wiggins (1972) points out that the diagonal elements of R alone can be used to characterize the resolution of the system. The resolution compactness about some point in the model is.indicated by the distance spread of adjacent model parameters whose associated diagonal elements in R sum to 1. This allows a concise characterization of the resolution. Assuming that the data errors are uncorrelated, have zero mean and have equal standard deviation of 1.0 (due to the assumed normalization), then the variance of the jth parameter estimate is ™=YMt=l 2.3.1 (2-34) i Regularization Because the particular solution lies in the row space of A , which is spanned by the columns of V , the solution can be thought of as a linear combination of the column r vectors V j . The weighting coefficient for the jth singiriar vector V j in the solution is (2.35) Chapter 2. 2 D Travel Time & Amplitude Tomography 28 1 Any errors in the data b will be amplified by Xj and thus the contribution to the solution of those singular vectors with corresponding small singular values will be erroneously large, and possibly of the wrong sign. To stabilize the solution when the data contain noise, the problem must be regularized. This can be accomplished by winnowing of singular vectors in the solution according to some cut-off criterion (Wiggins, 1972), or by damping the solution (Marquardt, 1970). In the latter case, (2.35) is replaced by w = = T^T^ U J • B 2 3 6 < - ) The damping factor u hmits the size of the solution T x = V A;U b, r (2.37) r and the resolution matrix becomes R = V A;V r where A' — A ( A T r 2 r 2 1 4- a !)' . T (2.38) r The solution (2.37) is the solution of the regularized normal equations ( A A + yu. I)x = A b T 2 T (2.39) which are obtained by minimizing T T </, = e - t V ( x x ) . e (2.40) It is referred to as the ridge regression estimate or the damped least-squares solution (Menke, 1984b). D Chapter 2. 2 Travel Time & Amplitude Tomography 29 More generally, we can minimize the objective function 2 T <t> = e e + / x ( W x ) W x T (2.41) where W is some weighting matrix. Formally it is referred to as the parameter covariance matrix (Wiggins, 1972), but its form can be arbitrarily chosen to obtain different classes of solution models. Postmultiplying A in (2.29) by W _ 1 W , the approximate solution (2.37) becomes x = W^VrA^U/b (2.42) where it is understood that U and V contain the left and right singular vectors of r r 1 AW" . For example, consider the linear system (2.16) relating interface depth changes to travel time changes. Instead of determining the smallest solution, we may want to find the flattest solution; i.e. the solution that minimizes ||6\.z||. This is conveniently done by using a backward finite difference approximation for the derivative, which results in the weighting matrix W and the inverse W w / 1 0 -1 1 = _ 1 0 / i w 1 -1 (2.43) 1 1 - 1 1 1 - 1 1 V o -1 1/ Vl 1 • 1 1 1/ It is important to realize that postmultiplying A will result in a product matrix which is no longer necessarily sparse. Even when W has a narrow bandwidth, as it does in the above example, the inverse W product matrix A W - 1 _ 1 may have a broad bandwidth which will result in a that is fuller than A . This is an important consideration as it will increase storage requirements for the matrix as well as requiring solution of a nonsparse D Chapter 2. 2 Travel Time Sz Amplitude Tomography 30 system. 2.3.2 Global Norms When system (2.11) is solved by minimizing the objective function in (2.41) the norm is applied to the current perturbation. In other words, the solution minimizes the perturbation relative to the previous solution. Instead, it may be more desirable to construct models which minimize a global norm. After Oldenburg (1983) this can be accomplished as follows. For the kth iteration the linear system of equations (2.11) is (2.44) St = Ak^vjj where £v k = v f c + 1 - (2-45) v . fc This can be rewritten as, (2-46) St' = AfciSvk where, St' = new data = St -f Ak(vk — v ) 0 (2-47) and, £v ' k = v f c + i - (2.48) v . 0 The objective function to minimize is then < M e ' V + /i (wK) wK 2 r (2.49) where e' = St' — A^Sv^'. The solution so constructed will be referred to as the globally smallest or flattest deviatoric model as opposed to the smallest or flattest deviatoric model as is appropriate. The globally smallest solution minimizes the deviation of the Chapter 2. 2 Travel Time &: Amplitude Tomography 31 solution from the starting model. Systems (2.14) and (2.16) can be recast in a similar fashion. 2.3.3 Solution Method. Singular value decomposition facilitates construction of the generalized inverse as well as resolution and variance estimates for the resulting solution. It will be used to study resolution characteristics for small problems where the matrix A has dimensions of 2 2 2 (3(10 x 10 ). But, SVD requires computing time of 0(MN ) and memory space of O(MN) (van der Sluis & van van der Vorsi, 1987), and thus is unsuitable for solution of larger systems. To avoid this problem, iterative solution techniques are commonly employed to solve linear systems in tomography. The most common of these methods are referred to as row action methods of which Stationary Iterative Reconstruction Techniques (SIRT) and Algebraic Reconstruction Techniques (ART) are examples (see Censor (1981) for a review). Conjugate gradient methods provide an alternative to row action methods, at the cost of a meagre increase in storage requirement. Ivansson (1986) compared ART, SIRT and a conjugate gradient algorithm and found that the conjugate gradient method converged more rapidly. To solve the regularized normal equations (2.29) we employ a modified conjugate gradient algorithm (subroutine LSQR) developed by Paige & Saunders (1982). It is particularly suited to solving large sparse asymmetric systems. Their method is analytically equivalent to other conjugate gradient methods except that it is more stable when the system is ill-conditioned. Test comparisons of this algorithm and SIRT methods done by Nolet (1985) and van der Sluis & van der Vorst (1987) found that the Paige-Saunders algorithm converged more rapidly for the examples tested. Unlike SVD analysis, this method does not provide resolution and variance estimates for the resulting solution. D Chapter2. 2 Travel Time &; Amplitude Tomography 2.4 32 Iterative Inversion Scheme The complete inversion scheme is depicted in Figure 2.3. To begin, an initial model estimate v must be supplied (where v represents v , u or zo). This might be ar0 0 0 0 D rived at by applying some conventional 1 analysis technique (e.g. plane layer analysis (Ewing, 1963), slant stacking and wavefield migration (Clayton & McMechan, 1981), extremal bounds inversion (Bessonova et al., 1974)) o r simple forward modelling of a subset of the data. Ray tracing is done for the current model to determine the travel time residuals 8t and the Jacobian matrix A . If the data misfit is unacceptable, the resulting linear system is solved using LSQR with an appropriate damping factor \L. In this procedure, during a particular iteration either a velocity perturbation or an interface depth perturbation is determined. Obviously these could be solved for simultaneously, but we have chosen to keep them separate to allow the user to control the trade-off between velocity versus interface depth geometry in poorly constrained regions of the model. The model is updated and the procedure is repeated until a satisfactory data misfit is achieved. A common misfit measure is M (2.50) J = I If it is assumed that the data uncertainties are uncorrelated and that each is due to a 2 2 zero mean, Gaussian process with standard deviation crj, then X is % -distributed and 2 the expected value of % is M, the number of data. The misfit parameter will be referred 2 to as the chi-squared misfit or simply % . At the kth iteration, having obtained a solution Xfc we can obtain a linear approximation of the misfit using the difference Ax*. — b. This 2 approximate misfit value is referred to as % . p Once the velocity structure of the medium has been determined, amplitudes for the velocity model can be calculated by ray tracing and used to construct the linear system (2.25). This system is then solved directly for the effective attenuation a(x). D Chapter 2. 2 Travel Time & Amplitude Tomography Starting model v o Trace rays through current model: (5t=A(5v f 2 X ~M ? Yes Amplitude inversion. No Solve linear system for <5v k JL Update model: F i g u r e 2.3. T h e complete inversion scheme. Chapter 2. 2 Travel Time &c Amplitude Tomography 2.4.1 34 Grid Cell Size The choice of the grid cell dimensions (or equivalently the number of model parameters) is based on several considerations. Ideally the objective is to resolve as much structural detail in the model as possible and thus a finely divided grid is preferable. Likewise, a fine grid is required to permit a realistic representation of the model allowing accurate solution of the forward problem. Realistically the data contain only a finite degree of constraint on the model and thus an increase in the number of parameters is only useful up to a point, past which the solution becomes disjointed. A coarse grid relative to the ray spacing prevents this and maintains some smoothness in the solution. Also, a coarse grid accommodates more efficient calculation of the forward problem and results in a smaller linear system to be solved in each iteration of the inversion. The choice of an appropriate grid will be a compromise among these factors. Typically a grid is chosen such that the number of parameters exceeds the number of data resulting in an underdetermined system. The effect of the grid cell size will be investigated in examples that follow. 2.4.2 Damping Factor In each iteration of the velocity inversion an appropriate damping factor must be chosen. 2 As the ultimate goal of the inversion procedure is to find a model such that x ~ M, it might seem reasonable to choose the damping factor which reduces the linear approximation to the misfit so x to the final solution, xl 2 p ~ m a P But, unless the kth iteration solution is linearly close n 2 y ° t be a good estimate of the true misfit % . Consta- ble, Parker & Constable (1987), applying a similar inversion scheme to E M sounding data, suggest a method in which solutions for several values of \i are calculated at each iteration. The true misfit for each of these models is calculated by solving the forward problem in each case, and the solution corresponding to the damping factor which reduces 2 X maximally is chosen. This method removes the subjective choice of /x, but multiple Chapter 2. 2 Travel Time Si Amplitude Tomography 35 solution of the forward problem at each iteration is extremely time consuming. If the required computing capabilities are available this method is very attractive. Our approach is more subjective. In each iteration velocity solutions and misfit 2 estimates (% ) are determined for several damping factors. A solution is chosen which p doesn't overfit the data, and which changes the model significantly while satisfying the user's notion of what is realistic. Typically the approximate misfit x\ underestimates p 2 2 2 the true misfit % , and as the misfit is reduced % becomes a better estimate of % . In p most cases this allows an appropriate choice of p. to be made based on the value of x%The true misfit is provided when the forward problem is solved in the next iteration, and thus a check is always provided on the chosen velocity solution. For the attenuation inversion X X% since the amplitude-attenuation relationship is 2= linear. Thus the true misfit reduction achieved for a particular choice of damping factor 2 is accurately provided by % , without having to repeat the forward problem solution. p 2.5 Resolution Resolution of velocity or attenuation structure in a tomographic inversion relies on independent sampling of each region by various rays (McMechan, 1983). This requires not only dense ray coverage over the model, but also sampling of each region by rays at different angles. The result of limited angular sampling is smearing of the solution in the predominant direction of ray alignment in that zone. The resolution in a particular refraction experiment is determined by the number and geometry of ray paths, which are controlled by the source-receiver geometry employed and the velocity structure of the medium. To characterize the associated resolution, the linear system (2.11) can be constructed for the source-receiver geometry of interest using a representative velocity model, and the resolution matrix R can be calculated via SVD. In so doing, it is assumed that the true velocity is linearly close to the velocity model used (i.e. so that the ray paths are close to the true ray paths). In most practical situations Chapter 2. 2 Travel Time & Amplitude Tomography 36 (2.11) is very large and thus this analysis is restricted to examples with relatively few model parameters and shot-receiver pairs. Resolution characteristics of larger systems can be inferred by analogy with these examples. In the examples that follow, the resolution analysis is performed on the travel time-velocity Jacobian of (2.11), but because of the similarity of the Jacobian expressions for the velocity and attenuation formulations (compare (2.10) and (2.26)), the conclusions are valid for the attenuation problem as well. Figure 2.4 shows the ray diagram for a velocity model in which velocity increases _1 linearly with depth at 0.5 s , for a deployment of 6 sources (receivers) and 20 receivers (sources). This is representative of most refraction surveys where there is typically a high ratio of number of receivers/number of sources or vice versa. The model consists of 126 parameters. Inspection of the ray diagram reveals that the upper and central parts of the model are well covered by rays both in density and angle and thus should be well resolved, while toward the edges of the covered region the density decreases and the sampling rays have very similar inclination. This suggests that the resolution will be poor at depth and toward the edges of the model. The ordered singular values (Figure 2.5) range over 22 orders of magnitude, and there are 15 identically zero singular values. Due to the extreme ill-conditioning of the system the data are very insensitive to model structures represented by the singular vectors associated with very small singular values. Examples of the singular vectors from which the solution model is constructed are shown in Figure 2.6. There is generally an increase in singular vector structure with depth and toward the edge of the model as the associated singular values decrease in size. The singular vectors associated with the very small singular values (10 -16 21 to 10~ ) contain structure at the edge of ray coverage in the model, and the singular vectors associated with zero values are delta-like functions centred over model grid points unsampled by the rays. In a damped solution, these very small singular vectors will be excluded. D Chapter 2. 2 Travel Time & Amplitude Tomography 0 5 10 37 15 20 Distance (km) _1 Figure 2.4. Ray diagram for velocity model with constant vertical gradient (0.5 s ) . There are 120 source-receiver pairs with 6 sources (large triangles) located at 4 km intervals with receivers every 1 km. The grid cells are 1 km x 1 km in size for a total of 126 velocity parameters. D Chapter 2. 2 Travel Time & Amplitude Tomography o in CN_| 0 38 , , , , 25 50 75 100 125 Singular Value Number Figure 2.5. Ordered singular values determined for the ray geometry of Figure 2.4. There are 15 identically zero values which are plotted on the lower axis. D Chapter 2. 2 Travel Time & Amplitude Tomography b 39 Singular Vector 13 A=0.16 Singular Vector 87 A=0.17xl0 -16 F i g u r e 2.6. Singular vectors for the ray geometry of Figure 2.4. The contour intervals are 0.05, 0.05 and 0.10 km/s for figures (a), (b) and (c), respectively. Positive contours are solid whereas dashed contours are negative. Chapter 2. 2 Travel Time &: Amplitude Tomography 40 Given accurate enough data each of the model parameters sampled by the data (i.e. excluding those with associated zero singular values) can be perfectly resolved. But, due to the ill-conditioned nature of this linear system introduction of even the smallest level of noise leads to a severe decrease in resolution. This is demonstrated in Figure 2.7 where the diagonals of the resolution matrix are contoured for two different damping 8 factors. In the first case (Figure 2.7 a) where a very small damping factor of ^,=10~ has been applied, the model is perfectly resolved to 2 km depth, but resolution deteriorates at greater depths and toward the sides of the model. When heavier damping is applied -2 (/i=10 ) the resolution decreases with depth starting at the surface of the model. Thus, the effect of noise in the data will be to reduce the resolution of the model, particularly at depth. This is further demonstrated in Figure 2.8 where the resolving kernels in the centre of the model at increasing depths are shown for fi=10~ & -2 (in a,b,c) and ^=10 (in d,e,f). It is immediately obvious that damping results in a decrease of resolution with depth, but closer inspection of the kernels for fi=10~ 2 reveals that the horizontal reso- lution is more adversely affected than the vertical resolution. This horizontal smearing results from the predominantly horizontal ray angles sampling the deeper parts of the model. Resolution in the upper part of the model is maintained due to the greater variety of sampling ray angles. To demonstrate the effect that the medium can have on resolution the velocity gra_1 dient in the medium was reduced to 0.25 s . This reduces the curvature of the rays resulting in flatter rays and shallower penetration. The resolving kernels for [i=10~ 2 are shown in Figure 2.9, and as can be seen they are even broader than the corresponding kernels for the previous model. The more severe horizontal smearing is a result of the reduced range of ray angles sampling the model. The effect of the source-receiver geometry on resolution is demonstrated in the following example. Initially 3 sources are deployed at 10 km intervals across the original model with 20 receivers at 1 km intervals. In this instance the 126 model parameters are to be Chapter 2. 2 Travel Time <fe Amplitude Tomography 41 O O 0.0 10.0 20.0 Distance (km) o 0.0 10.0 20.0 Distance (km) Figure 2.7. Resolution matrix diagonals for the ray geometry of Figure 2.4. In (a), the damping factor is ^ = 10~ , and in (b) /i = 10~ . The contour interval in both cases is 0.1. 8 2 D Chapter 2. 2 Travel Time & Amplitude Tomography F i g u r e 2.8. Resolution kernels for the ray geometry of Figure 2.4. (a), (b.) and (c) show resolving kernels for p.= 10~ at i = 10 km and depth z as indicated, (d), (e) and (f) show the resolving kernels for fi.=10~ at the same locations. The maximum value of the kernel in each case is indicated. The kernels are plotted for a 20 k m x 5 km grid. To enhance the visual presentation, the kernels have been plotted using a grid which is twice as dense as that used for the velocity parameterization. The delta-like functions in (a) and (b) represent perfect resolution of the model parameter at the associated grid node. 8 2 42 Chapter 2. 2 P Travel Time & Amplitude Tomography 43 _1 2.9. Resolution kernels for a medium with a reduced vertical gradient (0.25 s ). A damping factor of /x=10~ has been used. The positions correspond to those shown in Figure 2.8, and they are plotted on the same grid. Figure 2 Chapter 2. 2 Travel Time &: AmpEtude Tomography 44 determined from 60 ray paths. The resolution matrix diagonals for fi=10~ 2 are contoured in Figure 2.10 and the 60 nonzero singular values are shown in Figure 2.11. To improve the resolution characteristics of the system the number of sources or receivers might be increased to increase the number of data constraining the model. The resulting resolution matrix diagonals are compared to the 3 x 20 case for 6 sources and 20 receivers in Figure 2.10 (a), and for 3 sources and 39 receivers in Figure 2.10 (b). The ordered singular values for the 3 cases are compared in Figure 2.11 . Referring to Figure 2.10, the overall resolution of the model has been improved in both cases as indicated by the depression of the contours relative to those corresponding to the 3 x 20 case. But, overall resolution has clearly been increased more by doubling the number of sources (i.e. that component which there are least of) as opposed to doubling the number of receivers. Also, although the overall resolution is improved by increasing the number of source-receiver pairs, resolution at depth remains poor. Both of these results are related to the angular coverage of the model. First, whereas both of the augmented source-receiver geometries result in increased ray density over the model, the 6 x 20 geometry introduces rays which sample the model at different angles and thus they contain more independent information. The 3 x 39 geometry introduces rays which have very similar geometry to those of the 3 X 20 deployment, and thus the amount of additional independent information is limited. Secondly, resolution at depth in the model is fundamentally limited by the geometry of the ray paths in this region. Regardless of the number of source-receiver pairs deployed over the limited aperture of the experiment (0-20 km), turning rays in this region will have very similar inclinations resulting in limited, angular coverage. Angular coverage in this region can only be improved by increasing the aperture of the experiment or by including other ray paths such as reflections which sample this region at steeper inclinations. 2 The superiority of the 6 x 20 geometry (for /i=10~ ) is demonstrated in Figure 2.11 by the more gradual decline of singular values for X/XMAX > -2 10 . Further inspection of Figure 2.11 indicates that the relative drop-off of the 3 x 20 singular values relative D Chapter 2. 2 Travel Time & Amplitude Tomography 0.0 10.0 45 20.0 Distance (km) Figure 2.10. Resolution matrix diagonals for various source-receiver geometries, (a) compares the resolution diagonals for 3 sources x 20 receivers (solid contours) vei^sus 6 sources x 20 receivers (dashed contours), and (b) compares the diagonals for 3 sources x 20 receivers (solid contours) versus 3 sources x 39 receivers (dashed contours). A damping factor of /i=10 " was used in all cases, and the velocity model was that used in Figure 2.4. The contour interval in (a) and (b) is 0.1. - Chapter 2. 2 D 46 Travel Time Sz Amplitude Tomography Figure 2.11. Ordered singular values for the 3 different source-receiver geometries compared in Figure 2.10. Only the singular values for X/X x > 10~ are shown. 5 MA Chapter 2. 2 Travel Time &: Amplitude Tomography to those of the other two geometries increases below X/XMAX 47 -2 < 1 0 . This suggests that the improvement in resolution of the augmented geometries relative to the 3 X 20 geometry will further increase if noise levels in the data allow use of a smaller value of u. Conversely, the relative resolution improvements to be obtained by increasing the number of shot-receiver pairs will be limited by the noise level of the data. These simple examples provide useful illustrations of how resolution is related to the geometry of the rays sampling the medium. In more complicated models the sampling of the medium will be less uniform, and the placement of sources and receivers relative to regions of anomalous velocity structure will be important in determining the resolution characteristics. Velocity anomalies will result in ray focusing and shadow zones. Velocity layering will also have a large effect on the ray geometry. Resolution in these situations will have to be qualitatively inferred from the geometry of the ray coverage. 2.6 2.6.1 Inversion Examples Velocity Inversion To demonstrate the velocity inversion method, travel times generated for the test anomaly shown in Figure 2.12 (a) were inverted with the results shown in Figures 2.12 (c)-(f). The test model, specified on a 20 x 5 km grid with l x l km cells, consists of a circular low D velocity anomaly superimposed on a 1 background velocity field. The anomaly has a magnitude of —1.0 km/s at its centre, and the background field increases from 4.0 km/s at the surface to 6.5 km/s at 5 km depth. The same source-receiver geometry as in Figure 2.4 was used resulting in 120 travel times. Travel times were calculated using the same ray tracing algorithm used for the inversion procedure, and inversions were performed on the same grid as the test model using the background field as the starting model. There are 126 parameters for this case. The effect of the velocity anomaly on the ray paths is shown in Figure 2.12 (b), where the rays avoid the low velocity zone. The presence of such low velocity anomalies can result in regions of a model being unsampled. Iteration 10 8 d -.74 Y=96 /u=.003 0 X =9.9xl0 ;U=.8 J Iteration 6 -.41 \ o ) —-.1 4 '/x=.05 / f Figure 2.12. (a) Test model anomaly, (b) Ray diagram for test model, (c) Globally smallest deviatoric solution (method A), (d) Smallest deviatoric solution (method A), (e) Globally smallest deviatoric solution (method B). (f) Globally smallest deviatoric solution (method A) for <r=5 ms noise added to the travel times. The model dimensions in all cases are 20 km x 5 km, and the contour interval in each case is 0.1 km/s. The number of iterations, anomaly amplitude, % and \i for each solution are indicated. The expected value of x is 120. 2 2 Chapter 2. 2 Travel Time & Amplitude Tomography 49 The globally smallest deviatoric solution and the smallest deviatoric solution obtained using method A are shown in Figures 2.12 (c) and 2.12 (d), respectively, and the globally smallest deviatoric solution for method B is shown in Figure 2.12 (e). No noise was added to the test model travel times for these examples, and a standard deviation of o"=0.1 ms was assigned to the data, based on how accurately the rays were traced. Both of the method A solutions are reasonably good images of the original velocity anomaly, although the global solution is more symmetric and comes closer to reproducing the correct central magnitude of the anomaly. This indicates that the global solution may produce better results when a good estimate of the background velocity field is available. The inversion using method B diverged after the 3rd iteration. Although the 3rd iteration solution (Figure 2.12 (e)) represents the gross features of the anomaly, this method was unable to resolve the finer details of the anomaly. Method B is based on the premise that velocities within individual cells of the model can be well approximated by a constant velocity. In this example where the anomaly has a maximum dimension of 4 cells, across which there are significant velocity variations, this restriction is violated. This, as well as smoothing introduced by the required interpolation at each iteration, results in the poor performance of method B for this example. Method B is better suited to situations in which the anomaly dimensions are large relative to the individual cell dimensions. In all examples that follow, velocity and attenuation solutions are calculated using method A. Figure 2.12 (f) shows the global solution for method A when Gaussian distributed noise with a=5 ms is added to the test model travel times. A damping factor one order of magnitude greater than that used in the previous examples is required to stabilize the solution. This results in increased horizontal smearing of the anomaly and a reduced magnitude. As expected, the horizontal smearing is more severe at depth. As will be seen in a later example, some of this smearing can be removed using an appropriate weighting scheme. Chapter 2. 2 Travel Time & Amplitude Tomography 50 2.6.2 Effect of Damping The effect of the damping factor p is demonstrated in Figure 2.13 , where the first iteration solution using method A is shown for three values of p ranging over two orders of magnitude. In this example, the receiver spacing is reduced to 0.25 km and the solution is specified on a 20 x 5 km grid with 0.5 x 0.5 km cells. Gaussian noise with <r=5 ms was added to the travel times. This results in 470 data and 451 model parameters. The solution for p,=1.0 (Figure 2.13 (a)) is overdamped. As a consequence, the minimization of the perturbation size dominates over the attempt to reduce the travel time residuals, and the solution is a broad, low magnitude anomaly spread over the region of densest ray 2 coverage. The linear estimate of the % misfit is Xa =9643, compared with the expected P value of 470. As a is decreased to 0.1 (Figure 2.13 (b)) the solution becomes more localized and the magnitude increases by a factor of about 10. In this case, X =586 2 p is closer to the expected value. The underdamped case is shown in Figure 2.13 (c) for /x=0.01. Here the data are overfit on the first iteration with x =335. Although the 2 p magnitude of the anomaly is increased, the most significant result of underdamping is the introduction of spurious small scale structure in the solution. In this example, the solution for /z=0.1 would be chosen over the other two solutions as it doesn't overfit the data, it appears realistic and it changes the model significantly. 2.6.3 Streak Removal Refraction source-receiver geometries, like those employed in other seismic methods, yield data for limited projection angles. The result of this is a smearing of the solution parallel to the most prominent ray orientations. But, most refraction geometries are also hampered by an imbalance of the number of sources to the number of receivers producing an inhomogeneous distribution of rays sampling the medium. These can produce streaking artifacts in the solution. The introduction of these artifacts when the model is overparameterized and their removal by an appropriate weighting scheme is outlined in Chapter2. 2 Travel Time & Amplitude Tomography 51 Figure 2.13. First iteration solutions demonstrating the effect of damping. Damping factors of 1.0, 0.1 and 0.01 were used for the solutions, and the contour intervals are 0.007, 0.05 and 0.15 km/s for (a), (b) and (c), respectively. The test model is that of Figure 2.12 (a) with 6 sources and receivers at 0.25 km intervals. The expected value of X is 470. 2 Chapter 2. 2 Travel Time &; Amplitude Tomography 52 Figure 2.14, where we consider the test model travel times used in the previous example. Figure 2.14 (a) shows the first iteration solution determined on a grid with l x l km cells so that there are 126 model parameters to be determined using 470 data. The anomaly is localized with a magnitude of —0.412 km/s. The first iteration solution was recalculated on a grid with 0.5 x 0.25 km cells resulting in 861 parameters. A damping factor was chosen so that the x\ w a s v the same for these two solutions. The solution on the finer grid is smeared out, particularly in the regions of the model where the ray paths from individual sources are highly colhmated. The streaks are a direct result of the nonuniform ray distribution. The damped least-squares solution tries to minimize the size of the model perturbation in a least-squares sense. Minimization of the size of the perturbation means that the model will be modified in the most efficient manner, i.e. the smallest perturbation which produces the required change to the travel times. As a consequence, when the ray coverage within the model is nonuniform, the solution tends to concentrate in those regions of the model most densely sampled by rays. Normally introduction of the perturbation in regions of dense ray coverage is desirable, but near the sources the angular distribution of rays sampling the model is very limited and thus will lead to smearing of the anomaly along the pencil of rays. The least-squares feature of the model gives preference to broad small amplitude solutions as opposed to localized large amplitude solutions. As the number of model parameters is increased resulting in less constraint on individual model parameters, the solution will tend to spread out. Thus, in the overparameterized case (Figure 2.14 (b)) the streak effects occur as the perturbation is smeared out into the regions where rays focus near the sources. The removal of these streaks can be approached in two ways. First, we might impose a smoothness constraint on the model which would discriminate against these narrow elongate features. This can be done using the parameter covariance matrix in (2.41) or by applying post-solution smoothing. Of these two possibilities the latter is preferred as use of a non-diagonal weighting matrix may require solution of a nonsparse system. a d Figure 2.14. Streaking effects, (a) 1st iteration solution on a grid of 1 km x 1 km cells, (b) 1st iteration solution on a grid of 0.5 km x 0.25 km cells, (c) Solution (b) after 5 X 3-point moving average smoothing applied, (d) 1st iteration solution on the same grid as (b), but with streak removal weighting and a 3 x 3 smoothing operator applied, (e) Final solution on a grid of 0.5 km x 0.25 km cells using streak removal weighting and 3 x 3 smoothing at each iteration, (f) Final solution on a grid of 1 km x 1 km cells with streak removal weighting applied. Each model has dimensions of 20 x 5 km, and a contour interval of 0.05 km/s is used for each plot. D Chapter 2. 2 Travel Time & AmpHtude Tomography 54 Figure 2.14 (c) shows the solution from Figure 2.14 (b) after application of a 5 x 3-point moving average operator. The dimensions of the operator refer to the horizontal and vertical dimensions, respectively. Some reduction in the streaking has been achieved, but the anomaly remains severely smeared. An alternative approach is to try to compensate for the nonuniform ray sampling directly. This can be done by weighting the model parameters so as to remove the bias introduced by the ray geometry. Since streaking in the solution occurs in regions where rays converge near sources, the bias toward introducing perturbations in these regions due to the high ray density may be counterbalanced by weighting each parameter inversely by the maximum number of rays from a single source affecting that parameter. The weighting matrix W in this case is diagonal so the matrix A W - 1 remains sparse. The iteration 1 solution was recalculated on the finer grid with this weighting scheme applied, and the result after applying a 3 X 3-point moving average operator to the solution is shown in Figure 2.14 (d). As desired, most of the streak effects have been eliminated and the magnitude of the central anomaly is only slightly reduced from the iteration 1 solution on the coarser grid. However, there is a broad low level anomaly that persists. The final solutions for the fine grid parameterization and the coarse grid parameterization are shown in Figures 2.14 (e) and 2.14 (f), respectively. Both solutions were calculated using the weighting scheme mentioned, and in the former case, a 3 x 3-point moving average operator was applied to the calculated perturbation at each iteration. The overparameterized solution suffers more from horizontal smearing, resulting in an anomaly with central value of —0.560 km/s as compared to —0.723 km/s for the overconstrained case. The weighting scheme improved the solution for the overconstrained case even though streaking effects were not apparent in the absence of weighting. The solution for the overconstrained case without weighting (not shown) was very similar to the first iteration solution shown in Figure 2.14 (a). D Chapter 2. 2 Travel Time Si Amplitude Tomography 55 2.6.4 Interface Inversion The fourth example (see Figure 2.15) demonstrates the inversion procedure for determination of interface depth and geometry. The test model and the source-receiver geometry deployed are shown in Figure 2.15 (a). The interface is sinusoidal with a period of 40/3 km and an amplitude of 0.5 km. It separates the upper medium, which has a _1 velocity of 2.5 km/s at the surface and a vertical gradient of 0.5 s , from the lower medium which has a velocity of 5.0 km/s at 1 km depth (at position x=0) and a gradient _1 of 1.0 s . Sources are located at 4 km intervals, with receivers at 0.50 km intervals extending to 6.75 km from the source. Gaussian distributed noise with <r=5 ms has been added to the calculated travel times. The interface was parameterized using 40, 0.5 km segments, and the starting model (i.e. the 0th iteration model) in this case used the correct velocities above and below the planar interface at 1.0 km depth. Figure 2.15 (b) shows the 0th to 3rd iterations as well as the true model. A 3-point smoothing operation was applied to the interface perturbation at each iteration except for the final iteration. By the 5th iteration the solution is visually indistinguishable from the true model except at the endpoints of the model where it is not as well constrained. Figure 2.15 (c) demonstrates the effect of using an incorrect estimate of the velocity in the lower medium. The velocity used was 0.25 km/s too low. As a result, a solution fitting the data was not found. The 5th iteration solution fits the data to an RMS value of 17 ms. As can be seen, the solution represents the flanking highs of the interface reasonably well, but underestimates the amplitude of the central low. Again the solution poorly represents the true interface at the edges due to lack of constraint there. As in.the previous case, a 3-point smoothing operation was applied to the interface perturbation for all but the final iteration. D Chapter 2. 2 Travel Time Sz Amplitude Tomography Distance (km) Figure 2.15. Interface inversion, (a) Ray diagram for test model used to generate travel times, (b) Iterations 0-3 for a starting model with the correct velocities above and below the interface. The 5th iteration solution overlies the true model (solid curve), except near the edges of the model, (c) Iterations 0,1,3 and 5 for a starting model using a velocity for the lower medium which is 0.25 km/s too low. For both cases (b) and (c) a 3-point smoothing operator was applied to the perturbation solution at each iteration except for the final iteration. 56 D Chapter 2. 2 Travel Time & Amplitude Tomography 2.6.5 57 Combined Velocity-Attenuation Inversion The final example (Figure 2.16) demonstrates inversion for both a velocity-attenuation anomaly and interface geometry. The test model (Figure 2.16 (b)) contains a sinusoidal interface centred at 1 km depth (the same as in example 4) separating the upper medium, _1 with a velocity of 2.5 km/s at the surface and a vertical gradient of 0.5 s , from the lower medium where the velocity outside of the anomalous region is 4.0 km/s at 1 km -1 depth, and the vertical gradient is 0.5 s . In the lower medium is a low velocity-high attenuation zone which is 4 km wide and dips at an angle of 65° from horizontal. The velocity in the anomalous zone is reduced by 1.0 km/s relative to the background velocity field at the same depth, and the attenuation in the anomalous zone is 0.020 whereas it is 0.0 elsewhere. The velocity model was specified on a 101 x 50 grid with cell dimensions -1 of 0.2 X 0.1 km. Travel times and amplitudes (incorporating attenuation for / =10.0 s ) c were generated for this model using the ray tracing algorithm for sources at 4.0 km intervals and receivers at 0.25 km intervals, for a total of 448 source-receiver pairs. The ray diagram for the test model is shown in Figure 2.16 (a). Gaussian distributed noise with cr=5 ms was added to the calculated travel times. For the inversion, the model was specified on a 41 x 20 grid with cell dimensions of 0.5 x 0.25 km. The starting model consisted of a plane interface at 1 km depth separating media with velocities as given above, excluding the anomalous zone. The true RMS misfit for this starting model was 107 ms (% =1.8 x 10 ). Initially a subset of 130 2 5 travel times were chosen corresponding to rays turning below the interface, to determine only the geometry of the interface. A total of 4 iterations were performed (with 3-point smoothing applied to the interface perturbation in each iteration) which reduced the RMS misfit for these 130 travel times from an initial value of 95.3 ms (x =4.7 x 10 ) to 2 4 12.1 ms (% =765). The resulting interface geometry is that shown by the solid line in 2 Figure 2.16 (c). The complete set of travel times was then included, and a smallest deviatoric velocity D Chapter 2. 2 Travel Time & Amplitude Tomography 58 Combined velocity and interface inversion, (a) Ray diagram for test model used to generate the travel times. Gaussian distributed noise with <7=5 ms was added to the calculated travel times, (b) Test model with a sinusoidal interface and a dipping low velocity body. The anomalous zone has a velocity reduced by 1 km/s relative to the surrounding region of the model. The velocity (km/s) and vertical gradient (s ) in each layer are indicated, (c) The reconstructed velocity anomaly is shown (contour interval of 0.2 km/s) along with the resultant, interface (solid line). The true model interface and the outline of the low velocity zone are shown (dashed lines) for comparison. Figure 2.16. _1 D Chapter 2. 2 Travel Time <fe AmpHtude Tomography 59 inversion was performed with 6 iterations, after which the RMS misfit had been reduced to 5.8 ms (x =580). 2 Streak removal weighting and smoothing were applied at each iteration. The interface was kept fixed during these iterations, as was the velocity above the interface. The final velocity solution (shown in Figure 2.16 (c)) well represents the true interface, although there is some distortion to the left of the anomaly and directly above it. The anomalous body is quite well resolved in the upper half of the model and the dip of the body is apparent, although the reconstructed velocity anomaly doesn't have the magnitude of the true anomaly over its width due to horizontal smearing. The solution suffers from severe smearing near the bottom of the model as expected. Following the velocity inversion, amplitudes (without attenuation) were calculated for a smoothed version of the final velocity model to construct the amplitude-attenuation system (2.25). The need for smoothing the model before amplitude calculation is discussed in Section 3.2.1. To produce a smoothed model, the residual velocity was removed from the background velocity and smoothed using a 3 x 5-point (1 k m x l km) moving average operator before it was added back to the background velocity. Only amplitudes for rays turning below the interface were included when constructing the linear system, and rays arriving near the critical distance were excluded as calculated amplitudes are inaccurate in this region and thus would not be used when applying the inversion to field data. The observed amplitudes and amplitudes calculated for the reconstructed velocity model (labelled Initial) are shown in Figures 2.17 and 2.18. Both curves have been smoothed using a 9-point (2 km) moving average operator. The initial misfit values for the normalized amplitudes U' are % =3262 and an RMS 2 misfit of 0.018 for 239 data, where an uncertainty of ± 2 0 % has been assigned to U' to account for uncertainties in the observed and calculated amplitudes. The expected RMS value for uncertainties of ± 2 0 % is 0.0054. The attenuation solution calculated for ^,=0.007 and smoothed using a 5 x 3-point operator is shown in Figure 2.19, and the resultant amplitudes are shown in Figures 2.17 and 2.18. The misfit values are reduced Chapter 2. 2 D Travel Time & Amplitude Tomography 60 Source 1 €> "D 3 ~cL E < 10 15 20 Source 2 3 \\ ~Q_ *»% i l l E < —' \/ \\ \ ' 1 i 5 10 15 20 15 20 Source 3 o 3 ~Q_ E < i 5 10 Distance (km) Figure 2.17. Calculated versus observed amplitudes. Amplitudes calculated for the reconstructed velocity model without attenuation (initial), and including the attenuation solution (final) are compared to the observed amplitudes. Amplitudes within the stippled region were excluded in the inversion as calculated amplitudes near the critical distance are inaccurate. Triangles mark the source positions. D Chapter 2. 2 Travel Time & Amplitude Tomography Source A 4) "NT" \ \ \ \ \ V N V E < 61 1 1 //> a I • 10 - 15 20 Source 5 \J u. E «' ' "* < * " " / / b T I 10 15 20 15 20 Source 6 o "a. E < 1^ 10 Distance (km) Calculated versus observed amplitudes. Amplitudes calculated for the reconstructed velocity model without attenuation (initial), and including the attenuation solution (final) are compared to the observed amplitudes. Amplitudes within the stippled region were excluded in the inversion as calculated amplitudes near the critical distance are inaccurate. Triangles mark the source positions. Figure 2.18. Chapter 2. 2 D Travel Time & Amplitude Tomography 62 o Distance (km) F i g u r e 2.19. Attenuation solution. The reconstructed attenuation anomaly is shown, with the true model inteface and the outline of the true attenuation anomaly indicated by the dashed lines. Chapter 2. 2 Travel Time h Amplitude Tomography 63 2 to % =210 and an RMS value of 0.0051. The attenuation above the interface was kept fixed for the inversion. As for the velocity solution, the edges of the anomalous body are reasonably well determined, although the effect of horizontal smearing is apparent. Definition of the anomaly just below the interface is not as good as in the velocity solution since rays turning just below the interface were excluded because of their proximity to the critical ray. The dip of the body is not well defined in the attenuation solution. The regions of negative attenuation in the solution occur where the amplitude reduction due to elastic effects is greater for the velocity solution than for the true model. Thus, improper compensation of the amplitudes for elastic effects produces erroneous structure in the attenuation solution. In applying the amplitude inversion to field data additional complications will arise. A stable method for estimating the observed amplitudes from the data is required, and an appropriate scaling factor must be determined in order to compare the calculated and observed amplitudes. Also, there will be instances where there is more than a single arrival contributing to the observed amplitude. All of these complications are encountered in the amplitude inversion in Chapter 4. 2.7 Summary D The examples of the previous section demonstrate that 2 velocity and attenuation structure can be determined from surface-surface refraction first arrival travel times and amplitudes using an iterative tomographic scheme. This method is particularly suited to situations where a good estimate of the background velocity field of a medium is available (as is usually the case for seismic refraction experiments) to be used as a starting model D for determination of 2 anomalous structure. As the complexity of the model increases, a better estimate of the starting model is required and the resultant model will be more reliant on the starting model. D Chapter 2. 2 Travel Time Si AmpHtude Tomography 64 The resolution achieved in such experiments can be improved by increasing the number of sources and receivers, but as is the case in all seismic tomography the resolution is fundamentally determined by the velocity structure of a medium. Apart from regions of the model which may remain unsampled due to the effect of velocity variations on ray paths, the resolution is strongly affected by the background velocity gradients and layering of a medium. When the vertical gradient is small, turning rays are predominantly horizontal leading to a decrease in horizontal resolution. Velocity layering results in regions of a model (i.e. the lower regions of each layer) which are sampled exclusively by steeply oriented rays when only first arrival times are used. This results in vertical smearing in these regions or, for interface inversion, superposition of false structure on the underlying interface. Resolution can be improved by including travel times and amplitudes for secondary arrivals (reflections and wide angle refractions), but it is often difficult to determine these with the required accuracy. In the absence of sufficient constraint in different regions of a model, bias may be introduced using the weighting matrix W. Seismic tomographic methods which utilize straight ray paths require no iteration (i.e. the effect of velocity perturbations on the ray paths is not considered) and thus the existence of a velocity solution which fits the data depends only on the consistency of the linear system of equations. There is no assurance that the solution is physically reasonable, for instance, that it doesn't introduce shadow zones. In the iterative scheme used here, there is no guarantee that the method will converge, but if it does converge, then because rays are traced through the updated model at each iteration, we are assured that the solution model is capable of producing the observed travel times. This is some small consolation for having to deal with the nonlinear problem. Chapter 3 Dynamic Ray Tracing 3.1 Introduction Construction of the linear systems of equations used for inversion (equations (2.11), (2.14), (2.16) or (2.25) of Section 2.2) requires two-point ray tracing to determine either travel times or amplitudes, as well as ray paths for the set of source-receiver pairs. The actual ray paths are required to calculate the travel time—velocity Jacobian matrix or the amplitude-attenuation matrix in method A and the partial path length matrix in method B (see Section 2.2.1). In general, a large number of source-receiver pairs are involved in tomographic reconstruction and the iterative nature of the velocity inversion scheme requires that the entire suite of rays be traced through the model each time it is updated (Section 2.4). Thus it is essential that the method of ray tracing be highly efficient. Many ray tracing schemes have been outlined in the literature. These can generally be categorized as numerical methods (e.g. Andersen & Kak, 1982; Smith, Goldberg & Liu, 1980; Julian & Gubbms, 1977; Pereyra et al, 1980; Chander, 1975; Julian, 1970; Cerveny et al, 1977; Lytle & Dines, 1980 ) or analytic methods (e.g. Will, 1976; Whitt & Clowes, 1979; Marks, 1980; Cassell, 1982; Langan et al, 1985). Numerical metho tend to be computationally more intensive than analytic methods which utilize analytic formulae locally valid within subregions of the model. The two-point ray tracing problem can be formulated as an initial value problem or as a boundary value problem. Ray tracing methods based on these two approaches are referred to as shooting methods or bending methods, respectively, and have been compared by Julian & Gubbins (1977). They found 65 Chapter 3. Dynamic Ray Tracing 66 the bending method was more efficient than the shooting method when determining a single ray between fixed endpoints. We have chosen to use a shooting method which utilizes analytic expressions within subcells of the velocity model. For each source (receiver) we require many ray paths from a line of closely spaced receivers (sources). In this situation a shooting method becomes more efficient as each ray traced provides useful information for determining neighbouring rays. This method also avoids the need for providing an initial estimate of the desired ray path. The velocity-attenuation model is defined by specifying velocity and attenuation values at the nodes of a rectangular grid. For the purpose of ray tracing the rectangular cells of the grid can be divided into triangular (tetrahedral) subcells within which the velocity and attenuation gradients are constant and are determined uniquely from their values at the triangle (tetrahedron) vertices. This is a convenient parameterization as the velocity field so defined is continuous across the subcell boundaries in contrast to use of rectangular cells where artificial first-order discontinuities are introduced. A parameterization utilizing triangular cells was used by Marks (1980) and Chapman & Drummond (1982). Specification of the model using constant velocity gradient cells allows use of analytic formulae for both kinematic (i.e. ray paths and travel times) and dynamic (i.e. amplitudes) ray properties as described by asymptotic ray theory. Also, approximate attenuation as well as travel time-velocity and amplitude-attenuation inversion formulae can be expressed analytically. Section 3.2 briefly reviews asymptotic ray theory and its restrictions, followed by kinematic and dynamic formulae for ray tracing (Sections 3.3 and 3.4). The calculation of attenuation is described in Section 3.5 and inversion formulae for both travel time -velocity inversion and amplitude-attenuation inversion are given in Section 3.6. Chapter 3. Dynamic Ray Tracing 3.2 67 Asymptotic Ray Theory Asymptotic ray theory (ART) describes seismic wave propagation using a high frequency approximate solution to the elastodynamic equation for an isotropic medium (Cerveny & Ravindra, 1971). This solution is convenient as it allows separation of compressional and shear wave propagation in a heterogeneous medium. It is assumed that a harmonic solution for the particle displacement W of the form oo j W(cu,x) = exp{iu;(i-r(x))}^(^)- W (x) j (3.1) exists and is asymptotic to the exact solution as u —> oo. Wavefronts are defined as surfaces of constant phase (r(x) = t) where t is the travel time, and rays are defined as trajectories orthogonal to the wavefronts. The amplitude coefficients Wj(x) modulate the wave amplitude in space. In most seismic applications the zero-order term of the expansion provides a good approximation to the true solution. Substitution of (3.1) into the elastodynamic equation results in the eikonal equation 2 2 (Vr(x)) = v~ (x) (3.2a) and the zero-order transport equation VT-VU + ^ ( V r + V r • V(ln p)) = 0 2 (3.26) for compressional waves, where v(x) and p(x) are the compressional velocity and density of the medium, and the amplitude U in (3.2b) represents the particle displacement or some time derivative of the displacement (e.g. particle velocity). The operator V r - V represents a directional derivative along the ray. The kinematic properties of the wavefronts (or rays) are determined by (3.2a), whereas the dynamic properties (i.e. amplitudes U) are governed by (3.2b) in conjunction with (3.2a). For compressional waves, the direction of Chapter 3. Dynamic Ray Tracing 68 particle motion is parallel to the ray in the zero-order approximation. As the zero-order solution corresponds to a solution obtained using the principles of geometrical optics, the term geometrical will preface references to rays governed by (3.2a) and (3.2b). 3.2.1 Restrictions of ART That the ansatz (3.1) is an asymptotic expansion in a ; -1 implies that the error in re- taining only the zero-order term goes to zero as the frequency becomes infinite. finite frequencies, the zero-order solution will be in error. requires that the frequency f » At Minimization of this error | V u | and when interfaces are present the radius of curvature of the interface must be much greater than that of the wavefront. According to Cerveny & Ravindra, (1971, p.12) the expansion (3.1) is valid where the phase function T ( X ) is analytic. Regions where T ( X ) is non-analytic include shadow boundary rays, caustics and critical rays or grazing rays which correspond to end points, cusps or tangent branches of the travel time-distance curve. The size of the neighbourhood about these _1 points where A R T breaks down is proportional to u> . Except in the neighbourhood of these singular points, zero-order A R T provides a good approximation for direct, reflected and turning rays. First-order A R T also accurately describes the pure head wave away from the critical point. Other situations where A R T is inadequate include interference head waves, Fresnel shadows, edge and point diffractions, interface diffractions and gradient coupling (Chapman, 1985). At distances away from the critical distance, turning rays provide a good estimate of the interference head wave. The limitations mentioned above apply to A R T in general, irrespective of the numerical algorithm used to implement it. However, there are additional restrictions of A R T which are a result of the model parameterization that we have chosen to use. Although a continuous velocity field is defined by a model of triangular cells, the velocity gradient is generally discontinuous across cell boundaries. The effect of these discontinuities on the amplitudes of finite frequency waves may be insignificant, but because A R T is a high Chapter 3. Dynamic Ray Tracing 69 frequency approximation the calculated amplitudes are affected by gradient discontinuities. The resultant effect for 2 models consisting of many cells is to produce scatter D in the amplitude versus distance curves. Our approach to remedying this problem is an ad hoc one. For amplitude calculations the velocity model is arbitrarily smoothed in an attempt to minimize the size of gradient changes from cell to cell, and the resulting amplitude-distance curves may be smoothed as well. An alternative solution would be to modify the algorithm to implement either Maslov theory (Chapman & Drummond, 1982) or Gaussian beam theory (Cerveny, 1985) which include implicit smoothing. The analytic formulae used in Section 3.4 are readily incorporated into a Maslov-based algorithm (Chapman, 1985). 3.3 Ray Kinematics The eikonal equation (3.2a) can be solved using the method of characteristics to obtain the ray equations , Tr= V Vv du "' T r = ~ ~ 3 3 <'> x where u = V r is the wave slowness (ju| = v~ ) and the ray path is expressed parametrically as x = x(r). Alternatively, the ray equations can be written as i, ( ; t ) = v (;) <-> 3 4 where dl represents an increment in path length along the ray. Solution of these equations for a constant velocity gradient medium results in simple analytic expressions for the ray paths and travel times which can be used in the constant velocity gradient cells of the model. In later sections, the following notation conventions are adopted. Cartesian coordinates will be specified either as x,y,z or x x , Xz, according to which is more convenient. ly 2 Similarly, components of the slowness vector u in the x-z plane will be synonymously Chapter 3. Dynamic Ray Tracing 70 referred to as p and q, or u and u$. Vector components specified with respect to a local y coordinate system will generally be denoted by primes, whereas vector components specified with respect to the home coordinate system will be unmarked. Quantities defined at the start point of a ray (either the source, or the start point within a particular model cell) will be indicated by a (*), whereas quantities for subsequent points along the ray will be unmarked. 3.3.1 Analytic Ray Tracing Formulae Consider a ray emanating from x* in a medium with a constant velocity gradient Vv (see Figure 3.1). The home coordinate system is defined by the unit basis vectors e ly e and 2 e which form the x,y and z axes. The direction of the ray at any point x is denoted by 3 the unit vector n, with n* = n(x"). We can define a new coordinate sj'stem (x'-system) such that the z'-axis is aligned with Vv and the y'-axis lies perpendicular to the plane defined by n* and X?v. The a:'-axis is chosen to form a right-handed cartesian coordinate system. The unit basis vectors defining the x'-system are ~/ e l -i — e ±n*xe |rrxe | - / 2 X e 3 3> 3 (3.5) X7 V In 2 D we choose the sign in the expression for e such that e • e = 1. The ray normal 2 2 2 rh is a unit vector perpendicular to ii and rh* = rh(x*). The vector rh is defined by rh = nxe . 2 (3.6) These two systems are related by the orthogonal transformation x' = Tx, x = T-V, (3.7) Chapter 3. Dynamic Ray Tracing 71 F i g u r e 3.1. Circular ray geometry in region of constant velocity gradient. Symbols are denned in the text. Chapter 3. Dynamic where T _ 1 Ray Tracing 72 T = T , and Tn = K • S;. (3.8) In the x'-system the velocity is given by v(z') — v -f k(z' — z' ) where v is the velocity 0 Q 0 of the medium at z' and k = \Vv |. The wave slowness is 0 u' = ( p, 0, q), (3.9a) and the, local horizontal and vertical slowness are defined by Vv • rh P = L - Vt; • n 9 ' (3.96) 1 ' = It is well known that ray paths in a constant velocity gradient medium are circular arcs in the plane (x — x") • e' = 0. Telford et al. (1976) derive parametric expressions for the 2 ray path and travel time in terms of the ray direction. They can be expressed as (x' - x ) = - ^ ( V v • (rT - n), w anc i = 1 , f /v*\ (k InH — U k 0, V v {m - m")), + Vv • n")) (3.10) , 3.11 \ \ v J (k + Vv • n) J v ; The ray path formula (3.10) and conservation of the local horizontal slowness p along the ray segment are obtained by integrating (3.4). The radius r and centre of curvature x^. c of a ray are 1 \Vv • rh*| « - x*') = - ^ ( V T ; • n, \pk\ 0, - V v • rh*). (3.12a) (3.126) Chapter 3. Dynamic Ray Tracing 73 Because a ray follows a circular arc, the path length along the ray is simply . = (3.13) where A# is defined to increase with time along the ray. Applying T - 1 to (3.10) we obtain — - {Vv[(n • Vv)Vv • (n - n") + (m* • Vv)Vv • (rh - riV)] ¥ (x - x") = \Vv\*(Vv z • m') 2 -n*{\Vv\ Vv • (n - ii*)]} , (3.14a) with the additional restriction (x-x*)-e' = 0 (3.146) 2 which confines the ray to the x-z plane. 3.3.2 Cell Boundary Intersections D C In either 2 or 3 the exit point of a ray from a cell is determined by the intersection of the ray with a straight line A'x' + C'z' + D' = 0. D (3.15) D In 2 where the cell boundaries are straight line segments, this is obvious. In 3 the cell boundaries are planes which can be represented as A'x' + B'y' + C'z' + D' = 0. Because the ray is restricted to the sagittal plane y' — 0, it will intersect the cell boundary along the line formed by the intersection of the sagittal plane and the boundary plane. Using (3.15) and (3.10) we obtain a quadratic in either x' or z' resulting in two solutions. That solution corresponding to the shortest path length in the direction of propagation is chosen. The ray direction vector n', partial traveltime and partial path length at the intersection point are determined using (3.10), (3.11) and (3.13), respectively. The Chapter 3. Dynamic Ray Tracing 74 direction vector and intersection point relative to the x-system are given by applying the rotation T - 1 . 3.3.3 Interface Intersections and Snell's Law A ray is refracted/reflected at an interface according to Snell's Law = -!—, V (n;xt)-n = r 0, (3.16) Vi r where the subscripts i and r refer respectively to the incident and refracted/reflected ray, and t is the unit interface tangent at the intersection point. The second condition requires that the refracted/reflected ray lie in the plane formed by the incident ray direction vector (at. the interface) and the interface tangent. If | n • t| > 1 in (3.16) for the refracted ray, r then we obtain only the reflected ray. D Interfaces within the 2 model are represented by cubic splines as 2 3 z = a + bx + cx + dx . (3.17) The intersection of the ray with the interface is a solution of the 6£/i-order polynomial in x formed from (3.14a) and (3.17). The appropriate root is real, it corresponds to coordinates located within the cell and the associated path length is positive. If the path length to the intersection point is less than the path length to the nearest cell boundary intersection, then the ray intersects the interface within the cell. The unit tangent and normal to the interface at x are t = (1 + E')->fy 2 where E = b + 2cx + 3dx . , 1 = (1 + E )-> ( " f ) , 2 (3-18) Chapter 3. Dynamic Ray Tracing 3.3.4 75 Two-point Ray Tracing The 2 D source-receiver geometry of interest is one in which shots from a sequence of sources located along a subhorizontal interface are recorded by a series of recorders on another subhorizontal interface. The problem as formulated earlier requires that we perform ray tracing between specific source-receiver pairs. This is done using a shooting method. For each of the sources the range versus take-off angle function X(9) is defined. r The appropriate take-off angle 9 for a receiver at distance A can then be determined from the inverse function 9(X). In general X(9) will be multivalued over the domain of 9, and the different branches of the function may be discontinuous. As a result, 9(X) will also be multivalued, discontinuous and will contain gaps (i.e. shadow zones). A simple example of X(9) corresponding to the occurrence of a triplication in the travel time curve, is shown in Figure 3.2. To define X(6), rays are traced from the source for a range of angles 6MIN < 9 < 9MAX using a systematic method to ensure coverage over the distance range XMAX • XMIN < X(9) The resultant set of ranges obtained is then divided into groups in which < X(9) is a monotonic function for 9^ < 9 < 9f and X? < X{{9) < X*. The set of inverse functions 6i(X), 9* < 9i < 9f, X * < A" < A t r + ; are then single valued. Then for each receiver, the range of each of the inverse functions is checked for containment of the corresponding distance, and the appropriate take-off angle 9 for each group is obtained from 9i(X). Each angle corresponds to a different arrival at the receiver, from which the appropriate arrival must be determined. 3.4 3.4.1 R a y Dynamics Connection Formula for Geometrical Rays For compressional seismic waves propagating in an isotropic heterogeneous medium, a connection formula relating amplitudes at different points along a ray can be obtained Figure 3.2. Distance versus take-off angle curve for the occurrence of a triplication. For two-point ray tracing, the multi-valued inverse function 0(X) is defined in terms of several single-valued functions O^X), X* < X < Xf. Chapter 3. Dynamic Ray Tracing 77 from the transport equation (3.2b) (Cerveny & Ravindra, 1971, p.74). If U(x) and U(x ) 0 are the amplitudes at positions x and x , then 0 Zj ff(x) = G(x,x ) (3.19a) 0 where G(x, x ) is the geometrical spreading function defined by, 0 1/2 (3.196) and the various quantities are defined as ./V = number of velocity discontinuity interfaces along ray, X j = position where ray intersects j-th interface, Zj — Zoeppritz coefficient at the j-th interface, dtr(x) = cross-sectional area of the ray tube at x, p(x) = density of the medium at X , v(x) = velocity of the medium at x, /c(x,x ) 0 = K M A H index. The quantities marked with tildes are evaluated on the side of the interface where the ray emerges, while unmarked quantities are evaluated on the side of the interface where the ray is incident. The K M A H index /t(x,Xn) (Chapman & Drurnmond, 1982), which is initially set to zero, is incremented each time the ray passes through a caustic producing an appropriate phase retardation in (3.19a). The amount by which the KMAH index is incremented is equal to the number of ray tube dimensions that vanish at the caustic (i.e. either 1 or 2). Zoeppritz coefficients for the refracted or reflected wave amplitudes at an interface are calculated using a program described by Young & Braile (1976) which Chapter 3. Dynamic Ray Tracing 78 is based on expressions from Cerveny & Ravindra, (1971, p.63) for plane waves incident at a plane boundary. This same routine also allows calculation of surface conversion coefficients and free-surface reflection coefficients. It is assumed that Poisson's ratio for the medium is 0.25 except for water. The densities required in (3.19a) are estimated from the compressional wave velocities v using the empirical formula p = 1.733 + 0.1644?; (3.20) obtained by a straight line fit to the velocity and density values compiled by Ludwig, Nafe & Drake (1970, p.74). The geometrical spreading function G(x, x ) is the product of two factors. The second 0 factor (G ) accounts for the effect of interfaces on the ray tube which can be expressed 2 as 1/2 N n 1/2 m, (3.21) m„ (Cerveny & Ravindra, 1971, p.79). The first factor (Gi) represents the spreading of the ray tube between interfaces and can be determined as follows. The ray path can be specified as x = x(7i,7 ,r) using the curvilinear ray coordinates 71, 72 and r, where 2 the first two coordinates identify a particular ray and r is the travel time along the ray. For r = constant the coordinates 72 and 72 define the wavefront. Then in general, an element of area on the wavefront at x is by definition aV(x) = <9x dx 071072 (3.22) which simplifies to aV(x) = dx_ dx dj2 572 (3.23) Chapter 3. Dynamic Ray Tracing if the ray coordinates are orthogonal. 79 The quantities in (3.23) are governed by the following system 1 d (duj\ _ 2 * d Inv Idx \ k which is obtained by differentiating the ray equations (3.3) with respect to -jj. For a 2 D velocity medium (i.e. v = v(x,z)), we choose the ray coordinate axes 71 and 72 such that 71 lies in the x-z plane and 72 is normal to the x-z plane. Then (3.23) can be written as, dx dcr(x) = dy 071 d~ti 1 (3.25) d~]2- = dcr^da 1 where d*W and der - are the parallel and normal cross-sections of the ray tube with respect to the x-z plane. The connection formula (3.19a) requires knowledge of a reference amplitude (7(x ) 0 and ray tube cross-section dcr(x ) in order to determine the amplitude at some point 0 x further along the ray. To provide an arbitrary reference, velocity variations near the source are ignored so that the wavefront arbitrarily close to the point source is a sphere of radius r. For simplicity it is assumed that the source is omni-directional and the amplitude over the wavefront (at r) is set to 1. The ray tube cross-section components on the wavefront for a ray in the x-z plane are d<T 0 = dx dix where d-y-L rdp Z~ 1 9i do-Q dy 572 0 d-y = rpv\d^ 2 2 (3.26a, b) and q\ are the velocity and vertical slowness at the source. Because we are considering relative amplitudes, the factor r in (3.26) can be set to 1, but it will be retained to preserve proper dimensions in equations to follow. Chapter 3. Dynamic Ray Tracing 80 Out-of-plane Geometrical Spreading Equation (3.24) for the out-of-plane spreading factor d ( dy \ ~,du reduces to d Idu \ 2 2 for a 2 velocity medium. Integrating (3.27b) and (3.27a) we obtain D dj Jv dr. (3.28) =p. (3.29) 2 2 The constant term is du \ 2 2 4 V<W 0 Taking the ratio of (3.28) and (3.26b), we obtain do~ (x) 1 r ,, ,; \ = / u dr = L 2 J air^^xoj r u " Jo 1 rv^ r „ / vdl. , 3.30 JL An analytic form for this integral is easily obtained for our model parameterization. Within a given model cell (a;' — x' )/p for a circular ray path, 0 J vdl = < v (z' - z' ) + k(z' - z' ) /2 for rh • Vv = 0 and 2 0 . vlj, if Q 0 ^ 0, (3.31) Vt> = 0. The quantities in (3.31) are defined in Section 3.3.1. Thus, if the ray is a circular arc within each cell, where is the velocity at the starting point of the ray, p and x are the horizontal s 3 slowness and range in the sth cell, both specified in the local coordinate system. A similar result was obtained by Marks (1980) and was implemented by Spence et al. (1984). Chapter 3. Dynamic Ray Tracing 81 In-plane Geometrical Spreading Previous implementations of ART for amplitude determination (e.g. Marks, 1980; Spence et al, 1984) have used analytic expressions to calculate out-of-plane geometrical spreading, but numerical methods have generally been used for calculation of in-plane spreading. An exception to this is Cassell (1982), who used analytic expressions for constant velocity blocks. An analytic expression for the in-plane ray tube cross-section can be obtained (after Chapman 1985) by considering the propagation of a neighbouring ray within and across the boundaries of constant velocity gradient cells (see Figure 3.3). In the development that follows, the subscript s enumerates model cells and boundaries along a ray rather that specific cells. The in-plane cross-section of the ray tube is da — derm, where for the present da^ is denoted simply by da. Within the sth cell of the model a ray is identified by the local ray parameter (or horizontal slowness) p , and a neighbouring ray 3 has ray parameter p„ -+• dp . The horizontal projection of the ray tube cross-section at s x" is da~(m* • e' ) and at x it is da (m • e' ). The change in the horizontal projection le s e ls of the cross-section across the cell will be due to the change in orientation of the central ray and the divergence of the neighbouring ray. The latter contribution is associated with a change in the ra}' tube cross-section and can be expressed as | ^ d p . This term 5 is obtained by retaining onfy first order terms in an expansion of the range x about the s central ray parameter p . Removing the effects of the ray orientation on the cross-section s projection by dividing by the projection of the ray normal we obtain ox da,(rh • e ^ ) B 1 1 = oV;(rh; • e ^ ) " + -^d Pe (3.33) op s relating da and da*. The term | £ i is obtained by differentiating (3.10), s dx, dp* 1 Ps x s (3.34) (q>*)(q v )' 3 s Chapter 3. Dynamic Ray Tracing 82 Figure 3.3. Neighbouring rays propagating across the sth cell and refracted at the sth boundary. Symbols are defined in the text. Chapter 3. Dynamic Ray Tracing 83 T Denning a — (da,dp) , then for transmission across the sth cell a = Ta; (3.35) s where T T = u '-, 12 <z,X T = 0, T = l. 21 3.36 22 Vsq v e t In general, the local ray parameter p and the cross-section da will be discontinuous across cell boundaries, and thus to determine these quantities along the ray, we require connection formulae for dp and da. The rays are connected across the boundary by Snell's law (3.16). A connection formula for the cross-sections is obtained from the geometry of the ray tube at the boundary. From Figure 3.3, da. ~ da*,, dt. = -^-4r-t. = — 2 ± ^ t , m, • t, m; • t„ (3.37) +1 This relation in conjunction with Snell's law expressed by (3.16) allows determination of da* from da„ across the boundary. +l A connection formula for dp across the boundary can be obtained by perturbing both the ray parameter equation (3.9b) and Snell's law (3.16). The algebra required is straight forward but tedious, and thus we proceed directly to the result: 8+1 1 <+A +i (t. • m; )(t, • rh ) 2 +1 2 kv s s q 3 v s (t, • rh; )(t • m ) J s +1 s s m + '^ '^ qv s qv s }' • - dp,. t • rh* s +1 (3.38) The connection formula for da and dp can be written as &s l = R S Q , < T , + fl s (3.39) Chapter 3. Dynamic Ray Tracing 84 where the 2 x 2 matrices R„, S„ and Q are, s K n - 1, K 1 2 - 0, K 21 - - 2 ^s+i^s+i 5n = VM;+1, S = S = 0, 12 . It, • m , g + 1 ^ 2 - j , f t„ • m S = 1; 21 8 + 1 (3.40) 22 (t, • m ) ' <7 u. s a Note that Q is calculated in the sth cell, S at the sth boundary, and R„ in the (s -\-l)th s a cell. Combining this with the expression for transmission across the cell, propagation through a sequence of model cells is given by * = T R _ S _ Q _ ...R S Q T ov n n 1 n 1 n 1 1 1 1 (3.41) a T For a point source da^ = 0 at the origin, and thus a\ = (dai,dpi) = (0, \) dp. Referring the cross-section da to the spherical wavefront at r (see (3.26)) where da = rdp/q^, 0 ( £ ) H is given by (3.41) for <r\ = (0,^/r). 1 Caustics occur where one dimension of the ray tube vanishes (i.e. either da - or da^). D In a 2 velocity medium where v = v(x, z) only the in-plane dimension of the ray tube can vanish, and thus da^ can be monitored within each model cell for sign changes indicating whether the ray has passed through a caustic. The KMAH index can then be updated accordingly. 3.4.2 Amplitudes for Reversed Rays Often, for purposes of efficiency, rays are traced in a direction opposite to the physical direction of ray propagation (i.e. from receiver-source instead of source-receiver). A l though the zero-order ART kinematic ray properties are unaffected by the direction of propagation, the dynamic properties depend on the direction of propagation. To determine the amplitude for a ray traced in the opposite direction, the start and endpoint Chapter 3. Dynamic Ray Tracing 85 quantities as well as the incident and emergent quantities in the amplitude expression for U(x) (i.e. in (3.19a)) must be switched. The incident and emergent quantities used to calculate the Zoeppritz coefficients at each interface must also be switched, and the surface conversion coefficient (if appropriate) must be determined at the start point of the ray. The geometrical spreading factors along a ray traced in opposite directions between points x and x are related by v(x )G(x,x ) = v(x)G(x ,x) where G(x, x ) is the 0 0 0 0 0 spreading factor going from x to x and G(x ,x) is the spreading factor for the opposite 0 0 direction (Richards, 1971). Head Waves 3.4.3 Head wave amplitudes are also easily calculated for the model parameterization that we employ. The compressional head wave amplitude for a stack of homogeneous layers with dipping interfaces according to first-order ART is n-l v T (n-t ) H kH kH W3/2(^ /^) /2 H 1 G l G 2 ( . n l f c H ) n (3.42a) (3.426) (3.42c) where ku =number of head wave boundary, n =total number of ray segments tfe ,ljt =unit tangent and normal to interface kjj at ray incident point, H H nii,rh = incident and refracted unit ray normal at interface, r Vj =constant velocity between boundaries j — 1 and j Chapter 3. Dynamic Ray Tracing 86 Ij =length of the jth ray segment, Vjj =velocity below the head wave boundary, I =length of the ray path along head wave boundary, UJ = dominant angular frequency of the head wave, C Vk =velocity above the head wave interface, H Zj — Zoeppritz refraction/reflection coefficients. These are equations (5.22) and (5.29) of Cerveny & Ravindra, (1971). The head wave coefficient l \ is H (3.43) 2 2 where P = (1 -v Q ) 3 H 1/2 1 and 0 = (n-t*.*)^ {Cerveny & Ravindra, 1971, pp.108-109). Zf and Z31 are the Zoeppritz coefficients at the head wave boundary for a ray incident 3 at the critical angle and for the emergent ray from the head wave boundary, respectively. Spence et al., (1984) implemented these formulae by reparameterizing each constant velocity gradient model block into a series of thin homogeneous layers. This can be avoided by considering a constant velocity gradient as the limiting case of a stack of homogeneous layers as the layers become infinitely thin. Then, within the sth constant velocity gradient cell, the product in G\ reduces to (3.44) and the sum in G reduces to / vdl which is given by (3.31). 2 3.5 Attenuation The attenuated amplitude of a wave produced by a narrow band wave source with dominant frequency u> can be approximated by c Chapter 3. Dynamic Ray Tracing t /(x) = C / £ 87 ( x ) e x p j - | / t (3.45) ^ * } It is assumed that a(x) and v(x) are frequency independent, and the attenuation calculated at UJ is applied to the signal at all frequencies about u> . The calculation of C c the amplitude /7E(X) for propagation in an elastic medium has already been described in Section 3.4. Evaluation of the attenuation factor in (3.45) requires calculation of the integral term. Like the velocity, attenuation a(x) is defined so that V a is constant within each model cell. In a medium where a(x) = a + V a • (x — x ) we have 0 0 (3.46) Within each triangular (tetrahedral) cell of the model the attenuation gradient is determined using the values at the vertices of the triangle (tetrahedron). Thus to evaluate (3.46) it remains to evaluate each of the integral terms. In order they are designated as I , Ix, Iy and I - Defining 0 z (3.47) so that the integrals are expressed in the local x'-system, the required quantities J , I x y and J are obtained by 2 = T- l 0 Evaluating the integrals in (3.47), for a circular ray path h =t, — k - — k (3.48) Chapter 3. Dynamic Ray Tracing 88 If the initial ray direction n* is parallel to Vv and Vv ^ 0 then I' = [x*' — aj')io, or if x Vv = 0, then the integral expressions in the x-system are given directly by (3.50o,'fe) Thus the amplitude including attenuation can be written as (3.51) where the integral is evaluated in each of the n cells using the above expressions. 3.6 Inversion Formulae 3.6.1 Travel time-Velocity Jacobian Formulae To form the linear systems (2.11) or (2.14) we require the travel time-velocity Jacobian matrix or partial path length matrix. In the latter case the path length of a ray in each cell can be obtained directly from (3.13). In the former case, Bregman (1986) has derived analytic expressions for the Jacobian term (2.10) for a constant velocity gradient medium. In a medium where v(x) — v + Vi> • (x — x ) we have 0 N ( dv . dV v (X dvj 8VJ 3 =1 0 x 0 - X) 0 + 8V v v (y - yo) + (z - Z ) 0 SVj. (3.52) Comparing this with (2.7) and (2.8), (2.10) can be written as (3.53) Chapter 3. Dynamic Ray Tracing 89 Within each triangular (tetrahedral) cell of the model the velocity gradient is determined using the values at the vertices of the triangle (tetrahedron). Thus each of the factors preceding the integral expressions in (3.53) are constant for a particular triangle (tetrahedron), and they will generally be nonzero for only 3 (4) values of Vj. Thus to evaluate (3.53) it remains to evaluate each of the integral terms. Designating the integrals as I , 0 I , I and I , and defining x y z the required quantities I , I and I are obtained by (3.48). In this manner, Bregman x y z (1986) derived the following results. For the case of a circular ray path Jo = j { ? " - » } , t'^-j'v (3.55a, 6) +l The above expressions are simplified somewhat from those given by Bregman (1986). If the initial ray direction is such that rh • Vv = 0 but \Vv\ ^ 0 then the above expressions are replaced by i / „ \ i 1i _ 1i V I' = (»"' - x' )I , I' = ~-^I x 0 0 Z (3.56a,6) V* (3.56c, d) 0 and if the ray travels in a constant velocity medium then the integral expressions in the x-system are given directly by t = — J(x - z*) + (z - z*) , 2 2 v 0 ^—^ J = — 0 (3.57a, fe) v 0 + [x - x ) 0 I, 0 I = Z—P- + (z* - z ) z 0 J . 0 (3.57c, d) Chapter 3. Dynamic Ray Tracing 90 Using these expressions we can calculate the contribution to from each cell that the ray encounters. 3.6.2 Amplitude-Attenuation Inversion Formulae To form the linear system (2.25) we must be able to calculate (2.26). Again we consider the medium in which a(x) = a + V a • (x — x ) . This can be written as 0 ^ /doi 0 dV a, 0 . x -w = g (sj + -g^(' <9V„a „ - -) + . ^(y dV,a • .\ - y°) + , (^8) - * ) ) «>• Comparing this to (2.24), (2.26) can be written as A _ Ooo 0Vaa rj (x[x da f j dld£ 8V - -x ) daj JL, v daj v dVya j (y-Vo) daj ^Li v x lJ 0 x 0 (3.59) JLi + dl + d^a j (z-_ daj JLi v The integrals are the same as those in (3.46) and thus can be calculated using (3.49)(3.50). Because attenuation and velocity are interpolated in the same manner, the factors preceding the integrals in (3.59) are equivalent to the analagous quantities in (3.53). 3.7 Summary Parameterization of the velocity model in terms of constant velocity gradient cells allows use of simple analytic expressions for kinematic (i.e. travel times and ray paths) and dynamic (i.e. amplitudes) ray properties, as well as for the travel time-velocity Jacobian. Similarly, constant attenuation gradient cells accommodate use of analytic formulae for amplitude attenuation. The formulae presented in this chapter compose the basis of the two-point ray tracing algorithm employed to construct the linear systems of equations used in the inversion scheme of Chapter 2. This inversion scheme is applied to data from the cross-ridge line in Chapter 4. Chapter 4 Tomography Applied to the Cross-Ridge Line 4.1 The Cross-Ridge Line The cross-ridge line is a single refraction profile, coincident with the multichannel reflection line of Rohr et al. (1988), which is centred on and normal to Endeavour Ridge; see Figure 1.1. Eight analog direct recording OBSs (Heffler & Barrett, 1979) with a passband of 3.5-35 Hz (3 db down) were deployed over a distance of 20 km along this line 3 to record shots from a 32 litre (2000 in ) airgun. The airgun was deployed at a nominal depth of 30 m and was fired at a pressure of 13.79 MPa (2000 psi) at intervals of ~0.2 km. The closely spaced array of receivers was designed to allow definition of lateral structural variations across the ridge, and to provide multipath ray coverage particularly at depths of 2-3 km beneath the ridge where a crustal magma chamber might exist. Continuous bathymetry records, which were subsequently digitized at distance intervals equivalent to the shot spacing, were collected along the line using a 12 kHz depth sounder. Water depths at the OBS deployment sites were determined from these records. Little or no sediment cover exists along the line. Navigation during OBS deployment and airgun shooting was provided by the BIONAV system (Grant, 1980) which integrates passive ranging Loran-C and Transit satellite (SatNav) navigation. The BIONAV positions, logged at 1 minute intervals, were postprocessed (i.e. edited, missing points interpolated and smoothed using a 3-point moving average) and corrected for drift using an edited subset of the recorded satellite fixes, and the OBS positions along the fine were adjusted using the direct water wave arrival recorded by the OBS. Relative shot positions for this line are accurate to ±15 m, while 91 Chapter 4. Tomography Apphed to the Cross-Ridge Line 92 relative OBS positions (after water wave adjustment) are accurate to ±20 m, resulting in shot-receiver distance uncertainties of ±25 m. The OBS analog data were digitized at 120 samples per second using a 10 Hz time signal from the data tape to control strobing of the digitizer. Such a procedure eliminates the effects of long period tape speed variations and tape stretching, so that timing errors during digitization are normally less than a sample interval (8 ms). The shot clock was corrected for drift using the WWVB radio time signal, resulting in origin time uncertainties of ± 3 ms. OBS clock ratings made at deployment and recovery were used to correct the OBS clock times for drift assuming a linear drift rate; accepting this assumption, the drift corrections introduce negligible errors. A timing correction was also made to account for the advance of the time channel with respect to the three data channels on the analog tapes caused by skew of the recording head. This correction accounts for the cumulative head skew of the complete recording and playback system. Further details of the digitization process, timing corrections and OBS repositioning can be found in Appendix A. 4.1.1 Projection onto a Line All shots and OBS locations over the 30 km length of the cross-ridge profile lie within ±300 m of a straight line fit to the line of shots (see Figure 4.1). To accommodate a D 2 analysis, the OBS positions are projected normally onto the straight line whereas the shot positions are projected onto the line at a position corresponding to the actual shot-receiver distance for each OBS. Thus a single shot will project onto the line at a slightly different position for each OBS. This procedure preserves the travel time-distance relationship for each shot-receiver pair. The bathymetry collected for the profile is also projected normally onto the line, and as a result, the bathymetry will be shifted relative to the projected shot positions. For shot-receiver distances of > 3 km, the relative shift is < 40 m and in most cases it is << 40 m. This shift is small relative to the sampling interval used for the bathymetry and thus in general it will introduce negligible Chapter 4. Tomography Applied to the Cross-Ridge Line Figure 4.1. Projection lines. The straight line approximation to the cross-ridge profile is shown, with the sub-parallel line to the north used for OBS 8. Open circles indicate OBS positions and dashed lines represent the actual shot lines. Reference positions (mentioned in the text) for the lines are indicated by the solid circles. 93 Chapter 4. Tomography Applied to the Cross-Ridge Line 94 errors. The straight-line fit to the line of shots passes through the points (47° 56.595' N , 128° 53.247' W) and (48° 1.151' N,129° 15.409' W); see Figure 4.1. The positions denoted in the figures to follow are measured relative to the latter point. Because the taperecorder for OBS 8 had stopped recording by the time the cross-ridge line described above was shot, data from an earlier cross-ridge line, offset from the above line by 0.5-1.5 km, were used for OBS 8. As described above, the shot positions and bathymetry for this line were projected onto a straight line passing through (47° 57.656' N, 128° 56.025' W) and (48° 1.774' N, 129° 14.676' W), (Figure 4.1). The latter point is the reference point from which positions along this line are measured. The reference point was chosen so as to match the projected bathymetry profiles along the twofinesas closely as possible (see Figure 4.2). 4.1.2 The Data Vertical component data sections for the OBSs are shown in Figures 4.3, 4.4, 4.5 and 4.6. The most prominent features observed on these sections are the direct water wave arrival, and the first arrivals. There is a large amount of secondary energy observed, although it is mainly restricted to events which mimic the first arrival, suggesting that these events are a variety of multiple. Secondary events following the water wave arrival near the source will not likely be observed, due to the persistent ringing caused by the direct water wave arrival and the limited dynamic range of these OBSs, which introduces a soft clipping (see Section A.4). The linear, high apparent velocity events which stand out on the OBS 3 section (Figure 4.4 (a)) are due to a small airgun that was being used as a source for continuous seismic profiling (CSP). The system gains of OBSs 2, 3, 6 and 7 exceed the gains of the other OBSs by approximately 6 db resulting in a somewhat different appearance for the two sets of seismic sections. This difference in appearance occurs because signals which are clipped due to the limited dynamic range of the OBS will be improperly scaled when the amplitudes are adjusted for differences in Chapter 4. Tomography Applied to the Cross-Ridge Line Figure 4.2. Cross-ridge bathymetry. The bathymetry for the cross-ridge line (solid) and the sub-parallel line (dashed) to the north (used for OBS 8) are compared. The sub-parallel bathymetry was used in the model when ray tracing for OBS 8. 95 Chapter 4. Tomography Applied to the Cross-Ridge Line 96 O CO jttk'l'l k 1*1'' i 1' > 'Hk>»J J k b j k k W * * CO CN CN I i'-'l : „i|,|. I.'I'I q I.J CO IMI,II,,in. 77 „ V-/ CN CN I iiK .'i»)i\N'!»!»!'» :! !.!» •I'll'' II ,i •in nil ,;, b 0 i i,, • n l . . . - 5 "i 10 r"'f 15 20 Position (km) Figure 4.3. Vertical component seismic sections for OBSs 1 and 2. Amplitudes are scaled proportional to distance and a zero-phase Butterworth filter has been applied to the data (passband 0-20 Hz for OBS 1, and 5-30 Hz for OBS 2). Only positive amplitudes exceeding a threshold level are plotted. The distances indicated are measured from the western end of the cross-ridge line. Triangles denote the OBS positions, and slopes associated with different apparent velocities (in km/s) are indicated in (b). Chapter 4. Tomography Applied to the Cross-Ridge Line 97 q 4) o °i - CN - 1 I q CN CN 10 15 20 Position (km) Figure 4.4. Vertical component seismic sections for OBSs 3 and 4. Amplitudes are scaled proportional to distance and a zero-phase Butterworth filter has been applied to the data (passband 0-30 Hz for OBS 3, and 0-20 Hz for OBS 4). Only positive amplitudes exceeding a threshold level are plotted. The distances indicated are measured from the western end of the cross-ridge line. Triangles denote the OBS positions, and slopes associated with different apparent velocities (in km/s) are indicated in (b). Chapter 4. Tomography Applied to the Cross Ridge Line 98 q oo rt^PMi 1 CN CN iiv^i;,';)>. "li''i';">! •! •'' 'V v •!•».»'•' i, l ' - , l M , " l l | l ; "a::::" > . 1 1 .' 1—r 00 <»>>>>• •r;^«^i^:^.W't^»i»^»»i III'IIII'' 1 i. 'IIIUMI'I' 1 'i 'in'nii'i »i 'iinniiu! CN CN : ; i: i;i:;!:;;!'; ;:. , M i i • i Ii I r—i—i—T-^T—i—i 0 10 1 1 ••..7*71 P l ' T " ! i . ' i r " i ' I i . ' , I , L H * ...1>»» r .-,> .. i " r—i—r—T* 15 20 25 30 Position (km) Figure 4.5. Vertical component seismic sections for OBSs 5 and 6. Amplitudes are scaled proportional to distance and a 0—30 Hz zero-phase Butterworth filter has been applied to the data. Only positive amplitudes exceeding a threshold level are plotted. The distances indicated are measured from the western end of the cross-ridge line. Triangles denote the OBS positions, and the arrow in (a) indicates the travel time delay associated with the summit depression. Slopes associated with different apparent velocities (in km/s) are indicated in (b). Chapter 4. Tomography Applied to the Cross-Ridge Line 99 q CO « ASfrSM^ >>V i <D «/> 1 O t 'in CN I'I'>> V 1 CN , i), ; I, 1 .»'|> 'M» I ; •l-i'-'l.T •;',), 1 , 'll-'l,,. ! :a l v i" ,k , • '•.••»i :>;: ">!•'•'' 'K •: 1 . I'll l'l'il'-, y "....",V,: ;'' i ! '1)1 . | l l, ll tt' I •Ii •iiu; "I ,1). I •• O B S > '•., • 7:;! !•• Iii' 1) q CO 5P U CN CN I B S 8 ''>'• i^^-^W'ti* J_» til 0 10 kJ 15 ••l.,'.iii 20 11 V j ..I'll 25 , 1' . i',., J''i' 30 Position (km) Figure 4.6. Vertical component seismic sections for OBSs 7 and 8. Amplitudes are scaled proportional to distance and a 0—30 Hz zero-phase Butterworth filter has been applied to the data. Only positive amplitudes exceeding a threshold level are plotted. The distances indicated are measured from the western end of the cross-ridge line. Triangles denote the OBS positions, and slopes associated with different apparent velocities (in km/s) are also indicated in (b). The shot interval of the data in (b) is half that of the other data sections. Chapter 4. Tomography Apphed to the Cross-Ridge Line 100 gain among the different OBSs. Special processing was applied to the data for OBS 2 in an attempt to remove signal noise caused by cross-talk with the time code channel. Also, relative to the other sections, the trace spacing for the OBS 8 data is decreased by one half as these data were recorded along a sub-parallel line with a doubled shot frequency. Timing problems with OBS 4 prevented the complete set of data from being used for the travel time inversion; data in the range 17-30 km have been excluded. The OBS 1 site (Figure 4.3 (a)) was very noisy, restricting observation of first arrivals on the OBS 1 data set to a range of < 12 km, and first arrivals were only clearly observed for positions > 6 km on the OBS 6 data section (Figure 4.5 (b)). In this analysis, we only consider first arrivals. The wavelet associated with the first arrivals remains remarkably consistent over the length of the profile, facilitating identification of this phase. The effect of seafioor topography (Figure 1.2 (a)) is particularly clear on the OBS 8 section (Figure 4.6 (b)) where bumps along the first arrival event are observed near 6 and 13 km. In fact, the time delay associated with the summit depression can be clearly seen on the OBS 5 section at 13 km (marked by arrow on Figure 4.5 (a)). Although variations in the seafioor topography complicate the shape of the first arrival curve, it can be characterized using three apparent velocities (if the direct water wave is excluded). Considering the data section for OBS 5 as a reference (Figure 4.5 (a)), then moving east from the OBS location at 16.57 km, the first arrival curve can be described as follows. From 16.6 km to 18.5 km, the direct water wave is the first arrival. A slightly higher apparent velocity (~2.5 km/s) is defined by the arrivals to 19.3 km, at which point the first arrival curve bends sharply to an apparent velocity of ~4.8 km/s. From 22.5 km to the end of the section, the first arrival curve defines an apparent velocity of ~6.0 km/s. The increase in apparent velocity after 26 km is observed on all of the data sections indicating that it is due to the change in the slope of the seafioor (Figure 1.2 (a)) or to shallow structure. The low velocity branch (~2.5 km/s) is best defined on the data section for OBS 7 which is shown in Figure 4.7. Chapter 4. Tomography Applied to the Cross-Ridge Line 0 1 2 101 3 4 Position (km) Figure 4.7. Low velocity arrival. T h e near-receiver data for O B S 7 (located at 0.0 km) demonstrate a low velocity arrival (i>=2.5 k m / s ) emerging from the direct water wave ( A ) between 2.0 and 2.8 km. Curves B and C are travel time curves calculated for layer 1 - 1 (with Vi=2.5 k m / s and V,u=0.5 s ) overlying layer 2 with u =4.8 k m / s , whereas D 2 _ 1 is the travel time curve for Vi=2.5 k m / s , but with V,t>=3.0 s . T h e dots indicate the travel time picks. Because of the sharp breakover indicated by the picked travel times, a layer with a low velocity gradient is favoured. T h e positions indicated are relative to the O B S position. Chapter 4. Tomography Apphed to the Cross-Ridge Line 102 The three velocities defined above were used to separate the first arrival picks for all eight seismic sections into three groups. First arrival travel times were determined on the data sections with the aid of a template. The template consisted of a transparent overlay on which was drawn a typicalfirstarrival wavelet chosen from the section. This was done to assure that consistent onset points were picked for the entire data set. Travel time uncertainties were estimated using the RMS travel time misfits obtained when the direct water wave arrivals were used to reposition the OBSs. This misfit should include timing, positioning and picking errors for each individual OBS, although it will not include systematic picking errors between different OBSs. Where the direct water wave occurs as the first arrival, the onset times are generally clearer than for any other arrival observed on the section. Thus, the errors in picking these arrivals should represent the minimum uncertainties involved in travel time picking. Estimated uncertainties for the complete set of travel times range from 8-25 ms. The RMS uncertainty for the entire set of travel times is 15 ms, based on these estimated uncertainties. 4.2 Travel Time Inversion The tomographic travel time inversion method described in Chapter 2 has been applied to the set of first arrival travel times determined for data from the eight OBSs. For a velocity model characterized by layering, the solution process is most effectively applied by starting with a subset of the travel times which define the shallow structure (i.e. those recorded at small distances from the OBS) and progressively expanding the number of travel times used to include those recorded at greater distances in order to define the deeper structure. The following sections describe the steps involved in applying this procedure to the data. Chapter 4. Tomography Apphed to the Cross-Ridge Line OBS 1 2 4 5 6 7 N o . of arrivals 4 5 2 3 4 6 103 v (km/s) 2.5 2.6 2.4 2.7 2.5 2.5 p Table 4.1. Seafioor velocities. 4.2.1 Starting Model Preliminary observation of the eight data sets identified three apparent velocities which generally characterize the first arrival curves. This suggests a starting model consisting of 3 layers. As with other deep water refraction studies using near-surface shots, observation of first arrivals from the first few hundred metres beneath the sea floor is limited. First arrivals associated with the shallowmost structure are observed at six of the eight sites across the ridge. The number of such arrivals observed at each site ranges from two to six, and the maximum distance range over which they are observed is ~0.8 km. Ray tracing was performed with the seafioor topography included to match the travel times of these arrivals, and a velocity of 2.5 ± 0.2 km/s includes the best-fit seafioor velocities for all six sites. A vertical velocity gradient of 0.5 s _1 was used in the model. This information is summarized in Table 4.1. Ray diagrams for these arrivals show that this velocity information is restricted to the upper 200 m beneath the seafioor, and the sparse receiver spacing (2-3 km) relative to the distances over which these first arrivals are observed (0.2-0.8 km) indicates that these velocities are independent samples along the fine. To actually define a layer characterized by v=2.5 km/s, we are extrapolating this information horizontally based on the consistency of the observed sample velocities. Note also that the data provide no independent determination of the shallow velocity on the axial ridge. Due to the general lack of ray coverage for shallow depths (and thus lack of constraint on lateral velocity variations), the structure of this layer was fixed in the Chapter 4. Tomography AppHed to the Cross-Ridge Line 104 inversion that followed. If this were not done, velocity perturbations which actually belong at depth would be introduced in this region due to the lack of constraint. To establish the deeper structure for a starting model, the travel times from OBS 5, for shots to the east of the OBS, were used to define a 3-layer structure. A trial-anderror method was used in which layer velocities and thicknesses were modified until the calculated travel times (determined by ray tracing) generally matched the observed travel times. Data for shots to the east of OBS 5 were chosen as the seafioor topography is more subdued there (Figure 1.2 (a)). Velocities determined in this manner for the surfaces of the three respective layers are 2.5, 4.8 and 5.8 km/s, and the vertical velocity -1 gradients within each of the layers are 0.5, 1.0 and 0.5 s . A comparison of observed and calculated travel times for a starting model in which the base of layer 1 was horizontal and a starting model in which it was a smoothed version of the seafioor, revealed that a better fit was obtained in the latter case. Thus, a bottom following interface was used in the starting model. The thicknesses of layers 1 and 2 for the starting model are ~0.45 km and ~0.80 km, respectively. A relatively low velocity gradient was assigned to layer 1 based on the rapid slope change of the travel time curve from an apparent velocity of 2.5 km/s to a velocity of 4.8 km/s (see Figure 4.7), as well as some preliminary amplitude modelling. Amplitude analysis of 15 seismic sections recorded during the same experiment (Cudrak, 1988) supports the existence of a low gradient in this layer. The reduction in gradient in layer 3 was required to reduce the first arrival amplitudes at distance. 4.2.2 Inversion for Shallow Velocity Structure The starting model described above was specified on a grid with x and z cell dimensions of 0.50 km and 0.20 km. The horizontal cell size was chosen to be of the same order as the shot spacing (typically ~0.2 km) and the smaller z cell dimension was chosen as vertical resolution normally exceeds horizontal resolution for surface-surface refraction Chapter 4: Tomography Apphed to the Cross-Ridge Line 105 studies (cf. Section 2.5). As we were initially concerned with the shallow structure, only the seafioor (interface 1) and interface 2 were included in the model (see Figure 4.8 (c)). Depths of these interfaces were specified every 0.50 km resulting in 63 depth parameters for interface 2. Although not shown, a separate seafioor interface was used for tracing rays to OBS 8 as the bathymetry along the parallel cross-ridge lines was slightly different (see Figure 4.2). A subset of 193 travel times (those identified with the intermediate velocity of 4.8 km/s) were used during this stage. The outstanding feature of the travel time residuals for the starting model (Figure 4.9 (a)) are the late calculated arrivals for OBSs 2, 3 and 4. Two iterations were performed to determine the geometry of interface 2. In each iteration, the interface depth perturbation vector 8z was smoothed by applying a 3-point moving average operator before updating the current interface. As the system was over determined, a damping factor fi—0.0 was used. The main result of the first two iterations is a 100-200 m shallowing of the interface to the west of the ridge axis to advance the calculated arrival times for sites to the west of the ridge (compare Figure 4.9 (a) and (b)). Also, the depth of the interface was increased slightly east of 20 km to delay the arrivals east of sites 6 and 7. The interface east of 24 km is unconstrained by these data. The RMS misfit for this set of travel times is reduced from 53 ms to 21 ms. Further reduction of the misfit could be obtained by continued iteration, but the degree of smoothing would have to be reduced resulting in a very rough interface. Determination of the interface geometry was done before considering velocity variations within layer 2 since very large velocity variations within the layer would be required to reproduce travel time changes accounted for by small changes in the interface geometry. This is due to the large velocity contrast across interface 2. Three iterations (iterations 3-5) were performed to determine perturbations in the velocity structure of layer 2, using a damping factor /J.=0.015. In each of these iterations the velocities in layer 1 and the depth to interface 2 were kept fixed, and the velocity perturbation was smoothed horizontally using a 3-point moving average operator. In Chapter 4. Tomography Applied to the Cross-Ridge Line 106 o CN Position (km) Figure 4.8. Shallow inversion results, (a) Ray diagram for iteration 2 model, (b) Ray diagram for iteration 5 model. Note how the rays are affected by the low velocity zone, (c) The form of interface 2 is shown for iterations 0, 2 and 7, along with the residual velocity for iteration 5. The residual velocity is (current model)— (background layered model). The intermittent reflector of Rohr & Milkercit (1988), with depths determined using velocities from the refraction velocity model,is also shown. The positions indicated are distances from the western end of the cross-ridge line. Chapter 4. Tomography Applied to the Cross-Ridge Line 8 7 6 5 • 11. mill T . 11 .11 • ..Hi •III Ill.....llll..ll A illinium ilfii" 3 2 107 1 2 0 0 ms • 1 ii'" •PL w 'i-'irin v* r "— 8 7 6 5 b .1. _ ill. -V 2 0 0 ms I l l in' •) I l „ ..In • | A • 1 tllnlln.li II nu | II • ' • • H I •Mill .1 1)1 • 'II..' 1 • 3 2 1 •• 1 ..Hi. ml .1 •••III T 8 7 6 5 "'1 . .1. 1 . .1 . -w •|i> •! 1 . . I . . -w • • • 1 1 • • M ' II. A .1 3 2 1 ...1 2 0 0 ms • , 11. Hi .11. " ,, , || " • ili^ii.. . iti.ll V ..1 "•III 1 1 1 1 1 1 1 1 1 1 10 1 1 1 1 1 15 1 1 1 1 1 1 20 1 1 1 25 Position (km) Figure 4.9. Shallow inversion residual map. Travel time residuals (t f, — t i ) are shown for (a) the starting model (RMS=53 ms), (b) the iteration 2 model (RMS=21 ms) and (c) the iteration 5 model (RMS=16 ms). OBS numbers are shown on the left and OBS positions are indicated by triangles. A subset of 193 travel times was included for the shallow inversion stage. Positive excursions on the diagram indicate that the calculated arrivals are early relative to the observed arrivals. The positions indicated are distances from the western end of the cross-ridge line. 0 s ca c Chapter 4. Tomography AppHed to the Cross-Ridge Line 108 iterations 3 and 4, smoothing was applied once to the perturbation in each iteration, whereas smoothing was applied twice in iteration 5. The main effect of these iterations is the introduction of a low velocity zone beneath the ridge which has an amplitude of —0.47 km/s and a width of ~3.0 km, measured at the —0.2 km/s contour (see Figure 4.8 (c)). The low velocity zone extends from 0.7-1.8 km depth sub-bottom. This feature is included to delay arrivals to sites 3, 4 and 5, reducing their RMS misfits from 22, 22 and 20 ms to to 11, 21 and 13 ms, respectively. Less significant velocity anomalies with magnitudes of < 0.1 km/s are also introduced to either side of the low velocity zone. The overall RMS misfit is reduced from 21 ms to 16 ms. The ray diagrams for the 2nd and 5th iteration solutions are shown in Figure 4.8 (a) and 4.8 (b). The convergence of rays from several OBSs in the anomalous zone beneath the ridge suggests that the velocity structure in this region should be well defined. To evaluate the resolution in this region, travel times were generated for a test anomaly, which were then used in the inversion procedure to reconstruct the velocity anomaly. The results for low level noise data (RMS< 2 ms) and noisy data (16 ms Gaussian noise added to the travel times) are shown in Figure 4.10 along with the test anomaly. The results indicate that the boundaries of the anomaly are well represented in the solutions, although the amplitude of the anomaly is reduced and smearing toward the bottom of the layer is more severe for the noisy data case. The degradation of resolution at depth due to noise is discussed in Section 2.5. This suggests that the amplitude of the sub-axial low velocity zone (as shown in Figure 4.8 (c)) represents a conservative estimate of the true amplitude. 4.2.3 Inversion for Deep Velocity Structure At this point the remainder of the travel times were included in the inversion (a total of 676 travel times). The velocity model was extended to 7 km depth by including an interface at 3.5 km depth (at x—0 km), with a velocity of 5.8 km/s at the interface 10 \A Position (km) 18 10 14 Position (km) 18 10 14 18 Position (km) Figure 4.10. Sub-axial resolution, (a) The velocity anomaly shown was added to a layered model and travel times generated for this model (employing the same source-receiver geometry as was used in the experiment) were used to reconstruct the anomaly. The results of inversion for low level noise data (<2 ms) and for noisy data (16 ms Gaussian noise added to the data) are shown in (b) and (c), respectively. The result in (b) was obtained after 3 iterations, whereas the result in (c) was obtained after 2 iterations. Chapter 4. Tomography AppHed to the Cross-Ridge Line and a vertical velocity gradient of 0.5 s _1 110 below, as determined for the starting model. The interface used was a smoothed version of interface 2. The travel time residuals for this model are shown in Figure 4.11 (a). The prominent features observed on the residual map are the negative residuals (i.e. calculated travel times are late relative to the observed) at distances > 23 km for all OBSs, as well as negative residuals in the position range 0-4 km for all OBSs, except 1 and 2, for which data are available in these distance ranges. As these trends are common to all of the OBSs, they are likely caused by shallow structure. A single iteration was performed for the geometry of interface 3. The solution produced little change to the interface except toward the edges of the model where there is less constraint. In these regions, large perturbations occurred. Because the velocity contrast is greater across interface 2 than across interface 3, the same effect on the travel times' can be achieved using smaller depth changes applied to the shallower interface. As this latter solution represents a smaller deviation from the starting model, the original interface 3 was restored and two iterations (iterations 6 and 7) were performed to determine changes to interface 2 (see Figure 4.8 (c)) using the entire set of travel times. As before, the perturbation was smoothed in each iteration using a 3-point operator and a damping factor /z=0.0 was used. The result is a shallowing of interface 2 in the 21-30 km range, and a reduction of the downward bulge near 4 km. These changes reduced the total RMS misfit from 35 ms to 24 ms. The effect of these modifications is demonstrated by comparing Figures 4.11 (a) and (b). The shallowing of the interface east of 23 km, advances the calculated travel times in this range for all receivers. The downward bulge in interface 2 near 4 km was introduced in iterations 1 and 2 based on a small set of arrivals from shallow depth. It is clear (see Figure 4.11 (a)) that the bulge delays the westernmost travel times for all sources, and thus it is reduced in size by iterations 6 and 7. Finally, two iterations (iterations 8 and 9) were performed for the velocity structure using the complete travel time set, reducing the RMS misfit from 24 ms to 17 ms. Chapter 4. Tomography Applied to the Cross-Ridge Line 8 7 6 5 A 3 2 f'J llljll ]lM i | " 'll'llll mi l "••IB .•I. Illii M 8 7 6 5 A 3 2 1 • -•IIHI' ••'•••Il W|||l'' - IP " ••" 1 ^ Ii • .... 1 ll •••ill. .11. .. i l •••in 200 ms 'l|l||H|| ~""H||llH|||l|p| . II • |||IM ..i.in.Mi... ~w •ll 11*1 ""^IPIll '"'|ll||| ""jl IIB •IIH V 1 IIHIIII . L i . 111 * " |l iriii|iir • i . n 'i II i ' " " i" . . . ii... '1 11' I 1 III..ill "1 —•••fp|| ..... | » - ' l | | | - • M - | | | i | | | •w 1 ••« • i | . IT ,| 1 nri'" A • -•••••-•II f I - I I l" ' - |l 1 1 1 1 , |l||||||l||l|||l ••• • • | l l | " 111 •• • III..ill IIIL •I | l I"" I.I. lllu. 200 ms u III.! - 1 ill b ,| .I111.1..111I nil.. ..in •'" "I'll"'"'!!^!!!!!! .1.. l l . T . 1" • III'" 11 1. . • . J l . . . . ...ill If 1 •' "' ^tiriii |iii|iiiiiny I . ,11. ... • - i ... Hi •"Ill • • |H III • " • | l .mini.. ..i i III' •1. •••••i|iti| 11 ...II. Il»ll| .(Ill I i IHII|||l||l|||i i'••III1 ill. 8 7 6 5 3 2 "'1 1 1 \ H«l "|llll|||l||l||l M | I"" l" ' . lim. • iiI M "i i i m i i m i i M i . . I"'*, ilillllll.. ..••In.,.. 1 IjII111"• " -• -| • ||riT" 1! |i 1 111 11 • ill 11 200 ms 'IT Ill -rp .Jl. . . 3E. 5 • .-,.„ ." „,| ,||| || .illii. I -U • I-II I ••I-. • I i II . 111 III |I-|'I||UI | | —"Iff • • " 10 till* —I 15 i . . i | 20 i 25 30 Position (km) Figure 4.11. Deeper inversion residual map. Travel time residuals (t i, — t a k ) shown for the complete set of 676 travel times for (a) the extended iteration 5 model (RMS=35 ms), (b) the iteration 7 model (RMS=24 ms) and (c) for the final model (RMS=17 ms). OBS numbers are shown on the left and OBS positions are indicated by triangles. Positive excursions on the diagram indicate that the calculated arrivals are early relative to the observed arrivals. The positions indicated are distances from the western end of the refraction line. a 0 s C r e Chapter 4. Tomography Apphed to the Cross-Ridge Line 112 Interfaces 2 and 3 were kept fixed as was the velocity in layer 1, and /i.=0.010 was used. The velocity perturbation in each iteration was smoothed using a 3 x 3-point moving average operator. The resultant total residual velocity is shown in Figure 4.12 (b). The residual velocity is obtained by subtracting the layered background velocity model from the final velocity model. The two most prominent features of the residual velocity model are the low velocity zone beneath the ridge which has been maintained, and a large high velocity zone (with amplitude 0.35 km/s) in the range 3-10 km. The high velocity zone lies in the depth range 1.5-3.2 km sub-bottom and advances calculated travel times in the range 0-6 km for sites 3-8 by as much as 50 ms (compare Figure 4.11 (b) and (c)). A smaller high velocity anomaly (amplitude < 0.15 km/s) is found in layer 3 in the range 23-28 km. For sites 3, 5 and 6, this anomaly advances travel times by up to 50 ms at ranges > 24 km. There are several smaller amplitude, smaller scale anomalies in layer 2. Note the small, large amplitude low velocity zone in layer 2 in the range 4-6 km. This feature which delays arrivals' for sites 1 and 2 has replaced the excessive downward bulge in interface 2 that was reduced when the full set of travel times was included. The inversion for deep structure was repeated using a horizontal interface as a replacement for interface 3, resulting in a solution very similar to the previous solution. Either model fits the travel time data equally well. 4.2.4 Travel Time Inversion Results To summarize, the velocity model resulting from the inversion procedure is essentially a 3-layer model. The uppermost layer has been assigned a velocity of 2.5 km/s with -1 a vertical velocity gradient of 0.5 s . This velocity is based on apparent velocities of 2.5 ± 0.2 km/s observed at six of the eight OBS sites across the ridge. Layer 1 varies in thickness from 250-650 m, but because the layer 1 velocity is fixed in the inversion, velocity variations within the layer will be represented as changes in layer thickness. Also, travel time bias errors (i.e. static shifts) for individual OBSs may show up as very local Chapter 4. Tomography Applied to the Cross-Ridge Line 0 10 20 113 30 Position (km) Figure 4.12. Travel time inversion results, (a) Ray diagram for the velocity model, (b) Residual velocity and the geometries for interfaces 2 and 3 for the complete model. The residual velocity is (current model) — (background layered model). Contours at the edges of the model are not closed as there is no constraint on the structure outside of the region of ray coverage. Also shown are the sub-axial (Rl) and mid-crustal (R2) reflectors of Rohr et al. (1988), converted to depth using velocities from the refraction velocity model. The positions indicated are distances from the western end of the cross-ridge line. Vertical exaggeration is 4:1. Chapter 4. Tomography Applied to the Cross-Ridge Line 114 variations in the interface geometry. The geometry of interface 2 is obtained primarily from shallow turning rays in layer 2, and thus the thickness of layer 1 will depend on having correct velocities in layers -1 and 2. Using the velocity range of ±0.2 km/s observed for the velocities sampled in layer 1, and assuming that the velocities used for layer 2 are correct within ±0.2 km/s, thickness variations of ±15-20% are expected. An incorrect gradient used for layer 1 will result in a scaling of the layer thickness. For instance, if the vertical velocity gradient of 0.5 s _1 is replaced with 4.6 s _1 (i.e. the gradient required for a velocity increase from 2.5 km/s to 4.8 km/s over 500 m) then a layer thickness of 500 m calculated using the former gradient will increase by 53% to 765 m. If instead, a zero gradient is used, then the layer thickness is reduced by 4% to 478 m. Thus, using a velocity gradient of 0.5 s -1 provides a conservative estimate of the layer thickness. Layer 2 is ~800 m thick, and is characterized by a velocity of ~4.8 km/s. A vertical velocity gradient of 1.0 s _1 was initially assigned to this layer. The lower boundary (interface 3) was chosen to follow interface 2, although the travel time data are fit equally well by a horizontal interface at 3.5 km depth. Velocity perturbations of —0.4 km/s to 0.1 km/s are superposed on the background velocity structure of layer 2, the most prominent of which is a 3 km wide low velocity zone beneath the ridge. This low velocity zone extends through layer 2 into the upper part of layer 3. Layer 3 has a background velocity of 5.8 km/s, upon which are superposed velocity perturbations of —0.2 km/s to 0.35 km/s. A velocity gradient of 0.5 s _1 was assigned to this layer. The maximum depth of ray coverage is ~4.0 km sub-bottom. The most prominent feature associated with layer 3 is a high velocity zone west of the ridge in the range 3-10 km. This velocity anomaly has an amplitude of 0.35 km/s and extends over a depth range of 1.5-3.2 km sub-bottom. There is a smaller high velocity zone east of the ridge in the range 23-28 km. This anomaly extends over a depth of 0.5 km just below interface 3 and has an amplitude of < 0.15 km/s. Note also, that there is no significant velocity anomaly beneath the ridge in the depth range 1.5 — 3.0 km sub-bottom, and Chapter 4. Tomography Applied to the Cross-Ridge Line 115 this region is well covered by rays travelling to various OBS sites (see Figure 4.12 (a)). The velocity perturbation of —0.2 km/s is an extension into uppermost layer 3 of the sub-axial low velocity zone from layer 2. 4.3 Amplitude Inversion The airgun source utilized in this experiment has a repeatable source signature and the OBS receivers all have similar response characteristics, making the data obtained suitable for quantitative amplitude comparison. However, the amplitude data are deficient in that the source and receivers have a narrow bandwidth, and more importantly, the recording instruments have a limited dynamic range which results in clipping of arrivals from nearreceiver shots. The high background noise level at some of the recording sites results in a further reduction of the dynamic range of the observed signals. The narrow bandwidth of the source-receiver system (in particular, the source) makes application of a spectral ratio method impractical. The power spectrum of a typical first arrival is shown in Figure 4.13. As an alternative, amplitudes are estimated using the seismogram envelope, and the amplitude-attenuation inversion scheme described in Section 2.2.3 is applied. As the problem is formulated, it is assumed that discrepancies between the observed amplitudes and amplitudes calculated for the velocity model are due to attenuation within the model, and the attenuation occurs along the ray path of the first arrival. The observed seismograms for the eight OBSs have been compensated for clipping as described in Appendix A. Without this compensation, meaningful results would not be obtained from the amplitude inversion. In this amplitude analysis, the objective is to account for the gross features of the observed first arrival amplitudes, rather than detailed trace-to-trace waveform matching. In light of this, the approximate method used for amplitude calculation (which includes smoothing of the model) is acceptable. 116 Chapter 4. Tomography Applied to the Cross-Ridge Line ° 0 10 Frequency 20 30 (Hz) Figure 4.13. Power spectrum for a typical first arrival. The power spectrum calculated for a 500 ms window is dominated by a peak near ~11 Hz which is associated with the bubble pulse frequency of the airgun source. Chapter 4. Tomography Apphed to the Cross-Ridge Line 4.3.1 117 Amplitude Estimation for Inversion Amplitude comparison requires some method of determining amplitudes of the first arrivals from both the observed and calculated seismograms. Amplitudes are estimated in one of the two following ways. The first method estimates the amplitude as the mean of the magnitudes of the maximum and minimum values of the seismogram within some time window following the arrival; Vest = [\s in\ + | « m o » | ) / 2 . (4.1) m The second method estimates the amplitude as the mean value of the seismogram envelope within the time window following the arrival. That is, tf- = 4 i > ; l 4 2 ( - ) where 7Y is the number of data samples within the window, and Sj are the data values. Obviously the amplitude estimates from these two methods will not be the same. It doesn't matter which method is used for amplitude estimation as long as the same method is used for both sets of amplitudes which are to be compared. For the inversion procedure, amplitudes for both the observed and calculated seismograms were determined using the envelope estimate with a 150 ms time window following the first arrival. The calculated seismograms were formed using a method similar to that employed by Spence et al. (1984). The impulse response of the model at a particular receiver is constructed using the calculated arrivals for that receiver and the final seismogram is formed by convolution of the impulse response with a chosen source wavelet. A single source wavelet (275 ms duration) chosen from the observed data was used for all synthetic seismograms. Chapter 4. Tomography Applied to the Cross-Ridge Line 4.3.2 118 Scaling of the Calculated Amplitudes Before forming the normalized amplitudes U ' in (2.23), an appropriate scaling factor is required to relate the relative amplitudes calculated by ray tracing to the observed amplitudes. An initial attempt was made to determine this scaling factor by comparing the observed amplitudes of the direct water wave arrival to amplitudes for this arrival calculated by ray tracing. A 250 ms time window following the arrival time of the direct water wave was used to calculate the observed amplitudes with (4.1). As can be seen in Figure 4.14, the comparison is valid only at large distances from the receiver as the arrivals from the near-receiver shots are clipped. In this comparison a scale factor of ~ 15000 (for receiver gains normalized to OBS 1) is required to match the observed and calculated arrivals. However, when the observed first arrival amplitudes are compared to the calculated amplitudes (both calculated using (4.1)), a scaling factor of ~29000 is required. The cause of this discrepancy is unknown. One possibility is that the source is directional. Source directionality may be caused in part by the superposition of the sea surface reflection with the airgun bubble pulse. For an airgun depth of 30 m, a time delay of ~40 ms exists between the sea surface reflection and the primary pulse travelling vertically downward. This delay corresponds to 0.44 of the bubble pulse period Tf, (Tj,=91 ms for /=11 Hz). Accounting for the phase reversal (i.e. a phase delay of Tb/2) of the sea surface reflection, the reflection and the bubble pulse will constructively interfere. As the shot position is horizontally offset from the bottom-receiver, the reflection-bubble pulse total phase delay gradually approaches Tb/2 (i.e. just the phase reversal due to reflection), at which point the reflection and bubble pulse destructively interfere. This directionality is consistent with the discrepancy in scaling factors because the first arrivals generally correspond to energy leaving the source at steep inclinations, whereas the direct water wave arrivals recorded at distance leave the source at shallower angles. Instead of using the observed water wave, the first arrival amplitudes were used to Chapter 4. Tomography Applied to the Cross-Ridge Line 119 o o. o « B o — o E < a OBS 1 e OBS 5 b OBS f OBS 6 c OBS 3 g OBS 7 o o. o o •»- O — o CL«o E < o o o T3 »• O — o E < Position (km) Position (km) Figure 4.14. Observed versus calculated water wave amplitudes. Calculated amplitudes (dashed) scaled by 15000, are compared to the observed amplitudes (solid) determined using (4.1). The curves match reasonably well at distance where the observed amplitudes are not clipped. No compensation is applied to the observed amplitudes plotted, as the clipping of the water wave is too severe. Chapter 4. Tomography Applied to the Cross-Ridge Line 120 determine a scaling factor. A suitable scaling factor was chosen which was minimal and which produced calculated amplitudes generally exceeding the observed amplitudes. These two criteria were used in choosing a scale factor so that a minimum amount of attenuation would be required in the attenuation solution, and so that the attenuation would generally be positive. A scale factor of 29000 was applied to the calculated seismograms prior to determining envelope amplitudes for the inversion. The relative scale factors determined (using the direct water wave arrival) for the different OBS vertical component seismometers agree very well with the relative gains of the OBSs suggesting that the bottom coupling at the different sites is consistent from site to site. 4.3.3 Amplitude Calculation for the Velocity Model Before calculation of amplitudes, the final velocity model from the travel time inversion was smoothed (see Section 3.2.1) in the following way. The residual velocity was subtracted from the background velocity field, smoothed using a 7 x 5-point ( 3 x 1 km, x X z) moving average operator, and added back to the background velocity field. Interface 2 was also smoothed using a 5-point (2 km) operator. The smoothed model is shown in Figure 4.15 (a). Figure 4.16 shows the OBS 7 travel time-distance curves and the amplitude-distance curves calculated for the smoothed velocity model. The following observations can be made. First, there are other travel time branches within a 150 ms window following the first arrival curve. Outside of the stippled region, there is a triplication near 15 km where the reflection from the base of layer 2 (r2) arrives following the turning ray in layer 2 or layer 3 (2 or 3). Thus, energy associated with this secondary arrival will contribute to the envelope-amplitude calculated for the first arrivals. Referring to the amplitude-distance plot, it is apparent that the amplitude associated with the turning ray (2 or 3, which is normally the first arrival) is larger than the amplitude of r2 away from the critical distance (near 15.5 km) where the turning ray amplitude goes Chapter 4. Tomography Applied to the Cross-Ridge Line 121 km/s E33 -0.10—0.20 VZZ 0.10-0.20 0.20-0.30 T O 1 1 r km/s E 2 1 -0.10—0.20 YZA 0.10-0.20 10 20 30 Position (km) Figure 4.15. Smoothed velocity models for amplitude calculation, (a) Residual velocity resulting from application of a 3 X 1 km (x x z) smoothing operator to the residual velocity of Figure 4.12. Interface 2 has also been smoothed using a 2 km operator, (b) Residual velocity resulting from smoothing (as in (a)) of the residual velocity in Figure 4.21. Chapter 4. Tomography Applied to the Cross-Ridge Line 122 Position (km) Figure 4.16. Example curves, (a) T h e travel time versus distance curve calculated for O B S 7 using the smoothed model of Figure 4.15 (a). T h e labelled branches correspond to: turning rays in layers 1, 2 and 3 (1, 2 and 3); reflections from the base of layer 1 and layer 2 ( r l and r2). The dashed line marks a 150 ms window following the first, arrival curve, and stippling identifies the region where the r l arrival and the water wave contribute to the envelope-amplitude, (b) Amplitude ing to (a). length. versus distance curves correspond- T h e amplitude-distance curves have been smoothed using a 2 k m operator Chapter 4. Tomography Applied to the Cross-Ridge Line 123 to zero, and thus will be the principal contributor to the envelope-amplitude. Hence, at most distances it is a reasonable approximation to assign variations of the first arrival envelope-amplitude to the earliest arriving ray which is the turning ray. This approximation is made in the inversion that follows. An exception to this occurs at small distances (< 4 km, marked by stippling in Figure 4^16) from the receiver where the direct water wave arrival (not shown) and the reflection from the base of layer 1 (rl) dominate the envelope-amplitude. Because of this and the severe clipping of the observed amplitudes in this region, these amplitudes are excluded from the inversion procedure. Also, note that the region about the critical distance (near 15.5 km) where the turning ray amplitude goes to zero is relatively small (1-2 km wide), and the largest amplitude arrival (with a similar travel time) in this region is the turning ray in the layer above. Thus, the inaccuracy of the ART amplitudes in this region should not significantly affect the envelope-amplitude calculated for this distance. However, observed versus calculated amplitude differences in this region may be improperly assigned during inversion to the lower turning ray rather than to the turning ray in the layer above, depending on which ray arrives first. 4.3.4 Application to the Data The amplitude-distance curves calculated for the smoothed velocity model (labelled Initial) are compared to the observed amplitudes in Figures 4.17 and 4.18. The noise levels shown in these figures were determined by calculating the envelope-amplitude for a 150 ms window preceding the first arrival. Travel times calculated for the final velocity model have been used as onset times for determination of observed amplitudes in regions where first arrival picks were not available. An exception to this was made for OBS 1 because of the generally higher noise level at this site. Although observed amplitudes. for OBS 1 are shown over the entire distance range in Figure 4.17, only those amplitudes for positions < 12 km were used in the inversion. Also, data have been included for OBS 4 Figure 4.17. Observed versus calculated amplitudes for OBSs 1-4. Observed and noise amplitudes were determined from the data using (4.2) with a 150 ms window following and preceding the first arrival, respectively. Amplitudes within the stippled zone were excluded during inversion. Initial and final amplitudes were calculated for the velocity models in Figure 4.15 (a) and (b), respectively. Calculation of the final amplitudes also included the attenuation model (Figure 4.19). In (a), only amplitude values for positions < 12 km were used. The observed amplitudes have units of 4.79 x l O km/s. - 1 2 o o. o 3 o o. OBS 5 OBS 7 a ::: :::::::: - /'A ^^::;:;:v:v:-:;\»s ' ^ \ /it s —•—• // * / / / A E < \y- * I/ ; ^ *^«,. + j • c \ j sjl / 3 ::: ::^::::::::!t-:::: : : : : : '•: /fe : : : >x-x< :: >-«-. •i ^^v.\•v.v.y.^^•.: : : 1 / If ^s^^ /. ' \ : : : : : : : : : : :7s :->;-: : : :-: ki: : : •' : •x:x : : :$; ;;::;*•i^ r . l , ~r-i i l l l 1 ' ' l l l i i l 1 10 20 a' 30 OBS 8 o Obtsrvod Vs o o 0) i 10 Position (km) 20 30 Position (km) Figure 4.18. Observed versus calculated amplitudes for OBSs 5-8. Observed and noise amplitudes were determined from the data using (4.2) with a 150 ms window following and preceding the first arrival, respectively. Amplitudes within the stippled zone were excluded during inversion. Initial and final amplitudes were calculated for the velocity models in Figure 4.15 (a) and (b), respectively. Calculation of the final amplitudes also included the attenuation model (Figure 4.19). The observed amplitudes have units of 4.79 x l O km/s. - 1 2 to en Chapter 4. Tomography AppHed to the Cross-Ridge Line 126 (positions > 17 km) which were excluded in the travel time inversion because of timing problems. This results in a total of 768 amplitude data as compared to the 676 data used for the travel time inversion. Amplitude data in the distance ranges marked with stippling in Figures 4.17 and 4.18. have been excluded for the amplitude inversion as they are dominated by the direct water wave arrival. Both the calculated and observed amplitudes have been smoothed using a 9-point (~2 km) moving average operator. The normalized amplitudes U ' -1 have an associated % =3360 and an RMS misfit of 13.8 ms-rad , where 2 an uncertainty of ± 2 5 % has been assigned to each of the 768 normalized amplitudes. This error estimate is designed to include uncertainties both in the observed amplitudes U and in the calculated amplitudes U J J . The expected RMS misfit value for uncertainties -1 of ± 2 5 % is 6.6 ms-rad . The normalized amplitudes are displayed in Figure 4.20 (a). In the inversion that followed, the attenuation in the uppermost layer was set to zero as the attenuation in this part of the model is poorly constrained. A centre frequency of / =12 Hz was used in the inversion. c The attenuation solution for /i=0.020 is shown in Figure 4.19, with the accompanying ray diagram. The solution was smoothed using a 5 X 3-point operator, with a resultant -1 % =1185 and an RMS misfit of 8.2 ms-rad . The most prominent feature of the at2 tenuation solution is the broad attenuative zone introduced below 4.3 km depth which is introduced to reduce the amplitudes of the most distant arrivals which are generally too large for each of the OBSs (see Figure 4.20 (a)). As an alternative to introducing this zone of attenuation, the vertical velocity gradient at depth can be decreased. This modification is preferred as the velocity gradient in Layer 3 of the oceanic crust is usually inferred to be less than that of the shallow crust (e.g. Gettrust, Furukawa & Kempner, 1982), and a reduction of velocity gradient at depth is consistent with the model of Cu- drak (1988). The attenuation anomaly at the very bottom of the model is based on rays from a single receiver, and thus it can be ignored. To include a fourth interface beneath which the vertical velocity gradient was reduced Chapter 4. Tomography Applied to the Cross-Ridge Line 0 10 20 127 30 Position (km) Figure 4.19. Initial attenuation solution, (a) Rays used to calculate attenuation solution in (b). (b) Attenuation solution for /z=0.020. Contours at the edges of the model are not closed as there is no constraint on the structure outside of the region of ray coverage. Also shown are the sub-axial (Rl) and mid-crustal (R2) reflectors of Rohr et al. (1988), converted to depth using velocities from the refraction velocity model. The positions indicated are distances from the western end of the cross-ridge line. Chapter 4. Tomography Applied to the Cross-Ridge Line 8 7 6 5 Him iiiiiiiiiniinnDi.....illih. niiinHllllllllllllllliilllini min 8 7 6 5 A ^ • • " •••ll« IIIHIIIIIIIIIIH II l l l l l l l l l l l l l l . l|| ^ llli"* ...^illlliii.. m - n..illlllll lllllllllllllllllhiiH. Mm,. " • UIHIIIM'' b ....III •iilllllllllllli. •• ^ .,„„ -ill.. ..- 1 III -...a II1 l l l t l 'Hill" ...I • M|||B||l'i|«PII||H . H.riMHlllllllllBlllllnillllllll.. ••• •w 1 .......iiimiuiiilinittiiiiiM... ...II..,. 'W .•llllllllllii III"" .iiiiiiiiiiiiin.. 1 T i l l II Illlllllllillllln. ill *llll« - V " i •||| 120 ''IIIIIIIIUIIMIII III" mnpM ..illllHlllllllllllna. ||ll, ••>•"••• r M | l . . ! • ! ( 111.. ..mi... • — ..irtlMfllllB.. . r ' .•llllllllllii. .IJ -w " " ' M l T 8 7 6 5 . IMMIUI ' ••"'"•I" •••'•II •i|||i' " U||||I ••ii-- -uiui|||||iii .••"I" •lulu' •••.••in v "Ml Mil 1 a. , U ..Illlll Nlllll- •l|||||pi 120 H|| ...llllllllllltll.il,.. III" Mil 3 2 1 120 inuumim • 3 2 1 A ull b llllllllllii tlllHIIIItii.ii A 3 2 1 _ ••II 128 III'f " " — n w m ' - ••IIIPI' " II •w II" ...HMII. "IIIIIIP" "III' -rftll t.Illlll ll«- 1 10 ,,i|P "iimiiiiiiiuiifii"••|liilIHI " 15 20 " —••-"•"^ii"" 25 30 Position (km) Figure 4.20. Normalized amplitude map. The normalized amplitudes U' calculated using (2.23) with / =12 Hz are shown for the complete set of 768 amplitudes for (a) the travel time-velocity solution in Figure 4.12 (RMS=13.8 ms-rad ), (b) the revised velocity model in Figure 4.21 (RMS=10.3 msrad ), and (c) for the final velocity-attenuation model (RMS=7.7 msrad ). OBS numbers are shown on the left and OBS positions are indicated by triangles. Positive excursions on the diagram indicate that the calculated amplitudes exceed the observed amplitudes. The positions indicated are distances from the western end of the cross-ridge line. c -1 -1 -1 Chapter 4. Tomography Applied to the Cross-Ridge Line to 0.3 s -1 from 0.5 s _1 129 above the interface, the iteration 7 velocity model (i.e. the model prior to inversion for deep structure in Section 4.2.3) was extended to include the extra interface. The additional interface was a smoothed version of interface 3 placed 0.80 km beneath interface 3. The travel time inversion for deep structure was repeated, with the resultant residual velocity shown in Figure 4.21. The travel time misfit values for the 2 2 extended model of % =2115 and an RMS misfit of 28 ms, were reduced to % =865 and an RMS misfit of 17 ms for the final model. A damping factor of /i=0.0075 was used in both iterations, and the velocity perturbation solution was smoothed using a 3 x 3-point operator in the first iteration, whereas in the second iteration a similar operator was used above interface 3, but a 3 x 5-point operator was used below interface 3. Comparison of this solution with the original velocity solution (Figure 4.12) reveals only minor changes. The pattern and magnitude of the shallow velocity anomalies in the two solutions are verj' similar as would be expected. In the revised velocity model, the deeper velocity anomalies are reduced somewhat in magnitude and are horizontally smeared due to the flatter attitude of the rays resulting from the decreased velocity gradient. The pattern of anomalies at depth, as in the original velocity solution, is asymmetric, but the asymmetry is less pronounced. Amplitudes for this revised velocity model were calculated as for the original velocity model and the amplitude—attenuation inversion was repeated using the same damping factor and smoothing operator that had been used in the original amplitude inversion. The resultant attenuation solution is shown in Figure 4.22, and the normalized amplitudes before and after the inversion are shown in Figure 4.20 (b) and (c). Reduction of the 2 vertical velocity gradient at depth reduced the amplitude misfit values to % =1846 and an -1 RMS misfit of 10.3 ms-rad , and for the attenuation solution they are further reduced 2 -1 to % =1032 and an RMS misfit of 7.7 ms-rad . As mentioned in Section 4.3.3, to accommodate the inversion process, amplitude differences at a particular receiver are assigned to the earliest arriving ray when in reality they may be associated with several Chapter 4. Tomography Applied to the Cross-Ridge Line to CN E o _C UO LO Figure 4.21. Revised velocity solution, (a) Ray diagram for revised velocity model in (b). (b) Revised velocity model with a 4th layer included. Contours at the edges of the model are not closed as there is no constraint on the structure outside of the region of ray coverage. Also shown are the sub-axial (Rl) and mid-crustal (R2) reflectors of Rohr et al. (1988), converted to depth using velocities from the refraction velocity model. The positions indicated are distances from the western end of the cross-ridge line. 130 Chapter 4. Tomography Applied to the Cross-Ridge Line 0 10 20 131 30 Position (km) Figure 4.22. Final attenuation solution, (a) Rays used to calculate attenuation solution in (b). (b) Attenuation solution for ^,=0.020. Contours at the edges of the model are not closed as there is no constraint on the structure outside of the region of ray coverage. Also shown are the sub-axial (Rl) and mid-crustal (R2) reflectors of Rohr et al. (1988), converted to depth using velocities from the refraction velocity model. The positions indicated are distances from the western end of the cross-ridge line. Chapter 4. Tomography Apphed to the Cross-Ridge Line 132 rays arriving in the time window used for amplitude determination. To test the validity of this approximation, theoretical seismograms were recalculated for the velocity model with the attenuation solution included. The complete set of arrivals was used to construct the seismograms, and then envelope amplitudes were generated for these seismograms. The resultant misfit was x =1086 and an RMS misfit of 7.9 ms-rad 2 -1 which compares favourably with the results obtained using the single ray approximation. The attenuation solution (Figure 4.22) is characterized by two zones of relatively high attenuation beneath the ridge axis in layer 2 and layer 3 as well as attenuative zones in layer 2 at 8 km distance to either side of the ridge. These attenuative zones correlate remarkably well with the shallow low velocity anomalies observed in the final velocity model (Figure 4.21). Attenuation values of 0.01-0.02 correspond to Q values of 50-100 which are of the correct order of magnitude for the shallow ocean crust (c/. Jacobson & Lewis, 1988). As well, there are zones of negative attenuation at either edge of the model where the velocity model is not as well constrained. Negative attenuation zones occur when the amplitudes calculated for elastic wave propagation through the model are less than the observed amplitudes. This may result, for instance, because the velocity gradient in these regions is too low. The remainder of the model has attenuation values of < 0.01. Due to the limitations of the observed amplitudes and the approximate method by which the theoretical amplitudes are calculated we refrain from assigning significance to attenuation variations with magnitude < 0.01. The scale factor relating the observed amplitudes to the calculated amplitudes was determined in an arbitrary manner. In light of this, the results obtained might be questioned. To determine the effect that the scale factor has on the attenuation solution, the solution was recalculated using amplitudes scaled by factors of 0.75-1.50 relative to the original scale factor. The patterns of the resultant attenuation solutions were very similar to the solution in Figure 4.22, but the magnitude of the anomalies varied. A scale factor Chapter 4. Tomography Applied to the Cross-Ridge Line 133 of 0.75 reduced the maximum positive attenuation values to 0.011 whereas the magnitude of the negative attenuation increased to —0.023. A scale factor of 1.50 increased the maximum positive attenuation to 0.047 (Q ~ 20) and decreased the magnitude of the negative attenuation to —0.018. Thus the pattern of attenuative regions observed in Figure 4.22 is valid, although the magnitude of the anomalies is not well constrained. 4.3.5 Amplitude Inversion Results The initial attenuation solution suggested a revised velocity model (Figure 4.21) in which the original layer 3 was subdivided into 2 layers; layer 3 which is ~800 m thick, with a -1 background velocity of 5.8 km/s at its surface and a vertical velocity gradient of 0.5 s , and layer 4 which has a background velocity of 6.3 km/s at its surface and a velocity _1 gradient of 0.3 s . Introduction of a fourth layer is based strictly on the amplitude inversion, as the travel times had previously been satisfied using a 3-layer model (see Section 4.2.4). The geometry of interface 4 is not well constrained, but was arbitrarily specified using a smoothed version of interface 3. Interface 4 marks a change in the vertical velocity gradient with little or no velocity change across the interface. The structure of layers 1 and 2 in the final velocity model are essentially unchanged from the original velocity solution (compare Figures 4.12 and 4.21), although the velocity anomalies below layer 2 have been changed somewhat. The magnitude of the high velocity anomaly to the west of the ridge in layers 3 and 4 has been reduced from 0.35 km/s to 0.25 km/s in the final velocity model. Also, the smaller high velocity anomaly to the east of the ridge which originally extended from 23-28 km, has been smeared out at depth. This is a result of the flattening of ray paths in layer 4 due to a lower velocity gradient, reducing the horizontal resolution of the solution in this region. The maximum depth of ray penetration is reduced from 4.0 km sub-bottom for the original velocity solution to 3.5 km for the final velocity solution. Chapter 4. Tomography Applied to the Cross-Ridge Line 134 The attenuation solution (Figure 4.22) is characterized by zones of relatively high attenuation (a ~ 0.01-0.02 or Q ~ 50-100) in layers 2 and 3 beneath the ridge axis, as well as in layer 2 at 8 km distance to either side of the ridge. The attenuation values associated with these anomalies should be considered as order of magnitude estimates, because of the arbitrary scaling applied to the amplitudes. The positions of these anomalies correlate well with low velocity zones in the velocity solution. Zones of negative attenuation occur at either edge of the model where the velocity model is not as well constrained, indicating that the velocity model does not successfully reproduce the observed amplitudes for rays propagating in these areas. The rest of the model is characterized by attenuation values of a < 0.01 or Q > 100. The attenuation solution was obtained with a=0 in layer 1. As a result, attenuation occurring along a ray in layer 1 will appear in layer 2 of the solution, and thus the attenuative zones in layer 2 of the solution may in reality extend into layer 1. One shortcoming of the solution method is that we have no means of assessing the effects of arbitrarily smoothing the velocity model as required for amplitude calculation. It is possible that the good correlation between the attenuation anomalies and the low velocity zones occurs because the amplitudes calculated for the smoothed model don't sufficiently account for the effects of a more localized, high amplitude velocity anomaly. This seems unlikely as there is not a general correspondence between velocity anomalies and attenuation anomalies in the solution, but elimination of this possibihty requires calculation of amplitudes using an alternative method such as Maslov theory which introduces smoothing in a more physically meaningful way. Chapter 5 Discussion of Inversion Results The final velocity and attenuation models resulting from travel time and amplitude inversion applied to data from the cross-ridge line are shown in Figure 5.1. These results and their significance are discussed below. 5.1 Layered Nature of the Velocity Model The final velocity model is essentially a 4-layer model. Layer 1 is generally thinner on the Pacific plate (refer to Figure 5.1 (a)) than on the Juan de Fuca plate (east of the ridge). The observed thinning of layer 1 toward the eastern edge of the model (on the Juan de Fuca plate) is likely significant, but the shallow structure is not as well constrained east of 24 km. Because the velocity of layer 1 was fixed during inversion, velocity variations within the layer will be represented as changes in layer thickness. To obtain apparent thicknesses of 250-650 m for a 400 m thick layer would require velocities ranging from 1.7-3.3 km/s. This wide range of velocities contrasts with the consistency of the shallow velocities sampled (see Table 4.1), suggesting that the variation in thickness is true rather than apparent. An exception to this may be in the vicinity of the ridge, where no independent determination of the shallow velocity exists. Rohr & Milkereit (1988) report a shallow intermittent reflector along a coincident multichannel reflection line (see Figure 1.1 for location). They obtain interval thicknesses of 600-1000 m and interval velocities which, except for a few outliers, cluster about 3.2 ± 0.5 km/s. These interval velocities and the refraction model RMS velocities of 2.6 ±0.2 km/s agree within uncertainties, although the reflection interval velocities are 135 Chapter 5. Discussion of Inversion Results i 0 i i i | i 10 i i 136 i | i 20 i i i 30 Position (km) Figure 5.1. Final velocity and attenuation models. Contours at the edges of the model are not closed as there is no constraint on the structure outside of the region of ray coverage. Also shown are the sub-axial (Rl) and mid-crustal (R2) reflectors of Rohr et al. (1988), converted to depth using velocities from the refraction velocity model. The velocity (km/s) and vertical velocity gradient ( s ) for each layer of the background velocity model are shown in (a) along the right side of the figure. -1 Chapter 5. Discussion of Inversion Results 137 systematically higher. To better compare the refraction and reflection models, the interval reflection times were converted to depths using velocities from the refraction model (see Figure 5.2). Interval times for the shallow reflector segments were determined on an unmigrated stacked section. Inspection of the common depth point gathers indicated that the onset time of the event was ~40 ms earlier than the stack enhanced event, and this was taken into account when picking. Uncertainties of ±15 ms have been assigned to the interval times determined from the reflection section. As shown in Figure 5.2, the profile defined by the sequence of reflectors follows interface 2 of the refraction model, although the reflection depths are systematically greater. The mean layer 1 thickness difference between the reflection and refraction results is 26 ± 20%. It is not clear whether this difference in thickness is significant or whether it is due to a constant bias in determining interval times from the processed reflection data. We assume the latter; thus the reflection interface and the layer 1- layer 2 interface from the refraction study are presumed to represent the same geologic feature. The segmented appearance of the reflector observed by Rohr & Milkereit, suggests that the abruptness of this transition zone varies along the line. Shallow layering is observed near other spreading centres. Most notably, Herron (1982) observed a shallow reflector on multichannel seismic (MCS) line 17 across the East Pacific Rise at 9° N, with an interval thickness of 400 m and interval velocity of 2.7 km/s. They interpreted the layer as a porous lava flow overlying the sheeted dike complex. McClain, Orcutt & Burnett (1985) have reported a 400 m thick low velocity gradient layer determined from a refraction line across the East Pacific Rise at 13° N. The velocity at the seafloor is 2.5-3.0 km/s and the layer is underlain by a rapid transition zone. Similarly, Harding et al. (1988) report a 100-200 m thick low velocity layer (2.4-2.6 km/s) in this region. Other refraction surveys near the East Pacific Rise (Ewmg & Purdy, 1982; Bratt & Purdy, 1984; Fischer & Purdy, 1986) ha reported shallow layers similarly characterized by low seafloor velocities (2.1-2.5 km/s), Chapter 5. Discussion of Inversion Results 138 F i g u r e 5.2. The shallow segmented reflector of Rohr & Milkereit (1988) is compared to interface 2 of the refraction model. The reflector depths are determined using velocities from the refraction velocity model. Chapter Discussion 5. of Inversion Results 139 _1 but having high gradients (3.5-4.6 s ) to sub-bottom depths of 700-800 m. Except for one refraction experiment on the Mid-Atlantic Ridge near 23° N (Purdy, 1987), such high gradients have been estimated using an extrapolation technique (Ewing & Purdy, 1982) assuming a constant vertical velocity gradient in the upper crust. Such a method is required as first arrivals from the shallowmost crust are not observed. Application of this method in a region where a low velocity gradient layer exists, would result in a large estimated gradient. For instance, application of Ewing & Purdy's method to travel time data calculated for a layer 500 m thick (seafioor at 2.5 km depth), with v=2.5 km/s and V i>=0.5 s 2 _1 _1 overlying a second layer with v=4.8 km/s and V v = 1.0 s , would result z in an estimated surface velocity of 1.7 km/s and a gradient of 4.75 s _1 for the layer. Thus, it is possible that shallow layering, as observed across Endeavour Ridge, may exist elsewhere but has not been recognized due to the difficulty of obtaining first arrivals from this layer in deeper water. The importance of porosity in determining seismic velocity in the shallow oceanic crust has been noted in previous investigations (e.g. Spudich & Orcutt, 1980; Christensen, 1984; Salisbury et al., 1985). Velocities are strongly dependent on both the bulk porosity as well as the geometry of pores within the medium. Figure 5.3 shows velocity vs. porosity curves calculated for a two-phase composite of seawater and basalt. Curves calculated for spherical and thin-disc inclusions using the self-consistent scheme (SCS, Budiansky, 1965; Hill, 1965) are shown, as well as the Hashin-Shtrikman bounds (Hashin & Shtrikman, 1963) which provide upper and lower velocity bounds. At low porosities, the effect of pore shape on velocity is dramatic. The Hashin-Shtrikman bounds assume no knowledge of the geometry of the constituent phases within the composite; only their fractional volumes are required (Hashin, 1970). If the upper oceanic crust is represented as a suspension of seawater inclusions within a basalt matrix (i.e. the phases are identifiable based on geometry) then the result obtained for SCS spheres may be a reasonable upper estimate of porosity. Neither the Hashin-Shtrikman bounds nor the Chapter 5. Discussion of Inversion Results O 140 Velocity vs. Porosity HS", SCS discs O iri J o O P 5" o CN 20 AO 60 80 100 Total porosity (%) Figure 5.3. Velocity vs. porosity. Shown are the Hashin-Shtrikman upper and lower velocity bounds ( i / 5 and HS~) calculated for a 2-phase composite material, where elastic moduli for sea-water and basalt have been used. Also shown are results calculated using the self-consistent scheme (SCS) for spherical and thin-disc inclusions of seawater in basalt. For these calculations zero-porosity compressional and shear wave velocities of 6.4 km/s and 3.5 km/s and a density of 3.05 gm/cc were used for the basalt. The compressional velocity is that measured at 6 kbar for a water saturated basalt sample from the Juan de Fuca Ridge area (Christensen, 1970). The lone experimental data point is from another Juan de Fuca Ridge basalt sample (Christensen, 1984). The dots mark where velocities corresponding to layers 1 and 2 of the refraction velocity model (t>!=2.5 km/s and u =4.8 km/s) intersect the various curves. + 2 Chapter 5. Discussion of Inversion Results 141 SCS results account for inter-pore communication of the fluid filled inclusions. A theoretical analysis of saturated cracked solids by O'Connell &; Budiansky (1977) indicates that velocity reduction due to the presence of saturated cracks is enhanced when there is inter-crack communication. Other factors (such as confining pressure, pore pressure and fining of voids with mineral alteration products) are effective in producing velocity changes primarily because they affect the porosity of the medium. Lab measurements on core and dredge samples indicate an increase of velocity with confining pressure which is attributed to the closing of voids and microcracks (Christensen, 1970; Fox, Schreiber & Peterson, 1973). Figure 5.4 illustrates velocity vs. pressure measurements for a water-saturated basalt sample recovered from the Juan de Fuca Ridge area (Christensen, 1984), demonstrating the dependence of velocity on the confining pressure as well as the pore pressure. The lower curve (pore pressure =hydrostatic) represents the situation where the shallow crust is permeable. In this case pore pressure acts to keep pores open, thus maintaining porosity and reducing the velocity. Experimental results (Todd & Simmons, 1972; Christensen, 1984) indicate only a weak dependence of velocity on confining pressure when the differential pressure (P nfining CO — Ppore) is maintained. For a medium with in- clusions having a wide range of aspect ratios, we would expect the small aspect ratio inclusions (i.e. thin cracks) to close first as the confining pressure is increased. As the velocity is more sensitive to the presence of low aspect ratio inclusions, a larger velocity gradient with respect to pressure is expected at low pressures, as is observed in Figure 5.4. A further consequence of the tendency for low aspect ratio inclusions to close first is that increased pore pressure is likely to be most effective in reducing velocity by maintaining cracks and fractures. Porosity changes are also associated with changes in lithology (e.g. basalt pillows—dike complex boundary, Becker et al., 1982) and the filling of fractures and voids with mineral alteration products (Anderson et al., 1982). The low velocities (2.5 ± 0.2 km/s) observed at the surface of layer 1, suggest that Chapter 5. Discussion of Inversion Results 142 Sub-bottom depth (km) 0.2 0.4 0.6 0.8 1.0 1.2 Confining pressure (kbar) Velocity vs. confining pressure for different pore pressures are shown as determined by Christensen (1984) f ° sample from the Juan de Fuca Ridge. The sample had 4% porosity at P =atmospheric. F i g u r e 5.4. ra c Chapter 5. Discussion of Inversion Results 143 the shallow igneous crust is highly porous as has been observed for 6 m.y. old oceanic crust at DSDP hole 504B (Becker et al, 1982). Porosities of 10-48% are obtained from the Hashin-Shtrikman lower bound and the SCS spheres curve (Figure 5.3). Porosities determined for dredge and core basalt samples have average porosities of 6-8% (Hyndman & Drury, 1976; Salisbury et ai, 1985), whereas in situ bulk porosities for the shallow ocean crust have been measured as 12-14% (Becker et al,, 1982) and as high as 40%) (Kirkpartick, 1979). To obtain velocities as low as 2.5 km/s for bulk porosity values similar to the more conservative of these observed values, much of the matrix void must exist in the form of low aspect ratio inclusions (i.e. cracks and fractures). Based on the large estimated porosities for layer 1, it is reasonable to suggest that layer 1 may be permeable as well. In this case (Ppore—Phydrostatic), the velocity gradient of 0.5 s _1 in layer 1 can be entirely accounted for by an increase in differential pressure with depth (see Figure 5.4).. For a permeable crust, the differential pressure at some depth z in layer 1 is the difference between the overlying water column to depth z and the combined pressure of the water column to the seafioor and the overlying rock mass. The large velocity increase from 2.6-2.8 km/s to 4.8 km/s across interface 2 suggests an abrupt decrease in porosity occurs at this depth (assuming the crack geometry remains the same). Porosity bounds of 0-27% are obtained from the Hashin-Shtrikman lower bound and the SCS spheres curve (Figure 5.3). There also may be a reduction in permeability at this level as is observed in DSDP hole 504B where permeability is reduced by an order of magnitude over the upper 200 m (Anderson et ai, 1982). A reduction in permeability so that P poTe —> 0, would result in a velocity increase of ~0.4 km/s ( see Figure 5.4). The velocities observed in layer 2 (in this study) range from ~4.4-5.7 km/s which compare well with interval velocities for pillow basalts from DSDP hole 504B, whereas interval velocities for the dike unit at hole 504B range from 5.5-6.0 km/s (cf. Figure 6 of Salisbury et ai, 1985). Based on the thickness of layer 1 Chapter 5. Discussion of Inversion Results 144 and the velocities observed in layer 2, interface 2 likely represents a porosity-related transition within the pillow basalts as opposed to the transition from the extrusives to the dike complex. Rohr & Milkereit (1988) have suggested that the shallow reflector observed in their study (presumably interface 2 of our refraction model) is due to a metamorphic front within the pillow basalts. Analysis of hydrothermal alteration at DSDP hole 504B indicates that metamorphic fronts exist within the pillow basalt sequence associated with upwelhng hydrothermal fluids at the spreading axis (Alt et al, 1986). Axial alteration 3 4 occurs on a time scale of 10 — 10 yrs, and the geometry of the metamorphic front is controlled by patterns of hydrothermal circulation. Using this model, the increase in velocity across interface 2 is attributed to replacement of low velocity seawater with higher 3 4 velocity alteration products within the voids. The estimated time scale of 10 -10 yrs required for formation of such a front would allow sufficient time for this layer to develop 4 beneath the ridge, if the last volcanic episode occurred ~ 10 yrs ago as suggested by Davis et al. (1984). McClain et al. (1985) attribute a low velocity surface layer on the East Pacific Rise to high porosity which they suggest may be due to surface Assuring which occurs as the crust spreads. They observe high surface velocities on the axial rise suggesting that the layer is absent there (i.e. no Assuring has occurred). Such a velocity increase might be expected across Endeavour Ridge as submersible mapping indicates that fracturing on the ridge itself is restricted to the summit depression (Kappel & Ryan, 1986). But, in contrast to the model of McClain et al. (1985) we do not observe an increase in the velocity of layer 1 associated with the ridge. As the velocity in layer 1 was fixed during the inversion procedure, a region of high velocity associated with the ridge would be manifest as a thinning of layer 1 beneath the ridge. This is not observed, making this explanation for the shallow layering less attractive for the Endeavour Ridge region. Also, if the thickness of layer 1 represents the depth to which surface fissures penetrate as Chapter 5. Discussion of Inversion Results 145 suggested by McClain et al, it is not obvious what controls this depth. Our results are more consistent with interface 2 representing a metamorphic front within the pillow basalts. Interface 3, which represents a transition from the higher gradient layer 2 above to a lower gradient below, and an increase in velocities from 5.5-5.7 km/s to 5.8 km/s, is likely the extrusives —dike complex boundary. Layer 3 is ~800 m thick. The velocities in the upper part of the layer are consistent with the interval velocities logged for the upper section of the dike complex penetrated at DSDP hole 504B (as given above). Interface 3 has been modelled as a velocity discontinuity, but the small velocity increase may occur over several hundred metres. In DSDP hole 504B this transition occurs over 200 m and it is characterized by a decreasing ratio of pillow basalts to dikes with depth. _1 Interface 4, which marks a further decrease in the vertical velocity gradient to 0.3 s , can be associated with the oceanic Layer 2-Layer 3 boundary. The geometry of this interface has been arbitrarily chosen as it is only constrained by the amplitude data. Velocity versus depth profiles constructed for the Samail ophiolite (Christensen & Smewing, 1981; Kempner & Gettrust, 1982) show a reduction in velocity gradient (i.e. presumably the oceanic Layer 2-Layer 3 boundary) which occurs near the dike complex-massive gabbro transition. The gradient decrease may be associated with the change in lithology, or with a metamorphic boundary related to sea water penetration (Christensen & Smewing, 1982; Salisbury & Christensen, 1978). The Samail ophiolite is thought to represent a section of young ocean crust (< 5-8 m.y. pre-obduction age) formed at a spreading rate of 410 cm/yr (Kempner & Gettrust, 1982), suggesting that a valid comparison can be made with the results for the Juan de Fuca Ridge. The layered model described here is consistent with the 4-layer model obtained by Cudrak (1988), who analyzed 15 in-line airgun profiles from the SEISRIDG-85 experiment using a trial-and-error forward modelling method to match travel times and amplitudes. Of particular interest, on most profiles Cudrak observed the refracted arrival in layer 1 as Chapter 5. Discussion of Inversion Results 146 a secondary arrival to distances of ~ 5 km. This provides very strong evidence that layer 1 has a low velocity gradient. Also, with few exceptions, the seafloor velocities observed by Cudrak were within the range of 2.5±0.2 km/s, in agreement with the results for this analysis of the cross-ridge line. 5.2 Magma Chamber, Velocity and Attenuation Anomalies Following from the fact that no large travel time delays or shadow zones are observed in the expected pattern (see Section 1.3 entitled Magma chamber detection), and no velocity anomaly is observed in the appropriate depth range beneath the ridge, we conclude that a partially molten magma chamber of substantial size does not exist in the depth range 1.5-3.0 km beneath the ridge. The region of the model (Figure 5.1) to 2.5 km beneath the axial ridge has a high ray density as well as wide angular coverage of rays (see Figure 4.12), indicating that it should be well resolved. Based on the earlier discussion, a magma chamber with dimensions of ~ l x l km (using v magma =4.5 km/s) represents the minimum chamber size that would be clearly identified using travel time delays. In the depth range of 1.5-3.0 km beneath the ridge, velocities vary by < ±0.10 km/s from the background velocity field. This eliminates the possibility that a broader low velocity zone (with velocities reduced by > 0.5 km/s) associated with hot rock surrounding the magma chamber exists. Such a zone was found beneath the East Pacific Rise at 13° N (Harding et al, 1988). We cannot rule out the possibility that a magma chamber exists below 3.0 km depth sub-bottom. The 3-4 km wide low velocity zone in layer 2 beneath the ridge has velocities which are reduced by < 0.4 km/s relative to the background velocity field (Figure 5.1 (a)). As layer 2 has been identified as part of the extrusives sequence with porosities of 0 — 27%, this anomaly is readily accommodated by a porosity increase of a few percent in this region. The increase in porosity may be due to increased fracturing in this zone or the low velocity zone may represent a zone of increased permeability within layer 2, . Chapter 5. Discussion of Inversion Results 147 resulting in greater pore pressures. An increase in pore pressure to hydrostatic could result in a velocity decrease of up to 0.4 km/s (Figure 5.4). Either of these effects may be associated with a zone of sub-axial hydrothermal circulation. The other smaller scale velocity variations in layer 2 may be real, but as they are less important in fitting the observed travel times, less significance is attached to them. The correlation of positive attenuation anomalies and low velocity zones in layers 2 and 3 (Figure 5.1) is consistent with the interpretation of the low velocity zones in terms of increased porosity and/or permeability. Like velocity, the intrinsic attenuation of a medium is strongly affected by the presence of cracks. An increase in attenuation with porosity is inferred from lab measurements which document decreased attenuation as confining pressure is increased, presumably due to the closing of cracks (Winkler & Nur, 1982; Toksoz, Johnston & Timur, 1979). Attenuation is accentuated if the cracks are fluid filled. At seismic frequencies, the primary attenuation mechanism in a saturated cracked medium is energy dissipation resulting from wave-induced fluid flow between cracks (O'Connell & Budiansky, 1977). Also like velocity, the attenuation of a cracked medium appears to depend on the differential pressure (i.e. P nfining~PpoTe) CO rather than the confining pressure (Johnston, 1981). Thus, increased pore pressure, which is effective in maintaining low-aspect ratio cracks and thus lowering velocity, will also be effective in increasing attenuation. The correlation of the attenuative zones and low velocity zones in layer 2 to either side of the ridge provides independent evidence which suggests that these smaller amplitude low velocity zones are real. These anomalies may be associated with off-axis hydrothermal circulation in the shallow crust. Analysis of heat flow data from the Juan de Fuca Ridge system by Davis et al. (1980) suggests that ventilated convection of water in the young ocean crust is common with horizontal convection cell dimensions ranging from 1-3 km to 10-20 km. Large scale Assuring observed in the summit depression and at the base of the axial ridge flanks (Tivey & Johnson, 1987) would provide conduits for intake or discharge of convecting fluid. Based on Chapter 5. Discussion of Inversion Results 148 the spreading model of Kappel & Ryan (1986) (see Section 1.1), a similar venting point might be expected at the base of the ridge to the west of the axial ridge (near 5 km, Figure 5.1). The pattern of velocity anomalies in layers 3 and 4 is asymmetric. Although regions of increased velocity are found below interface 3 toward both edges of the model, most of layer 3 and layer 4 west of the ridge have an associated positive velocity anomaly where velocities are increased by < 0.25 km/s relative to velocities beneath the ridge and to the east. This pattern could of course be explained in terms of porosity differences if the voids were of the correct geometry (z.e. thin cracks), but another explanation is preferred since the results from DSDP hole 504B indicate that porosities in the dike complex are very small (Becker et al., 1982), and the attenuation solution determined here indicates generally low attenuation in layers 3 and 4. A zone of moderately elevated temperature at depth beneath the ridge axis might be expected as a remnant from the last episode of crustal magmatic injection, or from the generally greater vertical heat flux associated with spreading centres. The increased temperature beneath the ridge would reduce velocities relative to the surrounding crust, but the observed asymmetry would not be expected. Rohr et al. (1988) observed asymmetry in the deep crustal structure across the ridge which they attribute to the presence of a heat source in the lower crust or upper mantle to the east of the spreading axis. Such a source would act to increase temperatures in the overlying crust and thus reduce velocities east of the ridge, resulting in the asymmetric velocity pattern. Hughes & Maurette (1957) observed average temperature coefficients o of —0.20 to —0.40 km/s/100 C for gabbro samples at 1 kbar, for temperatures between 25° C and 400° C. For temperatures between 300-500° C, Christensen (1979) observed temperature coefficients of —0.08 to —0.13 km/s/100° C for a gabbro sample at 2-3 kbar. Using the two intermediate values of temperature coefficient as conservative estimates o ( — 0.13 and —0.20 km/s/100 C), temperatures in layers 3 and 4 beneath and to the east of the ridge may exceed temperatures to the west of the ridge by 125°-200° C. This Chapter 5. Discussion of Inversion Results 149 anomalous thermal structure may be caused by interaction of the spreading centre with seamount volcanism (Karsten et al, 1986). The origin of the sub-axial reflector (Rl, Figure 5.1 (a)) observed by Rohr et al. (1988) remains enigmatic. It is possible that it represents a reflection from the roof of a very narrow magma chamber that is unresolved by the refraction data, but this is highly unlikely as the velocity anomaly associated with the hot rock surrounding a zone of partial melt should be resolvable. Another possibility is that the reflector is related to the sub-axial low velocity-high attenuation zone which extends into layer 3. If this is a region of hydrothermal circulation as proposed, then the observed increase in velocity at the base of this zone will be accompanied by an increase in density (due to decreased porosity). The impedance contrast due to this combination may be sufficient to produce the resulting reflection. The bulk density of a water saturated basalt sample with 10% porosity is 2.85 gm/cc compared to 3.05 gm/cc for a zero-porosity sample. Using a velocity increase of 0.4 km/s (from 5.6 to 6.0 km/s) with the above densities results in a reflection coefficient of 0.07. If the porosity contrast and velocity contrast are considered separately, the reflection coefficient is reduced to 0.03 in either case. The reflection coefficient for the combined effect (0.07) is large enough to produce an observable event. Chapter 6 S u m m a r y & Conclusions The two primary objectives of this study have been achieved: (1) a tomographic inversion D procedure for determining 2 velocity and attenuation structure from seismic refraction travel times and amplitudes has been developed, and (2) the upper crustal structure beneath the Endeavour segment of the Juan de Fuca Ridge has been determined by applying the inversion method to refraction data recorded along a cross-ridge line. The travel time-velocity inversion procedure is based on an iterative solution of the linearized travel time-velocity problem and allows for determination of continuous velocity variations as well as for the geometry of interfaces across which the velocity is discontinuous. In a particular iteration, either a velocity perturbation or an interface perturbation is solved for. Amplitudes are inverted directly for attenuation (i.e. no iteration is required) once the velocity structure of the model has been determined. The inversion method for both travel times and amplitudes is designed for surface-to-surface source-receiver geometries and successful application of the method requires a good initial estimate of the background velocity structure. The number of iterations required to achieve convergence for a set of travel times, depends on the degree of accuracy with which they are to be fit and on the degree of damping applied in each iteration. In synthetic examples using low noise data, the first 3-4 iterations reduced the travel time misfit most significantly and provided a good representation of the velocity anomaly. Further iteration reduced the misfit slightly, and defined the detail of the anomaly. In applying the inversion procedure to travel times from the cross-ridge data, 2-3 iterations typically were required to determine either interface geometry or velocity structure. The velocity and attenuation model is parameterized in terms of constant gradient 150 Chapter 6. Summary & Conclusions 151 (velocity and attenuation) cells allowing the use of analytic expressions for travel times, ray paths, amplitudes, attenuation and inversion quantities. Test examples and resolution analysis indicate that velocity and attenuation inversion using only first arrivals (i.e. turning rays) produces solutions that suffer from horizontal smearing due to the predominantly horizontal attitude of turning rays. Horizontal smearing can be reduced by including reflection arrivals (if available) and by increasing the number of source-receiver pairs. With regard to the latter, resolution is improved most by increasing the number of elements in the source-receiver geometry of which there are fewest. Streaks in the velocity or attenuation solution caused by inhomogeneous ray distribution are removed by parameter weighting. Amplitudes, calculated using zero-order asymptotic ray theory for the velocity model obtained from the travel time inversion, are used with the observed amplitudes to invert directly for attenuation within the medium. Parameterization of the model in terms of constant velocity gradient cells introduces artificial second-order velocity discontinuities across cell boundaries, which results in scattering of the calculated amplitudes. An ad hoc approach in which the velocity model is smoothed prior to amplitude calculation has been adopted to compensate for this scattering effect. The resulting amplitude- distance curves are also smoothed. This approach is acceptable for determining the gross amplitude vs. distance characteristics of the data, but would not be suitable for more detailed amplitude analysis. However, the method used for amplitude calculation can be readily incorporated into a Maslov-based algorithm (Chapman, 1985), if more accurate amplitude determination is required. The primary objective of the SEISRIDG-85 experiment, as relates to the work presented here, was to define the shallow crustal structure beneath the Endeavour segment 0 of the Juan de Fuca Ridge near 48 N, 129° W; and in particular, to image a magma chamber if one exists. To accomplish this, eight OBSs were deployed over a distance of 20 km on a cross-ridge line to record shots fired at 200 m intervals. The multi-path ray Chapter 6. Summary & Conclusions 152 coverage provided by this source-receiver geometry makes these data suitable for application of the tomographic inversion scheme that was just described. Amplitude comparison for data from the various OBSs is valid due to the repeatability of the airgun source and the common response characteristics of the eight OBSs, although amplitudes for nearreceiver shots required compensation for clipping. A total of 676 first arrival travel times and 768 envelope-amplitudes were determined from the data to use in the travel time and amplitude inversion, respectively. A layered starting model was determined for the inversion procedure by trial-and-error forward modelling of first arrivals east of OBS 5. Velocity gradients for the starting model were determined by comparison of calculated and observed amplitudes for the data from OBS 5. The final model obtained from travel time and amplitude inversion applied to data from the cross-ridge line is essentially a four-layer model, with significant lateral variations of velocity and attenuation. The uppermost layer has been assigned a velocity of 2.5 km/s 1 with a vertical velocity gradient of 0.5 s" , based on apparent velocities of 2.5 ±0.2 km/s observed at six of the eight OBS sites across the ridge. Layer 1 varies in thickness from 250-650 m. The velocity and attenuation (a = 0.0) of this layer were kept fixed during inversion. Layer 2 is ~800 m thick, and is characterized by a velocity of ~4.8 km/s. A vertical velocity gradient of 1.0 s _1 was initially assigned to this layer. The lower boundary (interface 3) was chosen to be a smoothed version of interface 2, although the travel time data are fit equally well by a horizontal interface at 3.5 km depth. Velocity perturbations of —0.35 km/s to +0.10 km/s are superposed on the background velocity structure of layer 2, the most prominent of which is a 3-4 km wide low velocity zone beneath the ridge. This low velocity zone extends through layer 2 into the upper part of layer 3 and correlates well with zones of relatively high attenuation (a ~ 0.01-0.02 or Q ~ 50-100) located beneath the ridge. Smaller magnitude low velocity zones in layer 2, which also correlate well with attenuative zones, exist at a distance of ~ 8 km to either side of the ridge. Chapter 6. Summary & Conclusions 153 Layer 3 is ~ 800 m thick with a background velocity of 5.8 km/s, whereas layer 4 has a background velocity of 6.3 km/s. Layers 3 and 4 were initially assigned vertical velocity _1 gradients of 0.5 and 0.3 s , respectively. The geometry of interface 4, separating layers 3 and 4, was chosen arbitrarily as a smoothed version of interface 3, as it is constrained only by the amplitude data. There is little or no velocity change across this interface. Velocity perturbations of —0.30 km/s to 0.25 km/s are superposed on the background velocity structure of layers 3 and 4, the most prominent of which is a high velocity zone west of the ridge in the range 3-10 km. Except for an attenuative zone beneath the ridge in layer 3, layers 3 and 4 are characterized by attenuation values of a < 0.01 or Q > 100. There is no significant velocity anomaly beneath the ridge in the depth range 1.5 — 3.0 km sub-bottom. Layer 1 and layer 2 are inferred to comprise the extrusives sequence, whereas layer 3 and layer 4 represent the dike-complex and underlying massive gabbro sequence. Layer 1 and layer 2 are separated by a sharp transition zone (interface 2) across which bulk porosities are reduced from 10 — 48% to 0 — 27%, and which correlates well with a shallow segmented reflector observed on the coincident reflection line of Rohr & Milkere.it (1988). This abrupt transition may be associated with a metamorphic front within the pillow basalts (Rohr & Milkereit, 1988) or it may represent the depth to which surface fissures penetrate (McClain et al., 1985). Our results favour the former interpretation. Interfaces 3 and 4 represent transitional boundaries with small or zero velocity changes, separating regions of differing vertical velocity gradients. The data preclude the existence of a region of partial melt exceeding maximum dimensions of ~1 x 1 km in the depth range of 1.5-3.0 km beneath the ridge axis. Also, no broad low velocity zone associated with hot rock surrounding a magma chamber was observed, further limiting the size of a hypothetical magma chamber. However, the possibility of a magma chamber residing at greater depth cannot be excluded. The absence of a significant shallow level magma chamber suggests two possibilities. Either a large shallow crustal magma chamber is not required Chapter 6. Summary &: Conclusions 154 during the spreading process, or Endeavour Ridge is currently experiencing a period of diminished magma supply with the magma chamber reduced or ehminated by extensive hydrothermal cooling. The development of the summit depression and the observed lack of recent volcanism are consistent with the latter scenario. Crustal temperatures remain elevated as evidenced by the observed discharge from hydrothermal vents, and the lower velocities in layers 3 and 4 beneath the ridge. The 3-4 km wide low velocity zone (velocities reduced by < 0.4 km/s) observed in layer 2 beneath the ridge, and extending into layer 3, is consistent with zones of relatively high attenuation (Q ~ 50-100) observed beneath the ridge. The low velocityhigh attenuation anomaly can be accounted for by an increase in porosity in this region due to an abundance of fracturing and/or increased permeability resulting in increased pore pressure. Either of these effects may be associated with a zone of hydrothermal circulation beneath the spreading centre. The smaller amplitude low velocity anomalies 8 km to either side of the ridge, with corresponding attenuation anomalies, indicate that hydrothermal circulation is not restricted to the region directly below the ridge. Asymmetric distribution of velocity anomalies below layer 2 suggest crustal temperatures in layers 3 and 4 may be elevated by 125°-200° C beneath the ridge and to the east. The cause of the asymmetry may be a deep crustal or upper mantle melting anomaly east of the ridge. The results of this study demonstrate the potential of using multi-source—multiD receiver refraction profiles in conjunction with seismic tomography to determine 2 velocity and attenuation structure. The structural information obtained complements that provided by a coincident reflection survey. Ideally, refraction and reflection data should be combined during inversion to take advantage of the contrasting resolution characteristics of refracted and reflected ray paths. This combination would be most suitably implemented using the same source and receiver array to collect both refraction and reflection data, in order that discrepancies caused by differing source-receiver characteristics might Chapter 6. Summary h Conclusions 155 be minimized. The ray tracing formulae, which form the basis of the inversion procedure D presented here, are easily adapted for application in 3 . However, devising a dependable D method for two-point ray tracing in general 3 media presents a significant challenge. References Aki, K., Christoffersson, A. &: Husebye, E., 1977. Determination of the three-dimensional seismic structure of the lithosphere, J . geophys. Res., 82, 277-296. Alt, J.C., Honnerez, J., Laverne, C. &; Emmermann, R., 1986. Hydrothermal alteration of a 1 km section through the upper oceanic crust, Deep Sea Drilling Project Hole 504B: mineralogy, chemistry, and evolution of seawater-basalt interactions, J . geophys. Res., 91, 10,309-10,335. Andersen, A.H. &; Kak, A.C., 1982. Digital ray tracing iri two-dimensional refractive fields, J . Acoust. Soc. Am., 72, 1593-1606. Anderson, R.N., Honnerez, J., Becker, K., Adamson, A.C., Alt, J.C., Emmermann, R., Kempton, P.D., Kinoshita, H., Laverne, C , Mottl, M.J. & Newmark, R.L., 1982. DSDP Hole 504B, the first reference section over 1 km through Layer 2 of the oceanic crust, Nature, 300, 589-594. Backus, G. & Gilbert, F., 1969. Constructing p-velocity models to fit restricted sets of travel-time data, Bull, seism. Soc. A m . , 59, 1407-1414. Baker, E.T., &; Massoth, G.J., 1985. Heat and particle fluxes from active vent fields on the Juan de Fuca Ridge, Eos Trans. A G U , 66, 930. Bath, M., 1974. Spectral Analysis in Geophysics. Developments in Solid Earth Geophysics, 7, Elsevier, Amsterdam. Becker, K., Von Herzen, R.P., Francis, T.J.G., Anderson, R.N., Honnorez, J., Adamson, A.C., Alt, J.C., Emmermann, R., Kempton, P.D., Kinoshita, H., Laverne, C , Mottl, M.J. & Newmark, R.L., 1982. In situ electrical resistivity and bulk porosity of the oceanic crust Costa Rica Rift, Nature, 300, 594-598. Bessonova, E.N., Fishman, V.M., Ryaboyi, V.Z. & Sitnikova, G.A., (1974). The tau method for the inversion of travel times — I. Deep seismic sounding, Geophys. J.R. astr. Soc, 36, 377-398. Bishop, T.N., Bube, K.P., Cutler, R.T., Langan, R.T., Love, P.L., Resnick, J.R., Shuey, R.T., Spindler, D.A. & Wyld, H.W., 1985. Tomographic determination of velocity and depth in laterally varying media, Geophysics, 50, 903-923. 156 References 157 Bois, P., La Porte, M . , Lavergne, M. & Thomas, G., 1972. Well-to-well seismic measurements, Geophysics, 37, 471-480. Bratt, S.R. & Purdy, G.M., 1984. Structure and variability of oceanic crust on the flanks of the East Pacific Rise between 11° and 13° N, J . geophys. Res., 89, 6111-6125. Bregman, N., 1986. Tomographic inversion of crosshole seismic data, PhD thesis, Dept. of Physics, University of Toronto. Bryan, W.B. & Moore, J.G., 1977. Compositional variations of young basalts in the Mid-Atlantic Ridge rift valley near lat 36°49'N, Geol. Soc. A m . Bull., 88, 556-570. Budiansky. B., 1965. On the elastic moduli of some heterogeneous materials, J . Mech. Phys. Solids, 13, 223-227. Burnett, M.S., Orcutt, J.A. &; Olson, A.H., 1988. A finite element computation of seismic diffraction about the rise axis magma chamber, Geophys. Res. Lett., 15, 1487-1490. Cassell, B.R., 1982. A method for calculating synthetic seismograms in laterally varying media, Geophys. J.R. astr. Soc, 69, 339-354. Censor, Y., 1981. Row-action methods for huge and sparse systems and their applications, SIAM Rev., 23, 444-466. Censor, Y., 1983. Finite series-expansion reconstruction methods, Proc. I E E E , 71, 409-419. Cerveny, V., I., 1985. Gaussian beam synthetic seismograms, J . Geophys., 58, 44-72. Cerveny, V., I., &: Ravindra, 1971. Theory of Seismic Head Waves, University of Toronto, Toronto. Cerveny, V., I., Molotkov, A. & Psencik, I., 1977. Ray Method in Seismology, University of Karlova, Prague, Czechoslovakia. Chander, R., 1975. On tracing seismic rays with specified end points, J . Geophys., 41, 173-177. Chapman, C.H., 1985. Ray theory and its extensions: W K B J and Maslov seismograms, J. Geophys., 58, 27-43. Chapman, O H . k, Drummond, R., 1982. Body-wave seismograms in inhomogenous media using Maslov asymptotic theory, Bull, seism. Soc. A m . , 72, S277-S317. Chapman, C.H., 1987. The Radon transform and seismic tomography, in Seismic Tomography, pp. 25-47, ed. Nolet, G., Reidel, Dordrecht. References 158 Christensen, N.I., 1970. Compressional-wave velocities in basalts from the Juan de Fuca Ridge, J . geophys. Res., 75, 2773-2775. Christensen, N.I., 1979. Compressional wave velocities in rocks at high temperatures and pressures, critical thermal gradients, and crustal low-velocity zones, J . geophys. Res., 84, 6849-6857. Christensen, N.I., 1984. Pore pressure and oceanic crustal seismic structure, Geophys. J.R. astr. Soc, 79, 411-423. Christensen, N.I. & Smewing, J.D., 1981. Geology and seismic structure of the northern section of the Oman ophiolite, J . geophys. Res., 86, 2545-2555. Clayton, R.W. & McMechan, G.A., 1981. Inversion of refraction data by wave field continuation, Geophysics, 46, 860-868. Constable, S.C., Parker, R.L. & Constable, C.G., 1987. Occam's inversion: a practical algorithm for generating smooth models from electromagnetic sounding data, Geophysics, 52, 289-300. Crane, K., Aikman III, -F., Embley, R., Hammond, S., MalahofF, A. and Lupton, J., 1985. The distribution of geothermal fields on the Juan de Fuca Ridge, J . geophys. Res., 90, 727-744. Cudrak, C.F., 1988. Shallow crustal structure of the Endeavour Ridge segment, Juan de Fuca Ridge, from a detailed seismic refraction survey, M . S c thesis, University of British Columbia. Davis, E.E., Lister, C.R.B., Wade, U.S. k Hyndman, R.D., 1980. Detailed heat flow measurements over the Juan de Fuca Ridge system, J . geophys. Res., 85, 299-310. Davis, E.E., Currie, R.G., Riddihough, R.P., Ryan, W.B.F., Kastens, K., Hussong, D.M., Hammond, S.R. & Malahoff, A., 1984. Geological mapping of the northern Juan de Fuca and Explorer Ridges using SeaMARC I, SeaMARC II, and Seabeam acoustic imaging, Eos Trans. A G U , 65, 1110. Delaney, J.R., McDuff, R.E. &; Lupton, J.E., 1984. Hydrothermal fluid temperatures of 400° C on the Endeavour segment, Northern Juan de Fuca Ridge, Eos Trans. A G U , 65, 973. Delaney, J.R., Karsten, J.L. & Hammond, S.R., 1986. Petrology and tectonics of the Juan de Fuca Ridge, Pacific Section, Geol. Assoc. of Canada, Symposium on the Juan de Fuca Ridge System, programme and abstracts, 23-25. Dewey, F. & Kidd, S.F., 1977. Geometry of plate accretion, Geol. Soc. A m . Bull., 88, 960-968. References 159 Dziewonski, A.M. & Gilbert, F., 1976. The effect of small, aspherical perturbations on travel times and a re-examination of the corrections for ellipticity, Geophys. J.R. astr. Soc, 44, 7-17. Ewing, J.I., 1963. Elementary theory of seismic refraction and reflection measurements, in The Sea, 3, pp. 1-19, ed. Hill, M.N., Interscience, New York. Ewing, J.I. & Purdy, G.M., 1982. Upper crustal velocity structure in the ROSE Area of the East Pacific Rise, J . geophys. Res., 87, 8397-8402. Fawcett, J.A. <fe Clayton, R.W., 1984. Tomographic reconstruction of velocity anomalies, Bull, seism. Soc. A m . , 74, 2201-2219. Firbas, P., 1987. Tomography from seismic profiles, in Seismic Tomography, pp. 189-202, ed. Nolet, G., Reidel, Dordrecht. Fischer, K.M. &; Purdy, G.M., 1986. Seismic amplitude modeling and the shallow crustal structure of the East Pacific Rise at 12° N , J . geophys. Res., 91, 14,006-14,014. Fox, P.J., Schreiber, E. & Peterson, J.J., 1973. The geology of the ocean crust: compressional wave velocites of oceanic rocks, J . geophys. Res., 78, 5155-5172. Francheteau, J. & Ballard, R.D., 1983. The East Pacific Rise near 21°N, 13°N and 20° S: inferences for along-strike variability of axial processes of the Mid-Ocean Ridge, Earth Planet. Sci. Lett., 64, 93-116. Franklin, J.M. & Kappel, E.S., 1986. Compositional characteristics and volcanic setting of massive sulphide deposits on the Juan de Fuca and Explorer ridges, Pacific Section, Geol. Assoc. of Canada, Symposium on the Juan de Fuca Ridge System, programme and abstracts, 35-38. Geological Survey of Canada, 1987. Bathymetry, Northern Juan de Fuca Ridge, Marine Geophysical Maps of Western Canada, Map 6-1987, Scale 1:250,000. Gettrust, J.F., Furukawa, K. &; Kempner, W.C., 1982. Variation in young oceanic crust and upper mantle structure, J . geophys. Res., 87, 8435-8445. Grant, S.T., 1973. Rho-rho Loran-C combined with satellite navigation for offshore surveys, International Hydrographic Review, 50, 35-54. Grant, S.T., 1980. BIONAV - The Bedford Institute of Oceanography Integrated Navigation System, Proc. U.S. National Ocean Survey 7th Annual Hydrographic Survey Conference, Washington, pp. 19. Greenbaum, D., 1972. Magmatic processes at ocean ridges: evidence from the Troodos Massif, Cyprus, Nature, 238, 18-21. References 160 Gustavsson, M., Ivansson, S., Moren, P. & Jorgen, P., 1986. Seismic borehole tomography — measurement system and field studies, Proc. I E E E , 74, 339-346. Harding, A.J., Orcutt, J.A., Kappus, M.E., Vera, E.E., Mutter, J.C., Buhl, P., Detrick, R.S. & Brocher, T.M., 1988. The structure of young oceanic crust at 13° N on the East Pacific Rise from expanding spread profiles, submitted to J . geophys. Res.. Hashin, Z., 1970. Theory of composite materials, in Mechanics of Composite Materials, 5th Symposium on Naval Structural Mechanics, eds. Wendt, F.W. and others, Pergamon, New York. Hashin, Z., &: Shtrikman, S., 1963. A variational approach to the theory of the elastic behaviour of multiphase materials, J . Mech. Phys. Solids, 11, 127-140. Hearn, T.M. & Clayton, R.W., 1986. Lateral velocity variations in southern California. I. Results for the upper crust from P waves, Bull, seism. Soc. A m . , 76, 495-509. g Hefner, D.E. & Barrett, D.L., 1979. OBS development at Bedford Institute of Oceanography, Marine Geophys. Res., 4, 227-245. Herron, T.J., 1982. Lava flow layer-East Pacific Rise, Geophys. Res. Lett., 9, 17-20. Hill, R., 1965. A self-consistent mechanics of composite materials, J . Mech. Phys. Solids, 13, 213-222. Hughes, D.S. & Maurette, C , 1957. Variation of elastic wave velocities in basic igneous rocks with pressure and temperature, Geophysics, 22, 23-31. Humphreys, E., Clayton, R.W. &; Hager, B.H., 1984. A tomographic image of mantle structure beneath southern California, Geophys. Res. Lett., 11, 625-627. Hyndman, R.D. & Drury, M.J., 1976. The physical properties of oceanic basement rocks from deep drilling on the Mid-Atlantic Ridge, J . geophys. Res., 81, 4042-4052. Ivansson, S., 1985. A study of methods for tomographic velocity estimation in the presence of low-velocity zones, Geophysics, 50, 969-988. Ivansson, S., 1986. Seismic borehole tomography — theory and computational methods, Proc. I E E E , 74, 328-338. Jacobson, R.S. & Lewis, B.T.R., 1988. Attenuation results from the uppermost young oceanic crust, Eos Trans. A G U , 69, 406. Johnston, D.H., 1981. Attenuation: a state-of-the-art summary, in Seismic Wave Attenuation, Geophysics reprint series, ed. Toksoz, M.N. and Johnston, D.H., Society of Exploration Geophysicists, Tulsa. References 161 Julian, B.R. & Anderson, D.L., 1968. Travel times, apparent velocities and amplitudes of body waves, Bull, seism. Soc. A m . , 58, 339-366. Julian, B.R., 1970. Ray tracing in arbitrarily heterogeneous media, Technical Note 1970-45, Lincoln Laboratory, Massachusetts Institue of Technology. Julian, B.R. & Gubbins, D., 1977. Three-dimensional seismic ray tracing, J . Geophys., 43, 95-113. Kanasewich, E.R. & Chiu, S.K., 1985. Least-squares inversion of spatial seismic refraction data, Bull, seism. Soc. A m . , 75, 865-880. Kappel, E.S. & Ryan, W.B.F., 1986. Volcanic episodicity and a non-steady state rift valley along northeast Pacific spreading centers: evidence from Sea MARC I, J . geophys. Res., 91, 13,925-13940. Karsten, J.L., Hammond, S.R., Davis, E.E. & Currie, R.G., 1986. Detailed geomorphology and neotectonics of the Endeavour Segment, Juan de Fuca Ridge: new results from Seabeam swath mapping, Geol. Soc. A m . Bull., 97, 213-221. Karsten, J.L. & Delaney, J.R., 1987. Temporal and spatial variations in the petrology of lavas from the Endeavour segment, Juan de Fuca Ridge, Eos Trans. A G U , 68, 1540. Kempner, W.C. & Gettrust, J.F., 1982. Ophiolites, synthetic seismograms, and oceanic crustal structure 2. A comparison of synthetic seismograms of the Samail ophiolite, Oman, and the ROSE refraction data from the East Pacific Rise, J . geophys. Res., 87, 8463-8476. Kirkpatrick, R.J., 1979. The physical state of the oceanic crust: results of downhole geophysical logging in the Mid-Atlantic Ridge at 23° N, J . geophys. Res., 84, 178188. Lanczos, O , 1961. Linear Differential Operators, Van Nostrand, London. Langan, R.T., Lerche, I. &; Cutler, R.T., 1985. Tracing of rays through heterogeneous media: an accurate and efficient procedure, Geophysics, 50, 1456-1465. Ludwig, W.J., Nafe, J.E. & Drake, C.L., 1970. Seismic Refraction, in The Sea, 4, Part 1, New Concepts of Sea Floor Evolution, ed. Maxwell, A.E., Wiley- Interscience, New York. Lytle, R.J. &; Dines, K.A., 1980. Iterative ray tracing between boreholes for underground image reconstruction, I E E E Trans. Geosci. Remote Sensing, GE-18, 234-240. References 162 McClain, J.S. & Lewis, B.T.R., 1980. A seismic experiment at the axis of the East Pacific Rise, Marine Geol., 35, 147-169. McClain, J.S., Orcutt, J.A. Sz Burnett, M., 1985. The East Pacific Rise in cross section: a seismic model, J . geophys. Res., 90, 8627-8639. Macdonald, K.C., 1983. A geophysical comparison between fast and slow spreading centers: constraints on magma chamber formation and hydrothermal activity, in Hydrothermal Processes at Seafloor Spreading Centers, NATO Conference Series IV, 12, pp. 27 — 51,eds. Rona, P.A. and others, Plenum, New York. Macdonald, K.C., 1986. The crest of the Mid-Atlantic Ridge: models for crustal generation processes and tectonics, in The Geology of North America, The Western North Atlantic Region, M , pp. 51-68, ed. Vogt, P.R., and Tucholke, B.E., Geological Society of America, Boulder. Macdonald, K., Sempere, J.C. & Fox, P.J., 1984. East Pacific Rise from Siqueiros to Orozco Fracture Zones: along-strike continuity of axial neovolcanic zone and structure and evolution of overlapping spreading centres, J . geophys. Res., 89, 6049-6069. McMechan, G.A., 1983. Seismic tomography in boreholes, Geophys. J.R. astr. Soc, 74, 601-612. Marks, L.W., 1980. Computational topics in ray seismology, PhD thesis, Dept. of Physics, University of Alberta. Marquardt, D.W., 1970. Generalized inverses, ridge regression, biased linear estimation, and nonlinear estimation, Technometrics, 12, 591-612. Menke, W., 1984a. The resolving power of cross-borehole tomography, Geophys. Res. Lett., 11, 105-108. Menke, W., 1984b. Geophysical Data Analysis: Discrete Inverse Theory, Aca- demic, Orlando. MERGE group, 1984. Regional setting and local character of a hydrothermal field/sulfide deposit on the Endeavour segment of the Juan de Fuca Ridge, Eos Trans. A G U , 65, 1111. Morton, J.L., 1984. Oceanic spreading centers: axial magma chambers, thermal structure, and small scale ridge jumps, Ph.D. thesis, Stanford University. Murase, T. & McBirney, A.R., 1973. Properties of some common igneous rocks and their melts at high temperatures, Geol. Soc. A m . Bull., 84, 3563-3592. References 163 Nakanishi, I., 1985. Three-dimensional structure beneath the Hokkaido-Tohoku region as derived from a tomographic inversion of P-arrival times, J . Phys. Earth, 33, 241-256. Nolet, G., 1985. Solving or resolving inadequate and noisy tomographic systems, J . Comp. Phys., 61, 463-482. Normark, W.R., Delaney, J.R., Lichtman, G.S., Morton, J.L., Karsten, J.L. & Johnson, H.P., 1982. Geologic and morphologic constraints on crustal-accretion models for the Southern Juan de Fuca Ridge, Eos Trans. A G U , 63, 1146. Normark, W.R., Morton, J.L., Koski, R.A., Clague, D.A. & Delaney, J.R., 1983. Active hydrothermal vents and sulfide deposits on the southern Juan de Fuca Ridge, Geology, 11, 158-163. O'Connell, R.J. & Budiansky, B., 1977. Viscoelastic properties of fluid-saturated cracked solids, J . geophys. Res., 82, 5719-5735. Oldenburg, D.W., 1983. Funnel functions in linear and nonlinear appraisal, J . geophys. Res., 88, 7387-7398. Orcutt, J.A., 1987. Structure of the earth: oceanic crust and uppermost mantle, Reviews of Geophysics, 25, 1177-1196. Paige, C O & Saunders, M.A., 1982. LSQR: An algorithm for sparse linear equations and sparse least squares, A C M Trans. Math Softw., 8, 43-71. Pallister, J.S. & Hopson, C.A., 1981. Samail ophiolite plutonic suite: field relations, phase variation, cryptic variation and layering, and a model of a spreading ridge magma chamber, J . geophys. Res., 86, 2593-2644. Pavlis, G.L., 1986. Geotomography using refraction fan shots, J . geophys. Res., 91, 6522-6534. Pereyra, V., Lee, W.H.K. & Keller, H.B., 1980. Solving two-point seismic-ray tracing problems in a heterogeneous medium, Bull, seism. Soc. Am., 70, 79-99. Peterson, J.E., Paulsson, B.N.P. & McEvilly, T.V., 1985. Applications of algebraic reconstruction techniques to crosshole seismic data, Geophysics, 50, 1566-1580. Purdy, G.M., 1987. New observations of the shallow seismic structure of young oceanic crust, J . geophys. Res., 92, 9351-9362. Radcliff, R.D., Balanis, C A . & Hill, H.W. Jr., 1984. A Stable geotomography technique for refractive media, I E E E Trans. Geosci. Remote Sensing, GE-22, 698-703. Richards, P.G., 1971. An elasticity theorem for heterogeneous media, with an example References 164 of body wave dispersion in the earth, Geophys. J.R. astr. Soc, 22, 453-472. Richard, P.G. & Menke, W., 1983. The apparent attenuation of a scattering medium, Bull, seism. Soc. A m . , 73, 1005-1021. Riddihough, R., 1984. Recent movements of the Juan de Fuca plate system, J . geophys. Res., 89, 6980-6994. Rohr, K.M.M., Milkereit, B. & Yorath, C.J., 1988. Asymmetric deep crustal structure across the Juan de Fuca Ridge, Geology, 16, 533-537. Rohr, K.M.M. & Milkereit, B., 1988. A shallow crustal reflector across the Juan de Fuca Ridge, submitted to J . geophys. Res.. Salisbury, M.H. & 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., 81, 805-817. Salisbury, M.H., Christensen, N.I., Becker, K. & Moos, D., 1985. The velocity structure of Layer 2 at Deep Sea Drilling Project Site 504 from logging and laboratory experiments, Init. Reps. D S D P , 83, 529-539. Schoenberger, M. & Levin, F.K. 1974. Apparent attenuation due to intrabed multiples, Geophysics, 39, 278-291. Smith, A.M., Goldberg, M. &; Liu, E.S.K., 1980. Numerical ray tracing in media involving continuous and discrete refractive boundaries, Ultrasonic Imaging, 2, 291-301. Spence, G.D., Whittall, K.P., & Clowes, R.M., 1984. Practical synthetic seismograms for laterally varying media calculated by asymptotic ray theory, Bull, seism. Soc. Am., 74, 1209-1223. Spudich, P., &; 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. Stork, C. & Clayton, R.W., 1986. Analysis of the resolution between ambiguous velocity and reflector position for traveltime tomography, S E G Expanded Abstracts, 56th Annual International SEG Meeting. Tabata, S. & Peart, J., 1985. Canadian Data Report of Hydrography and Ocean Sciences, 38, 411-414. Telford, W.M., Geldart, L.P., Sheriff, R.E. & Keys, D A . , 1976. Applied Geophysics, Cambridge, Cambridge. Tivey, M.K., & Delaney, J.R., 1986. Growth of large sulfide structures on the Endeavour Segment of the Juan de Fuca Ridge, Earth Planet. Sci. Lett., 77, 303-317. References 165 Tivey, M.A. & Johnson, H.P., 1987. The central anomaly magnetic high: implications for ocean crust construction and evolution, J . geophys. Res., 92, 12,685-12,694. Todd, T., & Simmons, G., 1972. Effect of pore pressure on the velocity of compressional waves in low-porosity rocks, J . geophys. Res., 77, 3731-3743. Toksoz, M.N., Johnston, D.H. Sz Timur, A., 1979. Attenuation of seismic waves in dry and saturated rocks: 1. Laboratory measurements, Geophysics, 44, 681-690. Trappe, H., 1988. Seismic attenuation in the vicinity of the geothermal anomaly at Urach obtained from near-vertical reflection profiles, Geophys. Prospect., 36, 149-166. Van der Sluis, A. &; van der Vorst, H.A., 1987. Numerical solution of large, sparse linear algebraic systems arising from tomographic problems, in Seismic Tomography, pp. 49-83, ed. Nolet, G., Reidel, Dordrecht. Whittall, K.P. & Clowes, R.M., 1979. A simple, efficient method for the calculation of travel times and ray paths in laterally inhomogeneous media, J . Can. Soc. Exp. Geophys., 15, 21-29. Wiggins, R.A., 1972. The general linear inverse problem: implication of surface waves and free oscillations for earth structure, Ttev. Phys. Space Phys., 10, 251-285. Will, M . , 1976. Calculation of travel-times and ray paths for lateral inhomogeneous media, in Explosion seismology in Central Europe: data and results, pp. 168-177, ed. Giese, P., Prodehl, C. & Stein, A., Springer-Verlag, Berlin. Winkler, K.W. & Nur, A., 1982. Seismic attenuation: effects of pore fluids and frictional sliding, Geophysics, 47, 1-15. Wong, J., Hurley, P. &; West, G.F., 1983. Crosshole seismology and seismic imaging in crystalline rocks, Geophys. Res. Lett., 10, 686-689. Young, G.B. & Braile, L.W., 1976. A computer program for application of Zoeppritz's amplitude equations and Knott's energy equations, Bull, seism. Soc. A m . , 66, 1881-1885. Appendix A Data Reduction A.l OBS Timing Corrections The OBSs utilized in this experiment simultaneously record 4 data channels; internal time code, hydrophone, horizontal and vertical geophones. Absolute timing is established by rating the OBS internal clock with an external reference clock before deployment of the OBS and after recovery. Time drift corrections are made to the OBS data using linear drift rates established from these deployment-recovery clock ratings, assuming a constant drift rate for the intervening time period. Clock drift rates for the individual OBSs are provided in Table A . l . Timing for the individual data channels is provided by the parallel time channel. Timing corrections must be made for advance of the time channel relative to the data channels caused by tape head skew. In the record-playback system (see Section A.3) there are several different tape heads involved, and thus a cumulative skew value must be determined. This is done by comparison of a calibration signal recorded simultaneously on all 4 data channels. The comparison is. made at the output of the playback system so that a cumulative skew value is obtained. Skew corrections applied to the various OBSs are listed in Table A . l . OBSs numbered 1-8 in this thesis correspond to OBS sites 2, 11, 12, 3, 4, 13, 14 and 5, respectively, in the original experiment notes and in the thesis of Cudrak (1988). 166 Appendix A. Data Reduction OBS Drift Rate Initial Advance 1 2 3 4 5 6 7 8 (ms/day) -12.0 -92.1 -101.0 -137.5 -80.7 -37.9 -52.0 -48.2 (ms) -270 33 69 -149 -122 37 13 -254 167 Tape Head Skew Advance (ms) Vertical Hydrophone Horizontal 12 ± 5 50 ± 2 15 ± 1 59 ± 6 -9±5 6±4 -10± 1 19 ± 2 52 ± 2 -185 ± 8 -55 ± 2 106 ± 10 28 ± 2 -36 ± 2 21 ± 2 65 ± 2 60 ± 3 — 113 ± 4 -31 ± 2 151 ± 13 39 ± 2 -24 ± 3 18 ± 2 89 ± 4 Table A . l . OBS timing drift and skew corrections. The initial clock advance (column 3) refers to the OBS clock advance relative to the master clock at the pre-deployment rating, and the tape head skew values are advances of the data channels relative to the time channel. A.2 OBS Repositioning & Final Timing Correction OBS positions were initially determined using SATNAV adjusted p-p Loran-C coordinates with an absolute accuracy of ±180 m (Grant, 1973). The relative positions of the OBSs along the cross-ridge line are greatly improved by using the direct water wave travel times determined for shots passing directly over the OBS site. Relative shot positions were obtained using Loran-C for which random errors of ±30 m have been estimated (Grant, 1973). However, consistency checks applied to the smoothed navigation data for the cross-ridge line suggest that relative shot position errors are ±15 m. The water wave travel time data are very sensitive to position changes of the OBS parallel to the shot line, but are insensitive to position changes perpendicular to the line. Also, static time shifts are difficult to distinguish from timing shifts produced by moving the OBS perpendicular to the line. Because of this lack of discrimination, it is important to remove static time shifts from the travel times before adjusting the OBS positions perpendicular to the shot line. Static time shifts common to all OBSs are caused by uncertainties in the airgun depth and uncertainties in the average water velocity used for travel time calculations, whereas static shifts for individual OBSs are due to uncertainties Appendix A. Data Reduction 168 OBS AE AN 1 2 3 4 5 6 7 8 (m) 24 202 218 48 -50 -81 128 -12 (m) 23 93 -78 22 35 -3 199 34 R M S error (ms) 6.7 9.3 9.3 13.1 7.8 . 8.3 9.8 7.3 static shift (ms) -5 -8 0 -235 -26 -12 0 -38 Table A . 2 . OBS repositioning and static time shifts. The easting and northing position corrections and static time corrections are calculated for each OBS using the direct water wave arrival. The resulting RMS error (after repositioning) for the observed vs. calculated water wave travel times is given in column 4. in each of the following: water depth at deployment sites, OBS clock drift and tape head skew corrections. OBS repositioning is accomplished in two stages. First, the OBS position parallel to the shot line is adjusted (keeping the normal distance from the line constant) and a least squares best-fit static time correction is determined. During this stage, the shot-receiver distance parallel to the line is approximated by the actual shot-receiver distance. The static time shift is applied to the travel time data, and then the OBS position parallel and perpendicular to the shot line is adjusted. The position shifts and time corrections determined for each OBS in this manner are provided in Table A.2, and the algorithm used for repositioning is described in Section A.2.1. Figure A . l shows the direct water wave travel times for OBS 3 compared to the calculated travel times before and after repositioning of the OBS. A constant velocity of 1.485 km/s was used for the water column during OBS repositioning and for ray tracing in Chapter 4. This velocity was determined to be a suitable average velocity for a detailed water column sounding obtained during October (Tabata & Peart, 1985). Appendix A. Data. Reduction 169 CN CN Figure A . l . Travel time comparison for repositioning of O B S 3. Direct water wave travel times (solid line) calculated for an O B S positioned at 0.0 km, are compared with the observed travel times before (circles) and after (dots) repositioning of the O B S . T h e observed travel times are plotted at shot-receiver distances calculated using the corresponding O B S position estimate. T h e O B S was shifted 218 m E and —78 m N from its original position, with no static time shift required. T h e resulting R M S misfit for the travel times is 9.3 ms. \ Appendix A. Data Reduction A.2.1 170 OBS Repositioning Algorithm If an OBS is positioned at (x ,y , z ), then the travel time to the OBS through the water 0 0 0 column from a shot located at (x,y, 0) is t= 2 v - x y + (y- y f + zlf' o =~ , 0 v w (A.l) w where v is the constant water velocity and d is the shot-receiver distance. If the OBS w position is moved by 6x and 6y , then the travel time change is 0 0 St = dv - x )Sx + (y - y )Sy )0 0 0 (A.2) 0 w For a set of M shots, (A.2) can be written as, *-A(£) where 6t' = { (t b 0 eeTve d,i ~t ca i i cu atedti )d{V , w A itl — X{ (.,) — x and A^ = yi — yo- The linear 0 2 system (A.3) is overdetermined. It can be solved in a least squares sense for the position changes 8x and 8y which are then used to update the original estimate of the OBS 0 0 position. This process is repeated until the changes k~x and 8yo are insignificant. To 0 judge how well the OBS positions can be determined with this procedure, it was applied to sjnthetic travel times using an initial OBS position estimate which was 200 m along the shot line and 200 m off the shot line from the true OBS position. Gaussian noise was added to the shot positions (<r=±15 m) and to the travel times (cr=±15 ms), and a water velocity in error by 0.015 km/s was used during repositioning. The resulting OBS position estimate was 18 m from the true position parallel to the line. Based on this test, uncertainties of ±20 m have been assigned to the water wave adjusted positions. Appendix A. Data Reduction A.3 171 D a t a P l a y b a c k and D i g i t i z a t i o n The OBSs used in this experiment record data in analog form on audio cassettes, using the direct record method. Tape speeds of ~ 0.2 mm/s (0.008 in/s) were used to allow approximately 10 days of recording on a single C120 cassette. Prior to digitization, the cassette tapes were rerecorded onto 1/4 in audio tapes (FM) at a tape speed that ensured a data transmission rate that is compatible with the maximum throughput rate of the PDP-11 computer-based digitizer. The rerecord-playback process consisted of playback of OBS cassettes using a TEAC R-61D (direct record) recorder at l | in/s, and recording onto \ in audio tapes using a RACAL STORE-4D (FM) recorder at 30 in/s. The F M tapes were then played back on an HP 3960 Instrumentation recorder at 3 | in/s, with output of the 4 data channels to the digitizer. This sequence of playback steps results in an overall data speed-up factor of ~ 30. Time code from the OBS data tapes was utilized by a hardware-software system to digitize data according to a specified time table. This controlling system started the digitizer prior to a specified shot time, and digitized data within a time window of a specified length following the start time. A 10Hz signal from the time channel was used to control strobing of the digitizer, removing the effects of long-period tape speed fluctuations. The digitized data were subsequently demultiplexed, and reformatted. A.4 Amplitude Compensation Input vs. output voltage curves for 2 of the OBSs (including the playback system) used in this experiment are shown in Figure A.2. As can be seen, the system response is linear for input voltages of < 1.25 V after which the response flattens out, resulting in a soft clipping of amplitudes for higher voltages. It is not certain which part of the record-playback system is responsible for the flat response, although it is most likely due to saturation of either the tape head or the magnetic tape during the original OBS recording. Note that the response curve shapes for the two instruments are very similar, Appendix A. Data Reduction 172 but the maximum voltage levels are somewhat different. The output vs. input response of the OBS-playback system, consisting of a linear and non-linear regime, can be represented analytically using an empirical function of the form, V. WlBVnt _ W \n(l - ~ V /Vma )/A out 2 X W) (Wi + 2 where w = exp{A (V - V )}, x T and B, A, AT, V J and V T max w = exp{A (V; - V )}, out 2 T a n out out w , whereas for V 2 out (A.5) represents the linear response for output voltages < Vr, whereas — ln(l — V /V )/A out T are parameters which determine the shape of the response curve. Vi is weighted average of two functions; BV the response curve as V ut —* V represents the singular nature of max . For V max out « Vr, the weighting function W\ dominates >> Vr, w dominates w\. Thus the weighting in (A.4) produces 2 the appropriate response (linear or logarithmic) according to the output voltage range. Analytic functions determined using (A.4) and (A.5) are shown in Figure A.2. Knowing the response of the instruments, an attempt can be made to compensate the data for the clipping that occurs when the system is overdriven. Due to the singular nature of (A.4) as V out —> V , a voltage V u is specified in addition to the other parameters, so that if V out > V u , then V max c p c out p is set to V u . This stabilizes the compenc p sation procedure in the presence of noise, but we cannot expect to fully recover heavily clipped amplitudes. Calibration curves were only available for two of the OBSs. Due to the similarity of the curve shapes (differing only in their maximum values), analytic response curves of the same form were used for all of the OBSs, with the maximum output voltage levels (V ) determined individually for each OBS using the direct water max wave amplitude at near-receiver distances. At short distances from the OBS, the direct water wave strongly overdrives the recording system and thus provides a good estimate of V . The response curve parameters for each OBS used for amplitude compensation are max Appendix A. Data. Reduction 0 173 1 2 3 4 Input Volts (peak) Figure A . 2 . Input vs. output voltages for the OBS/playback system. The symbols represent electronically calibrated input vs. output voltages for 2 of the OBSs. The curves fit to these points were determined using the analytic formulae (A.4) and (A.5). Appendix A. Data Reduction OBS 1 2 3 4 5 6 7+ 78 174 B Knox A V 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.30 1.00 1.40 1.80 1.57 1.90 1.60 1.90 1.48 0.6 0.8 0.8 1.0 0.8 0.8 0.8 0.8 1.1 1.7 1.2 1.2 1.7 1.7 1.0 1.2 1.2 1.7 T A T 2.0 1.0 1.0 2.0 2.0 2.0 2.0 2.0 2.0 1.329 0.999 1.399 1.799 1.569 1.899 1.599 1.899 1.479 Table A . 3 . The amplitude compensation factors corresponding to (A.4) and (A.5) are shown for each OBS. Note that different parameters were used for positive and negative amplitudes for OBS 7. The units of V , Vr and V u are volts. max c p listed in Table A.3, and the compensated and uncompensated amplitudes are compared in Figure A.3. A.5 OBS-Playback Gain Factors The OBS-playback system gains as determined by electronic calibration are given in Table A.4, and the relative scale factors for each OBS as determined by comparison of calculated and observed direct water wave amplitudes are provided in Table A.5. Comparison of the normalization factors from the two tables for a single OBS, shows that they are in agreement, suggesting that the OBS-bottom coupling factors are similar from site-to-site. All amplitudes are normalized relative to OBS 1 for the amplitude analysis (Section 4.3). The amplitudes shown in amplitude-distance plots thoughout this thesis are related to ground motion in the following way. The 12-bit digitizer has an output range of 0-4096 (or ±2048 once the mean is removed) corresponding to an input voltage range of ± 5 V. Appendix A. Data Reduction 175 o o m 0 10 20 Position (km) 30 0 10 20 30 Position (km) Figure A . 3 . Compensated (dashed line) and uncompensated (solid line) first arrival envelope-amplitudes are compared for each of the OBSs. Envelope-amplitudes were determined using (4.2) with a 150 ms window following the first arrival time. Compensated amplitudes were determined in the same manner except that the seismograms were compensated using (A.4) and (A.5) with the appropriate parameter values from Table A.3, prior to amplitude calculation, OBS locations are identified by the triangles. Appendix A. Data Reduction OBS Amplifier Gain 1 2 3 4 5 6 7 8 (db) 95.4 ± 0.3 105.1 ± 0.2 105.2 ± 0.2 95.4 ± 0 . 3 95.7 ± 0.3 104.3 ± 0.2 103.4 ± 0.2 95.2 ± 0 . 3 176 Record/Playback Gain (db) -6.9 ± 0 . 5 -12.4 ± 0.5 -8.9 ± 0.5 -6.0 ± 0.5 -7.3 ± 0.5 -7.1 ± 0.5 -8.0 ± 0.5 -4.9 ± 0.5 Total Gain (db) 88.5 ± 0.6 92.7 ± 0.5 96.3 ± 0.5 89.4 ± 0.6 88.4 ± 0.6 97.2 ± 0.5 95.4 ± 0.5 90.3 ± 0.6 Normalization Factor 1.00 0.62 ±0.15 0.41 ± 0.18 0.90 ± 0.06 0.99 ± 0.01 0.37 ±0.19 0.45 ±0.18 0.81 ± 0.10 Table A.4. OBS/playback gain factors. The total gain of the OBS/playback system (neglecting the gain of the geophone) is the combined gain of the OBS amplifier and the recorder/playback sytem. These gains were determined electronically using a 10 Hz sinusoid applied to the amplifier inputs. The normalization factor listed in column 5 is that required for gain normalization relative to OBS 1. OBS Scale Factor 1 2 3 4 5 6 7 8 (xlO ) 15 30 35 17 15 37 30 19 3 Normalization Factor 1.00 0.50 0.43 0.88 1.00 0.41 0.50 0.79 Table A . 5 . Normalization factors from water wave matching. The factors required to normalize the amplitudes for OBSs 2-8 relative to OBS 1, as determined from the matching of observed water wave amplitudes to calculated water wave amplitudes, are in agreement with the normalization factors in Table A.4. Appendix A. Data Reduction 177 Thus the output voltage at the geophone is given by: Vgeophone = AoBS — - 12048 G.l where G is the gain of the record-playback system (see Table A.4), and AOBS is the A amplitude output by the digitizer (i.e. the amplitude values that are plotted). Once data from each of the OBSs have been normalized relative to OBS 1, a single gain of G =88.5 db can be used for all data. Then, 3 V GEOPHONE The = AQBS 8 8 5 2 O 8 x i - - / | = AOBS • 9.1 x 10~ volts. 0 velocity sensitivity of the Geospace HS-1 geophone (that is used in the OBSs) at 10 Hz, is 19 volts/m/s. Thus, the vertical component of the ground velocity is obtained from the observed vertical component amplitude by, V g r n m d = AQBS • 9.1 x 8 10~ • =As 0B 9 • 4.79 x 10~ m/s.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Shallow crustal structure beneath the Juan de Fuca...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Shallow crustal structure beneath the Juan de Fuca ridge from 2[sup D] seismic refraction tomography White, Donald John 1989
pdf
Page Metadata
Item Metadata
Title | Shallow crustal structure beneath the Juan de Fuca ridge from 2[sup D] seismic refraction tomography |
Creator |
White, Donald John |
Publisher | University of British Columbia |
Date Issued | 1989 |
Description | The formation of oceanic lithosphere along ocean ridges, and the role that crustal magma chambers play in the accretionary process, continue to be fundamental issues in plate tectonics. To address these issues, a multi-receiver airgun/ocean bottom seismograph refraction line, designed to allow definition of lateral velocity and attenuation variations within the shallow crust, was shot across the Endeavour segment of the Juan de Fuca Ridge near 48° N, 129° W. A tomographic inversion procedure has been developed to invert the first arrival travel times and amplitudes from this profile for 2[sup D] velocity and attenuation structure. The inversion method is suited to multi-source, multi-receiver refraction profiles where source/receiver spacings are denser than for conventional profiles. The travel time-velocity inversion scheme is based on an iterative solution of the linearized problem and allows for determination of continuous velocity variations as well as geometry of subhorizontal interfaces. The iterative procedure requires a good initial estimate of the velocity model. In each iteration, two-point ray tracing is performed to construct a linear system relating travel time residuals to velocity perturbations. A damped least-squares algorithm is used to solve this system for a velocity perturbation which is used to update the current velocity estimate. Once the final velocity structure of the model has been determined, amplitudes can be inverted directly for attenuation. Tests to ascertain resolution of the method reveal horizontal smearing of the solution due to ray geometry, drop-off in resolution with depth, and the effects of source-receiver geometry and velocity structure on resolution. Parameter weighting is important in removing streaking effects (caused by inhomogeneous ray coverage) from the solution. For the purposes of ray tracing, the model is parameterized in terms of constant gradient (velocity and attenuation) cells, which allow use of analytic expressions for kinematic and dynamic ray properties, attenuation and inversion quantities. This parameterization causes scatter in the amplitudes calculated using zero-order asymptotic ray theory, a problem which is remedied by smoothing the velocity models before amplitude calculation. Application of this 2[sup D] tomographic inversion scheme to first arrival travel times and amplitudes for the cross-ridge refraction line produced a 4-layer model for the shallow crust. Layer 1 is 250 — 650 m thick, with v₁ = 2.5 km/s and [Nabla, sub z]v₁ = 0.5 s⁻¹. Layer 2 is ~800 m thick, v₂ = 4.8 km/s and [Nabla, sub z]v₂ — 1.0 s⁻¹. Layer 1 and layer 2 likely represent the sequence of extrusives whereas layer 3 (~800 m thick, v₃=5.8 km/s, [Nabla, sub z]v₃=0.5 s-1) and layer 4 (v₄=6.3 km/s, [Nabla, sub z]v₄=0.3 s⁻¹) are associated with the dike complex and massive gabbro sequence, respectively. An abrupt velocity transition between layer 1 and layer 2 may be a metamorphic front within the pillow basalts. A low velocity-high attenuation anomaly (velocities decreased by < 0.4 km/s and Q ~20-100), which is interpreted as a zone of increased fracture porosity and/or permeability associated with axial hydrothermal circulation, exists beneath the ridge in layer 2 and upper layer 3. Smaller low velocity-attenuative zones in layer 2, located 8 km to either side of the ridge may be loci of off-axis hydrothermal circulation. No evidence is found for the existence of a crustal magma chamber in the depth range of 1.5 — 3.0 km below the seafloor. Tests indicate that a 1 X 1 km zone of partial melt represents the minimum dimension of such a feature that would be clearly detected by this refraction experiment. These results suggest that Endeavour Ridge may be experiencing a period of diminished magma supply with the magma chamber reduced or eliminated by hydrothermal circulation. Asymmetry of the velocity anomalies observed in layer 3 and layer 4 suggest that crustal temperatures are elevated by 125 — 200° C beneath the ridge and to the east relative to temperatures west of the ridge, indicating that a deep crustal or upper mantle melting anomaly may exist east of the ridge. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-10-18 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052620 |
URI | http://hdl.handle.net/2429/29317 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1989_A1 W44.pdf [ 11.23MB ]
- Metadata
- JSON: 831-1.0052620.json
- JSON-LD: 831-1.0052620-ld.json
- RDF/XML (Pretty): 831-1.0052620-rdf.xml
- RDF/JSON: 831-1.0052620-rdf.json
- Turtle: 831-1.0052620-turtle.txt
- N-Triples: 831-1.0052620-rdf-ntriples.txt
- Original Record: 831-1.0052620-source.json
- Full Text
- 831-1.0052620-fulltext.txt
- Citation
- 831-1.0052620.ris
Full Text
Cite
Citation Scheme:
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>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052620/manifest