A NUMERICAL INVESTIGATION OF EXPLOSIVE GENERATED NORMAL MODES AND LEAKING MODES IN AN UNSATURATED SURFACE LAYER OVERLYING A SATURATED HALF-SMGE Dissertation for the Degree of Ph. D. MICHIGAN. STATE UNIVERSITY ‘ CHING-IIAN MO 1974 This is to certify thgt the ‘1’ thesis entitled A T‘Uf-EEVRICAL TNVES‘TIGATION OF EXPLOEIVE GENERATE‘D NORI-Q‘tL MCDFS AND LEIAKING MCDES IN AN UNSATURATEL‘ SURFACE LAYER OVERLYING A SATURATFD HALF—SPACE presented by Chino-nan Kao has been accepted towards fulfillment of the requirements for Ph. D. Geology degree in " MHZ 2&M /Major professor i Date 2/25/74 0-7639 r' VL‘I :fIL'HSIZICAL 1-“. .- I r ‘ u-s—s Asa LEAH ”rel m. "ZRLYITIG A SAN \ Most of ti seismegrams are farm. The rea: near surface 1' r-ar.hquake sur Pattern, there the Cr‘JSt'Uppe to use the “9' have data frOr ABSTRACT A NUMERICAL INVESTIGATION OF EXPLOSIVE GENERATED NORMAL MODES AND LEAKING MODES IN AN UNSATURATED SURFACE LAYER OVERLYING A SATURATED HALF-SPACE BY Ching—nan Kao Most of the surface waves recorded on exploration seismograms are characterized by the poorly develOped wave form. The reasons are due to short recording distance, near surface inelastic material, etc. On the other hand, earthquake surface waves have a more complete dispersion pattern, therefore, they have long been used to interpret the crust-upper mantle structure. This study is designed to use the high speed computer to process low quality surface wave data from small explosions, in the hope that the geological information contained in the surface waves of exploration seismograms may be better utilized. A new moving window spectral analysis method successfully identified superimposed modes which included M11, M21, M12, three "pipe organ" type modes, and others. It also produced a smooth time-varying spectrum, from which the group velocity dispersion was determined. Some methods of computing phase velocities were found to give poor results using available data. Bloch and Hale's (1968) method1 and a method.deve10ped herein, that directly integrates the observational group velocity igperSiOD C ;:3'.'iding 9C recorded at By ass: layer overly dispersion c find the mod were shear w half-space, rigidity rat is to solve 1 The n nOn-l Observationa After set/era fimse Velocit atraging tin h’iOther appn of more than Sum of Squarr similar resuj The She¢' mind to be I Ching-nan Kao dispersion curve were superior to other methods in providing good phase velocity dispersion curves for data recorded at source distances of 1000 feet and beyond. By assuming a theoretical model--a single surface layer overlying a infinite half-space, the obtained dispersion curve of the M11 mode is then inverted to find the model parameters. The model parameters determined were shear wave velocities in both surface layer and half-space, thickness of the surface layer, and the rigidity ratio of the two layers. One method of inversion is to solve n non-linear equations in the n unknowns. The n non-linear equations are formed by substituting observational phase velocities into the period equation. After several trials using different sets of observational phase velocities, the final solutions are obtained by averaging the results for each trial of n unknowns. Another approach is to solve for n unknowns for a set of more than n non-linear equations by minimizing the sum of squared residues. Both methods have yielded similar results. The shear wave velocity in the surface layer was found to be 603 feet/second. It is in complete agreement with the result found by direct measurement. The depth of the surface layer of 22 feet is less than the depth to the water table at which the compressional wave velocity increases from 1150 feet/second to 5700 feet/second. This is because normal modes are several times more sensitive 1'— czmpression was anticip saute veloci amputed. ' A theo; organ type the water ta completely i multiple of laree source superior to dEpth. 15- Bloch Geterninatic Seism. Soc Ching-nan Kao sensitive to the shear wave velocity profile than to compressional wave velocity structure. This result was anticipated. A density ratio of unity and a shear wave velocity for the half—space of 904 feet/second were computed. They were all consistent with known information. A theoretical interpretation of the dispersive pipe organ type leaking modes was attempted. The depth to the water table, computed by using these modes, agreed completely with the results obtained by using the first multiple of the P head wave from the water table. At large source distances, both methods were found to be superior to the standard refraction technique in computing depth. lS. Bloch and A. L. Hales, New techniques for the determination of surface wave phase velocities, Bull. Seism. Soc. Am., vol. 58, pp. 1021-1034, 1968. A NUMERICA VOPQiAL MOD .SL'RFACE L in Part1, A NUMERICAL INVESTIGATION OF EXPLOSIVE GENERATED NORMAL MODES AND LEAKING MODES IN AN UNSATURATED SURFACE LAYER OVERLYING A SATURATED HALF-SPACE BY Ching-nan Kao A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Geology 1974 The auth testers, Drs. Carmichael, R A. Vogel for It is al dgscussions w of Physics an Department of ACKNOWLEDGMENTS The author is extremely grateful to his guidance committee members, Drs. Hugh F. Bennett (Chairman), Robert S. Carmichael, Robert Ehrlich, Charles M. Spooner, and Thomas A. Vogel for their guidance and encouragement. It is also a pleasure to acknowledge several helpful discussions with Dr. Wayne W. Repko of the Department of Physics and Astronomy and Dr. James F. Hannan of the Department of Statistics and Probability. ii List of Table: istof Figurw tepter I. Introdt II. Theory A. Pha 3. Per 1. 2. III. Spectra A. Met D l. 2. 3. B Int! IV' MGthOds ObSEIt A Peal B. Four 1. 2. 3. C' Comp by Grc List of List of Chapter I. II. III. IV. TABLE OF CONTENTS Tables . . . . . . . . . . . . . . . . . . . Figures . . . . . . . . . . . . . . . . . . Introduction . . . . . . . . . . . . . . . . Theory . . . . . . . . . . . . . . . . . . . A. Phase Velocity and Group Velocity . . . . B. Period Equation . . . . . . . . . . . . . 1. Normal Mode Roots of Period Equation. 2. Leaking Mode Roots of Period Equation Spectral Analysis and Data Interpretation . . A. Methods Used to Identify Signals or Derive Group Velocities . . . . . . . . 1. Moving Window Spectral Analysis Method . . . . . . . . . . . . . . 2. Fixed Window Fourier Transform . . . 3. Determination of Particle Motion . . B. Interpretation of Data . . . . . . . . . Methods of Computing Phase Velocities From Observational Data , , , A. Peak-and-trough Method . . . . . . . . . B. Fourier Transform Methods . . . . . . . . l. Fourier Phase Difference Method . . . 2. Crosscorrelation Method . . . . . . . 3. Fourier Sum-and-difference Method . . C. Computing Phase Velocity DiSpersion Curve by Using the Knowledge of Observational Group Velocities. . . . . . . . . . . . iii Page vi ll 13 15 20 27 27 28 38 38 48 55 55 55 57 57 58 66 :a;ter V. Invei A.1 l 2 BA VI . S umma A. I» l 2. 3. 8. TI C. Mi AppendiXA APPEndix B Appendix c aeferences Chapter V. Inversion of Data to Yield Model Parameters . A. B. Inversion of Normal Mode Data 1. Methods 2. Checking the Correctness of Inverted. Model Parameters by Reversing Process A Discussion of the Model Parameters Derived From Normal Mode, Leaking Mode, and Compressional Wave Data VI. Summary and Conclusions . . . A. C. Appendix A Appendix B Appendix C References Numerical Techniques . . l. Derivation of Time-varying Spectra and Group Velocities 2. Computation of Phase Velocities . . . 3. Inversion Methods . . The Use of Normal and Leaking Modes in Data Interpretation . . Miscellaneous tions . . . Conclusions iv Recommenda- Page 71 71 71 81 81 9O 9O 92 94 95 99 100 103 ;&le 4a. 4b. Body W Udel List c List c Phase Obsc Curt Phase Obse Curx Model of I Meti Model of n of l (Rec A COmE Obse DIS} No. The DI 4b. LIST OF TABLES Body Wave Velocity Distribution in Udell Hills Area . . . . . . . . . . . List of Refraction Records . . . . . . . List of Events . . . . . . . . . . . . Phase Velocities Obtained by Integrating Observational Group Velocity Dispersion Curve (Record No. 43, Second 8 trace). Phase Velocities Obtained by Integrating Observational Group Velocity Dispersion Curve (Record No. 43, Sixth Z trace) . Model Parameters Obtained by the Inversion of M1 Data Through the Use of the Exact Metho (Record No. 43, Second 8 trace) Model Parameters Obtained by the Inversion of M11 Data Through the Use of the Method of Minimizing Sum of Squared Residues (Record No. 43, 3 Component) . . . . . A Comparison of the Theoretical and the Observational Phase and Group Velocity Dispersions (M11 data taken from Record No. 43, Second 8 trace) . . . . . . . The Depths to Water-table . . . . . . . Page 51 52 69 70 74 80 83 87 L L P‘“ C,‘ o Layeret A Port: Time-v; (Recc Events (Recc Events (Recc Time-Va 43. s Observa AmPlitu AmPlitu Determi Phase V An EXan TraPeZc A Small Varyi Phfie V FOuri NO. 4 ObServe Mode LIST OF FIGURES Figure Page 1. Layered Model . . . . . . . . . . . . . . . . . lO 2. A Portion of Record No. 43 . . . . . . . . . . 23 3. Time—varying Spectra of Multiple Reflections (Record No. 43, Second 5 trace) . . . . . . . 37 4. Events with Distinct Time-varying Spectrum (Record No. 43, Second 3 trace) . . . . . . . 39 5. Events with Distinct Time-varying Spectrum (Record No. 43, Second 3 trace) . . . . . . . 40 6. Time-varying Spectrum of M11 Mode (Record No. 43, Second 3 trace) . . . . . . . . . . . . . 41 7. Observational Group Velocities of M11 Mode . . 45 8. Amplitude Spectrum (Fixed Window) . . . . . . . 46 9. Amplitude Spectrum (Fixed Window) . . . . . . . 47 10. Determination of Particle Motion . . . . . . . 49 11. Phase Velocities of M11 Mode (Record No. 43) . 56 12. An Example of Shifting in Time Domain . . . . . 62 13. Trapezoidal Weighting Function . . . . . . . . 62 14. A Small Segment of Seismic Trace and Its Time Varying Filtered Version . . . . . . . . . . 65 15. Phase Velocity Dispersion Curve Obtained by Fourier Sum—and-Difference Method (Record No. 43, Second and Sixth 8 traces) . . . . . 67 16. Observed and Theoretical Dispersion of M11 Mode . . . . . . . . . . °.' . . . . . . . . 82 vi Normal a: srsurface geu is: use the bi airple geomet seismic inter. leaking modes simple theori Flier, the su tremendously “in? the th ties of Crust ldSt two deca The surf been widgly u usually recor fCl‘e' the SUI addition, the the material assme a cor 1 Chapter I Introduction Normal and leaking mode seismic surface waves contain subsurface geological information as do the body waves. Because the body waves, in most cases, are predictable by simple geometrical ray theory, they have long been used in seismic interpretation. On the other hand, normal modes, leaking modes, diffraction, etc. cannot be described by simple theories. Since the advance of the high-speed com- puter, the surface wave solution computing time has been tremendously reduced. The use of surface waves in deter- mining the thickness of the crust and the physical proper- ties of crust and upper mantle has become popular in the last two decades. The surface waves from a small explosion have not yet been widely utilized in seismic interpretation. They are usually recorded at short distances from the source, there— fore, the surface waves are usually not well developed. In addition, the near surface heterogeneity and departure of the material from perfect elasticity makes it difficult to assume a correct theoretical model. These facts inhibit the exploration seismologist from making practical use of the surface wave. Our study aims at improving the available techniques and develOping new methods in order to make use of the surface wave normal and leaking modes from a small explosion. The grot .19 mOtiOn’ a sxface wave- :an determine Two meti city diSPersj city by diviC ngtmm 0f E byinsyectiol This is a trc’ persion. Am sis techniqw (1969), but I wand. The: the one used The pha BMW differe iiiles (1968) a Smoother c ‘nlocities ti Cities. 2 The group velocity and phase velocity dispersion, parti- cle motion, and amplitude Spectrum are characteristic of a surface wave. They must be correctly computed before one can determine the correSponding earth model. Two methods have been used to calculate the group velo- city dispersion curve. One is to determine the group velo— city by dividing the distance to the station by the arriv- ing time of a particular period. Measurement of the period by inspection is valid only when a single mode is present. This is a traditional method of deriving group velocity dis- persion. Another method is the moving window Spectral analy— sis technique. It is similar to Landisman, gt alfS method (1969), but was developed without prior knowledge of his method. There are some differences between his method and the one used here as will be explained in Chapter III. The phase velocity dispersion can be determined using many different approaches. A method first used by Bloch and Hales (1968) obtained a good result. As an attempt to find a smoother curve, a method was developed to derive the phase velocities by direct integration of the observed group velo— cities. The results obtained were in agreement with the previous method of Bloch and Hales. Several other methods we had tried did not yield usable results. It was apparently due to the short time duration of signals and noises. The details of various methods will be discussed in Chapter IV. The refraction records used in this study were kindly Supplied by Dr. H. F. Bennett of the Department of Geology. gtotal of 24 The first 8 t: :ction. The 1 :ansverse mo vertical moti' rent geOphone glories used w 2f the record ties. The na and the dampi response was used were rec Settings. Ba A digiti tile the reco frequency of aLTJCy WEI-e We riertz). The IQCQ 3 A total of 24 traces were recorded on photographic paper. The first 8 traces were X component or horizontal in—line motion. The second 8 traces were Y component or horizontal transverse motion. The third 8 traces were 8 component or vertical motion. At each detection position, a three compo— nent geOphone recorded the X, Y, and Z motions. The geo- phones used were of the moving coil type. The amplitudes of the recordings were proportional to the particle veloci— ties. The natural frequency of the geOphones was 4.5 Hertz and the damping was 62% of the critical frequency. The response was essentially flat above 7 Hertz. All the data used were recorded without Automatic Gain Control and filter settings. Bandpass was therefore about 7 to 125 Hertz. A digitizer owned by the University was used to digi— tize the records for use on digital computer. The digitizing interval was 0.005 sec. which corresponded to the Nyquist frequency of 100 Hertz. The frequencies encountered in this study were well below the Nyquist frequency (i.e., 10—50 Hertz). The records were recorded in the Udell Hills area of Manistee County, Michigan. The ground surface in the re- cording area is essentially horizontal. The topographic correction, therefore, was not necessary. This area contains approximately 500 ft. of Pleistocene glacial drift. It is composed mainly of sand with strips of Clay. The underlying formation contains sandstone with some shale layers which is of Mississippian age. The water table at the locati at a depth Of description c Rid (1971) . In study :crreSponding chcsen before shear wave ve selecting a n attained by c' distance curt 55 a marked C 7-79 critical rated laYer table and the energy trapp‘ TiEr f e Ore: WI «Ecords take' are the tram his ‘nfiwing. 4 at the location where the records were obtained was located at a depth of about 32 ft. below the surface. A detailed description of the geology of this area can be found in Todd (1971). In studying surface waves, a particular period equation correSponding to a particular theoretical model must be chosen before the study proceeds. The compressional and shear wave velocity distribution in this area is vital in selecting a model. Table l is a list of body wave velocities obtained by direct measurements and from refraction time- distance curves. It is observed that the water-table acts as a marked discontinuity for compressional wave velocity. The critical angle of the compressional wave in the unsatu— rated layer is about ll%°. If all shots were above water- table and the shear waves were not taken into account, the energy trapped in the surface layer would be about 98%. Therefore, we expect the prominent feature of the seismic records taken in this area to be composed of signals that are the transient reSponse of the surface unsaturated wave guide. Through this study, a single surface layer model will be assumed. The period equation for this model can be found in Ewing, gt al.(l957, pp. 193). A problem exists in assuming the water-table to be the interface of our single layer model. The leaking mode is more sensitive to the compressional velocity distribution, whereas the normal mode is more sensitive to the shear velo— City distribution (Su and Dorman, 1965, p. 1018). In general, an increase in compressional velocity accompanies an increase fable 1. Bod CmPress ionaj Have Veloc 1 t: (ft/Sec) shear have Vela; ity ift/SEC) Obtained Table 1. Body Wave Distribution in Udell Hills Area Unsaturated Saturated Glacial Glacial Drift Drift Compressional Wave Velocity (ft/sec) #1150 *5700 Shear Wave Velocity (ft/sec) # 600 ? # Obtained by direct measurements * Obtained by time-distance curves Mississippian 8.8. *13461 in shear ve10c sure that the ;e'.ow the watu snear velocit; level. This Cebrin, g a tan, l960) h lie experimen increase belo table is smal "'EIOCitY. I f l or. “995 due t I) 33513“ a sin .ace, fmmd b dies ..erent fOr surface-lay€r ChaPter II I A . ‘6”an mode 6 in shear velocity. In this particular case, we can be sure that the compressional velocity increases greatly below the water-table, but we hesitate to infer that the shear velocity also increases substantially at the same level. This is because several papers (Biot, 1956a, 1956b; Dobrin, gt 31, 1954; Ewing, gt a1, 1957; Kisslinger, 1959; Mann, 1960) have reported, using either the theoretical or the experimental basis, that the shear velocity does not increase below the water-table. There is no solid ground to accept or to reject their finding, but it may be safe to assume that the change in shear velocity at the water— table is small as compared to the change in compressional velocity. If this assumption is true, the shear velocity changes due to the clay layers is possibly of greater magni- tude than the change at the water—table. Because we have assumed a single surface layer model, the depth to the inter- face, found by the inversion of the observed data, may be different for normal and leaking modes. However, the single surface-layer model will still be used through the study. In Chapter II, we will make a theoretical treatment of the leaking mode based upon the above assumption. Compressional wave (P-wave) velocities, shear wave (S-wave) velocities, thickness of layers, and rigidity ratio of any two adjacent layers are important parameters in characterizing a horizontal stratified model. Body wave velo- cities and thicknesses are possibly deduced by a simple geo- metrical ray theory. The rigidity ratios and their dependent ;a:aneters (1 (3;,- be foun: Two metl ;‘:.ase veloci‘ :etizods are it. Chapter V deveIOp a le. iata. Theyt 5:: of squar. 2’2:- true pha. also lineari 539? Solve . C3-TP‘Jter. O. is: least-sq 7 parameters (i.e., density ratio and Poisson's ratios) can only be found by using surface wave data. Two methods were developed to inverse the observed phase velocity dispersion data. The results from both methods are in good agreement. The details will be given in Chapter V. Dorman and Ewing (1962) were the first to develop a least-square technique to invert the surface wave data. They derived the normal equations by minimizing the sum of squared errors which were the differences between the true phase velocity and the assumed phase velocity. They also linearized the entire process. The methods in this study solve a system of non-linear equations directly on the computer. One solves for exact solutions and another solves for least-square solutions by minimizing the sum of the squared residues of the period equation. As before, the development of these methods were accomplished without know— ledge of the one explained by Dorman and Ewing (1962). The methods of this study will not be compared to theirs, but the results are expected to be in agreement. The best method to check the correctness of the inver- sion methods described above is to reverse the process. By substituting the model parameters obtained into the period equation of the assumed theoretical model, one can check the correctness of the inversion methods by comparing the re- sulting phase and group velocities with the observed data. The results of the surface wave computation not only provide an independent check on the information extracted from the bOdY cerning the Ph which are usua 2.3.15 in 1211856 the high-speed siierably lowe use of the SU] :ion, uncerta :emine local still has som cosy waves on 8 from the body wave data, but also gives more information con— cerning the physical properties of the subsurface layers which are usually not provided by body waves. The calcula- tions in these processes are laborious, but with the aid of the high-Speed computer, the time and cost have been con- siderably lowered. Despite the intrinsic drawbacks of the use of the surface wave data, such as the limited penetra- tion, uncertainty in determining layering, inability to de— termine localized structure, etc., the use of surface waves still has some advantages that cannot be surmounted by using body waves only. A sin<_ Figure 1. denote the rescective her i with thick iégth B an ‘4 l ,L’ (D (D (1.1 . Maiu‘ltim Chapter II Theory A single layer over a half—space model is shown in Figure 1. The symbols (it , g; , and fl in the figure denote the P wave velocity, S wave velocity, and density respectively, where the subscript 1 corresponds to the layer 1 . The transient reSponse of the surface layer with thickness ii for any component of the motion 11, at depth 8 and distance r from an impulse point source at depth b is given by a double integral (Ewing, et a1, 1957) such as: oo 00 . 9 w k ZLb J kr dk .............(?.1) u = ReyexPIJthwj i L LPT\-;3k)£ ) ' -00 --00 where Re denotes the real part of the double integral w = angular frequency k = wave number = w/C C = Phase velocity Jo = the Bessel function of zero order a function of w, k, B, and b D II P = a function of w and k 347-7. Evaluation of this double integral can be performed in several 9 10 Layer 0 0‘.» 30, F0 Z =0 Layer 1 04., 5', P. 2 = h Layer 2 0&1, E2, P2 h 2 -—v 00 Figure l. Layered Model a‘h’s' one: Cathy's th 29‘3de of 5 signal pr tie integra solution is and the bra ehich may t variables ' able. The A. Pl“ The V6 velocity . as A epr. i,k,t, and mg without nent wt-kr vat ‘ Kdl‘ : For a Map velOI ~Uppose a I tuber of I V V» may be 1 11 ways. One method is to integrate the inner integral by Cauchy's theorem on the complex plane and the other by the method of stationary phase. Since we are interested in the spectral properties of the transient response, only one of the integrals has to be solved. The contribution to the solution is from the poles within the contour of integration and the branch lines. The poles are the zeros to P (w,k) which may be real or complex. Note that P (w,k) has two variables, w and k, the zeros being functions in one vari- able. The locus of a zero is the dispersion curve of a mode. A. Phase Velocity and Group Velocity The velocity of wave prOpagation is called the phase velocity. Assume that a propagating wave is expressible as A exp(j(wt-kr)) , where A is the amplitude and w,k,t, and r are defined as before. For this wave, travel— ing without continuously changing its wave form, the argu- ment wt-kr must be kept constant. This requirement leads to: wdt - kdr a O or C a dr/dt = W/K a Phase velocity....(2.z) For a diSpersive wave train, there is a velocity called group velocity U which is different from the phase velocity. Suppose a wave is made up of superposition of an infinite number of waves with continuously changing wave number k. It may be expressed as: 00 Mint) =I A'k) €Xp\j\wtk)t-kr)) dk...............(2.5) oO The conc 5:14) whi garticul scall fo the wave tax,“ [13 y expat: 372d negl 'fi'hEre f( me quar tud 80f l «E ' 5‘05 tc 12 The concept of group velocity applies only to those cases A(k) which have a significant value in the neighborhood of a particular wave number, say ka, and becomes vanishingly small for k outside a small range, denoted by kaiek . Then, the wave function can be written approximately as: 1“," Ak u(x,t) z/ik A(k) exp(j(W(k)t - kr)) dk ...........(2.4) “‘4 a By expanding w(k) about k = ka, W\k) = W\ka) + \dW/ak&=Ka(k - Ka) + 0000000000 and neglecting higher order terms, equation (2,4) becomes h+Ak u(x,t) ./’ A(k) exp(jf(k,r,t)) exptjgtr,t)) dk’...\2.b) kg -—Ak ow . Where f(k.r,t) =‘k-Ka)(dik§;r) and g(r,t) = Wtka)t Kar- The quantity A(k) exptjftk,r,t)) is the effective ampli- tude of the wave. The requirement of constant phase now leads to fkk,r,t) 3 O or \dW/dk)k-K 3 r/t 8 U oooooooooooo\zob/ a for the group velocity. Thus the energy of a wave, which is associated with the amplitude, travels with the group velocity 13 U. Equation (2.6) can be rewritten as: w d0 U 3 (C/(1- C E;))W=wa and U=Ca'°°°°°°.'°.°'.°°°°‘2'7) where lca _ wa/Ua . For a dispersive wave train, C = constant or dC/dw = 0 hence U = Ca° B. Period Equation The denominator of equation (2.1) is the period equa- tion which determines the diSpersion pattern of the model. A form convenient for computation is given below for refer- ence (Ewing, 3E.al, p. 193). We will use this form for computation throughout the study. all" i1“: 0 ...........................\2.8) l where 9- (2- 7;“ )(X cos r, h 4»?! sin r, h) ’ I + 2% T" sm s,h -% Z cos s,h) E 31"(2- 13'— )(é‘hfl cos r,h +11 Z sin r,h) + 2%(x sin s,h - 553! cos s‘h) V I 1 00(209) ‘1'=(2- 7;; )(—EW cos s.h +4.1. Z sin s,h) + 219(X sin r,h -—% Y cos r,h) man (2- {4 )(X cos s,h +§fY sin s,h) + 2+(éfw sin r,h -.'rf.z cos r,h) L z .. ILII Univ/u Hour Q. = rig II . Adiffers 1'. t n 3 2' z X - E3 116.} -2 [:2 -1) Y = LE! +2(E2-1) I k2 I“ kl “I 1 l .............. (2.10) z=29flz-fi;-2(Ei—1) w=2(/‘_‘}-1) ”I k1 It I ”I r? = k:'—k2 r: = k2 -k;1 ........................ (2.11) 82 = k2 2 2 2 rigidity of the layer 1 Hi\= u = frequency “1 = compressional wave velocity in the layer i éi = shear wave velocity in the layer 1 _ 2M kd‘- d; k = 217 I t 7.. A different form of the period equation of a wave guide is given by Tolstoy and Usdin (1953) as follows: | N-M = o ...... . ...... . ..... ... ..... . ..... .........(2.12) ‘l'IBICI and SS IESpect Efficients Comir ask'mstotic PP 193-1g Obtained ‘ Shown tha‘ ZEIQ' equé tiOn fOr . {. I Ike" hlt beCOmeS f. Rayleigh I ZEro Of t} intErfaCe As TI 15 where I I denotes the determinant = A exp(jr.h) C exp(jr.h) 2 13 N B exp(js' h) D exp(js'h) ooooooooo o oooooooo coco ( o ) = G exp(jr.h) J exp(jr.h) M (H exp(jS, h) I exp(j8,h) oooooooooooooooooooooo (2.14) ' . . -._ I exp(js h) J exp(Jr h)- - M - I;TI i-H exp(js1h) G exp(jr:h)’ inverse of M.. (2.15 A,B,C, and D are reflection coefficients of PP, PS, SP, and SS respectively at 3:0. G,H,J, and I are reflection co— efficients of PP, PS, SP, and SS respectively at B = h. 1. Normal Mode Roots of Period Equation Coming back to equation (2.8), let us discuss some asymptotic roots of this equation (Ewing, SE 31, 1957, pp. 193-196). Positive real values of Si’ ri are obtained when C > «05. and C 4 é< 0K2 . It can be shown that when the thickness of the layer h approaches zero, equation (2.8) becomes the simple Rayleigh wave equa— tion for the half-space. For c:o means trapping of P wave energy in the surface layer. The roots of this equation according to Gilbert (1964) FBI. are: S roots on the Riemann sheet (+,+) (i.e.,Re r2)>o, Re .sz>o)and (-,-) (i.e.,Re r2(0, Re 5240. and 15 roots on (+,-) and (-,+) sheets. S roots eventually become Ml branch of Rayleigh waves and P roots become M2 branch of Lg Rayleigh waves. The P +- and F-+ , compressional and shear pipe organ modes, respectively (not the true "pipe organ" modes) have been identified on earthquake seismograms. The leaking modes found in our data are similar to the normal mode propagation in a liquid layer overlying a liquid half-space in which only P waves enter the problem (liquid cannot sustain the shear motion) (Ewing, g5 gl, 1957, pp. 126-151; Officer, 1958, pp. 117-145). When C3 approaches (X1 , or equivalently the incident angle ap— proaches normal, the normal mode equation becomes the equa— tion for the simple pipe organ modes whose spectral band peaks at 1, 3, 5, --- times the fundamental frequency (Ewing, g2 gl, 1957, p. 185). This phenomenon has also been reported by Grant and West (1965, pp. 104-107). It is strange to see that the above phenomenon exists in the solid layer overlying a solid half-space, because in 2‘ .rv 1.. .H.« :~ 2‘ 22 solids, both P and S waves in the general case couple to each other at the interfaces. Two possible explanations are discussed here. One explanation may be decoupling of P waves and S waves due to their high angle of emergence in the surface layer. The critical angle of the P wave in the layer is 11° 34' and the shear wave generated by the com— pressional wave impinging at the boundaries has an incident raga angle of 6° 00'. Both are small. Another explanation is based upon the assumption that the reflected S wave from water-table has small amplitude because of the small S wave velocity discontinuity at the saturated-unsaturated boundary, L”; which has already been discussed in some detail in Chapter I. It is difficult to determine which explanation is correct. At near normal incidence, impinging SV waves generate neglegible amplitudes of the reflected P waves at a boundary, but the impinging P waves generate about equal amplitudes of reflected P waves and S waves (Muskat and Meres, 1940). This infers that X component amplitudes must be similar in magnitude to the 8 component after the arrival of (X1 . This is not the case we have observed on our records. It is observed that the X component traces have smaller ampli- tudes than the 5 component traces after the arrival of the refracted C£. A portion of record No. 43 is shown in Figure 2. Howeveq,we still have grounds to argue that since S wave energy is more easily dissipated in the near surface loose material, the amplitudes of S wave shown on the seismograms are more severely attenuated than the amplitudes of P wave. 23 me .02 950mm .45 onHmom 4V. .N 55on C..m OHIV AIHHE O.m V > :(b '1?! .54. .r‘V “ ‘Qmufifimw? hush. $34...» ...Vsé V. k... V. gnaw . 9 «a , nigh); .. f. ’ ’g ‘V’. V V’ . ’1’“) V I; r? kWh... AV»... ... «4.... 2.“...me mm. «waffle , . ,. .V. l 7% «......4... f a #&.w 34......) «numwfigV ......Vme ...V ......V was” 4% V... .V.. {“4 wen E... .: 2...... ...... . 4.745....«3... 2...... 4V. V... «V “V... VVVVVJ V7. % . V a: _ (INS/W um HNE ON 0.. Inommvmsfle 9‘- (f (_I) (7' Ft IN 24 If the S wave velocity change at the water-table is assumed to be small and the density change is also assumed to be not substantial, an interesting result can be derived from the period equation of the single surface layer over a half-space model. The existence of the dispersive "pipe organ" modes in the assumed earth model is then explainable. The derivation is made possible by simplifying the period . equation under the assumptions: [9" =fizand f, = fl . r.” In the following, we want to show that the normal mode equa- tion of the liquid layers is a possible equation governing the dispersion pattern of the leaking modes propagating in i an unsaturated layer overlying a saturated half-space if the above assumptions are true. By substituting £54.: w» 'and F. in equation (2. 8), it becomes (2 - )[Cos (r h) + 3%. sin mm] 2": + 43'/C:_I My: [—3 Cos (r h) - ssn (r'h)] = o ....... (2.30) Following Gilbert's (1964) method, we assume f= ir+’j +i , thenrh=k’%_: -1h=‘3_“_f g} -1h= d2 C 211(frijj-i) ’Q —1 h = R + jI. Equation (2.30) can be C (7‘ I written as: 25 1T. 2‘ + “[222 = o .................................. (2.31) where TT = (2- 9—2)2 Cosh (I) + 4/93-1 _C__z—1 sinh (1).... (2.32a) fi‘ (3* ) (if 5:): C05 (R) + J 1-Cz/dz’ sin (R) .................. (2.3213) P‘“ ; JCz/af~| V2: (2" E- ) sinh (I) + 4/9-1-4 C... Cosh (I)... (2.32c) .‘ )5: a: :2: WW Cos (R) - sin (R) .................... (2.32a) C3&:-| I For equation (2.31)to be satisfied, four possible branches of roots can be solved from the four sets of equations listed below: Th = O and W} = O . . . . . . . . (2.33a) 1T,=o and 22:0 ........(2.33b) X,= O and W2: 0 . (2.33c) XI: 0 and 22: 0 . . . . . . . . (2.33d) These equations may or may not have solutions. Equation (2.33c) can be explicity written as tanh( gjLLLh-Cos e) = _4-94fq C/Rffl ................. (2.34a) 0‘) (2-Cz/@:)1 up f! L) r“? .L‘... ."n‘ “a. 26 and ‘/C‘/o(:: tan (m C059) =- ‘ ......... (2.34b) OI, /'-—Cl/o(: where 6 = incident angle or P wave in the surface layer. Equation (2.34b) involves only fr which is related to the dispersion pattern, whereas equation (2.346) involves only fi which controls the attenuation of the amplitude. Note that equation (2.34b) is exactly the normal mode equation in mgj a liquid layer overlying a liquid half-space with a density ratio of unity. If the above analysis is correct, the existence of the pipe organ type modes on land will be an indication of the existence of water-table. More observa- tional evidence is needed to verify it. Chapter III Spectral Analysis and Data Interpretation Spectrum, dispersion, arrival time, time duration, ampli- tude, particle motion, cross-spread velocity, etc. are useful in identifying a signal and in deducing model parameters. A body wave usually has a definite arriving time and short signal length. In contrast, the arrival of a normal or a leaking mode is not as sharp as the arrival of a body wave and the time duration is always much longer. Therefore, in studying a body wave, the determination of the arrival time ‘1 {Artur . is important. On the other hand, in studying the normal or ~ leaking mode, the derivation of the dispersion pattern is stressed. In this chapter, we will discuss methods of iden— tifying signals or deriving group velocities, and present the properties of signals identified on our records. A. Methods Used to Identify Signals or Derive Group Velocities Various methods used in identifying signals or deriving group velocities are introduced below. The discrete Fourier transform (DFT) and discrete inverse transform which will be used in these methods are given as follows: (Gold and Rader, 1969) N-l F(m(2) = Z: f(nA) exp(—jmnAfl). ..... (3.1) “:0 l N-( f(nA) vii-é; F(mf2) exp(jmnAQ) (3.2) 27 28 where FfrnC?) = discrete form of the Fourier transform {(n A) = discrete form of the time function P4 = total data points in the time window ‘A-= sampling time interval J 2: .— = 2 W/(NA) m = O, l, 2, ...., N-l n O, l, 2, 0000’ N-lo The same formulas will also be used in other chapters. Note ' that for a particular data points N and a particular sampling : interval 43 , there is a particular set of frequencies corres— E ponding to them. For example, N=4OO and A§=0.005 sec., the flea corresponding frequencies are m£V(2fi) = 0, 1/2, 1, 3/2,..,399/2. 1. Moving Window Spectral Analysis Method A typical feature of a surface wave record is that many modes are superimposed on each other and the spectral proper- ties of each mode varies with time at different rates. If one is able to determine the exact amplitude and phase spec— tra at each time instance the separation of modes will easily be made and the group velocity dispersion curve of each mode will readily be found. Unfortunately, this particular method is not a feasible task. The uncertainty principle in spectral analysis states that the product of the spectral bandwidth (a measure of the bandwidth of the signal) and the time duration of a signal cannot be less than a certain minimum value (HSu, 1970, p. 229). That is, the shorter the time window the poorer the frequency resolution, and the converse is also true. A compromise must be made between the length of the time window and the frequency resolution. 29 a. Parameters Involved in Use of the Method The parameters considered in designing this method are: 1) the minimum time window for which the discrete Fourier transform is applicable; 2) the time window that is required to resolve the spectrum of our available records to the ex- tent that each signal of interest is visually distinguish- able; 3) the weighting function; and 4) the increment of window shift along the time axis. To determine the minimum time window required, such that the spectrum, corresponding to a particular frequency mC2/2TT is derivable by using the discrete Fourier trans- form, consider the Fourier transform of a time function f(n.A ). Let us rewrite equation (3.1) as follows: "1:. 2___Tm) Wm" 3L“) F(m{2) = >_ f(n9) exp ('3 Nfin = Z: f(n‘° ex? ( j“ Nfin w“ nzo n=o 2’ " N ‘+ 2 f(n43) exp (- jr 3) + ..+'"bh%'f(nA) exp(- j: %r)..(3. 3) =Nfi" n:(m¢uMhO Note that the exponential term exp (-j 4%%%;L ) is periodic. For some i = 1,2, -—- , or m, an arbitrary summation in equation (3.3) can be written as: iMhfl-I znn NAn—I f(n4) eXP(-j W737.) =2 f[(1+(i-1) (N/m)) A] nsu-uxfim) 211 (:0 exp [ amour-1) (N/m))) Nfin-I = f[(l+(i-l)(N/m))é 1 - exp [- 313% + 21r(i-l))]. . . L=o Nfin-l 1:: f [(1 + (1-1) (N/mHAl-exrk" 52% L=0 l I... (3.4) 30 because i is an integer. Define fi(14) = f [(l+(i—1) (N/m))A]. The time origin of fi(lA) differs from the time origin of f[(1+(i-l) (N/m))4] by ('1-1) (N/mM. i.e., fi (14) F».- is a time-shifted function of f(1A) . Also define }§({23 to be the Fourier transform of fi(LA) for the data points n=(i-l)(N/m) to n=i(N/m)-l assuming that N/m=integer We can write LWJ N/m-) 1 Pi ( Q ) = 2: fi (1A) exp(rjg-l) L=O Wm lh-' 2 = fi (14) exp(-jm1A.J-’) . . . (3.5) NA [=0 N ‘4 i = fi (1A) exp(-jmiAQ ). L=O / It is concluded that Q is equal m0 . Therefore both Fourier transforms, F(mQ) and Fi(m(2 ) , are defined for the same frequency ’%%2= {14‘ , except that F(m(2 ) is defined for the data points from n=o to n=N-l and Fi(m{2) is defined for the data points from n=(i-1)(N/m)to n=i(N/m)-l only. Note that the number of data points in Fi(m(2) is N/m which corresponds to a time interval of (N/m)43 31 This time interval is exact the period of the frequency m/(NA) defined by both Fi(mfl) and F(m(2) . Thus Fi(mzs) can be termed One-period Fourier transform. Equation (3.3) can therefore be written: F(m(2) = F|(mQ) + F2(mQ) + + mem, .......... (3.6) namely, a m-period Fourier transform is a sum of m one- period Fourier transforms. Note that for N/m gt integer, equation (3.6) is not valid. However, from the practical point of view, it is possible to approximate it. For example, N=100, m=3, and A=l, then 0 = 2Tr/1oo and !§_' 3“L| F(mQ) = F(3%) 2’: f(n) exp(— BjRTg-g) + f(n) exp(-3j '2—0—0) ... M ”"239 a... + Z: f(n) exp( 33.17%: [Z f(n) exp(- -3j’—’3-;7) + 1/3f(33)exp(- -3j§--6W)] +[2/3f(33)exp(-3jffy')+ ff: f(n)exp(- -3jf2¥) + 2/3f(66)exp(- 3j7;;)] +[1/3f(66)exp(- -3j'—-——3:")+ Z? f(n)exp(- 3%)] "=67 = Fl(3L 11') +-F2(3% E) +-F3(3fl ), I00 [00 32 where F,, F2, and F3 are the sums in the brackets [ ]. This approximation was tested on an exponential function. The maximum phase error observed was i 0.003 radians, or 10.05%. This is negligible in practice. The percentage amplitude error corresponding to this phase error can be shown to be [l-Cos(0.0005d)]XlOO , where c1 is the phase angle. It is a very small amount. The one-period Fourier transform described above is a DFT that is defined for the shortest time window, but it does not guarantee a satisfactory frequency resolution. It has been verified experimentally in this study that the five-period Fourier transform is a satisfactory one for our particular data. The programming method is as follows: Firstly, N data points are divided into groups in subsequent order. For a particular frequency m/(Nzi), the number of points in each group is N/m, which corresponds to a time duration of N.A./m, i.e., one period long. For example, N=400 and 43:0.005 seconds, the number of points in each group are 400/m= 00 , 400, 200, 133 1/3, ..., and 2 for m=0, 1, 2, ..., 200 (the maximum possible number in the first group is equal the maximum number of the available data points, i.e., 400). The corresponding periods are (400x0.005)/m=2/m=co , 2, 1, 2/3, ..., and 2/200, and the corresponding frequencies are m/2=0, 1/2, 1, 1 1/2, ..., and 100. Secondly, the one- period complex Fourier transform is performed on each group of data points. Note that if the number of points in a group is not an integer, the approximation method mentioned 33 previously must be used. A five—period Fourier transform is then obtained by summing five one-period complex trans- forms to represent the complex spectrum of the central period of the five. To find the complex spectrum of the next time section, the five-period Fourier transform is shifted one period to the right (i.e., direction of increasing time). The well known weighting functions such as Bartlett (triangular), Tukey (hanning), Hamming, Parzen, etc. are not easily programmed in the above process. A convenient method is to weight the five one-period transforms by 1/9, 2/9, 3/9, 2/9, 1/9, respectively and then sum them up to H“ 222 form the five-period transform. This kind of weighting is a rough version of the triangular weighting. b. A Special Method for Improving Time Resolution Since the shifting of the time window is in increments of one period, a further improvement in time resolution is desirable. The concept of Page's (1952) instantaneous power spectrum was used to accomplish this goal. The intuitive meaning of the so-called "instantaneous power spectrum" is the contribution of a single data point to the power Spectrum of a time function which extends from minus infinity to a certain prescribed time. In application, it is assumed that the beginning time of a one-period time section, defined previously, is the time of appearance of a one—sided, (i.e., zero displacement before time of arrival) band-limited signal and the complex spectrum defined for the same time section is contributed mainly by this signal. 34 The following is the definition of Page's (1952) instan— taneous power spectrum: Parseval's theorem for the energy of a continuous signal may be written as: oo 00 2 Jg’mdt =J [F(f)] df ...................... (3.7) ~00 ~00 where g (t) is the time function, F(f) is the frequency function, u t = time, f = frequency, I i ( denotes the absolute value. " Let /7(t,f) denote the instantaneous power of a particular frequency f and at certain time t. The energy expanded from time minus infinity up to time t is t 2 so I A f9 (1,)d/t =f df f(‘t, f)d’t .............. (3.8) —-oo "°0 -oO by using the definition of‘F(f,f) and the identity, equa- tion (3.7). The instantaneous power of all frequency at time t is found to be 00 g.2(t) =j f(t,f) df...... ......... .................(3.9) —c» by differentiating equation (3.8). Define a causal time function 9(1)=.[g(’13) °<¢§t .......................... (3.10) t 0 otherwise 35 The Fourier transform.of gt(t) is 00 t , Ft(f) =j gtVt) eXP(-jW1‘,) d ’L =f g<1>exp(-jwr)df.. (3.11) ..00 o ' The total energy of gt(f) is 0041 2 f o0 gt (1)dt=j1 g (1)d’[ =f d’lj f(f,f)df.... ........... (3.12) 0 - 0 ~00 .00 by using equation (3.9). It can also be written as 00.2 no 2 , Lifl’t’d’t =L|Fuf>| df ........................ (3.13) by using Parseval's theorem. By equating equations (3.12) and (3.13) we get fde F(1,£) dT=fiFt(f)'z df ................... (3.14) .... o w . . _ From equation (3.14), one can write t , - ~- \ ' [f(t pf) d1 = lFt(f) '1. 0000' 000000000000000000000 (3015) o . ‘ Differentiating equation (3.15), obtain» f(t,f) = _2f(t) I Ft-(f), Re [exp(jwt) exp(jut,f)] = 2f(t) ‘ Ft(f) 1 cos (wt +0(t,f) ......................... (3.15) where Cit"; denotes the phase of Ft(f) . 36 Equation (3.16) was used in computing the instantaneous power Spectrum. Note that F(t,f) may be negative, but the cumulative energy iff(1,f)df is always positive at all frequencies by equation (3.15). A negative f(t,f) indicates a temporary readjustment to the total energy. Although the positive fKt,f) does not necessarily indicate the presence of a particular frequency at time t, a successive "*“3 building up of the positive power is no doubt a good indica- i tion. With this in mind, the positive densities are reduced by an amount prOportional to the ratio of the negative sum to the positive sum and the negative densities are reduced by _ W; an amount proportional to the ratio of the negative sum to the positive sum and the negative densities are arbitrarily assigned zero values because they are not associated with the major signals in the time window. The isolated positive densities are again eliminated because they are apparently not associated with long duration normal or leaking modes. C. Uses of Moving Window Spectral Analysis Method The moving window spectral analysis method may be used for two purposes. One is for separating superimposed modes and for revealing time-varying spectra. Another is for de— riving group velocity dispersion curves. Figure 3 shows the second and third "pipe organ" modes, on a portion of the 2nd 8 trace of record No. 43, being separated by this method. This figure also indicates that the first mode (frequency about 10 HTB) does not appear right after the refracted (X2 arrival. The observed points in 3'7 Awom<>lmz; .n umber. mo inzmT D momT lo? 8 $071? mm .d> .5 (v ommtnmomv Km... $3 womn 3mm. enmn mmnn $3 I. u . . . 0 no 0. m: E omm mm? om_¢ 8.... 0mm 9% own 2% on. 3. HUD mm on 3 U U D NEE: 2 D rezmaqaux... D 38 Figures 4 and 5 are the peaks of the time—varying spectra of the subsequent time sections of the same record. Two superimposed modes are identified in Figure 5. The events, La, Lb, Lc, etc. are defined in Table 3, p. 52. By the definition of group velocity, introduced in Chapter II, a curve obtained by connecting the peaks of a time-varying spectrum is the observed group velocity dis- persion curve of a particular mode. Figure 6 shows the time- varying Spectrum of M11 mode. Due to the limited space, this spectrum is drawn separately on four pages in original sequence. The circles indicated are the peaks of the Spec- trum. The group velocity dispersion curve is obtained by connecting these circles. Figure 7 shows the group velocity curves of M11 mode found on several records by this method and also by hand picking. 2. Fixed Window Fourier Transform The fixed window Fourier transform is used as an auxi- liary means to the moving window Spectral analysis method. Occasionally, the latter does not reveal clear spectrum due to its limited length of time window, the fixed window Fourier transform must be applied using a longer time window. The triangular weighting function is used to reduce the side lobes resulting from truncation. Figures 8 and 9 are the Spectrum of "pipe organ" modes. 3. Determination of Particle Motion The particle motion is often a useful means in identi- '59 cull! Shana. >ozuaoug .32.... u 233» .36; 2.82: 2:53..» ez_>¢<>.ua: 52.53 :20 3:26 ... 232... auuo\ht 00: >h.004u> taoco o. a. v. n. u. ._ 0. 3:3... ozicgéa.» ..o 3.3.. u o.— o.. r x 4o 0...: .upcu3. >ozuaau¢a Auo<¢ h N 080006 .0 v .0! occuuc. lachoua» oz_>¢(>nul_.—. box—baa 3:3 u;u>u <¢houtn 02_»¢(>om2.h no nxh.ooau> aaoce .n acaci DOVE 1 (5| 3223533312; flwamn ' D h 4 o . u_ 5- . 2 T o W 4 . . I] J J) W1... fl/2/Mi%//S /22/¢.2%%/W1m o FEIL E R 43 FR!OUENCY(HTZ) F a I o (a 20 20 s o ' I g g I «our veL. ' ll Irv/sac U ”WU" -IO TO ~400 ucono) 1:, U“ ~4ou AID nova an (3.99) —-x — LINE counccrmc rams or Tm:- vmvma spam-nun 507 - (4.04) 496 _ (4J4) :r noun: 0. mom” a) 44 UENCYCHTZ) r") , % (31:53; It , E] 0000000 a - I X Ll 45 FREQENCY‘HT21—0 8 I0 )7. I4 I6 GROUP VEL. I l I I | FT/SEC 0 600 -' 500- X X‘pOX O 400- MOVING WINDOW SPECTRAL ANALYSIS METHOD SYMBOL RECORD NO. COMPONENT TRACE (NO. 45 z 2 -—— 43 Z ‘ DIRECT MEASUREMENT OF PERIODS AND ARRIVAL TIMES SYMBOL RECORD NO. COMPONENT TRACE NO. x 33 Z 2 0 33 z 3 o 43 z 2 1 43 2 8 FIGURE 7. OBSERVATIONAL GROUP VELOCITIES OF I)” MODE 46 LC RECORD NO. 33 COMPONENT: z TRACE NO. 3 DISTANCEzgag FT WINDOW LENGTl-IIO.75 SEC, PROM! NENT EVENTS SI-IOWN. THREE “PIPE ORGAN" MODES Lb AMPLITUDE SPECTRUM (ARBITRARY LINEAR SCALE) La /\ / \/ / /\/ V \/ \ 11 ‘2 ‘ 11 L2 4 '1 ‘2 55 103 16 21-5 263 32 373 423 FREQUENCY(I-ITZ) FIGURE 8. AMPLITUDE SPECTRUM(FIXED WINDOW) AMPLITUDE SPECTRUM (ARBITRARY LINEAR SCALE) 47 LC RECORD NO. 43 COMPONENT: z DIQSATCATNC'E'QéOSO FT WINDOW LENGTH: 1.2 SEC. PROMINENT EVENTS SHOWN: THREE *PIPE ORGAN” MODES LCI Lb I A A I j . A 1 ./\/— . 25 30 35 IO 15 20 FREQUENCY(HTZ) FIGURE 9. AMPLITUDE SPECTRUM(FIXED WINDOW) 48 fying a signal. The determination may be accomplished by superimposing the X trace on the Z trace and inspecting the phase lags between two traces. This method is illus- trated in Figure 10. The particle motion is one of the prOperties of events identified in Table 3. B. Interpretation of Data The refraction records that have been studied are F—T: listed in Table 2. The spacing between any two adjacent geophones, charge size, depth of Shot hole, and the dis- tance between shot and the spread are indicated. These records were recorded without automatic gain control (AGC), suppression, and filter setting. Because no part of the spectral band is negligible under theoretical consideration, any distortion, caused by the artificial means, Should be avoided. The high frequency noises such as those due to instrument, wind motion, etc. can be filtered out after a visual inSpection of the Spectrum is made. The data used for the study were digitized from the paper records. The errors introduced by the digitizing process are not easy to evaluate. To test for errors in digitization a trace was digitized twice and the average was determined. There was no notiCeable improvement in the result, thus we conclude that digitization errors were neg- ligible. Most of the data used was thus digitized only once. The characteristic seismic events having been identi— fied on seismic records by methods described previously are PARTICLE MOTION SEISMIC TRACES z INCREASING TIME ——+ l\l J / N\ / I ' I \ x I - ‘ ( ./ \ // / k (\x, . , \J ’ UV/ IV I I c a! I — Z TRACE ---- x TRACE FIGURE IO. DETERMINATION OF PARTICLE MOTION PARTICLE MOTION C\ ,. \3 /\ . EM 2’) . C/ O) 50 SEISMIC TRACES INCREASING Tl ME ——F Z TRACE —--- x TRACE 51 Table 2. List of Refraction Records Record Spacing Charge Hole Depth Distance From End Number (ft) (lb) (ft) of Spread (ft) 10 25 2 5 725 ll 25 2 5 725 33 10 16 10 969 34 10 4 10 600 42 50 8 10 2000 43* 50 8 10 2000 *Repeated version of record No. 42, except that No. 43 has lower gain setting than No. 42. onuauwuou O» wandsam OUOE HMEHOZ mummw Iona III mvv I mmw III o.na I o.v Hmfiuoz HHS wvoE Hmfiuoz m III own I How III N I m.ha m NH: OUOE Hmfiuoz .mouumm 0» .moum III ENS I osm III m.HH I m.h Husuoz MN: ucoc IomEoo x :A .Hmsu N mom mom uson< on I HA IIII IIII «ML 50H” Am UDJJUHHDK Rwy 0A opmum IIII com I oooa III m.ma I m.HH Hmfiuoz CA No :Owumscfiucoo Iona ..mmh umv mvma mucoqomeoo m .w .x ca occuflflm 1.0mm Isa cmflm .va A umo ooaa uoouwp ou onmsoo wwmum CODE mcflxqu IONA A.ooom oooa I OMMH III m.HH I o.ma wmum>om on any ombH .mmmuo mafia. and muumm IIII wsna I coma III m.HH I o.m Hmeuoz on Iouumm :OHumuso uuonm m IIII mwwm I momv III A.ooo~ paw oo.mm m NAI .cmmuo mafia. cum Rama umo new.~4 coflumuap uwonm m IIII mNNM I momv III A.ooo~ umv mm.m~ a n4 .cmmuo «mam. 6am imam um. pmm.o~ Hw>mH min—munimm m cahm IIII 0>Onm a ON Illl Illl NVO Eoum m wmuomummm xuou own Scum m vaMH IIII w>onm a on IIII IIII mxu m Uwuomummm meEmm :oHuoz Aowm\uwv.aw> A00m\umv Ammov Ammov :OHmuwmmwo uao>m OHONuHmm umoumm mmouu .Hw> .ND pcmm .voum .voum xdom nuqo>m mo umflq .n wands 53 listed in Table 3 (events on Y component are not included). La, Lb, ...Le are leaking modes arriving before the refracted shear wave fa . La, Lb, and Lo are actually mul- tiple reflections of refracted P wave (x3 . We have dis- cussed their Spectral properties and the underlying theory in Chapter II. Ld is a peculiar event. It is dispersive and has a cross Spread velocity similar to the P wave velo- city in the surface layer, i.e., 1150 ft/sec. It may be interpreted that when the direct P wave couples to Lc, it results in a low frequency, high amplitude disturbance which has a cross spread velocity of the P wave in the layer and also exhibits the properties of a leaking mode. Le, be- ing interpreted as the continuation of Lc, is solely based upon the similarity in Spectral bands. That 6; is being related with the refracted S, can be due to the following evidences: 1) high amplitudes appear in X component (at distances, 989' and 2050') indicating a disturbance with strong shear motion; 2) cross spread velocity and group velocity (more correctly, shot—receiver distance/arriving time) are both about 900 to 905 ft/sec, indicating a re- fracted event from the half-space; moreover, the S velocity in the half-space is found to be in this range by the inver— sion of observational phase velocities (discussed in Chapter V; 3) all normal modes appear after this event (the highest velocity of normal modes is S velocity in the half-space). The normal modes M21, M12, and M11 are identified on the basis of their time-varying spectra. The particle motion of A /—;I 54 M21 and M11 are found to be consistant with the theoreti- cal calculations made by Mooney and Bolt (1966) for the alluvium-case. There are many other events that are possible of identification on the records, for example, the one identi- fied by Bennett (1973) to be due to P—S conversion. The study was limited to those events that are easily detecta- ble by moving or fixed window spectral analysis, and those that have apparent relations to our single surface layer model. Further discussion concerning the events and their implications in the assumed model are found in Chapter VI. Chapter IV Methods of Computing Phase Velocities From Observational Data Various methods of computing phase velocities using ob- servational data are discussed in this chapter. The computed phase velocities will be inverted to yield model parameters. The details of the inversion is given in Chapter V. A. Peak-and—trough Method Having been mentioned in Chapter II, the phase velocity is the velocity with which the wave form prOpagates. On a seismic record, it is possible to correlate a particular wave form from one trace to another, provided that: 1) the distance between two stations is not excessive; 2) the noise level is not excessive; 3) no superimposed events are present. The phase velocity is computed by dividing the distance be- tween two stations by the time shift of this wave form on two traces and the corresponding frequency is determined by mea— suring the period of this wave form (Officer, 1958). The phase velocities of M11 mode obtained by this method, is Shown in Figure 11. This is the oldest and the simplist method of estimating phase velocities. It is observed on Figure 11 that the results are poor. B. Fourier Transform Methods The underlying theory is the same for all the methods that will be discussed in this section. Consider that a dis- persive event observed at station A is fA(t), and the same 55 56 PHASE VEL CC! 'I' Y 750 FT/SEC \ 1 740 - 730 - 720 *- 7|O - 700 - 690 L l l I 8.5 NTz 9.0 9.5 I0.0 IO.5 FREQUENCY -'- ‘ X PEAK-AND -T ROUGH METHOD (2ND ANO ETH Z TRACES) --- FOURIER SUM-AND-DIFFERENCE METHOD (2th AND GTH z TRACESI -- INTEGRATION OF THE OBSERVATIONAL GROUP VELOCITY DISPERSION CURVE (2ND Z TRACE) FIGURE ll. PHASE VELOCITIES OF M” MODE (RECORD NO. 43) 57 event observed at station B is fB(t). The Fourier transform of fA(t) and fB(t) are FA(w) and FB(w), respectively. The difference in amplitude spectra between FA(w) and FB(w) may be due to geometric Spreading, and material dissipation, and the difference between phase spectra is due mainly to the difference between arriving times, i.e., dispersion. The phase velocity can be found by the formula given by Sato (1960): Ad: 2m“: w At = znf (rB' rA)/C(f) or C(f) = 2TIf (rB-rA)/(AII:2mII) ........................ (4.1) where rA and rB are Shot-receiver distances at station A and station B, reSpectively,A°(= “A - 0(5 = phase difference, and m=0, l, 2, 3, ---, .At: time shift of the same wave form in two traces. l. Fourier Phase Difference Method This method finds phase difference by subtracting the phase spectrum of FB(w) from the phase spectrum of FA(w) directly, and then uses equation (4.1) to find phase velocity. It was applied to the data of this study, but the results ob- tained were very poor. The reason is probably due to the short time duration of signals and to rapid change of spec~ tral properties. 2. Crosscorrelation Method The crosscorrelation functions of two functions fA(t) and 58 fB(t) is defined as: (Hsu, 1967) 00 RAB(‘E) I...) fA(t) fB(t-’£) dt and 00 RBA(1) = 100135“) fA(t-’£)dt the Fourier transform of RAB(T) is FA(W)'FB('W) = 'FA(W)' . ' FB(w)' exp (jam) and the Fourier transform of RBA(T) is PA(-W) . FB(W):1 FA(W)‘. IFB(w)\ exp (-jA°() . Theoretically, the crosscorrelation can be performed either in frequency or time domain. Landisman, gt al.(l969) claims that the phase velocity dispersion curve obtained first by crosscorrelating two functions in the time domain and then taking the Fourier transform are smoother than those found by the Fourier phase difference method described previously. This method was also used on the data of this study by means of computer program provided by Dr. R. S. Carmichael of the Department of Geology. The results did not seem to .be much better than the results using the Fourier phase difference method. 3. Fourier Sum-and-difference Method This method, instead of making use of the phase differ- ence directly, calculates the amplitude Spectrum of the sum 59 or difference of the two functions fA(t) and fB(t +zxt), where zit = time Shift with respect to the time origin of fA(t). The Fourier transform of these two time functions are FA(f) and FB(f) exp (jw13t) , respectively. The phase velocity is found by C= (rB-rA)/At ...................................... (4.3) and the correSponding frequency is the frequency that ex- hibits maximum or minimum amplitude spectrum, depending on the sum or the difference being used. The details of this method can be found in Bloch and Hales (1968). A brief explanation of how this method works is given below for reference: Let S(f)= F5(f) exp (jWAt) + FA(f). D(f) = F5(f) exp(jwat)-FA(f). Assume the amplitude Spectra of FB(f) and FA(f) are the same (in application, this can be accomplished by normalization), i.e., both equal some constant A. S(f) can also be written as: S(f) = FB(f) exp(jwat) + FA(f) = AIeXPdeBI-eprijt) + exp(jo‘AH = A [(Cos(o2 + (sin(°(8+wAt) '.._sin (dwwow—Ismm ] + sin (13);]y2exp [j tan‘ COS (0% + WAN—ICosO‘A =A 2&[1 + Cos (0(8- %+wnt)]y2 . -I SI‘M‘D’OTW‘WflsmdA] """" ”(4.4) o QX ' “Ha“ co5(d3+wet)+cosd/I - 60 When the absolute value of S(f) is a maxima, Cos (UE- 0(A+ wot) = l or 046- aA12mn=Adt 2m11=wAt ....................... (4.5) It is the same as equation (4.1). Similarly D(f) can be written as: D(f) = A25 [l-Cos(0§3 'd/I + w At)]yz Sin (013 + wat) ~Sin0lA C05(olB +w4t)- CosdA . ‘exp [j tan“ When the absolute value of D(f) is a minimum, the same con- clusion results. In application, the ratio [D(f)] /|S(f)) is used so that the minima is more pronounced. A display of this ratio vs. time shift .At' is automatically printed out by the computer. The phase velocity dispersion curve is de- termined by correlating the minimum points of amplitude spectra correSponding to different time shifts. Care must be exercised in determining the dispersion curve because for any time Shift, the corresponding spectral minima is repeating for each period along the time-shift (or phase velocity) axis. This can be easily explained by using equation (4.1): 61 C(f) = ZUf (rb- rA)/(Aolj'2m11) = (r8 - rA)/(Atjm-period)o For m=0,l,2, ---, there are infinite values of phase velo— cities.’ This ambiguity can be resolved by checking with the results obtained by the peak—and-trough method described previously. There are some additional procedures which were used to improve the results and to reduce the computing time: a. Shifting Phase Spectrum The Fourier transform is a time consuming process. The method used here was designed to perform the time shift- ing in frequency domain without resorting to repeating the Fourier transform. It can be shown that the Fourier trans— form of fB(t + At) is exp (jwat) FB(w). For a time limited function, the multiplication of FB(w) by exp (jwat) is equivalent to the shifting of the time function. It is schematically illustrated in Figure 12. This method has two advantages. Firstly, one needs to perform Fourier transform only once for the two time functions, fA(t) and fB(t), respectively. The Fourier transforms of fB(t +At.), fB(t +413) , fB(t +At2) , ...etc. are replaced by multiplying exp (ijt.), exp (ijt,), exp (jwatz), ... etc. to FB(w). The saving in computer time is considerable. Secondly, the Shifting of the windowed time function will not introduce any additional noise as will the actual shifting in 62 Shifted Position of Samples Time Increments in Time Domain 0 O 1 2 3 4 5 6 7 1 7 0 1 2 3 4 5 6 2 6 7 O 1 2 3 4 5 3 5 6 7 0 1 2 3 4 4 4 5 6 7 0 1 2 3 5 3 4 5 6 7 0 1 2 Figure 12. An Example of Shifting in Time Domain T J , 1.0 -'H<- L 0.1 sec. Figure 13. Trapezoidal Weighting Function 63 time. Since most of the signals on our records are short in duration, the introducing of a few noises will distort the Spectrum substantially. b. Normalization of Time Functions This method of computing phase velocity is based upon the assumption that FA(w) and FB(w) both have the same am— plitude Spectrum. This is not true in practice. Two time functions, fA(t) and fB(t), were normalized with respect to the maximum amplitude in two functions. c. Trapezoidal Weighting Function A trapezoidal weighting function shown in Figure 13 is used to reduce the effects caused by trancations. Another reason for using this weighting function is to weight every sampling point in the window evenly, except for a few points near the truncations. This is done to avoid the possible distortion of the Spectrum. d. Time-varying Filter All surface waves possess time-varying spectra. To filter out noises, a filter with time-varying pass-band is desirable. The principles in designing this filter are the same as those described in Landisman, gt al.(1969). It was found that the phase velocity diSpersion curves, determined by using the filtered data, were better than the unfiltered ones . The method of applying this filter is given below: 'sl?‘h-I| IIIVZ—d ‘ 64 Firstly, a fast Fourier transform is performed on the time function. Using the discrete Fourier Spectrum obtained, one can compute the amplitudes and phases of the Fourier series which represents the time function. Secondly, by inspect— ing the time-varying spectra described in Chapter I, one is able to determine the pass-band of the signal at different time instances. Finally, the filtered signal is obtained by truncating each cosine function of the time series accord- ing to the shape of the time-varying spectrum and summing over all truncated cosine functions. The truncations of cosine functions are smoothed before summing by a half-cosine weighting function given below: L. W (t) = Cos (11 - L/2¢<>IulLP 9L.- ozt match ciguo no hamloua ...-(80 < .¢L mass... L L I L L x x L L L x L L L L L L L L L L L L L L L L L L L L L L L L L L L L . L L L L L L L . L _ L L I L L L L L L \ L — L L \ L L L L L L L L L L L L. L L L L L L L L L L L L L L L L L L L L L L LL L L L LL L L L L L L L L L L LL L L L L L L L L L L L LL L L L L L L L L L L L L L L L L /I\ L L L L L L 66 example of the filtered and the unfiltered M11 Signal. A phase velocity dispersion curve obtained by the Fourier sum-and-difference method is shown in Figure 15. The maximum Spectrum is normalized to 0 decibels. Those below -50 are set to zero and those higher or equal to '25 decibels are represented by X. The asterick (*) repre- sents the lowest and the next lowest points of the spec- y“ trum. The solid curve connecting the low points in the F—— trough is the phase velocity diSpersion curve. There are several possible curves in this diagram. This one is chosen because the phase velocities associated with this curve are similar to those found by the peak—and-trough method. C. Computing Phase Velocity Dispersion Curve by Using the Knowledge of Observational Group Velocities The use of this method is two-fold. Firstly, it may be used to check the correctness of the results obtained by other methods. Secondly, the results of this method is ex- ‘pected to be less scattered because the group velocity dis- persion curve, obtained by the moving window spectral analy- sis method, is relatively smooth. Equation (2.7), which relates the phase velocity to the group velocity, is a first order differential equation. The :numerical technique used to solve it is the Runge-Kutta Inethod of order 4 (Conte, 1965). The details of this method is given in Appendix A. The group velocity data are those from the moving window spectral analysis method. The disadvantage of this method is that an initial value 67 en‘s: sz. (FT/IE6) 512.821 X9XXXXXXYXXX9885*5566874k377776646468‘XXX99XXXXX 510.431 9Bx9xxxxxxx98775116124663378787767677*Xxxxxx9oxx 526.316 XXXXXXXXXYKXXXX96433RI9988XXXXXXXXXXX7XXXXXXXXXY 533.333 x9xxxxxxxxxxx7x87655539299xxxxxxxxxxxBxxxxxxxaxa 540.641 xxxxxxxxxxxxxxxxxxxxxxxXXXxxxxxxxxxxx9xxxxxxxxx7 547.045 97xaxquxxxxx9998777779xxxxxxxxxxxxxxse79XX9xa76 555.556 6494875xxxx9a7765555557xz9xxxxxxxxx9943244536715 563.380 xxxxxxxxxxxxxxxxxxxrxrXVXXxIXXXXxXXXX*7867B*XX6X 571.429 X9XBXX9XXXXXXXXXXXXYXXXYXXKKXXXXXXXXX‘794*569XBX 579.710 xxx9xxxxxxxyxvxxxyxnyx‘XXRXXXXXX9X1X~XX989XXXXX 588.235 xxx9xx9xxxxxxxxxxyxvxyxxxxxxx998877ex’XXXXXXXXXX 597.015 76x59359xxxxxxxvxxxvxrx99964433;i#*359xxK99X85X6 606.n61 55X37637X88XXXXXYYXYXV75653OOV*0011336X9XXXX7365 615.385 xxx7xx7xxxx8565435678899195X8XXX99067 634.921 xquxxaxxxxxxxxyxvxvx*27§{(3k l.90u363 0 *9 Q/ .9 is an example. Although the programming of this method is laborious, the computer time needed in computation is surprisingly small and the results are close to those obtained in preceeding sections. The model parameters computed by this method are listed in Table 6. The S wave velocity in the surface layer was not computed by this method. Because the value 603.3 ft/sec found by the exact method in previouu section is very close to the value 600 ft/sec determined by Bennett (1973), it was thought that these two values must be close enough to the true value. 80 Table 6. Model Parameters Obtained by the Inversion of M11 Data Through the Use of the Method of Minimi— zing Sum of Squared Residues (Record No. 43, 8 component) Second Trace Rigidity Ratio 2.26 Thickness (Surface Layer) 21.95 S velocity (Half—space) 902.00 Density Ratio 1.01 Poisson's Ratio (Surface Layer) 0.31 Poisson's Ratio (Half—space) 0 .49 Sixth Trace 21.34 902.47 81 2. Checking the Correctness of Inverted Model Parameters by Reversing Process The best method of checking the correctness of the inverted model parameters is to reverse the process and to check the differences between these dispersion curves ob- tained by the reverse process and the observational disper- sion curves. The reverse process is just the traditional method of deriving theoretical phase and group velocity dispersion curves (Ewing, gt gt, 1957). The details are given in Appendix C. Table 7 shows a comparison of the theoretical and the observational phase and group velocity dispersion curves. The same data is plotted in Figure 16. B. A Discussion of the Model Parameters Derived From Normal Mode, Leaking Mode, and Compressional Wave Data The P wave velocity in the surface layer has not been computed in the study due to the difficulty in determining the accurate times of the direct P wave arrivals. The value of 1150 ft/sec has been used throughout the study. It was determined by Bennett (1973) using short distance refraction surveys. In his study it was also found that a low velocity loose sand with variable thickness was present near the free surface. The P wave velocity in the saturated layer was computed by the standard refraction time-distance curves to be 5700 ft/sec. This velocity and the P wave velocity in the Mississippian sandstone, which was found by the same method, VELOCITY (FT/SEC) l TOO - 700- 630 " SOO ' 550- 9 SOO - 450 . . I . T S 9 PREOUENOY( HTZ) —-+ 82 O THEORETICAL OISERVED(UOVINO VIRDOW SPECTRAL ANALYSIS UETI-IOD) E OSSERVEDUNTEORATIOR OF OISERVED OROUP VELOCITIES) {/Ffldfléfli \flELCXZITY’ R I at X 0 (SFKDLWD \4ZLCXZFIY / + O O O 53 :2 ° I'D I'I 1'2 1‘3 1'4 1'6 FIGURE IO. OBSERVED AID THEORETICAL DISPERSION or u.,uoo: 83 Table 7. A Comparison of the Theoretical and the Observational Phase and Group Velocity Dispersion (Mll data taken from Record No. 43, second 3 trace) Group #Group Frequency Phase Velo- *Phase Velocity Velocity Velocity (HTS) city(Ft/sec) Difference (Ft/sec) Difference 8.50 740.6087 -0.1323 582.5392 13.0952 9.00 728.3312 0.7878 554.3113 3.9753 9.50 715.2862 0.4978 528.1900 -6.3600 10.00 701.7815 -0.2535 505.6809 -8.7481 10.50 688.2109 0.1815 487.8027 4.8817 ll.00 674.9856 0.3324 474.8996 -0.7384 11.50 662.4602 —0.0458 466.6908 -6.7502 12.00 650.8844 -0.5699 462.4858 -6.0852 12.50 640.3907 -0.8555 461.4311 -2.3699 13.00 631.0123 -0.7912 467.6952 -2.5382 13.50 622.7123 -0.2058 465.5661 13.0571 14.00 615.4122 0.9882 469.4813 23.8293 Theoretical phase velocity minus Observational Phase velocity (by integration of group velocitiesL Theoretical group velocity minus Observational Group velocity (by moving window spectral analysis method), ' 84 are believed to be reliable because both arriving times are accurately determined in the study. The S wave velocity in the surface layer has been de— termined to be 603 ft/sec by the inversion of the M11 data. It is surprisingly close to the value 550 to 600 ft/sec, which was measured by Bennett (1973) using hammer blows as sources and at short distances less than 100 ft. It was also computed, by using M11 data, that the S F‘ wave velocity in the half-space of the single surface layer model is 902 ft/Sec. It is different from the value of 1200 ft/sec for the saturated layer by Bennett (1973). In Chapter II other evidences that support the existence of the layer with 902 ft/Sec, have been discussed. Having been mentioned before, the clay layers in the Pleistocene glacial drift are possible of causing larger shear wave velocity changes than the water-table. The driller's well log indicates that the depth to the water-table is about 32 ft (Todd, 1971). By applying the standard refraction technique, it was found to be 46 ft in Record No. 43, 41 ft on Record No. 10. Todd (1971) identi— fied eight successive multiple reflections right after the refracted P wave from the saturated layer and found the average time lag between two sequential reflections to be 0.056 sec. The depth corresponding to this value is 33 ft. After examining his record, the first time lag, i.e., the one between the refracted P wave and the first reflection, has a stable value of about 0.060 sec. This is equivalent a .. 5.....- t. . m r _ 85 to a depth of 35 ft. It is felt that the first time lag rather than the average value must be used in the calcula- tion. However, the depth, based upon the first time lag, seems to be more in error as compared to the one based upon the average time lag. The low velocity material right near the surface and the elongated ray path due to velocity gra- dient may cause this discrepancy. By using the first time lag, the depths are determined to be 41 ft on Record No. 43, 36 ft on Record No. 33, and 33 ft on Record No. 10. The depth can also be determined from the Spectra of the "pipe organ" modes. The first three modes have peak amplitude Spectra at 10.000, 23.333, and 35.000 cps in Record No. 43, 10.667, 26.667, and 42.667 cps in Record No. 33, and 11.429, 23.143, and 44.286 cps in Record No. 10. The true "pipe organ" modes are supposed to peak at frequencies with ratios of l, 3, 5, 7, ... The departure of the observed ones from the true "pipe organ" modes has been discussed in some de— tail in Chapter II. Since these modes are dispersive, their frequencies actually vary with time. The frequencies of the true "pipe organ" modes as their phase velocities approach the P wave velocity of the saturated layer. The fundamental frequencies of the true "pipe organ" modes, derived under this requirement from the observed data, are 7.000 cps for Record No. 43, 8.533 CpS for Record No. 33, and 8.857 cps for Record No. 10. These frequencies are equivalent to depths of 41.9 ft for Record No. 43, 34.4 ft for record No. 33, and 32.5 ft for Record No. 10. The depths to the water- 86 table computed by various methods are summarized in Table 8. The shot-receiver distances and geOphone spacings of the records used can be found in Table 2, p. 51. It is observed that the depths based upon the first reflection and "pipe organ" mode data are very close. Al- though the exact depths at each station are not known, it is believed that these two sets of values are slightly higher than the true values. The small positive errors are due to the low velocity material near free surface and the elongated ray paths caused by the velocity gradient. The depths de- termined from the refracted P wave are obviously in error. The intercept times used in computing the depths are ob- tained by extending the time-distance curves to intercept the time axis. The range of a T-D curve, that is actually determined from the observed data, is limited by the Spread of the geOphone array, which is small as compared to the shot-receive distance for the records studied. The extension of a T-D curve to intercept the time axis is inevitably sub- ject to error. .The depths found by using record No. 43 data are greater than those found by using other records. This may be due to the descending of the water-table because the recording site of No. 43 is farther from the source of the ground water than the recording sites of other records (per— sonal communication with Dr. H. F. Bennett of the Department of Geology). The thickness of the surface layer of the single surface layer model, computed from the M11 data of record No. 43, is 87 Table 8. The Depths to the Water-table Data *Refracted *First # "Pipe Organ" I Based P Wave Reflection Mode ' Record 1 1 Number 10 33 43 10 33 43 10 33 43 fi-.& Depth (ft) 41.0 ? 48.0 33.0 35.7 41.0 32.5 34.4 41.9 ? Difficult to measure the slope of the time-distance curve All traces in 8 component were used in computation # 8 component traces used in computation are: Record No. 10---4th and 7th, Record No. 33—--lst, 2nd and 3rd, Record No. 43---2nd and 6th, Record No. 42 (results included in No. 43) ---1st, 4th and 7th 88 about 22 ft. This value is apparently not the depth to water-table. The original thinking is that because most seismic energy is trapped in the unsaturated layer above the water-table, the W-T may act as the interface of the assumed surface wave model. By re-examining the tech- niques used in studying surface waves, no problems that may introduce errors have been found. In the same evalua- tion, the 8 wave velocity of the surface layer which has been mentioned previously, is verified to be correct and the density ratio at the interface is about 1.0. This in- dicates that it is not an unsaturated-saturated boundary because the porosity of the glacial sand in this area is about 25% (Bennett, 1973) and the corresponding density ratio at the W-T is about 1.1. It has been well estab- lished that the normal modes are several times more sensi- tive to the S wave velocity distribution than the leaking modes. On the other hand, the leaking modes are more sen- sitive to the P wave velocity distribution (Su and Dorman, 1965, pp. 1008-1009). It is logical to accept the fact that a S wave velocity discontinuity exists at a depth of 22 ft. This discontinuity may be a result of lithology change (clay layer), compaction, or cementation. Unfor- tunately, this result does not contribute anything to the controversial problem discussed in Chapter II: Does the S wave velocity change at the water-table or not? One possible inference is that the S wave velocity change at W-T is of a magnitude less than the one at the level of 89 22 ft, where the velocity increases from 603 ft/sec to 902 ft/sec. If a multi-layered medium is assumed to be a single layer over a half-space model in surface wave inversion, the interface of the model will be found at the level where the maximum change in S wave velocity occurs. That is because the diSpersion pattern of the model, using this level as the interface, is close to the actual diSpersion pattern of the multi-layered medium. Chapter VI Summary and Conclusions This study has been centering around the problem of the utilization of leaking and normal modes from a small explosion. The difficulties of this problem result from the presence of a high level of noise, an immaturely de— veloped wave form, and a short time duration of the sig- nal. The emphasis of the study has been placed upon the use of numerical techniques to achieve better results. Some problems pertaining to the theory, and the particular data used, have also been studied. A. Numerical Techniques All the numerical techniques used in the study have been carefully examined. To check the correctness of the results of a newly developed or modified technique, one or more other methods were used to the same set of data, and a comparison was made. In some cases, the testing of techniques was made by reversing the process. Because of these careful considerations, the techniques used in the study are believed to be reliable. This study may be help- ful for other workers who wish to study the seismic data of similar nature by numerical techniques. 1. Derivation of Time—varying Spectra and Group Velocities A special version of the moving window spectral analy- 90 91 Sis method has been proved to be an effective means of separating superimposed modes and identifying signals. The group velocity dispersion curves derived by this method are smoother than those derived by the traditional, visual measurement of periods. To save computer time, an approxi- mation method of the one-period Fourier transform was de— velOped, and Page's (1952) instantaneous power spectrum was implemented in the moving window Spectral analysis method. 2. Computation of Phase Velocities Some methods were found to give poor results using the available data. The Fourier sum—and-difference method LLJE using the time—varying filtered data gives a good disper- sion curve. This method uses the shift of the time origin of one seismic trace with respect to another. By using a trapezoidal window, and doing the time shifting in the _' frequency domain instead of the time domain, two advanta- ges result. One is the reduction of time needed in per- 411...; forming Fourier transform. Another is that no additional noises are introduced. The derivation of phase velocities from the observational group velocity dispersion curves yields smooth phase velocity diSpersion curves suitable for the inversion computation. 3. Inversion Methods I The exact method for the inversion of the phase velo— cities to yield model parameters is easy to apply. The method of minimizing the sum of squared residues is difficult 92 in programming, but the computer time needed is small. The results of these two methods were compared and a good agreement was found. To further assure their correctness, the inverted model parameters were used to compute the theoretical dispersion curves, and then compared to the observational ones. The differences between the theoreti- cal, and the observational ones are small. There are lar- ger deviations in group velocities at low and high fre- quencies. The higher degree of uncertainty in determining the correct values at the two ends of the observational group velocity dispersion curve may be the source of errors. B. The Use of Normal and Leaking Modes in Data Interpretation In all the records studied, M11 mode has the most well develOped wave form. For shot-receiver distance less than 1000 ft, the time duration of this mode is still too short for the accurate phase velocity dispersion curve to be easily determined. This is why only the results, computed from the M11 data in record No. 43 are shown in previous chapters. The M12 and M21 data may be possible to utilize when the shot-receiver distance is longer than 2000 ft. The use of normal modes not only gives an independent check on the results obtained by using body wave data, but also gives lithological information. The dispersive "pipe organ" modes found in the data studied, can be used to find the depth to the water—table. The results agree with those obtained by using the time 93 lags between the refracted P wave and the first reflec- tion. The depth determined from the "pipe organ" modes may be more reliable than the ones determined from the refracted P wave arrivals when shot—receiver distance is large as compared to the spread of the geophone array. A theoretical derivation was made in Chapter II to Show the possible existence of the dispersive "pipe organ" modes in the unsaturated layer overlying the saturated half-space. This theoretical treatment was not done in rigorous manner. An actual numerical calculation must be made and more observational data must be collected to verify its correctness. If the dispersive "pipe organ" modes are verified to appear only in this particular situation, their appearance will be a good indication of the existence of a water-table. In the assumed model, the boundary with large S wave velocity change becomes the interface of the model for normal modes, and the boundary with large P wave velo- city change becomes the interface of the model for leak- ing modes. This is in agreement with the conclusions of Sn and Dorman (1965). In this study the S wave velocity change at the water-table was not detectable. This does not exclude the existence of this velocity discontinuity. In future study, the Haskell's matrix must be used in the surface wave computation. It is more flexible in shifting from the single-layered half-space model to the Imulti-layered model. A more accurate result may be H LJ 94 obtained. C. Miscellaneous Conclusions and Recommendations In addition to the conclusions of Tolstoy and Usdin (1953) regarding the existence of the symmetric and anti- symmetric modes in a solid layer overlying a solid half— space model, a theoretical derivation was made to prove that they exist under all conditions. Only X and 8 components of the records have been carefully studied. The spectra of the Y component was only lightly studied. It was found that the spectral properties of the Y component are different from X and Z components. This is because the theory underlying the SH motions of the Y component is independent from the one for the other two components. Additional information may be obtained through the study of this component. Although linear filters are not suitable for studying dispersive wave trains, they are recommended for use in the detection of direct, refracted, and reflected body wave arrivals and to add more information to the study. APPENDICES III-F— Appendix A. Runge-Kutta Method of Order 4 Used in Computing Phase Velocities The problem here is to solve a first order ordinary differential equation of the form dC(f) ___ C(f) df fu(f) [U(f) - C(f)] = F(f,C) ....... . ..... (A1) with initial condition C we) a where C = Phase velocity U = Group velocity (an implicit function of f) f = Frequency fo = Initial value of f Co = Initial value of C. In computer computation a value observed on the observa- tional group velocity dispersion curve is substituted for U. Co is a phase velocity chosen from the results of Fourier sum-and-difference method. Note that equation (A-l) is the same as equation (2.7). Here, w = wa and C = Ca are understood. The recursion formula for numerical computation is: (Conte, 1965, p. 223) Cn+|= Cn + 1/5 (Kl + 2K2 + 2K3 + K4) 95 96 where K1 h F (fn, Cn) K2 = h F (fn + h/2, Cn + l/2 K1) K3 = h F (fn + h/2, Cn + 1/2 K2) K4 = h F (fn + h, Cn + K3) n = 0, l, 2, ... = Iteration step h = Step size of frequency. Appendix B. Newton's Method for Simultaneous Non-Linear Equations The simultaneous non-linear equations Fl (X1, X2, 000 xn) = o 0 F2 (X17 X27 ... Xn) II 0 PD (X1, X2, ... Xn) are to be solved. Assuming that F1, F2, ... F and all their n derivatives through second order are continuous and bounded in a region containing the true solutions (Al, A2, ..., An) and the initial approximation (a1, a2, a3, ..., an) is chosen sufficiently close to (A1, A2, ..., An). By expand— ing F1, F2, ... Fn about (a1, a2, ..., an), one gets Fl ( X1, X2, ..., Xn) = F1(a1, a2, ..., an) + eF1(al’ 32, 000, an) ( Xl-al) €9X1 + .9 F1(a1, a2, ..., an) (Xz-az)+ ;;Fl(al, a2, ..., an) (Xn-an) ia’xz E3 Xn + Higher-order Terms F2 (X1, X2, .00, Xn) = F2(al,a2' 0.0, an) +aF2(a1'a21 '0'! an) axl "Xl‘a1I 97 98 + .e‘FZ(al,a2, ..., an) (XZ-az) EDXi +- IaF2(al, a2, ..., an) (Xn-an) 9X1: + Higher-order Terms + LLLLLL Fn (X1, X2, ..., Xn) = Fn(al,a2,..., an) + F - 9 n (a1,a2,...,an) (Xl a1) 6X. +9Fn(al, a2,...,an) + 9F (a1,32,..., an) ex“ (Xn - an) + Higher-order Terms. By neglecting the higher-order terms, the above equations can be written as: 3 . , GE (Xla')+ 2:- (X1. a).)"\’ -— — '+ :—%(XH an):-'F- F2 1 z 5 .(X’a') * 3%(Xz’azh -- — +.‘;L,Fzg(x.-an):- F. I I I I I a 6 n 5%(XI‘QI) + 5%(X1—QL)T s —— ... TéiF—(XVI—an):’ Fr] 99 Therefore, a recursion formula is: (Conte, 1965, pp. 45-46) - _ _. .. .. 9_£I_ —FI a )(2 9X" I I .' I l I I : I I: 9E" _____ Q Fn " ex. 9 X 0 . L a . 1 _1____---__------.. “ >00.» III+I "t j02 or °"