UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Structure of the upper mantle below soouth-central Saskatchewan from teleseismic travel-time inversion Bank, Carl-Georg 1997

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

Item Metadata


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

Full Text

S t r u c t u r e o f t h e u p p e r m a n t l e b e l o w s o u t h - c e n t r a l S a s k a t c h e w a n f r o m t e l e s e i s m i c t r a v e l - t i m e i n v e r s i o n By Carl-Georg Bank Diplom (Geophysik), Ludwig-Maximilians-Universitat Miinchen, Germany A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF SCIENCE in T H E FACULTY OF GRADUATE STUDIES DEPARTMENT OF EARTH AND OCEAN SCIENCES We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA February 1997 © Carl-Georg Bank, 1997 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Earth and Ocean Sciences The University of British Columbia 129-2219 Main Mall Vancouver, BC, Canada V6T 1Z4 Date: fctt I f f ? 11 A b s t r a c t An array of three-component seismographs was deployed over the central Trans-Hudson Orogen from October 1994 to July 1996 with the objective of characterizing subcrustal lithospheric structure below an area with diamondiferous kimberlite occurance. The two-dimensional array consisted of nine broad-band and eight short-period stations at approx-imately 100 km spacing. Most stations were situated on the Phanerozoic cover of southern Saskatchewan and located above the Glennie domain and adjacent domains as inferred from aeromagnetic anomalies. Additional data came from a feasibility test conducted in the fall of 1991 and from one permanent station run by the Incorporated Research Institutions for Seismology (IRIS). Relative travel-time residuals for 331 teleseismic events were obtained and inverted for slowness variations using a tomography code by VanDecar (1991). The ray coverage allows for depth resolution from 60 km to 400 km below the array. Our results show heterogeneities in mantle slowness above 220 to 270 km that deviate by ±2% from the iasp91 radial Earth model. The maximum slowness anomaly is a quasi-cylindrical vertical body of 120 km diameter having its maximum at (52.8°N, 105.9°W, 100km depth) and extending from 60 km depth to 220 km depth. The minimum slowness anomaly is kidney shaped and has its peak value at (51.8°N, 103.9°W, 170km depth). Cretaceous diamondiferous kimberlites and high concentrations of kimberlitic minerals in glacial tills occur above or near the rims of the low velocity anomalies. These are interpreted as remnants of a mantle plume conduit or Rayleigh-Taylor-like instability in the uppermost mantle. The correlation of the fast anomaly with a long-wavelength gravity low suggests a lithospheric keel, consistent with the presence of the proposed Saskatchewan craton. In the south, a correlation between low mantle velocities and high heat flow corroborates the interpretation of the low velocity I l l anomalies as being primarily thermal in origin. Low levels of heterogeneity below 220 to 270 km depth are interpreted as being representative of the convecting asthenosphere. Our preferred interpretation of the imaged anomalies is the thermomechanical erosion of the lithospheric keel of the Saskatchewan craton during the Cretaceous. Diamondiferous kimberlites are a direct consequence of this process. Material from deeper levels of the mantle now trapped in the lithosphere explains the low velocity anomalies, whereas the remnants of the lithospheric keel give rise to high velocity anomalies. Our analyses indicate a lithospheric thickness of 220 to about 270 km below the Trans-Hudson Orogen. iv Contents A b s t r a c t i i L i s t o f F i g u r e s v i i i A c k n o w l e d g e m e n t i x D e d i c a t i o n x 1 I n t r o d u c t i o n 1 1.1 Tectonic overview 1 1.2 Previous geophysical studies 5 1.2.1 Crustal studies 5 1.2.2 Teleseismic studies 8 1.3 Evolution and signature of velocity anomalies below Archean cratons . . . 10 1.4 Diamonds in the Trans-Hudson Orogen 12 1.5 Objective of this study 14 2 T e l e s e i s m i c D a t a S e t 16 2.1 Sites and Instrumentation 16 2.2 Event coverage 19 2.3 Preprocessing of data 21 CONTENTS v 3 M e t h o d 25 3.1 Relative delay times 26 3.1.1 Preparation of seismic records 26 3.1.2 Multi-channel cross-correlation 27 3.1.3 Minimization of error 29 3.2 Tomographic inversion 30 3.2.1 Rationale 30 3.2.2 VanDecar's inversion method 32 The forward problem 33 Model construction by linear inversion 35 Model appraisal 36 Non-linear inversion 37 4 R e s u l t s 3 8 4.1 Comparison of residual travel times 38 4.2 Tomographic models 40 4.2.1 Models from linear inversion 40 4.2.2 Resolution test 43 4.2.3 Models from non-linear inversion 46 4.2.4 Station and Event corrections 48 4.3 Summary of inversion results 51 CONTENTS vi 5 I n t e r p r e t a t i o n 52 5.1 Possible causes for the seismic anomalies 52 5.2 Comparison of tomographic results to other data sets 54 5.2.1 Surface geology and inferred subsedimentary geology 54 5.2.2 Seismic reflection and refraction and shear-wave splitting 57 5.2.3 Gravity 60 5.2.4 Heat flow 61 5.2.5 Summary of comparison 64 5.3 Hypothesis: Thermomechanical erosion of a lithospheric keel by a mantle plume conduit 64 6 C o n c l u s i o n s 72 R e f e r e n c e s 7 4 A S t a t i o n i n f o r m a t i o n 8 1 B E v e n t l i s t 8 2 C T i m i n g c o r r e c t i o n f o r P R S - 4 a n d W I T S 8 7 D P r a c t i c a l n o t e s o n V a n D e c a r ' s c o d e 8 9 E C o m p a r i s o n o f h e a t d i f f u s i o n f o r v a r i o u s g e o m e t r i e s 0 4 F Z u s a m m e n f a s s u n g - A b s t r a c t i n G e r m a n 9 7 vii List of Figures 1.1 Simplified tectonic map of the study area 2 1.2 PreUminary model of Proterozoic tectonic history 4 2.1 Topographic map of study area showing stations and equipment used . . . 17 2.2 Events used in tomographic inversion 20 2.3 Flowchart showing preprocessing stream 22 2.4 Timing corrections over a period of four deployments at PRS-4 station BIR 23 3.1 Simplified flow chart showing the steps in VanDecar's codes 25 3.2 Seismic records before and after cross-correlation 28 3.3 P-wave velocity structure for iasp91 31 3.4 Rationale of a code for inversion of teleseismic delay times 31 3.5 Taut spline parameterization of model region 34 4.1 Contour plots of four events arriving from different azimuths 39 4.2 Linear inversion with different levels of regularization : . . 41 4.3 Tradeoff curve for linear inversions 42 4.4 Synthetic spike model with alternating positive and negative anomalies . . 44 4.5 Recovered spike model 45 4.6 Model after non-linear inversion 47 4.7 Plot of station corrections 48 4.8 Plot of event time corrections 50 LIST OF FIGURES viii 5.1 Geotherm of Archean lithosphere compared to mantle solidus 53 5.2 Graphical comparison of tomographic results with other data sets 55 5.3 Results of studies relating to upper mantle anisotropy 59 5.4 Results of heat flow studies 63 5.5 Schematic sequence illustrating proposed hypothesis 65 5.6 Thermomechanical erosion of a thick continental lithosphere 67 5.7 Trace of a hypothetical hotspot across North American Plate 69 C . l WWVB time code as recorded by a PRS-4 station 87 C. 2 Extract from a WWVB time trace showing digitizing points 88 D. l Flowchart showing implementation of VanDecar's code 89 E. l Comparison of heat conduction for three geometrical cases 94 E.2 Effect of size of anomaly on heat diffusion 95 ix Acknowledgement This work would not have been possible without the help and support from many indi-viduals, institutions and companies. The people at the Department of Earth and Ocean Sciences (former Department of Geophysics and Astronomy) offered an atmosphere fruitful for learning and discussion. My supervisors Robert (Bob) Ellis and Michael Bostock gave advise, direction and encourage-ment and arranged for my funding. Selection of sites and maintenance of the equipment was in the hands of Bob Meldrum and Scott Dodd who did a great job in the field. John Amor I want to thank for all his programming and other help with the computer. Further codes used in preprocessing were provided by Michael Bostock, Costas Caravas, John Cas-sidy and Andrew Frederiksen. Ron Clowes and Maya Kopylova assisted in literature search and served on my committee. Bruce Buffett helped solving the heat conduction problem and Doug Oldenburg opened the door to inverse theory. My fellow students Denise Long, Stephane Rondenay and Andrew Gorman and the post-docs Dave Baird, Phil Hammer, Holger Mandler and Barry Zelt were open to discussions. John VanDecar kindly provided his inversion code and was available for answering questions. Preprints from William Griffin, Jean-Claude Mareschal, Zoli Hajnal and Suzan van der Lee as well as data from Jacek Majorowicz are kindly acknowledged. Funding for this project was provided by Cameco Corp., Monopros Ltd., Uranerz Exploration and Mining Ltd. with matching funds provided by NSERC. Anglo-American Corporation of SA Ltd. provided the WITS data-loggers. Maps in this thesis were created using GMT by Wessel & Smith (1995) and initial comparison of data was done using Xmap8 by Lees (1995). A big hug to the many unnamed people who shared thoughts and offered me a smile. Fur Katharina 1 Introduction 1 Our understanding of the tectonic architecture and evolution of the Trans-Hudson Oro-gen (THO) in northern Saskatchewan and Manitoba has improved significantly over the past decade. Geoscientists working in diverse fields and largely through the auspices of LlTHOPROBE, Canada's National Geoscience Project, have combined their research results to synthesize an interpretation of the present structure and geologic history of the THO. Some highlights of this cooperative effort are presented here with a focus on the tectonic outline and those results relevant to our study. Our present knowledge of the upper mantle in the area from previous teleseismic studies will also be presented, followed by a summary of the significance and implications of diamondiferous kimberlite discoveries in the region. We conclude this chapter with an overview of craton formation and evolution in the context of the THO and an outline of the objectives of this research. 1.1 Tectonic overview The THO is part of a series of Proterozoic belts that weld together the Archean provinces of the North American craton (Hoffman, 1988). The THO can be traced from South Dakota, where it lies beneath Phanerozoic cover, northwards through northern Saskatchewan and Manitoba and then west across Hudson Bay to the Cape Smith Belt of northern Quebec (see insert in figure 1.1). Its total length spans some 5000 km. It is generally interpreted to be a subduction-related orogenic belt containing juvenile crust sandwiched between the Superior province to the east and the Rae-Hearne provinces to the west (Hoffman, 1990a). From east to west across central Saskatchewan and Manitoba it comprises the following four sections: the Thompson belt, the Reindeer zone, the Wathaman batholith and the Cree Lake Zone, as shown in figure 1.1. The Thompson belt consists of metasediments 1 INTRODUCTION 2 110° 108 ~JC cities, towns v stations used in this study ^ diamondiferous kimberlites F: Fort a la Come S: Sturgeon Lk„ C: Candle Lk. 54' ^ kimberlitic mineral concentration J^ J Williston Basin (Paleozoic and Mesozoic sedimentary rocks) | | Western Canada Sedimentary Basin (Phanerozoic sedimentary rocks) 52 I | Athabasca Basin (Proterozoic sedimentary rocks) | Continental arc plutonic rocks r—i Marginal basin/collisional sedimentary and plutonic rocks | Arc plutons, mixed gneisses Continental margin deposits/ reworked archean basement I | Archean cratons and windows: m S: SahU Granite H: Hunter Bay, I: Iskwatikan Figure 1.1: Simplified tectonic map of the study area Dashed lines on the Phanerozoic cover denote tectonic boundaries as inferred from aeromagnetics (redrawn after Green et al., 1985, Lewry et al., 1994, and Ellis et al., 1996). Inset shows extent of the Trans-Hudson Orogen following Hoffman (1988). and mafic rocks which are interpreted as continental margin deposits that form a narrow foreland belt (Lucas et al, 1993). The Reindeer zone is a 400 km wide collage of arc-related volcanic and plutonic rocks, volcanogenic elastics, arkosic molasse and crustal melt rocks (Lewry & Collerson, 1990). This zone can be divided into a number of lithostratigraphic domains (figure 1.1). Also shown are the localities of three Archean basement windows that structurally underlie the Reindeer zone, the Sahli Granite in the Hanson Lake block and the Iskwatikan window and the Hunter Bay window in the western Glennie domain 1 INTRODUCTION 3 (Lewry & Collerson, 1990). The Wathaman batholith is interpreted as an Andean-type calc-alkaline magmatic arc (Hoffman, 1990a). The Cree Lake zone comprises mainly high-grade orthogneisses and migmatitic paragneisses and is believed to be a broad, reworked hinterland (Lewry & Collerson, 1990). A preliminary model of the Proterozoic tectonic history was presented by Lewry et al. (1996) and schematic tectonic scenarios illustrating this history are published in Ansdell et al. (1995). The model combines data from various geological and geophysical investigations of the past decade. According to this model (figure 1.2), subduction (believed to be towards the southeast) gave rise to oceanic arc magmatism (La Ronge-Lynn Lake arcs and Flin Flon-Snow Lake arcs) within the Manikewan ocean about 1.92 to 1.87 Ga ago. By 1.865 Ga this ocean was about 5000 km wide and similar to the present day Pacific, both in size and tectonic setting. Within this ocean the Glennie-Flin Flon protocontinent (now comprising Flin Flon belt, Hanson Lake block and Glennie domain) developed by amalgamation of arc and oceanic crustal elements. At that time the La Ronge-Lynn Lake arc(s) were still outboard of the passive margin, which followed the southeast side of the Hearne conti-nent. At 1.860 Ga the arcs collided with the continent by subduction of the Hearne plate beneath the incoming arc terrane. The Wathaman batholith was emplaced at 1.855 Ga and is possibly related to a reversal in the sense of subduction after collision of the oceanic arcs. The Kisseynew backarc basin was initiated at that time and was filled during 1.850 to 1.835 Ga with clastic sediments derived from the uplifted Wathaman arc and the uplifted Glennie-Flin Flon protocontinent. During that time span, the tectonic setting resembled a Mediterranean-type microplate-subduction zone. At 1.835 Ga the Saskatchewan craton arrived from the south and collided with the protocontinent, which gave rise to crustal thickening in the Reindeer zone. During the period 1.820 to 1.805 Ga the terminal collision of the Superior craton with the other components resulted in peak metamorphism and an 1 INTRODUCTION 4 Figure 1.2: Preliminary model of Proterozoic tectonic history (a) 1920 - 1870 Ma: subduction related oceanic arcs and oceanic closure (b) ca. 1865 Ma: accretion of Glennie-Flin Flon-Snow Lake protocontinent (c) ca. 1855 Ma: emplacement of Wathaman Batholith (d) 1850 - 1835 Ma: successor arcs and basins (e) 1835 - 1820 Ma: start of main orogeny and destruction of Kisseynew basin (f) 1820 - 1805 Ma: collision of Superior craton (g) ca. 1800 Ma: terminal closure and transpression (h) 1800 - 1780 Ma: late to post-orogenic events redrawn after Lewry et al. (1996) and Ansdell et al. (1995). The sketches are not to scale. 1 INTRODUCTION 5 oroclinal rotation of the Hearne province. By 1.800 Ga the Manikewan ocean was fully closed and transpressive structures evolved in the Thompson belt. Post-orogenic events, including cooling, emplacement of granites and orogen-parallel escape wedging lasted un-til 1.78 Ga, by which time the tectonic history of the Trans-Hudson Orogen was largely complete. During the Phanerozoic, the Williston Basin evolved in Saskatchewan, North Dakota and Montana (Fowler & Nisbet, 1985), creating >3500 m of sediment thickness at 48°N, 104°W (Baird et al, 1996). In the early Cretaceous, diamondiferous kimberlites of the Sturgeon Lake (53.4°N, 106.0°W) and Fort a la Corne (53.3°N, 105.0°W) fields were emplaced (Gent, 1992; Strnad, 1993). In addition, large quantities of kimberlitic diamond-indicator minerals surveyed in glacial and fluvial sediments south-west of Assini-boia prompted Swanson & Gent (1993) to postulate at least two diatreme sources in that vicinity. Glaciation during the ice ages removed the Phanerozoic cover in the North and rafted kimberlite bodies (Strnad, 1993), e.g. those found at Sturgeon Lake. 1.2 Previous geophysical studies Geophysical studies of the THO have focused in large part on the crust. Here we provide a brief summary of results. In addition we report on a number of teleseismic studies which place the Trans-Hudson Orogen in the context of the large-scale upper mantle structure of the Canadian Shield. 1.2.1 Crustal studies Green et al (1985) compiled an aeromagnetic map and inferred the tectonic boundaries beneath the Phanerozoic cover in the light of excellent correlation between magnetic pat-terns and geological domains in the North. Their results are incorporated in figure 1.1. 1 INTRODUCTION 6 Detailed geologic and tectonic mapping performed as part of the LlTHOPROBE project of the exposed THO was accompanied by diverse geophysical studies. These studies have improved understanding of the orogen's three-dimensional structure. Results from import-ant geophysical studies, including reflection and refraction surveys as well as heat flow and electromagnetic surveys, are summarized below. Ending 1981, COCRUST (Consortium for Crustal Reconnaissance Using Seismic Tech-niques) had surveyed an area of 300 km x 600 km in south-east Saskatchewan and south-west Manitoba with eight reflection and refraction lines totalling 2250 km in length. Morel-a-l'Hussier et al. (1987) noted a high-velocity layer at the base of the crust that increases in thickness toward the center of the Williston basin and a thin (41 km) and laterally homogeneous crust below the Williston basin. In addition, they noticed evidence for upper mantle anisotropy from both Moho reflected and refracted phases. These phases display higher velocities (8.4 km/s) parallel to the THO and lower velocities (8.0 km/s) across the THO. LlTHOPROBE'in 1991 and 1994 recorded seismic reflection and in 1993 seismic refraction lines both across and along strike of the THO. Reflection lines 2,3 and 9 are together 800 km long and cut across the entire orogen following closely the edge of the Phanerozoic cover. Analysis of the data by Lucas et al. (1993) and Lewry et al. (1994) revealed both a broad antiform within the crust and a narrow crustal root underlying the Glennie domain. This is interpreted as an exotic microcontinental block (the Saskatchewan craton) that underthrust the THO and crops out in the three basement windows. Hajnal et al. (1997) suggest that the interpreted root below the Glennie domain is actually an older Moho. Interpretation of the refraction data is currently in progress. Initial results suggest large lateral fluctuations in P-wave velocity (8.1 km/s to 8.45 km/s) directly below the Moho (Hajnal et al, 1997). 1 INTRODUCTION 7 A 400 km seismic reflection transect along 48°N was performed by COCORP (Consor-tium for Continental Reflection Profiling) in 1990 (Nelson et al, 1993). Baird et al. (1996) compare results from this survey to the aforementioned LlTHOPROBE survey to the north. Both transects revealed antiformal structures beneath the THO. Apart from these similar-ities a number of differences were noted. In the north, based on geology it was suggested that this core is made up of Archean material and overlain by juvenile rocks. The Moho is defined by prominent subhorizontal reflections. In the south, reworked Archean rocks have been recovered from drill holes; no juvenile crust has been found. In addition, the Moho is characterized by the cessation of dipping reflectors rather than subhorizontal ones. The authors suggested that this difference is a consequence of eclogitization of a crustal root beneath the Williston basin which initiated the subsidence of the Williston basin during the breakup of the supercontinent Rodinia in the Precambrian. Results from magnetotelluric studies on the North American central plains conduct-ivity anomaly were described by Jones & Craven (1990) and Jones et al (1993). They model this highly conductive anomaly that runs along the western La Ronge belt as two distinct westward dipping bodies. Its location in southern Saskatchewan correlates with a gravity high, a magnetic quiet zone and a geothermal anomaly. They conclude that the conductivity is a result of sulphide bodies and suggest that it originated during subduction. Heat flow programs targeted both the exposed part of the THO (e.g. Drury, 1985; Guillou-Frottier et al, 1996) and the platform (Majorowicz et al, 1986; Jessop & Vigrass, 1989). These studies in general found low heat-flow values for the THO as expected for a Proterozoic region. A high value in the Thompson Belt was explained by heat refraction effects (Guillou-Frottier et al, 1996) and high values within the Mesozoic sediments by the redistribution of heat by groundwater circulation (Majorowicz et al, 1986). We note that a maximum in the Weyburn area (SE Saskatchewan) could not be explained satisfactorily 1 INTRODUCTION 8 by the latter mechanism and was found to coincide with hydrocarbon occurences in the Williston Basin. 1.2.2 Teleseismic studies A number of teleseismic studies have attempted to elucidate the nature of the upper mantle below the Canadian Shield, both from a global and regional perspective. A complete review of these studies lies beyond the scope of this section. However, a brief discussion of some fundamental papers and their main results is useful. Over three decades ago Brune & Dorman (1963) inverted surface wave phase velocities and travel times of refracted waves for a layered structure underlying the Canadian Shield. They were perhaps the first to recognize that the underlying mantle was characterized by generally high velocities. Buchbinder & Poupinet (1977) obtained P-wave residuals, defined as observed minus expected arrival times, for Canadian stations after applying station corrections and distance corrections. They found that residuals vary from -0.6 to +1.1 s across Canada with negative residuals corresponding to regions of greater age and lower heat flow. A subsequent study of 5-waves (Wickens & Buchbinder, 1980) showed negative residuals to be centred around station FCC (Churchill, Manitoba) on the Shield, with a value of -1.5 s. These results imply increased lithospheric thickness in this area. Furthermore, Wickens &; Buchbinder (1980) inferred a rapid thinning of the lithosphere to the south. LeFevre & Helmberger (1989) determined upper mantle P-wave velocity structure by forward modelling long-period waveforms of the Moho-refracted phase P n , the direct phase P, and the surface reflected phase PP and matched these to seismograms representing paths largely confined to the Canadian Shield. Their preferred ID model exhibits a high-velocity 1 INTRODUCTION 9 (8.32 km/s) upper-mantle lid to depths of 165 km, a low-velocity zone around 200 km and the transition zone discontinuities marked by a 5% velocity increase at 405 km and a 4% velocity increase at 660 km. They found only slight lateral variations within the shield region. A global S-wave velocity model for the upper 500 km of the Earth was obtained by Zhang & Tanimoto (1993) by using surface wave data. The equal-area block parameteri-zation allows for a horizontal resolution of 1000 km and a depth resolution of 60 to 250 km. Below the Canadian Shield they find a velocity 4% faster than the global average to depths of about 400 km. Grand (1994) inverted travel times of horizontally polarised direct and surface reflected S-phases. The region of interest included North America, South Amer-ica and the Atlantic, and extended to the core-mantle boundary. Horizontal resolution is on the order of 275 km. His results suggest that lateral variations in the upper 400 km correspond to surface tectonic regimes. Deviations from a ID starting model under the Canadian Shield are as follows: +4% at depths of 0-175 km, +3% at 175 to 250 km, and +2% at 250 to 325 km. Higher than average velocities can also be seen below 325 km, although the author notes that resolution is poor in the transition zone. The area of in-terest in our study is located at the western margin of the mantle root underlying the Canadian Shield based on the deviations mentioned above. Polet & Anderson (1995) fur-ther investigated the depth extent of high velocity zones under cratons using the models obtained by Zhang & Tanimoto (1993) and Grand (1994). They observe similar anomalies beneath the Early Proterozoic and Archean parts of the Canadian Shield, but not below the Middle Proterozoic parts. They offer two hypotheses to explain these high velocity zones beneath the older regions, one involving thick roots and another involving thinner roots. In the former hypothesis, thick roots developed as a result of either extraction of komatiites, or shallow subduction combined with underplating of the lithosphere. In the 1 INTRODUCTION 10 latter hypothesis, a thinner lithospheric root induces downwelling of cold material. Other tomographic studies, either on a global scale (Pulliam et al., 1993) or confined to North America (van der Lee & Nolet, 1996), support Grand's model revealing the central THO to overly the southwest margin of an extensive high velocity cratonic root. During the second half of 1991 a feasibility study was conducted in the area of interest in the present work to test equipment and to assess the potential of teleseismic techniques to reveal lateral variations in lithospheric structure (Ellis & Hajnal, 1992; Ellis et al., 1996). Results from this investigation include variations in travel time residuals of several 100 ms across their array and a significant variation in residual patterns for events arriving from different azimuths. The data from this feasibility study are incorporated in the present work. Silver & Kaneshima (1993) measured shear-wave splitting parameters for five earth-quakes recorded at 25 stations along a 1500 km long SW-NE trending line that crossed the THO from the Wyoming province in the SW to the Superior province in the E. They noted that the fast polarization direction is generally parallel to exposed geological features outside the THO but changes abruptly to an E-W azimuth within the THO. They suggest fossil anisotropy to be responsible for the observed splitting of shear waves. 1.3 Evolution and signature of velocity anomalies below Archean cratons The Canadian Shield was one of the first regions where a high velocity mantle root was noted. High velocity roots have since been confirmed beneath other Archean shields, like the Baltic Shield and the West African Shield (e.g. Polet & Anderson, 1995). At present we do not clearly understand the cause of these roots. 1 INTRODUCTION 11 When considering mantle roots, it is important to clarify the terminology used. We ad-opt the one suggested by Jordan (1978) which distinguishes three regions: the lithosphere, the tectosphere, and the asthenosphere. The lithosphere is (by its original definition) the mechanical boundary layer which is capable of supporting large deviatoric stresses. The tectosphere is defined to be the region that translates coherently during plate motion; it is equivalent to the thermal boundary layer where heat flow occurs by conduction, and the lithosphere forms part of it. The 1200°C isotherm marks the boundary between the tectosphere and the underlying asthenosphere. In the asthenosphere, heat is transported by convection and mass movements account for isostatic adjustments. In addition to a mechanical division into rigid and viscous, or a thermal division into conductive and con-vective, a distinction can also be made on geochemical grounds. A convecting homogeneous mantle of MORB composition is distinct from the chemical boundary layer that consists of a depleted upper mantle and a felsic crust. Although the general concept of a tectosphere is now widely accepted, debate still exists on what may have caused thick roots to accumulate preferentially beneath Archean cratons. In the literature one can find three competing hypotheses: (i) depleted mantle, (ii) buoyant underplating and (iii) transient root. Jordan (1978) when introducing his concept suggested mantle roots to be a residual, depleted mantle. Extraction of basaltic melt from garnet lherzolite would leave a 0.9 to 1.5 % less dense residual than the original source. This density contrast would be sufficient to compensate for the negative thermal buoyancy. Helmstaedt & Schulze (1989) showed that buoyant underplating by subducted slabs could explain inclusions of Archean eclogites in South African kimberlites. Support for this mechanism was offered by Hoffman (1990b) on the basis of geological observations. Comparing crustal ages for the shield and the platform (the latter inferred from isotope data and drill cores) he concluded that the mantle underlying the Canadian Shield must 1 INTRODUCTION 12 have formed during the Archean, most likely by buoyant subduction. Abbott (1991) also favours buoyant subduction using geophysical arguments like thermal evolution. Polet & Anderson (1995) offered yet another explanation. A small "physical" root may well induce cold downwellings, causing a "transient" root to develop. These three hypotheses have different implications on the nature of high velocity roots. If depletion formed the mantle root, the anomalies would be due in part to chemical contrasts. If buoyant subduction caused the root, the anomalies would arise from this chemical difference between subducted slabs and upper mantle material. If a transient root is evoked as cause, a cooler temperature would suffice in describing the anomalies. 1.4 Diamonds in the Trans-Hudson Orogen Diamonds within the THO may have been discovered as early as 1948 (Lehnert-Thiel et al, 1992), but it was not until summer 1988 that vigorous aeromagnetic exploration was initiated following the discovery of a diamondiferous kimberlite cluster at Sturgeon Lake, 30 km NW of Prince Albert (Gent, 1992) (see fig 1.1). A second kimberlite field was found near Fort a la Come, 80km E of Prince Albert, where diamond grades are estimated to be as large as 1 carat per 5 ton of kimberlite ore (Kjarsgaard, 1996). More recently, exploration has spread to Candle Lake, 70 km NNE of Prince Albert due to the discovery of other kimberlite pipes (Kjarsgaard, 1996) which by now have been confirmed to be also diamondiferous (Bruce Kjarsgaard, 1996, personal communication). In addition to these areas which are being explored Swanson & Gent (1993) point to an area to the SW of Assiniboia where anomalously high quantities of kimberlitic minerals were surveyed in glacial tills. They postulate at least two diatreme sources in the vicinity of that location. The presence of diamondiferous kimberlites beneath the central THO has important 1 INTRODUCTION 13 implications for the structure and physical state of the upper mantle in the region. Diamond is the high-pressure phase of carbon with a stability field which is denned in pressure-temperature-space by a line extending from (50kbar, 1000°C) to (60kbar, 1600°C) (see Kirkley et al., 1991). Within the Earth, these pressures correspond to depths of 150 and 200 km below the surface; which are minimum depths for diamond genesis. The presence of diamonds at the Earth's surface indicates that they have been brought up sufficiently fast (a few hours ascent time) such that the phase transition to graphite does not occur. The transport media are kimberlites and lamproites, which are hybrid1, volatile rich, mafic or ultramafic igneous rocks, thought to originate between 150 to 350 km depth. Both occur in carrot or champagne-glass-shaped pipes and although widespread on all continents it was recognized by Clifford (1966) that diamondiferous kimberlites are mostly confined to the platforms of old cratons. Their emplacement appears to have occured in global pulses throughout the Earth's history (Haggerty, 1994). The age of an individual diamond crystal can be determined by dating cogenetic mineral inclusions, and diamond geochronology suggests that diamonds have been forming continuously through the Earth's history. Further it is known that diamonds are usually much older than the kimberlites that host them. As a consequence, a stable storage mechanism for diamonds must exist below old cratons. The thick cool cratonic keels, isolated from convective mixing within the mantle, are within the diamond stability field and thus are thought to act as effective repositories for diamonds. Occasionally they are intruded by kimberlitic magmas that bring the diamonds to the surface. A decade ago, the THO was believed to consist mainly of juvenile Proterozoic crust, and thus the discovery of diamonds within the central part of the THO was not expected 1 The term hybrid describes rocks containing xenoliths and xenocrysts that did not crystallize from the same magma; volatiles are mainly CO2 for kimberlites and F for lamproites 1 INTRODUCTION 14 (e.g. Helmstaedt & Gurney, 1995). As mentioned in section 1.2, seismic reflection studies (Lucas et al, 1993) have revealed a root zone beneath the Glennie domain which was interpreted as an "exotic" Archean block. Comparison of this result with those from a southern line and from drillcores in the South led Baird et al. (1996) to suggest that this crustal fragment extends as far south as North Dakota. Subsequently it was termed the Saskatchewan craton or abbreviated the Sask craton (Lewry et al., 1996; Chiarenzelli et al, 1996). In light of this new tectonic interpretation, diamonds may very well be associated with a relic cratonic keel. However, two questions remain unanswered: (i) is cool lithosphere present at ca. 200 km depth that allowed for storage of diamonds between the Archean and the Cretaceous? (ii) if kimberlites are, in fact, associated with mantle plumes as proposed by Haggerty (1994), can we image in the upper mantle those plumes that triggered kimberlite emplacement? 1.5 Objective of this study Seismology is one of the few tools that can be used to interrogate the structure of the upper mantle and to infer its composition and evolution. This study will add to the understanding of the tectonic history of the THO, the structure of the upper mantle below the area, velocity anomalies beneath Archean cratons and diamond finds in the study area. Weaving these four areas of research together in the context of the THO leads to the following working hypothesis: We adopt the current view arising from geological history and reflection seismology that the Glennie domain overlies the Archean Saskatchewan craton. If this is the case, then the Saskatchewan craton, like other Archean blocks (including the much larger, adjacent Superior province) may be underlain by a lithospheric root a few 100 km thick. Support-1 INTRODUCTION 15 ing evidence for this hypothesis is the presence of diamondiferous kimberlites within the Glennie domain. Geophysical studies discussed above have presented either a detailed pic-ture of the crust (e.g. Hajnal et al, 1997), or provided low resolution, incomplete views of the upper mantle below the central THO (e.g. Grand, 1994). The primary objective of this study is thus to establish the presence, and extent of a mantle root beneath the Saskatchewan craton; and more generally to understand what role the upper mantle has played in the tectonic evolution of the region. Before testing our hypothesis we will discuss the method we have used and show the results we have obtained. 16 2 Teleseismic Data Set In a teleseismic experiment, natural earthquakes or large man-made explosions at epicentral distances greater than 30° serve as sources of elastic waves. These waves when travelling through the Earth to a station encode a variety of information pertinent to the source itself, the Earth structure along the path, and the equipment used for recording. In seismology we analyse these recordings to obtain an understanding of the Earth's interior. An important parameter, and the information used in this thesis, is the time of arrival of seismic energy. Historically, travel times have given rise to our picture of the Earth. In this chapter we describe the acquisition of a data set collected to elucidate structure below south-central Saskatchewan. 2.1 Sites and Instrumentation The data used in this experiment were obtained from October 1994 to August 1996 using a temporary array of 17 portable seismographs (figure 2.1). To this data set, records from IRIS (Incorporated Research Institutions for Seismology) station FFC (Flin Flon, Manitoba) and data obtained during the feasibility study conducted in 1991-1992 (Ellis et al., 1996) were added. The array was designed with the objective of obtaining a station spacing of approxim-ately 100 km, based on spatial variations of observed travel-time residuals (Ellis & Hajnal, 1992). In some instances this was not possible due to the accessibility of suitable sites (figure 2.1). In late 1995 some stations were moved to improve station coverage. Appendix A lists information concerning station locations and recording times. Note that all stations with the exception of DLK and FFC are located on Phanerozoic cover. Sites were selected during the fall of 1994 and most instruments were deployed in cabins or abandoned farm 2 TELESEISMIC DATA SET 17 1 1 0 ° 100* • MCR/NARS station • WITS station • PRS-4 station • PRS-4 station (1991) o broad-band seismometer • short period seismometer A IRIS broad-band station Figure 2.1: Topographic map of study area showing stations and equipment used houses. Main criteria for the selection of sites were their distance from busy roads or rail-ways and access to electrical power to run the data loggers. Care was taken to set up the seismometers in concrete basements, or - where this was not possible - to bury them (see appendix A for details regarding each station). A variety of equipment was used for the experiment. Four types of recorders were deployed (see also figure 2.1): - N A R S CSD20 data loggers are manufactured by Interay B. V. in The Netherlands. TELESEISMIC DATA SET 18 Five were deployed in summer 1995. Recording continuously at 20 Hz, they write SEED (Standard for the Exchange of Earthquake Data) files to a Digital Audio Tape (DAT). The storage capacity of a tape allows for ca. 190 days of continuous recording at a single station. Internal batteries can power the system for about one day if external power fails. A GPS Global Positioning System (GPS) is used to maintain accurate timing. - W I T S data recorders, developed by the University of Witswatersrand, also record continuously at 20 Hz (50 Hz at USK) sampling rate and provide recording capacity to about 1 month (2 weeks at USK) worth of data. The east channel switched every 4 hours for 90 seconds to WWVB to allow for monitoring of internal clock drift (see section 2.3). - Sc iNTREX 4 channel Portable Recording Seismographs (PRS-4) were developed for applications in both seismicity and crustal scale refraction experiments. Eight of these instruments were used in this study in trigger mode. Data were sampled at 50 Hz with 30 s pre- and 60 s post-event time. Up to 128 events can be stored in R A M and data can be held for up to three days in the event of external power loss. During each trigger WWVB is recorded on the fourth channel to perform timing corrections later. Programming and data retrieval in the field is performed using LithoSEIS, a software package specifically developed for this instrument (Asudeh et al, 1993). - As an interim measure prior to the delivery of the NARS instruments, five T E L E D Y N E G E O T E C H Microcorder MCR-600 ( M C R ) were deployed.These instruments were run in trigger mode at 20 Hz sampling rate recording 60 sec of data per trigger. The field tapes can hold up to 60 events at these parameter settings. Data from these recorders are not included in this study since accurate timing is not retained. They 2 TELESEISMIC DATA SET 19 axe mentioned here for completeness, and because the data were converted and stored with the other data. This data is suitable for waveform studies, where absolute timing is unnecessary (e.g. receiver function analysis). In addition to different recording systems, three kinds of seismometers were employed: - GiJRALP CMG-3, bandpass 0.05 - 5 Hz - M A R K P R O D U C T S L4-C, bandpass from 1-5 Hz, - W l L L M O R E Mark II, 0.5 - 5 Hz The high cut is provided by filters in the data loggers. The latter two types of seis-mometer are short-period and useful primarily for timing whereas the CMG-3 are broad-band and allow also waveform analysis. 2.2 Event coverage Teleseismic studies rely on large sources of seismic energy, namely earthquakes and, to a lesser extent, nuclear explosions. Rather than being uniformly distributed over the Earth's surface and interior, most earthquakes occur along plate boundaries and all occur at depths shallower than 700 km depth. This is one of the key differences which distinguishes seismic tomography from medical tomography. In the latter case, the location of sources outside the human body is completely within the control of the investigator, which leads to a much higher resolution. A non-uniform earthquake distribution will cause uneven ray coverage within the region to be mapped and hence lead to lower resolution. In the present study, coverage is quite good as illustrated in figure 2.2, where a total of 342 teleseismic events is shown (a complete 2 TELESEISMIC DATA SET 20 magnitude O 7.0 Mb O 6.5 Mb O 6.0 Mb o 5.5 Mb c 5.0 Mb . 4.5 Mb Figure 2.2: Events used in tomographic inversion. A total of 331 events in the magnitude range 4.6 to 6.7 were used in the present study. Note that red colour indicates nuclear explosions. Numbers on circles denote epicentral distance to array center in degrees. list of the events is provided in appendix B). Four of these events were explosions from nuclear tests conducted by France and China during the operation of our array. All but one event occured at epicentral distances between 30° and 105° from the array. The lower limit is set to avoid regional phases that are strongly affected by lithospheric structures and triplications at the upper mantle discontinuities. The upper limit is a consequence of the Earth's structure. Direct P- or S-wave are not recorded between 105° and 142° because waves are refracted as P-waves into the core, and thus produce a shadow zone. (For the 2 TELESEISMIC DATA SET 21 event in the Sumbawa, Indonesia, region on Oct. 19, 1991 records from the PKiKP phase were used.) The events shown in figure 2.2 occured between August and December 1991 and between November 1994 and April 1996 and were recorded by at least four stations. The earthquakes are located along the major subduction zones of the Pacific ocean (South and Central America, Aleutian Islands, Kurile Islands, Japan, Mariana Trench and Tonga/Kermadec Trench). A smaller number of earthquakes from the Mid-Atlantic Ridge and the Eastern Mediterranean as well as nuclear test explosions in the Tuamuto Islands Region of the South Pacific and in China improve the event coverage. 2.3 Preprocessing of data A number of steps are necessary before processing of data can start. These include con-version to a common format, timing corrections and selection of seismograms with high signal-to-noise ratio. The key steps in preprocessing the different instrument types are summarized in the flowchart of figure 2.3. The SAC (Seismic Analysis Code of the University of California) format was chosen as the common format for all records. Estimation of relative travel times to within ±10 ms is essential for tomographic imag-ing. Over a few days or weeks internal quartz crystal clocks of most instruments drift too much to be reliable, and other methods are required to produce this accuracy. Accurate timing does not necessarily imply absolute timing but it does require that times at a particular station are known to within a constant offset. Correction of instrument delay times, for example, is not required in this study. Any clock variation with time, however, must be corrected before data processing can be undertaken. Any contribution 2 TELESEISMIC DATA SET 22 IRIS jequest^ NARS WITS DATtap?) ( h ^ d r i v e ) PRS4 MCR J Cfidid laptop^) play to system I 1 I I convert to SAC I I I I change filename to standard correct headers if necessary I I 1 J extract hourlong windows select triggers to view I 1 correct times with WWVB I I quality control and archive Figure 2.3: Flowchart showing preprocessing stream to the delay times that is constant for a particular station will be absorbed into the station correction term in the inversion process. This in turn implies that stations where the instrument setup changed must be treated separately. The PRS-4 instruments record WWVB radio signal on the fourth channel at every trigger. This record is used to correct for instrument drift - an example of how this is accomplished is outlined in appendix C. Figure 2.4 shows an example of the drift as a function of Julian day at PRS-4 site BIR. Every trigger is denoted by a point and it is apparent that the offset at the beginning of each deployment is not zero as one would expect but f«-90ms. This offset reflects an instrument delay and is thus absorbed by the 2 TELESEISMIC DATA SET 23 —J 1 1 : 1— 150 200 250 300 Julian day In 19S5 Figure 2.4: Timing corrections over a period of four deployments at PRS-4 station BIR station correction. WITS dataloggers record continuously on 3 channels. In order to monitor clock drift, the systems were programmed to switch the east channel from the seismometer to a WWVB receiver every seven hours and record the time signal for 90 seconds. (A seven hour interval was chosen to avoid recording at the same time of day since WWVB-reception may exhibit diurnal variations.) A correction procedure similar to that described for the PRS-4 was used to correct for instrument drift. However, since time traces are seldom recorded with an event, corrections have to be interpolated and hence, errors in timing for the WITS stations are higher than for PRS-4 stations, or about ±10 ms. Due to technical problems, WITS units at stations LAN and TSD recorded no time signal. Data from these seismographs were thus not incorporated into this study. NARS recorders adjust their internal clock at intervals of 50 s by recording GPS signal. For time intervals when GPS lock was lost, clock drift was found to be comparable to that of PRS-4 recorders. Occasionally major inconsistencies appeared in the time stamps of adjacent records, a problem which is probably software related. Events recorded during times when the NARS header indicated no GPS lock were not considered for the inversion. 2 TELESEISMIC DATA SET 24 However, this problem affected very few seismograms. Seismograms to archive were selected primarily by visual inspection with the main criterion being a good signal-to-noise ratio. Extracted time series from the continuously recording stations underwent a similar quality control before archival. A constraint put forward by the inversion process is that at least four stations must have recorded an event. This constraint reduced the data set from a total of 918 recorded events to 331 usable events. The reason for this constraint will be outlined in the next chapter. 25 3 Method The objective of this chapter is to describe (i) the extraction of accurate delay times from the seismic records, and (ii) the methodology employed to invert these delay times for a velocity model of the underlying lithosphere and upper mantle. Both the codes to obtain relative delay times and to perform the tomographic inversion were developed by VanDecar (1991). A simplified flow chart including both codes is shown in figure 3.1 with the steps described in the following sections. Practical notes on the code seismograms multi-channel cross-correlation (reference model • i ~ r (relative delay-times 1 O rays )+"{new ray paths i t Frechet derivatives raytracing inversion new model Figure 3.1: Simplified flow chart showing the steps in VanDecar's codes (Van-Decar, 1991). Rectangular boxes denote software, oval boxes denote files. Note that the "rays" file associates each delay time with a ray path through the reference model. A non-linear inversion (dashed arrows) is done by iteratively per-forming "inversion" and "raytracing" until the model norm no longer changes significantly at which step the process is interrupted. 3 METHOD are provided in appendix D. 26 3.1 Relative delay times Accurate estimation of relative delay times is essential for tomographic inversion. VanDe-car (1991) and VanDecar &; Crosson (1990) developed an efficient code to extract relative arrival times and estimates of timing uncertainties for teleseismic traces recorded by a regional array. The semi-automated program, MCCC (multi-channel cross-correlation), makes use of the coherency of teleseismic waveforms between stations of an array. It determines relative delay times for all possible pairs of traces and minimizes the inconsist-encies between these by a least squares method. The main advantage of a cross-correlation method over manual picking is that windowed waveforms comprising the first arrival are quantitatively compared to one another, thus mitigating biases due to low signal-to-noise ratio and emergent onsets. The following is summarized following VanDecar (1991) and VanDecar & Crosson (1990), and is divided into three parts: (i) preparation of seismic records, (ii) multi-channel cross-correlation, and (iii) minimization of error. 3 . 1 . 1 P r e p a r a t i o n o f s e i s m i c r e c o r d s The first criterion in the selection of records is imposed by the inversion code (see discus-sion in section 3.2), namely that a particular event be recorded by at least four stations. Vertical component recordings are preferred for traveltime determination since P-wave signal-to-noise ratio is generally highest and there is less interference from crustal phases. In cases where the vertical component seismogram shows low signal-to-noise ratio or was not recorded the horizontal components can be used since the P-waveforms on vertical and radial components are generally quite similar though some biassing may be expected from 3 METHOD 27 crustal phases. Station LBL, in particular, recorded no vertical component between July and November 1995 and thus radial waveforms were used. The different response of seismometer types to the ground displacement or velocity results in significantly different appearances between broad-band and short-period wave-forms (see figure 3.2) of a seismic phase. A seismic phase is denned by a particular raypath between source and receiver and the particle motion (longitudinal for P and transverse for S) along this raypath or segments of this raypath. This variation in appearance is primarily the result of significantly lower frequency content in the broad-band recordings. A narrow bandpass filter (generally 0.4 to 1.6 Hz) was therefore applied to all records. The phase delay can be considered constant over that narrow frequency band. In addition to filtering, the short period records were inverted to obtain similar waveforms that could be utilized for automated cross-correlation. The effect of this preparation is illustrated in figure 3.2. Prior to cross-correlation the user must pick the approximate onset time of the phase, if, for each trace i. This prepick is necessary to avoid cycle-skipping of the filtered traces. In all cases the first distinct feature (i.e. amplitude maximum or minimum) was chosen as the pick. 3.1.2 Multi-channel cross-correlation A relative delay time between two traces can be obtained by determining the maximum of their cross-correlation function. The cross-correlation function between two time functions Xi(t) and Xj(t) can be expressed as (e.g. Dobrin, 1976) CO #y(r)= J Xi{T-t)'Xi{T)'dT (3.1) —oo where r is the lag time by which the second trace is shifted relative to the first. In our case, data have been digitized at sampling frequencies of 20 or 50 Hz. MCCC uses a METHOD 28 9 5 / 0 9 / 2 3 266 22:31:58.3 10.529S 078.697W 073 5.9Mb NEAR COAST OF PERU SOU NEV 625 630 635 640 645 650 655 660 665 670 675 680 time in seconds after event time 10 12 relative time in seconds Figure 3.2: Seismic records before and after cross-correlation. Top: Original seismograms, vertical bar denotes iasp91 expected P-wave arrival. Note the difference in waveform between the broad-band records at ERW, USK, FFC and the other records. Bottom: Same traces aligned as to maximise their cross-correlation functions. Bold bars denote the expected arrival times, as above. Note that the time scale has been expanded. The window length for cross-correlation in this case was set to 4.5 s starting 0.5 s before the prepicked onset of the phase denoted by dashed lines. A bandpass filter of 0.4 to 1.6 Hz was applied, the short-period records were polarity reversed and all traces were resampled to 100 Hz before cross-correlation took place. 3 METHOD 29 weighted average to interpolate all traces to a common sampling interval (see Wiggins, 1976). For our study the time interval was chosen to be At =0.01 s . In addition, a finite window length T, between 3 and 5 s chosen by the user, was used for cross-correlation (the optimum window length is determined by the source function which generally depends on the magnitude of the event). Thus, s = ^ sample points fall within the correlation window. The window starts at a time interval to before the visual pick if. to has been set to 0.5 s . Equation(3.1) is rewritten for discretized time traces as = - E as*(*F + *o + kSt + T) • + t0 + kSt) (3.2) 3 k=i with the normalization 1/s performed so as to easily compare the results of cross-correlating events with differing window lengths. The maximum of the cross-correlation function (3.2) is determined by an efficient searching algorithm. It cascades over two levels of coarse estimates to a final pass at full resolution which shifts by only five samples. Because the user has prepicked an arrival time, a maximum first shift of 0.35 s is allowed. The three level approach reduces computation time considerably. The delay time between traces z and j is obtained by Aty = *? " £ - <3-3) with rt™ax being the lag time that maximises the value of the cross-correlation function. 3 . 1 . 3 M i n i m i z a t i o n o f e r r o r Delay times as defined in (3.3) for any selection of three traces will, in general, not be consistent with one another, i . e. the sum of delay times for two pairs will not equal the 3 METHOD 30 delay time for the third pair. A least squares method is implemented in MCCC to minimize these inconsistencies. A total of n traces yields pairs of traces to cross-correlate and a similar number of linear equations: U-t^AUj • (3.4) To this system a zero mean for the A i y is imposed to render the system non-singular. Least squares inversion of this system results in the residual time of a particular trace as the arithmetic mean of all delay times associated with that trace U=-J2 A t a • (3-5) MCCC also calculates an uncertainty estimate which will be used to weight delay times in the inversion. 3 . 2 T o m o g r a p h i c i n v e r s i o n In this section we explain the inversion of the delay times determined in 3.1 for upper mantle velocity structure. 3.2.1 Rationale Non-zero travel-time residuals indicate the presence of topography and, more important, deviations of the subsurface velocity structure from an assumed earth model. The radial iasp91 (Kennett & Engdahl, 1991) model was chosen as reference earth model for this study. Figure 3.3 shows the P-wave velocity according to iasp91 from the surface to the 3 METHOD 31 Figure 3.3: P-wave velocity structure for iasp91 center of the earth as well as an enlargement for the upper 670 km, which is the region of interest in this study. Figure 3.4 shows the rationale of teleseismic inversion. The total travel-time delay is an accumulation of delays along the raypath traversing regions of variable velocities. model: - deviation from iasp91 - minimum structure event corrections station corrections 3|C event • station Figure 3.4: Rationale of a code for inversion of teleseismic delay times 3 METHOD 32 In seismic tomography one seeks to exploit redundancy of information provided by many rays crossing at many points. An array of stations is used in teleseismic tomography to delineate the structure of the lithosphere below where rays cross. Rays do not cross in the crust; illumination is poor and any slowness variation below a station will be pulled into the station correction term along with instrumental delay times. Baseline shifts (like the choice of pick time in cross-correlation), lower mantle structure and near source structure (e.g. descending slab in subduction zones) are largely accommodated by event timing corrections and event relocations. Because the inversion cannot produce a unique model, the philosophy behind the inversion is to obtain a minimum structure model, i.e. to retrieve only structure necessitated by the data. As much of the delay times as possible will be assigned to station and event corrections as mentioned before. A minimum of four records per event is necessary since delay times between three or two stations can be fully explained by earthquake relocations and thus do not impose any constraint on the slowness (reciprocal of velocity) model. 3.2.2 VanDecar's inversion method The code used for inverting the travel times was developed by VanDecar (1991) for an elastically isotropic earth. Bostock and VanDecar (1995) and VanDecar et al. (1995) explain later modifications of the source code. The following is a summary of these works. For a detailed description the reader is referred to VanDecar (1991). To understand the inversion, three questions have to be posed (Oldenburg, 1994): - How is the forward problem solved? - How are models being constructed? - What do these models have in common? 3 METHOD 33 These questions are answered in sections to below. Section explains how the non-linearity of the problem is taken into account. The forward problem Seismic tomography is based on ray theory which in essence states that energy travelling from source to receiver has travelled along narrow paths within the earth. Using the ray theoretical approximation we can relate the total travel time t required for the wave to pass from source to receiver with the velocity distribution at each point along the ray path where s0(r) is the slowness value at point r. The value of the slowness is the reciprocal of the seismic wave velocity. The usage of the slowness is preferred in this context as this provides a linear relationship between the physical Earth and travel time. Fermat's principle states that the travel-time between two points is stationary, i . e. the ray along which the wave travels from one point to the other corresponds to a stationary time path, which in many cases is a minimum. As a consequence, the travel time along a perturbed ray path equals to first order the travel time along the true ray path. This allows us to linearly relate a small travel-time perturbation At to the slowness perturbation A-s0 along the reference raypath as To solve equation (3.7) numerically the model space has to be discretized. VanDecar (1991) has chosen to parameterize the model region over a regular grid using splines under tension to smooth variations between knots. In our study 36,540 knots spaced at approximately I by (3.6) L (3.7) 3 METHOD 34 Figure 3.5: Taut spline parameterization of model region. The model region is described by 36,540 knots (36 knots in latitude, 35 knots in longitude and 29 knots in depth). Note that the region outside the bold lines will not be shown later since anomalies within that region are damped to zero to smoothly match the radial earth model. Values in curly brackets denote knot spacing within the inner region. Small squares denote stations. 35 km intervals in the central region describe the model (figure 3.5). Discretization of equation (3.7) leads to the description of the total travel-time perturbation being the sum 3 METHOD of travel-time perturbations associated with each knot 35 Ati = ,£PijAsj (3.8) 3 with index i denoting the ray, index j the slowness knot and Pij — | jv being the matrix of partial derivatives (Frechet derivatives). Initial ray paths through the target region are computed using the iasp91 reference earth model. The linear equations (3.8) are weighted by the uncertainty estimates output by MCCC. Additional equations account for station time corrections c and event correc-tions e as mentioned in section 3.2.1. The whole system of linear equations becomes A x = b (3.9) T where A = ^ W P W C W E ) and x = ( A s c e ) and b = ( W A t ) and i ray number P,j matrix of partial derivatives A s slowness perturbations j model knot C i r c station corrections r station number E,-g e event corrections q event number Wjj weighting matrix 3 . 2 . 2 . 2 M o d e l c o n s t r u c t i o n b y l i n e a r i n v e r s i o n In inverting equation (3.9), the objective is to find a minimum structure solution for As. By minimum structure solution we mean the simplest model required to explain the data. In fact, VanDecar's code allows for a combination of smoothing and flattening to take place. The objective becomes minimize the model objective function <f>m — A s ' | | 2 + A|| A s " | | 2 (3.10) subject to data misfit <f>d = | |A* x — b | | 2 3 METHOD 36 where || || 2 denotes the L2 norm. Smoothing and flattening are described by centred and forward differences of the slowness values at each knot, respectively. Equation (3.10) is solved with an augmented matrix technique (Oldenburg, 1994), by which the smoothing and flattening equations form additional rows in equation (3.9). The augmented system is solved using a conjugate gradient method. Equations with large residuals are down-weighted by means of an L i norm (Bostock & VanDecar, 1995). At the five open sides of the model region (surface excluded) the slowness perturbations and their derivatives are gradually damped to zero to merge the model with the outside radial earth model. Model appraisal T Since we are not attempting to calculate the inverse of the matrix A A in equation (3.9) we are unable to calculate model covariance and resolution matrices, and can only indirectly assess how well a model may represent the Earth structure that gave rise to the travel-time perturbations. A technique commonly used to check the resolution that can be obtained by the present ray distribution is called a spike test. A synthetic model with isolated velocity anomalies is created and synthetic data calculated using the same ray geometry pertinent as for the experiment. These data are then inverted and the result compared to the original. One spike represents the smallest perturbation we can expect to resolve and it corresponds to one column of the resolution matrix. A way to explore the model space is to run a number of inversions using different values of regularizations A and \i (see equation (3.10) ). This results in a tradeoff curve between data misfit and model norm: a smooth (and/or flat) model gives rise to a large misfit whereas one fitting the data produces models with too much structure. The construction of the tradeoff curve allows (i) the selection of a preferred solution where tradeoff is not biased in favour of either misfit or structure, and (ii) the assessment of which structures 3 METHOD are stable and common to all models. 37 3 . 2 . 2 . 4 N o n - l i n e a r i n v e r s i o n At this point it should be stressed that the ray paths used to calculate the Frechet de-rivatives (equations (3.8) and (3.9) ) for the linear inversion are those for a ID earth as described by iasp91. This assumption provides us with a first estimate of earth structure. The non-linearity of the problem stems from the fact that the raypath along which a wave travels from source to receiver depends on the velocity structure. Because the reference model raypaths do not represent the true raypaths any linear inversion will be in error. In the non-linear approach we recalculate the ray paths and the partial derivatives of slowness with respect to time using the model derived in the linear inversion. On the basis that the inverse problem is only weakly non-linear, i . e. the difference in raypaths' position is small compared to the size of the heterogeneities, VanDecar's method approximates the solution by iteratively performing ray tracing and linear inversion as outlined in sections and We stop this iterative process when the model norm does no longer change significantly from the previous iteration. In the linear case, we calculated the effect slowness perturbations have on the travel-time along segments of the ray path of a spherical earth model. Following Snell's law, rays will preferably travel through faster regions. Travel-time perturbations therefore tend to enlarge regions with fast velocity during linear inversion. The non-linear inversion focuses these anomalies. 38 4 Results The last chapter introduced VanDecar's method and we now show the results obtained from applying it to our data set. A brief examination of travel-time residuals is followed by results from tomographic inversion. 4 . 1 C o m p a r i s o n o f r e s i d u a l t r a v e l t i m e s After obtaining relative delay times for several seismograms from one event these may be compared to the expected iasp91 travel-times. The difference between the observed time and the expected arrival time is termed a travel-time residual. After subtraction of the mean for each event, positive residuals denote slow arrivals and negative residuals denote fast arrivals. The residuals may be plotted and contoured on a map, thereby allowing a comparison of residual patterns. Although static station terms are ignored, diagrams of this kind are a helpful qualitative aid for identifying major features. Figure 4.1 shows four examples. The most obvious feature for all subplots in figure 4.1 is a trend to fast arrivals towards the north. Both a decrease in site elevation and a decrease in thickness of the sedimentary cover with slow velocities contribute to this effect (compare to figures 1.1 and 2.1). It should be noted that the trend is usually SE-NW for events arriving from southern azimuths. The trend is generally SW-NE for events from northern azimuths, which is the expected trend, however, for events arriving from southern directions this changes more to a SE-NW trend, implying large scale changes in the underlying mantle. In all subplots, station BIR is associated with a slower arrival than areas to the W, E and N . This may be due to peculiarities of a particular instrument or very local crustal structure. 4 RESULTS 39 95/11/13 08:43:14.7 5.9Mb East of Lake Baykal, Russia 96/02/16 09:44:58.1 6.0Mb North of Ascension Island Figure 4.1: Contour plots of four events arriving from different azimuths. Contour interval is 0.2 s; dashed lines are negative contours. Arrows point in azimuthal directions to source. Ellis & Hajnal (1992) chose to correct travel-time residuals for sedimentary thickness, elevation and ellipticity. Furthermore they averaged residuals for sets of earthquakes in the Kurile Islands and from South America. Their results show significant differences between these two sets, with the main differences apparent at the south-eastern stations of the 4 RESULTS 40 1991 array. It should be stressed here that their data-set will be incorporated into our tomography study. No attempt has been made here to correct for site elevation, crustal structure or different instrument type as these will be absorbed by the station correction terms during inversion. 4 . 2 T o m o g r a p h i c m o d e l s We now compare models derived by a linear inversion of the data set at different levels of regularization to identify a model which displays the desired balance between travel time misfit and model roughness. Following this, we will show results of a resolution test which indicates which anomalies in our model region can be reliably interpreted. Finally, we will show the result of non-linear inversion of the travel-time residuals. 4.2.1 Models from linear inversion Results of the linear inversion of teleseismic travel-time delays are shown in figure 4.2. Black shades in the plots denote areas which are sampled by fewer than four rays. The inversion was run at several levels of regularization in order to explore the model space for consistent anomalies, i.e. to distinguish between anomalies required by the data and spurious anomalies resulting from poor ray coverage. Anomalies are in general small, varying between ±2% for the least regularized to ±0.5% for the most regularized model. Stable features are a fast (negative slowness perturbation) kidney shaped anomaly centred at (51.5N, 104W, 200km) and slow (positive slowness perturbations) localized anomalies at (52.5N, 106W, 100km) and (50.5N, 106W, 100-200km). Figure 4.3 shows the tradeoff curve corresponding to the various inversions performed. It plots the root mean square of all remaining travel-time residuals as function of the model roughness (defined as second RESULTS Figure 4.2: Linear inversion with different levels of regularization. Regularization decreasing by factors of two from top to bottom. A is smoothing weight and u is flattening weight. Left column: depth sections at 100km, center column: depth section at 200 km, right column: NW-SE cross section. The preferred model corresponds to a smoothing weight of A=17,500. 4 RESULTS 42 «2 •a 0 0.005 0.01 0.015 0.02 0.025 0.03 -3i 0.035 0.04 model roughness [s km"J] Figure 4.3: Tradeoff curve for linear inversions. Open circles denote models shown in figure 4.2. Values denote smoothing weight, flattening weight was changed by factors of two also. The preferred model corresponds to A =17500 at the elbow of the tradeoff curve. order derivative of the slowness perturbation model) and, in a qualitative sense, allows us to determine the preferred level of regularization corresponding to a balance between data misfit and model roughness. The optimum level seems to be given at a smoothing weight of 15,000 to 40,000. Since the anomalies are small we tend to be in favour of a lower level of regularization and our preferred model corresponds to a smoothing weight of 17,500 shown in the center row in figure 4.2. 4 RESULTS 43 4 . 2 . 2 R e s o l u t i o n t e s t A spike test was performed to assess where within the model region structures are well re-solved with the given source-receiver configuration. To perform the spike test, a synthetic model consisting of alternating positive and negative spikes was created (see figure 4.4). Spikes were chosen to have either ±5% slowness anomaly and separated by 5 knots (170 km) in the horizontal directions and every 4 knots (130 km) in the vertical direction, respect-ively. Synthetic travel times were calculated by ray tracing along each of the ID earth ray paths defined by the source-receiver pairs used for the actual inversion. These times were contaminated with noise (0.1 s standard deviation) and then inverted for slowness structure using the same values of regularization parameters, A and ft, used for the preferred model in linear inversion (see figure 4.2). The result is shown in figure 4.5. A comparison of the recovered model (figure 4.5) with the synthetic model (figure 4.4) shows that the magnitudes of the anomalies were not well recovered. This was to be ex-pected as a direct consequence of regularizing the problem using first and second derivative damping. Within a volume defined by a pseudo-rectangular prism with opposing corners at (50°N, 108°W, 0 km) and (54°N, 102°W, 400 km) and sides running along gridlines, the locations and dimensions of the spike anomalies were recovered quite well. In particular, reasonable depth resolution to 400 km was achieved as can be seen along the column of three anomalies on top of one another at (52.33°N, 105.0°W) and the recovered altern-ating pattern at all depths shown in the figure. Outside this volume anomalies were less well resolved, note the smearing in the southern part of the cross section along 102.5°W. However, the fact that the stable features identified in section 4.2.1 he within the volume provides confidence in their veracity. RESULTS 44 cross section aloug 50.66N 110W 105W 100W cross section along 107.5W cross section along 52.33N cross section along 54.QN now losw IOOW now losw IOOW cross section along 105.0W cross section along 1G2.5W P-wave % slowness anomaly Figure 4.4: Synthetic spike model with alternating positive and negative anomalies. Top row: E-W sections; middle row: depth sections; bottom row: N-S sections. Little squares in middle row panels indicate station locations. RESULTS cross section along 50.66N cross section along 52.33N cross section along 54.0N now i05w IOOW now IOSW IOOW now iosw IOOW P-wave % slowness anomaly Figure 4.5: Recovered spike model. Top row: E-W sections; middle row: depth sections; bottom row: N-S sections. Note that the scaling differs from that for figure 4.4. 4 RESULTS 4.2.3 Models from non-linear inversion 46 Having determined a preferred preliminary model and regularization parameters after lin-ear inversion (see figure 4.2) ray paths were recalculated and the inversion rerun at the same level of regularization. The result is shown in figure 4.6. The data misfit did not change significantly (85% reduction in data misfit in both cases), therefore no further itera-tions were performed. The maximum slowness perturbation is +1.40% at (52.8N, 105.9W, 100km) and the minimum value is -1.33% at (51.8N, 103.9W, 170km). A well resolved, extensive region showing little variation exists in the NE of the study area. In comparing the results after non-linear inversion (fig 4.6) with those after linear inversion (figure 4.2, middle row) one can see that the slow anomaly increased in amplitude from +1.12% to +1.40% and that the fast anomaly became more focussed. RESULTS depth 100km Long i tude W e s t '" p . w a v e % s | o W n e s s anomaly Figure 4.6: Model after non-linear inversion. 4 RESULTS 48 4 . 2 . 4 S t a t i o n a n d E v e n t c o r r e c t i o n s Station correction terms for the stations used in inversion are shown with their equipment in figure 4.7. Not all stations (see figure 2.1) were used because VanDecar's code requires a minimum of ten records per station. Since we have used two totally independent data v short period PRS-4 (1991) • broad-band PRS-4 (1991) O broad-band NARS O broad-band WITS O broad-band IRIS station Figure 4.7: Plot of station corrections. Symbols denote equipment used and colours the station corrections. Contours were plotted using station corrections for the "short-period PRS-4" stations denoted by both types of triangles, contour intervals are 0.2 s. 4 RESULTS 49 sets (from 1991 and from 1994-96) station corrections between these two sets will have a shift. Station GRY was equipped with the same equipment during both deployments but was treated as two separate stations. The corrections output by the inversion differed by as much as 3.352 s for the two GRY stations; this value indicates a baseline shift in station correction between the two independent data sets. Corrections of similar magnitude were found for the station pairs SHE and HLB as well as SND and CHO. Corrections plotted for the 1991 stations in figure 4.7 were shifted by this amount and a new zero mean calculated. A strong correlation exists between station equipment and station correction as can be seen by comparing stations RAY and SEM or TOB and CHO, or even at the multi-equipped stations USK and LBL. This in retrospect demonstrates the necessity of treating different equipment at the same location as individual stations. These differences could be used in a similar way as described above for the two data sets to account for different delay times of different instrument times. This was not attempted here because the PRS-4 short-period stations are highest in number and contour lines of their correction terms show the underlying trend across the array quite clearly (figure 4.7). We note that the decrease in station correction reflects the general trend we observed in section 4.1 (figure 4.1), in particular, we note that the difference in times of ~0.6 s across the array is reproduced as well as the anomaly at station BIR noted previously. A comparison of residual times and station corrections using these stations is sufficient to conclude that most of the traveltime residuals can be attributed to crustal and local structure and do not reflect upper mantle structure. Examination of event corrections is useful in order to assess the effect of mantle hetero-geneity outside the region of study since systematic correlation in event corrections would reflect projection of such structure. It was mentioned in section 3.2.1 that lower mantle structure and near source structure find their way into the event corrections, the most 4 RESULTS 50 obvious example being rays that run faster along a descending slab. Figure 4.8 shows that most event corrections are close to zero after applying two baseline shifts to the data: the first arises from the fact that we have used two independent data sets, the correction applied was the negative value of the one applied to the station terms as described above; the second shift was done to compensate the corrections for a systematic error due to pick-Figure 4.8: Plot of event time corrections. A mean of 0.250 s was added to all values to compensate for pick bias. Note that large positive and negative correction were applied to core-grazing rays. Scalebar is in s, circle denotes 104° epicentral distance. 4 RESULTS 51 ing the first maximum of the phase rather than its onset (see figure 3.2), this correction was exactly 0.250 s as expected for the time between onset and maximum of a 1 Hz sine wave. (Note that a narrow bandpass filter around 1 Hz was applied to the data prior to cross-correlation, see chapter 3.) Large positive and negative values are found for earth-quakes with epicentral distances approaching 104°. These correspond to core grazing rays for which the actual raypath may differ from the calculated raypath. 4.3 Summary of inversion results The anomalies found by non-linear inversion deviate by less than ±2% from the iasp91 radial Earth model and this magnitude does not take into consideration a possible baseline shift due to the position of the model region at the edge of a large scale positive mantle anomaly. The maximum slowness perturbation found at (52.8°N, 105.9°W, 100km depth) is quasi-cylindrical with a radius of 50 km and extending from 60 to 220 km depth. The minimum slowness perturbation at (51.8°N, 103.9°W, 170km depth) is somewhat kidney shaped and its bottom at 250-280km is not as well defined as for the maximum due to smearing at the model edges. Several smaller features are seen: two other slow anomalies of smaller magnitude along 106°W are well resolved at the top but smear at the bottom, and one fast anomaly at (51.6°N, 106°W) connects to the larger anomaly between 100 and 200 km depth. The mantle below 270 km shows very low levels of heterogeneity (±0.5%). A spike test showed the two main features (maximum and minimum anomalies) to be well resolved and the mantle below those features (between 220 to 270 km) and 400 km also to be well resolved. These observations are here classified as the "robust observations" and the interpretation will largely focus on them. 5 Interpretation 52 In the previous chapter we have identified the location of anomalies that deviate from the average seismic velocities for the iasp91 radial Earth model. Our interpretation will focus mainly on the largest positive and negative anomalies and the homogeneous mantle below. We will examine possible causes of the anomalies, compare their locations to other geological and geophysical information we have for the region and then offer a hypothesis which explains the observations. 5.1 Possible causes for the seismic anomalies Velocity anomalies found in tomographic inversions may be caused by temperature, anelas-ticity, anharmonicity, partial melt, composition or anisotropy (e.g. Humphreys & Dueker, 1994; Sobolev et al, 1996), or a combination of these. This section will examine their effects on seismic velocities. The effect of temperature on seismic P-wave velocities is usually given as ^ « —0.5ms _ 1 K _ 1 (e.g. Anderson & Bass, 1984). This relationship was calculated at room temperature for a pyrolite model composition assuming a linear temperature dependence and is often extrapolated to high temperatures. Using this relationship a 250°K increase in temperature is necessary to explain a 1.5% increase in slowness such as that observed. A total temperature variation of 500°K between high and low velocity regions would bring temperatures close to or above the solidus (see figure 5.1). The temperature derivative of velocity mentioned above does not take into consideration non-linear effects like min-eral reactions, dissipation effects and partial melting that will occur at high temperature. Sobolev et al. (1996) interpreted 3-D seismic tomographic models in terms of temperat-ure, degree of partial melt and rock composition. They pointed out that contributions to 5 INTERPRETATION 53 depth (km) 0 156 304 443 573 I 1 1 1 1 pressure (GPa) Figure 5.1: Geotherm of Archean lithosphere compared to mantle solidus solid line: geotherm for Archean lithosphere after Morgan (1984), dashed line: dry mantle solidus after Herzberg et al. (1990), conversion from depth to pressure after Dziewonski & Anderson (1981). seismic velocities from mineral reactions, anelasticity and anharmonicity1 are significant and all act to lower seismic velocities with increased temperature. In one example cal-culation, they show that these effects will decrease the temperature difference necessary to explain a velocity perturbation of -3% from 432°K to 225°K. Following their argument we may assume that the temperature contrast between high and low velocity anomalies is significantly less than the beforementioned 500° K in our case and probably closer to 200° K. The most important effect of composition on seismic velocities is the iron content in 1 I n a harmonic solid the potential energy is proportional to the square of the displacement of atoms and all changes of the elastic moduli with pressure and temperature are controlled by changes in volume. Anharmonicity accounts for deviations from this model (S. Sobolev, 1997, personal communication). 5 INTERPRETATION 54 olivine and garnet, two major constituents of the upper mantle. An increase in iron at the expense of magnesium in these minerals by small amounts will decrease P-wave velocity and increase density, for example changing the iron content in olivine from 8% to 10% decreases the compressional velocity by -1% and increases the density by +1% (see data compilations by Gebrande, 1982, and Anderson et al, 1992). This relationship between seismic velocities, density and iron content is one of the main arguments in Jordan's (1978, 1979) tectosphere hypothesis. He argues that lithospheric roots evolved from a garnet lherzolite by extraction of basalt, i.e. depletion in iron with respect to the surrounding mantle, leaving harzburgite, a less dense and seismically faster residuum. Anisotropy, though certainly present, is believed to be a subordinate cause for P-wave slowness anomalies in our case as will be discussed in the next section. Physical parameters like preferred mineral orientation have to exhibit lateral variations to explain the observed seismic anomalies but this would not alter other observables like the gravity signature. Therefore a combination of temperature effects and, to a lesser extent, compositional effects is thought to be the most likely cause for the observed perturbations. 5.2 Comparison of tomographic results to other data sets The Trans-Hudson Orogen has been a focus of research over the past decade. A variety of studies have been outlined in chapter 1. We now compare our results to those from several of these studies before introducing an explanation for the tomographic image. 5 . 2 . 1 S u r f a c e g e o l o g y a n d i n f e r r e d s u b s e d i m e n t a r y g e o l o g y No obvious correlation between the locations of anomalies and the outlines of tectonic terrains inferred from aeromagnetics is possible (figure 5.2). However, it is now well es-5 INTERPRETATION 55 a) tomographic result (100 km depth) b) surface tectonics and aeromagnetics c) depth to Moho d) filtered Bouguer gravity Figure 5.2: Graphical comparison of tomographic results with other data sets. a) Result of non-linear tomographic inversion at 100 km depth. The +0.7% and -0.7% contours are highlighted in all plots. b) Tectonic units of the Trans-Hudson Orogen and inferred boundaries un-der the platform. Locations of diamondiferous kimberlites and high kim-berlitic mineral concentrations are indicated (for legend see figure 1.1). c) Depth to Moho from seismic reflection, points along seismic lines used for digitizing are shown as squares. d) Filtered Bouguer gravity data (passband centered on 200 km). tabhshed from reflection and refraction seismic surveys and drillcore sampling that the Glennie domain and the Hanson Lake block are underlain by the Saskatchewan craton 5 INTERPRETATION 56 (Ansdell et al, 1995; Baird et al, 1996; Chiarenzelli et al, 1996) which is not limited to the aeromagnetic outline of the Glennie domain but extends southward to Montana and South Dakota. This would place the high seismic velocities at depth between 107°W and 102°W, in particular the kidney shaped anomaly, beneath the surface boundaries of the Saskatchewan craton. We interpret the high velocity anomalies as being associated with a lithospheric keel of the Saskatchewan craton. This interpretation is consistent with the hypothesis that three Archean basement windows outcropping within the Glennie domain and Hanson Lake block are a tiny surface expression of a very large tectonic unit within the THO. The locations of diamond occurences and high concentrations of diamond-indicator minerals are shown on figure 5.2b. Diamondiferous kimberlites at Sturgeon Lake, Fort a la Come and Candle Lake and high concentrations of indicator minerals south-west of Assiniboia are all situated above or near the rims of low velocity anomalies. A complicating factor is that the two kimberlite bodies at Sturgeon Lake which are known to have been ice-rafted to their present position (Kjarsgaard, 1996) and it is not known how far they have journeyed. However, it is known that the Laurentide Ice Sheet moved from the NNE into Saskatchewan (Dyke & Prest, 1987) and therefore it is most likely that these kimberlites followed a path above low mantle velocities (see figure 5.2a and b). Thus the correlation between locations of diamond and diamond indicator minerals and low velocities in the upper mantle is a strong indication that the genesis of Saskatchewan kimberlites is correlated to these anomalies. In addition, the fact that these kimberlites are diamondiferous allows us to infer that thick lithosphere was present at the time of their eruption. Either the diamonds had been preserved in a lithospheric keel and were picked up by ascending kimberlites or they were brought up from below the asthenosphere and then erupted from the base of the lithosphere. In both cases a thick lithosphere that reaches 5 INTERPRETATION 57 into the diamond stability field is necessary to prohibit the retrogression of diamond to graphite (Morgan, 1995). 5 . 2 . 2 S e i s m i c r e f l e c t i o n a n d r e f r a c t i o n a n d s h e a r - w a v e s p l i t t i n g Values for the depth to the Mohorovicic discontinuity (Moho) shown in figure 5.2c were di-gitized along seismic refraction lines from Nemeth et al. (1996) and then contoured. Data points in the south mark the thickening of the crust towards the centre of the Williston basin. Seismic refraction data points are sparse in the western central part of the area, in particular above our low velocity anomaly. The seismic line along 104°W ends above the centre of our high velocity anomaly and shows extremely variable Moho topography with depth values ranging between 37 and 51 km. Hajnal et al. (1997) speculate that the deep Moho together with underlying high velocity and a layered uppermost mantle are a consequence of the tectonic collision between Hearne-Rae, Superior and Saskatchewan cratons in the Proterozoic. (Note that the high velocity may be a manifestation of aniso-tropy, see below.) With respect to the thickness of the lithosphere below orogens, Kincaid & Silver (1996) argue that strain energy generated by orogeny destroys the subcontinental lithosphere. Dissipation of energy heats up the downward advecting material which will therefore not show up as a high velocity zone in tomographic studies. The deformation of the lithosphere during orogeny does however give rise to anisotropy as evidenced by splitting of shear waves (Silver & Chan, 1991). Following this idea we expect thickened Moho under orogenic belts to be underlain by thinner lithosphere than in the adjacent Archean regions. The shallowest Moho lies at 104°W, 54.2°N within the inferred extension of the Glennie domain and is underlain by low velocities in the upper mantle. On the other hand the deepest Moho at 104°W, 52.6°N lies above a transition from low velocity to high velocity. It is clear that depth to Moho does not correlate to our tomographic results. 5 INTERPRETATION 58 Since we expected such a correlation our interpretation must explain why our results do not conform to this. Seismic anisotropy for the upper mantle in the study area has been demonstrated by seismic refraction studies and teleseismic studies (e.g. Hajnal et al, 1997; Ellis et al., 1996). Seismic refraction information on the upper mantle from the velocity of the Pn phase shows high velocity blocks (8.4 km/s) in the 50 to 90 km depth range (Hajnal et al, 1997). An east-west crossing line failed to see high velocity block but instead imaged a velocity of 8.0 km/s, for which anisotropy is a possible explanation (Nemeth & Hajnal, 1996). This interpretation would be consistent with results of Morel-a-l'Hussier et al. (1987) further south where faster upper mantle velocities were also determined along north-south trending lines, i.e. along strike of the orogen. The segments showing high seismic Pn velocities are represented in figure 5.3. The existing teleseismic studies in this region investigating the splitting of SKS waves indicate that anisotropy characterizes the whole lithospheric column over the whole region. Silver &; Kaneshima (1993) report on an experiment that crossed the THO south of the study area. Their results show the polarization of the fast direction to change from an E-W direction within the THO to SW-NE directions in adjacent regions. This observation was not supported by Ellis et al. (1996) for stations further north where this direction is still SW-NE. Both studies, though disagreeing in splitting direction, render travel time differences between fast and slow split wave of f«l sec for stations within the inferred boundaries of the THO (fig. 5.3). According to Silver & Chan (1991) this delay time corresponds to an anisotropic layer having for example an inherent anisotropy of 4% and a vertical extent of 115 km. Thus the absence of high velocities in the tomographic image corresponding to the up-permost mantle beneath 54°N, 104°W may be simply the result of anisotropy. In controlled source studies rays travel subhorizontal below the Moho but when earthquakes between 5 INTERPRETATION 59 Figure 5.3: Results of studies relating to upper mantle anisotropy: Splitting results taken from Ellis et at. (1996) and Silver & Kaneshima (1993). Circles denote time difference between the fast and the slow split waves and lines show azimuth of fast polarization direction. Refraction data from Morel-a-l'Hussier et al. (1987), Nemeth et al. (1996) and B. Nemeth, 1996, personal communication. Small dots indicate location of refraction lines and bold lines indicate high Pn velocities (8.4 km/s). Bold area and contours are our slowness model at 100 km depth. 30° and 100° epicentral distance are used as sources they pass at angles of 20° to 40° from the vertical. If we assume a horizontal axis of hexagonal symmetry for the upper mantle then Pn rays from controlled sources run within the plane containing axes of maximum 5 INTERPRETATION 60 and minimum velocity whereas steep rays from earthquakes do not sense anisotropy in P velocities. A further point of note is that the high velocities are shallower than 90 km depth where tomographic resolution becomes poorer since rays become steeper and inversion may attribute travel time variations to station corrections rather than slowness perturbations within the model. Refraction seismology and SKS studies indicate the presence of anisotropy, most likely 3D variations in anisotropy are present. The inversion code used in our study assumes an elastically isotropic earth. A recent study by Gresillaud & Cara (1996) attempted to invert for an anisotropic structure using P-wave travel-time delays. They used initial models comprising a combination of layered hexagonal anisotropy and P-wave velocity perturbations. This model was recovered reasonably well for an ideal ray distribution. However, when applying the inversion technique to data from experiments in the Rhine Graben and the Pyrenees, no reliable inferences on anisotropic structure could be made since the results were similar to those obtained using an isotropic inversion technique. The main complicating factor is inadequate ray sampling, and a priori assumptions would be necessary to constrain the direction of the fast axis. They concluded that 3D velocity dis-tribution recovered by inversion codes assuming an isotropic earth give a reasonable picture even if anisotropy is present in the lithosphere. We thus expect that we are extracting a representative isotropic component. 5 . 2 . 3 G r a v i t y Gravity values shown in figure 5.2d come from the Western Canadian Whole File provided by the Geological Survey of Canada (Natural Resources Canada, Ottawa) and were filtered using a bandpass filter with a peak wavelength at 200 km. This choice of passband is justified since a spherical anomaly situated at 100 to 200 km depth (i.e. the depth region 5 INTERPRETATION 61 we are most interested in) gives rise to a vertical gravity signatures with a half width of 150 to 300 km, respectively. The assumption of a spherical anomaly is a reasonable first order estimate for upper mantle density anomalies although near surface structures like sedimentary basins can also cause anomalies of similar wavelength provided they have sufficient lateral dimension. Morgan (1995) suggested that attention be focussed on features with wavelengths of ~200 km when looking for lithospheric roots using potential surveying. The gravity map shows a general trend from low values in the west to high values in the east which reflects the transition from the Western Canada Sedimentary Basin to the Canadian Shield. A comparison of filtered Bouguer gravity with locations of slowness anomalies suggests the correlation of the high velocity region with a gravity low at about 51.5°N, 103.5°W which is consistent with the notion of a keel that has a high velocity and low density (Jordan 1988). We also note that the low velocity region coincides with a relative gravity high. This is consistent with an increased iron content, but a higher temperature would lower the density and thus compensate for compositionally increased density. If the correlation between seismic structure and gravity is meaningful this would imply compositional contrasts in the lithospheric mantle. 5.2.4 Heat flow Correlating surface heat flow patterns to upper mantle structure is difficult since mantle contributions to heat flow may be masked for a variety of reasons (Nyblade & Pollack, 1993, and references therein), such as variations in crustal heat production, groundwater circulation, heat refraction in shallow subsurface structures, blanketing effects of sediment formations, surface erosion and local climatic variations. Nevertheless variations in heat 5 INTERPRETATION 62 flow may provide constraints on the origin of our anomalies. Majorowicz et al. (1986) obtained thermal data from boreholes drilled at 30 km spac-ings across Alberta, southern Saskatchewan (south of 52°N) and southwest Manitoba. They distinguished between heat flow in the Mesozoic strata and below the Paleozoic surface and found heat in the Mesozoic strata being redistributed by water movement. This effect is not expected for the Paleozoic sediments which therefore would give a better picture of the heat flow input from below. Data below the Paleozoic surface for Saskatchewan (figure 5.4) document large variations in heat flow. A prominent heat flow high (>100 mW/m 2) in the SE (Weyburn area) is underlain by low mantle velocities (although resolution of our model is poor in this area; figure 4.6) and low heat flow values (40 mW/m 2) at 51°N 103°W are underlain by high velocities. Measurements taken at a deep well (2,200 m) into the Precambrian basement at Regina yielded a heat flow of 51 mW/m 2 for this basement (Jessop & Vigrass, 1989). This value is lower than data by Majorowicz et al. (1986) for that area (see figure 5.4) because thermal conductivities were measured directly on the cutting, and not estimated on the basis of geological logs (T. Lewis, 1996, personal communication). Studies by Drury (1985) and Guillou-Frottier et al. (1996) in the exposed part of the THO show average low heat flow values with the variability attributed mainly to surface effects. A distinction between the Flin Flon Belt (42±3 mW/m 2) and the Thompson Belt (54±8 mW/m 2) was attributed to differing radiogenic heat production within the upper 13 km of the crust which adds upon a uniform heat flow of 38 mW/m 2 from below. Very local high temperatures in the Thompson Belt (81 mW/m 2) were attributed to refraction of heat in highly-conductive quartz-rich metasedimentary formations. Since no heat flow measurements have been made in the central part of our study area, 5 INTERPRETATION 63 no comparison of the largest low velocity anomaly with heat flow data can be attempted. The largest high velocity anomaly shows correlation mainly with low heat flow values. In the south where dense heat flow measurements were made we can make some correlation between high heat flow and low-velocity mantle structure. Low heat flow average in the north relates to mantle velocities that do not show large heterogeneities. Figure 5.4: Results of heat flow studies: Results in the south show heat flow below the Paleozoic surface (Majorowicz et al, 1986) and results in the north are from Drury (1985) and Guillou-Frottier et al. (1996). Also shown is a value for heat flow from the Precambrian basement from a deep well at (50.4°N, 104.6°W) in Regina (Jessop & Vigrass, 1989). Our study area and contours for slowness model at 100 km depth are indicated. 5 INTERPRETATION 64 5 .2 .5 S u m m a r y o f c o m p a r i s o n The comparisons done in this section can be summarized as follows: - high and low upper mantle velocities fall within the outline of the proposed Saskatchewan craton - diamondiferous kimberlites are located directly above or close to low velocity anom-alies - there exists no obvious correlation between Moho depth and upper mantle structures - seismic anisotropy may have some effect but cannot alone be responsible for the observed anomalies - high seismic velocities correlate to a gravity lows and lower seismic velocities to gravity highs - there is a suggestion that high heat flow in the south relates to low velocity anomalies The following hypothesis is based on these observations. 5.3 Hypothesis: Thermomechanical erosion of a lithospheric keel by a mantle plume conduit In this section we want to describe a hypothesis that is consistent with our results and the other data sets presented in the previous section and the basis for this proposition. Our hypothesis evokes a hotspot in the form of either a mantle plume conduit or small scale convective instability that eroded into the lithosphere of the Saskatchewan craton 100 Ma ago leading to the eruption of kimberlitic magmas. The zone of erosion is now trapped 5 INTERPRETATION 65 Figure 5.5: Schematic sequence illustrating proposed hypothesis. a) Mantle plume or asthenospheric instability impinging on lithosphere of Saskatchewan craton in the Cretaceous. b) Thermomechanical erosion of lithosphere triggers diamondiferous kim-berlites and erodes to 100 km depth within 20 Million years. cj Hot material from plume or instability trapped in lithosphere has heated the lithosphere until the present. in the lithosphere and is imaged as a low velocity anomaly, whereas the remnants of a lithospheric keel of the Saskatchewan craton show up as a high velocity anomaly. Figure 5.5 shows our hypothesis in a schematic sequence. A relationship between kimberlites and mantle plumes was first proposed by Green & Gueguen (1974). Their geochemical analysis of xenoliths from South African kimberlites led to the conclusion that kimberlitic magmas were a consequence of mantle diapiric up-welling. Crough et al. (1980) relocated kimberlites that were emplaced in North and South America and in Africa during the last 150 Ma to their place of origin in the hotspot reference frame. Relocation was done using the paleopoles of the continents derived from paleomagnetic studies. Kimberlites cluster at present locations of hotspots when rotated to their position at the time of eruption. The conclusion was that the hotspot concept predicts ages of kimberlites quite well. Haggerty (1994) reviewed the relationship between kimberlites, diamonds and mantle plumes. Inclusions of majorite in diamond point to a 5 INTERPRETATION 66 depth of origin of kimberlites at or below that of the transition zone. Since no heat is being generated in the transition zone a deep heat source is necessary. The D" layer could serve as source for sulphur, potassium and other trace elements as well as carbon, which in dia-monds show primordial composition (i.e. comparable to meteorites). Haggerty suggested that protokimberlitic magma rises from D" as highly reduced plumes with the reducing environment being necessary to preserve diamonds. The shape of inferred mantle plumes depends to a large degree on the viscosity contrast between plume material and mantle. Lab experiments (e.g. Olson &; Singer, 1985) have shown that a plume rises as a cylindrical diapir with a mushroom shaped head or a spherical chamber with a thin conduit. For the Earth's mantle they favour the latter shape because the higher temperature within the plume (sd300°K temperature difference) causes a sharp viscosity contrast. Their experiments also showed that a plume will break into individual plumes as the source moves steadily through a highly viscous fluid and that convection of the background fluid does alter the rise path but the position where the plume will impinge on the base of the lithosphere remains fixed. The interaction of a plume with the lithosphere was modelled by Davies (1994). The process involves softening of the lithosphere by heating and its removal by convection and is called thermomechanical erosion. Above the plume conduit, thermomechanical erosion causes rapid local thinning but the plume head itself is unlikely to erode into the lithosphere due to insufficient thermal inertia (see figure 5.5a). Davies calculates the effect a plume, 250°K hotter than the surrounding mantle, has on the lithosphere in several cases: (i) Oceanic lithosphere (initial age 90 Ma, thickness 120 km) has within 5 Ma penetrated to a depth of 100 km from the surface by a hot plume conduit. After 2 Ma the geotherm crosses the dry solidus causing basaltic volcanism. An example is the Hawaii islands for which the thickness of the plume conduit is estimated to be 70 km. (ii) For erosion into a thin 5 INTERPRETATION 67 Figure 5.6: Thermomechanical erosion of a thick continental lithosphere. Shown are the geotherm after 2, 5, 10 and 25 Ma, temperature is plotted against depth. Note that the plume erodes from 200 km to 100 km depth and that the solidus is just touched (figure taken from Davies, 1994). (120 km) continental lithosphere melting within the plume starts after 5 Ma and melting at the base of the crust is approached within 25 Ma. Erosion progresses to a depth of 50 km from the surface. An example is the Snake River Plain as the track of the Yellowstone plume, (iii) The effect of a plume conduit on thick (200 km) continental lithosphere is most applicable to the present case and is shown in figure 5.6. Only minor melting occurs in the plume after 20 Ma, no crustal melting is initiated and the erosion reaches a depth of 100 km after 25 Ma. This context allows the following interpretation of our anomalies. A plume conduit eroded into the thick lithosphere and caused minor melting of kimberlitic magma. Volatiles necessary to create this magma were brought up from deeper within the mantle and the diamonds were brought up from deep and/or picked up in the lithospheric 5 INTERPRETATION 68 keel (figure 5.5b). Having established the possibility that a plume conduit eroded into the lithosphere creating a cylindrical anomaly in the lithosphere we need to assess the thermal diffusion of this anomaly. Let us assume the thermal conductivity of the lithosphere to be K = 1.1 • 1 0 - 6 m 2 s - 1 (Stacey, 1992) and .4795AT to be a reasonable length scale (this value corresponds to in the case of a heating half space). Using the results outlined in appendix E, heat will have diffused 40 km away from a 50 km thick cylinder after 100 Ma i.e. the total diameter of anomaly has become 130 km, (or 160 km for a 70 km thick conduit as in the case of Hawaii) which corresponds to the scale of our low velocity anomaly (figure 5.5c). The sharp lower boundary of the slow anomaly and the relative homogeneity of the mantle below leads to the question of what has happened to the plume head and bottom part of the plume conduit. The answer is that the top part of the conduit that has eroded into the lithosphere is trapped there now and that the plate has moved away from atop the hotspot source. We interpret the lack of heterogeneity below 220 to 280 km as evidence for the convecting asthenosphere. This interpretation is in agreement with the interpretation of electromagnetic results to the west below central Alberta that become more laterally homogeneous around 250±50 km depth (D. Boerner, 1997, personal communication). We will test whether or not a plume gave rise to the anomaly by trying to trace this plume. Hotspot traces are usually associated with volcanic activity. It is possible that a plume breaks up into separate blobs during its ascent through the highly viscous mantle (Olson & Singer, 1985) in which case we may find an interrupted trace. Let us assume that the Fort a la Corne kimberlite field was caused by the passage of a hotspot 100 Ma ago. The trace of this hypothetical hotspot across the North American continent 5 INTERPRETATION 69 Figure 5.7: Trace of a hypothetical hotspot across North American Plate. Numbers denote time in Million years. 100 Ma ago this hotspot would have initiated the Fort a. la Come kimberlite field. An error circle for the uncertainty in location for a time span 100 Ma and the location of a present day low velocity region (van der Lee & Nolet, 1996) is added. and its estimated present-day location today are shown in figure 5.7 following the plate reconstructions done by Engebretson et al (1985). The error circle of 900 km diameter is that for a 100 Ma reconstruction and it encloses a V-shaped low anomaly imaged by van der Lee &; Nolet (1996) from Rayleigh waveform inversion where S-wave velocities are lowered by 2% with respect to the mean value. This anomaly coincides with changes in S-wave splitting (Wysession et al, 1996) and low P-velocities extending down to the transition zone inverted using P-wave travel-times (Taylor &; Toksoz, 1979). To explain these low velocities, van der Lee & Nolet (1996) follow an argument used by Nolet & Zielhuis (1994) who regard low seismic velocities in this area as evidence for the presence of water brought down by subduction of the Iapetus plate 400 Ma ago. The area was passed by the Great Meteor hotspot 180 Ma ago and by the Monteregian hotspot 100 Ma 5 INTERPRETATION 70 ago; both hotspots are apparently unrelated to the Saskatchewan lithosphere on the basis of timing considerations. Heating of hydrated upper mantle is supposed to have triggered crustal magmatism and the low seismic velocities are an indirect effect of water content (water content lowers the solidus, when heated the material will create partial melt at lower temperatures which in turn lowers seismic velocities). In the argument used by van der Lee & Nolet (1996) no convective mixing in the upper mantle but a coherent translation of the upper mantle with the North American plate is necessary to preserve the low P-wave anomaly over 100 Ma to a depth of 600 km. VanDecar et al. (1995) interpret a low-velocity anomaly extending to depths of 600 km below the Parana basin as remnant of the Tristan de Cunha plume that has translated coherently with South America for the past 130 Ma; if their interpretation is right the same argument could hold for the deep low velocities below the New England states. However, a present day mantle plume would enhance the seismic signature created by a hydrated upper mantle. If we interpret the anomalies seen at the East Coast as anomalies partly thermal in origin we would add evidence to our hypothesis that the Saskatchewan kimberlites were created by a mantle plume. Asthenospheric instabilities have been proposed as an alternative to mantle plumes by Tackley & Stevenson (1993). They modeled mantle rock being close to its solidus in the asthenosphere and found that an infinitesimal upward velocity could lead to a positive feedback situation between melt creation and buoyancy. The instability is similar to a Rayleigh-Taylor instability and capable of producing large volumes of magma within a short time. Implications of this mechanism, including lithospheric erosion, are currently being tested (Tackley, 1996). The conclusions drawn here for a plume would equally hold for a such short lived instability with the exception that no track could be observed. Because it is less restrictive than the plume hypothesis the mechanism of Rayleigh-Taylor-like instabilities is appealing. 5 INTERPRETATION 71 Our preferred explanation of the anomalies is the thermomechanical erosion of mantle material into the lithosphere. The seismic anomalies imaged by non-linear inversion of travel-time residuals are caused by a combination of thermal and compositional hetero-geneities in the upper mantle. The fast and slow anomalies originated separately: the fast anomaly is the remnant of an Archean keel which was eroded by a plume conduit or an asthenospheric instability, the slow anomaly, in the Cretaceous. The lithosphere is 220 to 280 km thick as manifested by low levels of heterogeneity below. 6 CONCLUSIONS 6 Conclusions 72 We now come back to evaluate the objective outlined in chapter 1. We identified four areas of research to which our study could contribute: the tectonic history of the Trans-Hudson Orogen, the structure of the upper mantle in the area, velocity anomalies beneath Archean cratons, and diamond finds in the study area. We have found well resolved seismic anomalies that may be explained as remnants of a lithospheric keel and of a mantle plume or instability that has eroded parts of this keel. This hypothesis is consistent with the following: - the Saskatchewan craton is in part underlain by a lithospheric keel, - the upper mantle below southern Saskatchewan is heterogenous at our level of sampling, - the lithosphere in this area is 220 to about 270 km thick, - the eruption of diamondiferous kimberlites may have been triggered by the ther-momechanical erosion of the lithosphere underlying the Saskatchewan craton via a mantle plume or a Rayleigh-Taylor like instability. A valuable test for our interpretation of the seismic anomalies was proposed by O'Reilly & Griffin (1996). Suites of xenoliths brought up from different depths define a paleogeo-therm at the time of eruption, allow one to infer the depth to the crust-mantle boundary and show the depth distribution of rock types and geochemical processes occuring at depth. Petrophysical parameters determined on xenoliths can constrain geophysical observations and linking geophysical studies to xenolith studies allows improvement of a one-dimensional picture of the upper mantle and crust to its three-dimensional understanding. If kimberlites 6 CONCLUSIONS 73 of different ages are present in the area the evolution of the thermal state and composi-tion of the lithosphere with time can be mapped and thus a four-dimensional image of the lithosphere achieved. In our case, a combination of our results with kimberlite studies could prove or disprove that the lithospheric keel remained unaltered until the emplace-ment of the kimberlites. Unfortunately, no xenolith data for the Saskatchewan kimberlites have been published and this method cannot be applied. Three other teleseismic experiments similar to ours are currently being conducted on Archean cratons actively mined for diamonds. In Canada, a two-dimensional array of a dozen broadband stations was installed in the fall of 1996 on the Slave province as part of the LlTHOPROBE SNORCLE transect with data acquisition expected to last about 9 months. In Australia, the SKIPPY project makes use of temporary broadband arrays consisting of 12 portable stations that move across the continent. In its final stage this array will be deployed on the western and northwestern part of the continent (Zielhuis & van der Hilst, 1996) where some of the world's economically most important diamond occurences are located. In South Africa, historically the leading supplier of gem-quality diamonds, an array of 50 broadband stations has been proposed to study the crust and upper mantle beneath the Kaapvaal craton (Carlson et al., 1996). All three projects focus on areas that are actively being mined for diamonds and for which a variety of geophysical, geological and penological data, including xenolith data, is available from other surveys or supporting experiments. The implications drawn from the teleseismic results of these projects will be a test for the validity of the conclusions we have drawn for the Trans-Hudson Orogen. 74 References Abbott, D., 1991. The case for accretion of the tectosphere by buoyant subduction, Geophys. Res. Lett, 1 8 , 585-588. Anderson, D. L. and J. D. Bass, 1984. Mineralogy and composition of the upper mantle, Geophys. Res. Lett., 1 1 , 637-640. Anderson, 0. L., D. Isaak and H. Oda, 1992. High-temperature elastic constant data on minerals relevant to geophysics, Rev. Geophys., 3 0 , 57-90. Ansdell, K. M . , S. B. Lucas, K. Connors and R. A. Stern, 1995. Kisseynew metasedi-mentary gneiss belt, Trans-Hudson Orogen (Canada): back-arc origin and collisional inversion, Geology, 2 3 , 1039-1043. Asudeh, I., R. Wetmiller and C. Spencer, 1993. LithoSEIS version 5.00 User Manual, Geological Survey of Canada Open File, N o . 2 6 7 9 , 105 p. Baird, D. J., K. D. Nelson, J. H. Knapp, J. J. Walters and L. D. Brown, 1996. Crustal structure and evolution of the Trans-Hudson orogen: Results from seismic reflection profiling, Tectonics, 1 5 , 416-426. Behler, R. E. and M . A. Lombardi, 1991. NIST - Time and Frequency Services, National Institute of Standards and Technology Special Publication, 4 3 2 , 23 p. Bostock, M . G. and J. C. VanDecar, 1995. Upper mantle structure of the northern Cascadia subduction zone, Can. J. Earth. Sci., 3 2 , 1-12. Brune, J. and J. Dorman, 1963. Seismic waves and Earth structure in the Canadian Shield, Bull. Seis. Soc. Am., 5 3 , 167-210. Buchbinder, G. and G. Poupinet, 1977. P-wave residuals in Canada, Can. J. Earth Sci., 14, 1292-1304. Carlson, R. W., T. L. Grove, M . J. de Wit and J. J. Gurney, 1996. Program to study crust and mantle of the Archean craton in Southern Africa, EOS, 7 7 , 273+277. Chiarenzelli, J. R., L. B. Aspler and M . Villeneuve, 1996. Characterization, origin, and Paleoproterozoic history of the Saskatchewan Craton and possible implications for Trans-Hudson Orogen, Report of Sixth Trans-Hudson Orogen Transect Meeting, Lithoprobe Report, 5 5 , 26-38. Clifford, T. N. , 1966. Tectono-metallogenic units and metallogenic provinces of Africa, Earth and Planet. Sci. Lett., 1 , 421-434. Crough, S. T., W. J. Morgan and R. B. Hargraves, 1980. Kimberlites: their relation to mantle hotspots, Earth and Planet. Sci. Lett., 5 0 , 260-274. Davies, G. F., 1994. Thermomechanical erosion of the lithosphere by mantle plumes, J. Geophys. Res., 9 9 , 15709-15722. REFERENCES 75 Dobrin, M . B., 1976. Introduction to geophysical prospecting, McGraw-Hill, New York, 3rd ed. Drury, M . J., 1985. Heat flow and heat generation in the Churchill province of the Canadian Shield, and their palaeotectonic significance, Tectonophysics, 115, 25-44. Dyke, A. S. and V. K. Prest, 1987. Late Wisconsinan and Holocene Retreat of the Lauren-tide Ice Sheet, Geological Survey of Canada, Map 1702A, scale 1:5,000,000. Dziewonski, A. M . and D. L. Anderson, 1981. PreUminary reference Earth model, Phys. Earth Plan. Int., 25, 297-356. Ellis, R. M . and Z. Hajnal, 1992. Investigation of the Properties of the Lithosphere using Teleseimic Waves, Report for Energy, Mines and Resources Canada and Cameco Corporation, 66 p. Ellis, R. M . , Z. Hajnal and M . G. Bostock, 1996. Seismic studies on the Trans-Hudson Orogen of western Canada, Tectonophysics, 262, 35-50. Engebretson, D. C , A. Cox and R. G. Gordon, 1985. Relative motions between oceanic and continental plates in the Pacific Basin, Geol. Soc. Am. Special Paper, 206, 59p. Fowler, C. M . R. and E. G. Nisbet, 1985. The subsidence of the Williston Basin, Can. J. Earth Sci., 22, 408-415. Gebrande, H., 1982. Elastic wave velocities and constants of elasticity of rocks and rock forming minerals, in: K.-H. Hellwege, Landolt-Bornstein - Zahlenwerte und Funk-tionen aus Naturwissenschaft und Technik, Gruppe V: Geophysik und Weltraum-forschung, Band lb, Springer-Verlag Berlin, 1-99. Gent, M . R., 1992. Diamond exploration in Saskatchewan, CIM Bulletin, 84, 64-71. Gresillaud, A. and M . Cara, 1996. Anisotropy and P-wave tomography: a new approach for inverting teleseismic data from a dense array of stations, Geophys. J. Int., 126, 77-91. Grand, S. P., 1994. Mantle shear structure beneath the Americas and surrounding oceans, J. Geophys. Res., 99, 11591-11621. Green, H. W. and Y. Gueguen, 1974. Origin of kimberlite pipes by diapiric upwelling in the upper mantle, Nature, 249, 617-620. Green, A.G. , Z. Hajnal and W. Weber, 1985. An evolutionary model of the western Churchill province and western margin of the Superior province in Canada and the north-central United States, Tectonophysics, 116, 281-322. Guillou-Frottier, L., C. Jaupart, J. C. Mareschal, C. Gariepy, G. Bienfait, L. Z. Cheng and R. Laponte, 1996. High heat flow in the Trans-Hudson Orogen, central Canadian Shield, Geophys. Res. Lett., 23, 3027-3030. REFERENCES 76 Haggerty, S. E., 1994. Superkimberlites: A geodynamic diamond window to the Earth's core, Earth and Planet. Sci. Lett, 122, 57-69. Hajnal, Z., B. Nemeth, R. M . Clowes, R. M . Ellis, G. D. Spence, M . J. A. Burianyk, I. Asudeh, D. J. White, D. A. Forsyth, 1997. Mantle involvement in lithospheric collision: seismic evidence from the Trans-Hudson Orogen, western Canada, Geophys. Res. Lett, in press. Helmstaedt, H. and D. J. Schulze, 1989. Southern African kimberlites and their mantle sample: implications for Archaean tectonics and lithosphere evolution, in: J. Ross, A. L. Jaques, J. Ferguson, D. H. Green, S. Y. O'Reilly, R. V. Danchin and A. J. A. Janse (eds.) Kimberlites and Related Rocks, Volume 1, GSA Special Publication, 14, 358-368. Helmstaedt, H. H. and J. J. Gurney, 1995. Geotectonic controls of primary diamond de-posits: implications for area selection, in: Griffin, W. L. (ed.), Diamond Exploration into the 21st Century, J. Geochem. Explor., 53, 125-144. Herzberg, C , T. Gasparik and H. Sawamoto, 1990. Origin of mantle peridotite: con-straints from melting experiments to 16.5 GPa, J. Geophys. Res., 95, 15779-15803. Hoffman, P. F., 1988. United plates of America, the birth of a craton: early Proterozoic assembly and growth of Laurentia, Ann. Rev. Earth Planet.Sci, 16, 543-603. Hoffman, P. F., 1990a. Subdivision of the Churchill province and extend of the Trans-Hudson Orogen, in: J. F. Lewry and M . R. Stauffer (eds.), The early Proterozoic Trans-Hudson Orogen of North America, Geol. Assoc. Can. Spec. Pap., 37, 15-39. Hoffman, P. F., 1990b. Geological constraints on the origin of the mantle root beneath the Canadian shield, Phil. Trans. R. Soc. Lond. A, 331, 523-532. Humphreys, E. D. and K. G. Dueker, 1994. Physical state of the western U.S. upper mantle, J. Geophys. Res., 99, 9635-9650. Jessop, A. M . and L. W. Vigrass, 1989. Geothermal measurements in a deep well at Regina, Sakatchewan, J. Volcanol. Geotherm. Res., 37, 151-166. Jones, A. G. and J. A Craven, 1990. The North American Central Plains conductivity anomaly and its correlation with gravity, magnetic, seismic, and heat flow data in Saskatchewan, Canada, Phys. Earth Planet. Int., 60, 169-194. Jones, A. G., J. A. Craven, G. W. McNeice, I. J . Ferguson, T. Boyce, C. Farquharson, R. G. Ellis, 1993. North American Central Plains conductivity anomaly within the Trans-Hudson orogen in northern Saskatchewan, Canada, Geology, 21, 1027-1030. Jordan, T. H., 1978. Composition and development of the continental tectosphere, Nature, 274, 544-548. REFERENCES 77 Jordan, T. H., 1979. Mineralogies, Densities and Seismic Velocities of Garnet Lherzolites and Their Geophysical Implications, in: F. R. Boyd and H. 0. A. Meyer (eds.), The Mantle Sample: Inclusions in Kimberlites and Other Volcanics - Proceedings of the Second International Kimberlite Conference, Volume 2, Am. Geophys. Union, 1-14. Jordan, T. H., 1988. Structure and Formation of the Continental Tectosphere, in: M . A. Menzies and K. G. Cox (eds.), Oceanic and continental lithosphere: similarities and differences, Journal of Petrology, Special Lithosphere Issue, 11-37. Kennett, B. L. N. and E. R. Engdahl, 1991. Traveltimes for global earthquake location and phase identification, Geophys. J. Int., 105, 429-465. Kincaid, C. and P. Silver, 1996. The role of viscous dissipation in the orogenic process, Earth and Planet. Sci. Lett, 142, 271-288. Kirkley, M . B., J. J. Gurney and A. A. Levinson, 1991. Age, origin and emplacement of diamonds: scientific advances in the last decade, Gems & Gemology, 27, 2-25. Kjarsgaard, B. A., 1996. Prairie kimberlites, in: LeCheminant, A. N. , D. G. Richardson, R. N . W. DiLabio and K. A. Richardson (eds.), Searching for diamonds in Canada, Geological Survey of Canada, open file 3228, 67-72. Lees, J. M . , 1995. Xmap8: a free program for three-dimensional GIS, Seis. Res. Lett., 66, 33-37. LeFevre, L. V. and D. V. Helmberger, 1989. Upper mantle P velocity structure of the Canadian Shield, J. Geophys. Res., 94, 17749-17765. Lehnert-Thiel, K., R. Loewer, R. G. Orr and P. Robertshaw, 1992. Diamond-bearing kimberlites in Saskatchewan, Canada: The Fort a la Come case history, Explor. Mining Geol, 1, 391-403. Lewry, J. F. and K. D. Collerson, 1990. The Trans-Hudson Orogen: extent, subdivision, and problems, in: J. F. Lewry and M . R. Stauffer (eds.), The early Proterozoic Trans-Hudson Orogen of North America, Geol. Assoc. Can. Spec. Pap., 37, 1-14. Lewry, J. F., Z. Hajnal, A. Green, S. B. Lucas, D. White, M . R. Stauffer, K. E. Ashton, W. Weber and R. Clowes, 1994. Structure of a Paleoproterozoic continent-continent collision zone: a LITHOPROBE seismic reflection profile across the Trans-Hudson Orogen, Canada, Tectonophysics, 232, 143-160. Lewry, J., S. Lucas, R. Stem, K. Ansdell and K. Ashton, 1996. Tectonic assembly and oro-genic closure in the Trans-Hudson Orogen, Geological Association of Canada Annual Meeting, Abstracts, p. A56. Lucas, S. B., A. Green, Z. Hajnal, D. White, J. Lewry, K. Ashton, W. Weber and R. Clowes, 1993. Deep seismic profile across a Proterozoic collision zone: surprises at depth, Nature, 363, 339-342. REFERENCES 78 Majorowicz, J. A., F. W. Jones and A. M . Jessop, 1986. Geothermics of the Williston basin in Canada in relation to hydrodynamics and hydrocarbon occurences, Geophys-ics, 51, 767-779. Morel-a-l'Hussier,P., A. G. Green and C. J. Pike, 1987. Crustal refraction surveys across the Trans-Hudson Orogen - Williston Basin of south central Canada, J. Geophys. Res., 92, 6403-6420. Morgan, P., 1984. The thermal structure and thermal evolution of the continental lithos-phere, in: Pollack, H. and V. R. Murthy (eds.) Structure and evolution of the continental lithosphere, Phys. Chem. Earth, 15, 107-193. Morgan, P., 1995. Diamond exploration from the bottom up: regional geophysical signa-tures of lithospheric conditions favorable for diamond exploration, in: W. L. Griffin (ed.): Diamond Exploration into the 21st Century, J. Geochem. Explor., 53, 145-165. Nelson, K. D., D. J. Baird, J. J. Walters, M . Hauck, L. D. Brown, J. E. Oliver, J. L. Ahem, Z. Hajnal, A. G. Jones and L. L. Sloss, 1993. Trans-Hudson orogen and Williston basin in Montana and North-Dakota: New COCORP deep-profiling results, Geology, 21, 447-450. Nemeth, B. and Z. Hajnal, 1996. Structure of the lithospheric mantle beneath the Trans-Hudson Orogen, Saskatchewan, Report of Sixth Trans-Hudson Orogen Tran-sect Meeting, Lithoprobe Report, 55, 190-196. Nemeth,B., Z. Hajnal and S. B. Lucas, 1996. Moho signature from wide-angle reflec-tions: Preliminary results of the 1993 Trans-Hudson Orogen refraction experiment, Tectonophysics, 264, 111-122. Nolet, G. and A. Zielhuis, 1994. Low 5 velocities under the Tornquist-Teisseyre zone: Evidence for water injection into the transition zone by subduction, J. Geophys. Res., 99, 15813-15820. Nyblade, A. A . and H. N . Pollack, 1993. A global analysis of heat flow from Precam-brian terrains: implications for the thermal structure of Archean and Proterozoic lithosphere, J. Geophys. Res., 98, 12207-12218. O'Reilly, S. Y . and W. L. Griffin, 1996. 4-D Lithosphere Mapping: methodology and examples, Tectonophysics, 262, 3-18. Oldenburg, D. W., 1994. Lecture notes in Linear Inverse Theory, University of British Columbia, unpublished. Olson, P. and H. Singer, 1985. Creeping plumes, J. Fluid Mech., 158, 511-531. Polet, J. and D. L. Anderson, 1995. Depth extent of cratons as inferred from tomographic studies, Geology, 23, 205-208. REFERENCES 79 Pulliam, R. J., D. W. Vasco and L. R. Johnson, 1993. Tomographic inversions for mantle P wave velocity structure based on the minimization of P and I1 norms of Interna-tional Seismological Centre travel time residuals, J. Geophys. Res., 98, 699-734. Silver, P. G. and W. W. Chan, 1991. Shear wave splitting and subcontinental mantle deformation, J. Geophys. Res., 96, 16429-16454. Silver, P. G. and S. Kaneshima, 1993. Constraints on mantle anisotropy beneath Pre-cambrian North America from a transportable teleseismic experiment, Geophys. Res. Lett, 20, 1127-1130. Sobolev, S. V., H. Zeyen, G. Stoll, F. Werling, R. Altherr and K. Fuchs, 1996. Upper mantle temperatures from teleseismic tomography of French Massif Central including effects of composition , mineral reactions, anharmonicity, anelasticity and partial melt, Earth and Planet. Sci. Lett, 139, 147-163. Stacey, F. D., 1992. Physics of the Earth, Brookfield Press, Kenmore, 3rd. ed., 513p. Strnad, J. G., 1993. Diamondiferous kimberlites in Saskatchewan, Canada: global, con-tinental, regional and local setting, in: K. P. E. Dunne and B. Grant (eds.), Mid-continent Diamonds, GAG-MAC Symposium Volume - Edmonton, Alberta - May 17-18, 1993, Mineral Deposit Research Unit, Dept. of Geological Sciences, The Uni-versity of British Columbia, Vancouver, 11-17. Swanson, F. J. and M . R. Gent, 1993. Results of reconnaisance diamond indicator mineral sampling, Saskatchewan, in: K. P. E. Dunne and B. Grant (eds.), Mid-continent Diamonds, GAC-MAC Symposium Volume - Edmonton, Alberta - May 17-18, 1993, Mineral Deposit Research Unit, Dept. of Geological Sciences, The University of British Columbia, Vancouver, 113-119. Tackley, P. J., 1996. Spontaneous melt-driven asthenospheric instabilities: an alternative to plumes? EOS, 77(46 Suppl.), F768-F769. Tackley, P. J. and D. J. Stevenson, 1993. A mechanism for spontaneous self-perpetuating volcanism on the terrestrial planets, in: D. B. Stone and S. K. Runcorn (eds.), Flow and Creep in the Solar System: Observations, Modeling and Theory, Kluwer, Amsterdam, 307-321. Taylor, S. R. and M . N. Toksoz, 1979. Three-dimensional crust and upper mantle struc-ture of the northeastern United States, J. Geophys. Res., 84, 7627-7644. Turcotte, D. L. and G. Schubert, 1982. Geodnamics applications of continuum physics to geological problems, Wiley, New York, 450 p. van der Lee, S. and G. Nolet, 1996. The upper-mantle S-velocity structure of North America, J. Geophys. Res., submitted. VanDecar, J. C. and R. S. Crosson, 1990. Determination of teleseismic relative phase arrival times using multi-channel cross-correlation and least squares, Bull. Seis. Soc. Am., 80, 150-169. REFERENCES 80 VanDecar, J. C , D. E. James and M . Assumpgao, 1995. Seismic evidence for coherent flow of the crust and upper mantle below South America since the breakup of Gondwana, Nature, 378, 25-31. VanDecar, J. C., 1991. Upper-mantle Structure of the Cascadia Subduction Zone from Non-linear Teleseismic Travel-time Inversion, PhD. thesis, University of Washington, Seattle. Wessel, P., and W. H. F. Smith, 1995.. The Generic Mapping Tools (GMT) version 3.0 Technical Reference & Cookbook, SOEST/NOAA. Wickens, A. J. and G. G. R. Buchbinder, 1980. 5-wave residuals in Canada, Bull. Seis. Soc. Am., 70, 809-822. Wiggins, R. A., 1976. Interpolation of digitized curves, Bull. Seis. Soc. Am., 66, 2077-2081. Wysession, M . E., K. M . Fischer, T. J. Clarke, G. I. Al-eqabi, M . J. Fouch, P. J. Shore, R. W. Valenzuela, A. Li and J. M . Zaslow, 1996. Slicing into the Earth, EOS, 77, 477-482. Young, J. B., B. W. Presgrave, H. Aichele, D. A. Wiens and E. A. Flinn, 1996. The Flinn-Engdahl Regionalisation Scheme: the 1995 revision, Phys. Earth Plan. Int., 96, 223-297. Zhang, Y.-S. and T. Tanimoto, 1993. High-resolution global upper mantle structure and plate tectonics, J. Geophys. Res., 98, 9793-9823. Zielhuis, A. and R. D. van der Hilst, 1996. Upper-mantle shear velocity beneath eastern Australia from inversion of waveforms from SKIPPY portable arrays, Geophys. J. Int., 127, 1-16. 81 A Station information 1991-1992 station latitude longitude elev. area operational equipment footing [m] from-to date DLK E4.9164 -103.3760 374 Deschambeault Lake 910806-920109 BB PRS4 Guralp buried FLX SI.3800 -109.6812 704 Flaxcomb 910703-920104 BB PRS4 Guralp buried GRY SO.7088 -102.5736 605 Grayson 910910-920107 SP PRS4 Hark buried HLB S3.3022 -106.1309 506 Holbein 910917-910110 SP PRS4 Mark basement RAY 51.2969 -104.4512 623 Raymore 910902-920106 BB PRS4 Guralp floor SID 53.5777 -104.6771 472 Snosden 910913-920111 SP PRS4 Hark buried TSD 52.9088 -103.9498 462 Tisdale 910703-920112 BB PRS4 Guralp fl o o r USK 52.1961 -106.3982 605 Bergheim 910702-920112 BB PRS4 Guralp vault 1994-1996 station latitude longitude BIR 53.1622 -105.3931 CHO 53.4254 -104.584S DAV 51.1112 -106.0266 DLK 54.9165 -103.3766 DLK 54.9165 -103.3766 ERV 52.8167 -102.1987 ERW 52.8167 -102.1987 FFC 54.7250 -101.9783 GRY SO.7087 -102.5742 HUD 52.7404 -102.2721 LAN 51.8943 -104.8977 LAB 51.8943 -104.8977 LBL 54.3029 -104.6778 LBL 54.3029 -104.6778 LBL 54.3029 -104.6778 HID 52.4726 -105.4473 NEV 50.0161 -107.5965 BOR 51.8160 -102.2167 qui 52.2324 -104.1041 SEH 51.3298 -104.6443 SHE 53.3204 -106.3186 SOU 50.8475 -104.3462 TQB 53.5684 -103.6199 TOB 53.5684 -103.6199 TSD 52.9065 -103.9173 TSD 52.9065 -103.9173 USK 52.1962 -106.3982 USK 52.1962 -106.3982 WAS 53.8633 -106.0952 WAS 53.8633 -106.0952 HAS 53.8633 -106.0952 WEL 53.0647 -105.2950 elev. area Cm] 463 Birch H i l l s 364 Choiceland 705 Davidson 374 Deschambeault La 374 368 Ervood 368 338 F l i n Flon 562 Grayson 421 Hudson Bay 545 Lanigan 545 it 570 L i t t l e Bear Lake 570 M 570 M 523 Middle Lake 875 Seville 496 Borquay 526 Q u i l l Lake 521 Semans 460 Shellbrook 585 Southey 305 Tobin Lake 305 it 448 Tisdale 448 Tisdale 605 Bergheim 605 n 671 Lake Haskesiu 671 •i 671 M 438 Veldon operational equipment from-to date 950619--960603 SP PRS4 951210--960608 SP PRS4 941206--960601 SP PRS4 941015--950824 BB Hicro 950911--960610 BB BARS 950317--950707 BB Hicro 950707--960728 BB BARS sinee 930828 BB IRIS 941126--960602 SP PRS4 941105--950319 BB Hicro 950323--960612 BB HITS 960612--960727 SP PRS4 941010--950705 BB Hicro 950324--950705 BB PRS4 950705--960402 BB BARS 951203--960603 SP PRS4 941122--951203 SP PRS4 941111--960611 SP PRS4 941114--951209 SP PRS4 951209--960602 SP PRS4 941031--960603 SP PRS4 941116--951210 SP PRS4 941013-950706 BB Hicro 950706--960728 BB BARS 950518--960611 BB HITS 960611--960827 BB BARS 941007--950321 SP PRS4 941028--960730 BB HITS 941029--950911 BB Hicro 950704--960610 BB BARS 960610 -960726 SP PRS4 941029--950619 SP PRS4 footing Hillmore basement Hillmore basement Hark floor/buried Guralp buried Guralp " Guralp basement Guralp " SE vault Hark buried Guralp basement Guralp basement Hark ii Guralp buried Guralp ii Guralp II Hark flo o r Hark basement Hillmore fl o o r Hark flo o r Hark f l o o r Hark basement Hillmore fl o o r Guralp basement Guralp ii Guralp basement Guralp ii Hark (z) vault Guralp II Guralp buried Guralp II Hark II Hillmore basement SP : short period seismometers - Hark Products L-4C (eigenfrequency lsec) Hillmore Hark II (eigenfrequency 2sec) BB : broad band seismometers - Guralp CHG-3 (eigenfrequency 20sec) SE: Streckeisen STS-1/VBB (eigenfrequency SOsec) PRS4 : EDA PRS-4 data logger Hicro: Geotech HCR-600 recorder HITS : Hits designed Anglo American data logger BARS : BARS CSD20 data logger IRIS : station of the Incorporated Research Institutions for Seismology vault: concrete pier in underground vault buried: buried vault basement: concrete floor of basement floor: concrete floor 82 B Eventlist This appendix contains a list of all 331 events used in the tomographic inversion. The list was compiled using the Preliminary Determination of Epicenters (PDE) weekly listings of the National Earthquake Information Service (NEIS) of the United States Geological Survey. The geographical regions were converted from numbers to names using Young et al. (1996). A couple of events did not appear in the PDE listing but were found in the near-real-time Earthquake Bulletin ("quake") provided by the NEIS. Locations for these events have been kept in capital letters to distinguish between the two source listings. date location depth mag= geographical region [km] nitude 91/08/13 91/08/14 91/08/17 91/08/26 91/09/15 91/09/18 91/09/21 91/09/30 91/10/06 91/10/08 91/10/19 91/10/19 91/10/20 91/10/25 91/10/26 91/10/27 91/10/30 91/11/10 91/11/11 91/11/11 91/11/13 91/11/13 91/11/13 91/11/16 91/11/19 91/11/23 91/11/24 91/11/25 91/11/26 91/11/26 91/11/27 91/11/27 91/12/07 91/12/08 91/12/11 91/12/13 91/12/13 91/12/13 91/12/13 91/12/14 91/12/15 91/12/15 91/12/17 91/12/19 91/12/22 91/12/23 94/11/09 94/11/18 94/11/22 94/12/05 94/12/06 94/12/07 94/12/08 02:11: 12:53 06:18 14:59 06:39 09:48 15:19 00:21 16:48 03:31 08:48 21:23 01:17 10:39 02:27 16:17 10:35 08:20 17:46 22:34 13:21 19:20 20:27 11:12 22:28 20:00 03:47 14:15 10:41 19:40 05:03 11:45 16:34 14:09 17:03 02:33 05:45 08:00 18:59 00:07 18:56 21:47 06:38 01:33 08:43 13:10 18:21 16:47 11:11 07:09 09:06 03:37 03:29 :16.3 :28.3 :20.1 :46.6 :12.5 :13.0 :50.6 :47.5 :21.2 :15.6 :18.6 :15.5 :02.1 :02.5 :31.5 :29.8 :46.8 :20.0 :00.6 :40.6 :47.6 :04.2 :35.7 -.38.2 :50.8 :41.3 :10.0 :44.7 :35.8 :48.2 :32.2 :19.1 :55.3 :54.0 :09.3 :52.2 :29.9 :08,4 :11.0 :57.2 05.2 00.2 :16.2 42.2 13.3 :04.8 -.03.4 :25.2 :55.7 :06.5 :08.8 :55.S :12.0 52.8071 34 54.3345 169 10.0521 69 42.0251 144 17.804S 115 14.5831 91 16.162S 173 20.973S 178 7.349S 74 45.5571 149 8.987S 117 30.7381 78 63.7081 167 43.2121 144 18.5811 145 22.356S 67 15.183S 173 9.9731 83 17.9221 105 24.7101 142 2.8981 76 33.4561 137 3.1681 78 5.335S 100 4.5571 77 26.804S 114 16.6001 97 8.766S 74 51.8261 176 42.0191 142 48.2321 154 4.740H 77 45.5091 151 45.4951 151 17.819S 116 45.5821 151 45.538B 151 45.3391 151 45.5541 161 23.108S 66 17.512S 70 45.2531 151 47.3471 151 46.1721 151 45.4751 161 45.8731 151 43.5231 147 43.5721 147 43.9121 147 48.0341 146 15.332S 75 23.458S 66 52.9851 171 .977V .343V .862V • 698E .933V ,042V .128V .608V .855V .039E .209E .792E .038V .386E .540E .591V .258V .436V .608V .549E .552V .870E .564V .726V .485V .899V .790V .422V .120V .555E .864E .626V .605E .596E .028V .629E .531E .807E .718E .555V .541V .954E .525E .198E .048E .947E .189E .489E .295E .824E .359V .736V .105E 10 299 10 45 10 5 33 680 144 147 90 19 33 111 193 178 58 33 33 23 148 308 29 10 21 10 30 147 69 55 35 39 45 33 10 32 33 33 52 211 102 47 149 43 26 23 61 33 33 450 40 243 33 5.3Mb North Atlantic Ocean 5.7Mb Fox Islands, Aleutian Islands, United States 5.1Mb Venezuela 5.8Mb Hokkaido, Japan, region 5.8Mb Southern East Pacific Rise 5.7Hb Guatemala 5.9Mb Tonga Islands 6.0Hb F i j i Islands region 5.4Mb Peru-Brazil border region 6.0Hb Kuril Islands, Russia 5.7Mb Sunbava, Indonesia, region 6.5Hb Horthern India 5.0Hb Fox Islands, Aleutian Islands, United States 5.7Hb Hokkaido, Japan, region 5.8Mb Mariana Islands 5.8Mb Chile-Bolivia border region 5.7Hb Tonga Islands 5.1Mb Costa Rica 5.5Hb Off coast of Jalisco, Mexico 5.9Mb Volcano Islands, Japan, region 5.0Mb Colombia 5.4Mb Hear south coast of eastern Honshu, Japan 5.3Mb South of Panama 5.3Mb Central East Pacific Rise 6.5Mb Bear vest coast of Colombia 5.2Mb Easter Island region 5.1Mb Oaxaca, Mexico 5.3Mb Peru-Brazil border region 5.8Mb Andreanof Islands, Aleutian Islands, United States 6.2Mb Hokkaido, Japan, region 5.9Mb Kuril Islands, Russia 5.1Hb Hear vest coast of Colombia 5.4Mb Kuril Islands, Russia 5.3Mb Kuril Islands, Russia 5.9Mb Southern East Pacific Rise 6.2Mb Kuril Islands, Russia 6.0Mb Kuril Islands, Russia 5.5Mb Kuril Islands, Russia 6.2Mb Kuril Islands, Russia 5.0Mb Jujuy Province, Argentina 5.6Mb Hear coast of Peru 5.1Mb Kuril Islands, Russia 5.9Mb Kuril Islands, Russia 6.1Mb Kuril Islands, Russia 6.3Mb Kuril Islands, Russia 5.9Mb Kuril Islands, Russia 6.1Hb Kuril Islands, Russia 5.4Mb Ku r i l Islands, Russia 5.7Mb Kuril Islands, Russia 4.7Mb Sea of Okhotsk 5.3Mb Hear coast of Peru 5.8Mb Jujuy Province, Argentina 5.0Mb Hear Islands, Aleutian Islands, United States B EVENTLIST 94/12/10 94/12/10 94/12/10 94/12/12 94/12/13 94/12/14 94/12/18 94/12/18 94/12/2E 94/12/26 94/12/27 94/12/28 94/12/30 94/12/31 95/01/04 95/01/06 95/01/07 96/01/08 95/01/10 95/01/11 95/01/11 96/01/12 95/01/13 95/01/15 96/01/16 95/01/19 95/01/19 95/01/20 95/01/20 95/01/21 95/01/22 95/01/29 95/02/02 95/02/04 95/02/05 95/02/06 95/02/08 95/02/10 95/02/11 95/02/12 95/02/14 95/02/14 95/02/14 95/02/16 95/02/18 95/02/19 95/02/21 95/02/22 95/02/23 95/02/23 95/02/23 95/02/23 95/02/24 95/02/28 95/03/01 95/03/03 95/03/03 96/03/06 95/03/08 95/03/08 95/03/10 96/03/10 95/03/11 95/03/14 95/03/14 95/03/18 95/03/18 95/03/19 95/03/19 96/03/20 95/03/21 96/03/23 95/03/23 95/03/25 95/03/26 95/03/26 95/03/26 95/03/30 95/03/31 95/04/01 95/04/01 95/04/01 95/04/04 95/04/07 95/04/08 03:39: 15:24 16:17: 07:41 08:23: 07:28 11:07 20:38: 02:39: 03:08: 20:42: 12:19: 15:12: 13:50: 02:47: 22:37: 02:13: 09:22: 13:18: 07:48: 10:26: 10:26: 03:12: 02:40: 20:46: 03:00: 15:05: 03:35: 13:59: 08:47: 10:41: 04:53: 12:53: 19:07: 20:37: 13:51: 18:40: 20:27: 22:45: 01:02: 07:35: 15:53: 20:47: 09:36: 13:29: 04:46: 02:09: 10:43: 06:19: 20:06: 21:03: 21:40: 20:09: 05:03: 02:04: 12:02: 21:12 18:44 03:45 07:10 05:22 08:22 16:21 12:56 17:33 09:27 18:02 11:43 15:59 14:49 09:10 02:08 17:51 10:42 05:21 06:57 15:12 22:15 14:01 03:49 05:60 17:54 07:10 22:06 17:45 31.4 18.6 41.0 56.4 48.3 55.8 34.7 32.6 27.8 17.7 59.3 23.6 26.9 22.6 27.1 37.9 :28.2 :19.0 49.9 :20.4 25.3 46.9 :59.8 :18.9 51.1 23.1 03.6 46.1 20.2 29.9 27.7 37.6 53.1 04.4 10.9 35.6 25.1 03.3 33.2 07.0 22.1 56.9 41.1 17.9 06.8 50.3 50.9 24.6 02.3 06.7 02.0 31.0 07.8 03.2 23.6 10.7 :37.6 :01.3 :59.7 :35.5 :22.4 :37.8 :10.7 :30.6 :50.7 :19.2 -.37.0 30.2 :19.8 18.1 :35.7 34.5 :07.9 :18.7 :19.0 :11.6 :31.7 :52.0 :40.8 :33.4 20.6 02.3 15.3 58.0 18.1 23.579S 6.8251 18.2431 17.504S 16.3411 9.461S 37.3171 17.858S 11.1771 53.6881 44.8661 40.451B 18.595B 40.186B 51.5901 40.227B 1.529S 8.444S 5.429B 41.965B 7.960S 43.986B 43.102B 27.619B 34.5491 43.326B 5.075B 43.261B 5.2141 43.335B 5.1261 29.234B 10.755B 19.559S 6.849B 41.124B 4.162B 19.770S 12.607B 5.796S 46.158B 23.290S 43.991B E.778S 46.6671 43.098B 45.942B 52.7331 24.136B 18.992B 35.0391 35.006B 27.185B 11.5041 55.662B 6.472S 14.617S 15.572B 16.555B 22.144S 46.075B 11.538S 44.008B 15.168S 54.793B 29.279B 42.502B 43.824B 8.675S 13.309B 22.799S 7.405B 15.518B 42.9891 54.8211 54.8561 2.063S 44.821B 38.1SOB 37.924B 52.279B 18.365B 33.713B 15.187S 21.804B 70.522V 72.969V 101.354V 69.650V 98.410V 1S9.308E 139.855E 178.690V 86.592V 164.480V 149.558E 143.491E 14S.266E 142.658E 173.661V 142.242E 78.014V 74.289V 82.630V 142.460E 73.946V 147.088E 147.070E 128.443E 135.002E 147.010E 72.918V 146.821E 72.893V 146.717E 72.966V 141.202E 42.557V 70.407V 82.675V 142.188E 76.644V 68.544V 81.603V 76.136V 146.822E 67.702V 148.098E 76.130V 146.894E 146.842E 151.553E 171.372E 121.600E 146.101E 32.266E 32.270E 127.389E 85.936V 161.273E 154.987E 176.627V 92.016V 59.S74V 68.369V 143.540E 130.829V 148.132E 64.794V 161.295V 140.691E 87.179E 147.910E 108.959V 89.243V 68.506V 76.696V 46.734V 147.541E 161.413V 161.363V 79.497V 137.572E 135.068E 139.177E 159.133E 146.802E 38.593V 173.S94V 142.632E 37 5.7Hb 158 5.0Hb 67 6.5Hb 151 5.8Kb 17 S.3Hb 33 5.7Mb 25 5.2Mb 551 5.6Mb 33 4.8Mb 33 5.3Hb 33 5.6Hb 33 6.4Mb 235 5.6Kb 33 5.8Kb 33 5.4Mb 57 6.7Hb 166 5.4Mb 149 5.1Hb 10 4.6Kb 33 5.3Mb 174 6.3Mb 33 6.2Mb 33 5.8Mb 47 5.7Mb 16 6.4Mb 40 B.BHb 18 6.4Mb 61 5.7Mb 33 5.2Mb 63 6.5Hb 22 5.6Kb 67 5.6Kb 10 5.6Mb 108 5.0Mb 10 5.8Mb 71 5.6Mb 69 6.3Mb 164 5.4Mb 11 5.3Hb 22 5.7Mb 400 4.9Kb 156 5.7Mb 37 5.9Mb 18 5.3Mb 355 5.6Mb 33 5.8Mb 33 6.8Mb 33 4.8Mb 44 5.8Mb 121 5.3Mb 15 5.8Mb 10 5.2Mb 92 5.4Mb 131 5.1Mb 92 5.3Kb 35 5.8Mb 22 5.6Mb 174 5.0Mb 15 6.3Mb 114 4.9Mb 350 5.4Mb 10 5.3Hb 33 6.0Mb 588 5.1Mb 33 6.1Hb 104 5.5Hb 21 5.2Mb 49 5.3Kb 10 5.0Mb 82 4.9Mb 103 5.0Mb 10 5.1Hb 33 5.1Mb 53 5.2Mb 33 5.8Mb 33 5.3Kb 91 5.1Mb 319 5.4Mb 365 6.0Mb 11 5.8Mb 47 5.9Mb 57 5.2Mb 10 5.2Mb 31 6.7Mb 319 6.3Mb Bear coast of northern Chile florthern Colombia Guerrero, Mexico Peru-Bolivia border region Bear coast of Guerrero, Mexico Bougainville - Solomon Islands region Eastern Honshu, Japan F i j i Islands region Bear coast of Bicaragua Unimak Islands, Alaska, United States, region Kuril Islands, Russia Off east coast of Honshu, Japan Mariana Islands Bear east coast of eastern Honshu, Japan Andreanof Islands, Aleutian Islands, United State Bear east coast of eastern Honshu, Japan Ecuador Peru-Brazil border region South of Panama Hokkaido, Japan, region Peru-Brazil border region Kuril Islands, Russia K u r i l Islands, Russia Ryukyu Islands, Japan Bear south coast of western Honshu, Japan Kuril Islands, Russia Colombia Kuril Islands, Russia Colombia Kuril Islands, Russia Colombia Southeast of Honshu, Japan florthern Mid-Atlantic Ridge Bear coast of northern Chile South of Panama Hokkaido, Japan, region Colombia Chile-Bolivia border region Caribbean Sea florthern Peru florthsest of K u r i l Islands, Russia Chile-Argentina border region East of Ku r i l Islands, Russia florthern Peru Sea of Okhotsk Kuril Islands, Russia Kuril Islands, Russia Bear Islands, Aleutian Islands, United States Taiwan Mariana Islands Cyprus region Cyprus region Ryukyu Islands, Japan Bicaragua Bear east coast of Kamchatka Peninsula,Russia Bougainville - Solomon Islands region Samoa Islands region Mexico-Guatemala border region Leeward Islands florthern Chile Sakhalin Island, Russia South Pacific Ocean Kuril Islands, Russia Central Bolivia Alaska Peninsula, United States Southeast of Honshu, Japan Borthern Xinjiang, China Kuril Islands, Russia Central East Pacific Rise E l Salvador florthern Chile florthern Colombia florthern Mid-Atlantic Ridge Off southeast coast of Hokkaido, Japan Alaska Peninsula, United States Alaska Peninsula, United States Bear coast of Ecuador Eastern Sea of Japan Sea of Japan Eastern Honshu, Japan Off east coast of Kamchatka Peninsula,Russia Mariana Islands florthern Mid-Atlantic Ridge Tonga Islands Mariana Islands region B EVENTLIST 84 95/04/15 95/04/17 95/04/17 96/04/17 95/04/18 95/04/19 95/04/20 95/05/13 95/05/15 95/05/16 95/05/18 95/05/19 95/05/20 95/05/23 95/05/25 95/05/27 95/05/29 95/05/29 95/05/29 95/05/29 95/05/30 95/05/31 95/06/02 95/06/09 95/06/09 95/06/12 96/06/21 95/06/24 95/06/25 95/06/25 95/06/27 95/06/29 95/07/07 96/07/08 95/07/08 95/07/21 95/07/25 95/07/29 95/07/30 95/07/30 95/07/30 96/07/30 95/07/31 95/08/02 95/08/03 95/08/03 95/08/06 95/08/17 95/08/19 95/08/23 95/08/24 95/09/06 95/09/08 96/09/11 95/09/14 96/09/17 95/09/19 95/09/23 95/09/23 95/09/26 96/10/01 95/10/01 95/10/01 95/10/02 95/10/03 95/10/03 95/10/04 96/10/04 95/10/06 95/10/07 95/10/08 95/10/09 95/10/18 95/10/18 95/10/19 95/10/19 95/10/20 95/10/21 95/10/27 95/10/28 95/11/01 95/11/13 95/11/13 95/11/14 95/11/22 05:36: 07:14: 12:02: 23:28: 05:23: 03:50: 20:49: 08:47: 04:05: 03:35: 00:06: 09:31: 22:18: 10:01 04:59 13:03 04:58 07:29 09:06: 10:21: 04:12: 16:08: 19:07: 05:35: 08:10: 03:35: 16:33: 06:58: 06:59: 12:25: 10:09: 07:45: 21:15: 05:42: 17:16: 22:44: 15:13: 16:18: 05:11: 05:47: 10:35: 21:05: 08:48: 00:14: 01:57: 08:18: 11:59: 00:59: 21:43: 07:06: 01:55: 22:48: 17:25: 04:22: 14:04: 07:25: 03:31: 20:56: 22:31: 07:14: 12:60: 15:67: 17:06 01:35 01:51 12:45 09:17 15:12 11:39 21:28 08:55 15:35 09:30 10:37 00:32 02:41 19:21 02:38 21:59 14:38 00:35 07:38 08:43 15:14 04:15 32.3 35.2 30.6 08.3 58.7 05.6 10.3 12.6 58.0 03.0 26.7 46.0 57.9 30.6 51.5 55.6 32.4 46.2 22.7 34.2 43.9 40.0 22.6 49.3 46.9 48.8 05.9 :06.5 :05.6 37.8 58.0 :09.7 18.5 56.9 28.1 07.6 26.7 44.7 23.5 02.1 42.5 60.6 34.7 09.6 21.7 -.53.5 34.6 57.8 32.4 02.6 34.6 52.6 49.5 52.8 31.6 30.7 53.6 03.7 58.3 37.7 15.6 16.0 :03.0 :48.4 :24.1 :00.1 :30.2 :34.9 :36.0 :06.1 :49.9 :55.7 :38.6 :26.3 :06.4 :37.9 :28.6 :57.5 :57.0 :33.0 :32.3 :45.9 :14.7 :03.5 :11.7 88. 70. 21. 66. 73. 43.8881 147. 33.778B 38. 11.3231 86. 45.904B 151. 45.806H 151. 44.0271 148. 45.901H 151. 40.1441 21. 41.6651 36.4851 0.950S 23.811S 5.641B 43.634B 141 43.912B 147. 52.563B 142 35.056B 32 10.144S 163. 10.176B 104 52.681B 142. 29.468B 138. 18.949B 107 31.760S 71. 21.664S 68. 43.957B 147. 8.308S 75. 11.546S 77. 3.979S 153. 24.599B 121. 26.1311 124. 18.8471 81. 48.784B 154. 33.947B 137. 39.640B 143. 53.6461 163. 36.443B 103. 10.694B 41. 30.365B 138. 23.364S 70. 23.290S 24.392S 23.317S 10.396S 23.152S 23.132S 28.346S 44.370B 147 41.587B 88 5.096B 75 18.857B 145 18.920B 144 14.9821 94 14.938B 94 0.998B 101 16.844B 98 35.557S 21.228S 24.345S 128 10.529S 78 41.795B 143 31.377S 71 38.099B 30 29.314B 138 66.739B 179 2.774S 77 2.818S 77 75.983B 6 52.7131 172 20.021S 175 2.783S 77 40.9891 72 19.147B 104 36.441B 70 27.920B 130 28.145B 130 28.075B 130 18.775B 145 16.890B 93 21.950S 139 6.330S 154 28.968S 71 15.1O0S 173 56.067B 114 18.8011 145 28.818B 34 70 70 70 78 70 70 69 74 68 299E 600V 616V 288E 484E 204E 253E 684E 821E 896E 985V 496V 818V 783E 374E 814E 264E 734E 007V 825E 425E 434V 297V 029V 450E 913V 564V 945E 713E 760E 730V 459E 122E 353E 546V 105E 211V 292E 312V 325V 691V 590V 316V 578V 602V 198V 263E 782E 690V 186E 951E 251V 216V 441V 599V 052V 740V 143V 897V 387E 266V 176E 954E 223E 884V 897V 869E 463E 922V 840V 118E 220V .428E .337E .206E .309E .391E .451V .190V .330E .503V .602V .477E .195E .861E 48 5.4Kb 10 5.8Kb 32 4.9Hb 34 6.1Kb 33 5.7Hb 33 5.9Mb 33 5.7Hb 13 6.2Kb 0 6.1Mb 190 5.7Hb 10 6.2Kb 230 5.0Mb 142 5.0Mb 33 5.4Mb 77 5.5Mb 33 6.8Mb 10 5.3Mb 33 5.8Mb 33 5.0Hb 33 5.3Hb 468 4.9Mb 33 5.4Mb 70 5.4Mb 132 5.2Mb 95 5.1Mb 34 5.7Mb 70 5.5Mb 386 6.2Mb 47 5.8Mb 171 5.1Mb 10 5.7Mb 62 5.9Mb 324 5.8Mb 40 5.7Mb 33 5.8Mb 33 5.7Mb 10 5.4Mb 436 5.5Mb 47 6.6Mb 33 5.8Mb 33 5.8Mb 33 5.8Mb 93 5.8Mb 33 5.4Mb 33 5.4Mb 104 5.9Mb 74 5.3Mb 0 6.0Mb 125 6.1Mb 596 6.3Mb 589 5.9Mb 33 5.8Mb 33 5.5Mb 10 5.3Mb 21 6.4Mb 33 5.8Mb 110 5.7Mb 10 5.3Mb 73 5.9Mb 36 5.7Mb 64 5.4Mb 33 5.7Mb 427 5.5Mb 10 5.1Mb 27 6.5Mb 33 6.0Hb 10 5.1Hb 33 5.6Mb 209 5.7Mb 33 5.8Mb 48 5.9Mb 49 6.8Mb 226 5.4Mb 27 6.5Mb 33 5.9Mb 31 6.3Mb 225 5.3Mb 161 6.2Mb 0 5.6Mb 33 5.8Mb 20 6.3Mb 33 5.7Mb 24 5.9Mb 600 5.0Hb 10 6.2Mb Kuril Islands, Russia northern Mid-Atlantic Ridge Bear coast of Bicaragua Kuril Islands, Russia Kuril Islands, Russia Kuril Islands, Russia Kuril Islands, Russia Greece Southern Xinjiang, China Hindu Kush, Afghanistan, region Central Mid-Atlantic Ridge Jujuy Province, Argentina Colombia Hokkaido, Japan, region Kuril Islands, Russia Sakhalin Island, Russia Cyprus region Bougainville - Solomon Islands region Off coast of Mexico Sakhalin Island, Russia Southeast of Honshu, Japan Off coast of Jalisco, Mexico Bear coast of central Chile Chile-Bolivia border region Kuril Islands, Russia Central Peru Bear coast of Peru les Ireland, Papua lew Guinea, region Taiwan Hortheast of Taiwan Caribbean Sea Kuril Islands, Russia Bear south coast of eastern Honshu, Japan Off east coast of Honshu, Japan Unimak Islands, Alaska, United States, region Gansu, China Borthern Mid-Atlantic Ridge Southeast of Honshu, Japan Bear coast of northern Chile Hear coast of northern Chile Bear coast of northern Chile Hear coast of northern Chile Bear coast of Peru Hear coast of northern Chile Bear coast of northern Chile Chile-Argentina border region Kuril Islands, Russia Southern Xinjiang, China Colombia Mariana Islands Mariana Islands Off coast of Chiapas, Mexico Off coast of Chiapas, Mexico Eastcentral Pacific Ocean Bear coast of Guerrero, Mexico Off coast of central Chile Chile-Bolivia border region South Pacific Ocean Bear coast of Peru Hokkaido, Japan, region Bear coast of central Chile Turkey Southeast of Honshu, Japan Eastern Siberia, Russia Peru-Ecuador border region Peru-Ecuador border region Greenland Sea Bear Islands, Aleutian Islands, United States Tonga Islands Peru-Ecuador border region Kyrgyzstan Bear coast of Jalisco, Mexico Hindu Kush, Afghanistan, region Ryukyu Islands, Japan Ryukyu Islands, Japan Ryukyu Islands, Japan Mariana Islands Chiapas, Mexico TUAMOTU ARCHIPELAGO REGIOH S0L0H0H ISLABDS Bear coast of central Chile Tonga Islands East of Lake Baykal, Russia Mariana Islands Bgyp* B EVENTLIST 95/11/24 95/11/27 95/11/27 95/11/30 95/11/30 96/11/30 96/12/01 95/12/02 95/12/02 96/12/03 95/12/06 95/12/07 95/12/07 95/12/07 95/12/10 95/12/10 95/12/11 95/12/11 95/12/17 95/12/20 95/12/25 95/12/25 95/12/29 95/12/30 95/12/30 95/12/30 95/12/31 95/12/31 96/01/01 96/01/02 96/01/02 96/01/08 96/01/08 96/01/11 96/01/14 96/01/19 96/01/22 96/01/22 96/01/25 96/01/27 96/01/29 96/01/30 96/01/30 96/01/31 96/01/31 96/02/01 96/02/01 96/02/02 96/02/07 96/02/11 96/02/12 96/02/14 96/02/14 96/02/16 96/02/16 96/02/16 96/02/17 96/02/18 96/02/19 96/02/19 96/02/21 96/02/22 96/02/22 96/02/22 96/02/22 96/02/25 96/02/25 96/02/25 96/02/25 96/02/26 96/02/26 96/02/29 96/03/01 96/03/01 96/03/03 96/03/03 96/03/04 96/03/06 96/03/05 96/03/06 96/03/07 96/03/09 96/03/10 96/03/12 96/03/13 17:24: 05:24: 15:52: 13:15: 16:09: 23:37: 05:20: 17:13: 19:40: 21:38: 23:17: 03:22: 05:12: 19:30: 03:27: 22:23: 14:09: 19:44: 23:48: 11:39: 03:19: 18:18: 13:01: 03:26: 12:11: 16:15: 07:26: 15:12: 09:57: 06:41: 09:04: 10:04: 11:52: 07:19: 06:28: 19:01: 00:10: 13:14: 12:45: 21:29: 10:06: 02:28: 21:14: 19:21: 20:30: 07:18: 17:57: 18:36: 21:36: 09:28: 02:58: 20:31: 21:26: 09:44: 11:34: 15:22: 13:25: 23:49: 02:28: 07:10: 12:51: 06:05 08:38 13:40 14:59 03:08 04:17 09:18 14:17 01:37 07:17 16:58 02:27 06:48 14:55 16:37 05:66 14:62 17:32 01:35 08:38 16:15 08:56 18:43 21:04 12.5 26.0 58.0 37.0 23 36 27 21 10 38.0 20.0 03.0 22.0 24.0 49.0 14.7 24.3 11.0 31.2 18.2 44.8 :17.0 :40.9 :08.2 :07.4 :31.7 :13.7 33.3 :51.7 :03.4 :09.7 51.3 09.2 :09.9 20.7 58.3 23.0 56.0 02.0 57.0 56.7 42.5 55.5 27.0 47.0 05.3 55.5 11.3 45.1 49.3 53.4 06.3 :56.2 :58.1 :30.3 :57.8 :37.7 :28.0 :31.8 :10.3 :04.3 :42.2 :36.5 :53.5 :09.7 :18.8 :09.8 :00.1 :21.0 :36.5 :26.5 :43.9 :24.1 :53.1 :11.8 :31.5 :38.6 :28.6 -.10.1 :03.0 :57.5 :36.2 :22.3 :42.8 :17.9 44.542B 149 22.840S 70 44.620B 149 27.700S 70 44.2201 145. 44.1701 149. 10.060B 104. 44.810H 149. 44.320E 149. 44.697B 150. 44.230B 149. 44.360B 149. 44.430B 149. 44.680B 149. 34.780B 24. 44.331B 149. 18.829B 105. 18.6501 105. 62.566B 32. 27.379B 128. 36.432B 70. 17.512B 82. 10.041B 70. 4.708S 104. 40.790B 143. 31.024B 140. 53.877B 160. 10.080B 70. 53.917B 159. 18.917S 69. 54.054B 159. 63.233B 142. 20.439S 68. 45.083B 150. 44.557B 149. 10.424S 78. 44.560B 147. 40. HOB 142. 18.350B 101. 22.300S 138. 14.322B 92. 47.081B 151. 36.156B 135. 44.506B 149. 44.456B 149. 44.863B 146. 37.757B 19. 11.439B 141. 45.321B 149 4S.291B 150 45.147B 150 45.386B 150 29.220B 140 1.524S 15 15.106S 173 37.343B 142 6.160S 154 1.283S 14 1.223S 42.091S 9.620S 20.195S 8.620B 33.693S 45.208B 148 16.204E 97 16.039B 16.1681 13.022B 16.068B 28.647B 4.154B 23.479S 34.198B 11.657B 11.9051 11.6371 24.092B 122 24.033B 122 18.6S6S 174 23.267S 70 43.4251 148 12.969S 69 48.442B 88 16.7041 98 14 75 79 68 83 71 97 97 91 97 36 84 70 26 86 86 86 091E 080H 270E 660V 660E 370E 040V 220E 570E .299E .390E .380E .380E .S70E .HOE .784E .483V • 360V .169V .22SE .260E .872V .089V .783V .382E .143E .614E .050V .496E .201V .307E 732E .904V .025E .010E .772V .910E .180E .940V .940V 841V .929E 329E .417E .372E .293E .860E .631E .909E .471E .326E .360E .383E .222V .587V .474E .737E 274V 227V 085V 568V .921V .116V .706V . SS7E .963V .875V .909V .050V .807V . 181E .167V .358V .248E 856V .772V 978V .215E .241E • 864V .285V .009E .425V .143E .870V 33 6.1Hb 33 5.3Hb 33 6.0Hb 33 5.5Kb 145 6.1Kb 33 6.0Mb 10 6.2Mb 33 6.5Mb 33 5.4Mb 33 5.7Hb 33 5.6Mb 33 5.7Hb 33 5.8Mb 33 5.9Mb 33 5.3Mb 33 5.6Mb 33 5.7Mb 33 5.6Hb 10 5.2Mb 33 5.2Mb 230 5.4Mb 10 5.2Mb 33 5.4Mb 10 5.2Mb 33 5.7Mb 106 5.4Mb 55 6.0Hb 40 5.1Hb 33 5.8Mb 108 5.1Mb 33 4.9Hb 33 5.5Mb 110 4.8Mb 13 5.5Mb 33 5.5Mb 36 5.6Mb 33 6.2Mb 33 6.1Mb 33 5.2Mb 0 5.4Mb 33 4.8Mb 132 5.3Hb 358 5.0Mb 33 5.5Mb 58 5.8Mb 180 5.7Hb 10 5.2Mb 33 S.3Hb 33 6.3Mb 33 5.3Hb 33 5.5Mb 33 5.9Mb 142 5.9Mb 10 6.0Mb 33 5.1Mb 33 6.2Mb 56 5.5Mb 10 6.3Hb 12 5.5Mb 33 5.9Mb 33 5.8Mb 118 4.9Kb 33 5.0Mb 44 5.9Mb 133 6.2Mb 33 5.9Mb 33 5.1Mb 33 5.6Mb 33 5.5Mb 33 5.1Mb 10 4.9Hb 33 5.1Hb 27 5.1Mb 33 4.9Mb 33 5.7Mb 33 5.7Mb 84 4.8Mb 30 6.1Hb 33 5.6Mb 134 5.4Kb 31 5.3Mb 33 5.6Mb 33 5.8Mb 17 5.8Mb 25 5.2Mb Kuril Islands, Russia BEAR COAST OF BORTHERB CHILE KURIL ISLABDS BEAR COAST OF BORTHERB CHILE HOKKAIDO, JAPAB REGIOB KURIL ISLABDS OFF COAST OF MEXICO KURIL ISLABDS KURIL ISLABDS East of Kuril Islands, Russia KURIL ISLABDS KURIL ISLABDS KURIL ISLABDS KURIL ISLABDS CRETE Kuril Islands, Russia Off coast of Jalisco, Mexico OFF COAST OF JALISCO, MEXICO Borth Atlantic Ocean Ryukyu Islands, Japan Hindu Kush, Afghanistan, region Caribbean Sea Venezuela Central East Pacific Rise Off east coast of Honshu, Japan Southeast of Honshu, Japan Hear east coast of Kamchatka Peninsula,Russia Venezuela Bear east coast of Kamchatka Peninsula,Russia Borthern Chile Bear east coast of Kamchatka Peninsula,Russia Sakhalin Island, Russia Chile-Bolivia border region Kuril Islands, Russia Kuril Islands, Russia Bear coast of Peru KURIL ISLABDS BEAR E COAST OF HOHSHU, JAPAB GUERRERO, MEXICO TUAMOTU ARCHIPELAGO REGIOB Bear coast of Chiapas, Mexico Kuril Islands, Russia Sea of Japan Kuril Islands, Russia Kuril Islands, Russia Kuril Islands, Russia Ionian Sea Vestern Caroline Islands, Micronesia Kuril Islands, Russia Kuril Islands, Russia Kuril Islands, Russia Kuril Islands, Russia Southeast of Honshu, Japan Borth of Ascension Island Tonga Islands Off east coast of Honshu, Japan . Bougainville - Solomon Islands region Borth of Ascension Island Borth of Ascension Island Off coast of southern Chile Off coast of northern Peru Chile-Bolivia border region Costa Rica Bear coast of central Chile Kuri l Islands, Russia Oaxaca, Mexico Oaxaca, Mexico Oaxaca, Mexico Bear coast of Guatemala Oaxaca, Mexico Vestern Arabian Peninsula Off coast fo central America Bear coast of northern Chile Crete, Greece Hear coast of Bicaragua Bear coast of Bicaragua Bear coast of Bicaragua Taiwan region Taiwan region Tonga Islands Bear coast of northern Chile East of Kuril Islands, Russia Central Peru Mongolia Bear coast of Guerrero, Mexico B EVENTLIST 86 96/03/18 03:16:16.7 24.027S 67.972V 149 4.9Kb Chile-Argentina border region 96/03/19 17:12:43.0 16.850B 97.306V 33 6.8Kb Hear coast of Oaxaca, Mexico 96/03/20 04:53:27.0 16.960H 97.258V 33 5.3Mb Hear coast of Oaxaca, Mexico 96/03/22 03:24:20.0 51.221H 178.695E 20 5.7Mb Rat Islands, Aleutian Islands, United States 96/03/28 23:03:49.8 1.036S 78.737V 33 5.8Kb Ecuador 96/04/02 07:59:26.1 37.827H 26.995E 33 4.9Mb Dodecanese Islands, Greece 96/04/03 23:00:47.4 14.493H 93.490V 33 5.0Hb Hear coast of Chiapas, Mexico 96/04/07 14:18:59.0 53.260H 159.781E 52 5.1Mb Hear east coast of Kamchatka Peninsula,Russia 96/04/08 16:55:08.4 15.043H 61.609V 181 4.7Hb Leeward Islands 96/04/11 10:51:14.9 56.986H 33.745V 10 4.9Mb Horth Atlantic Ocean 96/04/11 11:24:26.9 10.803S 161.543E 43 5.9Hb Bougainville - Solomon Islands region 96/04/16 00:30:54.6 24.061S 177.036V 111 6.4Mb South of F i j i Islands 96/04/18 17:08:26.1 47.346H 154.171E 33 5.3Mb Kuril Islands, Russia 96/04/19 00:19:31.1 23.944S 70.093V 50 6.0Mb Hear coast of northern Chile 96/04/20 19:17:06.1 24.071S 66.786V 197 5.2Mb Salta Province, Argentina 96/04/21 15:50:04.2 51.437H 178.420V 28 5.0Mb Andreanof Islands, Aleutian Islands, United States 96/04/22 11:27:55.9 29.797H 129.051E 197 5.2Mb Ryukyu Islands, Japan 96/04/23 06:53:35.3 17.280H 101.363V 33 5.3Mb Hear coast of Guerrero, Mexico 96/04/24 17:06:36.3 8.128S 74.362V 151 5.6Mb Peru-Brazil border region 96/04/26 07:01:27.5 36.366H 28.042E 75 5.2Mb Dodecanese Islands, Greece 96/04/27 08:40:41.8 2.368H 79.341V 10 4.8Mb South of Panama 96/04/28 08:18:42.0 52.830H 164.684V 33 5.1Hb South of Alaska 96/05/03 03:32:47.1 40.774H 109.661E 26 5.5Mb Vestern Hei Mongol, China 87 C Timing correction for PRS-4 and WITS WWVB is broadcast from Fort Collins, Colorado. At the start of each second the carrier power is reduced 10 dB. Figure C . l shows an extract from a time trace as recorded by one of the PRS-4 stations. As can be seen, the instrument transforms the drop in power as a rise in amplitude. The pulses are not equal in length, but rather three different pulse widths can be found. These denote reference markers, ones and zeroes by which information is encoded within each minute sequence. Thus each minute has a unique sequence of pulses. How each second is used and an example of a full minute is shown in figure C . l (For minutes: 28 hours: 14 Julian day: 98 UT1 correction: +400 ms 11 a | | | year: 96 '2 g « l 88"~&2°i2 OOTCCN— 8888 © g o o ooTt-ts-H ^ § ' • 0 SIS 0OT^'sl'-H a £ S OOTI-CM-H OOTJ-CJ— WTCN™ >££^-I o t o o » o o o i o o o « o o « o o i o o o o o » o o « i « o o o o o » o » i o « o o o « o o » i o » » o o « o » o I -8 3 10 20 30 40 50 seconds Figure C . l : WWVB time code as recorded by a PRS-4 station. This particular WWVB trace was recorded at station CHO on Easter Sunday 1996. Three different pulse widths can be distinguished. They are labelled above the trace: vertical bar - reference marker, solid circle - 1, open circle - 0. The binary coded decimal system is explained at the top. Decoding the sequence tells us that the first pulse (denoted by an arrow) was broadcast at 14:28 UTC on Julian day 98 (i.e. April 07) in 1996. C TIMING CORRECTION FOR PRS-4 AND WITS 88 Figure C.2: Extract from a WWVB time trace showing digitizing points detailed information on the WWVB time code see Behler & Lombardi, 1991). The first five seconds of the same trace are shown in figure C.2. Superimposed on the trace are the sampling points. A decimation filter in the PRS-4 ensures that the leading edge is «30ms wide. As a consequence, a digitizing frequency of 50 Hz will always sample this edge (see appendix A of Ellis & Hajnal, 1992). A code written by John Amor at UBC was used to obtain the timing corrections. The code works two-fold: a stacking procedure obtains deviations of the internal clock from a full second and an error estimate for this deviation. Cross-correlation of the recorded time trace with a synthetic trace for that minute finds second skips. Only a handful of traces had error uncertainties larger than ±5ms . These were rejected for the tomographic studies. D PRACTICAL NOTES ON VANDECAR'S CODE D Practical notes on VanDecar's code 89 The source code used in this study is that of (VanDecar, 1991) in Fortran 77. The elements that make up the code are shown in figure D . l and notes that were written while working with that code are reproduced on the following pages. The intention of these notes is to aid other students at the Department of Earth and Ocean Sciences who may use this code. seismograms ref. model 1 mccc i 1 1 Q x-files X loc.sta )^ / | procx ref.t* rays sta eqs cals cale rays.total O partials 1 C part-file ^^[grabigl raytrc (jelin.su dtelin Q new.model (jresidualsj) (new.cals, new.cale Figure D . l : Flowchart showing implementation of VanDecar's code Compare with fig. 3.1. Rectangular boxes denote programs, oval boxes denote files. D PRACTICAL NOTES ON VANDECAR'S CODE This i s a brief description of how to run John VanDecar's code for teleseismic inversion. Oftentimes you w i l l see explanations intermixed with unix commands. This i s my say of documenting thoughts and working procedures. The text i s intended as a supplement to the README f i l e s compiled by the gifted programmer John himself. programs are located i n (include in your path) \ivana\vandecar\bin \snorcle\inversion\bin \earth2\bank\tomo\bin note that in a l l of John's code cartesian coordinates are t (theta) latitude p (phi) longitude r (radius) depth 1) Multi-channel crosscorrelation (mccc) note that you do not have to f i l t e r or resample before. Seed to cross correlate more than 3 stations (else a l l residual times can be explained by event mislocations) I use: mcee -w 3 - i 0.5 -s 0.4 - f BP 0.4 1.6 2 -ipick A -wpick T7 *.SAC options: v correlation window length i time substracted from i n i t i a l pick to start of window s cross-correlation s h i f t (here 0.4 s to avoid cycle skip) f f i l t e r (here bandpass from 0.4 to 1.6 Hz with 2 passes) ipick read i n i t i a l pick from time stamp in sac header wpick write mccc times to specified time stamp in sac header check man page mccc output f i l e Column Description 1 Station name. 2 Relative delay time (with zero mean). 3 Sample standard deviation. 4 Average cross-correlation coefficient. 5 Sample standard deviation of cross-correlation coefficient in Fischer-transform space. To obtain error bounds on cross-correlation coefficient, one must transform the average value to Fischer space, add and subtract a standard deviation and then transform those values back to to the cross-correlation domain, (see VanDecar and Crosson, 1990, BSSA, vol. 80, pp. 150-169) 6 Flag to indicate i f polarity of station i s reversed. This is not a problem with the portable broadband data, but often with the short-period network data some stations are known to have reversed polarity. 0 - normal polarity, 1 = reversed. 7 F i l e name. 2) set up model parameterization and zero reference model for inversion see folder /earth2/bank/tomo/src/conmod similar f i l e s to be created for other regions. Knot locations are specified in conhud.f note: i f changes are made you - need to re-edit f i l e /earth2/bank/tomo/include/mod3d.h - keep old versions of changed f i l e s , otherwise binary model f i l e s cannot be read can use chkmod to check model (see /earth2/bank/tomo/src/utils ) remember to copy ref.model into ../partials ( w i l l become mod.zero) 3) create reference f i l e s for inversion see /snorcle/inversion/src/procx/README a l l mccc output f i l e s have to be in l o c a l directory (named e.g. 96031.20.30.47px ) plus f i l e loc.sta (coordinates of a l l stations) awk '-Cprintf " %-4s X8.5f XlO.Sf KS.3f \n" ,$1 ,$2,$3,*4/1000>» \ /earth2/bank/1 ables/st at ion*dat run (on ivana): Is *px I procx then: sh cut.sh D PRACTICAL NOTES ON VANDECAR'S CODE 91 reference f i l e s are: ref.trays info on every ray (source-receiver pair) column description 1 ray number 2 event 3 station 4 r e l . delay time } 5 sample standard deviation } from mccc 6 cross-corr. coefficient } 7 take-off angle (radians) ] Id earth now (•) 8 azimuth (radians) ] 9 travel time (sec) 10 0.00000 (•) 11 phase (•) s i l l be changed for non-linear case 7 and 8 are then for the 3D earth model 10 then denotes by how many degrees eq location was missed (value should f a l l below 0.1 degree) ref.teqs event, PDE entry (time, location, magnitude) ref . t s t a station name, l a t , Ion, elevation ref.tcale event, event mistiming (0.000) * ref.teals station name, station correction(O.OOO) *, no. of rays (*: compare to new.teal* after the inversion, see E) 4) calculate matrix of p a r t i a l derivatives read /earth2/bank/tomo/partials/READHE ( needs f i l e s r ef.tsta, ref.trays, ref.model, iasp91.vel and partial.sh ) copy iasp91.vel and partial.sh into l o c a l directory get r e f . t s t a and ref.trays (copy from procx or link i t ) (in non-linear inversion use the new.rays) create soft l i n k to ref.model: In -s mod.zero ref.model (note: l i n k to ../../sre/conmod/ref.model when you are starting, or to output of last linear inversion in non-linear case) now run s p l i t -S00 ref.trays ref.rays. ; compress ref.rays.?? run partial.sh a b (c... i f you have more ref.rays.??) p a r t i a l s i l l run a long time in ivana creating part.out.??.Z t i t ' s a good idea to nice partial.sh :) = smiling supervisor] nice partial.sh a b c d > partial.log (note that partial.sh s i l l run the sequence setpar ; sort r h i t s > rhits.sort ; p a r t i a l for each of the ref.rays, run grabig on each of the part.out.??.Z f i l e s , l i k e zcat part.out.ad.Z I grabig 1.0 > tmp sort tmp > part.sort_ad ; compress part.sort.ad ; rm tmp and combine sorted p a r t i a l derivative f i l e s zcat part.sort.??.Z I compress > part.modXX.Z where XX i s number of your inversion 5) linear inversion read /earth2/bank/tomo/models/READHE.dtelin (inversion needs the following 8 f i l e s ref.trays ref.tsta ref.teqs ref.tcale ref.teals ref.model part.modXX.Z telin.su ) cp ../partials/part.modXI.Z . cp ../proex/ref.t* . In -s ../partials/mod.zero ref.model or l i n k from sre/conmod as ref.model and copy telin.su from earlier model, change entries as desired check i f r e f . t * are sorted run (on ivana) nice zcat part.modXX.Z I tho.dtelin > dtelin.log ft nice zcat part_mod03.Z I /earth2/bank/tomo/src/dtelin3/tho_dtelin > dtelin.log ft nice nightshell.sh > nightshell.log ft i f problems occur check number of stations, events, rays in /earth2/bank/tomo/src/dtelin/dtelcm.h and matrix dimensions etc. in /earth2/bank/tomo/src/dtelin/dscgls.h then recompile dtelin (run "make dt e l i n " in /earth2/bank/tomo/src/dtelin remove ».o f i l e s before recompiling) output f i l e s (after 12 robust iterations and 2000 CG iterations) - new.bmod.12.2000 binary model (knot locations and slowness perturbations in s/km : to convert these to a s c i i see comments on Xmap8 below) - new.tcale.12.2000 earthquake relocations D PRACTICAL NOTES ON VANDECAR'S CODE 1 •vent 2 timing correction (s) 3 latitude correction (in km) 4 longitude " E depth " - new.teals.12_2000 station corrections 1 station 2 timing correction (s) 3 number of rays for that station - resid_12 residual information (as written to screen also) 1 robust iteration 2 conjugate gradient iteration 3 rms error of A transpose residual 4 rms error of travel-time equations 5 rms error of amplitude 6 rms of damping 7 rms of flattening (1st derivative) 8 rms of smoothing (Laplacian) 9 rms of jerking (3rd derivative) 10 t o t a l rms of a l l equations (values decreasing monotonically) - resids_all.012 residuals for each ray ordered by ray number 6) analysis check for possible cycle skipping by comparing residuals for each ray paste resids.all.001 ref.trays I sort -n > checkOl rather high, or rather low, values usually indicate cycle skipping copy residuals above 0.2 and below -0.2 into separate f i l e "too.high" and sort by event sort +2 too.high > too_high_events usually one event with a large positive residual w i l l have a large negative residual also => cycle skip redo crosscorrelation i f you notice cycle skip REMEMBER: ray paths do not change unless you matched seismograms to a wrong earthquakes. Therefore do not have to recalculate p a r t i a l s ! However, need to update reference and partials f i l e : (i) new crosscorrelation: - change entry in ref.trays (make sure numbers match, crosscorrelation may have changed order) ( i i ) single rays to leave out: - remove in ref.trays - remove in part.modXX (ray number is the f i r s t number) ( i i i ) whole event to leave out: - remove from ref.teqs - remove from ref.tcale - remove entries in ref.trays - remove a l l corresponding rays i n part.modXX Bote: program w i l l not complain about missing ray numbers but rays should be sorted in r i s i n g order a) create postscript sections through binary model read /earth2/bank/tomo/models/README.plot to mask " c e l l s " with few rays crossing run zcat part.modXX.Z I tho.hitcnt don't forget to copy loc.sta into l o c a l directory to have current stations on depth plots (else depth.ps w i l l look for an older default f i l e ) be on ivana depth.csh 200 new.bmod.12.2000 > depth_200.12_2000.ps slice.A.csh new.bmod.12.2000 > slice_A.12_2000.ps b) plot relocation vectors note that ref.teqs has PDE lines in same order as relocations appear in new.tcale* (ref.teqs has to be in directory) cat new.tcale.01.2000 I awk '{print $3, $4, *6}' |\ plotre -o S3 -10E -c > tmp.ps 'origin "colour c) to u t i l i s e Xmap8 see /earth2/bank/xmap8/READHE (to convert binary model to a s c i i and use this as input for matlab code which projects lat/long knots to km grid) 7) to obtain tradeoff curve - run inversion with smoothing turned off how well does data f i t ? D PRACTICAL NOTES ON VANDECAR'S CODE - run inversion with flattening turned off which value do you have to choose to get similar/same lev e l of regularization? (in practise Cflat * Csmooth / grid spacing ) - now change flattening and smoothing weights just found by factors of 2 to work through the tradeoff curve... - rms travel-time residual reduction i s 4th column in resid_XX ?? - rms model roughness and smoothness are 7th and 8th column 8) resolution test a) create synthetic model see /earth2/bank/tomo/src/synth/README synth.f s i l l create one spike at a time (current version of synth.f; John has also implemented slabs, cratons and f o l d belts) I have written a shell script that w i l l add any number of spikes as defined in spikes.dat to come up with a f i n a l model P.synth.model (see /earth2/bank/tomo/inv_02/Resolution/createsynth.sh) copy loc.sta into l o c a l directory check output with depth.csh and slice.A.csh as described in (6). b) forward model to get synthetic data cp .,/partials/part_modXX.Z . cp .,/procx/ref.trays orig.ref.trays synthetic travel times are produced by running (very fast) zcat part_modXX.Z I tho_forwrd -d orig.ref.trays -m P_synth.model\ -p - > synth.ref.trays now add random noise to the synthetic times (try using 0.1 s for a start - i t ' s a bit large compared to the standard error of the data, but the real noise probably has a l l sorts of biases and correlations in i t also). addnois -e 0.1 < synth.ref.trays > synth.noise.ref.trays c) run inversion In -s synth.noise.ref.trays ref.trays copy mod.zero, ref.teal?, ref.teqs, ref.tsta, telin.su into l o c a l directory In -s mod.zero ref.model and run inversion as shown in (5) 9) non-linear inversion a) ray trace through 3D model obtained before see /earth2/bank/tomo/raytrc create directories mkdir raytrc ray.tables ; cd raytrc and copy ref.trays ref.tsta ref.teqs into raytrc need ID model: iasp_91.vel and 3D model: In -s ../model.l/new.bmod.12.2000 mod.P96al (important: mod.P96al i s name implemented in raytracing code) sort +2 ref.trays I splref w i l l create ray f i l e for each station by s p l i t t i n g ref.trays (see /earth2/bank/tomo/raytrc/Iew_rays/ and /Old.rays ) and calculate table from delta 30 to 100 degrees over a l l azimuths (10 degree spacing), see /earth2/bank/tomo/ray_tables get new ray paths (picks 3 closest rays from tables calculated above and uses lewton's method to compute new ray paths) mkdir Old.rays Bew.rays copy ray.sh into raytrc and run ray.sh rays.* this s h e l l c a l l s raytrc and creates a new reference f i l e Bew.rays/rays.total . Entries in Hew_rays/rays.total are similar to those in ref.trays (3); changes: column 7 new take-off angle (radians) 8 new azimuth (radians) 10 by how many degrees was eq location missed (value should f a l l below 0.1 degree for a l l rays) ( i f not a l l rays have been processed: may need to change no of events in /earth2/bank/tomo/src/raytrc/raytrc.f and recompile) b) recalculate partials (see 4) remember that ref.trays i s output Bew.rays/rays.total (see a) and ref.model i s new 3D slowness model c) run inversion as in the linear case (see 5), using new ref.trays and 3D slowness model 94 E Comparison of heat diffusion for various geometries This appendix shows results of heat conduction in different symmetry systems. The cal-culations are done using similarity solution for the following cases: - cartesian symmetry (i.e. homogeneous half space heating from a plane boundary), - cylindrical symmetry (infinite cylinder of radius a heating a homogeneous space), and - spherical symmetry (sphere of radius a heating a homogeneous space). Boundary conditions applied are the same in all three cases, namely that the temperature within the anomaly and at infinite distance from the anomaly remain constant (i.e. no cooling of the anomaly in these models). Mathematical solutions for the three cases are outlined on page 96. The results for the cylindrical and spherical symmetries differ from the one for cartesian symmetry by factors of y/r/a and r/a, respectively, and all results involve an error function. cartesian cylindrical 1 O O 2 0 0 t i m e ( M a ) T3 2 0 1 O O 2 0 0 t i m e (!N<Ia) spherical 1 O O 2 0 0 t i m e (!N<Ia) Figure E . l : Comparison of heat conduction for three geometrical cases. Shown are cartesian, cylindrical and spherical geometries. The parameter is given in fraction of A T and the radius of the infinite cylinder or the sphere is fixed at 25 km. \//c£ is drawn into all three plots to ease comparison. E COMPARISON OF HEAT DIFFUSION FOR VARIOUS GEOMETRIES 95 Graphical representations of these results are shown in figures E . l and E.2. K =1.110~6 m 2 s - 1 was chosen as the thermal conductivity of lithosphere (Stacey, 1992). Figure E . l shows that at an arbitrary distance away from the heat source the temperature rises fastest in the spherical case followed closely by the cylindrical case. A length scale of y/ni corresponds to .4795AT, where A T = TS-T0 ' m * n e c a rtesian case. The effect of the size of the anomaly on the diffusion distance is shown in figure E.2. The lines shown correspond to .4795AT with A T as defined above. After 100 Ma .4795AT is at nearly constant distance from sphere and does not extend further with time for the cylinder (radii < 50 km). After 100 Ma a sphere with 50 km diameter has heated its surroundings to 30 km, a cylinder has heated to 40 km and the half space has heated to 60 km away from source. cylinder sphere 100 200 time (Ma) 100 time (Ma) 200 Figure E.2: Effect of size of anomaly on heat diffusion. Shown are lines of .4795AT for various radii of a cylinder or sphere, s/rti is added in both panels. E COMPARISON OF HEAT DIFFUSION FOR VARIO US GEOMETRIES 96 to a > c o '•§ 3 cr 4) a o ed 4) D X +» bO o IT •c v •a COP d § V M tig V T3 O • ^ ^ _ ^ «R ^ « 8 ft ft ft E ? E ? EH II II o EH EH EH 8 T I *T o 8 <^c> ft ft ft E ? EH E ? II II o f EH EH M ii V l ° S o H O O cc> <^ > ft ft ft EH5 EH° £ II £ S, i EH EH J l a o •i-t d o o o •a 9 .|co <0F 2 8 V u d 8 .a II Cr-•c o o o ^> s o d _o v» J9 d o u o. 3 d d o d i O II II <*> •8 II t i « V eli. E ? « v 4 •8 d v •c J* <o u a .2 d -3 97 F Zusammenfassung - Abstract in German Zwischen Oktober 1994 und Juli 1996 zeichnete ein Array von Drei-Komponenten Seis-mographen uber dem zentralen Abschnitt des Trans-Hudson Orogens Erdbeben auf. Ziel dieser Studie war es, Aufschlufi uber die Struktur des oberen Erdmantels unter einem Vorkommen diamantfuhrender Kimberlite zu gewinnen. Das temporare Stationsnetz be-i stand aus neun Breitband und acht kurzperiodischen Rekordern mit einem Stationsabstand von etwa 100 km. Es iiberdeckte das Glennie Domain und seine Nachbargebiete, deren Grenzen unter den phanerozoischen Sedimenten der Plattform mit Hilfe aeromagnetischer Daten festgelegt worden waren. Erganzt wurde der gewonnene Datensatz durch Aufzeich-nungen sowohl einer permanenten Station des IRIS Consortiums als auch einer im Jahre 1991 durchgefuhrten Pilotstudie. Aus Laufzeitresiduen von 331 Erdbeben wurden Slownessvariationen bestimmt. Dazu wurde ein Tomographie-Programm von VanDecar (1991) eingesetzt. Die Strahliiberdeckung erlaubt Aussagen liber einen Tiefenbereich zwischen 60 und 400 km. Unsere Analysen zeigen Abweichungen von ±2% vom radialsymmetrischen iasp91 Erd-modell im Erdmantel oberhalb 220 bis 270 km. Das Modell zeigt die grofite Slowness (niedrigste Geschwindigkeit) in einer vertikalen zyhnderformigen Anomalie von 120 km Durchmesser zwischen 60 und 220 km Tiefe mit Maximum bei (52.8°N, 105.9°W, 100km Tiefe). Die Anomalie niedrigster Slowness (hochste Geschwindigkeit) lafit sich als nie-renformig beschreiben. Das Anomaliemaximum liegt bei (51.8°N, 103.9°W, 170km Tiefe). Uber oder nahe den Randern von Niedrig-Geschwindigkeits-Anomalien befinden sich kretazische diamantfiihrende Kimberlite sowie Orte mit erhohten Konzentrationen kim-berhtischer Minerale in quartarem Till. Wir deuten diese Anomalien als Uberreste eines Mantel-Plumes oder einer Rayleigh-Taylor-artigen Instabilitat. Die Korrelation F ZUSAMMENFASSUNG - ABSTRACT IN GERMAN 98 der Hoch-Geschwindigkeits-Anomalie mit einem langwelligen Schweretief spricht fur eine Interpretation dieser Struktur als Lithospharenkiel. Diese Deutung bestarkt Vermu-tungen, wonach der vor kurzem aus reflexions- und refraktionsseismischen Ergebnissen pos-tulierte Saskatchewan Kraton eine wesentliche Struktureinheit des Trans-Hudson Orogens darstellt. Niedrig-Geschwindigkeits-Anomalien im Siiden des Arbeitsgebietes korrelieren mit erhohten Warmeflufiwerten. Daraus leiten wir einen thermischen Ursprung dieser Ano-malien ab. Eine konvektierende Asthenosphare wird im Tiefenbereich von 220 zu 270 km bis 400 km als vergleichsweise homogene Zone abgebildet. Wir interpretieren die Slownessanomalien als Reminiszenz einer thermomechanischen Erosion der Lithosphare des Saskatchewan Kratons wahrend der Kreidezeit. Die kretazi-schen diamantfuhrenden Kimberlite sind eine direkte Folge dieses Prozesses. Uberbleibsel des Lithospharenkiels werden als Hoch-Geschwindigkeits-Anomalien abgebildet und Ein-lagerungen von Gesteinen aus tieferen Regionen des Erdmantels in die Lithosphare als Niedrig-Geschwindigkeits-Anomalien. Unsere Studie deutet auf eine Machtigkeit der Lithosphare von 220 bis 270 km unter dem Trans-Hudson Orogen hin. 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items