Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Linear systems theory applied to a horizontally layered crust Jensen, Oliver George 1970

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

Item Metadata

Download

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

Full Text

LINEAR SYSTEMS THEORY APPLIED TO A HORIZONTALLY LAYERED CRUST by OLIVER GEORGE JENSEN B.Sc. U n i v e r s i t y of B r i t i s h Columbia, 1964 M.Sc. U n i v e r s i t y of B r i t i s h Columbia, 1966 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in the Department of GEOPHYSICS We accept t h i s t h e s i s as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December, 1970 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of Brit ish Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Depa rtment The University of Brit ish Columbia Vancouver 8, Canada ABSTRACT E l a s t i c wave propagation i n a m u l t i l a y e r e d crust with causally or acausally attenuating layers i s formulated d i r e c t l y i n terms of l i n e a r systems theory. The s o l u t i o n of the l i n e a r systems analog determines the wave motions at the free surface, motions of a l l i n t e r n a l boundaries and the waveforms propagated i n t o the mantle i n response to plane P or S waves generated w i t h i n the l a y e r i n g or entering i n t o the cr u s t . The p a r t i c u l a r problem solved i s the determination of the response of an n-layered crust to t e l e s e i s m i c P or S waves i n c i d e n t with a r b i t r a r y angle at the c r u s t a l base. T r i v i a l extensions of t h i s problem would allow r e f l e c t i o n s o l u t i o n s . Numerical s o l u t i o n s have been accomplished i n both the frequency domain and the time domain. The F o u r i e r s o l u t i o n r e s t r i c t e d to a non-attenuating crust i s equivalent to the standard H a s k e l l matrix s o l u t i o n of e l a s t i c wave equations. D i r e c t time domain s o l u t i o n s allow the syntheses of seismograms considering i n t e r n a l c r u s t a l absorption. For demonstration of the u t i l i t y and advantages of the theory, the l i n e a r systems formulation has been applied to studies of the e l a s t i c wave response of the c e n t r a l A l b e r t a crust using P waves from 6 t e l e -seismic events. Frequency domain comparisons between the t h e o r e t i c a l and experimental s p e c t r a l amplitude V/H r a t i o s have shown that, although the t h e o r e t i c a l e f f e c t of attenuation w i t h i n the crust can be considerable, l i t t l e improvement i n c o r r e l a t i o n s between theory and experiment has been achieved by considering p l a u s i b l e c r u s t a l absorption models. Although s i g n i f i c a n t s i m i l a r i t i e s between the t h e o r e t i c a l and experimental V/H r a t i o s were found below 2 Hz, l i t t l e i i c o r r e l a t i o n was apparent at higher frequencies. Background and s c a t t e r i n g noise p a r t l y contributed to t h i s e f f e c t and i t i s also probable that i n s u f f i c i e n t d e t a i l and accuracy was a v a i l a b l e f o r the model c r u s t a l sections. Time domain synthetic seismograms have been determined which w e l l correspond to the e a r l y P coda of several of the experimental records. Assumptions on the event source motions and the mantle pro p e r t i e s are required to determine i n c i d e n t P waveforms f o r these s o l u t i o n s . Causal attenuation w i t h i n the c r u s t a l l a y e r i n g was included. Correlations between the s y n t h e t i c seismograms and the experimental records has been found to decrease r a p i d l y with time following the P onset. I t i s suggested that t h i s e f f e c t i s p r i m a r i l y due to background noise and p o s s i b l e s c a t t e r i n g of the waves w i t h i n the c r u s t . Furthermore, i t i s probable that the waveforms used f o r these s o l u t i o n s based on the source motion and mantle attenuation assumptions were not s u f f i c i e n t . A major apparent advantage of the new formulation i s that causal and acausal attenuation s o l u t i o n s are permitted i n both the frequency and time domains. Al s o , the large body of communications theory mathematics can now be applied d i r e c t l y to the propagation problem and could prove u s e f u l i n attempts at the s o l u t i o n of the non-normal incidence inverse problem. i i i TABLE OF CONTENTS CHAPTER I GENERAL INTRODUCTION 1.1 Introduction 1 1.2 Outline of Thesis 4 CHAPTER II THE LINEAR SYSTEM THEORY FOR A MULTILAYER CRUST 2.1 D e s c r i p t i o n of the Problem 6 2.2 Linear Systems Approach 6 2.3 The Linear System Analog for Normal Incidence 9 2.4 The Frequency Domain Desc r i p t i o n 11 2.5 The Linear System Analog for Non-normal Incidence 18 2.6 Time Domain Formulation 27 2.7 Attenuation and Dispersion i n the System Model 34 CHAPTER I I I ANALYSIS OF THE ALBERTA SEISMOGRAMS 3.1 Introduction 38 3.2 The Analysis Problem 38 3.3 The Experiment " 3 9 3.4 The Crust Model 41 3.5 Data Preparation 48 3.6 The Seismograms f 50 3.7 Frequency Analyses - T h e o r e t i c a l S p e c t r a l Amplitudes 59 3.8 Frequency Analysis - Spectral Ratio Technique ^ 63 3.9 T h e o r e t i c a l Spectral Ratios , 6 8 3.10 Experimental S p e c t r a l Ratios 71 i v 3.11 General Features of Spectral Ratios 79 3.12 D e t a i l s of the Sp e c t r a l Ratios 81 3.13 Time Domain Syntheses 85 3.14 Syntheses with Layer Attenuation 88 3.15 Impulse Synthetic C r u s t a l Responses 92 3.16 Incident Waveform Models 97 3.17 Seismogram Syntheses 106 3.18 Comparisons of the Syntheses and Records f o r 107 55° Epicenter 3.19 Comparisons of the Syntheses and Records f o r 110 85° Epicenter 3.20 Synthesis Using an Input Wavelet Obtained by 110 Deconvolution CHAPTER IV CONCLUSIONS 4.1 Summary of Analyses 113 4.2 Concluding Remarks and Suggestions 115 REFERENCES 118 APPENDIX I IMPULSE RESPONSE FUNCTION FOR NORMAL INCIDENCE 121 APPENDIX II PHASE ADVANCE TRANSFER FUNCTIONS 128 APPENDIX I I I EXAMPLE BLOCK REDUCTION PROCEDURE 130 APPENDIX IV FORTRAN COMPUTER PROGRAMS FOR TIME DOMAIN 136 SEISMOGRAM SYNTHESES AND FOURIER COMPONENT SOLUTIONS BASED ON THE LINEAR SYSTEMS THEORY V LIST OF TABLES TABLE I C r u s t a l Sections 44 TABLE I I Formation Q's 47 TABLE I I I Events Used i n Analyses 57 TABLE IV Attenuator Weights f or the 0.05 Second 90 Sampling Increment - • LIST OF FIGURES FIGURE 1 The h o r i z o n t a l l y multilayered medium with n l a y e r s o v e r l y i n g a basement halfspace. FIGURE 2 The input-output diagram for a l i n e a r system cr u s t . FIGURE 3 A representation of the possible input and output waveforms i n l a y e r i under normal incidence. FIGURE 4 A ray path representation of the i n t e r n a l reverberations w i t h i n l a y e r i g i v i n g r i s e to output transform s i g n a l s U 1 ( s ) and D i + l < S > -FIGURE 5 A f u n c t i o n a l block representation of the l i n e a r system analog for layer i under normal incidence. FIGURE 6 A representation of the "pseudo-coordinate" transformation f o r l a y e r i . FIGURE 7 A diagram showing the non-coincidence of wave constant phase p o s i t i o n s and a h o r i z o n t a l l a y e r boundary. FIGURE 8 A f u n c t i o n a l block representation f o r non-normal incidence f o r layer i . FIGURE 9 The free surface r e f l e c t i o n f o r non-normal incidence. FIGURE 10 A comparison of a normal incidence l a y e r impulse response f o r i n f i n i t e Q and ca u s a l l y absorbing conditions. FIGURE.11 The time-domain representation of the system analog of one l a y e r at normal incidence. FIGURE 12 The time-domain representation of the system analog of one layer f o r non-normal incidence. FIGURE 13 A map of the seismic s t a t i o n and o i l w e l l l o c a t i o n s by which the c r u s t a l s e c t i o n was determined. v i i FIGURE 14 The sedimentary c r u s t a l section obtained by i n t e r p o l a t i o n between o i l w e l l geophysical logs. FIGURE 15 The t h e o r e t i c a l standard seismograph response compared to a response c a l i b r a t i o n obtained by the tra n s i e n t method. FIGURE 16 The v e r t i c a l , r a d i a l h o r i z o n t a l and transverse h o r i z o n t a l components of the experimental seismograms f or the events studied i n t h i s work. FIGURE 17 The s p e c t r a l amplitude response of the t h e o r e t i c a l non-attenuating Leduc c r u s t a l s e c t i o n to a normal incide'nce white amplitude spectrum P-wave incident at the Precambrian-Cambrian boundary. FIGURE 18 The s p e c t r a l amplitude response of the t h e o r e t i c a l non-attenuating Leduc c r u s t a l s e c t i o n to a 30° i n c i d e n t P wave at the Precambrian-Cambrian boundary. FIGURE 19 The s p e c t r a l amplitude response of the t h e o r e t i c a l c r u s t a l s e c t i o n to a 30° incidence P wave at the Precambrian-Cambrian boundary. FIGURE 20 S p e c t r a l V/H r a t i o s f o r various attenuation conditions compared to the non-attenuating crust r a t i o s . 42 49 51 60 61 62 66 FIGURE 21 V e r t i c a l - t o - v e r t i c a l s p e c t r a l r a t i o s f o r non-attenuating crust models. FIGURE 22 Experimental V/H s p e c t r a l r a t i o s obtained from the s u i t e of seismograms. FIGURE 23 Sp e c t r a l s i g n a l - t o - n o i s e r a t i o s f o r the v e r t i c a l and h o r i z o n t a l f o r three events. 70 72 78 FIGURE 24 T h e o r e t i c a l v e r t i c a l - t o - v e r t i c a l r a t i o s are compared to those obtained from the experimental seismograms f o r 30° incidence angle. FIGURE 25 Comparison of non-normal and normal incidence s y n t h e t i c seismograms. 84 87 FIGURE 26 The impulse response of l a y e r 1 of the c r u s t a l s e c t i o n assuming causal attenuation with Q=23. 89 v i i i FIGURE 27 The t h e o r e t i c a l frequency response of the crust f o r one passage of an impulse for average Q=200 compared to the frequency response due to the sum of the time domain attenuation operators of Table IV. FIGURE 28 Normal incidence s y n t h e t i c seismograms dependence on Q. FIGURE 29 The v e r t i c a l and r a d i a l h o r i z o n t a l s y n t h e t i c impulse responses of the three c r u s t a l sections representing the three s t a t i o n s . FIGURE 30 The normal incidence s y n t h e t i c impulse response f o r the LED c r u s t a l s e c t i o n . - • FIGURE 31 The basement r e f l e c t e d wave sy n t h e t i c impulse response f o r the LED c r u s t a l s e c t i o n , FIGURE 32 The e f f e c t of various mantle path average Q values on an impulse source displacement motion which has t r a v e l l e d 10 minutes at the undispersed phase v e l o c i t y through the mantle. FIGURE 33 Wavelets constructed by convolution of the seismograph response with absorbing mantle responses f o r various T/Q parameters f o r both impulse and step sources. FIGURE 34 V e r t i c a l and h o r i z o n t a l component s y n t h e t i c seismograms f o r various T/Q wavelets of the form shown i n Figure 32 f o r impulse source motions. FIGURE 35 V e r t i c a l and h o r i z o n t a l component s y n t h e t i c seismograms f o r various T/Q wavelets f o r step source motions. FIGURE 36 Impulse source motion s y n t h e t i c seismograms f o r 55° and 85° e p i c e n t r a l distances. FIGURE 37 A comparison of the v e r t i c a l s y n t h e t i c seismograms with the experimental seismo-grams of Figure 16 for s t a t i o n LED. FIGURE 38 A comparison of the v e r t i c a l component LED seismogram f o r the Venezuela event and the re-synthesis using the Ulrych wavelet. FIGURE 39 Pole maps of the one l a y e r t r a n s f e r f u n c t i o n f o r normal incidence: Case I. 91 93 94 95 96 100 101 102 103 105 108 111 125 ix FIGURE 40 Pole maps of the one l a y e r t r a n s f e r function f o r normal incidence: Case I I . 126 FIGURE 41 A representation of phase advance time c a l c u l a t i o n s necessary f o r non-normal condition analyses. 128 FIGURE 42 The desired block diagram representation i n terms of t r a n s f e r function and transformed s i g n a l s . 130 FIGURE 43 Block diagram of p a r a l l e l summing paths. FIGURE 44 Block diagram of the l i n e a r p o s i t i v e feedback system. FIGURE 45 The reduced system block diagram having the la y e r response to inputs U^(s) and V_^(s). 130 131 134 X ACKNOWLEDGEMENTS I wish to thank Dr. R. M. E l l i s for his invaluable aid and encouragement during the development of this research. I am also grateful to my advisory committee with Drs. G. K. C. Clarke, R. M. Clowes and R. D. Russell for useful discussions and criticisms of various phases of the work. I would also like to acknowledge the contribution of the Department of Geophysics staff whose preparation of the seismograph f i e l d equipment, recording of the experimental data and preparation of the digitization equipment was so necessary to much of this study. During this experiment, equipment was made available by the Seismology Division, Earth Physics Branch and both equipment and f a c i l i t i e s were provided by the University of Alberta. I would also like to thank Judy McLean for a fine job typing the manuscript. This research was supported by the National Research Council (Grant A-2617) and the Defence Research Board of Canada (Grant 9511-76). The author was partly supported during this research by a University of B r i t i s h Columbia Graduate Fellowship. CHAPTER I GENERAL INTRODUCTION .• .1.1 Introduction Seismologists have been i n v e s t i g a t i n g the propagation of e l a s t i c waveforms i n h o r i z o n t a l l y multilayered media since the early development of seismic exploration methods i n order to determine the c r u s t a l structure-and the e f f e c t s of the structure on seismic waves. In addition to the obvious p r a c t i c a l applications i n o i l and mineral exploration, a knowledge of c r u s t a l propagation i s necessary to the continuing advance of earthquake seismology. For example, much of the character of a seismogram i s gener-ated by the e l a s t i c properties of the crust underlying a seismograph s t a t i o n . Thus i t i s useful to remove these complicating e f f e c t s i n order to study the d e t a i l s of f o c a l mechanisms and to r e f i n e our knowledge of the earth's i n t e r i o r . Furthermore, since c r u s t a l information i s contained i n earthquake seismograms, analysis procedures based upon a thorough understanding of e l a s t i c propagation could be u s e f u l to obtain c r u s t a l structure from these seismograms. The f i r s t s a t i s f a c t o r y mathematical analysis of e l a s t i c wave propagation i n h o r i z o n t a l l y layered structures was that of Thomson (1950) which was subsequently generalized to include both body waves and surface waves by Haskell (1953). As w e l l as surface wave a p p l i c a t i o n s , several i n v e s t i g a t o r s (e.g., Phinney, 1964; Ibrahim, 1969) have applied the Haskell matrix method to the analysis of earthquake body waves. The Haskell-Thomson approach obtains the propagation s o l u t i o n i n terms of frequency component amplitudes; the more recent work of Wuenschel (1960) f o r normal incidence and F r a s i e r (19 70) f o r non-normal P and SV motions permit the determination of time domain c r u s t a l responses. For v e r t i c a l incidence P waves, extensions of these methods and other a n a l y t i c a l approaches to the inverse problem of determination of the la y e r paramenters from earthquake and exploration seismograms have received considerable a t t e n t i o n i n the recent l i t e r a t u r e (e.g., Goupillaud, 1961; Claerbout, 1968; Ware and A k i , 1969). This t h e s i s presents a new formulation f o r the d i r e c t propagation problem using l i n e a r systems theory which allows both P and SV motions of a r b i t r a r y angle with the v e r t i c a l w i t h i n the l a y e r i n g and s o l u t i o n i n e i t h e r the frequency or time domain. Linear systems and analog procedures have been used p r e v i o u s l y . For example, the more r e s t r i c t i v e one-dimensional problem of v e r t i c a l P-wave propagation has been studied i n systems (Lindsey, 1960), analog (Sherwood, 1962) and communications theory (Rob inson and T r e i t e l , 1966) terms. Also, surface wave normal d i s p e r s i o n has been described by Mueller and Ewing (1962) by means of a l i n e a r systems model. I t i s notable that when an e l a s t i c wave t r a v e l s through a str u c t u r e d m a t e r i a l , a complex record of the e l a s t i c p r o p e r t i e s along the t r a v e l path i s retained by the wave. Hence, analysis of the wave can i provide s t r u c t u r a l information. This i s the e s s e n t i a l goal of the inverse seismic problem. The r e l a t i v e l y simple geometry of the h o r i z o n t a l l y m u l t i l a y e r e d crust leads one to expect that inverse s o l u t i o n s may be f i r s t p o s s i b l e f o r t h i s problem. However, the complete inverse s o l u t i o n also demands that the propagation c h a r a c t e r i s t i c s due to a l l e l a s t i c e f f e c t s i n the crust are w e l l understood. This t h e s i s should help to extend t h i s understanding. 3 The i n c l u s i o n of attenuation, both causal and acausal i s new i n th i s work. However, extensions of other approaches may also allow the feature. For instance, Haskell (1953) suggested a possible means for including attenuation i n his formulation; derivatives of the Haskell method (e.g., Wuenschel, 1960; Fr a s i e r , 1970) could possibly be modified to permit c r u s t a l absorption. The l i n e a r systems approach has a further advantage that the dir e c t time and frequency domain problems are cast i n very s i m i l a r forms. Thus, the correspondence of time series and spectral response character-i s t i c s are p a r t i c u l a r l y apparent. The equivalence of the time and frequency solutions by other approaches i s not d i r e c t . One outstanding feature of a systems formulation i s that a great body of a n a l y t i c a l techniques already e x i s t s , p a r t i c u l a r l y i n automatic control engineering and communications engineering, which now can be applied to the problem. Possible further success i n determining c r u s t a l properties from earthquake and r e f l e c t i o n seismograms could be derived by using these mathematical tools. In the past, communications theory has been of great interest to.many geophysicists and much recent advance of geophysics has been gained using these methods. Thus, a dir e c t formu-l a t i o n of the propagation problem i n the systems form may even broaden the i n t e r e s t among these geophysicists already using communications theory methods. A l a s t consideration which may be of consequence i s the fact that the systems form of the problem can lead to d i r e c t mechanical and e l e c t r i c a l analog procedures. With the recent developments i n micro-e l e c t r o n i c d i g i t a l l o g i c , for example, small hard-wired computers could e f f i c i e n t l y solve propagation problems. The p o s s i b i l i t y of analog time domain approaches also e x i s t s . 4 The basic assumptions required by the linear systems approach are similar to those used in other methods. Distinct layers having isotropic velocities, uniform densities and welded boundaries are necessary. Plane rotational and compressional waves are considered. Attenuation with an arbitrary complex dependence on frequency i s also allowable in each layer. The general solution of the problem requires the determination o f the motions on the free surface and within the medium in response to source wave motions generated within the layering or entering into the medium. The choice of the source location of the plane wave generation determines the nature of the physical problem; either the reflection or transmission problem is soluble. 1.2 Outline of Thesis (a) The systems description of crustal transmission of elastic waves is f i r s t obtained in the frequency domain for normal incidence plane waves impinging on a horizontally multilayered crust from a basement halfspace. This restricted formulation is then generalized to non-normal conditions. (b) The analogous time domain formulation i s derived from the general frequency description. (c) The u t i l i t y of the theory i s demonstrated in comparison analyses with a set of earthquake and nuclear event seismograms recorded in central Alberta. Both time and frequency analyses as well as seismo-gram synthesis are obtained. (d) Extensions of the method are suggested which could allow analysis of propagation in crusts possessing complex linear and non-linear properties. The possibility of inverse solutions is also indicated. 5 The p r o b l e m has been s p e c i f i c a l l y f o r m u l a t e d to d e s c r i b e p r o p e r t i e s o f t e l e s e i s m i c seismograms: t h a t i s t h e c r u s t a l t r a n s m i s s i o n p roblem. The o b v i o u s e x t e n s i o n t o p e r m i t n o r m a l i n c i d e n c e r e f l e c t i o n s o l u t i o n s i s n o t e d . 6 CHAPTER I I THE LINEAR SYSTEM THEORY FOR A MULTILAYER CRUST 2.1 D e s c r i p t i o n of the Problem In t h i s work e l a s t i c propagation i s to be formulated i n two dimensions f o r waveforms i n c i d e n t on a h o r i z o n t a l l y l a y e r e d medium from a h a l f s p a c e beneath. This model represents a t e l e s e i s m i c a r r i v a l impinging on a l a y e r e d c r u s t from a basement. Each c r u s t a l l a y e r i i s s p e c i f i e d by i t s - compressional and shear wave v e l o c i t i e s a., 6.: d e n s i t y p. and t h i c k n e s s n.. The P and SV ray angles i n each l a y e r , designated e^ and f ^ r e s p e c t i v e l y , represent the angles between the ray path and the v e r t i c a l a x i s ( F i g u r e 1). The f r e e s u r f a c e displacement, the i n t e r m e d i a t e boundary displacements and the displacement waves propagated downwards i n t o the basement generated by a plane P or S wavefront i n c i d e n t on the c r u s t from the basement are sought. I n t h i s chapter the normal i n c i d e n c e wave propagation problem w i l l be formulated f i r s t and then g e n e r a l i z e d to the non-normal problem. This approach should most simply demonstrate the e s s e n t i a l nature of the f o r m u l a t i o n without obscuring the method i n the d e t a i l s of the ge n e r a l problem. Time domain system d e s c r i p t i o n s w i l l then be d e r i v e d from these frequency f o r m u l a t i o n s . 2.2 L i n e a r Systems Approach The l i n e a r system f o r m u l a t i o n d e s c r i b e s wave propagation i n terms of f i l t e r or communications theory. In t h i s manner, the c r u s t can be regarded as a f i l t e r system or communications channel which r e c e i v e s an 7 LAYER 1 LAYER 2 ° 2 P2 v2 • • • LAYER n fln 7 LAYER n + 1 y ( BASEMENT )y A e n + l / (^ n + 0 / / F i g . 1. T h e h o r i z o n t a l l y m u l t i l a y e r e d medium w i t h n l a y e r s o v e r l y i n g a b a s e m e n t h a l f s p a c e . E a c h l a y e r i i s s p e c i f i e d b y P - and S-wave v e l o c i t i e s ct^ a n d B i ; d e n s i t y p-j_ and t h i c k n e s s n^. T h e p l a n e wave f r o m t h e b a s e m e n t i s i n c i d e n t a t a n g l e e , , f o r P wave ( o r f , . f o r S w a v e ) . ° n+1 n+1 i n p u t s i g n a l a n d o p e r a t e s o n i t t o c r e a t e t h e o u t p u t s i g n a l ( o r t h e s i g n a l s e t i n t h e c a s e o f many o u t p u t s ) . F o r c r u s t a l t r a n s m i s s i o n t h e i n p u t i s a t e l e s e i s m i c P o r SV w a v e m o t i o n i n c i d e n t on t h e c r u s t a t t h e b a s e m e n t b o u n d a r y ; t h e o u t p u t s i g n a l s e t c o n t a i n s t h e c a r t e s i a n c o o r d i n a t e s o f t h e f r e e s u r f a c e m o t i o n s . T h e w a v e f o r m s r e - t r a n s m i t t e d i n t o t h e b a s e m e n t a n d a l l i n t e r n a l b o u n d a r y m o t i o n s a r e a l s o a v a i l a b l e . T h e s e i n p u t - o u t p u t r e l a t i o n s c a n b e b l o c k d i a g r a m m e d as i n F i g u r e 2. U n d e r t h e a s s u m p t i o n t h a t t h e c r u s t a l r e s p o n s e i s l i n e a r , t h e r e s p o n s e o p e r a t o r i s i n d e p e n d e n t o f t h e i n p u t w a v e f o r m . T h u s , i f two 8 i(t) INPUT SIGNAL OUTPUT S SIGNAL SET CRUST RESPONSE Fig. 2. The input-output diagram for a linear system crust. The output signal set O j ( t ) is obtained from the convolution of the input signal i(t) with the crustal response vector hj(t). among the input, the output set and the response operator are known, the third is uniquely determinable. In particular, i f the input signal and response operator are known, the output signals can be calculated. A f i l t e r system is often described in terms of the output signal set generated by specific input signals such as the Dirac delta function <5(t). This particular response is called the impulse-response function set of the system, h.(t). The output signal set o.(t) is J 3 obtained by the convolutions of the input signal i(t) with the impulse response function vector h^(t): i C t ' ^ C t - t ' J d t ' (2-1) = l(t)*h (t) For a l i n e a r system, an equivalent equation i n terms of complex frequency components may be obtained by the Laplace i n t e g r a l transformation such that: Cv.(s) = I(s).H (s) (2-2) where s = o+ito i s the complex frequency and the functions of the v a r i a b l e s designated with c a p i t a l s are the Laplace transformed time functions, o GO e.g. I(s) = L { i ( t ) } = | i ( t ) e ~ S t d t . (2-3) o T h i s • r e p r e s e n t a t i o n r e l a t e s the transformed input s i g n a l to the transformed output s i g n a l set by m u l t i p l i c a t i o n with the system t r a n s f e r function v e c t o r at each complex frequency s. Both equations completely characterize the system and the p a r a l l e l functions i n the time and frequency domains are generally obtainable from each other through the Laplace transform and i t s i n v erse. 2.3 The Linear System Analog f o r Normal Incidence The l i n e a r system analog f o r the simple case of normal incidence i s now constructed to show the e s s e n t i a l concepts. This r e s t r i c t i o n eliminates any conversions between compressional and r o t a t i o n a l motions at the l a y e r boundaries and l i m i t s the analog to two equivalent s i g n a l channels f o r each of the up and down-travelling P (or S) waveforms. The m u l t i l a y e r l i n e a r system model may be developed i n terms of i n d i v i d u a l layers as follows. On the boundaries of any one l a y e r i , four p o s s i b l e input-output s i g n a l s represent the up and down-travelling waveforms; the s i g n a l s u^(t) and d^(t) represent the input displacement waveforms to l a y e r i , while u . _ , ( t ) and d, . ( t ) represent the output 10 LAYER i-I Ui_ j ( t ) LAYER i • J - , ! l U j ( t ) dj (t) 2 = 0 LAYER i+l dj + |(t) Fig. 3. A representation of the possible input and output waveforms in layer i under normal incidence. waveforms from layer i and the equivalent input waveforms to layers i - l and i+l respectively (Figure 3). Between the boundaries, the up-and down-travelling waveforms are not coupled to each other because the propagation characteristics are assumed to be linear. However, these waveforms are linked at the boundaries by the standard displacement wave reflection and transmission coefficients. The four possible input-output signals are related through four impulse response functions, h. (t), h (t), h. (t) and h (t),representing the layer. These functions xl 12 x3 itt must be determined from the layer parameters. Convolution of the input 11 s i g n a l s and impulse response f u n c t i o n p a i r s y i e l d s the output s i g n a l s as f o l l o w s : u 1 _ 1 ( t ) = h l i ( t ) * u i ( t ) + h i 2 ( t ) A d i ( t ) (2-4) d (t)=h (t)*u.(t)+h. ( t ) * d f t ) . 2.4 The Frequency Domain D e s c r i p t i o n - ' I t i s convenient to determine these input-output r e l a t i o n s i n the frequency rather than the time domain by the Laplace (or Fourier) transform methods. The f i r s t equation above under the Laplace t r a n s -formation becomes: D i - l ( s ) = H i 1 ( s ) U i ( s ) + H i 2 ( s ) D i ( s ) (2-5) D (s)=H (s)U.(s)+H (s)D.(s). 1+1 13 .1 it* 1 In t h i s representation each impulse response function corresponds to a t r a n s f e r f u n c t i o n : H (s) = L{h. (t)} (2-6) k \ and each time s i g n a l corresponds to a transform s i g n a l : U ±(s) = L{u.(t)} . (2-7) ^ T h e existence of these convolution i n t e g r a l s demands that the t r a n s f e r functions H-^_(s) are s e c t i o n a l l y continuous. Any p h y s i c a l l y reasonable choice of r e f l e c t i o n and transmission c o e f f i c i e n t s and the attenuation functions A^(s) w i l l f u l f i l l these requirements. 12 F i g . 4. A ray path representation of the i n t e r n a l reverberations within layer i giving r i s e to output transform signals U ^ ^ C s ) a n d D i + i ^ s ^ ' T h e transmission and r e f l e c t i o n c o e f f i c i e n t s for the top and bottom boundaries are t ^ and r ^ and t i ' and r^' respectively. This diagram represents the normal incidence condition; the diagram i s expanded h o r i z o n t a l l y for c l a r i t y . Direct c a l c u l a t i o n for the layer transfer function i s possible on the basis of ray theory by the following procedure. Consider a l l reverberations due to an input P (or S) wave within a layer i (Figure 4). The transfer function H. (s) r e l a t i n g input s i g n a l U.(s) to output si g n a l ^ ^ _ ^ ( s ) m a y be determined as a series i n which the f i r s t term represents the d i r e c t path transmission and each advancing term represents an increasing order of i n t e r n a l reverberation (1,2,3,....): 13 -sT. -3sT H ± (s) = t 1 { A i ( s ) e 1+A J L 3(s)e- Vr.* (2-8) - 5 s T + A 5 ( s ) e r 2r ,2+...} where the displacement wave transmission and r e f l e c t i o n c o e f f i c i e n t s at the upper and lower boundaries are represented by t . . r. and t , ' and r . ' 1 1 X 1 r e s p e c t i v e l y . A^(s) represents an a r b i t r a r y attenuation function of the complex frequency per s i n g l e t r a n s i t through the l a y e r and-T i s the t r a v e l time through the la y e r determined by the wave phase v e l o c i t y and l a y e r thickness. The exact l i m i t i n g sum of t h i s s e r i e s i s : -sT. A,(s)e 1 • • H (s) = * — C i (2-9) ± 1 „ ~ 2 s T i 1-A 2 ( s ) e r r ' G(s) . t (2-10) l-F(s)G(s) 1 V -sT where G(s) = A i ( s ) e 1 and F(s) = G ^ r ^ r ' S i m i l a r l y , the three remaining t r a n s f e r functions f o r l a y e r i are: V> - • V i ' <2-U) H (s) = G 2(s) . t 'r ( 2 _ 1 2 )  ± 3 l-F(s)G(s) 1 1 U ' 14 H (s) = G(s) . t 1 . (2-13) X * l-F(s)G(s) 1 These i n f i n i t e s e r i e s l i m i t i n g sums consider a l l p o s s i b l e i n t e r n a l rever-b e r a t i o n and ray paths within a l a y e r for the normal incidence condition. The Laplace i n v e r s i o n of these t r a n s f e r functions to obtain an equivalent set of time domain impulse response functions i s discussed i n Appendix I . A p a r a l l e l formulation for S waves requires only the appropriate choice of the shear wave phase v e l o c i t y and r e f l e c t i o n and transmission c o e f f i c i e n t s . G(s) The evident analogy of the form -, \ to the simple l i n e a r b J l-F(s)G(s) r p o s i t i v e feedback network (Bohn, 1963) leads to a f u n c t i o n a l block representation of the l i n e a r system analog for the l a y e r i (Figure 5). This system exactly diagrams a l l four l a y e r t r a n s f e r functions. For example, the input s i g n a l U\ (s) enters the system at the lower boundary and c i r c u l a t e s throughout the system as i n d i c a t e d by the arrows. The r e s u l t i n g output response s i g n a l ^(s) becomes the product of t h i s input s i g n a l and the exact t r a n s f e r function ( s ) . S i m i l a r l y , the three remaining t r a n s f e r functions of l a y e r i are described i n the s i n g l e diagram. Any number of such networks, each representing one l a y e r can be combined to construct a t o t a l m u l t i l a y e r crust l i n e a r system network. To completely determine the problem, input s i g n a l s and free surface boundary conditions are required. The necessary input s i g n a l co n d i t i o n r e l a t e s the i n c i d e n t waveform coming from the basement to the s i g n a l u (t) i n the bottom layer n: n U n ( t ) = " I N P U T ^ ' V l • <2-14> 15 Uj(s) - 0 Aj(s)e-sTi D j + , (s ) Aj(s)e-sTi Uj_,(s) t - l Di ( s i y LAYER i+l LAYER i LAYER i -1 F i g . 5. A f u n c t i o n a l block representation of the l i n e a r system analog for la y e r i under normal incidence. The heavy v e r t i c a l l i n e s denote l a y e r boundaries. The corresponding transform s i g n a l s are also s i m i l a r l y r e l a t e d : V S ) = U I N P U T ( s ) , t : n + l ' (2-15) The free surface boundary condition i s included by the choice of the transmission c o e f f i c i e n t ti=0 and the s i g n a l di(t)=0 f o r la y e r 1. These 16 two conditions imply that no wavemotion energy i s l o s t or gained at the free surface boundary. For the case of n layers above a basement halfspace, a set of 2n complex l i n e a r simultaneous equations i n terms of the i n d i v i d u a l layer t r a n s f e r f u n c t i o n sets and the si g n a l s transforms represents the system: U!(s) = H 2 (s)U 2(s)+H 2 (s)D 2(s) 1 2 U 2(s) = H 3 (s)U 3(s)+H 3 (s)D 3(s) 1 2 D 2(s) = Hi (s)U!(s) 3 U 3(s) = H t t i(s)U I +(s)+H l t 2(s)D 4(s) (2-16) D (s) = H ,(s)U ,(s)+H ,(s)D (s) n n - l 3 n - l n-14 n - l U n ( s ) = W U I N P U T ( 8 ) V l ( s ) = V ( S ) D n ( s ) + H n 3 ( s ) U n ( s ) + r n + l - U I N P U T ( s ) The s o l u t i o n of the problem requires the evaluation of each U^(s) and D^(s). Then, the desired surface v e r t i c a l displacement amplitude trans-form W Q(s) i s r e l a t e d to the s i g n a l Uj(s) by W Q(s) = H o(s)U!(s) (2-17) 17 H q(S) may be obtained from the block diagram f o r l a y e r 1 ( r e f e r to Figure 5 ) noting that t i = 0 , r j = - l and D 1 ( s ) = 0 . The surface displacement s i g n a l , W Q(S), i s the sum of si g n a l s due to the i n c i d e n t and r e f l e c t e d waveforms at the free surface boundary. Since the r e f l e c t i o n c o e f f i c i e n t is -1, the two si g n a l s have opposite signs and d i r e c t i o n s so that the sum is twice the incident s i g n a l displacement amplitude. This i n c i d e n t -sTi s i g n a l i s formed by Uj(s) which has been delayed ( m u l t i p l i c a t i o n by e l ) , attenuated ( m u l t i p l i c a t i o n by A i ( s ) ) and reverberated (divided by —2sT 1+A 1 2(s)r 1e l) so that: 2 A i ( s ) e " s T l H0 < s ) = • ( 2 ~ 1 8 ) l + A 1 2 ( s ) r 1 e _ 2 s T In cases i n v o l v i n g only one or two l a y e r s , a s o l u t i o n by d i r e c t a l g e b r a i c manipulations on the matrix of the f u n c t i o n a l c o e f f i c i e n t s i s p o s s i b l e . However, f o r many l a y e r s , the a n a l y t i c a l s o l u t i o n s f o r the transform surface motion and r e f l e c t e d s i g n a l s and the subsequent Laplace i n v e r s i o n to give the output time domain response functions i s not f e a s i b l e . In l i e u of the a n a l y t i c a l approach, numerical s o l u t i o n s of the simultaneous equation set at incremented F o u r i e r frequencies (obtained from the Laplace frequencies by s e t t i n g s=ico) may be accomplished to give a steady s t a t e frequency response function for the multilayered medium. This more l i m i t e d a nalysis corresponds to the H a s k e l l matrix numerical computation f o r s p e c t r a l responses. Results of such numerical analyses w i l l be presented i n a l a t e r chapter. 18 2.5 The Linear System Analog for Non-normal Incidence The formulation f o r non-normal incidence plane waveforms proceeds analogously to that f or normal incidence. However, immediate complication a r i s e s i n the problem of handling conversions between compressional and r o t a t i o n a l motions which occur at each l a y e r boundary. Because such conversions occur, the independent P and SV motions must be considered separately r e q u i r i n g a four-channel rather than a two-channel analog approach (Figure 8). For normal incidence, the waveform ray paths d i r e c t l y correspond to the p o s i t i v e and negative c a r t e s i a n z directions., but for non-normal conditions, the ray paths are not coincident with the coordinate axes. Consequently, i t i s convenient to create a set of pseudo-coordinates i n which the separate P-and S-wave motions e f f e c t i v e l y comprise the orthogonal vector p a i r (Figure 6). In each l a y e r , the d i r e c t i o n s along the ray paths determine the required pseudo-coordinate set and are r e l a t e d to the c a r t e s i a n coordinates by the angles between the normal and the P and S rays. I t i s evident that a d i f f e r e n t transformation i s required f o r every l a y e r . With h o r i z o n t a l l a y e r i n g , the compressional and shear motions throughout the medium possess a constant apparent h o r i z o n t a l phase v e l o c i t y c determined by the input conditions. For example, f o r an input P-waveform from the basement: c = c x n + 1 / s i n ( e n + 1 ) (2-19) where i s the P ray normal angle and °< n +^ is the P phase v e l o c i t y i n the basement. S n e l l ' s law then determines the P and S ray angles with the normal i n any other l a y e r i : c = a i / s i n ( e i ) = Bj/sinU ) . (2-20) 19 Up- and down-travelling waveforms are thus represented by the positive and negative directions of the two axes. The four possible directions on the two pseudo-coordinate axes determine the four required s i g n a l channels. Fi g . 6 . A representation of the "pseudo-coordinate" transformation for layer i . The arrows on the rays of the cartesian representation show the p o s i t i v e displacement conventions used throughout this work. The pseudo-coordinate representation i s for layer i only. The second diagram shows the P and SV displacements i n layer i transformed into the pseudo-coordinate representation. 2 0 Z NODES OF REFRACTED WAVE NODES OF INCIDENT WAVE F i g . 7. A diagram showing the non-coincidence of wave constant phase p o s i t i o n s (e.g. nodes) and a h o r i z o n t a l l a y e r boundary. A f u r t h e r complication a r i s e s f o r non-normal incidence because the phase of a plane wave i s not constant along a boundary or along any surface of constant depth (Figure 7 ) . For each l a y e r i , the i n t r o d u c t i o n of phase advance (negative time delay) t r a n s f e r functions with the forms S T S . S T P . e f o r SV motions and e for P motions (Figure 8) to the l i n e a r system analog r e - e s t a b l i s h the corre c t phase r e l a t i o n s f o r the s i g n a l s at the boundaries. For a d e s c r i p t i o n and c a l c u l a t i o n of these phase advance t r a n s f e r functions see Appendix I I . Just as each l a y e r i n the normal incidence problem was character-i z e d by two equations i n layer t r a n s f e r functions and s i g n a l s , each layer /-vST*l-\ . e Pj A p . (s)e s T P j Fig. 8 . A functional block representation for non-normal incidence for layer i . Note that four channels are shown, two for the up- and down-travelling P-waves and two for the up- and down-travelling SV waves. 22 in the non-normal problem is represented by four such equations; the transfer function relations for a layer i in a matrix form i s : u i - l ( s ) " H (s)H (s)H (s)H (s)" l j i 2 1 3 lu U ±(s) v i - l ( s ) • # • • • V ±(s) D 1 + 1 ( a ) * a • • D i(s) E i + 1 ( s ) * • • H . (s) E i(s) (2-21) where the transforms U\(s) and D^(s) correspond to the up- and down-travelling P waves, the V^(s) and E^(s) correspond to the up- and down-travelling SV waves and the H (s) are the transfer function set for the •Me layer. A direct algebraic determination of these transfer functions i s impractical. However, i t is possible to construct directly the layer linear systems block diagram representation (Figure 8) by analogy to the normal incidence system and simplify i t using "block diagram reduction procedures" (Kuo, 1962). One example reduction i s performed in Appendix III. For frequency component analysis, such reductions are best accomplished during numerical analysis; for time domain analysis, reductions are unnecessary. The equivalent linear system block representation of one layer for non-normal ray conditions contains four channels interconnected at the boundaries by the conversion reflection and transmission coefficients calculable using the Zoeppritz' displacement relations (Richter, 1958). For ray angles less than the c r i t i c a l angle, the sign of the real valued Zoeppritz' coefficients represents phase changes at refractions and r e f l e c t i o n s . S u p e r c r i t i c a l angles between the ray d i r e c t i o n s and the boundary normal angles determine complex transmission and r e f l e c t i o n c o e f f i c i e n t s representing phase changes other than 0 and T T. The l i n e a r system approach described i n t h i s t h e s i s has not been extended to permit such s u p e r c r i t i c a l conditions. Thus, trapped wave phenomena are not described. In the f u n c t i o n a l block diagram of Figure 8, the su b s c r i p t i re f e r s to the l a y e r i . The s i n g l e s u b s c r i p t s p and s designate the separate channels f o r the P and SV waveforms; each f u n c t i o n a l block so i n d i c a t e d pertains only to the corresponding wave type. For example, A (s) i s the t r a n s f e r function representing the attenuation of a P wave during one t r a n s i t through l a y e r i . The double subscripts of the Zoeppritz' c o e f f i c i e n t s i n d i c a t e conversion from the waveform type i n d i -cated by the f i r s t s u b s c r i p t to the waveform type i n d i c a t e d by the second. For example, t i s the transmission c o e f f i c i e n t which determines the s p ± conversion of S waves from l a y e r i i n t o P waves i n l a y e r i - l ; r . i s PP i the r e f l e c t i o n c o e f f i c i e n t which determines the P wave to P wave r e f l e c t i o n at the lower boundary of the la y e r i . The waveform types f o r the separate channels are i n d i c a t e d by P and S symbols. A l l other symbols r e t a i n the d e f i n i t i o n s of the normal incidence case. An n-layered non-normal incidence problem requires the s o l u t i o n of 4n complex l i n e a r simultaneous equations formed s i m i l a r l y to those, discussed i n the normal incidence problem: 24 Ui(s) = H 2 U 2(s)+H 2 V 0(s)+H 2 D 2(s)+H 2 E 2(s) 1 2 3 h V1(s) = H 2 U 2(s)+H 2 V 2(s)+H 2 D 2(s)+H 2 E 2(s) 5 6 7 8 D 2(s) = H : U1(s)+H, Vj(s) 9 10 E 2(s) = Hx U 1(s)+H 1 V^s) 13 14 U 2(s) = H 3 U 2(s)+H 3 V 2(s)+H 3 D 2(s)+H 3 E 2(s) 1 2 3 V 2(s) = H 3 U 2(s)+H 3 V 2(s)+H 3 D 2(s)+H 3 E 2(s) 5 6 7 8 D (s) =H U(s)+H V (s)+H D (s)+H E (s) n n - l g n n - 1 l o n 1 n - ± l l n - i n ~ M 2 n _ 1 E (s) = H . U ,(s)+H V ,(s)+H . D (s)+H E (s) n n - 1 ^ 3 n - l n— 1 ^ n - l n—±15 n - l n - l ^ g n - l V s ) * %p n + 1 UINPUT ( s ) VS> = ' p s ^ I N P U T ^ (2-22) D .(s) = H U (s)+H V (s)+H D (s)+H E (s) n + l ng n n^g n . n l l n n l 2 n + rpp n + 1 UINPUT ( s ) E . (s) = H U (s)+H .V (s)+H D (s)+H E (s) n+1 iij3 n n ^ n 1115 n n 15 n + rps n + 1 UINPUT( S) 25 The necessary free surface boundary conditions are included by the free surface Zoeppritz' r e f l e c t i o n c o e f f i c i e n t s and the requirement that no wavemotions enter from the vacuum halfspace ( i . e . (s)=E 1(s)=0). The required input conditions r e l a t e the incident s i g n a l from the basement to the input s i g n a l s at the lower boundary of the bottom la y e r n. For example, the input condition transforms generated by an input SV wave are: ™ V S ) - S p - VINPUT ( S ) n+1 .... ( 2 - 2 3 ) V S ) • S s - VINPUT ( S ) n+1 A s o l u t i o n requires evaluation of a l l s i g n a l s U^(s), D^(s), V^(s) and E^(s) which are u n s p e c i f i e d by e i t h e r the boundary or input conditions. A f o l l o w i n g f u n c t i o n a l operation then r e l a t e s the P and SV motions incident on the free surface boundary to the s i g n a l s U^(s) and V^(s) obtained from the s o l u t i o n of the 4n simultaneous equations: U (s) = H (s)U 1(s)+H ( s ) V 1 ( s ) o o j ° 2 V (s) = H (s)U 1(s)+H r t ( s ) V 1 ( s ) o O3 1 Oi+ 1 (2-24) where U (s) i s the P and V (s) the SV surface incident motions. The set o o of H (s) are obtained by reduction procedures on the f u n c t i o n a l block T2] I f S waves with a r b i t r a r y p o l a r i z a t i o n determine the incident waveform, the"two dimensional analysis i s not s t r i c t l y v a l i d . However, any S motion can be separated i n t o SV ( v e r t i c a l p o l a r i z a t i o n ) and SH ( h o r i z o n t a l p o l a r i z a t i o n ) components. The response to each component can be c a l c u l a t e d separately. 26 F i g . 9. The free surface r e f l e c t i o n f o r non-normal incidence. U Q and V Q are the in c i d e n t P- and S-waves and and V R are the r e f l e c t e d P- and S-waves. representation of la y e r 1 considering the free surface boundary conditions i n a manner s i m i l a r to the normal incidence case. Then, the surface motion transforms, v e r t i c a l component W Q ( S ) and h o r i z o n t a l component Y Q ( S ) , are obtained by re-combination of the displacement components i n the c a r t e s i a n z and x d i r e c t i o n s due to the P and S V waveforms inc i d e n t on and r e f l e c t e d from the free surface (Figure 9). The r e f l e c t e d transform ' s i g n a l s , U (s) and V ( s ) , are r e l a t e d to the surface i n c i d e n t s i g n a l s by the Zoeppritz' free surface r e f l e c t i o n c o e f f i c i e n t s : 27 s p i (2-25) The surface displacement amplitude c o e f f i c i e n t s are then: W o ( s ) = [ U o ( s ) + U R ( s ) ] c o s ( e 1 ) + [ V o ( s ) + V R ( s ) ] s i n ( f 1 ) (2-26) Y o ( s ) = [ U o ( s ) - U R ( s ) ] s i n ( e 1 ) - [ V o ( s ) - V R ( s ) ] c o s ( f 1 ) Because the a l g e b r a i c forms of the t r a n s f e r functions are so complex, even the simple one l a y e r system i s not r e a d i l y soluble a n a l y t i c a l l y . Again, numerical methods with incremented r e a l frequencies are most e a s i l y used.to determine the s p e c t r a l response of the medium. 2.6 Time Domain Formulation A development of the time domain systems d e s c r i p t i o n i s now derived from the p a r a l l e l frequency d e s c r i p t i o n f o r normal incidence conditions. The g e n e r a l i z a t i o n to permit non-normal incidence i s not described e x p l i c i t l y , but the complete system analog i s presented. The transform approach described e a r l i e r i s o f t e n used i n the a n a l y s i s ^ o f system response c h a r a c t e r i s t i c s . However, f o r d i g i t a l analyses, transform methods can be l e s s e f f i c i e n t ( i n terms of computer time and storage requirements) i f l i t t l e convolution i s necessary f o r a d i r e c t time domain s o l u t i o n . e x i s t . F i r s t , the system of l i n e a r simultaneous equations i n terms of the s i g n a l transforms and l a y e r t r a n s f e r functions (Equations 2-16, 2-22) Two d i r e c t approaches to a time domain l i n e a r systems s o l u t i o n can be transformed into the equivalent system of l i n k e d convolution i n t e g r a l s . This system of equations i n convolutions may be solved by standard numerical procedures for s o l u t i o n of i n t e g r a l equations using d i g i t a l computers. (e.g., Carnahan et a l , 1969). However, the approach has some drawbacks as each l a y e r requires s i x t e e n equivalent convolutions i n the general non-normal case. For non-attenuating conditions, each of the s i x t e e n required impulse responses f o r each l a y e r i s an i n f i n i t e impulse t r a i n with decreasing impulse amplitudes with i n c r e a s i n g time (Figure 10). The rate of the decrease of amplitude depends upon the r e f l e c t i o n c o e f f i c i e n t s at the two l a y e r boundaries. E f f i c i e n t com-putation using d i g i t a l methods requires f i n i t e (and p r e f e r a b l y short) convolutions which necessitates the truncation of the i n f i n i t e duration impulse responses. However, f o r accuracy, the truncated impulse responses must be at l e a s t as long as the time required f o r s e v e r a l waveform t r a n s i t s through the l a y e r i n order that the reverberation properties of each l a y e r are w e l l characterized. With the i n c l u s i o n of attenuation the l a y e r impulse response functions are no longer simple impulse t r a i n s but rather each impulse i s broadened by the p r e f e r e n t i a l low-pass nature of attenuation and the attending d i s p e r s i o n (Figure 10). This broadening increases with time at a rate determined by the l a y e r Q. The second approach which i s computationally more e f f i c i e n t and used i n t h i s work accomplishes a time domain formulation i n terms of system block representations by analogy to the frequency d e s c r i p t i o n s of Figures 5 and 8. Each frequency domain system block t r a n s f e r function corresponds to an equivalent time function which can be d i r e c t l y obtained from the frequency function through the inverse Laplace transform. The time domain system d e s c r i p t i o n s of Figures 11 and 12 have been obtained 29 IMPULSE RESPONSE Q=00 Q<00 F i g . 10. A comparison of a normal incidence l a y e r impulse response f o r i n f i n i t e Q and causally .absorbing conditions. i n t h i s way. For example, f o r la y e r i at normal incidence the frequency _ s T i domain delay and absorption t r a n s f e r function a^(s)e (Figure 5) i s inverse Laplace transformed to give the equivalent time domain delay and absorption function 5(t-T )*a^(t) (Figure 11) where a^(t) describes the impulse response of the layer per s i n g l e t r a n s i t determined by the la y e r 30 LAYER i + l LAYER i LAYER i - l F i g . 11. The time-domain representation of the system analog of one l a y e r at normal incidence. Compare to the frequency representation of Figure 5. absorption p r o p e r t i e s . (The form of t h i s absorption f u n c t i o n i s discussed l a t e r i n t h i s chapter). To describe the e s s e n t i a l nature of the time domain d e s c r i p t i o n , the normal incidence formulation (Figure 11) i s now f u r t h e r demonstrated. At. the lower boundary of layer i , a waveform u ^ ( t ) enters the system. This waveform combines by add i t i o n to any waveform u (t) which r e s u l t s Fig. 12. The time-domain representation of the system analog of one layer for non-normal incidence. Compare to Figure 8. 32 from the r e f l e c t i o n of any down-travelling wave at the lower boundary. u^(t) i s the output s i g n a l from the system block r ' which represents the lower boundary r e f l e c t i on c o e f f i c i e n t . The sum waveform i s then delayed by the delay operator 6(t-T^) and absorbed according to the f u n c t i o n a^(t) during i t s t r a n s i t to the top boundary so that the equivalent waveform incident on the top boundary i s : y i ( t ) = ( u 1 ( t ) + u r ( t ) ) * 6 ( t - T 1 ) * a i ( t ) (2-27) = ( u i ( t - T i ) + u r ( t - T i ) ) * a i ( t ) . I f the l a y e r possesses no i n t e r n a l attenuation: a t ( t ) = 1 , (2-28) and y ± ( t ) = u ^ t - T ^ + u ^ t - T ^ . (2-29) The remaining system blocks are merely constants describing the other r e f l e c t i o n and transmission c o e f f i c i e n t s at the l a y e r boundaries and are i d e n t i c a l to those of the frequency domain system. For example, the waveform y i ( t ) incident at the top boundary of the layer i s p a r t i a l l y transmitted i n t o the next layer through the transmission c o e f f i c i e n t t^ to form u ^ _ 2 ^ ^ : U i - l ( t ) = V Y i ( t ) * ( 2 ~ 3 0 ) y ^ ( t ) i s also p a r t i a l l y r e f l e c t e d back i n t o l a y e r i to form d^(t) = r£ ,y ;j_( t) which combines by a d d i t i o n to any input down-travelling waveform d ^ ( t ) . 33 This combined waveform i s delayed and absorbed during i t s t r a n s i t to the lower boundary to form: x ± ( t ) = ( d i ' ( t - T 1 ) + d r ( t - T 1 ) ) * a i ( t ) . ~ (2-31) x^(t) i s p a r t i a l l y r e f l e c t e d at the lower boundary to form u^(t) = r^'»x^(t) which again combines by addition to u^(t) entering the l a y e r at the lower boundary. < For normal incidence, only two convolutions with the attenuation and delay operator are necessary f o r each l a y e r ; four are required for non-normal incidence. I f attenuation i s not present, the form of these convolutions s i m p l i f i e s to mere delays of the input waveforms. Then only the past h i s t o r i e s of the input function must be retained f o r a duration equal to the t r a n s i t time. Two l i n k e d i n t e g r a l equations thus describe the normal l a y e r responses: v U i - l ( t ) = t i [ u i ( t _ T i ) + , d i + 1 ( t ) ] * a . ( t ) t i ' (2-32) r i d . + 1 ( t ) = t . ' [ d i ( t - T . ) + u . _ 1 ( t ) ] * a i ( t ) ; a s i m i l a r set of four equations describe the non-normal case. For d i g i t a l time simulations, each operator and s i g n a l within the system model must be time incremented. Then, for time incremented sol u t i o n s (simulations) of the n-layered crust response to an input waveform, a set of 2n l i n k e d equations of the above form must be solved. 34 (4n such equations describe non-normal incidence). As i n the frequency formulation, the necessary input conditions are determined by the input s i g n a l at the basement boundary; the correct free surface r e f l e c t i o n c o e f f i c i e n t s determine the required boundary condit i o n s . For the d i g i t a l s i m u l a t i o n , i t i s important to choose the constant time increment to be small enough as to prevent a l i a s i n g of the system frequency response or the s p e c t r a l character of the input waveform due to sampling. 2.7 Attenuation and Dispersion i n the System Model E l a s t i c wave attenuation i n an absorbing medium i s u s u a l l y described i n terms of the m a t e r i a l Q f a c t o r (Knopoff, 1964). For a Q constant at a l l frequencies, a s p a t i a l absorption f a c t o r a(w) = £Q^" determines the decreasing amplitude of an harmonic wave with distance: U(x,co) = U(o,co)exp(-a(to)x) . (2-33) I f the absorption function i s r e a l and the phase v e l o c i t y c i s independent of frequency a non-physical acausal attenuation law r e s u l t s . Futterman (1962) extended t h i s d e s c r i p t i o n of attenuation to allow the p r i n c i p l e of c a u s a l i t y to hold while r e t a i n i n g the e s s e n t i a l l i n e a r nature of the absorption (Q remains, constant i n frequency) and wave motion i n order that s uperposition remains v a l i d . Subsequent analyses by Wuenschel (1965) demonstrated that Futterman's theory w e l l described attenuation i n r e a l m a t e r i a l s . Under the condition of causal attenuation, a plane wave having t r a v e l l e d through distance x can be represented by a superposition i n t e g r a l i n terms of the o r i g i n F o u r i e r c o e f f i c i e n t s : 35 f r 31 u(x,t) = U(o,w)exp(i[K(u>)x-ut])du> , 1 J (2-34) where K(oi) = k(a))+ia(u)) . K(oj) also describes the disp e r s i o n properties of the m a t e r i a l . For Futterman's model ( A 1 , D 1 ) which i s used i n l a t e r analyses K(w) has the form: K(w) = •< o to c tii £ CJ (2-35) where tu^ i s an a r b i t r a r y c u t o f f angular frequency below which Q can be considered to be i n f i n i t e . For frequencies w e l l above w^ , Futterman shows that the dispersed wave phase V(OJ) and group u(to) v e l o c i t i e s are r e l a t e d to the l i m i t i n g undispersed phase v e l o c i t y c by: v(«o) S u(u) S c(l- lnC^-))" 1 (2-36) c, thus describes a lower v e l o c i t y l i m i t on an inverse d i s p e r s i o n c o n d i t i o n . The absorption t r a n s f e r functions and impulse responses which are required by the l i n e a r systems formulation are described by t h i s attenuation model. For the time domain system: •[.3] The i n t e g r a l representation f o r Fourier superposition i n t h i s case i s that normally used i n wave propagation problems. Throughout the re s t of t h i s t h e s i s , the communications theory convention i s used. 36 a±(t) = | exp(i[K i((o)n i-ut])du (2-37) —CO where is the thickness of layer i and K^(OJ) the complex propagation function determined by the layer Q, phase velocity and the cutoff frequency O) q. This a^(t) i s the Futterman model impulse response (i.e. U(o,co)=l) of the layer. The equivalent frequency description (for the frequency domain linear system analog) in Fourier terms i s : . Ai(a>) = exp(iK i(o))n i) . <2"38> If the causal property for attenuation i s not required (i.e. no phase information is desired) then the Futterman absorption function can be simplified to the usual acausal form: -uT. Ai(co) - e x p ( - ^ ) (2-39) where i s the layer Q value and T^ the layer transit time. It should be noted that the Futterman analysis of attenuation requires cu^  and c. In practice, however, the dispersed phase velocity v(to) and group velocity u(to) are the observable quantities while O)q and c must be obtained by calculation. Because the form of the absorption impulse response function a^(t) is only weakly dependent upon the actual cutoff frequency, can be chosen a r b i t r a r i l y to be very small compared to the frequencies of interest. Futterman's choice of to =10"3 s e c - 1 has ^ o been used i n the computation of the absorption functions used later in this thesis. 37 In the preceding discussion, the absorption t r a n s f e r function, A_^(a)), and the absorption time function, a_^(t), can represent e i t h e r P-wave or S-wave attenuation. I t i s only necessary to speci f y the p a r t i c u l a r t r a n s i t time and Q corresponding to the wave type to determine . the frequency or time function desired. The t r a n s m i s s i o n - r e f l e c t i o n conditions at a f l a t boundary between two causally absorbing media are affe c t e d by the coincident dispersion of the e l a s t i c waves. Because the v e l o c i t y r a t i o s between the two media are not n e c e s s a r i l y constant at a l l frequencies, each Fourier component observes transmission and r e f l e c t i o n c o e f f i c i e n t s which are d i f f e r e n t from the non-dispersed case. The v e l o c i t y r a t i o s also determine the ray di r e c t i o n s of the transmitted frequency components by Snell's Law. I f the r a t i o s are not constant, the transmitted frequency component ray d i r e c t i o n s are not a l l coincident and a geometrical condition on.dispersion r e s u l t s . In a frequency domain a n a l y s i s , t h i s dispersion property i s not a major problem because transmission and r e f l e c t i o n c o e f f i c i e n t s and ray d i r e c t i o n s can be computed for each Fourier component frequency. However, i n a time domain approach, the separation of frequency components cannot be accomp- • l i s h e d and the method becomes an approximation of the true s o l u t i o n . Wuenschel's (1965) analysis of attenuation and dispersion in the P i e r r e Colorado shales suggests that the frequency dependence of the dispersed phase v e l o c i t y i s small over frequency ranges of several octaves.: Con-sequently, a constant t r a n s m i s s i o n - r e f l e c t i o n c o e f f i c i e n t and ray d i r e c t i o n approximation appears to be v a l i d and i s assumed i n the following analyses.' • 38 CHAPTER III ANALYSIS OF THE ALBERTA SEISMOGRAMS 3.1 Introduction In order to demonstrate the unique features and u t i l i t y of t h i s theory, the re s u l t s of analyses of seismograms recorded i n ce n t r a l A l b e r t a w i l l be compared to l i n e a r systems solutions i n both time and frequency domains. The e f f e c t s of c r u s t a l attenuation on s p e c t r a l r a t i o s w i l l be shown to be s i g n i f i c a n t f o r the low Q values which are thought to characterize c r u s t a l sediments. I t w i l l be observed that the c o r r e l a t i o n 4 between experiment and theory i s best at frequencies below 1 Hz. Time domain synthetic seismograms w i l l also be cal c u l a t e d by the systems approach and compared with experimental recordings. The e f f e c t s of c r u s t a l attenuation and the assumed form of the incident wavelet w i l l be shown. Although exact syntheses have not been obtained, s i g n i f i c a n t s i m i l a r i t y with experiment w i l l be evident. The primary aim of the thesis has not been to determine the Al b e r t a crust, but rather to develop the new l i n e a r systems technique and demonstrate i t s p a r t i c u l a r usefulness. Thus, the nature of each of the' following experimental analyses i s b a s i c a l l y one of example. 3.2 The Analysis Problem The character of the P coda of a seismogram depends upon: a) the source generated waveform, b) the p h y s i c a l properties of the source to receiver paths such as Q, mantle l a y e r i n g , and the near re c e i v e r c r u s t a l structure and c) the response c h a r a c t e r i s t i c s of the 39 recording instruments. This l a s t e f f e c t can be w e l l determined by the c a l i b r a t i o n of the seismograph system. The combination of the source waveforms and the mantle path e f f e c t s determines the waveform incident on the base of the crust which i s necessary for.the s o l u t i o n techniques presented i n t h i s work. The assumed l i n e a r nature of e l a s t i c propagation i n the crust permits inverse approaches. For example, i f a seismogram i s known and the c r u s t a l response i s l i n e a r and known, the corresponding input wave- . form can be d i r e c t l y computed. Also, from a known input s i g n a l and i t s r e s u l t i n g seismogram, the unique response of the crust to e l a s t i c waves can be obtained. However, since there i s seldom any comparison information on the nature of the source or incident t e l e s e i s m i c wave-, forms, these inverse approaches are not u s e f u l to demonstrate the u t i l i t y of the present theory without further assumptions. Rather, frequency and time domain propagation solutions w i l l be compared with analyses of the experimental seismograms. 3.3 The Experiment ' For analyses a chosen set of seismograms recorded south of Edmonton, Albe r t a during a f i e l d program of October and November 1968 . has been used. Three portable FM tape recording short period seismo-graph stat i o n s at the U n i v e r s i t y of Alberta Edmonton Seismological Obser-vatory (our seismograph i s designated LED i n t h i s t h e s i s ) , and at the nearby Chamulka (CHA) and Larsen (LAR) farms (Figure 13) recorded 16 us e f u l events from teleseismic distances. Six of these were chosen f o r a n a l y t i c a l comparisons with the theory (see l a t e r Figure 16, Table III). This region of A l b e r t a i s characterized by a f l a t l y bedded sedimentary 40 F i g . 13. A map of the seismic s t a t i o n and o i l w e l l l o c a t i o n s by which the c r u s t a l s e c t i o n was determined. The w e l l s i t e s are designated by l e t t e r s and correspond to those of Figure 14. The Dominion Land Survey Range and Township coordinates are shown. 41 s e c t i o n which can i n general be approximated by the h o r i z o n t a l l a y e r model. Maximum dip angles on the notably continuous formation boundaries are generally l e s s than 1° except i n anomalous areas such as the Leduc-Woodbend o i l f i e l d (which l i e s beneath s t a t i o n CHA) i n which o i l - b e a r i n g Devonian reef structures e x i s t . Therefore, these records are i d e a l l y s u i t e d to t h i s problem. Also, the sedimentary s e c t i o n i s w e l l known from r e f l e c t i o n seismic surveys and o i l w e l l borehole geophysics f o r s e v e r a l nearby d r i l l s i t e s . 3.4 The Crust Model The s t a t i o n l o c a t i o n s formed a t r i a n g u l a r pattern 15 kilometers on a s i d e . Over these r e l a t i v e l y short distances, v e l o c i t y and r e s i s t i -v i t y logs from deep o i l w e l l s i n the v i c i n i t y permitted easy and c e r t a i n c o r r e l a t i o n s of the major sedimentary formations to determine an i n t e r -polated p r o f i l e s e c t i o n f o r the upper crust (Figure 14, Table I ) . A r e g i o n a l downward dip normal to the Rocky Mountain s t r i k e i n a d i r e c t i o n N 130°W was assumed i n order to f a c i l i t a t e the i n t e r p o l a t i o n between the w e l l s i t e s to the s t a t i o n l o c a t i o n s . (e.g. E l l i s and Basham, 1968). In order to a f f e c t the character of a seismogram, a formation l a y e r must have a thickness comparable to the wavelength of the highest frequencies recorded by the seismograph. More exa c t l y , a l i n e a r systems s o l u t i o n shows that the fundamental reverberation frequency of a l a y e r i s determined by the imaginary component of the Laplace frequency s at the p r i n c i p a l pole of the t r a n s f e r function set (Appendix I ) . For a non-attenuating l a y e r under normal incidence, t h i s condition corresponds to"one of two frequencies depending upon the r e f l e c t i o n c o e f f i c i e n t s at the l a y e r boundaries: CRUST SECTION LAYER A PSURFACE -I KM 10 KM -2 KM -3 KM Fig 14. The sedimentary crustal section obtained by interpolation between oilwell geophysical logs. The layer numbers correspond to those in Table I. The dotted boundaries underlying station CHA show the effect of the Leduc-Woodbend o i l f i e l d . The control oilwell sites used are designated by letters which correspond to those of Figure 13. Extrapolations are indicated by dotted lines. 43 either f = ^ °i 2 T i or f (3-1) 1 °i 4 T i [4] where T. is the layer P-wave (or S-wave) transit time. Thus, f i o^ determines the minimum frequency component for which the layer rever-berates; frequencies significantly lower than this are unseen. The sedimentary crust models (Figure 14, Table I) include a l l distinct layers for which the transit time, T^, is greater than 0.05 seconds. Consequently, the layer modelling is sufficiently detailed to include the fundamental layer reverberation character for frequencies markedly less than 10 Hz. Furthermore, a l l other layer boundaries which showed -notable velocity contrasts were also considered i n the crustal models, . thus allowing layers for which T^ approximates only 0.025 seconds ' -( f Q = 20Hz). The frequencies recorded by the seismographs are limited to a band below 5 Hz. The crust model detail should be sufficient to describe reverberations in this band. The subsedimentary crustal section included i n the total crustal model of Table I follows the gross crust of Cumming and Kanasewich (1966) for southern Alberta. The total crust model i s desirable since the Mohorovicic discontinuity f i r s t causes a part of an incident P wave to be converted into SV motions and particularly affects much of the early character of the horizontal seismograms. [4] The f i r s t condition results i f the product of the two reflection coefficients i s positive; the second i f the product i s negative. Also for the positive product, there is a pole on the u) = 0 axis (i.e. f = 0). This pole determines only the zero frequency character of the i layer response and does not affect wave propagation (Appendix I ) . 1 TABLE I CRUSTAL SECTIONS LAYER FORMATION (TOP) THICKNESS (km) P-VEL0CITY (km/sec) DENSITY (gin/ cc) LED CHA LAR LED CHA LAR LED CHA LAR 1 Post Alberta 0.62 0.58 0.71 2.7 2.7 2.7 2.2 2.2 2.2 2 Alberta 0.46 0.49 0.49 3.0 . 3.0 3.1 2.3 2.3 2.3 3 Blalrmore 0.31 0.28 0.35 3.7 4.1 4.1 2.4 2.4 2.4 4 Wabumum 0.16 0.23 0.21 5.5 5.6 5.6 2.6 2.6 2.6 5 Ireton 0.17 0.22 0.19 4.0 4.2 4.5 2.4 2.4 2.5 6 Cooking Lake 0.26 0.27 0.30 5.3 5.4 5.5 2.6 2.6 2.6 7 Elk Point 0.19 0.19 0.18 4.6 4.7 4.7 2.5 2.5 2.5 8 Cambrian 0.42 0.40 0.37 4.1 4.2 4.2 2.4 2.4 2.4 9 Precambrian 10.0 10.0 10.0 6.1 6.1 6.1 2.7 2.7 2.7 10 Sub-Layer I 24.0 24.0 24.0 6.5 6.5 6.5 2.8 2.8 2.8 11 Sub-Layer II 11.0 11.0 11.0 7.2 7.2 7.2 2.9 2.9 2.9 12 Mantle Boundary - - - 8.2 8.2 8.2 3.0 3.0 3.0 45 These model considerations do not exactly represent non-normal incidence as the fundamental resonance frequencies are s l i g h t l y increased as the i n c i d e n t angle increases. But since the waveform in c i d e n t angles are generally small for the teleseismic source distances, only small e f f e c t s are to be expected. Independent data on the layer S v e l o c i t i e s , formation d e n s i t i e s and e l a s t i c wave attenuation properties are not a v a i l a b l e f o r the A l b e r t a c r u s t . However, reasonably probable values f o r the S-wave v e l o c i t i e s were determined assuming Poisson's r a t i o i s 1/4 so that: B i = a ±//3 (3-2) Formation d e n s i t i e s were obtained from the P-wave v e l o c i t y - d e n s i t y r e l a t i o n of Nafe and Drake (Grant and West, 1965). Few r e s u l t s have been reported f o r the attenuation properties of body waves w i t h i n the crust and among these l i t t l e consistency i s apparent. The studies of McDonal et a l (1958) on the thick shale formations near P i e r r e , Colorado are among the most widely quoted. How-ever, other i n s i t u Q measurements by subsequent workers (e.g., T u l l o s and Reid, 1969) have not shown attenuation c o e f f i c i e n t s which c l o s e l y agree with t h e i r r e s u l t s . Various authors (Press, 1964; Anderson et a l , 1965) have i n f e r r e d body wave Q's from surface wave analyses. But such studies do not provide the r e s o l u t i o n required i n the present a n a l y s i s . Clowes and Kanasewich (1970) experimentally obtained Q f o r the c e n t r a l A l b e r t a crust using deep r e f l e c t i o n seismograms. They determined average Q near 300 i n the sediments and approximately 1500 i n the lower c r u s t . T h e i r sediment Q appears to be s i g n i f i c a n t l y greater than that found by 46 McDonal et a l for similar materials. In fact, order of magnitude differences appear to be common among the several Q's reported for sedimentary crusts. For this thesis, material Q values have been obtained from the literature (Knopoff, 1964; McDonal et a l , 1958; Savage, 1965) and correlated with the known formation materials in the sedimentary portion of the Alberta crust. Within the lower crust, Q's have been estimated on the basis of the average crustal attenuation: (Press, 1964). Com-prehensive data on the composition of the Alberta sediments has been compiled by the Alberta Society of Petroleum Geologists (1964). Since most of the major formations are substantially homogeneous, the correlations between the formation rock material and published Q factors for the materials should give plausible values. Q has been reported to be approximately frequency independent over a wide range of frequencies between 1 Hz and 10 6 Hz (Knopoff and MacDonald, 1958 and 1960). Thus, Q factors which have been published for relatively high frequencies have been assumed to be constant to the Futterman cutoff angular frequency U Q . Furthermore, at frequencies below 1 Hz, r e a l i s t i c Q values have almost no effect on body wave propagation i n the crust and the assumption that Q i s constant with frequency to the cutoff presents no apparent d i f f i c u l t i e s . Shear and compressional Q's have been assumed to be identical. The crustal model Q factors for the individual layers which have been chosen for this work are tabulated and referenced i n Table II. Of these, only those from the Alberta and post Alberta formations can be expected to be accurate. These formations are basically uniform shales which appear to be contiguous with the Pierre, Colorado shales for which 4 7 TABLE II FORMATION Q's FORMATION MATERIAL CHOSEN. Q ' APPLICABLE REFERENCE Post A l b e r t a . sandy shale 2 3 2 3 (Pierre shale) Alberta shale 3 0 2 1 - 5 2 : 1 . Blairmore shale, sandstone 3 0 2 1 - 5 2 . 1 Wabumum .dolomite 1 0 0 V 4 5 - 1 9 0 ? (limestone) 1 ' ,: Ireton shale 3 0 2 1 - 5 2 • ; 1 Cooking Lake limestone 1 0 0 ' 4 5 - 1 9 0 ? Elk Point s a l t , evaporites 1 0 0 7 0 - 2 2 0 (hard chalk) ' . 3;"' Cambrian sandstones, . e l a s t i c s 5 0 5 2 (sandstone) • I PreCambrian granite? 2 0 0 5 7 - 2 0 0 I Sub-Layer I 5 0 0 estimated, 4 Sub-Layer II — . • 5 0 0 - estimated, 4 Q (average f o r whole crust) = 1 8 6 : Q (average f o r sediments) = 4 0 . The average Q1 's were computed considering the t r a v e l times through the i n d i v i d u a l layers, References: 1 ) Knopoff ( 1 9 6 4 ) 2 ) McDonal et a l ( 1 9 5 8 ) 3 ) Savage ( 1 9 6 5 ) 4 ) Press ( 1 9 6 4 ) The v a l i d i t y of the Q's f o r the sub-sedimentary layers i s questionable. The r e s u l t s of Press ( 1 9 6 4 ) apply p a r t i c u l a r l y to Pg waves i n the t e c t o n i c a l l y a c t i v e region of the Nevada Test S i t e . For the remaining chosen Q's, values obtained i n the lab must be assumed to apply i n s i t u . 48 attenuation properties are w e l l known (Williams and Burk, 1964). McDonal et a l (1958) obtained Q=23 f o r P waves and Q=10 f o r S waves with frequencies beyond 20 Hz i n t h i s formation. In order to conform to the Q assignments t o other layers the Q=23 value has been assumed to describe the S wave attenuation i n t h i s l a y e r a l s o . The average Q f o r the sedimentary portion of the crust above the Precambrian by the above model i s approximately 40 while the average Q f o r the whole crust i s near 200. 3.5 Data Preparation The FM tape recorded seismograms were d i g i t i z e d with a sampling rate of 20 per second for the numerical analyses. The r a p i d f a l l o f f of the seismograph s e n s i t i v i t y from 5 Hz ensures that there i s l i t t l e seismo-gram content beyond the 10 Hz Nyquist frequency so that a l i a s i n g e f f e c t s are minimal. Transient c a l i b r a t i o n pulses were recorded at the beginning of each 5-day tape record. From these pulses, the absolute ground motion amplitude frequency response of each s t a t i o n component seismograph was determined using the transient c a l i b r a t i o n procedure o u t l i n e d by Deas (1969). However, the high frequency instrumental response character i s poorly determined by t h i s method since the d i g i t i z a t i o n (quantization) noise swamps the high frequency components i n the c a l i b r a t i o n pulses (Figure 15). Consequently, the t h e o r e t i c a l response f o r frequencies beyond 2 Hz was assumed to be representative of the systems. The response c h a r a c t e r i s t i c s of the three components of any one s t a t i o n were designed t o be v i r t u a l l y i d e n t i c a l ; only the absolute s e n s i t i v i t i e s s i g n i f i c a n t l y v a r i e d between components. The r e s u l t s of the pulse c a l i b r a t i o n s were used to e s t a b l i s h absolute s e n s i t i v i t i e s . 4 9 >-0.1 6.2 05 [5 £O 5'.0 16. ' FREQUENCY (HZ) F i g . 15. The t h e o r e t i c a l standard seismograph response (dotted l i n e ) compared to a response c a l i b r a t i o n obtained by the transient method ( s o l i d l i n e ) . The nature of the propagation problem i s e s s e n t i a l l y two dimensional; motions only occur v e r t i c a l l y and h o r i z o n t a l l y along the wave ray paths. For the comparisons of the theory with the seismograms, i t i s then necessary to perform a transformation on the N and E h o r i z o n t a l components to e s t a b l i s h the equivalent r a d i a l component motion along the surface p r o j e c t i o n of the ray path and the transverse component normal to 50 I t . For transformation, the two h o r i z o n t a l seismogram components must have i d e n t i c a l responses. For t h i s purpose, a l l s t a t i o n component s e n s i -t i v i t i e s were normalized to the response c a l i b r a t i o n of the v e r t i c a l component of s t a t i o n LED on October 4, 1968. The absolute s e n s i t i v i t i e s o f the instruments v a r i e d with time. E f f e c t s of temperature on the t r a n s i s t o r i z e d a m p l i f i e r s and seismometers p a r t i a l l y account f o r t h i s e f f e c t . 3.6 The Seismograms The s i x events used i n these comparison analyses (Table I I I , Figures 16a-f) were chosen p r i m a r i l y for t h e i r short and apparently simple P coda. For these events, the wave motion i n c i d e n t on the base of the crust should comprise a b a s i c a l l y t r a n s i e n t P wave of short duration. I f the source to r e c e i v e r path i s not str u c t u r e d i n i t s e l a s t i c p r o p e r t i e s , P to S wave conversions do not occur along the path. However, i t should be noted that d i s t i n c t s t r u c t u r e i n the upper mantle and continuous v e l o c i t y and density v a r i a t i o n s (hence acoustic impedance v a r i a t i o n s ) along the source to r e c e i v e r path must cause some conversion o f energy from compressional to r o t a t i o n a l motions. For these analyses, source generated P and.S motions are assumed to a r r i v e at the r e c e i v e r s e v e r a l minutes apart with the incident waveform e n t i r e l y compressional. The P coda i s thus considered to be the r e s u l t of the purely compressional waveform impinging on the base of the crust at the Mohorovicic d i s c o n t i n u i t y . E l l i s and Basham (1968) have suggested that an increase i n the transverse component amplitude with time ( t h i s increase i s observed f o r 51 LED LAR EVENT I VENEZUELA R 5 SEC. F i g . 16a-f. The v e r t i c a l , r a d i a l h o r i z o n t a l and transverse h o r i z o n t a l components of the experimental seismograms for the events studied i n t h i s work. The P-coda onset f o r each event i s exactly 5 seconds from the beginning of each record. (a) Venezuela EVENT 2 FIJI ISLANDS T 5 S E C . (b) F i j i Islands 53 LED EVENT 3 CHILE I C H A L A R 5 SEC. ( c ) Chile 54 EVENT 4 KURILE IS. LED LAR 5 SEC. (d) K u r i l e Islands 55 EVENT 5 NOVAYA ZEMLYA LED CHA R LAR 5 SEC. (e) Novaya Zemlya 56 EVENT 6 RYUKYU IS. . LED LAR 5 SEC. (f) Ryukyu Islands TABLE III EVENTS USED IN ANALYSES EVENT LOCATION LAT LONG DEPTH MAG DELTA INCIDENCE AZIMUTH (S-E) 1 Venezuela 10.7°N 62.6°W 97 km 4.8 58.6° 30° 117° E of N 2 F i j i Is. 20.9°S 178.8°W 607 km 5.7 92.9° 19° 238° 3 Chile 19.6°S 68.9°W 107 km 5.3 82.1° 21° 138° 4 Kurile Is. 49.7°N 155.8°E 35 km 5.5 53.0° 31° 305° 5 Novaya Zemlya 73.4°N 54.9°E 0 km 6.0 53.4° 31° 4° 6 Ryukyu Is. 27.5°N 128.4°E 48 km 5.8 83.4° 21° 308° 58 example on the Venezuelan event recorded at s t a t i o n LAR i n Figure 16a) may be due to s i g n a l generation of noise by s c a t t e r i n g w i t h i n the c r u s t a l l a y e r i n g or topographic features. This e f f e c t may then i n d i c a t e the rate of degradation of the record f o r the purposes of t h i s analysis-. For most of the chosen events, the r a d i a l and transverse motion ampli-tudes become comparable wi t h i n 15 to 30 seconds of the P onset. Twenty-f i v e seconds of the record has been used i n the subsequent frequency domain analyses of the seismograms. Generally,- one would not expect the noise generated transverse amplitudes to exceed the r a d i a l amplitudes f o r h o r i z o n t a l l y layered s t r u c t u r e s . However, some records such as the CHA seismograms from the Chilean earthquake (Figure 16c, Table III) show a predominance of h o r i z o n t a l transverse motion. I t should be noted that the s t a t i o n CHA was located on the f l a n k of the Leduc-Woodbend o i l f i e l d , a region i n which there e x i s t s s i g n i f i c a n t departure from h o r i z o n t a l l a y e r i n g i n the sedimentary crust (see Figure 14). Indeed, i t appears that the o i l bearing Devonian reef s t r u c t u r e causes the o v e r l y i n g s e d i -mentation l a y e r s to be l i f t e d . The complicated subsurface geometry at the s t a t i o n CHA does not permit easy p r e d i c t i o n of the departures from the expected zero transverse motions. Also, strongly topographic features of the nearby North Saskatchewan r i v e r v a l l e y ( t y p i c a l l y 1 kilometer across and 100 meters deep) could a f f e c t - t h e CHA records. Furthermore, although the CHA record of the.Chilean event,is notably anomalous with respect to transverse motions, the Novaya Zemlya records are not. This may i n d i c a t e that there i s an azimuthal dependence of anomalous transverse motions. The e s s e n t i a l l y anomalous nature of the CHA l o c a t i o n was known before analyses and the two events were included 59 only f o r the purpose of comparing the d i s t i n c t l y h o r i z o n t a l l a y e r i n g at LED and LAR with the more structured geometry of CHA. 3.7 Frequency Analyses - T h e o r e t i c a l S p e c t r a l Amplitudes For the determination of t h e o r e t i c a l frequency amplitude response c h a r a c t e r i s t i c s of the model cr u s t s , impulse (white, zero phase spectrum) s i g n a l s representing P waveforms were used as inputs to the l i n e a r system analogs. Amplitude and phase spectra were c a l c u l a t e d f o r the surface motions, f o r a l l i n t e r n a l boundary motions and f o r the wave-, forms retransmitted i n t o the basement halfspace. For non-normal incidence conditions, both the v e r t i c a l and h o r i z o n t a l motion amplitude components have been obtained. A frequency sampling i n t e r v a l of 0.1 Hz between 0 and 5 Hz was chosen to reveal the major s p e c t r a l peaks i n t h i s range while l i m i t i n g the computational time. The surface motion amplitude response f o r the sedimentary p o r t i o n of the Leduc model crust (Table I I I ) above the Precambrian for normal incidence i s shown i n Figure 17. The phase spectra showed mainly a lagging trend with frequency due to the time delay between the input waveform and the surface motions. The phase spectra and i n t e r n a l boundary motion spectra w i l l not be presented i n t h i s a n a l y s i s . I t i s notable that reverberation resonances occur- at- se v e r a l frequencies (Figure 17). The f i r s t - resonance peak of the LED c r u s t a l response at approximately 0.4 Hz i s p r i m a r i l y due to the reverberation m u l t i p l e with the longest path bounded by the surface and the Precambrian basement. (Appendix I describes a procedure f o r the computation of the fundamental reverberation frequency and harmonics f o r one l a y e r ) . • 60 F i g . 17. The s p e c t r a l amplitude response of the t h e o r e t i c a l non-attenuating Leduc c r u s t a l s e c t i o n to a normal incidence white amplitude spectrum P-wave incident at the Precambrian-Cambrian boundary. S p e c t r a l peaks f o r higher frequencies are the r e s u l t of the combined e f f e c t s of higher harmonics of t h i s reverberation arid primary resonances f o r multiples with shorter paths (which cause higher frequency reverber-ation) along with t h e i r coincident harmonics. The s p e c t r a l amplitude response f o r a 30° (at the Precambrian boundary) in c i d e n t P wave consists of both v e r t i c a l and h o r i z o n t a l components (Figure 18). There i s a marked s i m i l a r i t y i n the form of the v e r t i c a l amplitude component with that for the normal incidence condition although i t s amplitudes are reduced. P a r t i c u l a r l y 61 F i g . 18. The s p e c t r a l amplitude response of the t h e o r e t i c a l non-attenuating Leduc c r u s t a l s e c t i o n to a 30° in c i d e n t P wave at the Precambrian-Cambrian boundary. WQ represents the v e r t i c a l component motion, Y Q the h o r i z o n t a l . the resonance of the fundamental reverberation and the peaks near 1.6 Hz, 2.3 Hz and 2.9 Hz are common to both responses. This i s expected because f o r small incidence angles, the P motions wi t h i n the l a y e r i n g dominate the v e r t i c a l component while the SV motions dominate the h o r i z o n t a l com-ponent. In t h i s example, the h o r i z o n t a l amplitudes are lower than the v e r t i c a l ; however, f o r S wave inputs or P wave input motions with much l a r g e r incident angles the opposite condition would occur. A 30° inc i d e n t angle at the Precambrian boundary corresponds to a 42° angle at the Mohorovicic d i s c o n t i n u i t y . This incidence angle corresponds to a source 62 at only 25 e p i c e n t r a l distance ( a f t e r Ichikawa, see Bashaun, 1967). For such nearby sources, signals from r e f r a c t i v e and r e f l e c t i v e paths • contaminate the early P. coda. Also, geometrical spreading of the . wavefront i n v a l i d a t e s the necessary plane wave assumption. Experimentally, source distances greater than approximately 35° should be used. ' . ; . The e f f e c t of c r u s t a l absorption on the s p e c t r a l amplitudes f o r the 30° incidence condition i s shown i n Figure 19. For the c a l -c u l a t i o n of these example spectra Q=50 was chosen to represent every l a y e r f o r both P and S motions. This Q i s near the computed average for the sediments (Table I I ) . Acausal, zero phase attenuation with no attending dispersion i s considered. S i g n i f i c a n t l y , there i s a decrease O U J 00 Q f -10 Q o CC o O oJ LEDUC SECTION 30° INCIDENT P WAVE W0  Yrs —-\ 7 \ ' / \ / 1.0 2.0 FREQUENCY 3.0 (Hz) 4.0 5.0 F i g . 19. The s p e c t r a l amplitude response of the t h e o r e t i c a l c r u s t a l section to a 30° incidence P wave at the Precambrian-Cambrian boundary. The c r u s t ' i s assumed to have a Q of 50 i n every layer. '•'•. o f t h e s p e c t r a l a m p l i t u d e s w i t h i n c r e a s i n g f r e q u e n c y f o r b o t h v e r t i c a l and h o r i z o n t a l components. The c o n s t a n t Q model g i v e s a n e g a t i v e e x p o n e n t i a l dependence o f a m p l i t u d e s w i t h f r e q u e n c y which d e s c r i b e s t h i s c o n d i t i o n . A f u r t h e r s i g n i f i c a n t e f f e c t o f t h e a b s o r p t i o n i s the a p p a r e n t r e d u c t i o n i n t h e q u a l i t y o f t h e r e s o n a n c e s , p a r t i c u l a r l y a t h i g h e r f r e q u e n c i e s . The resonance q u a l i t y i s d e t e r m i n e d by the energy a b s o r b e d p e r c y c l e a t the r e s o n a n t f r e q u e n c y harmonic and i s a d i r e c t measure o f t h e c r u s t a l a t t e n u a t i o n . I t was n o t e d p r e v i o u s l y t h a t p h y s i c a l a t t e n u a t i o n must be c a u s a l which demands a d i s p e r s i o n e f f e c t s u c h as i n t h e Futterman model. The d i s p e r s i o n o f phase v e l o c i t i e s would cause a s l i g h t s h i f t o f t h e r e s o n a n c e peaks b u t the form o f the c u r v e a m p l i t u d e dependence on f r e q u e n c y would remain n e a r l y unchanged. S i n c e t h e e f f e c t s o f d i s p e r s i o n on t h e s p e c t r a l a m p l i t u d e s a r e m i n i m a l , t h e a c a u s a l model r e a s o n a b l y r e p r e s e n t s the c r u s t a l a m p l i t u d e r e s p o n s e o f an a t t e n u a t i n g c r u s t . 3.8 F r e q u e n c y A n a l y s i s - S p e c t r a l R a t i o T e c h n i q u e To compare the s p e c t r a l a m p l i t u d e r e s p o n s e s w i t h ground d i s -placement a m p l i t u d e s due t o e a r t h q u a k e s , i t i s n e c e s s a r y t o remove the e f f e c t s o f t h e seismograph r e s p o n s e . T h i s r e q u i r e s .extreme enhancement o f low f r e q u e n c y components where m i c r o s e i s m i c and i n s t r u m e n t n o i s e may exceed s i g n a l l e v e l s and h i g h f r e q u e n c i e s where wind, water and c u l t u r a l n o i s e may swamp motions g e n e r a t e d by the t e l e s e i s m i c waves. A l s o , t h e s e seismograms a r e b a s i c a l l y a r e c o r d o f t h e ground motion v e l o c i t i e s r a t h e r t h a n d i s p l a c e m e n t s and an i n t e g r a t i o n o f the r e c o r d s w i t h r e s p e c t t o time would be n e c e s s a r y b e f o r e s p e c t r a l a m p l i t u d e comparisons c o u l d be made. R a t h e r than a p p l y such r e c o n s t r u c t i o n s to the d a t a , the t e c h n i q u e 64 using s p e c t r a l r a t i o s (Phinney, 1964) i s more u s e f u l . Since the frequency response character of the seismographs were normalized, a l l motion com-ponents can be assumed to have been i d e n t i c a l l y a f f e c t e d by the instruments, For a given P a r r i v a l , the v e r t i c a l and r a d i a l h o r i z o n t a l s e i s -mogram P coda are responses to the same input waveform which i s i n c i d e n t at the base of the c r u s t . Their F o u r i e r spectra may be described: Wo(U) = H v ( c ) . U I N p u T ( c o ) (3-3) VU) -V U )' UINPDT ( U ) * where H^Cw) and H^(OJ) are the t o t a l crust v e r t i c a l and r a d i a l h o r i z o n t a l t r a n s f e r functions and u'j.NPUT^1^ i-s t n e F o u r i e r spectrum of the a r b i t r a r y input waveform. I t i s evident that the v e r t i c a l - t o - h o r i z o n t a l s p e c t r a l amplitude r a t i o depends only upon the c r u s t a l t r a n s f e r function r a t i o and i s completely independent of the input waveform (Phinney, 1964): W (co) | H>) V / H = X (3-4) |Yo(U)| I lyoo | Furthermore, over small h o r i z o n t a l distances, one may expect that the low frequency content of a P a r r i v a l i n c i d e n t on the base of the crust due to a t e l e s e i s m i c source would not vary s i g n i f i c a n t l y . Therefore, the r a t i o s between the v e r t i c a l component s p e c t r a l amplitudes at two nearby s i t e s should describe another r a t i o which i s independent of the input waveform. For example: V u ) 1 2 = V U ) ] 2 ' U I N P U T ( W ) (3-5) and: |Wo<o))]i| |H v(o))] i| V g / V j = - . (3-6) |w (co)] | |H (co)] | o 2 v 2 These v e r t i c a l - t o - v e r t i c a l component spectral amplitude r a t i o s may be useful i n the comparisons of the cru s t a l sections at two stations. Also, because the v e r t i c a l components usually possess high signal-to-noise r a t i o s r e l a t i v e to the horizontal components, the accuracy obtainable from such analyses should be better. Phinney showed that a comparison of the V/H ra t i o s determined experimentally from long period P-wave seismograms and by dir e c t c a l -culation using the Haskell method could be useful for c r u s t a l interpre-tations. Although, h i s work yielded good comparisons between theory and experiment, other investigations have been less successful at short periods (Utsu, 1966; E l l i s and Basham, 1968). The work of E l l i s and Basham.included an area i n central Alberta near that studied i n these analyses. They ascribed part of t h e i r l i m i t e d success to scattering which generates noise from s i g n a l . However, thei r Haskell matrix analyses did not include attenuation. The present work w i l l attempt to determine the significance of c r u s t a l absorption to spectral r a t i o s . Later, i t w i l l be shown that for short periods and the low Q's of c r u s t a l sediments, the effect of absorption can be considerable. F i g . 20a. Spec t r a l V/H r a t i o s f o r various attenuation conditions i n the sedimentary crust compared to the non-attenuating r a t i o s . The designated Q's have been applied to every l a y e r ; the s o l i d l i n e i s f o r i n f i n i t e Q. F i g . 20b. S p e c t r a l V/H r a t i o s f o r Q values i n each layer assigned according to Table II compared to non-attenuating r a t i o s . The s o l i d l i n e i s f o r i n f i n i t e Q; the dashed l i n e f o r the chosen Q's. 66 SPECTRAL V/H RATIOS 0 1.0 2.0 3.0 4.0 5.C FREQUENCY (HZ) 0 1.0 ZO 3.0 4.0 5.0 FREQUENCY (HZ) Fig. 20a. SPECTRAL V/H RATIOS 1.0 2.0 3.0 4.0 5.0 FREQUENCY (HZ) 0 1.0 2.0 3.0 4.0 5.0 FREQUENCY (HZ) CHA o i.o 2:0 3.0 4.0 5.0 FREQUENCY (HZ) Fig. 20b. 68. The work of E l l i s and Basham in central Alberta included the analysis of Al events which allowed them to use a s t a t i s t i c a l approach i n their experimental comparisons. In this present work only 6 events are used and s t a t i s t i c a l l y valid inferences cannot be deduced. 3.9 Theoretical Spectral Ratios Theoretical V/H ratios were computed for the three station crust models for an incident P wave making a 22° angle with the Precambrian upper boundary (which is equivalent to a 30° incidence angle at the Mohorovicic boundary). The computed ratios are presented i n Figures 20a and 20b. Several different internal layer attenuation conditions are compared. Since the broad features of the V/H ratios are not greatly affected by neglecting the large scale layering of the sub-sedimentary Alberta crust ( E l l i s and Basham, 1 9 6 8 ) , only the portion of the crust above the Precambrian has been considered in these analyses. The major effect of neglecting the three deeper model layers is to remove closely spaced harmonic peaks due to the long path multiples involving the deep crustal boundaries. Smoothing of the computed ratios was not considered necessary but the frequency sampling increment was chosen sufficiently small to determine a l l of the major spectral ratio features. Small changes in the P wave incidence angle have l i t t l e effect on the location of the V/H peaks although the peak amplitudes are modified. Consequently, i f the absolute ratio amplitudes are not required, the one set of ratios for 22° incidence should be sufficient for comparisons with the ratios experimentally determined from the suite of earthquake events with incidence angles ranging from 14° to 23°. The theoretical ratios for 22° incidence should underestimate the ratio amplitudes for the smaller 69 angles of P-wave incidence. The curves designated by i n f i n i t e Q (Figures 20a and 20b) describe non-attenuating c r u s t a l conditions; designations with other Q's represent the V/H r a t i o s f o r the crust model with the s p e c i f i e d Q applied to each l a y e r . The value Q=100 i s thought to be a conservative estimate of the absorption i n the A l b e r t a sediments but the value Q=25 must be implausibly low. The Q values of Table I I average 45 f o r the sediments. In these r a t i o comparisons, i t i s apparent that at low frequencies the e f f e c t of c r u s t a l absorption i s v i r t u a l l y i n s i g n i f i c a n t . But, above 1 Hz absorption e f f e c t s increase n o t i c e a b l y . For example, the LED and LAR V/H r a t i o peaks near 1.2 Hz are increased sharply by lowering Q while the r a t i o peaks near 1.7 Hz are decreased. Furthermore, the LED V/H peak near 2.4 Hz i s more than doubled by decreasing the t h e o r e t i c a l Q from i n f i n i t y to 25. For a l l three c r u s t a l models, large increases i n the r a t i o amplitudes are p a r t i c u l a r l y evident beyond 4 Hz f o r attenuating conditions. Attenuation of a harmonic wave increases exponentially with frequency and the e f f e c t of absorption must become a dominant f a c t o r i n V/H r a t i o s at high frequencies. The p o s i t i o n s of s e v e r a l of the major V/H peaks are not a f f e c t e d by the a d d i t i o n of attenuation. Therefore, these features should be useable f o r comparisons of the b a s i c c r u s t a l models with experiment. S i m i l a r l y , the features which- strongly depend on attenuation could pro-vide separate information on the c r u s t a l absorption c h a r a c t e r i s t i c s . T h e o r e t i c a l r a t i o s of the v e r t i c a l component s p e c t r a l amplitudes from s t a t i o n p a i r s are shown i n Figure 21. V e r t i c a l - t o - v e r t i c a l r a t i o s have been computed only f o r non-attenuating crust conditions because the absorption c h a r a c t e r i s t i c s at the three close l o c a t i o n s should be very 70 o o UJ r 6 O => t s! P < SPECTRAL VERTICAL-VERTICAL RATIOS CHA/LED — LAR/LED 1.0 2.0 FREQUENCY 3.0 (HZ) 4.0 50 F i g . 21. V e r t i c a l - t o - v e r t i c a l s p e c t r a l r a t i o s f o r non-attenuating crust models. s i m i l a r and the e f f e c t on the r a t i o s must be small. S i g n i f i c a n t l y , the t h e o r e t i c a l LAR/LED v e r t i c a l r a t i o s are fea t u r e l e s s below 3 Hz while the CHA/LED r a t i o s possess more character. This i s not s u r p r i s i n g since the r a t i o describes the r e l a t i v e v e r t i c a l c r u s t a l responses and since the LED and LAR crust models are quite s i m i l a r t h e i r response c h a r a c t e r i s t i c s should be s i m i l a r . The CHA crust model describes an anomalous l o c a t i o n on the Leduc-Woodbend o i l f i e l d . The d i f f e r e n c e i n the CHA and LED models i s manifested by the greater v e r t i c a l - t o - v e r t i c a l s p e c t r a l r a t i o character. 71 3.10 Experimental S p e c t r a l Ratios For comparisons of experimental and t h e o r e t i c a l s p e c t r a l r a t i o s , d i s c r e t e F o u r i e r spectra have been computed using the algorithm of Cooley and Tukey (1965) on the f i r s t 25 seconds (500 sample points) of the normalized records following the P onset. V e r t i c a l - t o - h o r i z o n t a l experimental s p e c t r a l r a t i o s computed from these spectra are shown for the chosen events i n Figures 22a-f. Within the 25 second time window the t r a n s i e n t P-coda has generally decayed to ground noise l e v e l s (e.g. Figure 16a-f) so that i t may be assumed that the spectrum, c a l c u l a t e d from the truncated records generally represents the i n f i n i t e duration t r a n s i e n t . The choice of a longer record would decrease the average s i g n a l - t o - n o i s e r a t i o and reduce confidence i n the r a t i o c a l c u l a t i o n s . Also l a t e r a r r i v i n g phases could be admitted. The 25 second record length i s a reasonable compromise f o r the chosen event records. The sampling i n t e r v a l of 0.05 seconds from the d i g i t i z a t i o n and the record length of 25 seconds gives a frequency r e s o l u t i o n of 0.04 Hz using the Cooley-Tukey algorithm. This r e s o l u t i o n i s excessive f o r comparisons with the t h e o r e t i c a l r a t i o computations. Smoothing of the experimental s p e c t r a using 5 and 9 point F e j e r ( t r i a n g u l a r ) running operators was applied to reduce the h i g h l y resolved d e t a i l while r e t a i n i n g the broad s p e c t r a l features which should r e s u l t from reverberations i n the,upper sedimentary c r u s t a l l a y e r s . For three chosen event seismogram s e t s , smoothed r a t i o s of the Fourier spectra of the s i g n a l and 25 seconds of pre-event noise f o r both v e r t i c a l and h o r i z o n t a l components were compared (Figure 23a-c) to determine confidence l i m i t s on the computed V/H r a t i o s . For these three EVENT I VENEZUELA SPECTRAL V/H RATIOS 0 1.0 2.0 3.0 4.0 5.C FREQUENCY (HZ) F i g . 22a-f. Experimental V/H spectral ratios obtained from the suite seismograms. Error bars are included on 3 of the spectral ratios for comparison. These error conditions are thought to be typical. (a) Venezuela 73 EVENT 2 FIJI ISLANDS SPECTRAL V/H RATIOS 8 LED °-\ o % SI < 1.0 2.0 FREQUENCY 3.0 (HZ) 4.0 5.0 2.0 FREQUENCY 5.0 (b) F i j i Islands EVENT 3 CHILE SPECTRAL V/H RATIOS 74 8 LED U J « > | o => I-j I -(c). Chile EVENT 4 KURILE IS. SPECTRAL V/H RATIOS (d) Kurile. Islands 76 (e ) Novaya Zemlya 77 LAR 1.0 2.0 3 0 4.0 5.0 FREQUENCY (HZ) (f) Ryukyu I s l a n d s 78 F i g . 23a-c. S p e c t r a l s i g n a l - t o - n o i s e r a t i o s f o r the v e r t i c a l ( s o l i d l i n e ) and h o r i z o n t a l (dotted l i n e ) f o r three events. The n o i s e spectrum was obtained from the 25 seconds preceding the P onset of each event. 7 9 analyses, large s i g n a l - t o - n o i s e r a t i o s were found between 1 and 3 Hz. The expected root mean squared errors due to the background noise are desig-nated by the erro r bars on the corresponding s p e c t r a l r a t i o s (Figures 22a,d,e). The LED Venezuela recording showed t y p i c a l r e l a t i v e pre-event noise; only the F i j i Islands records at LED and LAR possessed lower apparent s i g n a l - t o - n o i s e r a t i o s . Thus the r e l a t i v e errors found f o r the LED Venezuela r a t i o i s considered to be t y p i c a l of (or worse than) the remaining r a t i o s f o r which e r r o r conditions were not obtained. (One further p o s s i b l e exception i s the LAR Venezuela recording f o r which there i s a p a r t i c u l a r l y high frequency noise contamination. However, the signal - t o - n o i s e r a t i o at useable frequencies does not appear to be poor). 3.11 General Features of Spectral Ratios The most notable general feature of the experimental V/H r a t i o s i s that the amplitudes are much lower than the theory p r e d i c t s . (the t h e o r e t i c a l computations f o r the large 22° incidence angle should under-estimate r a t i o amplitudes f o r the steeper angles). The only p o s s i b l e exception i s the CHA C h i l e seismograms (Figure 16c, Table III) f o r which high r a t i o s are evident between 0 and 2 Hz but return to the anomalous lower values beyond t h i s . This event recording i s also anomalous with respect to h o r i z o n t a l motions since the transverse amplitudes exceeded r a d i a l amplitudes. Reduced r a t i o amplitudes must r e s u l t from a r e l a t i v e enhance-ment of the h o r i z o n t a l ground motions compared to the v e r t i c a l motions. P o s s i b l e causes include background noise and s i g n a l generated noise by s c a t t e r i n g . I t i s s i g n i f i c a n t that the h o r i z o n t a l component background noise was generally higher than the v e r t i c a l . This e f f e c t would be 80 expected to reduce the r a t i o amplitudes towards values l e s s than 1 at frequencies where the s i g n a l - t o - n o i s e r a t i o i s small. However, records f o r which the average s i g n a l - t o - n o i s e r a t i o appears to be very high would then be expected to show s u b s t a n t i a l l y b e t t e r agreement with the theory. But even those with the l e a s t r e l a t i v e background noise also e x h i b i t reduced r a t i o amplitudes. For example, the LED seismograms from the C h i l e and Novaya Zemlya events show anomalously low r a t i o amplitudes. The'computed r a t i o s f o r which e r r o r conditions were c a l -culated (Figures 22a,d,e) determine that these reduced amplitudes are not due to background noise, at l e a s t i n the band from 1 to 3 Hz. In t h e i r work on V/H r a t i o comparisons, E l l i s and Basham (1968) suggested that s c a t t e r i n g noise generated from the s i g n a l could have caused part of t h e i r l a c k of comparison between t h e i r t h e o r e t i c a l and experimental r e s u l t s . S c a t t e r i n g can r e s u l t from both i n t e r n a l c r u s t a l i r r e g u l a r i t i e s and from surface topography which determines a n o n - f l a t f r e e surface boundary. To explain a reduction of r a t i o amplitudes by s c a t t e r i n g , i t i s necessary that the mechanism p r e f e r e n t i a l l y reduces v e r t i c a l component motions or p r e f e r e n t i a l l y converts v e r t i c a l motions to h o r i z o n t a l . Since a v e r t i c a l to h o r i z o n t a l conversion i s b a s i c a l l y a P to S conversion, a s c a t t e r i n g mechanism which increased shear motions r e l a t i v e to compressional motions could explain the reduced apparent experimental r a t i o s . The e f f e c t s of topographic s c a t t e r i n g should most a f f e c t the LED s t a t i o n which was located near a broad v a l l e y and the CHA s t a t i o n which was s i t u a t e d w i t h i n 2 km of the North Saskatchewan River v a l l e y which i s approximately a h a l f km wide and 100 m deep. S u r p r i s i n g l y , only the CHA s t a t i o n gave s p e c t r a l r a t i o s which were not anomalously low. 81 An attenuation condition which, p r e f e r e n t i a l l y absorbed P motions r e l a t i v e to S motions could also help to explain t h i s c o n d i t i o n . However, most experimental evidence (e.g. McDonal et a l , 1958) shows that S attenuation i s greater than P which would determine the inverse e f f e c t . 3.12 D e t a i l s of the S p e c t r a l Ratios The V/H'ratios from the recordings of the s i x events have a common s i g n i f i c a n t s p e c t r a l r a t i o maximum between 0.5 and 1 Hz. (The r a t i o s computed for the CHA recordings show l e s s character f o r t h i s peak and i t s v a l i d i t y i s questionable i n t h i s case). The t h e o r e t i c a l c a l c u l a t i o n s f i n d a major r a t i o maximum near 0.7 Hz which c o r r e l a t e s with the experimental feature. Most r a t i o s computed f o r the LAR and LED s t a t i o n events show a furth e r common maximum between 1.5 and 2 Hz. However, f o r the LED Novaya Zemlya r a t i o s , t h i s peak i s di s p l a c e d to approximately 1.3 Hz which i s suspicious although the record q u a l i t y and s i g n a l - t o - n o i s e r a t i o are e s p e c i a l l y good. T h e o r e t i c a l r a t i o s f o r LED and LAR show maxima at 1.6 and 1.7 Hz which apparently c o r r e l a t e with these experimental features. Low frequencies are p r i m a r i l y a f f e c t e d by the structures of the crust which have a s c a l e s i z e of the order of the harmonic wavelength. Since, the gross c r u s t a l features should be w e l l modelled, i t i s not s u r p r i s i n g that t h i s area of agree-ment between theory and experiment e x i s t s . However, beyond.approximately 2 Hz the l i t t l e c o r r e l a t i o n between the theory and experiment could i n d i c a t e that the small scale model structure inadequately describes the a c t u a l c r u s t s . Also, since P - a r r i v a l s i g n a l l e v e l s decrease r a p i d l y 82 with frequency, background noise begins to play an i n c r e a s i n g r o l e at higher frequencies. Beyond -3 Hz, the s i g n a l - t o - n o i s e r a t i o s become so low (Figure 23) that good c o r r e l a t i o n s cannot be expected. Where the s i g n a l - t o - n o i s e r a t i o i s near (or below) 1, background noise must l a r g e l y determine the V/H r a t i o s . But, even beyond 3 Hz the smoothed r a t i o s computed s o l e l y from the pre-event noise f o r the three events were found not to be s i m i l a r to the equivalent event r a t i o s . Con-sequently, i f the-background noise can be assumed to be a s t a t i o n a r y process, i t does not appear to be the whole cause of the high frequency discrepancies between theory and experiment. This lends weight to a s c a t t e r i n g noise hypothesis. The cause of the large LAR V/H r a t i o peak near 3.5 Hz i s not evident. However, s i g n a l - t o - n o i s e r a t i o s are very small (as evidenced by the large e r r o r bars (Figure 22d)) beyond 3 Hz and the peak i s probably not an i n d i c a t i o n of c r u s t a l features. Any s i m i l a r i t y i n the experimental and t h e o r e t i c a l r a t i o s f o r the CHA C h i l e event recordings i s s u r p r i s i n g because t h i s s t a t i o n was the most anomalous i n terms of subsurface l a y e r i n g and topography. Furthermore, the C h i l e event records showed p a r t i c u l a r l y large transverse h o r i z o n t a l motions. Apparent c o r r e l a t i o n s should be discounted as a c c i d e n t a l . The i n c l u s i o n of attenuation i n the t h e o r e t i c a l s o l u t i o n s did not improve the agreement between theory and experiment. Furthermore, the experimental d e t a i l s do not discriminate the various c r u s t a l Q models. For example, although the 1.7 Hz peak was noted as a p o s s i b l e i n d i c a t o r of c r u s t a l attenuation conditions i t i s apparent that there i s " no consistency i n i t s r e l a t i v e amplitude among the experimental r a t i o r e s u l t s . 83 For i d e a l h o r i z o n t a l l a y e r i n g , the V/H r a t i o i s independent of the azimuth of the P-wave ray approach. However, two p a i r s of e x p e r i -mental r a t i o s show strong azimuthal dependence. The LED r a t i o s from the Venezuela (azimuth 117° E of N) and C h i l e (azimuth 138° E of N) events give nearly i d e n t i c a l r a t i o character, within the expected e r r o r l i m i t s , throughout the strong s i g n a l band below 3 Hz. A s t r i k i n g s i m i l a r i t y i s not so noticeable on the corresponding LAR recordings. A l s o , the LED Ryukyu I s . (308° E of N) and K u r i l e I s . (305° E of N) r a t i o s are quite s i m i l a r to about 3 Hz; again, the LAR r a t i o s do not e x h i b i t much s i m i l a r i t y . For the LED recordings, such azimuthal c o r r e l a t i o n s are more notable than c o r r e l a t i o n s between events at s i m i l a r receiver-to-event e p i c e n t r a l distances. Admittedly t h i s sample o f events i s small, but the r e s u l t s seem to i n d i c a t e that the i d e a l h o r i z o n t a l l a y e r i n g condition does not e x p l i c i t l y apply to the crust beneath LED. Experimentally determined s p e c t r a l r a t i o s between s t a t i o n v e r t i c a l component pai r s are shown i n Figure 24a-c. The LAR/LED r a t i o s from the Venezuela, C h i l e , Novaya Zemlya and K u r i l e I s . events generally show the.expected r e l a t i v e l y f e a t u r e l e s s r a t i o s below 2.5 Hz. But, beyond 3 Hz the r a t i o s become much l a r g e r than the theory p r e d i c t s . The F i j i and Ryukyu r a t i o s also appear to be anomalous at the low frequencies. Both have more character and l a r g e r amplitudes than expected. L i t t l e character from the t h e o r e t i c a l r a t i o s c o r r e l a t e s with the experimental r a t i o s ; however, the K u r i l e I s . , Venezuela and Novaya Zemlya experimental r a t i o s a l l show s i m i l a r character below 3 Hz. Since the s i g n a l - t o - n o i s e r a t i o s are generally high at these frequencies, some of t h i s s i m i l a r i t y must be s i g n i f i c a n t although i t s i n t e r p r e t a t i o n SPECTRAL VERTICAL-VERTICAL RATIOS 84 FREQUENCY (HZ) i s not evident. The CHA/LED r a t i o s from the C h i l e and Novaya Zemlya events show l i t t l e c o r r e l a t i o n with theory. However, beyond about 3 Hz a r a t i o enhancement r e l a t i v e to theory i s also observed as f o r the LAR/LED r a t i o s . The enhancement suggests that the LED v e r t i c a l com-ponent la c k high frequencies r e l a t i v e to the other two s t a t i o n s . The -LED instruments were located on the Canadian seismic network p i e r on bedrock while the LAR and CHA instruments were les s i d e a l l y located on alluvium. High frequency noise i s e a s i l y generated by wind i n unconsoli-dated materials and probably contributes to the greater noise at s t a t i o n s LAR and CHA. 3 . 1 3 Time Domain Syntheses In p r i n c i p l e , i t i s p o s s i b l e to determine exactly one among the c r u s t a l transmission response, the i n c i d e n t waveform from the base-ment and the output seismogram set provided the other two are known e x p l i c i t l y . However, the nature of the propagation problem only allows d i r e c t knowledge of the received seismograms. In the following syntheses, the approach taken i s to obtain synthetic seismograms from the response of the system model of the crust to incident waveforms created on the b a s i s of simple and p l a u s i b l e assumptions. Then, com-parisons between the s y n t h e t i c and experimental seismograms are expected to demonstrate the usefulness of the l i n e a r systems approach and the v a l i d i t y ( i n part at l e a s t ) of the A l b e r t a crust models f o r propagation problems at low frequencies. I t i s arguable that the inverse approach which s t r i v e s to obtain the input waveform or the c r u s t a l response function from the experimental seismograms i s s c i e n t i f i c a l l y b e t t e r . However, i t i s evident that assumptions on both the nature of the waveform and the 8 6 c r u s t a l response function w i l l s t i l l be required f o r comparison. Thus, by the inverse approach, comparisons can be d i r e c t e d only between com-puted and assumed s o l u t i o n s . However, the creation of synthetic seismograms allows the d i r e c t comparison of solutions (which are dependent on the applied assumptions) with the r e a l i t y of the e x p e r i -mental seismograms. By simulating the continuous l i n e a r systems crust models using the methods discussed i n Chapter I I , time domain responses of the crus.t to any input displacement waveform can be c a l c u l a t e d . A FORTRAN simulation program f o r these syntheses i s l i s t e d i n Appendix IV. A requirement of t h i s d i g i t a l synthesis procedure i s that each l a y e r t r a n s i t time i s an integer m u l t i p l e of the sampling i n t e r v a l . However, i f the sampling i n t e r v a l i s chosen small enough to allow a Nyquist frequency beyond the frequencies of i n t e r e s t , the requirement i s not p a r t i c u l a r l y r e s t r i c t i v e . The experimental seismograms used i n t h i s work are band l i m i t e d to frequencies below about 5 Hz. The FORTRAN program of Appendix IV assigns a minimum layer t r a n s i t time of one time sample i n t e r v a l . For the f o l l o w i n g syntheses the time increment of 0.05 seconds has been used for comparisons with the d i g i t a l experimental data. V e r t i c a l component synthetic responses of the upper sedimentary p o r t i o n of the LED crust model to an example a r b i t r a r y incident P d i s -placement waveform at normal and 30° angl es are shown i n Figure 25. The reverberations of the primary motion l a s t i n g s e v e r a l seconds due to long path m u l t i p l e s w i t h i n the modelled crust are p a r t i c u l a r l y notable. The small s c a l e high frequency character r e s u l t s from shorter path multiples 87 0.5 SEC. Fig. 25. Comparison of non-normal (a) and normal (b) incidence synthetic seismograms. Wavelet (c) incident on the Precambrian-Cambrian boundary at 30° produces synthetic seismogram (a), at normal incidence produces synthetic seismogram (b). within the thinnest layers. For the 30° incidence synthetic response, the l a t e reverberations arrive e a r l i e r than the corresponding reverber-ations for normal incidence. In frequency solutions t h i s same eff e c t appears as a displacement of the reverberation spectral peaks toward higher frequencies as the incidence angle increases. (This effect can also be noted i n Figures 17 and 18). For non-normal incidence, the 88 h o r i z o n t a l motion response i s also a v a i l a b l e . The method permits the P and S waveforms retransmitted into the basement to be found but these were not generally used i n t h i s work. 3.14 Syntheses with Layer Attenuation Causal minimum phase attenuation within each l a y e r of the c r u s t Is allowed i n the time domain syntheses by the i n c l u s i o n of the attenu-a t i o n time operators a^(t) (see Chapter II) where each operator i s the impulse response of the attenuating layer f or one t r a n s i t . In t h i s work, the Futterman (1960) absorption model (A1,D1) described e a r l i e r has been used to compute the time incremented attenuation operators which represent each l a y e r f or i t s p a r t i c u l a r Q and t r a n s i t time. An example.of one such attenuator impulse response of the post A l b e r t a c r u s t a l l a y e r formation with Q=23 i s presented i n Figure 26. I t i s p a r t i c u l a r l y s i g n i f i c a n t that i t s response r a p i d l y decays to values near zero by 0.1 seconds following onset. The t r a n s i t time through the post A l b e r t a formation at i t s P-wave phase v e l o c i t y (undispersed) i s approximately 0.23 seconds which i s considerably longer than t h i s e f f e c t i v e 0.1 second "length" of the operator. The e f f e c t i v e "length" increases with decreasing. Q. Therefore, other c r u s t a l l a y e r formations which have higher Q's determine attenuation operators which are even shorter r e l a t i v e to t h e i r l a y e r t r a n s i t t i m e s . ^ ^ The FORTRAN simulation program of Appendix IV demands that the operator i s shorter than the layer t r a n s i t time which l i m i t s l a y e r Q's to values greater than approximately 5. Although such Q's are already u n r e a l i s t i c a l l y low, program modifications could be made to allow the very long operators which correspond to a r b i t r a r i l y low Q's. 89 0.06-0.04-0.02-0.20 0.25 0 30 SECONDS F i g . 26. The impulse response of l a y e r 1 of the c r u s t a l s e c t i o n (Table I) assuming causal attenuation with Q=23. Assuming the Futterman attenuation d e s c r i p t i o n the non-attenuated impulse would a r r i v e at 0.23 seconds. The amplitude s c a l e i s a r b i t r a r y . The time incremented l a y e r attenuators must be sampled and properly normalized to preserve the r e l a t i v e s p e c t r a l response. For the 0.05 second sampling i n t e r v a l , the layer attenuator operator weights f o r the P and S motions of a l l model layers are given i n Table IV. For each l a y e r , the Q values represented by these attenuator operators are approximately equivalent f o r both P and S motions to those given i n Table I I . However, s e v e r a l of the thinnest model layers produce such small attenuation i n one t r a n s i t that no broadening or d i s p e r s i o n of an 90 TABLE IV ATTENUATOR WEIGHTS FOR THE 0.05 SECOND SAMPLING INCREMENT FORMATION Q P WEIGHTS .S WEIGHTS Post Alberta 23 0.839,0.139,0.022 1 0.621,0.320,0.044,0.015 Alberta 30 0.932,0.06 0.831,0.154,0.015 Blairmore 30 0.953,0.047 0.947,0.053 Wabumum 100 1.0 1.0 Ireton 30 1.0 1.0 Cooking Lake 100 1.0 1.0 Elk Point 100 1.0 1.0 Cambrian 50 0.979,0.021 0.974,0.026 Precambrian 200 0.946,0.054 0.940,0.060 Sub-Layer I 500 0.915,0.085 0.815,0.168,0.017 Sub-Layer I I 500 0.946,0.054 0.940,0.60 91 F i g . 27. The t h e o r e t i c a l frequency response of the crust f o r one passage of an impulse f o r average Q=200 compared to the frequency response due to the sum of the time domain attenuation operators of Table IV. The value Q=200 i s near the average model Q value of 186. impulse can be seen on the 0.05 second time s c a l e . For these l a y e r s , the s i n g l e weight value unity which corresponds to i n f i n i t e Q or zero absorption must be applied. The s p e c t r a l response of the layer attenu-ators for one t r a n s i t of an impulse through the t o t a l c r u s t a l s e c t i o n i s compared to the t h e o r e t i c a l s p e c t r a l amplitude response f o r an average 92 Q equal to 200 (Figure 27). At the low frequencies, i t i s evident that the computed l a y e r attenuation operators underestimate the average Q while they s l i g h t l y overestimate Q at high frequencies. A r b i t r a r i l y b e t t e r agreement can be achieved by using smaller time sampling i n t e r v a l s at the expense of increased computational time. However, the chosen attenuators are considered to be an adequate d e s c r i p t i o n of attenuation i n the following syntheses. Normal incidence s y n t h e t i c crust responses for the upper s e d i -mentary p o r t i o n of the crust model and f o r the waveform used i n the syntheses of Figure 25 show the e f f e c t of i n c r e a s i n g the c r u s t a l attenuation (Figure 28). The i n f i n i t e Q response i s compared to the crust response with the appropriate attenuators of Table IV and to another response f o r which Q was a r b i t r a r i l y chosen to be approximately 15 i n every l a y e r . The s i g n i f i c a n t d i s p e r s i o n e f f e c t s on the i n i t i a l waveform and the p r e f e r e n t i a l removal of the high frequency character are evident. The apparent energy and the e f f e c t i v e length of the record are also reduced. 3. 15 Impulse Synthetic C r u s t a l Responses The Dirac d e l t a function i s often a convenient input s i g n a l f o r many continuous l i n e a r systems problems. However, numerical simulations with i n f i n i t e amplitude zero length impulses are not p o s s i b l e . Rather, i t is u s e f u l to represent the impulse area numerically and compute the areas of the response impulses. Synthetic crust responses to a 30° incidence P impulse at the Mohorovicic d i s c o n t i n u i t y f o r the three c r u s t a l models are presented i n Figure 29. S i g n i f i c a n t s i m i l a r i t y i n both the v e r t i c a l and h o r i z o n t a l component syntheses for the three s t a t i o n crust models i s evident. 0.5 SEC. F i g . 28. Normal incidence synthetic seismograms dependence on Q. I n f i n i t e Q (a) i s compared to the model Q values of Table IV (b) and to an average c r u s t a l Q=15 ( c ) . Note the dispersion of the primary a r r i v a l . The i n c i d e n t waveform of Figure 25 has been used. A l a r g e r scale P impulse normal incidence synthetic response f o r the LED crust model i s shown i n Figure 30. Several of the rever-beration multiples which cause delayed impulses are designated by integer s e t s . For these, each boundary of the crust model i s numbered from 0 f o r the surface to 11 f o r the crust-mantle boundary. The f i r s t i nteger represents the f i r s t r e f l e c t i o n ; the next represents the CRUST IMPULSE RESPONSE 94 LED CHA H — LAR H 5 SEG. Fig. 29. The vertical and radial, horizontal synthetic impulse responses of the three crustal sections representing the three stations. TRANSMISSION IMPULSE RESPONSE LEGEND B C D E F G 1 SECOND I 1 A: DIRECT B: 5-7-5-0 C: 0-4-0 D: 0-7-0 E: 0-3-0-2--0 & F: 0-3-0-3--0 G: 0-3-0-8--0 & H: 2-3-2-0 I: 0-2-0 J: 0-3-0 K: 0-5-0 L: 0-8-0 M: 0-3-0-4--0 & N: 0-7-0-8-•0 & 0: 0-9-0 P: 0-10-0 Q: 0-11-0 0-2-0-3-0 0-8-0-3-0 0-4-0-3-0 0-8-0-7-0 Fig. 30. The normal incidence'synthetic impulse response for the LED crustal section. Various multiples or reverberations are designated by letters. The legend describes the internal reverberation paths which cause these later arrivals. BASEMENT REFLECTION IMPULSE RESPONSE LEGEND D E F G 7 T A: B: C: D: E: F: G: H: 11-BASEMENT 10-BASEMENT 9-BASEMENT 8-BASEMENT 3-BASEMENT 1-BASEMENT SURFACE BASEMENT 0-8-O-BASEMENT 1 SECOND I 1 Fig. 31. The basement reflected wave synthetic impulse response for the LED crustal section. Reflections from various internal crust boundaries into the basement are designated by letters. One reverberated multiple, H, is also shown. The legend describes the internal reverberation paths causing the later arrivals. 97 second and so on. The l a s t integer represents the boundary on which the motion i s observed. Primary multiples which f i r s t r e f l e c t from the surface and then from an i n t e r n a l boundary before returning to the surface u s u a l l y d i s p l a y the l a r g e s t amplitudes. However, some high amplitude impulses from e n t i r e l y i n t e r n a l multiples are a l s o evident. For example, the multiples designated (2-3-2-0) and (5-7-5-0) are not r e f l e c t e d from the surface. The p a r t i c u l a r l y strong (0-3-0) m u l t i p l e due to the Wabumum upper boundary and surface and the (0-8-0) due to the Precambrian upper boundary should a f f e c t the e a r l y character of the experimental seismograms. The impulse response displacement wave propa-gated back i n t o the mantle shows s e v e r a l r e f l e c t i o n s from the major boundaries of the crust (Figure 31). These waveforms r e f l e c t e d back i n t o the mantle can r e t a i n much character created i n c r u s t a l reverber-ations . 3.16 Incident Waveform Models The computation of any reasonable s y n t h e t i c seismogram requires a r e a l i s t i c input waveform f o r the l i n e a r systems model which i s based upon the known properties of the mantle path and the nature of seismic source motions. But, because the form of P motions at sources are not generally known, p l a u s i b l e assumptions must be applied. For the Novaya Zemlya event, which r e s u l t e d from a Soviet underground nuclear t e s t of a device i n the 0.1 to 1.0 megaton TNT equivalence range, the r e s u l t s of Werth and Herbst (1967) f o r explosion waveforms have been considered. Their work on small nuclear explosions (Werth and Herbst, 1967; Figure 1) shows that the r a d i a l displacement time function near an explosion source has the form of a p o s i t i v e pulse 98 of short duration (les s than 0.5 seconds f o r t h e i r three examples). Based on t h e i r evidence, an impulse displacement time function source model has been assumed for the Novaya Zemlya event. Werth and Herbst suggested that waveforms at t e l e s e i s m i c distances are not h i g h l y dependent upon the exact nature of the explosion source motion pulse so that the impulse approximation should be an adequate d e s c r i p t i o n . I t should be noted that the f i e l d seismograph instruments observed only a band of frequencies below about 5 Hz. Assuming" a l i n e a r nature f o r the wavemotions from the source through the mantle, only those source frequencies w i t h i n the instrument bandpass can a f f e c t the seismogram. Consequently, the "impulse model" merely requires that the source wave-forms look l i k e an impulse below 5 Hz. The exact i n f i n i t e impulse of i n f i n i t e s m a l duration i s not required. An impulse model i s also p a r t i c u l a r l y a t t r a c t i v e because i t Imposes the minimum p o s s i b l e character on the source displacement motions. I t determines no s p e c t r a l amplitude or phase character and describes only an o r i g i n time; i t i s the maximum entropy (minimum Information) waveform. The chosen events i n t h i s work were l i m i t e d to those e x h i b i t i n g very short P coda. Since c r u s t a l reverberations can only lengthen the P coda, the i n c i d e n t P wave must have even shorter duration. Among the events, the LED recording of the Venezuela event has character which i s s u r p r i s i n g l y s i m i l a r to the Novaya Zemlya seismogram. This should i n d i -cate that the waveforms i n c i d e n t on the crust were a l s o quite s i m i l a r . Furthermore, since the two events occurred at nearly equal e p i c e n t r a l distances from the s t a t i o n s , the e f f e c t s of the mantle path should also be s i m i l a r . This implies s i m i l a r i t y i n the source generated waveforms, 99 at l e a s t as seen by the narrow band instruments. However, the remaining earthquakes do not appear to generate such s t r i k i n g l y s i m i l a r seismograms and the s t r i c t impulse model cannot be used as c o n f i d e n t l y . Some step displacement source models were also considered and-syntheses using these compared to impulse source syntheses. I t seems reasonable that a step model could b e t t e r represent e l a s t i c motions near the source of an earthquake f a u l t . Following the lead of Werth and Herbst (1967) and Carpenter (1967), the e f f e c t of the mantle path on the source waveform i s considered to be due only to absorption and attending d i s p e r s i o n . Carpenter shows the Futterman absorption model obeys a p r i n c i p l e of s i m i l a r i t y which he describes: F(t,T/Q) = (Q/T)F Q(tQ/T) (3-7) where F(t,T/Q) i s the waveform amplitude at time t f o r an event-to-station t r a v e l time of T and a medium with q u a l i t y Q. Thus, F Q(tQ/T) para-m e t r i c a l l y describes only the shape of the waveform and T/Q i s a c h a r a c t e r i s t i c time which i s the e s s e n t i a l parameter of F . In Figure 32, Futterman model responses to impulse displacement source motions are com-pared f o r various mantle Q's ranging from 500 to 3000 f o r a source to s t a t i o n t r a v e l time of 600 seconds. For these responses, the c h a r a c t e r i s t i c time parameter ranges from 1.2 seconds down to 0.2 seconds, the 10 minute t r a v e l time corresponds roughly to the e p i c e n t r a l distance to the Novaya Zemlya and Venezuela event sources or approximately 55°. The combined source motions and.mantle absorption waveforms may be f u r t h e r combined with the seismograph response to determine the 100 00 A 596 598 600 602 604 SECONDS F i g . 32. The e f f e c t of various mantle path average Q values on an impulse source displacement motion which has t r a v e l l e d 10 minutes at the undis-persed phase v e l o c i t y through the mantle. Q=1500 has been found to be the best representation for construction of synthetic seismograms i n t h i s work. input wavelet to the systems model. Since the source motions, mantle absorption and instrument response c h a r a c t e r i s t i c s are a l l l i n e a r t h e i r order of combination does not a f f e c t the r e s u l t a n t seismogram syntheses. Wavelets obtained by the convolution of the seismograph impulse response and the response of the mantle absorption model to impulse and step source motions f o r various T/Q are presented i n Figure 33. S i m i l a r wavelets were obtained by Carpenter using somewhat d i f f e r e n t instrument response c h a r a c t e r i s t i c s . 10.1 IMPULSE SOURCE T/Q= 1.2 SEC 0.6 0 .4 0.2 STEP SOURCE 0;6 0;3 5 SECONDS F i g . 33. Wavelets constructed by convolution of the seismograph response with absorbing mantle responses f o r various T/Q parameters f o r both impulse and step sources. T/Q=0.4 seconds has been found to be the best representation f o r seismogram synthesis at an e p i c e n t r a l distance of 55°. IMPULSE SOURCE AT 55° T/Q (SEC) 3 1-2 0.6 0.4 o to 5 SEC. Fig. 34. Vertical and horizontal component synthetic seismograms for various T/Q wavelets of the form shown in Figure 32 for impulse source motions. STEPfSOURCE AT 55e T/Q (SEC)-12 0.2 5 SEC. o Fig. 35. Vertical and horizontal component synthetic seismograms for various T/Q wavelets for step source motions. 104 Using these computed wavelets f o r the input source to the l i n e a r systems model of the LED crust f o r 30° incidence (which corres-ponds to the 55° t e l e s e i s m i c distance) a set of synthetic seismograms have been c a l c u l a t e d (Figure 34). Step source synthetic seismograms are also presented (Figure 35). Comparisons of these syntheses with the LED Novaya Zemlya and Venezuela records i n d i c a t e s that the impulse source model f o r T/Q =0.4 seconds best represents the frequency character o f the i n i t i a l motions. (These two p a r t i c u l a r events are used i n t h i s comparison because t h e i r e p i c e n t r a l distances r e s u l t i n the 30° incidence angle used i n t h i s s y n t h e s i s ) . Carpenter suggested that f o r t e l e s e i s m i c P waves, T/Q should approximate unity. However, those syntheses f o r T/Q = 1.2 seconds and 0.6 seconds lack high frequency character observed i n the experimental records. Short duration source motions which could contain r e l a t i v e l y more high frequency character than an impulse are p h y s i c a l l y unreasonable. I t i s b e t t e r to r e t a i n the impulse model and require a lower T/Q =0.4 seconds value which appears to best f i t experiment. This T/Q value determines a mantle Q=1500. Q's i n the mantle exceeding t h i s (to 6000) f o r P waves with frequencies near 1 Hz have been reported by se v e r a l i n v e s t i g a t o r s (Jackson and Anderson, 1970) so that t h i s mantle absorption model i s not u n r e a l i s t i c . For a l l subsequent syntheses Q has been held constant to 1500 while T/Q v a r i e s with t r a v e l time. Synthetic seismograms for impulsive displacement sources at 55° and 85° e p i c e n t r a l distances are shown i n Figure 36. In these com-putations, the i n t e r n a l c r u s t a l attenuation has been included using the attenuator weights of Table IV. For the 55° distance, the t r a v e l time determines T/Q =0.4 seconds; f o r 85°, T/Q = 0.6 seconds. Impulse 55° EPICENTER 85° EPICENTER LEDUC CHAMULKA LARSEN U J — \ A A _ A A A ^ — . — 5 SEC. Fig. 36. Impulse source motion synthetic seismograms for 55° and 85° epicentral distances. The mantle Q has been chosen to be 1500 for a l l syntheses. Internal crustal absorption (causal) of Table II has been included. 106 source wavelets with, these mantle absorption model parameters were used as systems inputs. A comparison between a v e r t i c a l component non-attenuating crust seismogram and the corresponding absorbing crust seismogram showed l i t t l e d i f f e r e n c e . Only a s l i g h t reduction i n the apparent high frequency content was evident. 3. 17 Seismogram Syntheses Most of the recorded v e r t i c a l component seismograms possess e a r l y character which i s quite s i m i l a r to the equivalent syntheses. In p a r t i c u l a r , the apparent frequency of the f i r s t few cycles of the s y n t h e t i c responses compares w e l l to the frequency nature of the e a r l y part of a l l of the equivalent experimental seismograms. This i s not t r i v i a l , since the mantle Q assumption which a c t u a l l y determines the in c i d e n t wave frequency content was based only on the Novaya Zemlya and Venezuela records. The frequency content of the syntheses f o r the l a r g e r 85° e p i c e n t r a l distance also favourably compares with the experimental records. The h o r i z o n t a l records are not so v i s i b l y com-parable. This p a r t i a l l y r e s u l t s from the greater r e l a t i v e noise contamination of the low amplitude h o r i z o n t a l motions. Noise also reduces the p o s s i b i l i t y of good comparisons with low amplitude v e r t i c a l motions i n the l a t e r portions of the v e r t i c a l component records. I t i s apparent that a l l the recorded seismograms, a n d - p a r t i c u l a r l y those from the Novaya Zemlya, K u r i l e Islands and Ryukyu events, do not decay as r a p i d l y as the equivalent syntheses. Furthermore, the K u r i l e Islands and Ryukyu events show d i s t i n c t secondary a r r i v a l s (probably pP) within the f i r s t 25 seconds which suggests that the i n c i d e n t waveforms are longer and more complex than these simple source and mantle absorption 107 models can describe. I t would be d i f f i c u l t to create the t h e o r e t i c a l l y p o s s i b l e but complicated source motion models which would exactly describe the experimental records. Therefore, i t i s reasonable to r e s t r i c t any comparisons to the i n i t i a l few seconds of record which are dominated by the beginning of the incident waveform r e s u l t i n g from the onset of source motions. However, even wit h i n these f i r s t few seconds, many discrepancies between the theory and experiment are s t i l l obvious. In p a r t i c u l a r , the synthetic seismograms appear to exaggerate the amplitude of the f i r s t pulse of onset. This primary character should r e s u l t from the instrument response to the surface motions generated by the very f i r s t passage of the incident waveform through the crust to the surface. I n t e r n a l c r u s t a l reverberations can only s l i g h t l y a f f e c t these very e a r l y mtoi ons. Thus , i t becomes obvious that the source motion or mantle attenuation models are not completely adequate. Most probably, the source displacements are not as abrupt as the impulse and step displacement models. But, since the experi-mental seismograms are contaminated by background noise, p r e c i s e c o r r e l a t i o n s could not have been expected even from exact source motion models. In the following more d e t a i l e d comparisons between theory and experiment, the LAR recordings from the Venezuela event are not used because of t h e i r high apparent noise l e v e l s . 3.18 Comparisons of the Syntheses and Records for 55° Epicenter Neglecting the noisy LAR Venezuela records, a l l v e r t i c a l com-ponent seismograms show some motions which are c o r r e l a t a b l e with the COMPARISON OF EXPERIMENTAL AND SYNTHETIC SEISMOGRAMS FROM STATION LED (a) Venezuela '(e) F i j i Islands ' (b) Kurile Islands (c) Novaya Zemlya (d) Synthesis for 55° epicenter (f) Chile (g) Ryukyu Islands (h) Synthesis for 85° epicenter 5 SEC. o CO Fig. 37. A comparison of the vertical synthetic seismograms with the experimental seismograms of Figure 16 for station LED. equivalent s y n t h e t i c motions to beyond 4 seconds following onset. One such notable c o r r e l a t a b l e feature i s a strong, short negative pulse at 2 seconds which corresponds to a consistent experimental feature (compare Figure 36 and Figure 16a,d,e; see also Figure 37). This pulse could r e s u l t from the strong reverberation of the i n i t i a l motions i n v o l v i n g a r e f l e c t i o n from the free surface and the upper Precambrian boundary. (See Figure 30). Another broad p o s i t i v e pulse beginning at 1 second appears to c o r r e l a t e with s i m i l a r but more compli'cated experimental features. The synthesized pulses p a r t i c u l a r l y w e l l represent those of the K u r i l e Islands events i n both general character and r e l a t i v e ampli-tude. The LED Novaya Zemlya seismogram contains more character i n t h i s area (near 1 second) than the equivalent syntheses. However, the CHA synthesis and record compare s u r p r i s i n g l y w e l l . The s e v e r a l strong m u l t i p l e s which a r r i v e between 1 and 2 seconds (Figure 29) must l a r g e l y contribute to the formation of t h i s broad c h a r a c t e r i s t i c pulse. Apart from the preceding comparisons, i t i s also p o s s i b l e to c o r r e l a t e the v e r t i c a l component syntheses and the equivalent e x p e r i -mental seismograms on a peak to peak basis to almost 5 seconds following onset (compare Figure 37). However, the theory p r e d i c t s amplitudes f o r these l a t e r motions which are lower than those observed. This i n d i c a t e s that these l a t e r motions are at l e a s t p a r t i a l l y generated by some con-t i n u i n g i n c i d e n t wavemotion at the lower c r u s t a l boundary. Also, noise could contribute. L i t t l e s i g n i f i c a n t c o r r e l a t i o n i s obvious between any of the h o r i z o n t a l seismograms and t h e i r equivalent syntheses. 110 3. 19 Comparisons of the Syntheses and Records f o r 85° Epicenter Only the F i j i event v e r t i c a l seismograms d i s p l a y i n i t i a l character which d i r e c t l y corresponds to the equivalent syntheses (compare Figure 37). However, the C h i l e event records can also be made compatible by i n v e r s i o n ( m u l t i p l i c a t i o n by -1). This suggests that the source motions f o r t h i s event are b e t t e r represented by a d i l a t a t i o n a l impulse model rather than the compressional impulse model used. The Ryukyu seismograms possess such very complex e a r l y character that c o r r e l a t i o n s are somewhat ambiguous. The broad pulse feature at 1 second followed by the narrow pulse at 2 seconds i s also evident on these synthetic s e i s -mograms and again corresponds to s i m i l a r features of the F i j i and Chile records. For a F i j i c o r r e l a t i o n , the small amplitude precursor to the strong p o s i t i v e motion of the e a r l y record has been neglected. The h o r i z o n t a l motions from t h i s event also appear to compare quite favourably with the theory. However, the high background noise l e v e l makes a d e t a i l e d comparison d i f f i c u l t . 3. 20 Synthesis Using an Input Wavelet Obtained by Deconvolution Ulrych (1970) has recently developed a method of separating the convolved incident wavelet from the c r u s t a l impulse response function i n a seismogram using homomorphic non-linear f i l t e r i n g - ' methods (Oppenheim et a l , 1968). The technique requires no assumptions on the character of the wavelet except a choice of i t s representative length and requires no assumptions of c r u s t a l properties except l i n e a r i t y . Using h i s homomorphic deconvolution procedure on the LED Venezuela seismogram, Ulrych obtained a short (30 sample points) input wavelet which was used f o r the re-synthesis of the seismogram by the I l l VENEZUELA-LEDUC RESYNTHESIS I 5 SEC. F i g . 38. A comparison of the v e r t i c a l component LED seismogram f o r the Venezuela event (Event 1, Figure 16a) and the re-synthesis using the Ulrych wavelet. The i n t e r n a l c r u s t a l causal absorption of Table II has been included. l i n e a r systems approach. (Figure 38). The seismogram synthesized with t h i s wavelet provides a good t e s t of the crust model and method. I t i s not s u r p r i s i n g that the synthesized seismogram i s not i d e n t i c a l to the experiments. Such a f o r t u i t o u s r e s u l t would demand an exact crust model and wavelet. However, the remarkable degree of s i m i l a r i t y between the theory and experiment for the v e r t i c a l component i s extremely encouraging. The e a s i l y c o r r e l a t e d character extends to almost 6 seconds following the onset. Since the input wavelet was only 1.5 seconds i n length, i t i s obvious that these l a t e r motions are 112 generated by the c r u s t a l response. The synthesis s t i l l poorly represents the h o r i z o n t a l component. Because the h o r i z o n t a l motions are mostly converted SV energy which arises at the i n t e r n a l c r u s t a l boundaries, i t is apparent that s i g n i f i c a n t model discrepancies must e x i s t . 113 CHAPTER IV CONCLUSIONS 4.1 Summary of Analyses Between 0 and 2 Hz most of. the experimental V/H s p e c t r a l r a t i o s e x h i b i t general features which were t h e o r e t i c a l l y p redicted by the l i n e a r systems s o l u t i o n s . In particular,, expected 0.7 and 1.7 Hz r a t i o peaks were found. However, except f o r such low frequencies, the s p e c t r a l r a t i o comparisons between theory and experiment have revealed l i t t l e f u r t h e r confirmation of the chosen c r u s t a l models. The observed c o r r e l a t i o n s below 2 Hz suggest that on a scale determined by the long wavelengths of the low frequencies, the crust models are almost s u f f i c i e n t , but on a small s c a l e , s i g n i f i c a n t model discrepancies must e x i s t . The e f f e c t of attenuation was v i r t u a l l y indeterminable experimentally since the s p e c t r a l peaks which theory i n d i c a t e d should reveal c r u s t a l absorption properties were not consistent enough to derive any conclusions. However, the t h e o r e t i c a l s o l u t i o n s c o n c l u s i v e l y demonstrated that attenuation e f f e c t s are considerable beyond 2 Hz f o r a r e a l i s t i c c r u s t a l absorption model and must be included i n any comprehensive r a t i o a n a l y s i s . The- v i s i b l e e f f e c t s of causal c r u s t a l attenuation of the time synthesized seismo-grams appeared to be s l i g h t . T h i s , however, does not imply that. attenuation cannot be important at higher frequencies. The v e r t i c a l - t o - v e r t i c a l experimental r a t i o s showed l i t t l e character as was t h e o r e t i c a l l y p redicted below 3 Hz. Unfortunately, the s p e c t r a l d e t a i l necessary to d i s t i n g u i s h and i n t e r p r e t r e l a t i v e c r u s t a l l a y e r i n g and absorption properties does not appear to be a v a i l a b l e . 114 Considering the extremely simple assumptions on source motions and mantle properties which were used f o r the time domain syntheses, the s i m i l a r i t y between the e a r l y portions of the t h e o r e t i c a l and most of the experimental v e r t i c a l component seismograms was very encouraging. However, there was also a d i s t u r b i n g lack of c o r r e l a t i o n between the equivalent h o r i z o n t a l records and syntheses. The Ulrych homomorphic deconvolved input wavelet determined the best comparison between the t h e o r e t i c a l syntheses and experiment. • With these reasonable successes already achieved using the simplest of models f o r source motions, one could expect that b e t t e r conceived input incident wavelets could provide some remarkably accurate syntheses. The p o s s i b l e causes of any moderation of the success are v a r i e d . Background noise, inexact c r u s t a l models, s c a t t e r i n g and p o s s i b l e non-l i n e a r e f f e c t s could a l l contribute to discrepancies between theory and experiment. Among these p o s s i b l e f a c t o r s , s c a t t e r i n g e f f e c t s due to inhomogenieties i n the l a y e r s , non-plane i n t e r l a y e r boundaries and surface (topography) could i n v a l i d a t e t h i s simple l i n e a r systems s o l u t i o n for a r e a l sedimentary bedded cru s t . Chernov (Grant and West, 1965) has shown that the e f f e c t of s c a t t e r i n g i n an inhomogeneous medium increases r a p i d l y with frequency. The c o r r e l a t i o n s obtained at low frequencies i s not i n c o n s i s t e n t with t h i s s c a t t e r i n g hypothesis. Although, the chosen events appeared to possess large average s i g n a l - t o - n o i s e r a t i o s , back-ground noise dominates the P coda outside the s i g n a l band from 0.5 Hz to 3 Hz. This undoubtedly contributed to the lessened success at high frequencies. I t i s p o s s i b l e that inaccuracy and i n s u f f i c i e n t d e t a i l i n the c r u s t a l models has been attained using the i n t e r p o l a t i o n between the 115 a v a i l a b l e o i l w e l l d r i l l s i t e s . C o r r elations between theory and experiment may be g r e a t l y improved at high frequencies by using more accurate models which would require more complete knowledge of the c r u s t . Although the frequency domain computations can provide t h e o r e t i c a l l y exact s o l u t i o n s for any model, f u r t h e r assumptions which generate model inaccuracies are necessary for d i g i t a l time domain syntheses. For example, the time domain s o l u t i o n s require that the l a y e r t r a n s i t time f o r both P and SV motions i s an exact i n t e g r a l m u l t i p l e of the sampling increment. For the 0.05 second i n t e r v a l used i n the t h e s i s , the thinnest layers would e x h i b i t both P and SV t r a n s i t times p r e c i s e l y equal to t h i s i n t e r v a l . Such d i f f i c u l t i e s could e a s i l y be overcome by reducing the sampling increment to obtain any d e s i r e d l e v e l of accuracy. However, p a r a l l e l increases i n computer time and storage requirements must be t o l e r a t e d . 4.2 Concluding Remarks and Suggestions The l i n e a r systems s o l u t i o n of the m u l t i l a y e r propagation problem has already shown some p a r t i c u l a r l y u s e f u l features. I t becomes * the f i r s t mathematical procedure to e x p l i c i t l y include l a y e r attenuation for both normal and non-normal incidence although t h i s extension has been suggested f o r the standard techniques. The f a c t that the theor-e t i c a l P coda V/H r a t i o analyses f o r the crust and attenuation models used h e r e i n showed the s i g n i f i c a n t dependence on attenuation beyond 2 Hz demonstrates that c r u s t a l absorption cannot be neglected i n any complete propagation a n a l y s i s using short period P-coda. Al s o , since attenuation showed l i t t l e e f f e c t below 1 Hz, the previous V/H r a t i o work of s e v e r a l authors at low frequencies has been v a l i d a t e d with respect to attenuation. 116 The l i n e a r systems approach has the p a r t i c u l a r advantage that the time and frequency formulations are completely p a r a l l e l . This should a i d i n the conception of the r e l a t i o n s between the c h a r a c t e r i s t i c s of frequency and time s o l u t i o n s . Furthermore, two independent viewpoints should prove to be very u s e f u l . The formulation admits analysis and s o l u t i o n using the large body of communications theory mathematics. Such methods have, i n the past, been invaluable f o r the reduction and analysis of geophysical data. T h e i r further a p p l i c a t i o n to the d i r e c t formulation of the propagation problem may prove as u s e f u l by providing new i n s i g h t s and producing g e n e r a l i z a t i o n s of s o l u t i o n s . I t should be p o s s i b l e to extend solutions to include c r u s t a l p r o p e r t i e s which are non-linear. For example, appropriate non-linear operators which could describe p o s s i b l e n o n - l i n e a r i t i e s of c r u s t a l propagation could be i n s e r t e d i n t o the l i n e a r l y r e s t r i c t e d formulation and solved by the time synthesis approach. Since superposition i s not s t r i c t l y v a l i d f o r non-linear systems, previous methods which require frequency incremented or modal solutions would not be s u f f i c i e n t i n t h i s regard. A l s o , the s o l u t i o n of sets of non-linear d i f f e r e n t i a l equations d e s c r i b i n g non-linear propagation would probably be much more d i f f i c u l t than the non-linear extension of the time domain systems s o l u t i o n . I t may be d i f f i c u l t to conceive n o n - l i n e a r i t i e s i n i d e a l i n f i n i t e s m a l s t r a i n propagation. But, i n r e a l c r u s t a l materials where s c a t t e r i n g and absorption are known to occur, non-linear d e s c r i p t i o n s have been hypothesized (e.g., Walsh, 1966). Although no solutions f o r incident S waves were attempted i n t h i s t h e s i s , the l i n e a r systems formulation i s d i r e c t l y a p p l i c a b l e . 117 Furthermore, r e f l e c t i o n seismology could be considered by a t r i v i a l extension of the normal incidence transmission problem. However, i t i s important to note that the systems formulation demands plane waves so that point sources w i t h i n the crust cannot be allowed. Most other r e f l e c t i o n seismic s o l u t i o n s also require plane waves, normal incidence, and neglect the e f f e c t s of geometrical spreading which a c t u a l l y r e s u l t s from explosion sources. The systems formulation may provide a new viewpoint which would permit a general s o l u t i o n of the important inverse problem. For i d e a l l y layered c r u s t s , the l a y e r t r a n s f e r functions (see Appendix I) appear to contain a l l the information required to determine layer t r a n s i t times and r e l a t i v e acoustic impedance v a r i a t i o n s f o r both normal and non-normal incidence conditions. I t should be p o s s i b l e to obtain t h i s information from the equivalent seismogram. This approach, however, cannot be so unambiguous as to determine l a y e r v e l o c i t i e s , d e n s i t i e s and thicknesses. Inverse s o l u t i o n s are now being attempted. The author f e e l s that studies using the l i n e a r systems formulation could prove to be extremely f r u i t f u l f o r both more generalized propagation and inverse s o l u t i o n s . 118 REFERENCES Alberta Society of Petroleum Geologists. 1964. Geological h i s t o r y of western Canada. Anderson, D.L., Ben-Menahem, A. and Archambeau, C.B. 1965. Attenuation of seismic energy i n the upper mantle. J . Geophys. Res. 70, pp. 1441-1448. Basham, P.W. 1967. Time domain studies'of short period teleseismic P phases, M.A.Sc. Thesis, U n i v e r s i t y of B r i t i s h Columbia. Bohn, E.V. 1963. The transform analysis of l i n e a r systems. Addison-Wesley, Reading, p. 100.' ' • . Carnahan, B., Luther, H.A., and Wilkes, J.O. 1969. Applied numerical, methods. Wiley, New York. p. 80. Carpenter, E.W. 1967. Teleseismic signals calculated for underground, underwater and atmospheric explosions. Geophysics, _32, pp. 17-32. Claerbout, J.F. 1968. Synthesis of a layered medium from i t s acoustic transmission response. Geophysics, 3_3, pp. 264-269. Clowes, R.M. and Kanasewich, E.R. 1970. Seismic attenuation and the nature of r e f l e c t i n g horizons within the crust. J . Geophys. Res. 75, pp. 6693-6705. Cooley, J.W. and Tukey, J.W. 1965. An algorithm f o r the machine c a l c u l a t i o n of complex Fourier s e r i e s . Math, of Comp. 19, pp. 297-301. Cumming, G.L. and Kanasewich, E.R. 1966. C r u s t a l structure i n western Canada, F i n a l Rept., Contract AF 19 (628)-2835, AFCLR, Bedford, Mass Deas, A.T. 1969. A c r u s t a l study using teleseismic P phases recorded near Port Arthur, Ontario. M.Sc. Thesis, U n i v e r s i t y of B r i t i s h Columbia. E l l i s , R.M. and Basham, P.W. 1968. C r u s t a l c h a r a c t e r i s t i c s from short period P waves. Bull-. Seism. Soc. Amer. 58_, pp. 1681-1700. F r a s i e r , CW. 1970. Discrete time s o l u t i o n of plane P-SV waves i n a plane layered medium. Geophysics, _35, pp. 197-219. Futterman, W.I. 1962. Dispersive body waves. J . Geophys. Res. 67, pp. 5279-5291. Goupillaud, P. 1961. An approach to inverse f i l t e r i n g of near-surface l a y e r e f f e c t s from seismic records. Geophysics, 26, pp. 754-760. Grant, F.S. and West, G.F. 1965. I n t e r p r e t a t i o n theory i n applied geophysics. McGraw-Hill, New York. pp. 48, 200. 119 Haskell, N.A. .1953. The dispersion of surface waves on multilayered . media. B u l l . Seism. Soc. Amer. _43, pp. 17-34. Ibrahim, A.K. 1969. Determination of c r u s t a l thickness from s p e c t r a l behaviour of SH waves. B u l l . Seism. Soc. Amer. 59, pp. 1247-1258. Jackson, D.D., and Anderson, D.L. 1970. P h y s i c a l mechanisms of seismic-wave attenuation. 1970. Rev. Geophys. and Space Phys. _8, pp. 1-63. Knopoff, L. 1956. The seismic pulse In material possessing s o l i d f r i c t i o n , I: Plane waves. B u l l . Seism. Soc. Amer. 46_, pp. 175-183. Knopoff, L. and MacDonald, G.J.F. 1958. Attenuation of small amplitude stress waves i n s o l i d s . Rev. Mod. Phys. 3_0, pp. 1178-1192. Knopoff, Li and MacDonald, G.J.F. 1960. Models for acoustic los s i n s o l i d s . J . Geophys. Res.^65, pp. 2191-2197. Knopoff, L. 1964. Q. Rev. Geophys. 2_, pp. 625-660. Kuo, B.C. 1962. Automatic control systems. Prentice H a l l , Englewood C l i f f s , New Jersey, p. 40. Lindsey, J.P. 1960. E l i m i n a t i o n of seismic ghost r e f l e c t i o n s by means of a l i n e a r f i l t e r . Geophysics, 25_, pp. 130-140. McDonal, I.J. , Angona,' F.A.., M i l l s , R.L., Sengbush, R.L. , Van Nostrand, R.G. and White, J.E. 1958. Attenuation of shear and compressional waves i n P i e r r e Shale. Geophysics, 23, pp. 421-439. Mueller, S. and Ewing, M. 1962. Synthesis of normally dispersed wave trai n s by means of l i n e a r systems theory. Lamont Geological Observatory, Columbia U n i v e r s i t y , Contribution No. 559. Oppenheim, A.V., Schafer, R.W. and Stockham, T.G., J r . 1968. Non-linear f i l t e r i n g of m u l t i p l i e d and convolved s i g n a l s . Proc. IEEE 56, pp. 1264-1291. Phinney, R.A. "1964.. Structure of the Earth's crust from s p e c t r a l behaviour of long period waves. J . Geophys. Res. 69, pp. 2997-3017. Press, F. 1964. Seismic wave attenuation in.the crust. J . Geophys. Res. 69_, pp. 4417-4418. Richter, C F . 1958. Elementary seismology. Freeman, San Francisco, p. 670. Robinson, E.A. and T r e i t e l , S. 1966. Seismic wave propagation i n layered media i n terms of communication theory. Geophysics, 31, pp. 17-32. Savage, J . C. 1965. Attenuation of e l a s t i c waves i n granular mediums. J . Geophys. Res. 7_0, pp. 3935-3942. 120 Sherwood, J.W.C. 1962. The seismoline, an analog computer of t h e o r e t i c a l seismograms. Geophysics, 27, pp. 19-34. Thomson, W.T. 1950. Transmission of e l a s t i c waves through a s t r a t i f i e d s o l i d medium. J . Appl. Phys. _21, pp. 89-93.' T u l l o s , F.N. and Reid, A.C. 1969. Seismic attenuation i n Gulf coast sediments. Geophysics 34, pp. 516-528. Ulrych, T.J. 1971. A p p l i c a t i o n of homomorphic deconvolution to seismology. Submitted f o r p u b l i c a t i o n . Walsh, J.B. 1966. Seismic wave attenuation i n rock due to f r i c t i o n . J . Geophys. Res. 71, pp. 2591-2599. Ware, J.A. and Aki, K. 1969. Continuous and d i s c r e t e i n v e r s e - s c a t t e r i n g problems i n s t r a t i f i e d e l a s t i c medium. J . Acoust. Soc. Amer. 45, pp. 911-921. Werth, G.C. and Herbst, R.F. 1963. Comparisons of amplitudes of seismic/.' waves from nuclear explosions i n four mediums. J . Geophys. Res. 68, pp. 1463-1475. Williams, G.D. and Burk, C.F., J r . 1964. Geological History of Western Canada. Alb e r t a Society of Petroleum Geologists, Calgary. Chapter 12. Wuenschel, P.C. 1965. Dispersive body waves - an experimental study. Geophysics 30, pp. 539-551. 1 2 1 APPENDIX I IMPULSE RESPONSE FUNCTION FOR NORMAL INCIDENCE A l . 1 Laplace Inversion of Transfer Functions The normal incidence t r a n s f e r function s e t : H. (S) = , g<?> N • t ± V l-F(s)G(s) U i 2 ^ S ; l-F(s)G(s) V i H ( S) ° 2 ( s ) • t 'r ( A 1 _ 1 ) a±3KS) l-F(s)G(s) t i r i V ( " l - F ( s ) G ( s ) fci -sT where G(s) = A j,(s)e and F(s) = G C s ) ' ^ ' must be inverse Laplace transformed to obtain the c h a r a c t e r i s t i c l a y e r impulse response s e t : h (t) = L _ 1 { H (s)} . k ^ Rather than d i r e c t l y attempt to i n v e r t the closed forms above, i t i s e a s i e r to return to the s e r i e s representation and accomplish i n v e r s i o n term by term. This approach i s used, f o r example, i n transmission l i n e theory (e.g., Bohn, 1963, Chapter 11) because the closed forms are gen e r a l l y unsuitable f o r Laplace i n v e r s i o n s . A l s o , with the i n c l u s i o n of a general absorption function A ^ ( s ) , the pole set required f o r d i r e c t 122 i n v e r s i o n cannot be cal c u l a t e d . Thus f o r : -sT -3sT -5sT H ± (s) = t j L { A i ( s ) e Kk^Cs)* 1 r i r i , + A i 5 ( s ) e \ 2 * ± ' 2 + ••*> (Al - 2 ) the impulse response s e r i e s i s e a s i l y determined by term by term i n v e r s i o n : h. (t) = t . { a . ( t ) * 5 ( t - T . ) + a , ( t ) * a , ( t ) * a , ( t ) * S ( t - 3 T )r r.1 ,1^ i i i l I I i l l + a ± ( t ) * . . . * a i ( t ) * 5 ( t - 5 T i ) r i 2 r i 2 + ...} where a (t) = L~ 1{A ±(s)} Noting th a t : and H (s) = -r- H (s) we may fu r t h e r determine: ( A l - 3 ) - s T H i 2 ( s ) = ( r i ' A i ( s ) e > Hi/ s> t ' r —sT H i 3 ( s ) = ("VSi06 1 ) H i x ( s ) <M-*> h. (t) = r ' a. ( t ) * 6 ( t - T . ) * h . (t> (Al-5) l 2 1 1 l ±i 123 t *r. h, (t) = - f - i -a.(t)*5(t-T.)*h. (t) (Al-5) 13 ' ti 1 i l l t,1 and h, (t) = *h. (t) . iu t ± 11 For non-attenuating conditions, a^(t)=l and the four response functions may be represented in simplified i n f i n i t e series form: >i (t) = t ± I (r.^'T X5(t-[2n-l]T i) 1 n=l h. (t) = t ± r ± I (r.r.'V^Vt^nT^ ±2 n=l (Al-6) and h i ( T ) = V r i ^ ( r i r i * ) n " 1 6 ( t - 2 n T i ) 3 n=l h. (t) = t.' I (r.r » ) n 16(t-[2n-l]T ) 1(* 1 n=l 1 1 1 Al.2 Poles and Zeros of the Normal Incidence Transfer Functions The four transfer functions (equation Al-1) each have a common denominator which determines the pole set at complex frequencies for which: -2sT l - ( r i r i , ) A . 2 ( s ) e 1 = 0 (Al-7) 124 Since the A^(s) can represent any general attenuation f u n c t i o n , no s p e c i f i c pole set can be determined without assumptions on A ^ ( s ) . For the non-attenuating case, A^(s) = 1 . A p h y s i c a l constraint (conservation of energy) demands that the r e f l e c t i o n c o e f f i c i e n t s at the two e l a s t i c boundaries of l a y e r i obey: < 1 (Al-8) s i s a complex number and the exponential must obey the p e r i o d i c r e l a t i o n s h i p : -2sT i -2sT i i2nir e = e • e (Al-9) where n i s any a r b i t r a r y i n t e g e r . Thus a s o l u t i o n of equation (Al-2) follows: 1 (Al-10) e Taking the n a t u r a l logarithm of t h i s equation: 2 ( 1 ^ - 3 ^ ) = - l n ( r r ^ ) (Al-11) and there are an i n f i n i t e number of s o l u t i o n s for s: S r = ^-{inir+ -|ln(r r ^ ) } (Al-12) Two conditions on the co e f f i c i e n t s y i e l d d i f f e r e n t general pole patterns: Case I - 0 < r . r . ' < 1 i l Thus l n ( r ^ r ^ ' ) i s negative r e a l and the pole positions are shown i n Figure 39. The poles are periodic with the rate to = TT/T. (f = ——) i 2T 1 1 the f i r s t non-zero frequency pole i s at f j = . + s - / LU 2 T £ ;lniVi'l Fig.. 39. Pole maps of the one layer transfer function for normal incidence: Case I. 126 Case II 1 < r . r . ' < 0 l l Set r ^ r . ' = [ r . r . ' l e * 7 1 . Thus l n ( r . r . ' ) i s complex with p r i n c i p a l value: l n ( r . r . ' ) = l n l r . r . ' l+iir . i i 1 l l 1 (Al-13) The poles = ~ {i(n+l/2)+l/21n| r | } are shown i n Figure 40. Again, the poles are p e r i o d i c with rate f =' 2 T ± ' pole has the frequency component f = 4T, However, the fundamental. These fundamental poles +s +s + S; +s 7 LU 2T £ F i g . 40. Pole maps of the one l a y e r t r a n s f e r function f o r normal incidence: .Case I I . 1 2 7 at fi - o r \ f o r Case I and f = — f o r Case I I represent the fundamental i ° i reverberation resonance frequencies due to the l a y e r ; the p e r i o d i c poles represent the harmonics. Each t r a n s f e r function has zeros at i n f i n i t e s. For non-normal incidence, the computation of poles and zeros i s much more d i f f i c u l t because the nature of the denominator of the tr a n s f e r function i s h i g h l y complicated (Appendix I I I ) . However, s i m i l a r general features should p r e v a i l f o r incidence angles which are not l a r g e . 128 APPENDIX II PHASE ADVANCE TRANSFER FUNCTIONS Consider any l a y e r i w i t h o r i g i n as shown i n F i g u r e 41. The d i r e c t P wave ray which e x i t s the upper boundary at, x=0, z=n^ enters the lower boundary at x=a, z=0. The phase at a given i n s t a n t at ( a,0) is time advanced r e l a t i v e to the phase at po i n t P and hence at the o r i g i n (0,0) according to the t r i g o n o m e t r i c r e l a t i o n : x = ( n . t a n ( e . ) s i n ( e . ) ) / a . . (A2-1) p. i l i i LAYER i (a ,0) (b,0) F i g . 41. A r e p r e s e n t a t i o n of phase advance time c a l c u l a t i o n s necessary f o r non-normal c o n d i t i o n analyses. 129 Similarly for the direct S wave entering at (b,0), the time advance between (b,0) and point Q i s : T = (n tan(f )sin(f.))/g. . (A2-2) s^ 1 1 1 1 These advance times represent any one transit from one boundary to the other for the P and S waves respectively. Higher orders of reflections are included through further multiplication by identical advance functions. The included advance transfer functions of layer in Figure 8 exactly establish the correct phase relations for any possible com-binations of internal reverberations. The phase advance transfer functions and the layer delay functions may be combined into single delay functions: -sT' -sT* e and e n i where T' = T -x = — c o s ( e J (A2-3) P ± P ± P ± a ± i n i and T' = T -x =-r icos(fJ . (A2-4) S i S i S i B i 1 1 3 0 APPENDIX I I I EXAMPLE BLOCK REDUCTION PROCEDURE I t i s de s i r a b l e to obtain a s i n g l e t r a n s f e r function H(s) to represent a system component (Figure 42) such that V (s) = H(s)V.(s). o 1 Vj (s) V0(s) F i g . 42. The desired block diagram representation i n terms of t r a n s f e r function and transformed s i g n a l s . Two common system block diagrams can be e a s i l y conformed to t h i s form using the standard block diagram reduction methods (Kuo, 1963). Case I - Summed p a r a l l e l systems (Figure 43): H2(s) F i g . 43. Block diagram of p a r a l l e l summing paths. 131 V (s) = (H 1(s)+H 2(s))V.(s) O X (A3-1) Set: H(s) = HiCsJ+^Cs) (A3-2) Case II - Continuous (feedback systems) loops (Figure 44) Hj(s) V S > = l - H 1 ( s ) H 2 ( s ) * V S > (A3-3) Set: H(s) = H!(s) l - H 1 ( s ) H 2 ( s ) (A3-4.) H[(s) r ' H2(s) V0(s) F i g . 44. Block diagram of the l i n e a r p o s i t i v e feedback system. 132 Example Reduction: Refer to Figure 8 during the following reductions. Given the equation: U, (s) = H. (s)U.(s) (A3-5) 1—JL X which relates an up travelling P wave entering layer i to the equivalent up travelling P wave entering layer i-1, determine H i ! ( s ) = U i _ l ( s ) / U i ( s ) ' ( A 3 _ 6 ) It Is first necessary to determine the P and SV motions incident on the upper boundary due to the P and SV motions at the lower boundary. For simplicity, let the forward absorption and delay transfer functions be represented by (see Appendix I I for definition of T 1 and T '): Gi(s) = Ap (s)e . (A3-7) G2(s) = A (s)e 1 (A3-8) S i and the return transfer functions: B!(s) = r r .G^s) B 2(s) = r r ,G2(s) (A3-9) ps ± sp i' z B 3(s) = r r ,Gi(s) S p i p p i 133 Bi*(s) = r r ,G 2(s) s s ± s P ± B 6 ( s ) = r r .G^s) (A3-9) PP ± P s ± B 7 ( s ) = r r ,G 2(s) s s ^ s s ^ Combining the p a r a l l e l t r a n s f e r functions (Case I) the s i m p l i f i e d block representation of Figure 45 i s p o s s i b l e . Further combining Gj^s) with i t s return path t r a n s f e r function and G 2(s) with i t s return path t r a n s f e r f u n c t i o n form: G X(s) F l ( s ) = l - G 1 ( s ) ( B 1 ( s ) + B 2 ( s ) ) (A3-10) G 2(s) ? 2 ( S ) = i - G 2 ( s ) ( B 7 ( s ) + B 8 ( s ) ) Since i t i s intended only to obtain the l a y e r response to the i n c i d e n t P wave, the only non-zero input s i g n a l i s U (s) . The t o t a l open loop t r a n s f e r f unction (obtained by breaking the continuous loop of Figure 44) i s F 1 ( s ) F 2 ( s ) ( B 3 ( s ) + B ^ ( s ) ) ( B 5 ( s ) + B 6 ( s ) ) . The open loop and forward t r a n s f e r functions combine to form (Case I I ) : 134 B 3 (s ) + B 4 (s) t : V , ( s ) & U'(s) V'(s) Fig. 45. The reduced system block diagram having the layer response to inputs U^(s) and V^(s). U'(s) V'(s) uTTsj" F i ( s ) l - F 1 ( s ) F 2 ( s ) ( B 5 ( s ) + B 6 ( s ) ) ( B 3 ( s ) + B 4 ( s ) ) (A3-11) F i ( s ) F 2 ( s ) ( B 5 ( s ) + B 6 ( s ) ) l - F 1 ( s ) F 2 ( s ) (B 5(s)+B 6(s)) (B 3(s)+B t t(s)) The. f i n a l transmission of these P and SV motions incident on the upper boundary through the upper boundary determines U\ ^ ( s ) : 135 U, ,(s) = t U*(s)+t V'(s) (A3-12) i- 1 pp ± s p 1 so that: U i - l ( s ) F l ( s ) ' t D D + F l ( s ) F 2 ( s ) ( B 5 ( s ) + B 6 ( s ) ) - t H, (s) = 1 1 U ±(s) l - F 1 ( s ) F 2 ( s ) ( B 5 ( s ) + B s ( s ) ) ( B 3 ( s ) + B l t ( s ) ) (A3-13) Each t r a n s f e r function H (s) can be'obtained i n a s i m i l a r manner, •k A l l 16 t r a n s f e r function possess the same denominator. 136 APPENDIX IV FORTRAN COMPUTER PROGRAMS FOR TIME DOMAIN SEISMOGRAM SYNTHESES AND FOURIER COMPONENT SOLUTIONS BASED ON THE LINEAR SYSTEMS THEORY The following FORTRAN IV-G programs were used f o r the numerical solutions of the transmission problems i n t h i s t h e s i s . I: NON-NORMAL INCIDENCE CRUSTAL FOURIER TRANSFER FUNCTION PROGRAM This program was used f o r the frequency domain t r a n s f e r function s o l u t i o n s . I I : NON-NORMAL INCIDENCE SYNTHETIC SEISMOGRAM FOR ANGLES LESS THAN CRITICAL This program was used f or the time domain seismogram syntheses f o r non-normal incidence conditions. On an IBM System 360 Model 67 computer operating under the Michigan Terminal System, the computational time was approximately 0.005 second per l a y e r per time increment f o r non-attenuating c r u s t models. The computational time was s l i g h t l y increased f o r attenuating l a y e r models. • I I I : SUBROUTINES A l l subprograms required by the above two programs are following. •»Uri"iMOKMAl I N C I D E N C E C R U S T A L F O U R I E R T R A N S F S * . - U N C T I O N P R O G R A M • * : ; ; * * * * * * * * * * * * * ; | : : ; : * * * : | : * * * : 1 : * * * * * * * * * * * i \ t * :|c * * * * * * * * * * 5;s * * *i\: * * * * * :\: PROGRAM DESCRIPTION TOE PROGRAM USES THE FPJHUJENCY DOMAIN LINEAR SYSTEMS CRUSTAL MODEL FOR COMPUTATION Oh SAMPLED F RE OU EM CY C O M P O N E N T FOUR t E R TRANSFER FONCTIONS PROGRAM LIMITATIONS •i"«X 11 H !i i MO. OF LAYERS-10 MAXIMUM MO. O F FREQUENCY POINTS C A LCULABL F=3 0 '\ CAUSAL 7. E RO'-'P HAS t ATTENUATION IS ASSUMED «•«••> THE PROGRAM MAY J -. GENERALIZED TO CAUSAL ATTENUATION BY ALLOWING COMPLEX ABSORPTION FNCS. 'INPUT CARDS . CARD 1: TITLE CARD 2 : NO. OF LAYERS,. NO. OF FREQUENCY POINTS* I N I T I A L FREQUENCY', FINAL FRF.OUF.NCY CARD 3: OUTPUT LIST AND ATTENUATOR FLAGS SKIP OUTPUT L I S T NFLAG 1 = 0 S <IP ATTENUATION NFL AG 2 = 0 CARD 4: FLAGS TO CHOOSE PROGRAM OPTIONS "••« . -S U P RAT 10 PROGRAM MR/\TIO = 0 SKIP DISPLACEMENT PROGRAM NOISPL=0 CARD 5 FF: LAYER THICKNESS, P VELOCITY, S VELOCITY, DENSITY, 0 FO< 0 FOR SV. IF THE ATTENUATOR FLAG = 0 S K I P THE it SPEC IF ICAT IONS CARD 6: P OR S-'WAVE INCIDENCE ANGLE VII TH NORMAL «•» THE INPUT NAVE TYPE IS DETERMINED BY WHICH ANGLE I S N U W . F . R O . FOR NORMAL IMC IDE N C F VERY SMALL INCIDENCE ANGLES MAY USED W I T H O U T DETECTABLE i=RR-') I 01 MEMS I Or! AL P ( 11 ) , B ET ( 1.1 ) ,RHO( 11 ) , II ( 10 ) , : ( 11 ) , I- ( 11 ) , T P P ( 2 1 ) ,TPS( 21 1) ,TSP( 21 ) ,TSS ( 21 ) ,RSP(21 ) ,RSS ( 21 ) ,RPS( 21 ) ,RPP( 21 ) , FR F0{ 50) , VP( 11 ) 2,VS(11) ,TP(10 ) , T S ( 1 0 ) , A P ( 1 0 ) , A S ( 10),E X P P( 10) ,E X P S( ID) , S T A T ( 1 5 ) , 0 X ( 311) ,W7. ( 11 ) , PANGLE ( 11 ) , S ANGL E ( 11 ) , A ( 4 0, 41 ) , ii { 40,41 ) , X( 40, 1 ) , S INE( 11 4) , S TMF( 11 ) , COSE ( 11 ) ,COSF { 11 ) ,TAUP( 10 ) , TAILS ( 10 ) , R A T ( SO ) ,OAMP( 11, 50 ) 5,0PHS(11,50) ,WAMP(1L,50) ,WPHS(11,50),0P(10) ,0S( 10) COMPLEX EXPP,EXPS ,G1 ,G2 ,A ,H1 ,H2,H3 ,H4 ,X ,DET, B, 0 1 , 02 , 03,04,0 5,06,07 1,D8,P1 ,S1 ,0X ,W/ DATA P 1/3.141 5 927/ READ TIT'. Ir CARD: BO CHARACTERS TO BE PRINTED ON OUTPUT RE A I'M 5,2) STAT •I :•! I T E ( 6 ,2 ) S TA T C READ INPUT: NUMBER 1.1 F LAYERS NUMBER OF FREQUENCY P f l l NTS » I NI T IAL FREOIJENCY C A NO FINAL RE DO EN I Y R E A I") ( 5 ,1 () ) N »I" M A X , E R E Q ( 1 ) , F R E Q ( (•'! M AX) i.|P=N+l l-INC = ( F R B P h i M A X ) - F R E O d ))/FLOAT (MHAX-1) OO I I=2,Mi-1/\X 1 . F R E 0 ( I ) =FRE'> ( I«l ) +F INC C READ I-.'PUT: FLAGS TO DETERMINE. OPTIONS KEAO(5,3) N FL A Gl , r-l FLA G2 REAI'X ii ,3 ) MR AT 111,WO IS PL IF(M R AT 10 .E0..0 .AND .NO IS PL .EO.O) STOP -R IT E ( 6 ,2 9 ) READ(5,11) ANGINP , ANG I MS IF(MFLAG5.EO .0) GO TO 13 C READ INPUT: LAYER THICKNESS, P VELOCITY, S VELOCITY, DENSITY, ATTENUATION C CONSTANTS FOR P AMD SV MOTIONS REA 0 ( 5 , 1 2 ) ( H ( 1 ) , A L P ( I ) » B E T ( I ) ,RHO(I ),OP(I ),DS( I ) ,I = 1,N ) GO TO 111 C READ INPUT: LAYER PARAMETERS WITHOUT ATTENUATION CONSTANTS "=> THIS OPTION C IS USED WHEN ATTENUATOR FLAG IS SET TO 0 13 READ(5,121) ( H ( I ) , A L P ( I ) , B E T ( I ) , RHO(I ),I = 1,N) 111 READ!5,121) DUMMY,ALP (MP.),BET(NP),RHO(M?) WRITE(6,14) ( I , A L P ( I ) ,BET(I ) ,RHO(I ) ,H(I ),I=1,M) WRITE(6,14) NP,ALP(MP),BET(MP),RHO(NP) EOO IVALENCE (VP,ALP ) ,(VS,8ET ) C REFLECTION AND TRANSMISSION ZOEPPRITZ' COEFFICIENTS ARE CALCULATED Ai-!GIi-]P=A!v'GINP=:--PI/180. ' ' ANGINS = ANGTMS-:;P 1/180 . CALL COEGEN(N,VP,VS,RHO,AMGINP,ANGINS,E,F,TPP,TPS,RPP,RPS,TSP,TSS, 1RSP,RSS) C LAYER TRANSIT TIMES AND PHASE ADVANCE TIMES ARE CALCULATED 0 0 5 1=1,N COS E( I ) = C O S ( E ( I ) ) S INE( I ) = S IN( E( t ).) ' COS F( I ) = C O S ( F ( I ) ) S INF( I ) = S IN( F'( t ) ) TPt I )=H( I ) /(V P ( I )*COSE(I )) TS( I )'=H( I ) / ( V S ( I ) * 0 O S r ( I )) TAOP( I) =H( I )*SINE ( I )*SINE (I ) / (VP { I )*COSE ( I ) ) T/MJS( I ) = h!( I )*S INF( I ) *S I N F { I ) / ( V S ( I )*COSF(I )) 5 CONTINUE C OS E ( N P) = C OS (E (N P ) ) S INE( NP ) =S INI E (NP ) ) 0 OSE( N P ) = COS (i- ( N P ) ) SINF(iMP) =SIN( F (IMP ) ) C CALCULATION OF TRANSFER MATRIX VS. F RF OUENC/ IF ( NFLAG2 .NE .0 ) GO TO 221 ' |J() 222 I = l , i M . C ABSORPTION FUNCTIONS ARE SET TO ONE FOR NO ABSORPTION A P ( I ) = 1. 222 A S( I)=1 . 221 CONTINUE 00 99 H=1,MMAX W = PI*i-RE(.)(M)*2. C CALCULATION OF ABSORPTION FUNCTIONS OF FREQUENCY 1 F( MFLAG2 .EC) .0 ) GO TO 24 00 22 1=1,N A S( I)=EXP(»0,5*W*TS(I )/QS(I )) 22 . AP( I)=EXP(«0.5*W*TP(I ) / P P ( I ) ) 24 CONTINUE C CALCULATION OF OF DELAY TRANSFER FUNCTIONS AT FREQUENCY COMPONENT 0 0 4 1 = 1,1-1 A1.= W * ( TP (I )™TAUP ( I ) ) A 2=U* ( TS ( I ) MTAUS ( I ) ) E X P P( I)=C M PL X(COS(A 1 ) , " S I N ( A 1 ) ) E XPS ( I)=CMPL X ( C O S ( A 2 ) , " S I N ( A 2 ) ) 4 CONTINUE <HAX = 4*N LMA XP =K MAX+ 1 0 0 2 0 K=l ,KMA'X DO 20 L = 1,LMAXP 20 A(K,L) =0. . M Ui C F G R i-i HA I N D I A G O N A L O F M A T R I X A L L V A L U E S = " 1 . D O 2 1 K=l , K M A X 2 1 A ( K , K ) = « 1 . C C O M P O T E R E M A I N I N G M A T R I X C O M P O N E N T S I F ( N . E O . l ) GO T O 2 3 0 0 A I = 2 , N G 1 = A P ( I ) * E X P P ( I ) G 2 = A S ( I ) * E X P S ( t ) J = I + N P C S O i i R O O T I N E .TRAMS. I . C O M P O T E S T H E C O M M O N P A R T S U F T H E L A Y E ' T R A N S F E R F U N ; i " I O N S C F O R O P " T R A V E L L I N G O U T P U T S I G N A L S C A L L T R A N S l ( R P J ( I ) , R P S ( I ) , R S P ( I ) , R S S ( I ) , * P P ( J ) , R P S ( J ) , R S P ( J ) , R S S ( J 1 ) , G l , G 2 , H 1 , H 2 , H 3 , H 4 ) 0 1 = T P P ( I ) * H L + T S P ( I ) * H 3 D2 = T P P ( I ) * H 2 + T S P ( I ) * H 4 D 3 = T P S ( I ) * H 1 + T S S ( I ) * H 3 D 4 = T P S ( I ) * H 2 + T S S ( I ) * H 4 D 5= R P P ( J ) * G 1 D 6 = R P S ( J ) * G 1 0 Y = R S P ( J ) * G 2 D 8 = R S S ( J ) * G 2 < = 4 * I - 7 L = K + 4 .<P=K + 1 L P = L + 1 C I N P U T - O U T P U T R E L A T I O N S A R E I N D I C A T E D BY T H E F O L L O W I N G C O M M E N T C A R D S C A ( l)( I'-l ) / U ( I ) ) A ( K , L ) = D 1 C A ( 0 ( I " 1 ) / V ( I ) ) A ( K , L P ) = 0 2 C A ( V ( I»l)/0( I ) ) A ( K P , L ) = D 3 C A ( V ( I " l ) / V ( I ) ) A ( K P , L P ) = D 4 L = K + 2 L P = L + 1 C A ( U . ( I » l ) / D ( I ) ) A ( K , L ) = 0 5 * 0 1 + 0 6 * 0 2 C A ( 0 ( T "1 ) / E ( I ) ) A ( K , L P )=D7*D1+08*02 C A ( V ( I ~ l ) / D ( I ) ) A(K P , L ) = D5*D3+06*D4 C A(V( I"1 ) / E ( I ) ) A(KP,LP)=D7*D3+D8*04 • ., C SUBROUTINE TRANS 1 COMPUTES THE COMMON PARTS OF THE LAYER TRANSFER FUNCTIONS C FOR DOWN-TRAVELLING OUTPUT SIGNALS CALL TRANS 1 ( R P P ( J ) , R P S ( J ) , R S P ( J ) t R S S ( J ) , R P P ( I ) , R P S ( I ) , R S P ( I ) , R S S ( I 1) ,G1,G2,H1,H2,H3 ,H4) 01=TPP(J)*H1+TSP(J)*H3 D2=TPP (J)*H2+TSP( J )*H4 D3=TPS( J) *H1+TSS ( J )*H3 D4=TPS(J)*H2+TSS(J)*H4 D5=RPP(I)*G1 Q6=RPS(I)*G1 •D7 = RSP( I ) *G2 08=R S S ( I ) * G 2 K=4*I-1 L = K=4 KP=K+1 LP=L+1 C A ( D ( I + l ) / 0 ( I) ) A(K,L)=D1 C INPUT«OUTPUT RELATIONS ARE INDICATED* BY THE FOLLOW ING COMMENT CARDS C A ( D{ I + l )/E ( I ) ) A(K,LP)=D2 C A ( E ( I +1 ) / 0 ( I )) A(KP,L)=D3 C A ( E ( I +1) / E ( I )) A(KP,LP)=D4 L = K™2 LP=L+1 C A(D( I + l )/U(I )) A( K ,L)=D5*D1.+ D6*02 ' C A ( 0 ( I + 1 ) / V ( I ) ) A( !< ,LP) =07*01 +08*02 C A ( E ( I + 1 ) / U ( I ) ) • • A ( K P,L)=05 *D3+06*04 C A ( E ( I + 1 ) / V ( I ) ) A l I \ K , L P J = U / * D 3 + D 8 * D 4 6 C O N T I N U E 2 3 C O N T I N U E C C A L C U L A T E T R A N S F E R M A T R I X C O M P O N E N T S F O R L A Y E R 1 1=1 J = 1 + N P G1 = A P ( 1 ) * E X P P ( 1 ) G 2 = A S ( 1 ) * E X P S ( I ) C A L L T R A N S 1 ( R P P ( J ) , R P S ( J ) , R S P ( J ) R S S ( J ) T R P P ( I ) , R P S ( I ) , R S P ( I ) , R S S ( I 1 ) , G 1 , G 2 , H 1 , H 2 , H 3 , H 4 ) 0 1 = T P P ( J ) * H 1 + T S P ( J ) * H 3 D 2 = T P P ( J ) * H 2 + T S P ( J ) * H 4 ' . D 3 = T P S < J ) * H 1 + T S S ( J ) * H 3 D 4 = T P S ( J ) * H 2 + T S S ( J ) * H 4 D 5 = R P P ( I ) * G 1 D 6 = R P S ( I ) * G 1 D 7 = R S P ( I ) * G 2 0 8 = R S S ( I ) * G 2 C A ( D ( 2 ) / 0 ( 1 ) ) A ( 3 , l ) = D 5 * D L + D 6 * D 2 C A ( 0 ( 2 ) / V ( l ) ) A ( 3 , 2 ) = 0 7 * 0 1 + 0 8 * 0 2 C A ( E ( 2 ) / 0 ( l ) ) A ( 4 , I ) = 0 5 * 0 3 + D 6 * D 4 C . A ( E ( 2 ) / V ( l ) ) A ( 4 , 2 ) = 0 7 * 0 3 + 0 8 * 0 4 K = 4 * N " 3 K P = K + 1 K P P = K P + 1 < P P P = K P P + 1 L = K + 4 C T H E C O N S T A N T V E C T O R D E T E R M I N I N G I N P U T C O N D I T I O N S I S C O M P U T E D I F ( A M G I N S . E O . O . ) G O T O 7 I F ( A N G I N P . 1 : 0 . 0 . ) GO T O 8 A ( U ( N ) / U ( N + 1 ) ) ^ A ( K » L ) = ™ T P P ( N P ) A ( V ( N ) / U ( N + 1 ) ) A ( K P t L ) = " T P S ( N P ) A ( D ( N + 1 ) / U ( N + 1 ) ) ' ' . A ( K P P , L ) = « R P P ( N P ) A ( E ( M + 1 ) / U ( N + l ) ) A ( K P P P , L ) = « R P S ( N P ) GO TO 9 A ( U ( N ) / V ( N + 1 ) ) A ( K ,1 , ) = « T S P ( N P ) A ( V ( N) / V ( N + 1 ) ) A ( K P , L ) = - T S S ( N P ) • • • A ( D( N + 1 ) / V ( N + l ) ) A ( K P P , L ) = " R S P ( N P ) A ( E ( N+ 1 ) / V ( N +1 ) ) A ( K P P P , L ) = « R S S ( N P ) C O N T I N U E DO 3 0 K = l f K M A X D O 3 0 1 = 1 , L M A X P B ( K , L ) = A ( K , L ) C A L L G A U C 0 M ( K M 4 X , 8 , X , D E T , 4 0 , 4 1 ) 1=1 J = 1 + N P C A L L T R A N S 1 ( R P P ( I ) ,RPS ( I ) ,RSP U ) , R S S ( I ) , R P P ( J ) ,RPS'.( J ) ,R.SP( J ) , R S S ( J 1 ) , G I , G 2 , H 1 , H 2 , H 3 , H 4 ) C O M P O T E T R A N S F E R F U N C T I O N T H R O U G H U P P E R t A Y - R T O F R E E S U R F A C E B O U N D A R Y P 1 = H 1 * X ( 1 t l ) + H 2 * X ( 2 , l ) S 1 = H 3 * X ( 1 , 1 ) + H 4 * X ( 2 , 1 ) C A L C U L A T I O N O F S U R F A C E D I S P L A C E M E N T S — C O M P L E X D 1 = P 1 * R P P ( 1 ) D 2 = S 1 * R S P ( 1 ) D3 = S 1 * R S S ( I ) D 4 = P 1 * R P S ( 1 ) 0 5 = ( P 1 + D 1 + D 2 ) * C O S E ( 1 ) ' , ' • D 6 = ( P 1 - D 1 - D 2 ) * S I N £ ( 1 ) D 7 = ( S 1 + D 3 + D 4 ) # S I N F ( 1 ) 0 8 = ( « S 1 + 0 3 + 0 4 ) * C O S F ( 1 ) U X ( 1 ) = D 6 + D 8 WZ{ 1 ) = D 5 + D 7 I F ( H D I S P L . E O . O ) G O T O 3 2 1 C C O H P U T A T I U N O F I N T E R I O R B O U N D A R Y D I S P L A C E M E N T S — C O M P L E X DO 3 2 1 = 2 » N P ' I ••!= I " l L 4 = 4 * I M L 3 = L 4 " 1 L 2 = L 3 ' : ' l L = L 2 " 1 D 1 = X ( L , 1 ) * C O S E ( I M ) + X ( L 3 , l ) * C O S E ( I ) D 2 = X ( L , 1 ) * S I N E ( I M ) - X ( L 3 , 1 ) * S I N E ( I ) D 3 = X ( L 2 , 1 ) * S I N F ( I M ) + X ( 1 4 , 1 ) * S I N F ( I ) D 4 = = X ( L 2 , 1 ) * C O S F ( I M ) + X ( L 4 , l ) * C O S F ( I ) OX'( I ) = D 2 + D 4 3 2 W Z ( I ) = 0 1 + 0 3 C C O M P U T E S U R F A C E A N D I N T E R N A L B O U N D A R Y D I S P L A C E M N T S V S . F R E Q U E N C Y I N T E R M S O F C A M P L I T U D E A N O P H A S E C O M P O N E N T S D O 3 3 1 = 1 , N P U X I = A I MA G ( U X ( I ) ) U X R = R E A L ( U X ( I ) ) ' • W Z I = A I MA G ( W Z ( I ) ) WZR = R E A L ( W Z ( I ) ) U A M P ( I , M ) = S Q RT ( U X R #U X R + U X I #U X I ) W A M P t I , M ) = S O R T ( W Z R * W Z R + W Z I * W Z I ) ' U P H S t I , M ) = A T A N 2 ( U X I , U X R ) * 33 W P H S ( I , M ) = A T A N 2 ( W Z I , W Z R ) 0 0 2 5 1=1 , N P P A N G L E ( I) = E U ) * 1 8 0 . / P I S A N G L E ( I ) = F ( I ) * 1 8 0 . / P I 2 5 C O N T I N U E ) I F ( N F L A G 1 . E O . O ) G O T O 2 6 W R I T E ( 6 , 1 8 ) F R E O ( M ) W R I T E ( 6 , 1 5 ) ( ( A ( K , L ) , L = 1 , L M A X P ) , K = 1 , K M A X ) W R I T E ( 6 , 1 5 ) ( X ( L , 1 ) , L = 1 , K M A X ) 2 6 C O N T I N U E . , " ' W R I T E ( 6 , 1 7 ) . M W R I T E ! 6 , 1 6 ) ( I , U X ( I ) , W Z ( I ) , U A M P { I , M ) , U P H S < I , M ) , W A M P t ( ,M> , W P H S ( I, M) 1 , P A ' N G L E t I ) , S A N G L E ( I ) , 1=1 , N P ) W R I T E ( 6 , 1 9 ) O E T 2 7 C O N T I N U E 3 2 1 I F ( N R A T I O . 6 O . 0 ) GO T O 3 2 2 C C O M P U T A T I O N O F S P E C T R A L A M P L I T U D E R A T I O S F O R S U R F A C E M O T I O N S - - I N T E R N A L C B O U N D A R Y S P E C T R A L R A T I O S C O U L D A L S O 8 E O B T I N E D C B O U N D A R Y S P E C T R A L R A T I O S C O U L D A L S O B E O B T A I N E D U X I = A I M A G ( U X (1 ) ) U X R = R E A L ( U X ( 1 ) ) W Z I = A I M A G ( W Z ( 1 ) ) W Z R = R E A L ( W Z ( 1 ) ) R A T ( M ) = S O R T ( W Z R * W Z R + W Z I * W Z I ) / S O R T ( U X R * U X R + U X I * U X I ) 3 2 2 C O N T I N U E 9 9 C O N T I N U E I F ( N R A T I O . N E . O ) W R I T E ( 6 , 1 1 0 ) I F ( N R A T I O . N E . 0 ) W R I T E ( 6 , 1 1 ) ( F R E Q ( M ) , R A T ( M ) , M = 1 , M M A X ) S T O P 2 F O R M A T ( 1 5 A 4 ) 3 F O R M A T ( 5 1 1 ) 1 0 F O R M A T ( 2 1 2 , 2 F 6 . 2 ) 1 1 F O R M A T ! 2 F 6 . 2 ) ' • -1 1 0 F O R M A T ( ' F R E Q R A T I O ' ) 1 2 F O R M A T ! 6 F 8 . 4 ) 1 2 1 F O R M A T ( 4 F 8 . 4 ) 1 4 F O R M A T ( 1 4 , 4 F 9 . 3 ) 1 5 F 0 R M A T ( T 8 ( 1 X , F 5 . 2 ) / ) 1 6 F 0 R M A T ( 3 X , I 2 , 1 X , 1 0 F 1 1 . 4 ) 1 7 F O R M A T ( / / 2 X , 5 H L A Y E R 8 X , 9 H U - C 0 M P L E X 1 3X , 9 H W - C 0 M P l. E X 9 X , 5 H U = A M P 6 X , 1 5 H U - P H S 6 X , 5 H W - A M P 6 X , 5 H W » P H S 6 X , 5 H P = A N G 6 X . , 5 H S - A N G / 1 2 X , 4 1 - I R E A L 7 X , 2 4 H I M A G 7 X , 4 H R E A L 7 X , 4 H I M A G / ) 1 8 F O R M A T ( / / 3 5 X , 1 0 H M A T R I X F O R F 6 . 2 , 1 7 H C Y C L E S P E R S E C O N D / / ) 1 9 F O R M A T ( / l 5 H D E T E R M I M A N T = 2 E 1 7 . 7 / / / ) 2 9 F O R M A T ! I X , 5 H L A Y E R 2 X , 5 H P - » V E L 4 X , 5 H S « V 6 L 3 X , 7 H D E N S I T Y 3 X , 5 H 0 E P T H / ) E N D 1 NON-NORMAL INCIDENCE SYNTHETIC SEISMOGRAM FOR ANGLES LESS THAN CRITICAL * * * * * * * * * ***** * * * * * * * * * * * ******** * * * * * * * ****** * * * * * * * * * * * * * * * ******* * * P R O G R A M D E S C R I P T I O N • . r HE P R O G R A M O S E S T H E T I M E D O M A I N L I N E A R S Y S T E M S C R U S T M O D E L F O R S Y N T H E S I S O F S A M P L E D T I M E S E R I E S S Y N T H E T I C S E I S M O G R A M S P R O G R A M L I M I T A T I O N S MA X I MUM N O . O F L A Y E R S = 2 0 M A X I M U M N O . O F H E L D P O I N T S P E R L A Y E R = 2 0 0 M A X I M U M N O . O F A T T E N U A T O R O P E R A T O R W E I G H T S P E R L A Y E R = 2 0 r HE A T T E N U A T O R S M O S T lit C A U S A L M A X I M U M N O . OF. O U T P U T T I M E P 0 I N T S = 1 0 0 ( ) I N P O T C A R D S : C A R D 1 : T I T L E C A R D 2 : F L A G F O R A T T E N U A T I O N M O D E L : F O R NO A T T E N U A T I O N E N T E R 0 0 ' C A R D 3 : N O . O F L A Y E R S , N O . O F T I M E S A M P L E P O I N T S , T I M E I N C R E M E N T P I N C I D E N C E A N G L E , S I N C I D E N C E A N G L E . T H E . N O N Z E R O A N G L E D E T E R M I N E S T H E F O R M O F T H E I N P U T W A V E F O R M . C A R D 4 F F : L A Y E R T H I C K N E S S , P V E L O C I T Y , S V E L O C I T Y , D E N S I T Y C O N T I N U E T O B A S E M E N T T H E F O L L O W I N G C A R D S A R E I N C L U D E D O N L Y I F N A T T E N I S N O T E Q U A L T O 0 0 C A R D 5 F F : P A N D S A T T E N U A T O R L E N G T H S F O R L A Y E R S 1 T O N C A R D 6 F F : P A T T E N W E I G H T S A N D S A T T E N W E I G H T S F O R L A Y E R S 1 TO N T H E I N P U T W A V E F O R M T I M E S E R I E S S A M P L E P O I N T S O F A N Y L E N G T H A R E R E A D F R O M D E V I C E 4 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 0 I M E N S I O N E P O L A Y ( 2 0 ) , E S D I A Y ( 2 0 ) , F P D L A Y ( 2 0 ) , F S 0 L A Y ( 2 0 ) , E P ( 2 0 , 2 0 0 ) , 4 E S < 2 0 , 2 0 0 ) , F P ( 2 0 , 2 0 0 ) , F S ( 2 0 , 2 0 0 ) , R P P ( 4 1 ) , R P S ( 4 1 ) , R S P ( 4 1 ) , R S S ( 4 1 ) , $ T P P ( 4 1 ) , T P S ( 4 1 ) , T S P ( 4 l ) , T S S ( 4 1 ) , N T A O P ( 2 0 ) , N T A 0 S ( 2 0 ) , S T A T ( 2 0 ) , A L P ( 2 $ 1) , B E T ( 2 1 ) , R H O ( 2 1 ) , H ( 2 1 ) , A I N ( 1 0 0 0 ) , P O U T ( 1 0 U 0 ) , S O O T ( 1 0 0 0 ) , U S U R F ( 1 0 0 S O ) , W S O R F ( 1 0 0 0 ) , E ( 2 1 ) , F ( 2 1 ) , A T T E N P ( 2 0 , 2 0 ) , A T T E N S ( 2 0 , 2 0 ) , J ' P ( 2 0 ) , J S ( Z S O ) D A T A P I / 3 . 1 4 1 5 9 2 7 / I N I T I A L I Z E A L L O U T P U T S 0 0 1 0 1 = 1 , 2 0 E P 0 L A Y ( I ) = 0 . F S D L A Y ( I ) = 0 . F S O L A Y t I ) = 0 . FPOLA Y ( I ) =0.' N T A U P ( I ) = 0 . N T A U S ( I ) = 0 . OU 10 J=1,?00 E P ( I , J ) = (). E S ( I , J ) = 0 . F P ( I , J ) = (). F S ( I , J ) = 0 . ;-lr'1IN = 0 READ!5,1) STAT W RIT E ( 6 ,1 ) S TA T REA D ( 5 , 9 ) MATTEN REAQ(5,2) i M , r W A X t ' r i i M C , ANG INP , AN"G INS i\jP = N + l i\J|-i= |\]L-I1 R E A 0 ( 5 , 3 ) ( H ( I ) , A L P ( I ) , B E T ( I ) ,RHO(I ) ,1 = 1 ,N ) R E A D ( 5 , 3 ) DUMMY,ALP(NP)»8ET(NP),RHQ(NP) W R I T E ( 6 , 4 ) W R I T E ( 6 , '5 ) ( I , A L P ( I ) , 3 ET ( I ) , RHO (I ) , H ( I ) , [ = 1, N ) WRITE(6,5) MP ,ALP(MP),BET(NP) ,RHO(NP) ANGINP=AM GIM P * P I / 1 8 0 . A N G IM S = A M G INS * P I / 1 80 . CALCULATE TRANSMISSION AND REFLECTION COEFFICIENTS . CALL CUEGEM I N , A L P,B E T , R H O , A N G I N P , A N G I N S , E , F , T P P , T P S , R P P , R P S , T S P , STSS ,RSP,RSS) CALCULATE TRANSIT TIMES AND DELAYS 0 0 50 1 = 1 ,N ... T A O P = H ( I ) * C O S ( E l l ) ) / A L P ( I ) TAOS = H (T ) *COS ( F ( I ) ) / 3 ET ( I ) NTAOPI I) = I F I X(TAUP/TINC + 0.5) NTAUS( I ) = I F I X(TAUS/TINC + 0.5) I F ( M T A U P ( I ) .EO.O ) NTAUP(I)=1 I F ( N T A U S ( I ) .EQ . 0 ) N T A U S ( I ) = 1 CONTINUE IF (NATTEN.EQ.O) GU TO 21 RE A 0 IIM MODEL ATTENUATION PARAMETERS READ IN ATTENUATION OPERATOR LENGTHS R E A 0 ( 5 , 9 ) ( J P ( I ) , J S ( I ),I=1,N) R E « 0 A TIE NO A TOR '..'EIGHTS WRITE(6,101) 00 21 1=1 ,N J P L I M = J P ( I ) .) SL I ;•'= JS ( I ) REAU(5,102) ( AT T E N P U ,K) ,K=1 , J P L I M ) R EA 0 ( 5 , 1 0 2 ) (ATT ENS(I,<),<=1,JSLIM) W R I T E R , 1 0 3 ) (ATTENPI I ,K) ,K=1,JPLIM) WRITE(6,103) (ATTEMS(I ,<) ,< = 1 , J S L I 1 ) 2 1 CONTINUE W R I T E ( 6 , a ) C REAO INPUT WAVEFORM FROM OEVICE 4 R EA 0 ( 4 , 7 , ENO =99 ) ( A IN ( I ) , I = 1, 1000 ) 9 9 M5N0=I 00 98 I = MENO ,MMAX A I N ( I ) = 0 .0 98 CONTINUE 00 100 H=1,MMAX 1 F ( ANG IMP .NE .0.0) GO TO 22 SIN PUT= A IN(M) PIN PUT = 0 .0 GO TO 2 3 22 PINPUT = AIN(M ) SIMPUT=().() 2 3 CUNT INOF C UPDATE ALL DELAY OUTPUTS C FOR ALL LAYERS IF(WATTEN.EO.O) GO TO 24 C ATTENUATING MODEL SYNTHESIS DELAY OPERATIONS FO* ALL LAYERS C DELAY 00 28 1 = 1 ,N NT P=NTAUP(I) NTS = NTAUS(I) MTPM = MTP«1 NTSM=NTS«1 EP INTP=EP(I ,NTP) F PIMTP = F P ( I , NT P ) ES INTS=ES(I ,NTS ) F S I NTS = FS( I , NTS ) C ATTENUATION OPERATOR CONVOLUTIONS F. P O L A Y ( I ) = E P I N T P * A T T E N P ( I » 1 ) E S U L A Y ( I ) = E S I N T S * A T T E N S ( 1 , 1 ) F P O L A Y U ) = F P I N T P * A T T E N P ( I , 1 ) F S O L A Y ( I ) = F S I N T S * A T T E N S ( I , 1 ) i )D 2 5 K = l , N T P M I F ( K . G E . J P ( I ) ) G O T O 2 6 E P ( I ,MTP"K + 1 ) = E P ( I , N T P ™ < ) + E P I N T P * A T T E N P ( [ , K + 1 ) F P ( I , N T P ™ K + 1 ) = F P ( I , N T P « K ) + F P I N T P * A f T E N P ( I , K + 1 ) GO TO 2 5 2 6 E P ( I , N T P « K + 1 ) = E P ( I , N T P « K ) F P ( I ,I\ITP-K + 1 ) = F P ( I , N T P " K ) 2 5 C O N T I N U E OO 2 8 K = 1 , N T S N I F ( K . G E . J S ( I ) ) G O T O 2 7 ES ( I , N T S - K + l ) = E S ( I , N T S » ; < ) + E S I N T S * A T T E N S ( I , K + 1 ) F S ( I , N T S - K + 1 ) = F S ( I , N T S " K ) + F S I N T S * A f " T E N S ( I , K + 1 ) GO TO 2 8 . 2 7 E S ( I , N T S « K + 1 ) = E S ( I , N T S - K ) FS ( I »N T S " K +1 ) = FS ( I , N T S r j K ) 2 8 C U N T I N O E GO TO 3 0 C M O N C 3ATT E M U A T I N G M O D E L S Y N T H E S I S D E L A Y O P E R A T I O N S • 2 4 C O N T I N U E C D E L A Y 0 0 3 0 1 = 1 , M N T P = N T A U P ( I ) N T P M = N T P « 1 N T S = N T A O S ( I ) N T S M = N T S « 1 -E P D L A Y ( I ) = E P ( I , IMTP ) E S O L A Y ( I )= E S ( I , N T S ) F S O L A Y ( I ) = F S ( I , N T S ) F P O L A Y I I ) = F P ( I , N T P ) 0 0 3 0 1 K = 1 , N T P M E P ( I , M T P » K + 1 ) = E P ( I , M T P « ' < ) 3 0 1 F P ( I , N T P , j K + 1) = F P ( I , N T P , ! , K ) 0 0 3 0 2 K = l , N T S M E S ( I , N T S - K + 1 ) = E S ( I , N T S « K ) 3 0 2 F S ( I , M T S - K +1 )= F S ( I , N T S - K ) CONTINUE 3 0I-1PUTE WAVE FORMS RETRANSMITTED INTO BASEMENT AT TIME SAMPLE <A POUTPT=TSP(M+NP)*FSDLAY(N)+TPP(N+NP)*FPDLAY( N)+RS° ( NP ) *SI NPUT + RPP S( NP ) *P INPUT SOUTPT=TSS(N+NP)*FSDLAY(N)+TPS(N+NP)*FPDLAY(N)+RSS(NP)*SINPUT+RPS $(NP)*PINPUT JUNCTION UPDATE TO TIME SAMPLE M BOTTOM LAYER E P ( M , 1 ) = TS P ( N P ) *S IN PUT +T PP ( NP ) *P IN PU T + RP P ( N+N P ) *FP DL AY ( N ) + RS P ( N + N P $)*FSULAY(N) ES(N,1)=TSS(NP )*SINPUT+TPS(NP)*P INPUT + RSS(N+NP)*FSDLAY(N)+RPS(N+NP $)*FPDLAY(N) FP(N,1)=TS P(N+NP»1)*FS DLAY(N"1)+TPP(N +NP"1)*F P D L A Y ( N - l ) + R S P ( N ) * E S D $ LA Y ( N ) + R P P ( N ) * E P D LA Y(N ) FS ( M , 1 ) = TS S ( N +N P» 1 ) * FS DL A Y ( N™ 1) +T P S ( N +N P •• 1 ) *F P DL A Y ( N-1 ) +R S S ( N ) * E S D $LAY(N)+RPS(N)*EPDL4Y(M) INTERMEDIATE LAYERS IF(N.EO.2 ) GO TO 41 DO 40 I=2,NM EP( 1,1 ) = TSP( I + L ) * E S D L A Y ( I + 1 ) + T P P ( I + 1 ) * E P D L A Y ( [ + 1 ) + R S P ( I+NP)*FSDL*Y $ ( I ) +RPP ( I +:•! P) *F PD.LA Y ( I ) . -ES( I ,1 ) = TSS( I + L )*ES OLA Y ( I + l ) + T P S ( I + l ) * E P D u A Y ( I + 1 ) + R S S ( I + NP)*FSDLAY S( I ) +RPS( I + N P ) * F P 0 L A Y ( I ) r P ( I , 1 ) = TSP( I+NP-1)*FSOLAY(I»l)+TPP(I+NP»1)#FPDLAY(I-l)+RSP( I ) * E S O $ L A Y ( I ) + R P P ( I ) * E P 0 L A Y ( I ) FS( I,1 ) = TSS(I+MP-1)*FS DLAY(I-l)+TPS(I+NP-l)*FPDLAY(I»1)+RSS(I)*ESD $ L A Y ( I J + R P S ( I ) * E P D L A Y ( I ) CONTINUE TOP LAYER E P ( 1 , 1 ) = T S P ( 2 ) * E S D L A Y ( 2 ) + T P P ( 2 ) # E P D L A Y ( 2 ) + R S P ( 1 + N P ) * F S D L A Y ( 1 ) + R P P ( S l + NP) *F P D L A Y ( 1 ) ESC 1 * 1 ) = T S S ( 2 ) * E S D L A Y ( 2 ) + T P S ( 2 ) * E P D L A Y ( 2 ) + R S S ( 1 + N P ) * F S D L A Y ( 1 ) + R P S ( $ 1 + NP) *FPDLAY( 1 ) F P ( 1 , 1 } = R S P { 1 ) * E S D L A Y ( 1 ) + R P P ( 1 ) # E P D L A Y ( 1 ) F S ( 1 , 1 ) = R S S ( 1 ) * E S D L A Y ( 1 ) + R P S ( 1)*EPDLAY(1) END OF JUNCTION UPDATE TO TIME SAMPLE M COMPOTE SURFACE MOTIONS FROM IMPINGING S AND P WAVEFORMS PZ=( E P U L A Y U )* ( 1 ,+RPP ( 1 ) ) +ESOL AY ( 1 ) *RSP( 1) ) *COS( E ( 1 ) ) PX = ( E POLAY ( 1) *C1..«RPP ( 1 ) )»ESOLAY ( 1 )*RSP ( 1 ) ) # SI N ( E ( 1 ) ) S Z = ! ES O L A Y ( 1 ) * (1,+RSS(1 ) ) +E P O L AY ( 1 ) * R P S ( 1 ) )*St M( F ( 1 ) . ) S X = { " E S O L A Y ( 1 ) * ( ! . » R S S ( 1 ) ) + E P O L A Y { 1 ) * R P S ( 1 ) ) * C O S ( F ( 1 ) ) 0 0 = P X + S X W O = P Z + S Z X M = F L O A T ( M » I ) * r I N : C W R I T E O U T F I R S T 2 0 0 T I M E S A M P L E S F O R C H E C K I F ( M . L T . 2 0 0 ) WR I T E ( 6 , 6 ) X M , P I N P U T , S I N P U T , P U U T 3 T , S O U T P T , W O , U O C S A V E O U T P U T S U R F A C E M O T I O N S A N O B A S E M E N T W A V E F O R M T I M E S E R I E S F O R C L A T E R P L O T T I N G O R W R I T E O O T P O U T ( M) = P O I J T P T S O U T ( M ) = S O U T PT U S U R F ( M ) = U ( ) W S U R F ( M ) = W O 1 0 0 C O N T I N U E S T O P 1 F O R M A T ! 2 0 A 4 ) 2 F O R M A T ( 1 2 , 1 4 , 3 F 8 . 4 ) 3 F O R M A T ! 4 F B . 4 ) 4 F O R M A T {» L A Y E R " P « V E L S - V E L D E N S I T Y D E P T H ' ) 5 F O R M A T ! 2 X , 1 2 , I X , 4 F 1 0 . 4 ) 6 F 0 R M A T ( F 8 . 4 , 6 E 1 7 . 7 ) 7 F O R M A T ! 1 0 F 8 . 5 ) 8 F O R M A T ( ' T I M E I N P U T P P U L S E I N P U T S P U L S E 8 A S E -1E N f P R E E L S B A S E M E N T S R E F L V E R T S U R F M O T I O N H O R Z S U R F M O T I O N ' ) 9 F O R M A T ( 2 0 1 2 ) . . ' 1 0 1 F O R M A T ) ' F I L T E R W E I G H T S F O R P A N D S M O T I O N S I N L A Y E R S 1 T O N 1 ) 1 ) 2 F O R M A T ( 1 6 F 5 . 3 ) 1 0 3 F O R M A T ! 1 6 F 6 . 3 ) . 1 ) 4 F O R M A T ( 1 0 E 1 0 . 3 ) E N D Z O E P P R I T Z C 0 6 F F I C I E N f G E N E R A T I O N P R O G R A M * * * * * * * * * * * * * * * * i'f * * * * * * * * * * * * * * * * * * sis * * * . S 0 8 R O O T I N E 0 I S C R.I P T I O N V A R I A B L E S T R A N S M I T T E D F R O M M A I N M : N U M B E R O F L A Y E R S A L P , B E T , R H O : P A N O S V E L O C I T I E S A N O L A Y E R D E N S I T I E S A N G I N p , A N G IN S : T H E I N C I D E N T P A N D S A N G L E S A T T H E B A S E M E N T B O O N O A R / E , F : T H E I N T E R N A L L A Y E R P A N D S A N G L E S W I T H T H E N O R M A L T P P , R P P E T C : Z O E P P R I T Z ' C O E F F I C I E N T S F O R T O P - B O U N D A <Y 1 TO N + 1 F O R B O T T O M B O U N D A R Y M+2 T O 2 N S U B R O U T I N E C 0 E G E N ( N , A L P , 8 E T , R H 0 , A N G I N P t A N G I N S , E , F , T P P , T P S , R P P , 1 R P S , T S P , T S S , R S P , R S S ) M C C A M Y N O T A T I O N F O R Z O E P P R I T Z ' E Q U A T I O N S 0 I M E M S I O N A L P I 1 ) , 8 E T ( 1 ) , R H O ( 1 ) , E ( 1 ) , F ( 1 ) , T P P ( 1 ) , T P S ( 1 ) , R P P ( 1 ) , R P S ( £ 1 ) , T S P ( 1) , T S S ( 1 ) , R S P ( 1 ) , R S S ( 1 ) , W ( 2 , 2 ) , Z ( 2 , 3 ) , C O F ( 2 , 1 ) N P = M + 1 I F ( A N G I M P . E O . O . ) G O T O 2 E ( N P ) = A N G I M P C = A L P ( N P ) / S I N ( E ( N P ) ) G O T O 3 F ( N P ) = A N G I N S C = 8 E T ( N P ) / S I N ( F ( N P ) ) C O N T I N U E 0 0 4 I = 1 , N P X = A L P U ) / C Y = B E T ( I ) / C I F I X . G T . l . ) G O T O 1 0 E ( I ) = A " f A N 2 ( X , S Q R T ( 1 . » X * X ) ) I F { Y . G T . 1 . ) G O T O 1 5 F ( I ) = A T A N 2 ( Y , S O R T ( l - . - Y * Y ) ) GO T O 4 W R I T E ( 6 , 1 1 ) I E ( I ) = 1 . 5 7 0 7 9 6 4 G O T O 1 4 W R I T E ( 6 , 1 1 ) I F ( I ) = 1 . 5 7 0 7 9 6 4 GO T O 4 11 F O R M A T ( 4 7 H C R I T I C A L R E F L E C T I O N ™ N O T R A N S M I S S I O N T O L A Y E R 1 2 ) 4 • C O N T I N O E C F R E E S U R F A C E R E F L E C T I O N C O E F F I C I E N T S T P P ( 1 ) = 0 . T P S ( 1) = 0 . T S P ( 1 ) = 0 . T S S ( 1 ) = 0 . W( 1 , 1 ) = « S I N ( 2 . * E ( 1 ) ) W ( 2 , 1) =+C O S ( 2 . * F ( L ) ) W ( 1 , 2 ) = + A L P ( 1 ) * C 0 S ( . 2 . * F ( 1 ) ) / B E T ( 1 ) W( 2 , 2 ) = + B E T ( 1 ) * S IN (2 - * F (1 ) ) / A I . P ( 1 ) 0 0 2 0 L = l , 2 0 0 2 0 M = l , 2 2 0 Z ( L , H ) = W ( L , M ) Z ( 1 , 3 ) = « W ( 1 , 1 ) . Z ( 2 , 3 ) = + W ( 2 , 1 ) C A L L G A O E L I ( 2 , Z , C O F , O E T , 2 , 3 ) R P P ( 1 ) = C O F ( 1 , 1 ) R P S ( 1) = C 0 F ( 2 ,1 ) 0 0 2 1 L = l , 2 ' • 0 0 2 1 M = l , 2 , 2 1 Z ( L , M ) = W (L , M ) Z ( 1 , 3 ) = - W ( 1 , 2 ) Z ( 2 , 3 ) = W (2 , 2 ) C A L L G A O E L I ( 2 , Z , C O F , D E T , 2 , 3 ) R S P ( 1 ) = C O F ( 1 , 1 ) R S S ( 1) = C O F (2 ,1 ) 0 0 1 2 I = 2 , N P J = I" 1 K= I 1 2 C A L L C O E F F ( A L P , 6 E T , R H 0 , E , F , T P P , T P S , R P P , R P S , T S P , T S S , R S P , R S S , I , J , K ) 0 0 1 3 1 = 1 , N J = I + l K N P = I + N P 1 3 C A L L C O E F F ( A L P , B E T , R H O , E , F , T P P , T P S , R P P , R P S , T S P , T S S , R S P , R S S , I , J , K N P 1 ) R E T U R N • E M U S E C O N D C O E F F I C I E N T S U B R O U T I N E ****** ******************* T H I S S U B R O U T I N E I S R E Q U I R E D B Y - C O E G E N * * : ; : * * * * * * * * * * =1: * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * S U B R O U T I N E C O E F F ( A L P • B E T , R H O , E , F , T P P , T P S , R P P , R P S , T $ P , T S S , R S P , R S S , I 1 , J , K ) 0 I i-i E N S I O N A L P ( I ) , B E T ( 1 ) , R H O ( 1 ) , E ( h) , F ( 1 ) , T P P ( 1 ) , T P S ( 1 ) , R P P ( 1 ) , R P S ( $ 1 ) , T S P ( 1) , T S S ( 1 ) , R S P ( 1 ) , R S S ( 1 ) , C ( 4 , 4 ) , Z ( 4 , 5 ) , C O F ( 4 , 1 ) C ( l , l ) = ™ S I N ( E ( [ ) ) ' C ( 1 , 2 ) = + C O S ( F ( I ) ) C ( 1 , 3 ) = - S I N ( E ( J ) ) C ( 1 , 4 ) = + C O S ( F ( J ) ) C ( 2 , l ) = + C O S ( E ( I ) ) C { 2 , 2 ) = + S I N ( F ( I ) ) C ( 2 , 3 ) = - C O S ( E ( J ) ) C ( 2 , 4 ) = « S I N ( F ( J ) ) C ( 3 , l ) = + C O S ( 2 . * F ( [ ) ) C ( 3 , 2 ) = + B E T ( I ) / A L P ( I ) * S I N ( 2 . * F ( I ) ) C ( 3 , 3 ) = + R H O ( J ) / R H O ( I ) * A L P ( J ) / A L P ( I ) *CO S { 2 . * F ( J ) ) C ( 3 , 4 ) = + R H O ( J ) / R H O { I ) * 3 E T ( J ) / A L P U ) * S I N C 2 . * F ( J ) ) C ( 4 , 1 ) = " S I N ( 2 . * E ( I ) ) C ( 4 , 2 ) = + A L P ( I ) / B E T ( I ) * C O S ( 2 . * F ( I ) ) C ( 4 , 3 ) = R H O ( J ) / R H O ( I ) * ( ( 8 E T ( J ) / B E T ( I ) ) * * 2 ) * A L P ( I ) / A L P l J ) * S I N ( 2 . * E ( J 1 ) ) C ( 4 , 4 ) = - R H 0 ( J J / R H J ( I ) * ( ( B E T ( J ) / B E T ( I ) ) * * 2 ) * A L P ( I ) / B E T ( J ) * C O S ( 2 . * F ( U ) ) D O 2 0 L = l , 4 D O 2 0 0 = 1 , 4 Z ( L , i ' l ) =C ( L , (A ) Z( 1 , 5 ) = " S I N ( E ( I ) ) Z ( 2 , 5 ) = « C O S ( E ( I ) ) Z ( 3 , 5 ) = + C O S ( 2 . * F ( I ) ) Z ( 4 , 5 ) = + S I N ( 2 . * E ( I ) ) C A L L G A O E L 1 ( 4 , 7 . , C O F ,OET , 4 , 5 ) R P P ( K ) = C O F ( l , l ) • R P S ( K ) = C 0 F ( 2 , 1 ) T P P ( :<) = C 0 F ( 3 ,1 ) T P S ( K ) = C ' O F ( 4 , 1 ) 0 0 2 1 L = l , 4 0 0 2 1 .14 = 1 , 4 21 Z ( L t\A) =C ( L ,M ) Z( 1 , 5 ) = + C 0 S ( F ( I ) ) . Z ( 2 , 5 ) = " S I N ( F ( ( ) ) Z ( 3 , 5 ) = + B F T ( I ) / A L P ( I ) * S I N ( 2 . * F ( I ) ) Z ( 4 , 5 ) = « A I. P ( I ) / B E T ( I ) * C O S ( 2 . * F ( I ) ) C A L L G A O E L I ( 4 , Z , C O F , O E T , 4 , 5 ) R S P ( K ) = C O F { I ,1 ) K S S ( K ) = C 0 F ( 2 , 1 ) T S P ( <) = C O F ( 3 ,1 ) T S S ( K ) = C U F ( 4 , 1 ) K E T O K M E N 0 C B A S I C L A Y E R T R A N S F E R F U N C T I O N G E N E R A T I O N P R O G R A M C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # * * * * * * * * * * * C T H I S S U B R O U T I N E I S C A L L E D B Y T H E F O U R I E R C O E F F I C I E N T M A I N P R O G R A M C # * * * * * * * * *:;= * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * S UB R O U T I N E T R A N S I ( R P P I , R P S I , R S P I , R S S I , R P P J , R P S J , R S P J , * S S J , G 1 , G 2 , H I 1 , H 2 , H 3 , H 4 ) C O M P L E X R P P I , R P S I , R S P I , R S S I , R P P J , R P S J , R S P J , R S S J , G l , G 2 , H I , H 2 , 1 H 3 . , H 4 , 8 1 , B 2 , B 3 , B 4 , B 5 , B 6 , B 7 , B 8 , D E N B 1 = G 1 * R P P J B 3 = I U * R S P I 3 1 = B 1 * R P P I 6 6 = G 1 * R P S J B 8 = B 6 * R S P I B 6 = B 6 * R P P I B 5 = G 2 * R S S J B 7 = B 5 * R S S I B 5 = B 5 * R P S I B 4 = G 2 * R S P J •3 2 = 8 4 * R P S I B 4 = B 4 * R S S I B 1 = B 1 + .B2 B 3 = B 3 + B 4 B 5 = B 5 + B 6 B 7 = B 7 + B 8 • B 2 = G 1 / ( 1 . - G l * S l ) B 4 = G 2 / ( 1 . » G 2 * . B 7 ) O E N = l . / ( l . ~ 8 2 * B 3 * B 4 * B 5 ) H1 = B 2 * U E N H 2 = B 4 * H 1 H3 = B 5 * H 2 H 2 = B 3 * H 2 H4 = 8 4 * U E N R E T U R N . • ' C GAUSS E L I M I N A T I O N SUBROUTINE FOR REAL S I M U L T A N E O U S E O O A T I O N S r * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * « * * * * * * * * * * ;|: * * * * * * s|: * * :1c * * * * * * C PROGRAM 01 S C R I P T ION C T R A N S M I T T E D V A R I A B L E S . C N: NUMBER OF E O O A T I O N S C A: I N P U T C O E F F I C I E N T S ( CONSTANT VECTOR [ N ROW N+1) C X: OUTPUT S O L U T I O N VECTOR C 0 ET : C HEC f< 0 ET E R -I INANT C 10: M A I N PROGRAM D I M E N S I O N OF F I R S T D I M E N S I O N OF A C J O : M A I N PROGRAM D I M E N S I O N OF SECOND D I M E N S I O N OF A C * * * * * * * * * -f * * * * * * * * * * * -:- * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * S U B R O U T I N E G A U E L I ( N , A , X ,OET , 10 , J O ) D I M E N S I O N A ( I D , J O ) , X ( I D , 1 ) NP=N+L NM=N«1 O E T = 1 . 0 C T R I A N G U L A R I Z A T I O N DO L2 K=1,NM KP=K+1 R = 1 . 0 / A K , < ) U E T = O E T * A ( K ,K ) DO 11 J = <P ,NP 11 A( K , J )=R*A ( K , J ) 0 0 12 I = K P , N S = A ( I , K ) . DO 12 J = K P , N P 1 2 A( I , J ) =A ( I , J ) ~S*A (K , J ) •DET = O E T * A ( N , N ) C BACK S U B S T I T U T I O N X ( N , I ) = A ( N , N P ) / A ( N , N ) DO 2 2 I= 1,NM K = N - I KP=K+1 S = A ( K , N P ) DO 21 J=KP,N 21 S = S » A ( K , J ) # X(J , 1 ) 2 2 X ( K , 1 ) = S RETURN 'END , G A U S S E L I M I N A T I O N S U B R O U T I N E F O R C O M P L E X S [ M U L T A N E O U S E O U A T I O N S * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * P R O G R A M I N S C R I P T I O N T R A N S M I T T E D V A R I A B L E S N : N U M B E R O F £ 0 0 A T I O N S A : I N P U T C O E F F I C I E N T S ( C O N S T A N T V E C T O R IN ROW N + 1 ) X : D O T P O T S O L U T I O N V E C T O R D E T : C H E C K D E T E R M I N A N T I D : M A I N P R O G R A M D I M E N S I O N O F F I R S T D I M E N S I O N O F A J O : M A I N P R O G R A M D I M E N S I O N O F S E C O N D D I M E N S I O N O F A * * * * * * :|c * * * * * * * * * ;\: * * * * * # « * * * * * * * * * * * * * * * * * * * * * * * s|s * * :;: * * S O B R O U T I N E G A O C O M ( N , A , X , D E T , I D , J O ) G A O S S E L I M I N A T I O N S U B R O U T I N E F O R C O M P L E X S I M U L T A N E O U S E Q U A T I O N S C O M P L E X A , X O E T , S 0 H E N S I O N A ( ID , J D ) , X ( 1 0 ,1 ) N P = M + 1 NM=N"1. 0 E T = 1 . 0 T R I A N G U L A R I Z A T I O N DO 1 2 K = 1 , N M K P = K + 1 R = 1 . 0 / A ( K , K ) D E T = D E T * A ( : < , K ) 0 0 11 J = K P , N P " ( A ( K , J ) = R * A ( K , J ) DO 1 2 I = K P , N S = A( I , K ) 0 0 1 2 J = K P , N P A ( I , J ) = A ( I , J ) - S * A ( K , J ) ' D E T = 0 E T A ( N , N ) B A C K S U B S T I T U T I O N X ( N , 1 ) = A ( N , N P ) / A ( N , N ) D O 2 2 1 = 1 , M H K = N " I K P = K + 1 S = A ( K , I M P ) 0 0 2 1 J = K P , N S = S " A ( K , J ) * X l J , l ) X ( K , 1 ) = S R E T U R N F MO 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items