University of Southern California Dissertations and Theses

Deterministic and stochastic modelling of the high frequency seismic wave generation and propagation in the lithosphere

(USC Thesis Other)

Deterministic and stochastic modelling of the high frequency seismic wave generation and propagation in the lithosphere

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content DETERM INISTIC AND STOCHASTIC M ODELLING OF THE HIGH FREQUENCY SEISMIC W AVE GENERATION AND PROPAGATION IN THE LITHOSPHERE by Yuehua Zeng A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillm ent of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Geological Sciences) December 1991 Copyright 1991 Yuehua Zeng UM I Number: DP28598 All rights reserved INFORMATION TO ALL USERS The quality of this reproduction is dependent upon the quality of the copy submitted. In the unlikely event that the author did not send a complete manuscript and there are missing pages, these will be noted. Also, if material had to be removed, a note will indicate the deletion. UMI Dissertation Publishing UMI DP28598 Published by ProQuest LLC (2014). Copyright in the Dissertation held by the Author. Microform Edition © ProQuest LLC. All rights reserved. This work is protected against unauthorized copying under Title 17, United States Code ProQuest ProQuest LLC. 789 East Eisenhower Parkway P.O. Box 1346 Ann Arbor, Ml 48106- 1346 UNIVERSITY OF SOUTHERN CALIFORNIA THE GRADUATE SCHOOL UNIVERSITY PARK LOS ANGELES, CALIFORNIA 90007 This dissertation, written by under the direction of h .l$ ...... Dissertation Committee, and approved by all its members, has been presented to and accepted by The Graduate School, in partial fulfillm ent of re quirements for the degree of D O C TO R OF PHILOSOPHY Dean of Graduate Studies Date SePtem^er i s » 1991 DISSERTATION COMMITTEE Chairperson Acknowledgement I feel very fortunate to have Professor Keiiti Aki as my advisor. Thank you very much, Kei, for your inspiration, encouragement and guidance throughout my entire study and stay at USC and for your stimulating ideas during this research. I have the highest respect for your scientific judgement and am always impressed by your enthusiasm in seeking the truth of nature and your efficiency in getting work well done. Equally I would like to thank Professor Ta-liang Teng for his innumerable guidance, and his generous help and support. I am grateful to Professors Charles Sammis, Charles Weber and Simeon Katz for their valuable suggestions and comments. I also thank Professors and Drs. Thomas Henyey, Terence Longdon, Niren Biswas, Peter Leary, Raul Madariaga, Peter Malin, Egill Hauksson, Kojiro Irikura, Haro Sato, Rushan Wu, Steve Lund, Yuantai Chen, Guoming Xu, Longshen Gao, Francisco Sanchez-Sesma and Kiyoshi Yomogida, from whom I broadened my horizon of my scientific understanding and benefited from discussions with them. I am also grateful to John McRaney for his assistance during our stay at USC. I will never forget my years passed at USC, all the happiness and difficulties, the excitement and dry spells. I am grateful to all my fellow students and scholars at USC who made my stay enjoyable. Stuart Koyanagi who proofread many of my thesis chapters and manuscripts deserves special acknowledgement. Thanks for Xiaofei Chen with whom I enjoyed our valuable discussions and friendship, for Tianqing Cao who helped me and my wife to get through our first registration at USC, for John Faulkner from whom I obtained the surface wave Gaussian Beam programs and started my first research project, for Valerie Ferrazzini who made us wonderful coffee in my final defense, for James Chin who provided my thesis figure 5.1, for my officemate Qiang Qu from whom I learnt more about computer, and for Kevin Mayeda, Shabei Li, Yonggang Li, Anshu Jin, Sandra Steacy, Geoffrey Saldivar, Michelle Roberson, Rafael Benites, Yin An, Joe Jackson and Calvin Lee of their friendship I will always remember. Thanks are extended to Hung-Chie Chiu, Jin Wang, Ouyang Hongping, Ron Biegel, Sun Ping, Yuhuda Benzion, Jenn-Yih Peng, Syhhong Chang, David Adams and Kenichi Kato. I also thank our computer manager Mustag Khan for his assistance with the VAX and SUN computers Finally, I would like to express my deepest gratitude to my parents for their endless love and encouragement, to my wife, Feng Su, for her love, support, and keeping my work in perspective, and to my brothers and sisters for their understanding and their assistance whenever I needed them. Without them, I would not have finished all my thesis work yet. This research was supported by NSF grant EAR-8904964 and DOE grant DE- FG03-87ER13807. iii Table of Contents C h a p te r P a g e AC KN OW LED GEM ENT ....................................................................................................... ii A B S T R A C T ................................................................................................................................... xvi I INTRODUCTION 1.1 Scope and Objective .................................................................... 1 1.2 R e v ie w .............................................................................................. 5 1.3 O utline of the Thesis .................................................................... 9 II REVIEW OF THE ISOCHRON METHOD FOR FORWARD SEISMOGRAM SYNTHESIS 2.1 Basic Assumptions and Equations ........................................... 13 2.2 A sym ptotic G reen’s Function .................................................... 19 2.3 Ray Tracing in 3-D Vertically Heterogeneous M edia................. 24 2.4 Other Numerical Treatments ..................................................... 26 2.5 N um erical Applications .............................................................. 30 2.5.1 Earthquake Nucleation and Stopping with a Hypothetic Iso ch ro n Velocity .......................................................................... 30 2.5.2 Comparison with Spudich and Frazer (1984) .............................. 34 2.5.3 Comparison with Discrete Wave Number Accelerograms 39 III RECURSIVE STOCHASTIC INVERSE ALGORITHM 3.1 Introduction ................................................................................... 44 iv 3.2 Recursive Stochastic Inverse Algorithm .................................... 45 3 .3 Model Covariance Matrix RA m o .................................................... 50 3.4 N oise C ovariance M atrix Rn ...................................................... 53 3.5 Applications to Seismic Inverse Problems ................................ 54 3.5.1 Seism ic T om ographic Problem .................................................. 54 3.5.2 Source Imaging Using the Isochron Method .............................. 57 3.5.3 Inversion of Site Effect Using Coda Power Spectra................... 61 3.6 D iscussion and Conclusion ......................................................... 63 IV APPLICATION TO THE WHITTIER-NARROWS EARTHQUAKE 4.1 D escription of the Problem ......................................................... 65 4.2 D ata Collection and Modeling .................................................... 65 4.3 Inversion Result and Analysis .................................................... 72 4.4 D iscussion and Conclusion ......................................................... 97 V APPLICATION TO THE LOMA PRIETA EARTHQUAKE 5.1 D escription of the Problem ......................................................... 103 5.2 Data Collection and Modeling .................................................. 103 5.3 Inverse R esult and Analysis ........................................................ I l l 5.4 D iscussion and Conclusion ...........................................................124 VI MAPPING OF THE HIGH FREQUENCY SOURCE RADIATION 6.1 I n tr o d u c tio n ...................................................................................... 130 6.2 B asic T h e o r y .....................................................................................131 6.3 Application to the Loma Prieta Earthquake ............................... 133 6.3.1 Data Collection and Modeling ................................................ 134 V 6.3.2 R esults and Analysis ......................................................................135 6.3.3 D iscussion ...................................................................................... 144 6.4 C o n c lu sio n ........................................................................................146 VII HIGH FREQUENCY SEISMIC WAVE SCATTERING 7.1 Introduction ..................................................................................... 148 7.2 B asic Theory ................................................................................... 149 7.3 Iterative Solution of the Scattered Wave Energy Equation 151 F irst-o rd er Scattering .................................................................. 153 S econd-order Scattering ............................................................... 155 H ig h -O rd er Scattering ................................................................. 156 7.4 Integral Solution of the Scattered Wave Energy Equation............ 158 7.5 E nergy Conservation ................................................................... 168 7 .6 Scattered Wave Energy Equation in a Non-uniform Random M edium ........................................................................................... 171 7.6.1 B asic Form ulations ...................................................................... 171 7.6.2 N u m erical Solutions .................................................................... 172 7.7 D iscussion and Conclusion ........................................................ 178 Appendix: First, Second and Third-order Scattering F irst-o rd er Scattering .................................................................. 180 S econd-order Scattering ............................................................... 181 T h ird -o rd er S c a tte rin g ....................................................................182 VIII SUMMARY AND CONCLUSION ......................................................... 185 RE FE R E N C E .......................................................................................................... 191 vi List of Figures F igu re Pa g-£ Chapter II 2.1 Geometry of source and receiver. X is the fault plane, xoy-plane is located on the free surface with z axis point downward, x is the observation point. Ray shot from x down to the fault plane at £,............. 13 2.2 Area element dX defined in terms of element isochrons dl and isochron velocity times time difference between two isochrons. q is a unit vector perpendicular to the isochron .................................................................... 16 2.3 Geometry of layered structure model with ray starts from the free surface down to the layer between n and (n+l)th interfaces ............................... 24 2.4a Triangular grid division used for the linear interpolation.......................... 27 2.4b An exam ple o f the triangular cell .............................................................. 27 2.5a Shear wave isochrons map for a constant shear wave isochron velocity c = 2 . 6 km /sec ................................................................................................ 30 2.5b Synthetic displacement and velocity seismograms calculated from the isochron model shown in figure 2.5a with constant slip over the fault 31 2.5c Synthetic displacement displacement seismograms calculated using equation (2 .1 0) and integrated from the result calculated using equation (2.11) for the isochron model shown in figure 2.5a with constant slip o v er the fault ............................................................................................... 32 2.6a Contours of the rupture time function interpolated from the digitized Spudich and Frazer's model with contour interval 0.05 sec .................... 34 2.6b Slip density models interpolated from the digitized Spudich and Frazer's vii model 34 2.7 a Geometry of the fault and receiver, and velocity structure adopted from Spudich and Frazer. The fault extends from 5 to 6 km in the vertical direction and -0.5 to 2.5 km in the horizontal direction.......................... 35 2.7b Shear wave isochrons at the receiver shown in figure 2.7a for the rupture in fig u re 2.6a 36 2.8a Velocity seismogram of x-component for the receiver shown in figure 2.7a. Results of our calculations (left) are compared with Spudich and F ra z e r (right) ................................................................................................ 37 2.8b Acceleration seismogram of x-component for the receiver shown in figure 2.7a. Results of our calculations (left) are compared with Spudich and F ra z e r (right) ................................................................................................ 37 2.9 Fault and receiver geometry adopted from Bouchon's (1982) rectangular dislocation model for the Gilroy 6 recording of Coyote lake earthquake of 1979 in California. The shear wave velocity is 2.6 km/sec. in the upper layer and 2.8 km/sec. in the lower half space ....................................... 40 2.10 Shear wave isochron for the model and receiver shown in figure 2 .9 ......... 40 2.11 The synthetic displacement seismogram from our calculation (left), and that from Bouchon (lower right) and Bernard and Madariaga (upper right) .. 41 2.12 The synthetic velocity seismogram from our calculation (left) and that from Bernard and Madariaga (right) ...................................................... 42 Chapter III 3.1 Exemplary segmentation of i-th ray passing through the square-area which is divided into N=m2 number of square sub-grids with each having a c o n sta n t velocity ......................................................................................... 54 3.2 (a) Map view of the fault trace (AA') and station distributions (triangles), (b) geometry of the fault Cross section along the fault strike, and (c) v i i i depth dependent P-wave and S-wave velocity structure of the m odel 57 3.3a Slip velocity density contours with value of the most inner line equals 1.0 m/sec and the interval value of any adjacent lines equals 0.1 m /sec 58 3.3b Fault rupture time contours with the starting point near the up-left comer and 0.025 sec contour interval ................................................................. 58 3.4a Inverse solution of the fault rupture time distribution obtained using C hebyshev sem i-iterative method ........................................................... 59 3.4b Inverse solution of the fault rupture time distribution obtained using our recursive stochastic inverse process ....................................................... 59 Chapter IV 4.1 Map view of the strong motion station distribution used for the Whittier- Narrows earthquake source inversion study .......................................... 67 4.2 Seismic velocity model for the Los Angeles basin a re a ......................... 69 4.3 Contours of rupture time model for the Whittier-Narrows earthquake with constant rupture velocity of 2.5 km/sec and contour interval of 0.5 se c 70 4.4a The fault slip (cm) distribution of the first inversion shaded in different intensity according to scales indicated below the fig u re .......................... 72 4.4b Contour plot of the same slip distribution as in Figure 4.4a with contour in terv al o f 15.21 cm ................................................................................. 72 4.5 Plot of seismograms for the first inversion with the observed traces on top of the synthetic ones for each station numbered along the vertical axis. Each seismogram is normalized to its peak values (cm), indicated to the right the trace followed by its epicentral distance (km). Seismograms to the left are components parallel to the fault strike and those to the right are perpendicular to the fault strike ...................................................... 73 4.6a The fault slip (cm) distribution of the second inversion shaded in different ix intensity according to scales indicated below the figure 76 4.6b Contour plot of the same slip distribution as in Figure 4.6a with contour in terv al of 17.93 cm ................................................................................... 76 4.7 Same as Figure 4.5 but for the second inversion resu lts.......................... 77 4.8a The fault slip (cm) distribution of the third inversion shaded in different intensity according to scales indicated below the fig u re......................... 80 4.8 b Contour plot of the same slip distribution as in Figure 4.8a with contour in terv al o f 19.50 cm ................................................................................... 80 4.9 Same as Figure 4.5 but for the third inversion results .............................. 81 4.10 Contours of rupture time model determined by the simultaneous inversion (forth) of both slip and rupture time from the Whittier-Narrows earthquake strong motion records with contour interval of 0.5 sec .......................... 84 4.1 la The fault slip (cm) distribution of the forth inversion shaded in different intensity according to scales indicated below the fig u re......................... 85 4.11 b Contour plot of the same slip distribution as in Figure 4.1 la with contour in te rv al of 17.50 cm ................................................................................... 85 4.12 Same as Figure 4.5 but for the forth inversion results ............................ 8 6 4 .13a The fault slip (cm) distribution of the fifth inversion shaded in different intensity according to scales indicated below the fig u re .......................... 90 4.13b Contour plot of the same slip distribution as in Figure 4.13a with contour in te rv al of 1 8 . 6 6 cm ................................................................................... 90 4.14 Same as Figure 4.5 but for the fifth inversion results ............................. 91 4 .15a The fault slip (cm) distribution of the sixth inversion shaded in different intensity according to scales indicated below the fig u re .......................... 93 X 4.15b Contour plot of the same slip distribution as in Figure 4.15a with contour in te rv al o f 22.91 cm ................................................................................... 93 4.16 Same as Figure 4.5 but for the sixth inversion results 94 4.17 A rep lo t of F igure 4.8 99 4.18 (a) Map, (b) south-north cross section, and (c) front view from the south of the mainshock focal mechanism and following aftershocks unitil 10:58, October 4, 1987, which was the time of the largest aftershock (reproduced Chapter V 5.1 Three component synthetic ground velocities (cm/sec) produced by a dislation point source at 12 km depth with the focal mechanism as that of the 1989 Loma Prieta earthquake and a seismic moment of 1020 dyne-cm. The source time function is a ramp function with 0.5 sec rise time. The hypocentral distance inside the parenthesis is in km. The reduction v elo city is 3.5 km/sec ............................................................................... 103 5.2 Map view of the strong motion station distribution used for the Loma Prieta earthquake source inversion study ................................................ 106 5.3 Seismic velocity model for the Santa Cruz mountain a re a ........................ 108 5.4 Intial fault slip distribution shaded in different intensity with a high central plateau that gradually tapers to zero at the fault ed g es............................ 109 5.5 Contours of rupture time model for the Loma Prieta earthquake with constant rupture velocity of 80 percent of the average local S-wave velocity and contour interval of 0.5 sec ................................................ I l l 5.6a The fault slip (cm) distribution of the first inversion shaded in different intensity according to scales indicated below the fig u re............................1 1 2 5.6b Contour plot of the same slip distribution as in Figure 5.6a with contour from Hauksson and Jones, 1989) 100 xi interval of 45.80 cm 112 5.7 Plot of seismograms for the first inversion with the observed traces on top of the synthetic ones for each station numbered along the vertical axis. Each seismogram is normalized to its peak values (cm), indicated to the right the trace followed by its epicentral distance (km). Seismograms to the left are components parallel to the fault strike and those to the right are perpendicular to the fault strike ................................................................ 113 5 .8 Contours of rupture time model for the Loma Prieta earthquake with constant rupture velocity of 70 percent of the average local S-wave velocity and contour interval of 0.5 sec ................................................ 115 5.9a The fault slip (cm) distribution of the second inversion shaded in different intensity according to scales indicated below the figure ................. 116 5.9b Contour plot of the same slip distribution as in Figure 5.9a with contour in terv al of 5 9 .4 0 cm .....................................................................................116 5.10 Same as Figure 5.7 but for the second inversion results ................. 117 5.11a The fault slip (cm) distribution of the third inversion shaded in different intensity according to scales indicated below the figure ................. 1 2 0 5.11b Contour plot of the same slip distribution as in Figure 5.1 la with contour in terv al o f 7 6 .3 2 cm .....................................................................................120 5.12 Same as Figure 5.7 but for the third inversion results ..............................121 5.13 Overlain of the projection of well-located aftershocks (U. S. Geological Survey staff, 1990) as well as the large early aftershocks (Simila, et a i, 1990) onto the fault slip contours of the first inversion ........................... 125 5.14 Same as Figure 5.13 but for the slip contours of the second inversion 126 5.15 Same as Figure 5.13 but for the slip contours of the third inversion 127 xii Chapter VI 6.1a High frequency energy intensity (km-cm) distribution of the first inver sion shaded in different intensity according to scales indicated below the figure .............................................................................................................. 135 6.1b Contour plot of the same energy intensity distribution as in Figure 6 . la with contour interval of 1 .0 0 (km-cm) .................................................. 135 6.2 Plot of envelope seismograms for the first inversion with the observed traces on top of the synthetic ones for each station numbered along the vertical axis. Each seismogram is normalized to its peak values (cm), indicated to the right the trace followed by its epicentral distance (km ) 136 6.3 a High frequency energy intensity (km-cm) distribution of the second inversion shaded in different intensity according to scales indicated below th e figure ....................................................................................................... 138 6.3b Contour plot of the same energy intensity distribution as in Figure 6.3a with contour interval of 1.60 (km-cm) ............................................................. 138 6.4 Same as Figure 6.2 but for the second inversion results ........................ 139 6.5 a High frequency energy intensity (km-cm) distribution of the third inver sion shaded in different intensity according to scales indicated below the figure .............................................................................................................. 141 6.5b Contour plot of the same energy intensity distribution as in Figure 6.5a with contour interval of 2.21 (km-cm) ..................................................... 141 6 . 6 Same as Figure 6.2 but for the third inversion results .............................142 Chapter VII 7.1 The Geometry for positions of the source, receiver and scatterers where S represents source, R represents receiver and numbers represent scatterers. (a) First-order scattering case, (b) Second-order scattering case, (c) Third-order scattering case ....................................................... 153 7.2 The contour of integration in complex k plane and location of branch point a ti(rj+ s/v ) .................................................................................................... 160 7.3a Scattered wave energy decay curves with the absorption coefficient r\f= 0 .0 1 km-1, the scattering coefficients r|s=0 .0 0 2 km-1 and source-receiver distance r=10 km. Solid curve represents our multiple scattered wave energy solution which is compared with the results of single scattering, energy-flux and diffusion models shown in dashed, dot-dashed and dot lin es, respectively ....................................................................................... 164 7.3b Same as 7.3a but for the scattering coefficients T|s=0.05 kirr1 .................. 164 7.4a Same as 7.3a but for the scattering coefficients rjs=0.01 km- 1 and the absorption coefficients tj ^ 0 . 0 0 2 k irr 1 .................................................... 166 7.4b Same as 7.4a but for the absorption coefficients ^=0.05 km-1 ................ 166 7.5 The Geometry for positions of the source and scatterers where S represents the source, numbers represent scatterers and the dashed circle represents the primary wave front, (a) First-order scattering case, (b) n* -order sc a tte rin g case ............................................................................................ 168 7.6 Scattered wave energy decay curves obtained at different source-receiver distances in the medium with absorption coefficient T j^ 0 .0 1 k m 1 and scattering coefficient rjs=0.04 km ' 1 ..........................................................172 7.7 S-wave Q' 1 and S-coda Q ' 1 vs frequency observed from different regions o v e r the world ............................................................................................. 173 7.8 Numerical solutions of scattered wave energy decay curves calculated using equation (7.25a) with the same absorption and scattering coefficients as fo r F ig u re 7.7 175 7.9 Same as Figure 7.8 but for a variable absorption coefficient that decreases with depth at the rate of 0.00018 km' 1 per km from a free surface value o f 0 .0 2 k m " 1 ................................................................................................ 176 xiv List of Tables Table Page 4.1 Strong motion station names and locations ............................................. 6 6 4.2 Site am p lificatio n factors ............................................................................ 79 4.3 Seismic moment of the inversion results ................................................ 8 8 5.1 Strong motion station names and locations .............................................. 105 5.2 Site am p lificatio n factors ............................................................................. 119 5.3 Seismic moment of the inversion results ................................................ 122 6.1 Site am p lificatio n factors ............................................................................. 140 7.1 Comparison of scattering energy coefficients ........................................... 157 xv Abstract The present thesis takes deterministic and stochastic modeling approaches for understanding earthquake source, station site and propagation path effects on short- period seismic waves, with special emphasis on earthquake source imaging using strong motion data and modeling of seismic scattering in the lithosphere. The deterministic studies of earthquake source processes have been done under the framework of the isochron synthetic method and a fast recursive stochastic inversion algorithm. We derived analytical expressions for the asymptotic Green's function and solutions for ray tracing in a vertically inhomogeneous medium. Using strong motion records, we determined the fault slip and rupture time for the 1987 Whittier-Narrows earthquake and the fault slip with different choices of constant rupture velocity for the 1989 Loma Prieta earthquake. We also included the site effect correction in our source inversions using the site amplification factors derived from the coda part of the strong motion records. Both earthquake inversion results indicate that slip is highly variable over the rupture plane and aftershocks are likely to occur at places deficient in mainshock slip. The above deterministic modeling is applicable to a relatively lower frequency range of the strong motion data (0.6-4.5 Hz). By assuming that the slip responsible for high frequency wave (> 5 Hz) generation is incoherent between neighboring patches on the fault, we formulated the isochron integration equation for the high frequency wave radiation from the fault. Using the recursive stochastic inversion scheme, we determined the high frequency energy radiation intensity for the 1989 Loma Prieta xv i earthquake from the mean squared strong motion seismograms. Our results indicate that the high frequency energy sources are located along or near the boundaries of large slip zones, which is consistent with the theoretical considerations that high frequencies are primarily generated from the rupture stopping or areas with large slip variation. Next, we studied the problem of the scattered wave energy propagation from a point source in a random isotropic scattering medium. We developed a scattered wave energy equation that unifies both the multiple ray-paths and energy transport approaches. This equation was formulated by extending the stationary energy transport theory to the time dependent case. Our numerical results indicated that single scattering and energy-flux models agree well with the multiple scattering model in weak scattering. For strong scattering case, diffusion model provides a good approximation. By generalizing our scattering wave energy equation to an inhomogeneous random scattering medium, we obtained numerical solutions of the equation in a medium with vertical increase in absorption that explains the observed coincidence in attenuation factor Q between direct S-waves and S-coda waves. Chapter I INTRODUCTION 1.1 Scope and Objective The subject of the present thesis is the generation and propagation of high- frequency seismic waves. In principle, seismic waves in an arbitrarily complicated earth model generated by an arbitrarily complex source model can be simulated by the finite difference or finite element method. We cannot, however, study high-frequency waves even with the most powerful supercomputer because of limitation in the memory size and speed. In the present thesis, we shall explore two approaches to deal with the generation and propagation of high-frequency seismic waves. One is the asymptotic method for delineating the smoother part of the source and medium deterministically. The other is the kinematic method for studying the energetics of wave generation and propagation by stochastic modelling. Let us first focus our attention on the seismic source study. In the early days of seismology, it was assumed that the earthquake was a point source. A finite rupture source model was not used to explain the data until the great Chilean earthquake of 1960. Benioff, Press and Smith (1961) were the first to apply a rupture propagation model to the Chilean earthquake of May 22, 1960. They found that the observed phase shift between the vertical and horizontal components of the earth's free oscillation could be explained by a finite rupture source. Now it is well recognized that earthquake faulting is a finite rupture process and moreover that earthquake faulting is a very inhomogeneous process with slip and rupture time highly variable over the fault. 1 In our deterministic source modelling, we use the asymptotic approach based on the high frequency approximation. The geometrical ray method is used for the wave propagation effect. The source radiation processes are also simplified using the asymptotic approach. Madariaga (1983) has studied the wave radiation from a dynamic earthquake fault model and found that the high frequency far field radiation comes from the rupture front. This result has been supported by Spudich and Cranswick (1984) from their observations of the 1979 Imperial Valley, California earthquake. Based on these studies, Bernard and Madariaga (1984) and Spudich and Frazer (1984) have developed independently the so-called isochron method to calculate synthetic seismograms from a complicated fault rupture model. For our theoretical preparation, we are going to follow their initial idea and re-derive these isochron formulations. Following the forward synthetic calculation, the next step is the inversion. Beroza and Spudich (1988) used the perturbation method to obtain a linearized inverse problem and calculated the solution of the problem using the Chebyshev semi-iterative method. This method provides a fast inversion by taking advantage of the sparseness in the matrix equations. However, the iterative solution usually does not give error and resolution estimates for model parameters and sometimes may suffer from convergence problems. As an alternative, we present another inversion method called recursive stochastic inversion. The recursive formula gives the same stochastic inverse obtained by the exact matrix inversion as well as the error and resolution matrices of the solution. Since the recursive formula can take advantage of the the sparseness of the matrix in the same way as the tomographic iterative method (Clayton and Comer, 1983), it leads to a fast and efficient inverse scheme for certain kinds of large sparse matrix equations that commonly arise in seismic tomographic problems, including the inversion for fault slip and rupture time distribution using the isochron method. 2 We shall apply our forward and inverse method to the Whittier-Narrows earthquake which occurred on October 1, 1987. The abundant strong motion data from the Whittier-Narrows earthquake provides a unique opportunity to study earthquake source, site and path effect. We will perform the inversion to determine fault slip distribution and rupture front propagation on the fault. We will also study the site effect on the waveform inversion solutions and try to minimize the propagation path effect in our inversion process. The results should provide us with a better understanding of the earthquake source and site effect. Of more importance is how this inversion result is correlated with other studies of earthquake faults, such as the aftershock distribution, and its implication for fault physics. We find, for example, that aftershocks occur in areas of low mainshock slip. Another interesting example is provided by the Loma Prieta earthquake. This earthquake of magnitude 7.1 occurred on October 17, 1989. It is the most significant event to strike the San Francisco Bay area since the 1906 San Francisco earthquake. The mainshock is well recorded by many strong motion stations. A surprise to everybody is that while strike slip faulting is observed for most earthquakes occurring along the San Andreas, this earthquake had a unexpectedly large dip slip component. Understanding the nature and extent of the faulting of the Loma Prieta earthquake is of great importance in studying the San Andreas fault activity and prediction of the strong motion for future earthquakes. This earthquake is recorded on many sites of various geology, thus providing seismologists with a better example for studying site effects on strong ground motion. We will also compare the mainshock slip distribution with the aftershock activity to further check the relation between the mainshock slip deficiency and the aftershock activity which is found for the Whittier-Narrows earthquake. 3 In order to study the high frequency source radiation, we will develop a new model of the earthquake fault which combines both deterministic and stochastic processes. By assuming that the slip responsible for high frequency wave generation is incoherent between neighboring patches on the fault, we can superpose the high frequency wave energies radiated from different patches on the fault at the receiver. Using the same isochron method we formulate the high frequency wave energy radiation equation. This equation is then applied to the Loma Prieta earthquake to determine the high frequency wave energy radiation intensity from the earthquake fault. We will compare this high frequency wave energy inversion result with the low frequency fault slip inversion solution and discuss physical implications. Next, we look at the high frequency scattered wave energy propagation in an isotropic scattering medium. In the wave scattering problem, the Bom approximation is usually made to obtain the first order solution. An important issue here is to what degree this approximation is adequate. To answer this question, we need to study the multiple scattering problem and compare the result with the single scattering solution. To obtain a complete description of the scattered wave energy propagation in a random isotropic medium, we will extend the stationary energy transport theory studied by Wu (1985) to the time dependent case and formulate an integral equation for the scattered wave energy propagation in a random isotropic scattering medium. By solving this equation, we should obtain the multiple scattering solution and compare the results with other studies, such as the single scattering, energy-flux and diffusion models. We will also be able to answer questions regarding the validity of the Bom approximation, and interpretations of the observations. 4 1.2 Review Strong motion seismology was initiated by earthquake engineers who developed and operated seismographs to record strong ground motion. The modem quantitative analysis of near source strong motion observations started from the publication of the now-famous station No. 2 record obtained from the 1966 Parkfield earthquake at a distance of only 80 meters from the fault break (Housner and Trifunac, 1967). This modest start has developed into an active field of earthquake source research during the past two decades. It was shown by Aki (1968) and Haskell (1969) using the unbounded homogeneous elastic medium model that the observed transverse component motion at station No. 2 (the parallel component of the instrument was not working properly) of the 1966 Parkfield earthquake was precisely what is expected for a right-lateral strike slip propagation from northwest to southeast. This simple model was later improved with better constraints from more observations. The first refinement was to divide the rupture plane into subfaults with each assigned a different slip vector (Trifunac and Udwadia, 1974). Nonuniform slip distribution over the rupture surface has been found to better explain observations for the Borrego Mountain earthquake of 1968 by Heaton and Helmberger (1979), and Shakal (1979), and for Imperial Valley earthquake by Hartzell and Helmberger (1982), and Olson (1982). A variable slip distribution over the fault plane has also been found from the analysis of teleseismic body wave data (e.g. Ebel and Helmberger, 1982; Langston, 1978; Kanamori and Stewart, 1976; Wallace et al., 1981). Waveform inversion methods were well developed for the study of the 1979 Imperial Valley, California earthquake. The extensive data set generated by the 1979 Imperial Valley earthquake has provided a unique opportunity for researchers to fully 5 explore the dynamic characteristics and the finiteness of the fault plane. It was found that the radiation responsible for strong ground motion mainly came from contributions from the deeper portion of the fault plane (Olson and Apsel, 1982; Archuleta, 1982, 1984; Hartzell and Helmberger, 1982; Hartzell and Heaton, 1983, 1986). By allowing the slip to occur more than once as the rupture front passes, Olson and Apsel (1982) used least squares inversion to obtain the slip on subfaults from the strong motion data. Using the same subfault division method, Hartzell and Heaton studied the earthquake source with a joint inversion of strong motion data and the teleseismic records. Archuleta (1984) derived a very detailed model of rupture for the Imperial Valley earthquake by a forward calculation to fit the velocity waveforms, resulting in a model with variable slip and rise time. Linearized inversions have been carried out recently by Fukuyama and Irikura (1986) and Takeo (1987) using nonlinear parameterization to solve both the slip and rupture time distribution on the fault for the Akita-oki and Naganoken-Seibu earthquakes. Olson and Anderson (1988) have proposed a frequency domain inversion method. By transforming the temporal convolution into a spatial multiplication, this method allows for simultaneous linear inversion of both fault slip and rupture time distribution on the fault. The method was later successfully applied to the study of the temporal and spatial evolution of the Michoacan earthquake September 19, 1988 (Mendez and Anderson, 1991). Spudich and Frazer (1984) have formulated an interesting and promising approach now known as the isochron method for short-period near-source radiations calculation. Independently, Bernard and Madariaga (1984) have made the same formulation and both have presented results of forward calculations that compare favorably with those derived from more elaborate numerical schemes such as finite difference, finite element, and discrete wavenumber methods. Later, Beroza and Spudich (1988) succeeded in 6 applying the isochron method to the strong motion data of the 1986 Morgan Hill earthquake, California. From the inversion results, they were able to postulate the possible existence of an asperity on the fault plane that may have impeded the rupture propagation, and to give an estimate of the earthquake fault's shear fracture energy. The waveform inversion studies reviewed above may be categorized as deterministic modelling. These methods provide us with direct physical descriptions of earthquake faulting, but their applications are limited to relatively low frequency waves which are coherent between stations. At high frequencies, the deterministic modelling becomes difficult to apply due to many unknown small scale inhomogeneities at the source, site and wave path. The stochastic models avoid this difficulty by characterizing the statistical properties of the source process using a small number of parameters (Blanford, 1975; Boore and Joyner, 1978; Hanks, 1979; Andrews, 1980, 1981; Papageorgiou and Aki, 1983a). An example of such a model is the specific barrier model of Papageorgiou and Aki (1983) constructed for a quantitative description of the earthquake faulting process. Papageorgiou and Aki (1983b) applied the model to the power spectrum of the high frequency accelerograms to estimate the source parameters for some major earthquakes in California. Let us now turn to the review of studies on high frequency seismic wave scattering and attenuation. One of the most prominent features of local earthquake seismograms is their long decaying train of waves after the direct S-wave arrival which puzzled seismologists for decades in the early stages of instrumentation seismology. It was not until the 1960's, aided by developments in seismic instrumentation, that seismologists started to solve the puzzle regarding this train of waves known as coda waves. Wesley (1965) first used the diffusion theory to explain these waves. But an acceptable 7 physical modelling was not obtained until Aki (1969), who recognized that the coda wave can be modelled as back-scattered waves from more or less uniformly distributed scatterers in the earth's lithosphere. By treating the earth as a random isotropic scattering medium, Aki (1969), Aki and Chouet (1975), and Sato (1977) studied the coda waves as single back-scattered waves. This single scattering model has been widely used to interpret observations and to calculate coda Q*1. To complement the single scattering wave description of the seismic coda, various attempts were made to investigate the multiple scattering process. For examples, Kopnichev (1977) and Gao et al (1983a, 1983b) studied the 2-D and 3-D multiple isotropic scattering problem by adding the higher order scattered wave contributions to the single scattering effect. Wu (1985), on the other hand, introduced the stationary energy transport theory for the multiple scattering problem. This theory can be used to separate the effects of scattering and intrinsic attenuation uniquely as demonstrated by Wu and Aki (1988) for the Hindu Kush region. Later, Frankel and Wennerberg (1987) proposed a 2-D and 3- D energy flux model for coda waves based on the finite difference simulation results obtained by Frankel and Clayton (1984, 1986) for 2-D random media. In their model, they assumed that the scattered wave energy is uniformly distributed immediately behind the direct wave front. More recently Hoshiba (1990) used a Monte-Carlo simulation to study the multiple isotropic scattering problem by incorporating the energy conservation law. By considering scattered wave energy propagation in a random isotropic scattering medium, Zeng et al. (1991) have shown that both the traditional single and multiple scattering formulas and the energy transport formulas can be unified in an integral equation. This integral equation is called the wave energy equation obtained by extending the stationary transport theory studied by Wu (1985) to the time dependent case. By solving this equation iteratively, they obtained a power 8 series in the scattering coefficient which corresponds to the exact order of the scattered waves discussed by Gao et al. (1983). By taking the Fourier transform of the equation, they found a compact closed form integral solution for the same problem, which is easy to evaluate numerically. The scattered wave energy equation was also generalized to a nonuniform random medium adequate for a more realistic earth model. 1.3 Outline of the Thesis This thesis uses the deterministic and stochastic modelling approaches to study earthquake source, station site and lithospheric wave propagation path effect on short- period seismic waves. The main body of the thesis is divided into six chapters. In Chapter II, we re-formulate all the isochron equations for the synthetic seismogram calculation based on a given fault rupture model. The analytical expressions for the asymptotic Green's function and for ray tracing in a vertically inhomogeneous medium are derived. We also include a detailed numerical procedure necessary for computing isochron integration. The formulations are applied to the synthetic seismogram calculations and the results are compared with the published ones. In Chapter III, a recursive stochastic inversion process is derived. It is shown that the recursive formula gives the same stochastic inverse solution obtained by the exact matrix inversion as well as the error and resolution matrices of the solution. We demonstrate that the recursive inverse method can take advantage of the sparseness of the matrix equations, thus leading to a fast and efficient inversion scheme. In this chapter, we will also discuss some useful techniques regarding the choice of the initial model covariance matrix and the noise covariance matrix. 9 In Chapter IV and Chapter V, the methods described in Chapter II and III are applied to the 1987 Whittier-Narrows earthquake and 1989 Loma Prieta earthquake, respectively. Extensive strong motion data are available for both earthquakes. Based on these strong motion data, we determine fault slip and rupture time for the Whittier- Narrows earthquake and fault slip distribution with different choices of constant rupture velocity for the Loma Prieta earthquake. To incorporate the site effect into our source inversion study, we determine the site amplification factors following Phillips and Aki (1986) using the coda part of the strong motion records. Based on these site amplification factors, we then correct strong motion data and perform the inversion to determine the earthquake source parameters. The site effect corrected results are compared to the ones without site correction and their differences are discussed. We then compare the mainshock slip distribution with the aftershock activity patterns. We find that the aftershocks tend to occur at places deficient in mainshock slip. This observation can be explained well by the barrier model (Aki, 1979) of fault heterogeneity. In Chapter VI, we developed a new model of the earthquake source based on a combination of the deterministic and stochastic processes. We assume that the slip, which is responsible for the high frequency generation, is incoherent between neighboring patches on the fault. According to the isochron method, we formulate the high frequency wave energy radiation equation. Using the recursive stochastic inversion process, we can determine the high frequency energy radiation intensity (see Chapter VI for definition) from the observed mean square seismograms. The method is then applied to the 1989 Loma Prieta earthquake. We performed inversions with and without site and path attenuation corrections. The inversion results of the high 10 frequency energy radiation intensity are analyzed and compared with the mainshock slip distribution obtained in Chapter V. In Chapter VII, we present a stochastic model for the scattering and attenuation of high frequency seismic waves. This model is applicable to a medium with randomly distributed isotropic scatterers. First, we formulate the scattered wave energy equation by extending the stationary energy transport theory studied by Wu (1985) to the time dependent case. Then we present two different expressions for the solution of the energy integral equation and discuss their physical implications and numerical implementation procedures. Examples of the solution are presented and compared with the single scattering, energy flux and diffusion models. We then discuss the energy conservation for our system by starting with our fundamental scattered wave energy equation and then show that our formulas satisfy the energy conservation law when contributions from all orders of scattering are summed. We also generalize our scattered wave energy equation to the case of nonuniformly distributed isotropic scattering and absorption coefficients. Using a discretization scheme, we numerically solve this equation in a vertically inhomogeneous random isotropic scattering medium and explain the observed coincidence in attenuation factor Q between direct S-waves and S-coda waves. In Chapter VIII, we summarize major conclusions from the preceding six chapters and propose future research directions for further study on the subject. 11 Chapter II REVIEW OF THE ISOCHRON METHOD FOR FORWARD SEISMOGRAM SYNTHESIS 2.1 Basic Assumptions and Equations We start with the representation theorem for the elastodynamic field due to discontinuity in displacement across an internal surface (Aki and Richards, 1980, p. 39). The displacement un(x,t) due to the slip discontinuity Aui(x,t) across fault plane X (see Figure 2.1) is given by where n is the normal of X, Cyki are the elastic coefficients of the medium, and is 3Gnk(x,t-x;^,0) --------------------- d2> (2. 1) the k1 * 1 component of displacement Green’ s function due to an impulsive point source acting on x in the xn-direction. In an isotropic medium, the traction due to the Green's function is Tn(x,t; 4,0) = CypqrijoG a^ x ,t-v ^ 0 ) = ^ (V |. Gn)n+|x[v 4 G n - Tn(x,t, ^,0 ) — Cijpqltj - ^(V^*Gn )n+p[V^Gn +(V^Gn ) ]*n Then, the n^ component of displacement can be expressed by (2.2) or in the frequency domain, the above expression becomes 12 Figure 2.1. Geometry of source and receiver. X is the fault plane, xoy-plane is located on the free surface with z axis point downward, x is the observation point. Ray shot from x down to the fault plane at £. 13 un(x,co)=J Au (5,co)Tn(x,co;^,0) dL (2.3) For high frequencies, following Spudich and Frazer (1984), we assume Gn of the following form, G n=Unexp[-ico£*r(£)/v] (2.4) where v=a for P waves and v=p for S waves. By considering only the highest term in frequency, we have Tn on the fault at ^ as T „=^[X .(r*U n)n+|i(rU„+Unr)«n]exp[-iti4*r(^)/v] (2.5) We further assume that U£ exp[-io^*r(^)/a]=rA(x,^)exp[-icotp(x,^)] URexp[-io^*r(^)/p]=[0B(x,^)+(lC(x,^)]exp[-ioxs(x,^)] where indices p and s are used for P and S waves expresions, respectively, and A(x,£), B(x,^), and C(x,^) are the corresponding wave amplitudes for P-, SV- and SH-waves, respectively. Consider the case in which Au(£,t) _ L n, we then get TP=-icoGP(x,^)exp[-icotP(x^)] (2.6) Tn= -icoG^(x,^)exp[-icots (x, £)] where Gg=^(r.n)rA(x.5) 14 G>n=— [X i B (x,£)+x2 C(x,£)] p and x i = (0 • n )r+ (r* n )0 /S X2 = (<(>*n)r + (r«n) < |> Now going back to equation (2.2), we assume Au(^,t) can be simply written as In this expression, we have assumed that point £ ruptures at time tr(^) and heals at time t ^ ) . fr(t) and fh(t) are the time functions describing the rupture part and the healing part of the slip, respectively, and are assumed to be independent of position. Aur(^) and Auh(^) are the spatial distribution of the slip intensity due to the rupture and the healing part of the slip, respectively. For a dynamic crack model, fr[t-tr(£)] in (2.7) can be approximated as (t-tr) 1/2 H(t- tr), and H(t-tr) for uniform dislocation model, where H(t-tr) is the step function. In the following, for the sake of simplicity, we will consider equation (2.3) for S waves with the rupture term only. The same formulas apply to the healing part of the slip. Substituting (2.6) and (2.7) with rupture term only into (2.3) gives Au(5,t) = Aur (5)fr[t-tr £)]+A uh £ )fh [t-th (5)] (2.7) Aur(^)*Gn(x,^) icofr(C D ) exp[ico(ts +tr)] d£ Taking the inverse Fourier transform, we get 15 ta(x, £)= t 2 C d t Figure 2.2. Area element dX defined in terms of element isochrons dl and isochron velocity times time difference between two isochrons. q is a unit vector perpendicular to the isochron. 16 Un ( X , t ) = - f r (t)* Aur(£)*Gg(x, £)8(t-t|)dZ (2 .8) where tl^ + tj. Noticing that the surface integral in equation (2.8) contains a delta function, we immediately conclude that the contribution of the integral is concentrated along the line defined by tl(x,^)=t. Its physical implication is, as pointed out by Madariaga (1983), that the high frequency wave is emitted from the rupture front. Taking this into consideration, we can convert the surface integral in equation (2.8) to a line integral over the isochron ti(x, £)=t and rewrite (2 .8 ) as where dE=c(x, ^)dt dL (Figure 2.2), and c(x, ^)=|(I-nn) • Vtf(x, ^)| 1 is the 'velocity’ of the isochron curves along the fault surface normal to their length (Spudich and Frazer, 1984). Equation (2.9) is the basis for the so-called isochron method (Spudich and Frazer, 1984). c(x, ^) in the line integral is called the isochron velocity at point % on the fault. The ground velocity and acceleration are of the form (2.9) Un ( X , t ) = - f r (t)* Aur (£)*G £(x,£)c(x, ^)dL (2. 10) and Aur(^)»Gn(x, £)c(x, ^)dL 17 = - f r(t)* I { ^ i * G ^ c ^ A u r* ^ c ^ A u r*G ^ d c d q d q d t +[-d-V 0 ^ / a L ) 2+ O n /3 L )2 JA ur»G„ c2)dL (2.11) dq respectively, where q=cVs ta . Equations (2.9), (2.10) and (2.11) are basically the same as that given by Spudich and Frazer (1984) except the term [^-V (0^/0L)2+(0t|/0L )2]Aur*GnC2 in equation (2 . 1 1) which they have probably missed when they converted the derivative with respect to t to a derivative with respect to q and a negative sign omitted in their formulation. The above term actually gives a non-negligible contribution to the ground acceleration. For example, when we integrate it over a circular isochron, its value will equal 27t[Au(£;)*G£c2 ], where the upper bar denotes an average over the isochron. Equations (2.10) and (2.11) directly relate the ground velocity and acceleration with the isochron velocity and acceleration. Their physical interpretations have been discussed in detail by Spudich and Frazer (1984). 2.2 Asymptotic G reen’s Function In this section, we are going to derive the asymptotic Green's function in a homogeneous half space. Following Aki and Richards (1980, p. 304), we can express the elastic displacement in a homogeneous media by three scalar potentials < p , \\f and % of the form G=V<p+VxVx(0,0,\j/)+Vx(0,0,x) These potentials can be expressed as a superposition of the terms 18 <p=D (k,m,co)Jm(kr)e“T t l > _ V \|/=E(k,m,oo)Jm(kr)e1 I T t : ) ^'0 ' X=F(k,rn,co)Jm(kr)eiml > - ' 1 where y=Vk2-k« , d=V k2 -kn , ka = ^ and kB = ^ P a p If we define T £ (r,0) = ^ ^ f - - L ^ 0 to 3 0 k 3r S J'(r,0)=£— + J -----— 0 k 3r ^ 3 0 R£’(r,0)=-Y C'z where Y ^ J n ^ r ^ e 1 1 1 1 ! * , then we can obtain the expression for the Green's function displacement G and traction T in the z-plane as f + ° ° G(r,0,z,co)=Z I dk[kFe4 )z TJ1 (r,0)+k(DeJ V z-u E e^ S ™ (r,0) J — oo + ( J)e'yz -k 2 Ee_ 'uz )R jj1 (r, 0 ] (2.12) and r+ °° T(r,0,z,co)=X I dk{ - |i tkFe- ll2 T |ik [ 2'pe'7 Z + (k?- 2k2)EeA ) z ] S ™ ml J — oo 19 +^[(kp-2k2 )De''*z+2vk2E e w] }S^ (2.13) For an impulsive point traction acting on the surface at x, we can write its expression in the Fourier transform domain as j»+oo -F 5 fr )= £ I kdkftfk.mJlY+fsCk.mJSE1 +fI (k,m)R£1 ] mJ- o o From the orthogonality relations of TJJ1 , SJJ1 and R™(e.g. Aki and Richards, 1980, pp. 308-309), we find that ft (k, m)=fs(k, m)=0, except for m =l, and fr (k, m)=0, except for m=0. For those non-zero values, we have f,(k, l)= i-(F y +iFx ) f,(k ,-l) = -l(-F y+iFx ) fs(k, l) = -±(-Fx +iFy) fs(k, -1) =r-(Fx+iFy ) and fr(k, 0) = FZ. Solving the equation T | z=o= — FS(x), we can obtain the solution for D(k,m,w), E(k,m,co) and F(k,m,co). By substituting these solutions into equation (2.12) and carrying out the integral according to the Steepest Descents method, we will obtain the first order asymptotic Green function displacement. The solutions for P waves are U i=2 V (a 2 /p2 — sin-^) costjsinGUp? U£=2V (a 2/p2— sin20) sin<|)sin0Upr (2.14) 20 U f= (a 2/p 2-2sin^0)U[? , i cosGe lkaR where Up = — --------------------------------------------- 2;i|iR ( - a 2/p 2_ 2 sjn2 0^2+2 sjn2 0 sin 0 V ( a 2/p 2— sin2 © ) The solutions for S wave are / N / N Uf =cos2 ©cos<tUsv© — sin<j)Ush < j > /S /N U | =cos20sin<|jUsvG + costjjUsh^ (2.15) U 3 = —2 V (p2/ a 2— sin2 © ) sin0Usv e u 1 cos0 e lkpR 1 _ik R where Usv — — ----------------------------------------------------- and 1]^ = — -±-— e p . 27t|iR cos2 2 0 +2 sin2 0 sin0 V (p2/ a 2 — sin2 © ) % For an inhomogeneous half space, we can use the ray method to trace this asymptotic Green function from the receiver down to the source and obtain the Green function expression in equations (2.9), (2.10) and (2.11) at the fault. By identifying 1/R as the geometrical spreading factor 1/RP(^, x) for P waves and 1/R S (^, x) for S waves, and 1/(2 ti(I) as factor l/(27tp1/2(x)p1/2(^)P3^ 2(x)P1 ^ 2(^)), we can complete the formulae for P waves in an inhomogeneous half space: U p =2V (a^ xyp ^ x^ sin ^ o) coscfosinGoUjf’ U P=2V (a 2(x)/p2(x)-sin 2©o) sin4bsin©oUpr (2.16) U ^(a-^xV P ^x)— 2sin2©o)UI i: cos0 oe ltp where Up = (a 2 (x)/p2 (x)-2sin2 0o)2 +2sin20osin0oV (a 2(x)/p2(x)-sin2 0o) _____________________E I ___________________ 27ip1/2 ( x ) p 1 /2 (^ )p 3/2(x )p 1/2(^ )R P (^ x ) and tp, F p are the P wave travel time and products of reflection and transmission coefficients for the ray passing through the interfaces, respectively. For S waves, we have U J =cos20ocos<j*)US v0 — sin<j*)Ush < t > U | =cos20osin<jbUsv0 + cos<JoUsh < f > (2.17) U | = -2 V (a 2(x)/p2(x)-sin 26 o) sin0 oU sv0 where Usv cos0 oe U s c o s 2 2 0 q + 2 s i n 2 0 o s i n 0 o V ( a 2 ( x ) / p 2 ( x ) —sin20 o 2 ?tp 1/2 ( x ) p 1/2(^ )p 3/2(x )p 1/2(^ )R s(^ ,x ) t~ sh Ush = ------------------------1--------------------------- , 27t p 1 /2 ( x ) p 1/2(O P ( x ) ( i1/2( O R s( t x ) and ts, F sv and F sh are the S wave travel time, and products of reflection and transmission coefficients for SV and SH waves across the interfaces, respectively. In 22 all the above formulations, 0 o and < |> o are the initial ray take-off angle and azimuth, r, 0 and < |) are the unit vectors along the P, S V and SH wave directions, respectively. 2.3 Ray Tracing in 3-D Vertically Heterogeneous Media As was discussed previously, we use the ray tracing method to simplify the Green function calculation in an inhomogeneous media. In this modeling scheme, we assume that the velocity changes only in the vertical direction. This assumption is acceptable since we only deal with the wave propagation at relatively short distances and we also correct for the travel time fluctuation caused by the lateral variation in the seismic velocity. We further assume that the velocity increases linearly within each layer. Thus for a N-layer media, in which the depth to the upper boundary of V * 1 layer is hj, seismic velocity in this layer v(z)=Vi+ai(z-hj), we can find the ray travel time for the slowness p and the corresponding horizontal travel distance across this layer as •h i+i /»h i+i A t i = [ ----------- d z = [ d z Jh. v(z)Vl-v2(z)p2 Jh [v i+ a i (z-h i)]Vl-[vi+ai(z-hi)]2p2 1 ln(v rt ~ai(hi+1-hi) Vl-v?p2+l_______ j Vl V 1— [v i+ a i (h i+i-hi)] 2p2+l and 'hi+i /* h i+ i [vi+ai(z-hi)]pdz J phi+ i r ' , v ( z ) p d z = [ hi V i-v2(z)P2 ; h h| V l- v 2(z)p2 )h i Vl-[Vi+a;(z-hi)]2p2 = ^ r { Vl-v?p2 - V 1— [vi+ai(hi+1-hi)]2p2 } 23 1 2 n n+ 1 h N N Figure 2.3. Geometry of layered structure model with ray starts from the free surface down to the layer between n and (n+1)1 * 1 interfaces. 24 respectively. The total travel time and horizontal distance are then t = y i - 1nrVi+aih‘ , V l~v?P2 + 1 ■ ] (2 . 1 8) i=1 & 1 Vl V l-(vi+aih*)2p2+ l and r = X V l-v?p2 - Vl-(vi+aih*)2p2 ] i=l dlP (2.19) respectively, where h*=hi+i-hi for z > hi+i and h£=z-hn for hn+i > z > hn (see Figure 2.3). The geometrical spreading term can be obtained analytically as p 2_ rcos 6 dr = rco s 6 ocos6 dr sinGo G G o vosinGo dp _ rV 1— [vn+an(z-hn)]2p2V l-v 2p 2 y i ^______ i___________ l j (2 20) v iP3 i=1 ai Vl-(vi+aih*)2 p2 V 1 -v?p2 where G o and G are the ray take off-angles defined at the initial and destination point, respectively, and h j is defined above. Before applying equations (2.18) and (2.20), we have performed two point ray tracing by fixing the depth of the ray destination point and adjusting its horizontal distance according to equation (2.19) using Newton's method to find the corresponding ray slowness p iteratively. The reflection and transmission coefficients are also calculated for each ray at each interface, and included in the evaluation of the amplitude. 2.4 Other Numerical Treatments In order to carry out the isochron integration, we subdivide the fault into a number of grids. The integrand values inside a grid are interpolated from its values at the node of the grid. In this study, instead of utilizing some sophisticated high order polynomial 25 interpolation algorithm, we used a simple triangular linear interpolation scheme (see Figure 2.4a). The interpolated values of the slip intensity s(x, y), the projection g(x, y) of the traction term along the fault slip due to an asymptotic Green's function of an impulsive point source acting on the receiver, and the isochron time distribution t(x, y) on the fault are given by s = si + (s3-si)y + (s2-si)x g = gi + (g3-gi)y + (g2 -g i)x t = ti + (t3-ti)y + (t2-ti)x respectively, where the indices 1, 2 and 3 correspond to the node points 1, 2 and 3 as shown in Figure 2.4b, x and y are the normalized coordinates by the grid interval along x and y directions, respectively. By solving y from the expression for t, we find t-ti t2-ti y = ----- -— -—-x t3-ti t3-ti Plugging this into the expressions for s and g, we get s = bi + b2x g = ai + a2x where bi = si (S3 - si), b2 = ( s2 - si) - ( s3 - sO, ai = gi +7 ^ (g3 - gi), and t 3 - ti t 3 - ti t 3 - ti a2= ( g2- gi) - ( g3 - gi). The isochron integral provided in equation (2.9) over the t3-ti triangular cell in Figure 2.4b can now be converted to 26 Figure 2.4a. Triangular grid division used for the linear interpolation. 3 1 Figure 2.4b. An example of the triangular cell. 27 / ss c d L = J g s [i^ V ( '^/t )2+1 ]'V ( ^ ) 2+ l dx where x2 > x l5 and x 1 and x2 are the x-coordinates of the two end points of the isochron across the triangular cell. The total value of this isochron integral is obtained by summing up the contributions from each triangular cell it passes through according to equation (2 .2 1 ). When we perform the ray tracing and calculate the reflection and transmission coefficients, we also encounter imaginary values if the incidence angle exceeds the critical incidence angle of the interface. In this case, our synthetic ground motion will be of the form u=a(t)+b(t)i. Following Achenbach and Harris (1978), we corrected the ground motion using the form u(t)=a(t)— H [b(t)], where H [b(t)] is the Hilbert transform: We used some smoothing procedure to remove singularities over the integral. Examples of such singularities in isochron velocity calculation are caused by extrema or saddle points in the arrival time function. The method we used here is the same as that the boxcar function and At is the half length of the smoothing interval. This smoothing will result in a more physically meaningful synthesis since observations have shown a H[b(t)] (2.22) used by Chapman (1978) in which we convolve our data with — B(— ) , where B(t) is At At 28 limited frequency band width (Hanks, 1982) for the accelaration spectra of the strong ground motions. 2.5 Numerical Applications The computer program for calculating 3-D high frequency radiation from a complicated earthquake rupture model has been written using the theory discussed in the preceding section. With this program, we will first calculate synthetic seismograms for some simple models to discuss the effects of nucleation and barriers on earthquake high frequency radiation. Then we will compare our results with published ones for more complicated models. 2.5.1 Earthquake Nucleation and Stopping with a Hypothetic Isochron V elocity In Figure 2.5a, we have constructed a dislocation model with a hypothetic shear wave isochron velocity c=2.6 km/sec and a constant slip velocity. The effects of the Green function and the rupture healing process have been neglected in this calculation. The rupture starts from the center of the bottom edge of the fault and gradually propagates to cover the entire fault. Rupture stops immediately when it arrives at the rectangular edges of the fault. The synthetic displacement and velocity have been calculated and plotted in Figure 2.5b. As expected, displacement increases linearly after the arrival of the rupture starting phase and drops down with the time dependence of - V t - tc when rupture stops, where ^ is the sum of the rupture time and the ray travel time from the point at which the rupture front is tangent to the fault edge. Figure 2.5c is a comparison between the result obtained by equation (2 .1 0) and integrated result calculated using equation (2 .1 1), where we have replaced the second derivative slip 29 Depth Shear Wave Isochron ro CM o Figure 2.5a. Shear wave isochrons map for a constant shear wave isochron velocity c=2 .6 km/sec. Distance Along Strike (KM ) 30 Seismic Displacement o o ■ 2 1 0 1 2 3 4 TIME (SEC) Seismic Velocity 8 < N 8 o o o T 8 w I 2 0 1 1 2 3 4 TIME (SEC) Figure 2.5b. Synthetic displacement and velocity seismograms calculated from the isochron model shown in figure 2.5a with constant slip over the fault 31 Seismic Displacement o M 1 s a < 5 -2 0 -1 1 2 3 4 TIME (SEC) Figure 2.5c. Synthetic displacement displacement seismograms calculated using equation (2 .1 0 ) and integrated from the result calculated using equation (2 .1 1 ) for the isochron model shown in figure 2.5a with constant slip over the fault. 32 time function with its first derivatives in equations (2 .1 0 ) and (2 .1 1) to obtain the corresponding displacement and velocity, respectively. The purpose for this comparison is to test the ability to simulate the velocity stopping phase using our program. For a dislocation model, the slip distribution across barriers is a step function and its spatial derivatives will be a delta function. But in our numerical calculation, instead of the step function, we make a sharp change for the slip distribution across barriers. The comparison shows that our computed result is quite satisfactory. 2.5.2 Comparison with Spudich and Frazer (1984) To compare the results obtained by our program with that by Spudich and Frazer (1984), we digitized their rupture and slip intensity models. The digitized data were mapped on rectangular grids using an interpolation program for irregularly distributed data points from the IMSLS Library. The contours for the rupture time tr(^) and slip intensity s(£) are plotted in Figures 2.6a and 2.6b. The S wave isochrons for the receiver shown in Figure 2.7a were calculated using our program and are plotted in Figure 2.7b. The x-component ground velocity and acceleration seismograms were calculated and plotted on the left hand side in Figures 2.8a and 2.8b, respectively. The corresponding result from Spudich and Frazer (1984) are plotted in the same figure on right hand side for comparison. From these figures, we find a general agreement between the two results. A possible cause of some differences, especially for acceleration, is that we directly carried out the integral along the isochron whereas Spudich and Frazer smoothed their result by calculating the integration over the strip between two nearby isochrons. They have not used equation (2.11) in their calculation, which will result in the error pointed out earlier. 33 DeP t h _ Depth Rupture Time in to — - 0.5 0.5 0 1 1.5 2 2.5 Distance Along Strike (KM) Figure 2.6a. Contours of the rupture time function interpolated from the digitized Spudich and Frazer's model with contour interval 0.05 sec. Strike Slip Velocity Density 0 .2 0 .4 in in < o 2.5 1.5 2 0.5 1 0 0.5 Distance Along Strike (KM) Figure 2.6b. Slip density models interpolated from the digitized Spudich and Frazer's model. 34 hypocenter z velocity km/s 1 0 5 0 0 E .x f s « • X3 10 Figure 2.7a. Geometry of the fault and receiver, and velocity structure adopted from Spudich and Frazer. The fault extends from 5 to 6 km in the vertical direction and -0.5 to 2.5 km in the horizontal direction. Depth 5 .5 Shear Wave Isochrons to -0.5 0.5 0 1.5 2 2.5 Distance Along Strike (KM) Figure 2.7b. Shear wave isochrons at the receiver shown in figure 2.7a for the rupture in figure 2 .6 a C T \ GROUND VELOCITY I ------------- 1 -------------1 _______ L J_______ L 3 4 TIME (S) 1.0 1.5 2.0 2.5 3.0 3.5 TIME (S) 4.0 Figure 2.8a. Velocity seismogram of x-component for the receiver shown in figure 2.7a. Results of our calculations (left) are compared with Spudich and Frazer (right). x GROUND ACCEL, 1 2 3 4 TIME (S) "J Figure 2.8b. Acceleration seismogram of x-component for the receiver shown in figure 2.7a. Results of our calculations (left) are compared with Spudich and Frazer (right). 2.5.3 Comparison with Discrete Wave Number Accelerograms Following Bernard and Madariaga (1984), we shall also compare the result from our asymptotic method with that of the result obtained by the more accurate discrete wave number method of Bouchon (1982). Bouchon used a dislocation model to explain the Gilroy 6 recording of the Coyote Lake earthquake of 1979 in California. His model is a vertical strike slip fault with rectangular edge (Figure 2.9). Rupture starts at 9.5 km depth and propagates at a constant velocity of 2.6 km/sec. Slip equals to 21 cm everywhere on the fault. The medium has a top layer 1.75 km thick, a shear wave velocity of 2.4 km/sec and a density of 2.6 gm/cm3 overlying a half space with a shear wave velocity of 3.5 km/sec and a density of 2.8 gm/cm3. Using Bouchon's model, we have computed the shear wave isochrons as shown in Figure 2.10. The synthetic displacement calculated using our program were compared with Bouchon's result as well as Bernard and Madariaga's result in Figure 2.11. For all results, the starting phase and the stopping phases can be clearly identified. The bump around the S wave first arrival in Bouchon's result and in our bandpassed asymptotic result are caused by a noncausal high pass filtering effect. The principle difference between these two results is the rapid increase of the displacement just before phase 2 which does not appear in Bouchon's result. As pointed out by Bernard and Madariaga (1984), this difference is probably due to the discretization of the source in Bouchon's model, where the top row of elementary sources in his discrete model is deeper than the nominal depth to the top of the fault, resulting in an absence from his model of the rapid change in reflection coefficient that causes the increase of the displacement in our result. In general, the agreement between the two results is excellent. Note that in this figure, there is a small difference immediately before phases 38 2 and 3 between our result and that of Bernard and Madariaga. By comparing velocity seismograms (Figure 2.12), we find a gradual increase by velocity in Bernard and Madariaga's calculation, whereas our calculation shows a gradual decrease in velocity. The S wave Green's function takes on negative values for all isochrons within the isochron which is tangent to the upper edge. Furthermore, these isochrons gradually approach the station which effectively increases the absolute integral value. The resultant velocity should keep decreasing except for the phase 4 arrival. So we believe that there were some errors in their calculation. In fact, our result fits that of the discrete wave number accelerograms better than Bernard and Madariaga's. 39 Epicenter Gilroy 6 1.75 km t 7.5 km 0.5 km 2 km 1 2 km Upper layer Rupture Front Sli p 4 ^ ■ » 0 km Figure 2.9. Fault and receiver geometry adopted from Bouchon's (1982) rectangular dislocation model for the Gilroy 6 recording o f Coyote lake earthquake of 1979 in California. The shear wave velocity is 2.6 km/sec. in the upper layer and 2.8 km/sec. in the lower half space. S hear Wave Isochrons Distance Along Strike (KM) Figure 2.10. Shear wave isochron for the model and receiver shown in figure 2.9. 40 DISPLACEM ENT - i L 4 6 8 TIME (SEC) ASYMPTOTIC METHOD T I M E BOUCHON' S METHOD 9 .3 cm Figure 2.11. The synthetic displacement seismogram from our calculation (left), and that from Bouchon (lower right) and Bernard and Madariaga (upper right). VELOCITY 5 0 cm/s t I 4 6 T I M E (SE C ) Figure 2.12. The synthetic velocity seismogram from our calculation (left) and that from Bernard and Madariaga (right). to Chanter IT T RECURSIVE STOCHASTIC INVERSE ALGORITHM 3.1 Introduction Consider a set of data d which is related to a set of model parameters m by a theoretical relation d=G(m). For example, the relationships between seismogram and earthquake faulting process described by equations (2.9), (2.10), and (2.11). The corresponding inverse solution can be obtained by applying a certain optimization criterion to the linearized system, provided that an adequate initial model estimate is available. For example, in seismological inverse studies (e.g. Aki and Richards, 1980) the generalized (Lanczos, 1961) and the stochastic (Franklin, 1970) inverse solutions have been used. These solutions, however, require a large amount of computational time and memory space when a large number of model parameters are involved. In order to save computer time and memory space, iterative solutions such as ART (Gordon et al., 1970) which do not require matrix inversion have been introduced to seismology (Mason, 1981; McMechan, 1983; Humphreys et al., 1984). However, these iterative solutions do not give error and resolution estimates for the model parameters and sometimes may suffer from convergence problems. In this chapter, we will show that the stochastic inverse process can be reduced to a recursive scheme which does not require matrix inversion. The basic idea is as follows: Suppose we obtained a stochastic inverse solution for a linearized inverse problem using a set of k observational datum. Our problem is to find how the above solution must be revised if we have an additional data. We will show that the solution for the (k+1) data can be obtained in a simple recursive manner using the solution for 43 the original k data plus the new one. This recursive formula gives the same stochastic inverse obtained by the exact matrix inversion as well as the error and resolution matrices of the solution. Since the recursive formula can take advantage of the sparseness of matrix does as the iterative method, it leads to a fast and efficient inverse scheme for certain kinds of large sparse matrix equations that commonly arise in seismic tomographic problems, including the inversion for fault slip distribution using the isochron method discussed in Chapter II. 3.2 Recursive Stochastic Inverse Algorithm Suppose we have a linearized physical system described by where m is a model vector, n is a noise vector, d is a data vector, and G is a function governed by the physical process of the model. Given a noise covariance matrix Rn, an a priori model covariance matrix Rm , and assuming zero means for both model and noise vector, the stochastic inverse solution for (3.1) is obtained (Franklin, 1970) as Gm + n = d. (3.1) m = R rn G T ( G R mG T + R p ) ^ d (3.2) with error covariance matrix and resolution matrix given by R Am = R m - R m G T ( G R m G T + R „ ) 1 G R m (3.3) and M = R m G T ( G R mG T + R „ ) ^ G , (3.4) respectively. 44 We call the procedure of obtaining an inverse solution by the above equations a batch stochastic inverse process. This process involves a direct matrix inversion. In this section we will derive a recursive stochastic inverse procedure which breaks the above inversion into a sequence of scalar inverse process. In this derivation, we assume that the model vector and noise vector obey the joint Gaussian distribution with zero means, and the noise vector is uncorrelated between its components and with the model vector. The result obtained under the joint Gaussian assumption applies exactly the same as to the linear stochastic inverse of the non-Gaussian case (Mendel, 1983). Let us first introduce the following notation; dk, nk : kth components of the data vector d and noise vector n, respectively. Dk : data vector given by ( di, • • • , dk )T. m : true model perturbation independent of data size k. irik : stochastic inverse solution of m given Dk. Arrik : solution error defined a sm - irik. Gk : kth row of the partial derivative matrix G. Adk+i • (k+ljth component of data residual defined as dk+i — Gk+imk- R^m, k • error covariance matrix given by E[Amk(Anik)T] in, k • noise variance for nk. Mk : resolution matrix of the stochastic inverse based on Dk. We now rewrite equation (3.1) for the component of the data vector as G km + nk= d k , k = 1,2, — , N. (3.5) In order to derive the recursive formula it is necessary to use the orthogonality of stochastic inverse solutions. According to the orthogonality principle (Mendel, 1983) any function of the data set Dk has a stochastic inverse solution error which is statistically orthogonal to f(Dk). We express the orthogonality principle as E tA m ^ D i) ] = 0 . (3.6) We assume that x, y and z are joint Gaussian random vectors with zero mean. The stochastic inverse solution of x for known y and z satisfies the following equality (Mendel, 1983) E(x I y, z) = E(x I y) + E(x I Az) (3.7) where Az = z - E(z I y). From the above equation we find m k+1=E (m lD k, dk+1)=E(m I Dk)+E[m l dk+1-E(dk+ 1 lDk)]. With equation (3.5) and the assumption that the noise is uncorrelated between itself and with model parameters, we find E (dk+i lDk)=Gk+iE(ml Dk)+E(nk+i lDk)=Gk+iink. This gives us fnk+i=m k+E(ml dk+i-G k+im k)=mk+E(inl Adk+i)=fnk+Kk+iAdk+i (3.8) where Kk+ 1 = Rm A dk„ . 46 Since m k is a function of Dk, the orthogonality principle (equation 3.6) gives E (A m km k) = 0 . Thus we obtain R m A 4 +i =E [m (G k+i Am k+nk+i )]=E(m Am J)G k+1+E(m nk+1) =E [(A m k+m k) A m k]G k+i= R A m kG k+1 , and likewise T fA dk.! = E [(G k+iA m k+nk+i)(G k+iA m k+nk+1) ] = G k+iE (A m kAm £ )G k+1+E(nk+1)= G k+iR A m kG k+1+rn, k+1. From the above results the expression for Kk+i in (3.8) is Kk+l = ^Am, k Gk+i/(G k+i ^Am, k^k+l"*"rn, k+l)* (3.9) We can also derive the recursive relation between Rahi, k and Rah i, k+i through the following steps A m k+i = m -m k+] = m -m k- K k+] Adk+i =A m k- K k+i Adk+!, R A m , k+i=E(Am k+1AmJ+1) =E (Am k Am k)+ Kk+k E (Adk+ iAdk+j) K k+r K k+1E (Adk+iAm k)-E (Am k Adk + {) K k+j — RA m , k +R k+ 1 R A d, k+lR k+l _Rk+l Rm A dk + i ~ Rm A d k + i R k + 1 = R A m , k-K k+iR L 4 ., = R A m , k-K k+iG k+1RA m > k . (3.10) 47 Rewriting equations (3.8), (3.9), and (3.10), we obtain the formulation for the recursive process of our stochastic inverse for k=0, 1, ..., N -l as m k+i=m k+K k+1Adk+i ! Kk+i = Ra h i, k^k+i/(Gk+iRA m > kG£+1+rn ? k+i) (3.11) ^ A m , k + l = ^ A m , k ' ^ k + l G j c + i R ^ k where Adk+i = dk+i — Gk+imk, mo=0 and Ra™, o * s defined a priori. The stochastic inverse solution given by the batch inverse process equals E(m I Dn). According to the above recursive process, it can be partitioned into a relation of the stochastic inverse solution E(m I D^-i) plus a revision based on d jM - Likewise E(m I Dn_i) equals E(m I Dn_2) plus a revision based on dN-i, and so on. Thus the original batch inverse process is converted into an equivalent scheme of a recursive scalar inverse process with each step giving a solution based on the data points up to it. This recursive inverse process was also called sequential estimation in early literature (Deutsch, 1965; Rodgers, 1976; Tarantola, 1987). But the method has been often considered as an alternative to the batch process without much computational advantage (Rodgers, 1976; Mendel, 1983). In the following discussion we will show that for seismological inverse problems involving large, sparse matrices, it provides a fast algorithm which could save a large amount of computer time. To complete our recursive process, we will also derive a recursive formula for the resolution matrix given in equation (3.2). A simple way to obtain the recursive formula is by comparing equation (3.2) with equation (3.4). If we replace each data component dk in equation (3.2) by Gk and the inverse solution m by the resolution matrix M, we arrive at the same expression in equation (3.4). Following this exact procedure and 48 replacing all the necessary terms in the recursive inverse solution expression in equation (3.11), we obtain where I is an identity matrix. Combining equations (3.11) and (3.12), we obtain our modified recursive stochastic inverse process for k=0,l, ..., N -l as where Adk+i = dk+i — Gk+imk, mo=0 and RA n ii o defined a priori. In the above recursive process, we see clearly that each accumulated data point contributes to the model revision only through its residual. Any additional observation, no matter how noisy it is, will always reduce the value of the error covariance matrix and improve the model estimate and resolution. 3.3 Model Covariance Matrix RA m o In this section, we shall discuss the issue of the model covariance matrix Rm=RAm , o» which brings a subjective element into our stochastic inverse process. A proper choice of Rm is very important because it constraints the deviation of our inverse solution from our initial guess. The simplest choice for Rm is to take it as a diagonal matrix with the same value for all the diagonal terms. This could lead to some disadvantages as we will explain in the following. For simplicity, we consider the equation (3.1) without the noise vector. The covariance matrices are then related by Mk+i =M k+Kk+i (Gk+i -Gk+i M k)=M k+Kk+i Gk+i (I-M k) (3.12) mk+1 =mk+Kk+1 Adk+1 K k + l = R aiii, k G k + l/(G k+lRAm, k G k + l+ rn, k+l) ^Am, k+ l- ^Am, k ”^ -k + lG k + lR Am k \ Mk +i=Mk+Kk+iGk+i(I-Mk ) (3.13) 49 G E (m m T)G T= E (d d T) or G R mG T= R d For R m = d m I, we can find ®m £G?k= r d i k where r^ is the data variance representing the signal energy and is the model variance. If the theoretical relation G is very sensitive to the model parameter mj, then the variation of this model parameter will control the variation of the data prediction. This becomes a problem when we deal with a nonlinear system where G gives a very unstable local mapping between model and observation data spaces. Furthermore, in some cases, we may already have certain control over these parameters in our initial model. For example, in our fault rupture time inversion, we are relatively certain about the location of hypocenter. In addition, we know that in the vicinity of the hypocenter, the fault will rupture in a way more or less close to the initial uniform model. For places further away from the hypocenter, the rupture time tends to more likely deviate away form initial model. So we have variable confidence over the choice of the initial model parameter which implies that the above choice of the same may not be adequate. One way to remedy this is to choose the initial model variance in such a way that it is inversely proportional to its corresponding information matrix values explained below. By assuming that the noise vector in (3.1) is Gaussian, we can find the probability density function as p ( d lm ) = (27to£) 'n/2exp[ - (d - Gm )T (d - Gm)] 2o* 50 where a 2 is the noise variance. The Fisher information matrix for this probability density distribution is obtained as I = ( 9 ^ 1 n £ ) = j _ G T G d m 2 Thus, we get o 2 C _ cog _ c* mt 1/OnX G?k X G fk S G ? k i i i We take the cross terms of the initial model covariance matrix as I-ij = C 7 m s O n,j exp{ - r [(^i- £J )2/a? + (T 1 i - T|j )2 /A 2] } where £ , and rj are the coordinates of the points on the fault plane. Such choices for the cross terms provide the same effect as the smoothing operator introduced by Hartzell and Heaton (1983). We can demonstrate this equivalence by a simple two model parameter example. For this case, the inverse solution of Hartzell and Heaton (1983) was obtained by minimizing the objective function minimize { - ~ [ d - G | )]T [ d - G ( J]+^i(m i-m 2)2+X2 (m?+m5)} The solution is given by T V c TG+og( Xl+^2 -x, V -i G T d (3.14) i m2 / I -*1 A ,i +X,2 /_ Our method gives ( W g t G+c s 2 R^>] '1 GTd I m2 I 51 Assuming Rm = | ^ ^ J , we then find GTG+og \m 2/ (3.15) where D = ri ir22 -r j:. The equivalence between our result and that of Hartzell and Heaton (1983) is shown by the same basic structure between equation (3.14) and (3.15), giving a physical meaning to Xx and X2 in terms of the correlation distance A. 3.4 Noise Covariance Matrix Rn Unlike the model covariance matrix, the data noise covariance matrix is a littile more straightforward to obtain. By assuming that the noise is uncorrelated, we get We further assume that the noise vector is uncorrelated with the model vector. According to equation (3.1), we can find By looking at the above equation from the energy viewpoint, the norm of GRm GT can then represent the signal energy and the norm of Rn represents the noise energy. Suppose the signal to noise ratio is r, the noise variance is then constrained by where Xmax is the maximum eigenvalue of GRm GT. In our current study, we take the equality sign in the above formula. Such a choice could be important in obtaining a stable inverse solution. Of course, we should be able to measure the noise variance Rn = og I E (ddT) = Rn+G R mGT a 2 < A m a x n r 52 directly from the data if there is no modeling error, which is usually not the case in reality. 3.5 Applications to Seismic Inverse Problems In this section, by showing some typical examples of seismic inverse problems, we will demonstrate that the recursive stochastic inverse process given by equation (3.13) provides a fast and efficient algorithm in dealing with inversion of large, sparse matrix systems. We divide the discussion of the applications into three topics; the seismic tomographic problem, source imaging based on the isochron method and Inversion of site effect using coda power spectra. 3.5.1 Seismic Tomographic Problem To illustrate the seismic tomographic problem, let us take a 2-D case of a square area (see Figure 3.1). The area is divided into N=m2 number of sub-grids with each having a constant wave velocity. All the sources and stations are deployed along the border of the square. We further postulate that rays travel along straight lines. Based on these assumptions, we write the travel time for ray i (see Figure 3.1) as I - k = I lijPj = ti, i = l L (3.16) j=l J j=l where lij is the length of the ray contained in grid j, is the travel time of the ray from the source to the station and Pj=l/vj with vj as the wave velocity of the grid. Putting G={ 1 ij}, P=(Pi, ..., Pn)T and t=(t1?..., tL)T, we can rewrite equation (3.16) in vector form as G P = t 53 i-th ray m m Figure 3.1. Exemplary segmentation of i-th ray passing through the square-area which is divided into N=m2 number of square sub-grids with each having a constant velocity. 54 Through examination of each ray path, we find the actual number of non-zero elements of 1 ij is about the order of m=VN, indicating that G is indeed a very sparse matrix. By applying the batch stochastic inverse method, we obtain P = RPGT (GRPGT+Rn)"1t (3.17) where the model covariance matrix Rp is defined a priori and R n is the noise covariance matrix. In this batch inverse process, we first compute (GRpGT+Rn), which becomes a non-sparse matrix. Using the Gaussian elimination method, we can find the inverse solution given in equation (3.17). By assuming each multiplication and addition as one computation counts, we can estimate that the total number of computation counts for the above process to be of the order of N3 (Press et al., 1986). Using our recursive inverse process provided by equation (3.13), we can compute the same inverse solution. In this recursive inverse process, the only time consuming calculation is to compute R ap, k^k+i and Gk+iM^. By tracing each non-zero element in Gk+i and neglecting the zero ones in performing the calculations, we find our number of computation counts to be of the order of NVN~=NL5. The total computation counts after all the recursive steps will be about N l5L, where L is the total number of data points. If we assume the total number of data points is of the same order as N, then our total number of computation counts will be around N2- 5. As a result, our recursive stochastic inverse scheme is about V N " times faster compared to the direct matrix inverse method. Although we may use iterative methods such as ART (Gordon et al., 1970) or tomographic back projection (Herman, 1980) to take advantage of the sparseness of matrix for a fast inversion of a large amount of data, these methods may suffer from the convergence problem and they usually do not provide the error covariance and resolution matrices of the solution which are important in judging the 55 quality of the result. Through a similar analysis for the 3-D seismic tomographic problem, we find that using the recursive algorithm can lead to N2/3 times reduction in computer time compared to the direct matrix inversion. 3.5.2 Source Imaging Using the Isochron Method Another example of applying the recursive stochastic inverse algorithm is for the source imaging problem using the isochron method (Spudich and Frazer, 1984). This method, based on the high frequency approximation, uses ray theory to calculate the near field seismogram. By assuming the same slip time function over the entire fault area, it reduces an area integration over the fault to an isochron line integration. In our study, we took a similar fault rupture model as Beroza and Spudich (1989). The corresponding station distribution, fault geometry, shear wave velocity model, slip velocity intensity distribution and rupture time function are shown in Figures 3.2a, 3.2b, and 3.2c, and Figures 3.3a and 3.3b, respectively. Using the isochron method, we calculated the vertical component synthetic seismograms for our stations as shown in Figure 3.2a. Then an inverse solution was obtained using these synthetic data to recover the rupture time function, which was parameterized into 800 equally-spaced points on the fault. A perturbation scheme was used to linearize this highly non-linear model. The seismograms were deconvolved with the fault slip time function which resulted in the data point depending only on an isochron integration on the fault. This leads to a very similar image problem as does 2-D seismic tomography discussed earlier. Thus we were dealing with an inverse system involving large sparse matrices. A direct matrix inverse method was first used to calculate the stochastic inverse by calling the matrix inverse routines provided in the IMSLS Library Package. However, due to the ill-posed nature of the problem and the large size (800x800) of the inverse 56 Station Distribution around the Fault (a) Velocity km/s Hypocenter Fault Strike Fault Depth (b) 0 S m J S 5 a a 10 (c) Figure 3.2. (a) Map view of the fault trace (AA') and station distributions (triangles), (b) geometry of the fault Cross section along the fault strike, and (c) depth dependent P-wave and S-wave velocity structure of the model. U i D epth (km ) D epth (km) Slip Velocity Density m in in C O -1 .5 - 1 -0 .5 0.5 Distance Along Strike (km) Figure 3.3a. Slip velocity density contours with value of the most inner line equals 1.0 m/sec and the interval value o f any adjacent lines equals 0.1 m/sec. Rupture Time m in in O C O - 1 -0 .5 0 0.5 t 1.5 Distance Along Strike (km) Figure 3.3b. Fault rupture time contours with the starting point near the up-left comer and 0.025 sec contour interval. 58 D epth (km ) D epth (km) R upture Time m 1.5 1 0.5 0 -0 .5 1 Distance Along Strike (km) Figure 3.4a. Inverse solution o f the fault rupture time distribution obtained using Chebyshev semi-iterative method. Rupture Time < o 0.5 0.5 - 1 Distance Along Strike (km) Figure 3.4b. Inverse solution o f the fault rupture time distribution obtained using our recursive stochastic inverse process. 59 matrix, we could not obtain a solution without an error message stating that our matrix was too ill-defined. Following the same Chebyshev semi-iterative procedure used by Beroza and Spudich (1989), we then computed the result shown in Figure 3.4a. This result is very similar to that obtained by Beroza and Spudich (1989). Figure 3.4b shows the result for the same inverse problem obtained using our recursive inverse scheme. The computer times for both methods are comparable. Of course, our inverse solution obtained through one iterative step was not exactly the same as the original model due to the non-linear nature of the problem. However, we have a result closer to the true solution plus the error covariance and resolution matrices, which could be very difficult to obtain using the Chebyshev-semi iterative scheme. Comparing the computer time of our recursive stochastic inverse with that of the batch stochastic inverse, our approach is about 20 times faster than the direct matrix inversion. 3.5.3 Inversion of Site Effect Using Coda Power Spectra We begin this section with a short review of the inversion of site amplification factors using coda power spectra. Phillips and Aki (1986) analyzed a large amount of seismic data from the USGS regional network in California, assuming that the coda wave power spectrum can be expressed as a product of site, source and average path effect in the following form ~2 InP ij(CDi, tk )= d ijk i= r i(coi)+Sj(coi)+c(to1, tk) (3.18) where rj(coi) is the site term, sj(coi) is the source term, and c(C0i, tk) is the coda propagation term which has been shown to be independent of source and receiver locations. Indices i, j, k and 1 represent the station, source, lapse time and frequency, respectively. By taking an average of dijki over all the available stations with fixed 60 indices j, k, and 1 and then subtracting it from the original dijki values, we arrive at the following equation dijki-dijk^n-fi (3.19) We can further write T[= £ 5imr m and f [ = 1/NjkiZ Imim, where Njki is the total number m m of usable stations for fixed indices j, k, and 1 and Im=l if station m is used, or Im=0 otherwise. Substituting these expression into (3.19), we obtain dijkrdijki=Z (Sim-Im/Njki)rm= S G(K)rm (3.20) m m where G(K) = 5im- Im/Njki and K is a function of indices j, k and 1 . For a given frequency (fixed index 1), equation (3.20) gives us a system of linear equations for different stations, sources and lapse time points which could be used to determine the relative site amplification factors. In practice the matrix provided by the linear equations could be very sparse. A singular value decomposition and generalized inverse technique used by Phillips and Aki (1986) is not only expensive when large number of model parameters and observation data points are involved but also inefficient. A fast and efficient algorithm could be obtained based on our recursive stochastic inverse. As addressed before, the most time consuming part of our recursive process is to calculate the product of Rr(K)GT (K +l), where R r(K) is the error covariance matrix in step K. According to equation (3.20), this can be obtained as R r(K)GT (K +1 )=ith row of Rr(K)- l/N jk,[Rr(K)( ■ • • Im ■ • • )T] where the computation of Rr(K)( • • • Im • • • )T equals the summation of all rows with Im=l. Since there are no multiplications involved in this computation and we ignore all 61 the calculations for Im=0, tremendous computer time is saved in this inverse process. In addition, with new data points we can always revise the solution according to the recursive inverse process and avoid the redundant computation once again over the whole data set. 3.6 Discussion and Conclusion The stochastic inverse method is one of the most commonly used algorithms in seismological inverse studies. But the conventional batch stochastic inverse process can be extremely time consuming for inverting large sparse matrix equations and sometimes we may not even be able to get any useful results. The recursive stochastic inverse process introduced in this paper not only leads to the same answer at its final step as that of the batch process, but also provides a fast and efficient solution. In addition, this recursive scheme provides some insight to the inversion process. From the recursive process given in equations (3.13) we see that each accumulated data point contributes to the model revision only through its residual. Any additional observation, no matter how noisy it is, will always reduce the value of the error covariance matrix and improve the model resolution and estimations. Thus these intermediate results indeed give us a good perspective of how each observation improves our inverse solution, helping us to design a better observational strategy. The procedure of our recursive process is a sequence of scalar inverses which requires only matrix and vector additions and multiplications. For the inversion of a r p large sparse matrix, we can compute the term R A m > k^k+i by tracing only the non-zero elements in G^+i ? which avoids a large amount of unnecessary calculations. Thus the computation time is greatly reduced. It is particularly effective in inverse problems with a large number of model parameters and observations. Seismic tomography, 62 earthquake fault rupture imaging and site inversion using coda power spectra are all typical examples of such problems. In application to the seismic tomographic problem, we found that the recursive algorithm is about Na times faster in computation compared with direct matrix inversion, where N is the size of the matrix equations and oc=0.5 for 2-D and 2/3 for 3-D problems. As a test, we also applied our method to the inversion of a fault rupture process using a synthetic data set analogous to that of Beroza and Spudich (1988). The result obtained by our recursive inversion process is better than the one obtained from the Chebyshev semi-iterative technique, and a nearly 20-fold reduction in computer time is achieved compared to the direct matrix inversion calculation. Our recursive method also leads to a simple rapid inversion procedure for inversion of site amplification factors using coda power spectra with a tremendous saving of computer time. 63 Chanter IV APPLICATION TO THE WHITTIER-NARROWS EARTHQUAKE 4.1 Description of the Problem With our theoretical preparation completed in Chapters II and III for the forward and inverse methods, we are now ready to apply them to actual data. We will first apply our methods to the Whittier-Narrows earthquake of 1987. The Whittier-Narrows earthquake occurred in the east Los Angeles metropolitan area at 1442 (UT) on October 1, 1987. The hypocenter of this earthquake was 34° 2.96' N, 118° 4.86’ W, at a depth of 14.6 km (Hauksson and Jones, 1989). Although it is a moderate-sized earthquake, it caused substantial damage to buildings and has become one of the best studied events in modem strong motion seismology. The mainshock occurred within one of the most densely instrumented areas and was recorded by a large number of strong motion instruments, providing us with a unique opportunity to study the earthquake source, site and path effects. The first motion source mechanism (Hauksson and Jones, 1989), geodetic modeling (Lin and Stein, 1989) and study of teleseismic body waves (Bent and Helmberger, 1989) all suggested that the mainshock was associated with a blind thrust on a fault plane striking approximately EW. The purpose of this study is to image the mainshock faulting process using the large number of available strong motion records and learn about the dynamics of the rupture propagation. 64 4.2 Data Collection and Modeling The mainshock of the Whittier-Narrows earthquake has been well recorded by the following three strong motion networks, the USC Los Angeles basin strong motion network, the California Division of Mines and Geology network and the U.S. Geological Survey network. There are over 100 three-component strong motion records available for analysis. Due to the complex stmcture of the Los Angeles basin, we try to minimize the wave propagation path effect by using only stations within 15 km of the epicenter. For such a short distance range seismograms may be dominated by the direct S-waves and suffer less from complication over wave paths through the heterogeneous velocity structure under the L.A. basin area. We eliminated the records at the station at 7215 Bright Ave, Whittier, which was strongly modified by the building structure response (Kawase and Aki, 1990). A total of 16 stations were used for our analysis. Table 4.1 lists the names and locations of these 16 stations. A map view of the station distribution is plotted in Figure 4.1 with the mainshock location marked by the star. We see that our station distribution provides good azimuthal coverage around the source, which should give good resolution for the source inversion (Hartzell and Iida, 1990). The instrument-response corrected ground acceleration records have been integrated twice to obtain the ground displacements. These displacement records were bandpass filtered from 0.6 to 4.5 Hz and then interpolated at a uniform time interval of 0.08 second. The filtering of the low frequencies is to limit the data within the frequency range of the validity of our ray method and the filtering of the high frequencies is due to our inability to model the detailed seismic velocity structure over the ray path as well as to account for the complicated site response, to which high 65 Table 4.1 Strong Motion Station Name and Location Station Latitude longitude 1 USC19 San Gabriel 600 E. Grand Ave. 34.091 118.093 2 USC66 El Monte 11338 Fairview Ave. 34.093 118.018 3 USC93 Arcadia 180 Campus Dr. 34.130 118.036 4 USC69 Baldwin Park 3699 N. Holly Ave. 34.100 117.974 5 USC94 South Gate 7420 Jaboneria 33.965 118.158 6 USC71 West Covina 1307 S. Orange 34.064 117.952 7 USC73 Hacienda Heights 16750 Colma 33.990 117.942 8 USGS Garvey Resevoir 34.050 118.110 9 USGS Whittier Narrows Dam-Upstream 34.030 118.050 10 USGS Alhambra 900 S. Fremont 34.090 118.150 11 USGS Los Angeles-Bulk Mail Center 33.990 118.160 12 USGS Vernon 4814 Loma Vista Ave. 34.000 118.200 13 CDMG Alhambra-Fremont School 34.070 118.150 14 CDMG San Marino-S Western Acad. 34.115 118.130 15 CDMG Los Angeles-Obregon Park 34.037 118.178 16 USGS Norwalk 12400 Imperial Hiway 33.920 118.070 66 J 1 I L . _ l____1 ____1 _ A 14 CDMi • 4 USC ♦ 10 USGS • 1 USC A 13 CDMG • 2 USC 6 USC ♦ 8 USGS ☆ 5 .9 MAINSHOCK A 15 CDMG ♦ 9 USGS 34° — ♦ 12 USGS ♦ 11 USGS • 5 USC • 7 USC 5 K M j—i i i i 1 1 1 r ♦ 16 USGS _ i --------1 --------1 --------1 --------1 --------1 --------1 --------1 --------1 --------1 --------1 --------,--------1 --------1 --------1 --------1 --------1 --------r 10’ 118' Figure 4.1. Map view of the strong motion station distribution used for the Whittier- Narrows earthquake source inversion study. frequency waves are very sensitive. The vertical components are not used in our study. One of the reasons for this is due to the triggering procedure of the instrument, some P- wave arrivals were lost from the records. Another reason is our incomplete knowledge of the seismic velocity structures for both P- and S-waves. The mainshock faulting is modeled as a 12 km by 12 km square planar region with the hypocenter located at the center of the fault plane. Following the work of Hauksson and Jones (1989), we take the hypocenter at 34° 2.96' N, 118° 4.86' W, at a depth of 14.6 km. The choice of the fault strike and dip angle, however, is not straight forward. Hauksson and Jones (1989) obtained a strike of 270° and dip of 25° from the local short-period first motion data. Bent and Helmberger (1989), on the other hand, modeled teleseismic body wave data and found the best fit to their model has a 280° strike and 40° dip. Lin and Stein (1989) modeled geodetic data, obtaining the same strike direction as Hauksson and Jones (1989) but preferred a dip angle of 30°. In our study, we shall vary the fault strike from 270° to 280° and dip from 25° to 40°. After many numerical trials, we find two preferred cases, namely, the fault strike direction and dip angle at 270° and 35°, and 280° and 30°, respectively. Based on these two different fault strike and dip combinations, we carried out an inversion study of the earthquake faulting process. The velocity structure (Figure 4.2) is adopted based on the works of Hartzell and Iida (1990), Wald et al. (1988), Hauksson (1987) and Apsel et al. (1981). Following Hartzell and Iida (1990), the high velocity gradient in the top 6 km is chosen to simulate the low velocity sediments in the Los Angeles Basin where most strong motion stations are located. The average rupture velocity is chosen as 2.5 km/sec and the corresponding rupture time contour is plotted in Figure 4.3 with 0.5 sec time interval 68 D epth (km) Velocity Structure for Whittier Narrows 0 2 4 6 8 o o o C V 2 o C O o T t< L 1 0 2 4 6 8 Velocity (km/sec) Figure 4.2. Seism ic velocity model for the Los Angeles basin area. 69 Fault D ip (km) Rupture Time Contour o in o 0 5 - 5 Fault Strike (km ) Figure 4.3. Contours o f rupture time model for the Whittier-Narrows earthquake with constant rupture velocity of 2.5 km/sec and contour interval of 0.5 sec. 70 between neighboring contour lines. The initial model has a uniformly distributed fault slip of 30 cm. For the sake of simplicity, we only take the dip slip component in our modeling. The slip time function is assumed to be [H(t-tr)-H(t-th)]/V t-tr, where the rise time defined by th-tr is chosen to be 0.2 sec. We also tried rise times of 0.25 and 0.3 sec and found 0.2 sec provided better results, though in general, the result was not very sensitive to the choice of rise times within this range. The synthetic ground motion displacement has been bandpass filtered in the same way as the observed records. 4.3 Inversion Result and Analysis As mentioned earlier, we tried several choices for the fault strike directions at 270° and 280°, and dip angles at 25°, 30°, 35° and 40°. Based on the fitting to the observational data, we found two preferred combinations of fault strike directions and dip angles of 270° and 35°, and 280° and 30°. The results present in this section will be based on these two fault geometries. We first discuss the result with fault strike of 270° and dip of 35°. For this fault geometry, we determined the fault slip intensity distribution from the strong motion displacement records of the 16 stations selected with constant rupture velocity of 2.5 km/sec. The resultant slip distribution on the fault is plotted in Figures 4.4a and 4.4b. In Figure 4.4a, we show the fault slip in cm with different shading intensity according to the scale indicated below the figure. Figure 4.4b is a contour plot of the same slip distribution as in Figure 4.4a. In these figures, we recognize two major slip regions below the earthquake hypocenter with peak slip of about 120 cm. The synthetic seismograms calculated for this fault slip model are compared with the observations in Figure 4.5, where the observed traces are plotted above the synthetic seismograms for each station. The vertical axis corresponds to the station number listed in Table 4.1. The peak displacement (in cm) and epicentral 71 Slip Amplitude Distribution - 5 0 5 Fault Strike (km) 0.000 21.30 42.80 83.90 83.20 108.3 Figure 4.4a. The fault slip (cm) distribution of the first inversion shaded in different intensity according to scales indicated below the figure. -J to Fault D ip (km) Slip Amplitude Distribution o O in o - 5 0 5 Fault Strike (km) Figure 4.4b. Contour plot of the same slip distribution as in Figure 4.4a with contour interval of 15.21 cm. U SC 2.71 4.91 use 1.54 7.61 use 1.47 9.82 1.70 use 0.82 11.41 1.03 use 2.53 11.62 1.69 U S C 0 .6 3 4.91 1 0.27 U S C 0.56 7.61 2 0.31 use 0 46 9 82 3 0.2 4 use 0 .2 8 11.41 4 0.33 use 2.26 11.62 5 usc 0.81 12.03 6 0.59 usc 0.61 14.57 7 8 9 0.51 usg 0.93 7.92 0 8 -2 0 10 2 4 6 usc 1.04 12.03 0.68 usc 0.52 14.57 usg 2.09 3.11 usg 0.97 3.11 usg 2.38 7.92 - 2 0 2 6 Figure 4.5. Plot of seismograms for the first inversion with the observed traces on top of the synthetic ones for each station numbered along the vertical axis. Each seismogram is normalized to its peak values (cm), indicated to the right the trace followed by its epicentral distance (km). Seismograms to the left are components parallel to the fault strike and those to the right are perpendicular to the fault strike. 73 11 . 12 . 13 . 14 . 15 . v A jV V 7^ usg 3.56 10.06 1.26 usg 1.95 12.42 3.76 cdm 1.26 6.94 0.52 cdm 0.33 8.50 0.22 cdm 1.65 9.05 1.57 usg 1.60 10.06 0.90 usg 1.69 12.42 cdm 1.95 6.94 0.90 cdm 1.65 8.50 cdm 2.10 9.05 1.78 - 2 10 -2 1 0 usg 0.74 14.40 0.56 17 2 0 . -2 usg 2.33 14.40 -2 Figure 4.5. Continue. 74 distance (in km) are shown to the right of each seismogram. For each station, the seismograms plotted on the left panel correspond to the horizontal component parallel to the fault strike and those on the right panel are the components perpendicular to the fault strike. In Figure 4.5, the amplitude of the seismogram is normalized to its peak value, which is indicated to the right of each seismogram. We also weighted observed seismograms by a factor between 1 and 2 with a increment of 0.25, by trial error, to get a better visual correlation between the observed data and the synthetic results in addition to a better fitting in terms of the Euclidean norm. The Euclidean norm is defined by the root of the sum of the squared residuals. The smaller the value the better the fitting. In order to quantify the correlation coefficients in our error analysis, we also introduced a Euclidean norm in which squared residuals are weighted by factors of one minus the corresponding correlation coefficients. The reason for this weighted Euclidean norm is to put additional emphasis on the similarity between shapes of the observed waveforms and its synthetic counterpart. The Euclidean norm and the weighted Euclidean norm of our inversion are 16.43 and 13.35, respectively. The above result is obtained from the inversion of 16 stations within 15 km from the epicenter. In order to see how the result depends on the number of stations, we delete 7 stations, leaving only the closest 9 stations to the epicenter. The fault slip distribution determined from records at these 9 stations is plotted in Figure 4.6a and 4.6b. The corresponding comparison between observed and synthetic seismograms is shown in Figure 4.7. Comparing the two results, we found that although there are some changes in the detail by deleting 7 stations, the main features are unchanged. The Euclidean norm and the weighted Euclidean norm for this case are 12.49 and 9.01, respectively. The Euclidean norm and the weighted Euclidean norm calculated for the same 9 stations but using the solution given in Figures 4.5a and 4.5b are 12.15 and 75 Slip Amplitude Distribution o E < r > J C Q . b - 5 0 Fault Strike (km) 0000 29.10 90-20 79.J0 100.4 129.9 Figure 4.6a. The fault slip (cm) distribution of the second inversion shaded in different intensity according to scales indicated below the figure. --j O n Fault D ip (km) Slip Amplitude Distribution o in o ■ 5 5 0 Fault Strike (km) Figure 4.6b. Contour plot of the same slip distribution as in Figure 4.6a with contour interval of 17.93 cm. 1 . 2 . U S C 0.63 4.91 0.47 usc 0.56 7.61 0.37 usc 0.46 9.62 0.24 usg 1.01 3.11 0.52 usg 1.10 3.11 0.55 - 2 0 2 4 6 8 10 - 2 0 2 4 6 8 10 use 2.71 4.91 1.29 usc 1.54 7.61 2.50 usc 1.47 9.82 1.79 usg 2.09 3.11 1.89 usg 0.97 3.11 1.08 usg 0.93 7.92 usg 2.38 7.92 1.67 cdm 1 26 6.94 cdm 1.95 6.94 0.58 1.59 cdm 0.33 6.50 cdm 1.65 6.50 0.21 cdm 1.65 9.05 cdm 2.10 9.05 15 1.36 1.47 — 1 1 0 -2 0 2 6 8 1 0 -2 4 Figure 4.7. Same as Figure 4.5 but for the second inversion results. 77 9.31, respectively. The increase in the weighted Euclidean norm indicates that there are some decreases in similarity between observed and synthetical seismograms, caused by including additional data. This exercise tells us that the overall picture obtained from the records of 9 stations is approximately the same with that obtained by adding 7 additional stations. We found, however, that we need at least 5 stations to keep generally the similar pattern. In the above inversion, the observed records were not corrected for possible recording site amplification effect. In order to apply the correction, we shall follow Phillips and Aki (1986) and use the coda part of the strong motion records themselves. For the valid application of the coda method, the measurement must be made at a lapse time greater than twice the travel time of the S-waves. In our case, the distance is less than 15 km, and the S-wave travel time is less than about 10 second. We chose the lapse time of 28 second for the coda measurement which satisfies the above condition even including a likely source duration. We calculate the root mean squares of the bandpassed displacement signal over a 10 sec time window centered at the lapse time of 28 second. Strictly speaking, the coda amplification factor is a directionally averaged amplification factor for S-waves. As a preliminary step, we assume that it applies to the whole seismogram. Figures 4.8a and 4.8b are the inversion results for the fault slip distribution from the site effect corrected data. Table 4.2 lists the site correction factor used for each station. In comparing Figure 4.8 and Figure 4.4, we find no significant differences between the new result and that without the site effect correction. This is probably because most of the stations are located on similar sediments with similar amplification factors as shown in Table 4.2. The synthetic seismograms based on this inversion are 78 Table 4.2 Site Amplification Factors STATION NAM E P arallel To F ault S trik e P erp en d icu lar To F au lt Strike 1 USC 0.8971 1.0121 2 USC 1.2393 2.0538 3 USC 1.0016 1.2532 4 USC 0.6821 0.7616 5 USC 1.2467 1.2181 6 USC 1.3931 1.5444 7 USC 1.2279 0.8155 8 USGS 0.9538 1.0319 9 USGS 1.1982 1.0099 10 USGS 0.7048 0.8162 11 USGS 1.4207 1.2361 12 USGS 1.7448 1.8316 13 CDMG 0.7881 0.7026 13 CDMG 0.7824 0.6393 14 CDMG 1.3361 0.9405 16 USGS 0.7831 0.8985 79 Slip Amplitude Distribution — 0.000 27.30 34,80 81.90 109.2 136.3 Figure 4.8a. The fault slip (cm) distribution of the third inversion shaded in different intensity according to scales indicated below the figure. o o o Foult D ip (km) Slip Amplitude Distribution o o in o - 5 0 5 Fault Strike (km) Figure 4.8b. Contour plot of the same slip distribution as in Figure 4.8a with contour interval of 19.50 cm. 1 . 2 . 3 . 4 . 5 . f ( k \ 7^ U SC 0.71 +.91 0.31 U SC 0.4 5 7.61 0.38 use 0.46 9 82 0.28 usc 0.42 11.41 0.42 usc 1.81 11.62 2.13 -2 0 2 4 6 8 10 - 2 0 2 4 6 8 10 U SC 2 .6 7 +.91 U SC 0 .7 5 7.61 1.25 usc 1.17 9.82 usc 1.08 11.41 usc 2 .0 7 11.62 3.08 usc 0.58 12.03 6 0.71 usc 0.5 0 14.57 7 0.98 usg 1.05 3.11 8 usg 0.92 3.11 9 use 1-33 7.92 1 0 0.62 -2 usc 0.67 12.03 usc 0.6 3 14.57 0 .39 usg 2.03 3.11 0.99 usg 0.96 3.11 0.64 usg 2.9 2 7.92 -2 Figure 4.9. Same as Figure 4.5 but for the third inversion results. 81 11 . 1 2 . 13 . 14 . 15 . usg 2.51 10.06 1.92 usg 1.12 12.42 3.62 cdm 1.61 6.94 0.66 cdm 0.42 6.50 0.33 cdm 1.24 9.05 2.26 - 2 0 2 4 5 8 1 0 1.16 usg 1.03 12.42 0.63 cdm 2.77 6.94 cdm Z 5 8 6.50 1.65 cdm Z 2 4 9.05 1.86 -2 0 2 6 8 1 0 4 usg Z 5 9 14.40 -2 2 0 4 6 8 1 0 usg 0.95 14.40 6 7 8 9 2 0 -2 Figure 4.9. Continue. plotted in Figure 4.9 in comparison with the site effect 'corrected' data. The Euclidean norm and the weighted Euclidean norm of this inversion are 18.96 and 15.33, respectively. Comparing with the ones without the site effect correction, we find the Euclidean norm and the weighted Euclidean norm have been increased. These could be partly because the coda site effect is the azimuthal average whereas the site response of the direct S-waves is incident angle dependent, or because the phase delay terms of the site responses are neglected. It could be also due to possible non-linear effect in the site response. With the same site effect corrected data, we performed a simultaneous inversion of the rupture time and slip intensity. The resultant rupture time contour is plotted in Figure 4.10 and the slip distribution is plotted in Figures 4.11a and 4.1 lb. This is a result of only one step iteration. The rupture process on the fault becomes rather irregular, but the main feature of the fault slip distributions with two large slip zones remains unchanged from the previous results. The synthetic seismograms based on this variable rupture time and slip model are plotted in Figure 4.12. Compared to the previous result (Figure 4.9), we found that a good fit is obtained for some of the phases which could not be fit by using the slip variable model only. This indicates that the actual faulting process may be an irregular rupture process caused by the geometrical or the strength inhomogeneity on the fault. The Euclidean norm and the weighted Euclidean norm are 16.38 and 12.26, respectively. The seismic moment from the above four different inversions are listed in Table 4.3. These moment values are adjustable with different choices of the initial model since we filtered out the low frequency information. Thus, the seismic moment for this source inversion should be a constrained parameter rather than a derived one. 83 Rupture Time Contour o E a . Q 3 a L u - 5 0 5 Fault Strike (km ) Figure 4.10. Contours of rupture time model determined by the simultaneous inversion (forth) of both slip and rupture time from the Whittier-Narrows earthquake strong motion records with contour interval of 0.5 sec. 84 Slip Amplitude Distribution ° r E « n a < 5 1 jL - i 1 < - j 1 ---- 1 ---- 1 ---- 4 . ; ----l - 5 Fault Strike (km) m m s m 0.000 24.90 49.00 75.90 90.00 122.9 Figure 4.11a. The fault slip (cm) distribution of the forth inversion shaded in different intensity according to scales indicated below the figure. o o t-n Fault D ip (km) Slip Amplitude Distribution o in o - 5 0 5 Fault Strike (km) Figure 4.11b. Contour plot of the same slip distribution as in Figure 4.1 la with contour interval of 17.50 cm. 2 . 3 . 4 . $ £ } & use 0.71 +.91 0.22 use 0.+5 7.61 0.30 use 0.+6 9.82 0 .2 + use 0.+2 11.+1 0.36 use 1.81 11.62 0.83 - 2 0 2 4 6 8 10 use 2.67 +.91 2.4+ use 0.7 5 7.61 0.91 use 1.17 9.82 1.35 use 1.08 1 1.+1 use 2.07 11.62 1.73 use 0.58 12.03 6 0 .6 + use 0.50 1+.57 7 8 usg 0.92 3.11 9 0.35 usg 1.33 7.92 1 0 0.50 0 2 4 -2 8 1 0 use 0.67 12.03 use 0.63 14.57 0.+2 usg 2.03 3 .11 usg 0.96 3.11 1.01 usg 2-92 7.92 -2 Figure 4.12. Same as Figure 4.5 but for the forth inversion results. 86 usg 1.45 10.06 0.90 usg 1.03 12.42 0.85 cdm 2 7 7 6.94 1.68 cdm 2 5 8 8.50 1.76 cdm 2 2 4 9.05 1.64 -2 0 1 0 2 6 8 4 usg 2 5 1 10.06 1.12 2 5 9 cdm 1.61 6.94 0.57 cdm 0.42 8.50 0.27 cdm 1.24 9.05 15 . 1.44 -2 usg 0.95 14.40 0.53 17 . 18 2 0 . -2 usg 2 .5 9 14.40 -2 0 8 1 0 2 6 4 Figure 4.12. Continue. 87 Table 4.3 Seismic Moment of the Inversion Results M ODEL M O M EN T (dyn cm) 1st INVERSION 9.69 E+24 2nd INVERSION 1.20 E+25 3rd INVERSION 1.19 E+25 4th INVERSION 1.05 E+25 5th INVERSION 9.86 E+24 6th INVERSION 1.01 E+25 The same inversion procedure was also carried out for the fault geometry with strike at 280° and dip at 30°. A detailed study for the same fault geometry has also been worked out by Hartzell and Iida (1990), using a similar strong motion data set but based on the finite element/discrete wavenumber Green's function calculation scheme. In our study, we first determined the slip intensity distribution using the records of 16 strong motion stations. The result of the fault slip is plotted in Figures 4.13a and 4.13b. A fault slip distribution with two distinct high slip zones were again found with slight differences in their positions and a higher peak slip value as compared to the previous results. The Euclidean norm and the weighted Euclidean norm are 16.34 and 13.08, respectively. The synthetic seismograms based on this result were compared with the observed ones in Figure 4.14. In this figure, we found a larger discrepancy in the peak displacement for the component parallel to the fault strike over several stations (e.g. station 1, 2, 3, 4, and 16, etc.) as compared to the previous case. For some stations, this component is mainly composed of the SH waves which appears to be very sensitive to the fault orientation. Since our previous fault strike direction is adopted from the first motion study, this may indicate that the short period first motion result provides a better control over the fault orientation than that derived from the long period data. We also inverted the site effect corrected data. The results of the determined fault slip distribution based on the ’ site corrected' seismograms are plotted in Figures 4.15a and 4.15b. The comparison between the synthetic results and the observations is plotted in Figure 4.16. The Euclidean norm and the weighted Euclidean norm are 17.18 and 13.30, respectively, showing a slightly better fitting than the corresponding inversion in our previous case. 89 o Slip Amplitude Distribution —i ------- » r 1 ------ j — Fault Strike (km) 0.000 20-12 52.24 78,55 104.4 150.6 Figure 4.13a. The fault slip (cm) distribution of the fifth inversion shaded in different intensity according to scales indicated below the figure. 's O o Fault D ip (km) Slip Amplitude Distribution o o - 5 0 5 Fault Strike (km) Figure 4.13b. Contour plot of the same slip distribution as in Figure 4.13a with contour interval of 18.66 cm. 0.27 0.18 use 0 .4 4 9.82 0.12 use 0.35 11.41 0.13 use 2.51 11.62 1.48 X -2 1.53 7.61 1.21 use 1.51 9.82 1.97 use 0 .8 4 11.41 0.90 use 2 .2 8 11.62 2 1 9 -2 0 2 4 6 8 1 0 use 0.69 1 2 0 3 use 0.59 14.57 0.31 usg 1.33 7.92 1 0 . 0.58 -2 use 1.12 1 2 0 3 use 0 .4 3 14.57 usg 2 1 0 3.11 0.85 usg 2 1 9 7.92 -2 Figure 4.14. Same as Figure 4.5 but for the fifth inversion results. 91 1 1 1 2 . 13 . 14 . 15 . ^/VAy^ywv usg 3.70 10.06 1.68 usg 1.74 12.42 3.57 cdm 1.56 6.94 0.69 cdm 0.49 8.50 0.28 cdm 1.48 9.05 -2 8 10 usg 1.29 10.06 1.24 usg 2 .1 3 12.42 1.03 cdm 1.78 6.94 1.42 cdm 1.63 8.50 1.93 cdm 2.19 9.05 1.84 -2 0 2 8 4 6 10 usg 1.08 14.40 6 7 18 9 20 8 6 - 2 0 2 4 10 usg 2.17 14.40 -2 0 2 4 6 8 Figure 4.14. Continue. Slip Amplitude Distribution i — . i i i i------ 1 ------ 1 ------ 1 ------ 1 i i - 5 0 5 Fault Strike (km) 0.000 32.08 8416 96.24 128 3 160.3 Figure 4.15a. The fault slip (cm) distribution of the sixth inversion shaded in different intensity according to scales indicated below the figure. v £ ) Fault D ip (km) Slip Amplitude Distribution o i f > o - 5 5 0 Fault Strike (km) Figure 4.15b. Contour plot of the same slip distribution as in Figure 4.15a with contour interval of 22.91 cm. 2 3 4 5 - 2 0 2 4 6 8 1 0 -2 0 2 4 6 8 10 x/vyv/\/A - w V W \ u sc a 3 6 7 6 1 A a A a ^ - u s c 0 + 1 9 8 2 ^ W V V A - usc 0.50 11.41 0.15 usc 1.81 11.62 1.34 use 2 .5 4 4.91 1.79 U SC 0 .8 0 7.61 1.18 \ A tf V w ^ u " c ,i+ 9 8 2 x A/Wa— use 1.08 11.41 1.02 use 2 .0 5 11.62 2 .0 3 use 0.50 1 2 0 3 6 usc 0.49 14.57 7 0.77 8 usg 0 .9 4 3.11 9 usg 1.88 7.92 10 0.70 - 2 usc 0.81 1 2 0 3 usc 0.51 14.57 usg 0 .9 5 3.11 0.58 usg 2 7 4 7.92 2 4 - 2 0 6 8 10 Figure 4.16. Same as Figure 4.5 but for the sixth inversion results. 94 11. 1 2 . 13 . 14 . 15 . - m — -2 usg 2.66 10.06 1.56 usg 1.06 12.42 3.36 cdm 2.02 6.94 0.81 cdm 0.59 6.50 0.40 cdm 1.12 9.05 1.55 J___ I ___ I 6 8 10 usg 1.05 10.06 1.25 usg 1.13 12.42 cdm 2.38 6.94 1.45 cdm 2.65 3.50 2.37 cdm 2.26 9.05 1.52 -2 usg 1.46 14.40 16 7 8 9 20 2 4 6 8 10 - 2 0 usg 2.32 14.40 2.28 — 2 Figure 4.16. Continue. 4.4 Discussion and Conclusion Two best-fit fault geometries, with strike direction and dip angle of 270° and 35°, and 280° and 30°, respectively, have been chosen from many fault geometry configurations. The fits provided by these two fault geometries with constant rupture velocity are about the same. The first model with strike of 270° and dip of 35° is preferred if we insist on a pure thrust fault. The fit can be improved, however, for the second model, if we introduce an additional strike slip component. Most of our study was based on the strong motion records from 16 stations selected within 15 km from the epicenter. The sufficiency of data was tested by eliminating 7 stations from the above 16 stations, and finding that the resultant slip distribution differs from the previous one only in the detail, and the overall picture of the fault slip distribution remained the same. This test tells us that the waveform inversion of a not so large number of records in an epicentral area well-constrains the earthquake faulting process. In order to obtain a rough picture of the source mechanism, a little more than 5 well-distributed stations may be sufficient. A site effect correction based on the coda waves of the same strong motion records has been performed. We did not find significant improvement or change between the results with and without this site effect correction, possibly because most stations are located on similar sediment sites and the amplification factor did not show a great variation from station to station. Further studies are needed to include the correction of the phase delay terms due to the site response, for example, using a simple one dimensional model to predict the possible phase response of the soft layer sediment site. In addition, the non-linear site response should also be taken into consideration and its effect on the strong ground motion records need to be clarified. 96 All the inversion results of this study indicate that there are two regions of major slip on the fault plane with large amount of seismic moment release. One is located immediately below the hypocenter. The other is near the lower left (east) comer of the fault. In addition, there are three or four other smaller localized source areas around the hypocenter. Together they form a highly heterogeneous fault slip distribution for the Whittier-Narrows mainshock. In fact, all our inversion results show a 10 km by 4 km regions of high seismic slip, from which the great majority of seismic energy is radiated. The inversion result with variable slip and rupture time model gives the best fit in terms of both the Euclidean norm and the weighted Euclidean norm indicating that the rupture velocity is not constant, and some directions faster than others. But we find no simple causal relation between these rupture process with the fault slip distributions. This conclusion is consistent with that of Hartzell and Iida (1990). Based on coseismic geodetic observations, Lin and Stein (1989) studied the coseismic folding and source mechanism of the Whittier-Narrows earthquake. They modeled 214 corrected geodetic observations with a simple half space model. Since their data have no control over the fault strike orientation, they adopted the geometry from Hauksson and Jones (1989). The best fitting solution they obtained without constraint from the aftershock data has a fault length of 12 km and width of 4 km, a thrust faulting mechanism with an average slip of 71 cm, and a dip angle of 34°. This result is in an excellent agreement with our inversion. The depth to the fault top obtained in their model is about 12 km. As pointed out by Xu and Su (1988) in their geodetic data modeling, a homogeneous half space model would result in a shallower source depth than that of a heterogeneous half space model where seismic velocity increases with depth, especially when a very soft sediment layer is present. By taking a more realistic earth model for the Los Angeles basin, the fault depth obtained by Lin 97 and Stein (1989) would be effectively increased, possibly resulting in an even better agreement between their geodetic inverse solution and our strong motion inverse results. An interesting result we have found is the negative correlation between the distribution of mainshock slip and the distribution of aftershocks. Figure 4.17 is a replot of Figure 4.8. Figure 4.18 is the aftershock distribution reproduced from Hauksson and Jones (1989). Figure 4.17 shows that the mainshock slip is largely concentrated below its hypocenter. On the other hand, Figure 4.18 indicates that most of the aftershocks occurred above the mainshock hypocenter. Furthermore, a large number of aftershocks with magnitude above 2 were occurred at the upper left (east) side and upper right (west) side of the hypocenter on the fault plane where there was small amount or even no slip during the mainshock faulting. This suggests that the area undergoing large amount of slip during the mainshock becomes a low stress zone and the stress concentration is shifted to the unbroken area, or the mainshock slip deficiency area. This observation is consistent with the fault barrier model proposed by Aki (1979). According to this barrier model, the mainshock rupture propagates, leaving behind unbroken patches call 'barriers' where aftershocks occur as the release of the stress concentration around barriers through static fatigue. The same negative correlation between mainshock slip and aftershock distribution has been observed by many other investigators (Mendoza and Hartzell, 1988; Beroza and Spudich, 1988; Schwartz et al., 1989; Engdahl et al., 1989; etc.). As pointed out by Aki (1984), identification of fault barriers and asperities may be useful in the interpretation of characteristic earthquakes and in quantitative prediction of the strong ground motion. Thus the understanding and characterizing of fault barriers and asperities from the strong motion data is a very important subject. 98 Fault D ip (km) ' 1 t; • v j . ; liilfll Slip Amplitude Distribution - 5 0 Fault Strike (km) Slip Amplitude Distribution - 5 0 Fault Strike (km) 0000 27.30 34.60 81.90 109.2 136.3 Figure 4.17. A replot of Figure 4.8. v O v O •< 6 ' 0.0 + 5* 2.0 + 3.0 + M l = 5.9 Mainshock October 1 4.0 + 4 ‘ 3' 2 ' 2 KM 3 4 * A ' B - 6 East North West South - 7 - 9 - 1 5 -1 8 1 1 DISTANCE (KM ) DISTANCE (KM) Figure 4.18. (a) Map, (b) south-north cross section, and (c) front view from the south of the mainshock focal mechanism and following aftershocks unitil 10:58, October 4, 1987, which was the time of the largest aftershock (reproduced from Hauksson and Jones, 1989). 100 Our inverse result of the fault slip distribution with and without the site effect correction is similar in most respects with that obtained by Hartzell and Iida (1990). Although we used quite different methods, we both found two major slip zones surrounded by several other small localized slip sources. The peak slip in both studies are about the same. By allowing a strike slip component in their study, Hartzell and Iida (1990) found a preferred fault geometry with strike of 280° and dip of 30°. We restricted our slip to a pure thrust, which resulted in a preferred fault geometry with strike of 270° and dip of 35°. In general, our results based on the ray method are consistent with that obtained by Hartzell and Iida (1990) based on the more time consuming discrete wavenumber method. Thus, our method should provide an economic way in both the forward synthetics of the strong ground motion and the inverse study of the earthquake source process as long as the ray method is applicable. 101 Chapter V APPLICATION TO THE LOMA PRIETA EARTHQUAKE 5.1 Description of the Problem The Loma Prieta earthquake struck the San Francisco Bay area on October 17, 1989 and was the largest earthquake to occur on the San Andreas fault since the 1906 San Francisco earthquake. This earthquake of magnitude 7.1 ruptured a 40 km segment of the San Andreas fault in the southern Santa Cruz Mountains of northern California. The mainshock hypocenter was at 37.04°N, 121.88°W, at a depth of 17.6 km (Dietz and Ellsworth, 1990). It resulted in 62 direct fatalities and over $6 billion in property damage over a widespread area from Monterey Bay to San Francisco Bay. Understanding the dynamic faulting of the Loma Prieta earthquake is of great importance in studying the San Andreas fault activity and prediction of the strong ground motion for future earthquakes. With the large amount of strong motion data available, it provides us another good opportunity to study the earthquake source process and its physical implication. 5.2 Data Collection and M odeling The Loma Prieta earthquake has been recorded by many strong motion accelerographs operated by the California Strong Motion Instrument Program of the California Division of Mines and Geology (CDMG). The seismograms from these strong motion accelerographs have been processed and are available to us for analysis. In order to minimize the propagation path effect, we choose to use only stations within about 34 km of the epicenter. The reason for choosing this distance range is the 102 Radial T r a n s v e r s e Vertical ( 13) 0.13E-03 ( 15) 0.13E-03 ( 19) 0.88E-04 ( 23) 0.69E-04 ( 27) 0.55E-04 ( 32) 0.32E— 04 ( 37) 0.JO E— 04 (41) 0.23E-04 ( 45) 0.1 B E — 04 ( 51) 0.12E— 04 ( 58) 0.12E-04 ( 61) 0.1 B E— 04 ( 66) 0.16E-04 ( 71) 0.22E-04 ( 75) 0.20E— 04 ( 80) 0.12E— 04 ( 85) 0.20E— 04 ( 90) 0.24E— 04 ( 95) 0.2IE-04 (100) 0.22E-04 (105) 0.17E-04 (110) 0.15E-04 (115) 0.1 IE-04 (120) 0.91E-05 (125) 0.93E-05 (130) 0.68E-05 (135) 0.87E-05 (140) 0.55E— 05 0 5 Time (s e c )-X /3 .5 s ScS A I -y -- SmS 0.14E-03 0.20E-04 0.29E-04 0.40E-04 0.39E-04 0.35E-04 0.24E-04 0.17E-04 0.1 B E — 04 0.22E-04 0.28E-04 0.18E-04 0.24E-04 0.3 IE-04 0.32E— 04 0 .3 7 E -0 4 0.2BE-04 0.28E-04 0.20E-04 0.13E-04 0.90E-05 0.88E-05 0.92E— 05 0.13E-04 0 10 5 SmS ScS 0.44E-04 0.61E-04 0.49E-04 0.17E-04 0.19E-04 0.13E-04 0.12E-04 0.7BE-05 0.5BE-05 0.75E-05 0.55E-05 0.51 E-05 0.81 E-05 0.7BE-05 0.61 E-05 0.73E-05 0.70E-05 0.87E-05 0.70E-05 0.52E-05 0.46E-05 0.40E-05 0.29E-05 0.33E-05 0 5 10 Time (s e c )-X /3 .5 Time (s e c ) -X /3 .5 Figure 5.1. Three component synthetic ground velocities (cm/sec) produced by a dislation point source at 12 km depth with the focal mechanism as that of the 1989 Loma Prieta earthquake and a seismic moment of 1020 dyne-cm. The source time function is a ramp function with 0.5 sec rise time. The hypocentral distance inside the parenthesis is in km. The reduction velocity is 3.5 km/sec. following. Somerville and Yoshimura (1990) and Chin and Aki (1991) calculated synthetic seismograms (see Figure 5.1) at various locations for a point dislocation source with the focal mechanism of the Loma Prieta earthquake buried at a depth of 12 km in a 11-layer crustal model, which was derived by averaging the velocities on both sides of the San Andreas fault. Three principle S-wave arrival time curves are shown in Figure 5.1; the direct S-wave (S), the reflected S-wave from the Conrad discontinuity (ScS) and the reflected S-wave from the Moho discontinuity (SmS). The hypocenter distance (in km) is plotted inside the parenthesis to the right of the radial components, followed by the peak amplitude of the trace. Their results demonstrated that for epicentral distances less than 30 km, the onset of the maximum peak velocity coincides with the arrival time of the direct S-waves. At distances between 30 and 50 km, the reflected waves from the Conrad discontinuity start to contribute to the direct S- waves. For distances greater than 50 km, the largest synthetic ground velocities come from the Moho reflection. Since our forward and inverse methods are limited to the direct S-waves, we are restricted to using stations closer than 50 km epicentral distance even if we ignore the weak Conrad reflected waves. By considering the extent of the mainshock faulting, we then choose stations with epicentral distances less than about 34 km. For these distances, records are available from 9 strong motion accelerograph stations. The names and locations of these stations are listed in Table 5.1. A map view of the station distribution is plotted in Figure 5.2 with the mainshock marked by a star. As demonstrated in Chapter IV, this quantity and distribution of stations should provide a relatively good image of the earthquake source process. Again we only modeled the horizontal components of the strong motion seismograms in our inverse process due to the reasons pointed out in Chapter IV. The instmment-response corrected ground acceleration records have been integrated twice to obtain the displacement seismograms. 104 Table 5.1 Strong Motion Station Name and Location Station Latitude longitude 1 CORRALITOS - Eureka Canyon Road 37.046 121.803 2 CA PITO LA - Fire Station 36.974 121.952 3 SANTA CRUZ - UCSC/Lick Lab. 37.001 122.060 4 W ATSONVILLE - 4-Story Commercial BLDG 36.909 121.756 5 LEX IN TO N Dam 37.202 121.949 6 SARATOGA - Aloha Ave. 37.255 122.031 7 GILRO Y #2 - HWY 101/Bolsa RD. Motel 36.982 121.556 8 COYOTE Lake Dam - Downstream 37.124 121.551 9 G ILRO Y #4 - San Ysidro School 37.005 121.522 105 S t r o n g M otion S t a t i o n s for L om a P r ie t a E a r t h q u a k e ♦ 6 SARATOGA ♦ 5 LEXINGTON 8 COYOTE 7.1 MAINSHOCK CORRALITOS ♦ 9 GILROY 4 ♦ 7 GILROY 2 ♦ 3 SANTA CRUZ PITOLA ♦ 4 YTATSONV1 10 KM 122 ° 50' 40' 30’ Figure 5.2. Map view of the strong motion station distribution used for the Loma Prieta earthquake source inversion study. 106 This displacement data were bandpass filtered from 0.8 to 4.5 Hz with a highpass transition band of 0.4 Hz to provide the desired frequency range valid for our asymptotic Green's function calculation and the results were interpolated at a uniform time interval of 0.08 second. Only few stations have their triggering times available which are used as references for timing of their S-wave arrivals. For stations without triggering times, we first pick the onset point of their S-wave arrivals and then calculate its arrival time according to our velocity model and the earthquake hypocentral distance. Our velocity model for the Santa Cruz mountain area is a smoothed version of the model used by Somerville and Yoshimura (1990). A profile of this model is plotted in Figure 5.3, where the solid line represents our smoothed velocity model and the dashed line corresponds to the original model of Somerville and Yoshimura (1990). We have removed the thin near surface sedimentary layer in the original model and simulate the average site response model using the strong motion coda waves. The mainshock fault is modeled as a 40 km by 14 km planar region with a strike of 130° and dip angle of 70° (USGS first motion mechanism, 1990). This planar region is then divided into 41 by 15 equally spaced grid points. We model the rupture process with two different constant rupture velocities, namely, 70 and 80 percent of the average local S-wave velocity . For simplicity, the fault slip is constrained to a rake direction of 130°. The initial model of fault slip distribution is plotted in Figure 5.4. It shows a high central plateau that gradually tapers to zero at the fault edges. Again, we take a uniform time function over the entire fault of the form [H(t-tr)-H(t-th)]/V t - tr, where tr is the fault rupture time and th is its healing time. The rise time defined by th-tr is set to be 0.2 sec. The synthetic seismograms are bandpass filtered between 0.8 Hz and 4.5 Hz, same as that for the observational data. 107 D e p th (km) Velocity Structure for L om a Prieta 2 4 6 8 o o C M o C O 8 6 2 4 Velocity (km/sec) Figure 5.3. Seismic velocity model for the Santa Cruz mountain area. 108 Slip Amplitude Distribution Fault Strike (km) 0.000 32.53 65.06 97.60 130.1 162-6 Figure 5.4. Intial fault slip distribution shaded in different intensity with a high central plateau that gradually tapers to zero at the fault edges. 5.3 Inverse Result and Analysis Three different inversions of the strong motion waveform data have been preformed. In the first inversion, we assume a bilateral rupture model with constant rupture velocity of 80 percent of the average local shear wave velocity. The rupture time contour is plotted in Figure 5.5. Our velocity structure may correspond to a rock site model. In order to simulate a velocity model appropriate for the "average" site, we multiply the Green's function calculated based on our velocity model by the site amplification factor averaged over all stations used for the analysis relative to the Santa Cruz site. The resultant fault slip distribution is plotted in Figure 5.6a. In this figure, we shaded the fault slip amplitude with different background intensities varying from black to white with dark regions representing high slip areas. Figure 5.6b is a conventional contour plot of the same slip distribution. These figures show that the slip on the fault is highly variable from place to place. There is one large high slip zone on the right (southeast) side of the fault plane with peak slip of 274.8 cm. On the left (northwest) side of the fault plane, there is also a large high slip zone segmented into a few regions of peak slip. The Euclidean norm and the weighted Euclidean norm, defined in our previous chapter, are 17.01 and 14.46, respectively. The synthetic seismograms based on this determined fault slip model are plotted in Figure 5.7 and are compared to the observations plotted on top of each synthetic counterpart. The vertical axis corresponds to the station number listed in Table 5.1. The peak displacement and epicentral distance are plotted to the right of each seismogram. Seismograms plotted on the left correspond to the horizontal component parallel to the fault strike and those to the right are for the component perpendicular to the fault strike. In this figure, we have normalized each seismogram by its peak value. We also weighted each observed seismograms by a factor between 1 and 2 with a increment of 0.25, by trial and error, 110 Fault D ip (km) Rupture Time Contour o ID o -20 -1 5 -1 0 5 15 20 0 5 10 Fault Strike (km) Figure 5.5. Contours of rupture time model for the Loma Prieta earthquake with constant rupture velocity of 80 percent of the average local S-wave velocity and contour interval of 0.5 sec. Fault D ip (km) fcuit D ip (km) Slip Amplitude Distribution O -----1 -----1 ----- 1 -----t-----1 -----1 -----« -----1 -----1 ------1 — I I 1 -----1 ----- 1 -----1 -----» ----- r---- 1 ----- 1 ----- 1 -----1 -----1 -----1 ----- 1 ----- 1 -----1 -----1 -----1 ----- 1 -----1 ----- 1 -----1 -----1 ----- X-----1 -----1 -----1 -----r --20 -1 5 --10 - 5 0 5 10 15 20 Fault Strike (km) 0.000 5*9C 109.9 104.8 219.8 2H.8 Figure 5.6a. The fault slip (cm) distribution of the first inversion shaded in different intensity according to scales indicated below the Figure. Slip Amplitude Distribut’on o to o 0 5 10 15 20 Fault Striko (km) F i g u r e 5 . 6 b . C o n t o u r p l o t o f t h e s a m e s l i p d i s t r i b u t i o n a s i n F i g u r e 5 . 6 a w i t h c o n t o u r i n t e r v a l o f 4 5 . 8 0 c m . 1 1 2 19.75 +.03 9.32 2.15 1.12 16.58 0.67 2.3+ 17.97 3.02 7.98 19.02 7.1 + - 3 0 3 6 9 12 15 18 21 +.07 6.9+ 1 6.96 +.57 9.32 2 3.02 1.31 16.58 3 0.88 +.7+ 17.97 4 1.52 2.6+ 19.02 5 0.73 - 3 0 9 12 15 18 21 6 3.93 27.25 7.22 3.58 0.61 30.58 0.58 +.07 31.97 2.37 - 3 0 9 12 15 18 21 3 6 2.82 27.25 6 2.71 29.+6 / +.+7 2.62 30.58 8 1.05 2.50 31.97 \JV J 9 2.57 10 9 12 15 18 21 - 3 0 3 6 Figure 5.7. Plot of seismograms for the first inversion with the observed traces on top of the synthetic ones for each station numbered along the vertical axis. Each seismogram is normalized to its peak values (cm), indicated to the right the trace followed by its epicentral distance (km). Seismograms to the left are components parallel to the fault strike and those to the right are perpendicular to the fault strike. 113 to get better correlations between the observed data and the synthetic results in addition to a better fitting in terms of the Euclidean norm. By examining the synthetic seismograms for each station, we found a large impulse like displacement response for station Corralitos (Station #1) where the observation shows multi-arrival phases. We increased the weighting parameter for this station by a factor of 2 and reperformed the inversion. But the result is not significantly improved. We then replaced the rupture model by a constant rupture velocity equaling 70 percent of the average local shear wave velocity. The rupture time contour for this model was plotted in Figure 5.8. The time step between neighboring lines of this contour map is 0.5 sec. The determined fault slip model is plotted in Figures 5.9a and 5.9b. Considerable differences are shown in the slip distribution between the new result and that of the preceding one. We still found two large high slip areas on the right (southeast) side and the left (northwest) side of the fault plane, respectively, with maximum slip of 356.3 cm. But the high slip area on the left (northwest) side of the fault plane has been segmented into three small zones of concentrated slip. The high slip area on the right side of the fault plane has been moved upward. The Euclidean norm and the weighted Euclidean norm are 14.90 and 11.20, respectively. Compared to the preceding case, we found improvements not only in terms of the residuals but also in terms of the waveform correlations between the synthetics and the observations. Figure 5.10 shows the comparison between the synthetic seismograms, based on this determined slip model, and the observed ones. In this figure, we see clearly the improvement in the synthetic seismogram for the Corralitos station. Although we are still unable to explain much of the observation at the Santa Cruz station, the similarities between the observations and the synthetics for the rest of the stations are quite encouraging. 114 Fault D ip (k m ) Rupture Time Contour o in o 20 10 15 -20 5 -1 5 - 5 0 -1 0 Fault Strike (km) Figure 5.8. Contours of rupture time model for the Loma Prieta earthquake with constant mpture velocity of 70 percent of the average local S-wave velocity and contour interval of 0.5 sec. F ault D i p (km ) Fault D i p (Km) Slip Amplitude Distribution i — i— i — i — |— i — i— i— i — |— i— r — - t — i — |— i — I i —'— r r|" P rF : • y * • : . : . . |i, /• iijjf • i j j l i • L|| • ; i j j • • , " J; - . ,1 - ' M b w m is-j:.-.!.: j: I I I ! ! ! : } : ' ■ o k - i t w f M ! w > i ? > "--'- r * ^ - ■ -- '• : v : ; j t I t I I 1 ----- (----- I -----1 ------1 ------1 ------1 I I I I I ------1 — L I ■ I ■ ■ I 1 t I L -1 5 -1 0 - 5 0 5 Fault Strike (km) 1 0 15 20 0 0O0 7t28 " « 5 2*3.8 285.! 31B.3 Figure 5.9a. The fault slip (cm) distribution of the second inversion shaded in different intensity according to scales indicated below the figure. S'ip Amplitude Distribution o o o 15 20 1 0 0 5 Fault Strike (km) F i g u r e 5 . 9 b . C o n t o u r p l o t o f t h e s a m e s l i p d i s t r i b u t i o n a s i n F i g u r e 5 . 9 a w i t h c o n t o u r i n t e r v a l o f 5 9 . 4 0 c m . 116 1 . 2 . 3 . 4 . 5 . w jl/w V \ [ \ / \ a ~ i/U > - 1 -----I -----» » ' +.07 6.9+ 2.7+ +.57 9.32 2.91 1.31 16.58 0.72 +.7+ 17.97 2.02 2.6+ 19.02 1.79 t t . 1 . +.+8 6.9+ 6.00 +.03 9.32 +.50 1.12 16.58 1.00 2.3+ 17.97 3.06 7.98 19.02 6.50 - 3 0 3 6 9 12 15 18 21 J -----------1 -------- 1 I I ______ I _____ I _____ I _____L -3 0 3 6 9 12 15 18 21 2.82 27.25 2.71 29.+6 +.79 2.62 30.58 0.86 2.50 31.97 2.79 10 . - 3 0 3 9 12 15 18 21 6 3.93 27.25 2.06 +.79 29.+S 3.95 0.61 30.58 +.07 31.97 1.82 - 3 0 3 9 12 15 18 21 6 Figure 5.10. Same as Figure 5.7 but for the second inversion results. 117 The observed ground motion data shows a relatively low peak displacement at the Santa Cruz station. As mentioned before, the Santa Cruz station is located on a well defined rock site. If we take the rock site as our reference site, the high peak displacement associated with other stations should be due to the sediment site amplification effect. In order to take into account the site effect, we perform a similar 'site effect correction' as described in Chapter IV for the Whittier-Narrows earthquake. In this case we take the lapse time for coda measurement at a 40 second for all the stations. We then calculate the root mean squares of the strong motion coda wave in a 10 second time window centered at the above lapse time. The amplification factors are listed in Table 5.2. Unlike the Whittier-Narrows earthquake, these amplification factors vary substantially among stations. We found an average amplification factor of about 4 to 5 for the sediment site stations we considered. We checked this average amplification factor in the same region and frequency range obtained from weak motion studies (Phillips and Aki, 1986; Su et al., 1991) and found that our result is consistent with the weak motion solution. We then corrected all the other seismograms based on the amplification factor assuming the Santa Cruz station as our reference rock site. Based on the corrected seismograms, we performed our recursive stochastic inverse for the mainshock fault slip distribution. Figures 5.1 la and 5.1 lb show the fault slip distribution of this inversion result. We can see that the previous large slip regions are now narrowed down to two concentrated high slip zones with peak value of 457.9 cm. On the left (northwest) side of the fault plane, we also note an upward shift in position of the high slip zones. The Euclidean norm and the weighted Euclidean norm for this case are 16.07 and 12.84, respectively. The synthetic seismograms of this inversion are plotted in Figure 5.12 along with the observations. The similarity between observed data and the synthetics is improved for stations Santa Cruz and Capitola 118 Table 5.2 Site Amplification Factors STA TIO N NAM E Parallel To Fault S trik e P erp en d icu lar To F au lt S trike 1 CORRALITOS 4.1406 1.7243 2 CAPITOL 5.8023 5.1368 3 SANTA CRUZ 1.0000 1.0000 4 WATSONVILL 5.6281 7.7695 5 LEXINTON 3.2157 2.8836 6 SARATOGA 4.3383 5.3108 7 GILROY #2 4.6063 2.1718 8 COYOTE 2.7125 3.1797 9 GILROY #4 11.738 6.6644 119 Fauit D i p (km ) Fault D i p (km) Slip Amplitude Distribution - 2 0 - 1 5 - 1 0 - 5 0 5 Fault Strike (km) Figure 5.1 la. The fault slip (cm) distribution of the third inversion shaded in different intensity according to scales indicated below the figure. Slip Amplitude Distribution o JO o - 2 0 -1 5 - 1 0 - 5 0 5 1 0 15 20 Fault Strike (km) F i g u r e 5 . 1 1 b . C o n t o u r p l o t o f t h e s a m e s l i p d i s t r i b u t i o n a s i n F i g u r e 5 . 1 1 a w i t h c o n t o u r i n t e r v a l o f 7 6 . 3 2 c m . 120 2. B O 6.94 2.45 0.78 9.32 1.12 16.58 0.58 0.30 17.97 0.45 2.77 19.02 2.59 - 3 Q 0.98 6.94 0.69 0.79 9.32 2 0.47 1.31 16.58 3 0.61 0.84 17.97 4 0.41 0.82 19.02 5 0.52 0.65 27.25 6 0.39 0.59 29.46 7 0.47 0.97 30.58 8 0.21 31.97 9 0.42 0 15 18 21 - 3 0 0.74 27.25 0.64 2.20 29.46 1.00 0.19 30.58 0.12 0.61 31.97 0.45 - 3 0 15 8 21 6 1 Figure 5.12. Same as Figure 5.7 but for the third inversion results. 121 Table 5.3 Seismic Moment of the Inversion Results M ODEL M O M EN T (dyn cm) 1st INVERSION 1.29 E+26 2nd INVERSION 1.67 E+26 3rd INVERSION 1.75 E+26 122 compared to the results without the site effect corrections. The fitting for the other stations are comparable to the results plotted in Figure 5.10 except for the first component in displacement at the Gilroy station, which is responsible for the increase in the Euclidean norm and the weighted Euclidean norm. The seismic moment for these three inversions are listed in Table 5.3. 5.4 Discussion and Conclusion We performed three different types of fault slip inversions with constant rupture velocity. In the first case, we assumed a constant rupture velocity of 80 percent of the average local shear wave velocity. In the other two cases, we modified this rupture velocity to be 70 percent of the average shear wave velocity. The inversion result showed smaller values of both the Euclidean norm and the weighted Euclidean norm for the 70 percent case. This result is somewhat different from the case of the Whittier- Narrows earthquake in which we found an average rupture velocity to be about 80 percent of the local shear wave velocity. The site response for the Loma Prieta earthquake inferred from the strong ground motion coda waves shows a strong variation across the Santa Cruz mountain area. Comparing the rock site underneath the Santa Cruz station and the sediment site underneath the Capitola station, for example, we found about a 4 to 5-fold amplification at the sediment site. Our site responses obtained from the strong motion coda are consistent with the result from weak motion coda waves (Phillips and Aki 1986; Su et a i, 1991). Using the coda amplification factors, we corrected the whole strong motion displacement seismograms for the site amplification effect. The inversion result based on these site effect corrected seismograms shows much improvement in terms of the fitting between the synthetics and the observations except for the first displacement 123 component at the station Gilroy #4. In comparison with the well located aftershock distribution, we found the fault slip distribution obtained from the site effect corrected data better correlates with the aftershock deficiency than that without the correction. The Loma Prieta earthquake provided us another opportunity to study the relation between the mainshock slip distribution and the aftershock patterns. As pointed out in Chapter IV, it has been noted by many investigators that aftershocks tend to occur in places with large mainshock slip deficiency. This negative correlation between mainshock slip distribution and the aftershock activity can be well explained by the fault barrier model proposed by Aki (1979). In Figures 5.13a, 5.13b, and 5.13c, we overlay the projection of well-located aftershocks (U. S. Geological Survey staff, 1990) as well as the large early aftershocks (Simila, et aL, 1990) onto our fault slip contours. This projection is obtained from Beroza (1991) who used only the well located events with magnitudes greater than 1.0 occurring within 5 km of the fault plane during the month of October, 1989. These figures provide additional evidence supporting our earlier conclusion about the negative correlation between the mainshock slip distribution and the aftershock activity pattern. Upon closer examination of these figures, we found stronger negative correlation for the inversion results using rupture velocity of 70 percent of the average local shear wave velocity. In addition, our result with site effect correction provides a slightly better correlation between high slip distribution during the mainshock and low aftershock activity. This may support the validity of our site effect correction method using the strong motion coda waves for the interpretation of displacement record. Chin and Aki (1991) found that the inversion of strong motion data corrected by the weak motion site effect provides poorer synthetic prediction. We think it could be due to the fact that they used the more high frequency 124 Slip Amplitude Distribution 20 -1 0 0 10 20 Fault Strike (km) Figure 5.13. Overlain of the projection of well-located aftershocks (U. S. Geological Survey staff, 1990) as well as the large early aftershocks (Simila, et al., 1990) onto the fault slip contours of the first inversion. Fault D ip (km ) Slip Amplitude Distribution o m o -2 0 10 0 Fault Strike (km) 10 20 Figure 5.14. Same as Figure 5.13 but for the slip contours of the second inversion. Is) Os Fault D ip (km ) Slip Amplitude Distribution o m o -2 0 -1 0 0 10 20 Fault Strike (km) Figure 5.15. Same as Figure 5.13 but for the slip contours of the third inversion. H -1 to acceleration data in their investigation, which is more likely to be subjected to the non linear site effect. The Loma Prieta earthquake occurred at the Santa Cruz mountain area, one of the slightly bending sections of the San Andreas where the fault plane has a non-vertical dip. A slight oblique slip could be expected at this area. But such a large thrust slip component for the mainshock was a surprise for seismologists. By allowing the rake to vary over the fault plane, Beroza (1991) found an oblique slip on the left (northwest) section of the fault and a strike slip on the right (southeast). The same phenomenon is also suggested by the preliminary aftershock study (Oppenheimener, 1990) in spite of a great diversity of mechanisms. Although we fixed the rake in our fault model, the slip pattern from our inversion results seems dividing the mainshock fault plane into two distinctly different parts, one on the left of the hypocenter and the other on the right. It appears that the Loma Prieta earthquake rupture was initiated at the offset of the fault at depth and propagated bilaterally along both sides of offset with different ways of slip. Detailed investigation of the aftershock focal mechanism may provide better evidence for this possibility. Similar phenomenon was also found for other earthquakes (Bakun, King and Cockermam, 1986) suggesting that earthquakes may initiate at offsets or bends and that fault heterogeneity controls the earthquake faulting dynamics. 128 Chapter VI MAPPING OF THE HIGH FREQUENCY SOURCE RADIATION 6.1 Introduction With the rapid development in both instrumentation and analysis methods, the study of earthquake source processes using strong motion seismograph records has received much attention from both seismologists and earthquake engineers. There are two distinct approaches pursued in the study of earthquake sources. One, which may be called deterministic modeling, utilizes the waveform inversion algorithm to determine earthquake source parameters by minimizing the difference between the observed seismogram and the predicted ones on the basis of model. The first successful simulation of strong ground motion was done by Aki (1968) for the station No. 2 record of the 1966 Parkfield earthquake using the homogeneous unbounded elastic model. Inversion procedure to map the fault slip and rupture time distribution for a fault in more complex media were later developed (Olson and Apsel, 1982; Hartzell and Heaton, 1983; Kikuchi and Fukao, 1985; Fukuyama and Irikura, 1986; Yoshida, 1986; Takeo, 1987; Matsu'ura and Hasegawa, 1987; Beroza and Spudich, 1988). The deterministic methods provide us with a direct physical description of the earthquake source process, but their application is limited to relatively low frequencies data (f< 4.0 Hz). In the second approach, which may be called stochastic modeling, small number of source parameters which characterize the statistical properties of the source process are estimated directly from high frequency accelerograms. An example of such a model is 129 the specific barrier model of Papageorgiou and Aki (1983) constructed for a quantitative description of the earthquake faulting. In this chapter, we will develop a model based on a combination of the deterministic and stochastic processes. We use the wave energy information to obtain the pattern of high frequency radiation on the fault. The wave energies are additive at the receiver if we assume that the slip responsible for the high frequency generation is incoherent between neighboring patches on the fault. Using the isochron method, we can determine the high frequency wave energy radiation intensity (see text for definition) from the observed mean squared (MS) seismograms. The method is then applied to the Loma Prieta earthquake. 6.2 Basic Theory We begin this section with equation (2.2) in Chapter II and rewrite it as where un(x,t) is the nth component of displacement at the observation p oint, Au(^,t) is the slip discontinuity across the fault plane £, and T n(x,t; ^,0) is the traction Green's function due to an impulsive unit point force located at x acting in the xn-direction. By squaring both side of the above equation, we get We assume the form of slip time function is the same over the entire fault. The MS seismogram defined by the expectation of the squared displacement is obtained by (6. 1) 130 E[u2(x,t)] = I I E[Au (£i)Au (£2)] [v*Tn(x,t; ^i,0)*f(t)] h ,2 h i [ v T n(x,t; ^,0)*f(t)] d Z id Z 2 (6.2) where Au(^j,t)=Au(^j)f(t)Vj, 1=1, 2. We further suppose that the slip intensity Au(^) is uncorrelated between neighboring patches on the fault and assume that E[A u(^) Au(^2)]=A(^1)E[Au2(^1)]5s (^1-^2)» where is a delta function, E[Au2(^)] is the high frequency mean squared slip intensity defined at point £ on the fault. A(^), centered at point £, is the statistical patch dimension whose slip is uncorrelated with that of its neighboring patches. Equation (6.2) then becomes We shall call the square root of A(^)E[Au2(^)] the high frequency wave energy radiation intensity. From equation (2.6) in Chapter II, we find for S-waves, respectively, where tp(x, £) and ts(x, £) are the P- and S-wave arrival times, respectively, tj-(^) is the rupture time at point £ on the fault and G £(x, £) and G j(x , £) are defined in equation (2.6) of Chapter II. For the sake of simplicity, we E[ug(x,t)] A © E [ A U 2 © ] [v*Tn(x,t; ^i,0)*f(t)]2 dL (6.3) v .T n(x,t; ^0)*f(t) = 4 )-tr(£,)]vG £(x, £) for P-waves and v•T„(x,t; £,0)*f(t) = -fit- ts (x, £ ) - tr( \ )]v.G 5(x, \ ) 131 now consider S-waves only. By substituting the above solution for S-waves into equation (6.3), we get E[u2(x ,t)]= j A(S)E[AU 2(S)] { v G |( x ,^ ) * f [ t - t f ( x ,^ ) ] } 2 dZ = f 2(t )* j A(q)E[Au2 (^)l |v*G n(x, q )l28ft- tf ( x , £,)dX (6.4) where t|= ts-i-tr. Equation (6.4) is of the same form as equation (2.8) in Chapter II. Following the same procedure in deriving equation (2.9), we find E[u2( x ,t ) ] = f 2(t)* A(£jE[Au2 (!;)] ( v « G ‘ n ( x , \ ) ]2c(x, q)d L (6.5) J t l (x ,£ )= t where c(x, £)=|(I-nn) • V t|(x, ^)|_ 1 . Equation (6.5) is the isochron integration equation for our high frequency earthquake source study. c(x, ^) is the isochron velocity at point § on the fault. In next section, we will apply equation (6.5) to map the high frequency earthquake source radiation. 6.3 Application to the Loma Prieta Earthquake In this section, we will apply our isochron integration equation (6.5) to the Loma Prieta earthquake to map the source of high frequency wave radiation. We will develop this section by starting with a description of the data collection and our model, followed by the presentation of our inversion results and discussion of their physical implication. 132 6.3.1 Data Collection and M odeling In this high frequency source inversion, we use the same stations selected in Chapter V. The station distribution is shown in Figure 5.2. The instrumentation- corrected ground acceleration records have been integrated twice to obtain the displacement seismograms. These displacement data were highpass filtered at 5.0 Hz. We square these filtered displacement data and average them over a 0.6 second time interval weighted by a cosine window function, resulting in the high frequency MS seismograms defined in section 6.2. We then summed the two horizontal components for each station and interpolated the data over a 0.2 second time interval. Our velocity model is the same as used in Chapter V. A vertical profile of this model is plotted in Figure 5.3 by a solid line. This velocity structure may corresponse to a rock site model. In order to simulate a velocity model appropriate for the "average" site, we multiply the Green's function calculated based on our velocity model by the site amplification factor averaged over all stations used for the analysis relative to the Santa Cruz site. The mainshock fault is modeled in the same way as before with a 40 km by 14 km planar region striking 130°, dipping 70° and raking 130° (USGS first motion mechanism, 1990). This planar region is divided into 41 by 15 equally spaced grid points. We take a rupture velocity of 70 percent of the local shear wave velocity. The initial energy radiation intensity is assumed to be zero everywhere on the fault. We take the time function of the slip velocity in the form f(t)=[H(t-tr+At/2)-H(t-tr- At/2)]/VAt, where At=0.2 sec. The value for At is chosen to be the same with the rise time used in Chapter V. 133 6.3.2 Results and Analysis In this section, we will present three different inversion results of the high frequency MS seismograms. The first result is the simplest and does not consider the path attenuation and station site responses. Using our recursive inversion procedure, we determined the high frequency wave energy radiation intensity on the fault from the MS seismograms. The resultant distribution of the energy radiation intensity is plotted in Figure 6.1a. We shaded the energy radiation intensity with different intensities varying from black to white with dark region representing high intensity areas. Figure 6.1 b is a contour plot of the same energy radiation intensity distribution. There are three large, high frequency energy sources and two smaller ones with peak intensity values ranging from 3.6 to 6.0 (cm*km). One localized source is centered immediately above the hypocenter and is surrounded by the other localized sources. The Euclidean norm for this inversion is 8.81. We plotted the synthetic envelope seismograms (square root of the MS seismogram) in Figure 6.2 along with the actual envelope plotted above each synthetic counterpart. The number indicated along the vertical axis corresponds to the station number listed in Table 5.1. Peak values of these envelope seismograms are plotted to the right of each envelope seismogram followed by the epicentral distance of each station. In this figure, we see a large, early energy arrival for station No. 2 which does not appear in the synthetic result. We checked the three componets of the original seismograms and found this arrival to be due to the later part of the direct P-wave energy. For the remaining stations, the fit between the synthetics and the observations is quite satisfactory. In the next inversion, we incorporate the attenuation effect by using the average coda Q"1 in central California (Peng, 1989) at 6 Hz as our attenuation factor and 134 F ault D i p (km ) F au lt D i p (km) Energy Intensity Distribution Fault Strike (km) 0.000 t.t95 2.390 3.560 *.78t 3.977 Figure 6.1a. High frequency energy intensity (km-cm) distribution of the first inversion shaded in different intensity according to scales indicated below the figure. Energy Intensity Distribution __/ 10 20 15 5 to - 5 0 -10 -1 5 -20 Fault Strike (km) Figure 6.1b. Contour plot of the same energy intensity distribution as in Figure 6.1a with contour interval of 1.00 (km-cm). 135 1 . 2 . 5 . H ill IIU IIU cdm 0.10 6.94 0.15 I cdm 0.08 9.32 0.09 in cdm 0.15 16.58 |||4 I H ..... 0. 1 1 |||[l||U ill||(iiiim iw "— cdm 0.10 17.97 0.10 7 . 8 . •"M ill -•••H ill lllllllllllllilllllilll cdm 0.10 27.25 0.12 cdm 0.08 29.46 0.06 cdm 0.07 30.58 0.05 cdm 0.06 31.97 0.06 10 . J_____ I _____ I ____ I _____ I _____ I _____ l l I -3 0 3 6 9 12 15 18 21 Figure 6.2. Plot of envelope seismograms for the first inversion with the observed traces on top of the synthetic ones for each station numbered along the vertical axis. Each seismogram is normalized to its peak values (cm), indicated to the right the trace followed by its epicentral distance (km). 74086950 6237061 063880311 corrected each ray path from the fault plane to the receiver. We then performed a recursive inversion to obtain the energy radiation intensity on the fault. The inversion result of this energy radiation intensity is plotted in Figures 6.3a and 6.3b. We found the intensity pattern distribution very similar to the one obtained in the first case except for a uniform amplification by 1.6. The synthetic envelopes are compared with the observations in Figure 6.4. The Euclidean norm of this inversion is 10.84. The fitting between the synthetics and the observation is about the same as the first case. Our third inversion is done by considering both the path attenuation effect and the station site responses. The site amplification factors are obtained from the highpass filtered strong motion coda waves. We first calculate the root mean squares of the strong motion coda in a 10 second time window centered at the lapse time of 40 second for all stations. By taking the Santa Cruz station as our reference rock site, we then compute the site amplification factors for other stations as the ratio of these RMS values to that for the Santa Cruz station. These site amplification factors are listed in Table 6.1. Comparing these to the site amplification factors (Table 5.2) obtained from low frequency strong motion coda, we found the amplification effects of the sediment sites at low frequency are drastically different from the effects at high frequency. The sediment sites show relatively lower amplification as compared to the rock site. Based on the site corrected MS seismograms, we performed our recursive stochastic inversion to obtain the high frequency energy radiation intensity over the fault. The result is plotted in Figures 6.5a and 6.5b. There are considerable differences between this energy intensity distribution on the fault and our previous results. There is one large area with high energy radiation intensity of 13.26 (cm*km) centered 10 km northwest of the hypocenter and two large areas with moderate energy radiation intensity centered at 2 km and 12 km southeast of the hypercenter, respectively. The Euclidean norm of 137 Fault D ip (km) Fault D io (km) Energy Intensity Distribution - 2 0 - 1 5 - 1 0 - 5 0 5 10 15 20 Fault Stnke (km) 0.000 1.917 3.6.V S 5.752 7.670 9.588 Figure 6.3a. High frequency energy intensity (km-cm) distribution of the second inversion shaded in different intensity according to scales indicated below the figure. Energy Intensity Distribution o o 15 20 0 5 10 - 1 0 - 5 -20 - 1 5 Fault Shike (km) Figure 6.3b. Contour plot of the. same energy intensity distribution as in Figure 6.3a with contour interval of 1.60 (km-cm). a n 't | ‘M g ': J ti 138 2 . 5 . n u i i i i i i i i i i n ' i r U lU lH U IIM cdm 0.10 6.94 0.19 cdm 0.08 9.32 0.10 I M cdm 0.15 16.58 0.12 cdm 0.10 17.97 0.12 iimiimiuiiih cdm 0.05 19.02 M IU IIIM 0.16 9 . 10 . I l l IIIII< IIIIIIII< <.... m iiii ..iiiiinii.i. in cdm 0.10 27.25 0.13 cdm 0.08 29.46 0.06 cdm 0.07 30.58 0.05 cdm 0.06 31.97 0.06 1 ____I ___ I ___ 1 ___ I ____I ___ i i i - 3 0 3 6 9 12 15 18 21 VO Figure 6.4. Same as Figure 6.2 but for the second inversion results. 82 3531 748575 C$9..++...^D 858754 9631 ^818110991 6119065 5415481099 999 15699899 54 Table 6.1 Site Amplification Factors STATION NAM E Parallel To F ault S trik e P erp en d icu lar To F ault Strike 1 CORRALITOS 0.4475 0.5609 2 CAPITOL 0.9319 0.9551 3 SANTA CRUZ 1.0000 1.0000 4 WATSONVILL 0.2546 0.3412 5 LEXINTON 0.2995 0.5132 6 SARATOGA 0.5000 0.6669 7 GILROY #2 0.4371 0.6333 8 COYOTE 0.3073 0.5415 9 GILROY #4 0.4485 0.5824 140 Fault D i p (km ) Fault D i p (km) Energy Intensity Distribution . i • ! s •••* * • • * V « k.. •:*: • • . « - 2 0 - 1 5 - 1 0 - 5 0 5 10 15 20 Fault Strike (km) CO O P 2-652 5 304 7.956 IP. 60 13.26 Figure 6.5a. High frequency energy intensity (km-cm) distribution of the third inversion shaded in different intensity according to scales indicated below the figure. Energy Intensity Distribution m o 15 5 10 20 0 Fault Strike (km) Figure 6.5b. Contour plot of the same energy intensity distribution as in Figure 6.5a with contour interval of 2.21 (km-cm). 141 142 4 . 5 . H U U U I I H i w u - m m iin cdm 0.11 6.94 0.29 m ii|[ll| luillllii . cdm 0.05 9.32 0.11 cdm 0.09 16.58 0.09 cdm 0.20 17.97 0.18 U iiiM iiw iN iiu i cdm 0.07 19.02 l i ll l l U i iH U M I I H I 0.11 —1 --------- 1 --------1 ---------1 --------1 _____ I ______I _____I ____L - 3 0 3 6 9 12 15 18 21 6 . 8 . 10 . cdm .• h u h cdm i n i - J L cdm cdm J I L - 3 0 3 6 9 12 1 Figure 6.6. Same as Figure 6.2 but for the third inversion results. 0.10 27.25 0.09 0.10 29.46 0.08 0.12 30.58 0.09 0.07 31.97 0.07 J______ I ____L 5 18 21 9673^9544343 65^488 698999999 579445144 52829999999 C94../...+^D 99999999 this inversion is 10.20. The synthetic envelop seismograms are plotted in Figure 6.6 along with the site corrected observations. On careful examination of the peak envelope values between synthetics and observation for each station, we found a nearly three fold over prediction of the peak values for station No. 1 and a two-fold over prediction of the peak value for station No. 2, both are sediment site stations located within 10 km of the epicentral distance. Synthetic results for the remaining stations agree well with the observations. 6.3.3 D iscussion We applied our isochron integral equation method to the high frequency (> 5 Hz) wave energy radiation of the Loma Prieta earthquake. Inversion results show that the high frequency energy radiation intensity is highly variable over the fault. All the inversion results indicate there are three localized major energy sources on the fault characterized by high intensity values and two or three minor ones. The distribution pattern of these energy sources is very different from the fault slip distribution inferred from the same strong motion data set at low frequencies (see Figures 5.1 la and 5.1 lb). On careful examination, we found these sources are located near the boundary areas of the large slip zones. This observation is consistent with the theoretical modeling results, which show the seismic high frequencies are mostly generated when rupture stops or the rupture front that passes through fault heterogeneity area where large slip variation occurs (e.g. Madariaga, 1976, 1977, 1983; Das and Aki, 1977; Bouchon, 1978; Papageorgiou and Aki, 1983). Chin and Aki (1991) studied the high frequency radiation of the Loma Prieta earthquake using the specific barrier model consisting of equal-size subevents proposed by Papageorgiou and Aki (1983). They obtained a total of five subevents on the fault 143 plane. By taking their subevents as the source of high frequency radiation, we found their number of subevents agrees with the total number of energy radiation sources (major and minor) obtained from our deterministic mapping. The effect of attenuation was included in one of the inversion results, based on the Q factor inferred from the coda wave decay (Peng, 1989). The resultant energy radiation intensity has the same distribution as our inversion without this attenuation correction, except for the increase by 60 percent. In addition to the path effect, we also considered the site effect correction, using the coda of strong motion resords. The inversion result based on these site corrected seismograms shows a significantly different energy radiation intensity distribution from the solutions without the site correction. The synthetic result shows large over prediction for the envelope seismograms at the nearest two stations. The sediment sites show realtively de amplification as compared to the rock site. These strong motion coda site responses appear to be different at high frequency from that of the weak motion results (Phillips and Aki, 1986; Su et al., 1991). We suspect that the property of the sediment may be changed during the strong motion. Thus, we cannot acertain the effectiveness of our site-corrected inversions since the observed strong motion seismograms may subject to possible non-linear effect. Our high frequency wave energy radiation intensity is characterized by the statistical length of the patch interval and the high frequency RMS slip. We may relate this RMS slip to the RMS stress drop by the following argument. We first write the stress drop following Brune (1970) as A o = = qu/(3 on 144 where p is the shear modules, p is the shear wave velocity, u and u are the particle 3u displacement and velocity on the fault, respectively, and is the derivative of u along the fault norm. We can approximate u by u = 27tf 4r~ = TcfAu 2 Thus we obtain A a = p7tfAu/p We use the above equation to estimate the RMS stress from the RMS slip. If we take the patch interval as that of the average energy source dimensions, we get roughly VA(^) » 5 km and Aum ax * 2 cm. This leads to ^ ® m ax — [J .7 tfA u m a x / P — 3 0 b a r s where we take f = 5 Hz, p. = 3*10n bars, p = 3.2* 105 cm/sec. This is the peak RMS stress drop. The average RMS stress drop over the fault plane would be about 10 bars. Compared to the estimated average stress drop over the fault plane (e.g. 119 bars, Chin and Aki, 1991), this is a very small fluctuation but mainly responsible for the high frequency wave radiation. 6.4 C onclusion A new method was developed to map the high frequency wave energy radiation on a fault. By assuming that the slip responsible for the high frequency generation is incoherent between neighboring patches on the fault, we formulated the isochron integration equation for the high frequency wave radiation from the fault. We applied our method to the Loma Prieta earthquake. The inversion results indicate that the high 145 frequency energy source distribution shows a very different pattern from that of the low frequency fault slip distribution. The high frequency sources are located along or near the boundaries of the large slip zones. This is consistent with theoretical considerations that high frequencies are primarily generated from the rupture stopping or areas of large slip variation (e.g. Madariaga, 1976, 1977, 1983; Das and Aki, 1977; Bouchon, 1978; Papageorgiou and Aki, 1983). The number of high frequency energy sources obtained in our inversion agrees with the number of subevents obtained by Chin and Aki (1991) using the specific barrier model of Papageorgiou and Aki (1983). The corresponding average RMS stress drop is estimated at about 10 bars. We performed inversions including the effect of the path attenuation and site response. The path attenuation does not affect the pattern of the energy source distribution, but only the absolute value of the solution. The pattern, however, changed significantly when we incorporated a site effect correction in the inversion. We cannot acertain the effectiveness of our site- corrected inversions since the observed strong motion seismograms may be subject to non-linear effect. A further investigation of site effect by considering possible non linear site responses is necessary. 146 Chapter VII HIGH FREQUENCY SEISMIC WAVE SCATTERING 7.1 Introduction The study of scattered seismic waves, which is an effective tool to investigate the inhomogeneity of the earth, has drawn increasing attention in seismology in recent years. By treating the earth as a randomly inhomogeneous medium, Aki (1969), Aki and Chouet (1975), and Sato (1977) studied the coda waves as single back-scattered waves. Since then, the single back-scattering model has been widely used to interpret the observations and to calculate coda Q 1 values. The single scattering model considered by Sato (1977) has included explicitly the scattering loss along the ray path and is referred in physics literature as " first-order multiple scattering" (Ishimaru, 1978). However, because this simple model ignores all other multiple scattered waves, its validity has always been questioned. In order to clarify this problem, many theories have been proposed and a number of numerical simulations have been conducted to investigate the multiple scattering process. For instance, Kopnichev (1977) and Gao et al (1983a, 1983b) studied the 2-D and 3-D multiple isotropic scattering problem by adding higher order scattered wave contributions. Using a different approach, Wu (1985) introduced stationary energy transport theory for the multiple scattering problem. This theory can be used to separate the effects of scattering and intrinsic attenuation uniquely as demonstrated by Wu and Aki (1988) for the Hindu Kush region. Later, Frankel and Wennerberg (1987) proposed a 2-D and 3-D energy flux model for coda waves based on the finite difference simulation results obtained by Frankel and Clayton (1984, 1986) for 2-D random media. In their model, they 147 assumed that the scattered wave energy is uniformly distributed immediately after the direct wave front. Recently Hoshiba (1990) used a Monte-Carlo simulation to study the multiple isotropic scattering problem by incorporating the energy conservation law. His results provided a numerical simulation proof of the energy conservation for the multiple scattering m odel. In this chapter, we will provide a complete analysis of scattered wave energy propagation in a random isotropic scattering medium. First, we will formulate the scattered wave energy equation by extending the stationary energy transport theory studied by Wu (1985) to the time dependent case. This wave energy equation is an integral equation applicable to a medium with random isotropic scatterers. By solving the equation iteratively, we will obtain a sum of Neumann series with each term characterized by the power of the scattering coefficient which corresponds to the exact order of the scattered waves discussed by Gao et al (1983b). Further simplifying these terms, we can derive a general solution for multiple isotropic scattering coefficients in the case of separate source and receiver locations. By applying the Fourier transform to the equation, we will obtain a compact closed form integral solution for the same problem which will be very easy to evaluate numerically. Examples of this solution will then be presented and compared with the single scattering, energy-flux and diffusion models. We will also show the energy conservation law of these multiple scattered waves theoretically. Finally, we will generalize the scattered wave energy equation to a non-uniform random medium. 7.2 Basic T heory Consider acoustic waves propagating in a 3-D elastic medium with unperturbed velocity v and randomly distributed scatterers which will generate incoherent scattered 148 waves upon an incident wave arrival. The scattered wave energies in such a medium are directly additive and we can write the energy equation in this medium at time t as E(r, t) = E in (r0, r, t - Es(n , r, t - ^ - ) (7.1) ri *r The first term on the right hand side of the above equation indicates the incident wave energy at the receiver point r, and the second term is a sum of scattered wave energies from all possible scatterer points iq to the receiver r. For a medium with uniformly distributed isotropic scatterers of density n0, the scattered wave energy can be expressed as a product of the scattering cross section a, wave energy density at the scattering point, geometric spreading factor and the scattering loss and intrinsic attenuation occurring over the path from the scattering point to the receiver, that is l r - r l I r - r l P -Tilri-rl Es (ri , r, t - J-) = aE (n , t - — ---------------- (7.2) v v 47t|r,-r|2 where Tj^s+Tk, rjs is the scattering coefficient which equals n0a and is the absorption coefficient. Likewise, we can obtain the incident wave energy expression „ . _ |r-r0k c /4 . |r-r0L e_nlr_ r °' _ E m (r0, r, t ) — E m (t ) ? (7.3) v v 4jqr-ro| Substituting (7.2) and (7.3) into equation (7.1) and rewriting it in a continuous form, we obtain |r-r0| e~nlr ro 1 ( v 47c|r-ro|2 Jv E (r,t) = Eul( t - ^ ) f ^ - + | n » E ( n , t - ^ ) f ^ ‘' ' dV, (7.4a) v 47c|ri-r| ^ We will call this the "scattered wave energy equation". The above equation is the most fundamental equation from which other scattering formulas published so far applicable for a medium with randomly distributed isotropic scatterers can be derived. By integrating the wave energy equation in (7.4a) with respect to time, we can find 149 e-ri!r-rol I g-rilri-r! E(r) = Ein ^ +1 T ls E(r1) f 1 ^ d V ) (7.4b) 47i|r-r0| Jy 4n\ri -r| This is exactly the integral equation derived in the energy transport theory (Wu, 1985) and is a special case o f our general scattered wave energy equation (7.4a). Equation (7.4a) has also been given in Shang and Gao (1988), but not as the fundamental equation from which both 3-D multiple scattering formulas and 3-D energy transfer formulas can be derived. In this paper, we w ill solve equation (7.4a) using two different approaches. First, we will derive an iterative solution for equation (7.4a) and then we will obtain a compact closed form integral solution o f the same equation by applying the Fourier transform technique. 7.3 Iterative Solution of the Scattered Wave Energy Equation A usual scheme in obtaining the solution o f the integral equation is by iterative manipulation. Suppose the k-th order iterative solution is E ^ (r, t), then the (k-i-l)-th order solution E^k+1)(r, t) can be expressed as t)=Ein ( t - t p l ) e ’ 1 'r ^ + f T\s E^Cnc, t - * * * * dVk (7.5) v 4ft|r-rc| ) v v 4ft|r-r0| I- - I p.-nlr- r01 with E( 0 ) (r, t) = Ei„ ( t - —------------ v 47t|r- r0|2 If we substitute the solution E ^ (r , t) in terms o f E ^k'^(r, t) into this equation we will obtain an integral relation between E ^ k+1^(r, t) and E ^ k _ 1 ^(r, t). Repeating this process o f substituting higher order solutions in terms o f lower order ones, we can write a relation between E^k+1^(r, t) and E^°^(r, t) as E< k + 1 ) (r> t)=Eil, ( t - ^ ) e T " r ^ + 1 T ls ( • • • f v 4jc|r- r0| i-i J v ) v 150 M + 1 |rj -rj i-i| Em (t ---------) - 5 ------- , -------' --' ----- dV, • • • dVi (7.6) 47c|r-rJ2 Ft 47i|rj -rj _i|2 j = i — T id r- ^ l+X lij-ij-iD The exact solution for E(r, t) may be obtained as E(r, t) = lim E< k + 1 ) (r, t) k — v 47c|r- rc|2 + n , | e , (1- .- ' ^ r '- r°l) e~1 1 (lr ~ r|Mri~r °l) d y 1 " 4it|r- rj 4it|iy r0| + I E" (' " “ j e. (, - rwiyv rH v s ) • , ; ,J‘ „ ■*. « J v J v v 47c|r- r j 2 4n\r2 - ri|2 47c|rr r0|2 + ................. k |r- rJ+X l rj ' rj-i| e n(lr r k l+ j z1l 1 3 r j_ ,l) + t£ | •••! E in (t---------- ^ ^ dVi • • • dVk 4n\r- rj2 n4:n:|rj - rj.J2 j . . . j * 'v * ^ v " I * j = 1 + ................. k k I r_rj+ X |ri _IVi| - n ( lr - rk l+I In- r,_jl) Ei, (t ^ dVi • • • dVk (7.7) 4jt|r-rJ2 Il4ji|ri -r, ) 2 i = 1 The above expression gives the scattered wave energy as a sum o f Neumann series with each term characterized by a power o f scattering coefficient rjs. By further examining each term in this expression, we see that the zero-order term gives the incident wave energy. The first-order term gives the wave energy scattered once from 151 all possible scatterer points rj to the observation point r. This corresponds exactly to the first-order scattered wave energy. The second-order term gives the incident waves scattered at scatterers located on r1 ? then scattered again at second scatterers on r2, and finally arrives at our solution point r. These are the second-order scattered waves. As shown below, the k-th order term in equation (7.7) can be identified exactly with the k- th order scattered waves discussed in the multiple scattering wave theory o f Kopnichev (1977) and Gao e t a l (1983b). By further simplifying these terms in equation (7.7), we can derive an explicit form o f the first-order, second-order, third-order, and then the general k-th order scattered wave expressions with separate source and receiver locations for an impulsive incident wave energy. For simplicity, we assume the source located at the origin r0=0 in the following development. First-order Scattering Consider a volume element dVj between two ellipsoidal surfaces with two foci located at source and receiver, respectively as shown in Figure 7.1a. W e can write dVj as dV : = 27tr1 sin |i 1 dX where dX is the area element shown in Figure 7.1a. Putting C j as half the distance between the two foci and a: is the semimajor axis of the ellipsoidal, we can express dX in terms o f dp-j and daj as a? + c f - 2aiCicospi , , dX = ----------------------------ridpidai (ar C i cospO then 152 2ci (a) 2c i (b) ' 2c- R S (c) Figure 7.1. The Geometry for positions of the source, receiver and scatterers where S represents source, R represents receiver and numbers represent scatterers. (a) First- order scattering case, (b) Second-order scattering case, (c) Third-order scattering case. 153 , , , a? + c ? -2aiCicos|Ui _ 9 . , , dVi = -------------------- -— 2 % r? sinqidqidai (7.8) (ai-ci cospi) Substituting the above volume element into the first-order scattering term in (7.7) for an impulsive incident wave energy, we obtain (see appendix for details) r r EoS(t- i d L n l jy J0 v 47trf 4 7 t|r-rj2(ar ci cosfii)2 „ / . | | EoSft ) e 1 1 ( 1 1 + 1 r ril) (af+cF2aiCiCosp.i) Ei(r, t) = ris | I -----------— ^--------- -— ——------ —----------------- 27trf sinjiidfiidai T ls H(t-r/v)Eo e~1 1 v l l+Vvt 4 jtr v t 1-r/^ This result agrees with the result obtained by Sato (1977) for the single scattered wave with separate source and receiver locations. Second-order Scattering In addition to the volume element dV\ considered previously between two ellipsoidal surfaces, which contains the first scatterer and foci located at source and receiver (Figure 7.1b), we introduce a second volume element dV2 between two other ellipsoidal surfaces, which contains the second scatterer and foci located at the first scatterer point and the receiver (Figure 7.1b). W e find d v 2 = a? ^ ~2a2C2cospa ^ ^ _^2 sin^jd>tedaj (7,l()j (a2'C2 C O SP2) j 2 9.1 C i cosL L i where 2c2 = ------------------------ . Parameters a^ and c2 have exactly the same meaning as ai— Cicosp. 1 aj and C j but for our current ellipsoidal. Substituting these two volume elements into the second-order scattering term in equation (7.7) for an impulsive incident wave energy, we obtain (see appendix for details) 154 rjJ vEo e-^1 dq2 . j 32?t(ai + d - 2a2 c2 ^2 )v 32=^1- Defining Ai = Q = ^ , i= l,2 , and * 1= ^ , then ltavt J c J , (A ,-C i^ ,)2 1)202 D2' ° 2 r a ? — C( ,, Ai +Ci -2 A , C, 5, where Q = — , ti = ----- 1 ——— , C2 = --------------------------, D 2= l — tj. vt 2(Ai-Ci %\) 2 (A ,-C i4 i) High-Order Scattering Following the same procedures as before, we can obtain the third-order scattering term in equation (7.7) as (see appendix for details) B ( r , t) = ^ f f A - + C ? - 2 A i G ^ P f 4 7 1 J c i J-1 (A 1 C l ) 2 J c 2 J . J A2 +C2 - 2 A2C2^2 , In ]?3 +° 3 d<;2dA 2 (7.12) (A2 -C2^2> 2 D 3 C3 D 3 -C 3 u r - A ?_i — G-i ^ +C^-i-2Ai_i Ci-i ^ -1 where Ci = — , D 2= l, and t,-i = ---— — -------, C, = -------------------------------- * 2(Ai_,-Ci-i^i-,) 2 (Ai-i— Ci-i 1) Di =D i_i-ti_i, i = 2, 3. The above result can be very easily extended to an arbitrary n-th order scattering term. The result is given as n _3 t i ! E o R ,( J - ) E i (r, t) = (risvt) -------------------------------------------(7,13) and 155 B„ (—) = —L ( D' [ AL+C'- 2 A ,C l g vt r - ' V c t J , (A 1 -C 1 ^ ,f A n — 2 + C ^ - 2 _ 2 A n -2 G i - 2 ^ n - 2 , . d^n-2dAn-2 II » Cn-2 7-1 n 7 Cn-l 7-1 (An— 2 — Gi — 2 ^n — 2 )2 (An-l-G ,-l^n_l)2 Dn_C» u r , A? i — C? i r ' Ai-i +Ci-i-2Ai_i Ci_i ^_ where C i = — , D 2= l , and ti-i = ---— — -------, C i = -------------------------------- 1 2 1 2(A i_i-C i-i^ i-i) 2 (A i-i-C i-i^ i-i) Di = D i_i-ti-i, i = 2, • • • , n. Equation (7.14) gives the general expression for multiple isotropic scattering coefficients with separate source and receiver locations. A sa special example o f the source and receiver coincident case, it gives a = J L rf dA ,P f ^ J c 1 J c 2 J 1 (A2 C 2^ 2 )2 . r f 1 A .ll e l - 2 A - l a - ^ - , ^ 1 ^ lev, J , (A^.-C-,^-,)2 Q -G , Dn - c where t 1= *1, C 2= D 2 = l - tl> and h-, = ~ C^ , G = M l +C" 1 ~2Ai-‘G " ^ Di =Di_i —ti-i, i = 3, • • • , n. Using our formulation (7.15), we calculated the scattering coefficients for the first, second and third-order scattered waves and obtained essentially the same result as that from the Monte-Carlo simulation done by Hoshiba (1990) (see Table 7.1). A similar 156 T ab le 7.1 C om p arison o f S catterin g E nergy C oefficien ts Hoshiba (1990) This study B LOO 1.22 ± 0 .0 2 5 0.71 ± 0 .0 4 3 1.000 1.234 0.715 study o f the multiple scattering coefficient for the source and receiver coincident case has been conducted by Gao e t a l (1983b). However, they have made some errors in their formulation by incorrectly assigning the distance across the last two scatterers to the receiver to be larger than one half of the total traveled distance, which significantly underestimated the multiple scattering coefficients for scattering orders higher than two. 7.4 Integral Solution of the Scattered Wave Energy Equation The above iterative solution gives us a basic physical insight to the scattering process inside the random scattering medium. By examining the physical meaning of each term in the solution, we can clearly identify their scattering orders and ways how wave energies are scattered in the medium. But the formula provided by the solution are difficult to evaluate numerically. Especially with increasing scattering order, it becomes extremely tedious and time consuming. A natural alternative is to make integral transform o f the scattered energy equation and solve the equation in the integral transform domain. Referring back to the integral equation (7.4a), if we take a Laplace transform o f both sides o f the equation with respect to time, we obtain 157 e _(n+v)r 1 e -(n+^)lr-r,l E(r, s) — Ein(s) v + 1 rjs E(ri, s) — —p---- [T “d v i (7 -16) 47tr2 / 47c|r-ri|2 where E(r, s) and Ein(s) is the Laplace transform o f E(r, t) and Ej^t). The above equation has almost the same structure as the stationary energy transport equation in an uniform random medium with isotropic scatterers ( Davison, 1957; Liu and Ishimaru, 1974; Ishimaru, 1978; Wu, 1985) except for the coefficient in the - (n+ A ) r P V exponential terms. Assuming G(r, s) = ---------- , we can rewrite equation (7.16) as 4?tr2 E(r, s) = Ein (s) G(r, s) +j rjs E(ri, s)G (r - r i, s) dVi (7.17) The second term in the right side o f the equation appears in the form o f a convolution integral. Following Ishimaru (1978), w e take a spatial Fourier transform o f equation (7.17) and considering that the Fourier transform o f a convolution integral is the product o f individual Fourier transform, we find E(k, s) = Ein (s) G(k, s) +TL E(k, s)G (k, s) (7.18) Where E(k, s) is the Fourier transform o f E(r, s) and G(k, s) is the Fourier transform o f G(r, s) which equals C r|+^ v )- From equation (7.18), we can solve for E(k, s) as E^ (s) G(k, s) E(k, s) 1- rjsG(k, s) N ow w e considering an impulsive incident wave energy with Em (s) = and taking the inverse spatial Fourier transform o f the above equation, we get 158 E(r, s) = — ^ — I E(k, s)elk rdk = — \— I E(k, s) sin(kr)kdk ( 2 k ) v J 2 n v T J 0 tan ) sin(kr) dk 2 k 2 v t J o l-ris/k an -1 i tan-1(— )sin(kr) - ^ - £ [ ^ t a n - ' ( - i - ) ] ndk 2 k 2v t k rj+s/v ' 2- l tan 1 (—~ ■ ) sin(kr) dk 2 k 2 vr I rj+s/v 1— -!!T ^ - [ tan"1( ^ r r T - ) ] n + ld k (7 . 19) ; vr J(| kn ri+s/v 27t The first term in the above equation represents the Laplace transform o f the incident 5(t- r/v) e- T 1 v t wave energy which has the time domain expression o f . The second term is a summation of all individual orders o f scattered wave energies. Suppose the Laplace transform of the n-th order scattered wave energy term is En(r, s), then :„(r, s ) = - ! f - f [ tan-1 ( - ^ - ) ] n+1 dk 2ji vr J o kn r i+ i [ta n -> (-^ -)]n+1e*rd k ( 7 ,2 0 ) c vr / TI+-X kn Jo 1 v 4k , jo 1 v W e can close the contour of this integral by adding two quarter circles and a branch cut as shown in Figure 7.2. Since the integral vanishes at infinity along the two quarter circles and the integrand is analytical inside the contour, we can convert the integral in equation (7.20) to an integral along the branch cut with branch point at k=i(r)+s/v). To 159 Im k Re k Figure 7.2. The contour of integration in complex k plane and location of branch point at i(r|+s/v). the right o f the branch cut along DC, we have k=i(rj+s/v)a and tan-1 (— - — ) =tanh 1 — rj+s/v a + 5*. Along BC on the other side o f the branch cut, tan-1 (— ^t -) =tanh_ 1 -1 — 5-. Thus 2 rj+s/v a 2 we find En (r, s ) = ^ - f | (J-mnh_ 1 J-+5-)"+1- (1 tajih~1-L -g-)n+l| e ^ ' * * l(h + s/v)da 4ji2vt I i a 2 ' i a 2 [i (n+s/v)a]n .n f *(n+-)ar I ---------- !------S ( - l ) n_k(2n k+ 1 i)(taiih-1J -)n-2k(S-)2k+ld^ vr Ji Cn+s/v )n_1 k=o a 2 a 2 k vr (rj+s/v )n_1 k=o where [n/2]=n/2 for n even, (n -l)/2 for n odd, and (iS+i) is the binomial coefficient. Defining t = we can convert the above integral to r v + s t [ n / 2 ] En (r, s) = - f - 6 (1 V + S > ‘ £ ( - l ) n- k(2 n k + + 1 l) (tanh-1 i ) " - 2k(a-)2k+1^ ^ 2 n v I h (Tl+s/v)"-1 ^ ) vt 2 (vt/r) [n/2] Tin / r n-2p-(nv+s)t _ , „ , = r u _ — £ ( - l ) n“ k (2k + + 1 i) ( t a n h - 1 J _ ) n - 2k(s . )2k d i Jx (riv+s)1 1 -1 k=o 1 This integral is clearly a Laplace transform o f the time domain integral n /•tn-2 rh V V V [n/2] S ( - l ) ^ k(2 T+ 1 i)(tanh-1 J -)" -2k( f ) 2kdt1 dt2---dtn-1 (7.21) k=0 Vt' 1 A s a special case for n = l, equation (7.21) gives n c e-T lv t 1 o ! rjsH (t-—)e _ T 1 v t l + ^r El(r,t)=-^— H(t-E )il( 2)anh-1 J L = , v In ^ 47tv v rt vt 47trvt 1-T_ vt 161 which is again the first-order scattered energy term shown in equation (7.9). For n=2, we obtain the second-order scattered energy term In the above equation, we have made an integral variable transform from q to a according to the relation a ^ l/q . We can further simplify the equation by carrying out the integration for the term n 2 in the integrand and noticing that However, for scattered orders higher than two, the time domain solution o f the scattered wave energy term provided by equation (7.21) still involves a multiple integral which is difficult to evaluate numerically. One way to remedy this is to keep the time domain solution for the direct, first, and second order scattered wave energy terms, and use the original spatial Fourier transform formula to calculate the sum of all the other scattered wave energy terms. This can be accomplished through the following steps. First, we sum up the scattered wave energy terms for the orders higher than two in v we get [ s L i (in _i+ffl-)2aa ] vt r 1 -a (7.22) 162 equation (7.19). Then we change the Laplace parameter s in the result to iQ which corresponds to its Fourier transform solution with respect to time. This gives us _ & 3[tan-'( k )]4sin(kr) X E n ( r , i i 2 ) = - i — I -------------- n+lQ/V------------------dk (7.23) r*3 2 l t v r J» k ) k rj+if^v The inverse Fourier transform o f the above formula leads to the time domain sum of all multiple scattered with orders higher than two. Thus a complete solution o f our wave energy equation (7.4a) is given by 2 E (r’ ^ = — —+X En(r’ $ 4t i vr2 ^ c o (— ) [tan_ 1 (— ^— )] sin (kr) I ei£2t I ^ T l+i^v + 1 — - d n l -------------------- -— dk (7.24) J ~ 7 1 Jo 27t2vr [ 1- — tarr‘(— - — )1 [ l - f W c k )] k T |+ i^ v The integral with respect to k in equation (7.24) can be computed using the method o f discrete wave number sum (Bouchon, 1979). The fast convergence rate o f the sum is guaranteed by the inverse cube k and the alternating sine function of k in the integrand. This equation gives a compact closed form integral solution to our scattered wave energy equation (7.4a) and could be easily calculated numerically. It provides us with the temporal energy decay of multiple isotropic scattered waves in a uniform random medium for an impulsive wave energy source. Numerical examples of equation (7.24) are shown in Figure 7.3a and b as solid lines with the same absorption term ^ = 0 .0 1 km-1 but for different scattering coefficients o f r|s=0.002 and 0.05 km-1, respectively. Comparing these scattered energy decay curves, we see that multiple scattering becomes very important as the scattering coefficient increases. We also 163 i n I O c o O t * - C j j O O 1 d C O M 'o T > ^ 0 3 I O m 2 o 0 20 4 0 60 8 0 tim e (sec) Figure 7.3a. Scattered wave energy decay curves with the absorption coefficient Tl^O.Ol km-1 , the scattering coefficients rjs=0.002 km-' and source-receiver distance r=10 km. Solid curve represents our multiple scattered wave energy solution which is compared with the results of single scattering, energy-flux and diffusion models shown in dashed, dot-dashed and dot lines, respectively. i n I O C O O 0 3 O o 20 4 0 0 6 0 8 0 tim e (sec) Figure 7.3b. Same as 7.3a but for the scattering coefficients " n s=0.05 km-1. 164 compared our multiple scattering results in these figures with that of the single scattering, energy-flux and diffusion models which are shown by dashed, dot-dashed and dotted lines, respectively. For the weak scattering case, results o f single scattering and energy-flux models are in very good agreement with our multiple scattering solution. However, for the strong scattering case, they both failed in different degree to predict the multiple scattered wave energy decay, especially at later travel times. For this case, the diffusion model can be used quite well to approximate the multiple scattering result as shown by Gusev and Abubakirov (1987). Figure 7.4a and b are similar plots as Figure 7.3a and b but with fixed scattering term ris=0.01 km-1 and different absorption coefficients o f rii=0.002 and 0.05 km-1, respectively. These figures show that the slope o f the scattered energy decay curve at later travel times is very sensitive to the change in the absorption coefficient. The increase in the absorption coefficient tends to cause the energy decay to become more steep. Another interesting observation is that the energy-flux model always overestimated the scattered energy for early travel times and underestimated the scattered energy at later travel times. This is because in the early time range the incident wave energy is strong and there is more scattered energy concentrated immediately behind the incident wave front, and in later travel times the incident wave energy has attenuated and there is a scattered energy deficiency behind it, while the energy-flux model has assumed that the scattered energies are uniformly distributed behind the incident wave front at all times. We will return to more detailed discussion about the physical meaning o f our results in comparison with other models in a separate study. 165 I 10 O CO G Q o> O 0 20 40 60 80 tim e (sec) Figure 7.4a. Same as 7.3a but for the scattering coefficients T is=0.01 km-1 and the absorption coefficients H j = 0 . 0 0 2 km-1. lO I c o O C h o o sq i 0 5 I O 80 40 60 0 20 tim e (sec) Figure 7.4b. Same as 7.4a but for the absorption coefficients ^ = 0.05 kirr1 . 166 7.5 Energy Conservation N ow w e proceed to investigate the energy conservation of our formulation. For any physical system, a complete theoretical description o f the system should satisfy the fundamental law o f energy conservation. In the single back-scattering model used for the usual coda Q-1 analysis, the Bom approximation has been made in the first place. So we should not expect the energy conservation law to hold for it. However the approach used by Kopnichev (1977) and Gao e t a l (1983a, 1983b) to describe the multiple scattering wave in a random medium should satisfy the energy conservation law if all possible scattering processes are considered. Unfortunately, a rigorous discussion o f the energy conservation for this theory has not been done and some confusion has appeared in geophysical literature (e.g. Richards and Menke, 1983; Frankel and Wennerberg, 1987; Frankel, 1989). In the following section, we will prove that the energy conservation law holds exactly for our theory based on the iterative solution o f equation (7.4a). 47tvr2 point of a medium with uniform random isotropic scatterers. For simplicity we assume there is no intrinsic attenuation. The total energy at t=0 will be At time t, the wave energy is propagated and scattered out into the medium with the total energy for the direct wave arrival as Consider an impulsive incident wave energy E o ^ t-C ) y— radiated from the original 167 < b ) Figure 7.5. The Geometry for positions o f the source and scatterers where S represents the source, numbers represent scatterers and the dashed circle represents the primary wave front, (a) First-order scattering case, (b) n1 1 1 -order scattering case. 168 The total energy for the first-order scattered wave (see Figure 7.5a) can be derived by . . T lsE05 ( t - ^ ± ^ ) E ,= | | ------ v— e-i. (lb+r.JdVodVi v 47trff 47trf 1 1 r r f r ^ E o s a - ^ = I I 27tr# sin|io d|iodr01 I -----------------^— e ^ ( 5 > + r > ) 2raf sinfii d |L L i dri I f I f v47tr6 47trf Jo Jo Jo Jo r = 1 r|sE0e_ T 1 ‘vt dr0 = (r|sv t ) e ^ 'E o Jo In general, the energy o f the n-th order scattered wave will be (see Figure 7.5b) n I r t ri?E08 ( t - ^ _ ) ° En= I • • • I ----------------- — e dV0 • • • dVn v v n 47cr£ k = 0 r r r r = I I 2 k t & sin|Ltod|Uodro I I 2;tr? sinfii dfii dri Jo Jo Jo Jo ff Jo Jo Zrk rj£E05(t- *=°_) “ Y — e 27Ci£sin|jndjiT 1 drn v IT 4 k t £ k = 0 - D O - O O » O O ^ = 1 dr01 dri • • • I r i^ E o ^ t-^ ^ — ) e ~ 7 1 \? irk ^ Jo Jo Jo - v t -v t - « = 'n?E0e-'’ -v,j dro! dr, • • • I Jo Jo Jo .Vt y»V t t) f V l - i_2 drn-i= • • • 169 - v t - v t - - v t - lo----------- lk-1 4 0 4 0 40 = ---=nsE0e n! So the total energy at this moment is O O O O ( v t )n Etotal I t = y ] n~th order scattering energy = V , ~^ s V , * — e_7 1 ,v t Eo = Eo “ n! n = 0 n = 0 The total energy at t=0 is identical with that at arbitrary time t. This concludes that the energy is conserved in our multiple scattering formulation. 7.6 Scattered Wave Energy Equation in a Non-uniform Random Medium So far, w e have studied the scattered wave energy propagation in an infinite medium with uniformly distributed random isotropic scatterers. In reality, we must deal with a finite medium with non-uniformly distributed random scatterers. In such a case, we can slightly modify our equations (7.4a) and (7.4b) by incorporating the position dependent scattering and absorption coefficients, plus the path dependent scattering loss and intrinsic attenuation. 7.6.1 Basic Form ulations Suppose the scattering coefficient is T |s(r) and absorption coefficient is T|i(r), the scattering loss and the intrinsic attenuation along path / will be equal to exp { - [ [Tls(r)+Tli(r)]d/ (r)} = exp[-j^ T|(r)d/ (r)]. Thus, w e can rewrite (7.4a) and (7.4b) for the non-uniform random medium as |r-roke • ,< v 47t|r-r„|2 Jv E(r, 0 = ^ ( 1 - ^ ) — -----r ^ — + 1 tis( r , ) E ( r , , t - £ ^ ) - — ----- — dV, (7.25a) 47t|ri-r| 170 for the time dependent scattered wave energy equation and -[ Tl^)df(^) f -[ T K ^)cf(^) for the stationary scattered energy transport equation. Equations (7.25a) and (7.25b) are very difficult to solve in general if we want an analytical or even a closed form solution. However, for practical purposes, we can discretize the above equations to form linear equations for the energy terms in (7.25a) and (7.25b). Using the iterative method, we can solve these linear equations simultaneously. Examples o f the solution will be present in our next section. 7.6.2 Numerical Solutions We start this section by presenting a theoretical result o f scattered wave energy decay curves in a uniformly random scattering medium (see Figure 7.6). These scattered wave energy decay curves are obtained at different source-receiver distances as indicated by their different onset time o f the direct arrivals. We assume that the absorption and the scattering coefficients equal rj^O.Ol k m 1 and rjs=0.04 km'1 for the medium, respectively. On examination of the these curves, we found that the direct wave arrivals decay in an obviously different slop compared to the later coda wave energy decay, indicating different apparent attenuation factor Q '1 between the direct wave energy and the coda wave energy. Figure 7.7 (reproduced from Wu, 1989) is a plot o f S-wave Q 1 and S-coda Q 1 vs frequency observed from many places over the world. In this figure, large scatters are shown for the direct S-wave Q 1 measurements. , — dVi (7.25b) n-r 171 o to O C O O 03 O o O rH 150 100 200 50 time (sec) F igu re 7 .6 . Scattered w a v e energy d ecay curves ob tained at different sou rce-receiver d ista n c e s in the m ed iu m w ith ab sorp tion c o e ffic ie n t rj^O .O l krrr1 and scattering c o e ffic ien t r|s= 0 .0 4 k m '1. 172 1 0 ' * 1 0 1 0 L 9 -: M i n a Ronpp W ttrarn. a S. A (6 ) ^ CoiiforwaN. ' (81 r C a lifo rn ia (10] (3 ) k k is C onfpctur* ( 1 9 8 0 ) (7 1 (I) Kuril itian a ( 7 ) ■ S o u i iw a tw * T U SA . ) I C a s n r n a ( 7 ) 0 Cknrrol U SA . J 0 t « ^ •^ K a n to . Jap an \ . Vx (4) • ' * J °R 'l -8 ' I | ^ " J to FA »r (12 C a n tra l U .S JL (9 ) Hindu Kutn ■ 1 i i i i ■ ■ ■ 0.1 o.s 10 20 IOO (6 ) ^ Hindu Kutn (7) Friuli. Italy PAC Gorm Cain (9 ) -*v i(2 ) O IS ,11) PAC Saul Horn Coiw O f mo (9) !<attuki. Japan 1 0 0 4 0 Figure 7.7. S-wave Q '1 and S-coda Q '1 vs frequency observed from different regions over the world. 173 In average, however, the values o f S-wave Q'1 are about the same as that of the S-coda Q"1. This observation is in contradicting to results shown in Figure 7.6 for uniform scattering medium. There two possible explanations. One reason could be that the scattering attenuation coefficient is very small inside the earth compared to its absorption coefficient. But it is unlikely since a large decay train of waves following the direct S-waves is observed on local earthquake seismograms. The other reason could be due to the inadequacy o f a simple uniform scattering model for the real earth. Using the numerical method, we are going to solve equation (7.25a) in a nonuniform random scattering half space medium. W e take a 220x220x110 km3 (Length x width x depth) finite half space model with about 105 scatterers distributed in the medium. According to equation (7.25a), we formulated linear equation for each scatterers and solved these equations in 11 iterative steps which corresponds to multiple scattering up to 11-th order. We compared the exact solution for the scattered wave energy equation in a uniform medium with the approximate solutions up to certain order of scattering, all were obtained based on equation (7.24) and found that up to the 11-th scattering order, the approximate solution well converges to the exact one. In order to check the program, we first calculated numerical solutions for a uniform scattering medium with the same absorption and scattering coefficients as that of Figure 7.7. We found the solution shown in Figure 7.8 are the same as that o f the exact ones. We then tested various models with absorption and scattering coefficients change with depth. Figure 7.9 shows the results for the model o f geometry and scattering coefficient the same as our previous uniform scattering model but with a variable absorption coefficient o f 0.02 k m 1 at the free surface decreasing with depth at the rate of 0.00018 k m 1 per km. This variable absorption coefficient is chosen to provide approximately the same average value as that of the uniform scattering medium. W e found that by 174 Scattered E nergy E(r,t) t o -9 i o -8 l o -7 i o -6 i o _ O v H 20 time (sec) Figure 7.8 Numerical solutions of scattered wave energy decay curves calculated using equation (7.25a) with the same absorption and scattering coefficients as for Figure 7.7. Scattered E nergy E(r,t) o m I o c o O O co O O 40 60 20 0 time (see) Figure 7.9 Same as Figure 7.8 but for a variable absorption coefficient that decreases with depth at the rate of 0.00018 km'1 per km from a free surface value of 0.02 km'1. allowing the absorption coefficient decreases with depth we were able to simulate approximately the same decay slops for direct wave energy and coda wave energy. Thus, the observed coincidence in attenuation factor Q between direct S-waves and S- coda waves may be the consequence o f the depth dependent intrinsic and scattering attenuation in the earth. We also studied the case in which the scattering coefficient deceases with depth, but found that the result is not as sensitive to this effect as to that o f the variable absorption coefficient. 7.7 Discussion and Conclusion Two different approaches for the theoretical study o f the seismic scattering problem have been developed in the past. The single scattering theory (Aki, 1969; Aki and Chouet, 1975; Sato, 1977) and the multiple scattering theory (Kopnichev, 1977; Gao e t a l , 1983a,b) are based on the ray-theoretical approach. The energy transport theory (Wu, 1985) on the other hand, is based on the balance o f energies carried by the primary and scattered waves. In the present paper, we extended the stationary energy transport theory to the time dependent case and obtained the scattered wave energy equation applicable to a random isotropic scattering medium. The iterative solution o f this equation resulted in formulas derived from the ray-theoretical approach. We obtained a general solution o f temporal variation o f scattered energy density for arbitrary source and receiver locations as a Neumann series expansion characterized by powers o f the scattering coefficient. The first term o f this series gave the first-order scattering formula obtained by Sato (1977). For the source and receiver coincident case, our solution gave the corrected version o f high-order formulas obtained by Gao e t a l (1983b). The calculated scattering terms up to the third-order have been confirmed by the Monte-Carlo simulation result o f Hoshiba (1990). These solutions are useful 177 only for physical interpretation but not for numerical calculations. Especially for high- order scattering processes, the ray-theoretical approach becomes extremely tedious and time consuming. Solving the scattered wave energy equation using the Fourier transform technique, we obtained a compact integral solution for the temporal decay in scattered wave energy which includes all multiple scattering contributions and can be easily computed numerically. Our numerical examples of the solution showed clearly that multiple scattering becomes very important as the scattering coefficient increases. We also found that the scattered wave energy decay is very sensitive to the changes in the absorption coefficient at later travel times. The increase in the absorption coefficient tends to cause the scattered energy decay to become more steep. We further compared our multiple scattering results with that of the single scattering, energy-flux and diffusion models. For the weak scattering case, both single scattering and energy-flux models were in good agreement with the multiple scattering result. But for the strong scattering case, both models underestimated the multiple scattering effect, especially at later travel times. For this case, the diffusion model could be used to approximate the multiple scattering process quite well. Then we discussed the energy conservation for our system. By starting with our fundamental scattered wave energy equation, we demonstrated analytically that our formulas satisfy the energy conservation when the contributions from all orders of scattering are summed up. This confirmed the result obtained by Hoshiba (1990) based on a Monte-Carlo simulation. We also generalized our scattering wave energy equations to a medium with non-uniformly distributed random isotropic scatterers. Using a discretization scheme, we numerically solved the scattering wave energy equation in a vertically inhomogeneous random isotropic scattering medium and explained the observed coincidence in attenuation factor Q between direct S-waves and S-coda waves. 178 Appendix: First, Second and Third-order Scattering F irst-o rd e r S cattering Consider a volume element dV\ between two ellipsoidal surfaces with two foci located at source and receiver, respectively. We can write dV\ as dVj = 2;cr1 sin p^ dX where dX is the area element shown in Figure 7.1a. Assuming cl is half the distance between the two foci and aj is the semimajor axis of the ellipsoidal, we can express dX in terms of dpj and daj as _ a? + c f- 2aiCicosfii dX = ------------------- — 1 — rid|Liidai (ar ci cosp.i) then a?+cf-2aiC icospi . 9 . dVi = ------------------- — 27trf sinpidpidai (ar C i cospi) Substituting the above volume element into the first-order scattering term in (7.7) for an impulsive incident wave energy, we obtain . f f E0 5 ( t - — — — -) e 1 1 ( 1 1 + 1 r ril) (a?+c?-2aiCiCosp.i) Ei (r, t) = r|s I I ---------------- ^------------------------------------------ — 27trf sinpidpidai Jy2 Jo 47tr?47i|r - r j 2(ai-ci cos|Lii) . , a?+cf- 2ai C i cospi , . Since r - rA = --------------------------and 2c, =r, so ai - C i cosp.i /■" r E o g ( t - ri+l - r- - rll)e - ’ i(r,-H r-r,l) Ei (r, t) = rjs I j -------------- Y sinp,i dpi dai )rf2 Jo 8;c(a?+c?-2ai ci cos^ti)v Eo 5(t- — )e~ "2 * =ris I sinpi dpi I ------------------ dai Jy 87c(a?+c?-2aiCicosp.i)v 179 H (t-^)E e n " = T l s | --------------- ^ ----------------------- <J4l { , 4 j i [(vt) +(2c,) — 2vt(2c,)^i] H(t- i)E n s EoH(t- •C -)e-’ lvl l+ r/ f v_o = J!— J_sr!-----ln _ [ 47t[(vt)2 +r2 -2vtr^i] 4jtrvt 1 r /vt S econd-order S cattering In addition to the volume element dV\ considered previously between two ellipsoidal surfaces, which contains the first scatterer and foci located at source and receiver (Figure 7.1b), we introduce a second volume element dV2 between two other ellipsoidal surfaces, which contains the second scatterer and foci located at the first scatterer point and the receiver (Figure 7.1b). We find dV2 = — 2a2C2COSlX 2 2t c |r2 - r j 2 sinji2 dji2 da2 (a2 - c2cos|i2 )2 2aiCicosLii where 2c2 = ---------------------- . Parameters a, and c? have exactly the same meaning as a r cicosjii aj and C j but for our current ellipsoidal. , ^ , al+c|-2 a2c2cosp2 We can also find r - ri =---------------------- . a2- c 2 c o s j i2 Substituting these two volume elements and expression for Ir — r 2l into the second- order scattering term in (7.7) for an impulsive incident wave energy, we obtain E2(r, t) = f f cosp.i ( ^ 2 h Jo ( a r c , COSH,) JC 2 J0 8(t_ri+l r2 - r,[+|r- r2 |)e_, ( r i+ l rr rl+l r. r! ^ C 2 cos^ ) 47tr?47t|r2- ri|247t|r- r2|2 (a2 - c2cosj^)2 2 k |r2 - rj sinji2diii2da2 180 r f a f J# > J-l (ar Ci^i)2 L J .i 16n(a$ +c$ - 2a2 c2 ^ )v 72 ^- 1 (ai_ C l ^ ^C 2 f '2 ( a? , +c?-2aiC,^,^ [dai f Ji2 J 1 (ar C i £1 )2 J. j 327c(al +c£-2a2 c2 £2 )v rjs vEo e_ T l v t d^2 a2 = X L _lL O O O t* Defining Ai = C i= - ^ i, i=l,2, and t\ = then E2(r,t) = f d fc d A .f J c , J , ( A r C i ^ i ) 2 J., T|?Eoe-,lvtd ^ 16Tcvt(A2 2+G -2A 2C2^2) • 1 # 1 t|2E0e-’ lv ' 167lVt f f A ? + c ? - 2 A ' y . ^ l n g ± G d ^ id A i J c , J , (A,-Ci ^1 ) DjC2 ^ C2 u r' r A? — C? o A? +Ci — 2Ai Ci -p . 1 where Cj = — , ti = ————— , C2 = ----------------------- , 02=1— 1!. * 2 ( A ,- C i \ x ) 2 ( A , - C , % x ) T h ird -o rd e r Scattering Following similar steps as before, we consider a third volume element dV3 between two ellipsoidal, which contains the third scatterer and foci located at the second scatterer and receiver (Figure 7.1c), we find ... a? +c?- 2a3C3cosp.3 , ,2 . , , dV3 = ------------------- -— 2 7 1 |r3 - v\ sinp3 dji3 da3 (a3 - c3cosfa3 ) , 0 a2+c£- 2a2 c2 cosp.2 . . , , , where 2c3 = --------------------- . Again parameters a, and c. have exactly the same a2 - c2cosp2 meaning as ax and C j but for our current ellipsoidal. ... . r- j 1 1 a^+cl-2a3c3cos|J3 We can also find r - rl=-----------------------. a3- c3cos|i3 181 Substituting the three volume elements and expression for Ir — r 21 into the third- order scattering term in (7.7) for an impulsive incident wave energy, we obtain x-; , x f [ a?+c?-2ai cicos)Lii E3 (r, t)= I I ------------------ — 2jtrf sinpi djiii dai h J o (arcicospi) f f J C2 Jo al+ci -2a2 c2cosp2 , . , , 2 7 tr2s1n p .2d p .2d a2 rc2 j n (a2-C2Cospi> E0 g(t-ri+^2 " ri* + ^ " r^ r" r3l )e~T 1 Mr2 -rHrrr^r-r3 |)(a^+ d -2a3 C 3 C O sp3 ) T }1------- ,C 3 j 0 47tr? 47t |r2 -rj2 47t |r3 -rj 247t|r-r3 | ^(arCs cospa) 27t |r3 -rj 2 sinp3 dp3 da3 r (' a ? ^ -2 a ic, (' a8 ^ 2 ^ ^ h l-i (ai-ci402 Jc l-i (a2 -c ^ )2 f “ f ] Eo8(t- r 1 +2a2 +2a3 'e-rlj1 + 2 a .+ 2 a1 ) 1 1 1 -------------- *-------------------- dq3 da3 Jc, J., 167c(a?+ci-2a3 c3 ^3)v = f ^ f a ? - K i? - 2 a ,c ^ , ^ i d a i r / 2 ~ r i (' ^ h l-i (arCi^i)2 )C 2 J.j (a2 -c2 ^ ) 2 /■ ' r^Eoe-^d^ J.j 647t(a^+c?-2a3 c3 ^3 ) a^vt-ri-2a2 2 Defining Aa= Ci= i= 1,2, 3, and t! = then E 3 ( r > t ) J^ f f A?+C?-2A,C1 ^i^ dAi ^ 1c, 1-, (A ,-C,^,)2 182 where Ci = - v D i = D i_ i — ti-i p ( ^ C ? - 2 A 3 C3^d^ JCl L (A 2 -C 2 & 2 J , (A3-C3 E,i)2 tllEoe ’T'f' [' A f+ C f-2A i Ci^i d t i(jA i I i 4 * JCl h ( A r C i ^ i ) D 2 n AfcC*-2A,Cfc 1 h g f o f r j A , (A rC ^ )2 15303 D3' C3 D2 = l, and tM= A?-i -C m ; C i=A.- +Q^-2A,-, a , , i = 2,3. 183 Chapter VIII SUM M ARY AND CO N CLU SIO N The present thesis took deterministic and stochastic modelling approaches for understanding earthquake source, station site and propagation path effects on short- period seismic waves, with special emphasis on earthquake source imaging using strong motion data and modelling of seismic scattering in the lithosphere. In Chapter II, we re-derived the isochron equations for seismic waves from a fault rupture model. Some errors in the original work (Spudich and Frazer, 1984; Bernard and Madariaga, 1984) have been corrected. The expressions for the asymptotic Green's function in a homogeneous half space medium were generalized to the inhomogeneous half space case using ray theory. The solutions for the kinematic and dynamic ray tracing in a vertically inhomogeneous medium are obtained in analytic forms. We also described a detailed numerical procedure necessary for computing isochron integration and synthesizing theoretical seismograms. The formulations are applied to the synthetic seismogram calculations for simple hypothetical model to discuss the effects of nucleation and barriers on earthquake high frequency radiation using the ray theory. Then we compared our results with the published ones for complicated models. In Chapter III, a fast stochastic inverse algorithm is derived in the form of a recursive formula. The method gives a relation between the solution of the inverse based on k observations and that based on (k+1). The final step of this process leads to the same solution as that obtained by the exact matrix inversion for the whole data set, together with the error and resolution matrices of the solution. Since the recursive 184 formula takes advantage of the sparseness of the matrix in the same way as the tomographic iterative method, it leads to a fast and efficient inverse algorithm for problems involving large sparse matrix equations that commonly arise in seismological inverse problems, such as seismic tomography, earthquake fault rupture imaging by the isochron method and site amplification factor inversion based on coda power spectrum analysis. For example, in applying to the seismic tomographic problem, we found that the recursive algorithm is about Na times faster in computation compared with the direct matrix inversion, where N is the size of the matrix and a=0.5 for 2-D and 2/3 for 3-D problems. When applied to a source rupture imaging problem, our recursive inversion process gives a better solution than that obtained from the Chebyshev-semi iterative technique, and nearly 20 times reduction of computer time is achieved in comparison with the direct matrix inversion calculation. Our recursive method also leads to a simple rapid inversion procedure for the site amplification factors using coda power spectra with a drastic saving of computer time. In Chapter IV, the methods described in Chapter II and III were applied to the 1987 Whittier-Narrows earthquake. Based on the strong motion records from 16 stations within 15 km from the epicenter, we determined the fault slip and rupture time with different fault geometries for the Whittier-Narrows earthquake. In our study, we have chosen two best-fit fault geometries, with strike direction and dip angle of 270° and 35°, and 280° and 30°, respectively. To find the minimum required data size for an effective inversion we eliminated 7 stations from the above 16 stations, and found that the overall picture of the fault slip distribution remained approximately the same as the one obtained from 16 stations. This test tells us that the waveform inversion of a not so large number of records in an 185 epicentral area well-constrains the earthquake faulting process. Using the site amplification factors derived from the coda part of the strong motion records, we corrected the site effect of the strong motion seismograms. However, we did not find significant improvement or change between the inversion results with and without this site effect correction, possibly because most stations are located on similar sediment sites and the amplification factor did not show a great variation from station to station. All the inversion results indicated two regions of major slip on the fault plane for the Whittier-Narrows mainshock, one located immediately below the hypocenter and the other near the lower left (east) comer of the fault. In addition, there are three or four other smaller localized source areas around the hypocenter. Together they form a highly heterogeneous fault slip distribution for the Whittier-Narrows mainshock. In fact, the great majority of seismic moment release occurs in a 10 km by 4 km region, which is consistent with the geodetic modelling results of Lin and Stein (1989). The inversion result for the model allowing variable slip and rupture time model gave the best fit in terms of both the Euclidean norm and the weighted Euclidean norm indicating that the rupture velocity is not constant, and is faster in some directions than others. But we find no simple causal relation between these rupture process with the fault slip distributions. On the other hand, by comparing the mainshock slip pattern of our inversion results with the aftershock (Hauksson and Jones, 1989), we found that the aftershocks tend to occur over patches on the rupture surface deficient in mainshock slip. This observation can be well explained by the barrier model (Aki, 1979) of fault heterogeneity. In Chapter V, we applied our inversion method to the 1989 Loma Prieta earthquake. We performed the following three different types of fault slip inversions 186 with constant rupture velocity. In the first case, we assumed a constant rupture velocity of 80 percent of the average local S-wave velocity. In the other two cases, we modified this rupture velocity to be 70 percent of the average S-wave velocity. The inversion result showed smaller values of both the Euclidean norm and the weighted Euclidean norm for the 70 percent case. The site response for the Loma Prieta earthquake inferred from the strong ground motion coda waves shows a strong variation across the epicentral area. Comparing the rock site at Santa Cruz and the sediment site at Capitola, for example, we found about a 4 to 5-fold amplification at the sediment site relative to the rock site. Using the coda amplification factors, we corrected the strong motion displacement seismograms for the site amplification effect. The inversion result based on these site effect corrected seismograms shows much improvement in terms of the fitting between the synthetics and the observations except for the first displacement component at station Gilroy #4. With regard to the correlation between high mainshock slip and low aftershock activity, we found that the correlation was stronger for the site effect corrected data than that without the correction. In Chapter VI, a new method was developed to map the high frequency wave energy radiation on a fault. By assuming that the slip responsible for the high frequency generation is incoherent between neighboring patches on the fault, we formulated the isochron integration equation for the high frequency wave radiation from the fault. Using the recursive stochastic inversion procedure, we can determine the high frequency energy radiation intensity from the observed mean square (MS) seismograms. We applied our method to the Loma Prieta earthquake. The inversion results indicate that the high frequency energy source distribution shows a very 187 different pattern from that of the low frequency fault slip distribution. The high frequency sources are located along or near the boundaries of the large slip zones, which is consistent with theoretical considerations that high frequencies are primarily generated from the rupture stopping or areas of large slip variation (e.g. Madariaga, 1976, 1977, 1983; D asandA ki, 1977; Bouchon, 1978; Papageorgiou and Aki, 1983). The number of high frequency energy sources obtained in our inversion agrees with the number of subevents obtained by Chin and Aki (1991) using the specific barrier model of Papageorgiou and Aki (1983). The corresponding average RMS stress drop is estimated at about 10 bars. We performed inversions including the effect of the path attenuation and site response. The path attenuation does not affect the pattern of the energy source distribution, but only the absolute value of the solution. The pattern, however, changed significantly when we incorporated a site effect correction in the inversion. Our coda site responses from strong motion appear to be different at high frequency from that of the weak motion results (Phillips and Aki, 1986; Su et al., 1991). We suspect that the property of the sediment may be changed during the strong motion. In Chapter VII, we presented a complete formulation of scattered wave energy propagation from a point source in a medium with randomly distributed isotropic scatterers. First, we formulated the scattered wave energy equation by extending the stationary energy transport theory studied by Wu (1985) to the time dependent case. The iterative solution of this equation gave us a general solution of temporal variation of scattering energy density at arbitrary source and receiver locations as a power series in the scattering coefficient. The first term of this series gives the first-order scattering formula obtained by Sato (1977). For the source and receiver coincident case, our solution gave the corrected version of higher-order multiple scattering formulas 188 obtained by Gao et al (1983b). Solving the scattered wave energy equation using the Fourier transform technique, we obtained a compact integral solution for the temporal decay in scattered wave energy which includes all multiple scattering contributions and can be easily computed numerically. We compared our multiple scattering results with that of the single scattering, energy-flux and diffusion models. For weak scattering case, both single scattering and energy-flux models were in good agreement with the multiple scattering result. But for strong scattering case, both models underestimated the multiple scattering effect, especially at long lapse time. In this case, the diffusion model could be used to approximate the scattering process quite well. We then discussed the energy conservation for our system starting with our fundamental scattering wave energy equation and found that our formulas satisfy the energy conservation when the contributions from all orders of scattering are summed up. We also generalized our scattering wave energy equations to a medium with non-uniformly distributed random isotropic scatterers. Using a discretization scheme, we numerically solved the scattering wave energy equation in a vertically inhomogeneous random isotropic scattering medium and explained the observed coincidence in attenuation factor Q between direct S-waves and S-coda waves. 189 R E FE R E N C E S Achenbach, J. D. and J. A., 1978. Harris, Ray methods for elastodynamic radiation from a slip zone of arbitrary shape, J. Geophys. Res. 83, 2283-2291. Aki, K., 1968. Seismic displacement near a fault, J. Geophys. Res., 73, 5359-5376. Aki, K., 1969. Analysis of Seismic Coda of Local Earthquakes as Scattered Waves, J. Geophys. Res. 74, 615-631. Aki, K., 1979. Characterization of barriers on an earthquake fault, J. Geophys. Res., 84, 6140-6148. Aki, K., 1984. Asperities, barriers, characteristic earthquakes and strong motion prediction, J. Geophys. Res., 89, 5867-5872. Aki, K., and B. Chouet, 1975. Origin of coda waves: source, attenuation and scattering effects, J. Geophys. Res. 80, 3322-3342. Aki, K., and P. G. Richards, 1980. Quantitative Seismology: Theory and Methods, vol. 2, pp. 804-805, W. H. Freeman, San Francisco. Andrews, D. J., 1980. A stochastic fault model, 1. Static Case, J. Geophys. Res., 85, 3867-3877. Andrews, D. J., 1981. A stochastic fault model, 2. Time-dependent case, J. Geophys. Res., 86, 10821-10834. Apsel, R. J., G. A. Frazier, A. Jurkevics and J. C. Fried, 1981. Ground motion prediction for the Los Angeles basin from a major San Andreas earthquake, USGS Open-File Report 81-276, 217p. Archuleta, R., 1982. Analysis of the near-source static and dynamic measurements, from the 1979 Imperial Valley earthquake, Bull.Seism. Soc. Am., 71, 1927-1956. 190 Archuleta, R., 1984. A faulting model for the 1979 Imperial Valley, California earthquake, J. Geophys. Res., 89, 4559-4585. Bakun, W. H., G. C. P. King and R. S. Cockerham, 1986. Seismic slip, aseismic slip, and the mechanics of repeating earthquakes on the Calaveras fault, California. In Earthquake Source Mechanics. AGU Geophys. Mono. 37, ed. S. Das, C. Scholz, and J. Boatwright. Washington, D. C.: American Geophysical Union, pp. 195-207. Benioff, H., F. Press and S. W. Smith, 1961. Excitation of the free oscillations of the earth by earthquakes, J. Geophys. Res., 66, 605-618. Bent, A. L. andD. V. Helmberger, 1989. Source complexity of the 1 October 1987 Whittier Narrows Earthquake, J. Geophys. Res., 94, 9548-9556. Bernard P. and R. Madariaga, 1984. A new asymptotic method for the modeling of near-field accelerograms, Bull. Seism. Soc. Am., 74, 539-557. Beroza, G. C. and P. Spudich, 1988. Linearized inversion for fault rupture behavior: application to the 1984, Morgan Hill, California Earthquake, J. Geophys. Res., 93, 6275-6296. Beroza, G. C., 1991. Near-source modelling of the Loma Prieta earthquake: evidence for heterogeneous slip and implications for seismic harzard, Bull. Seism.. Soc. Am., 81, submitted. Blanford, R. R., 1975. A source theory for complex earthquakes, Bull.Seism. Soc. Am., 65, 1385-1405. Boore, D. M. and W. B. Joyner, 1978. The influence of rupture incoherence on seismic directivity, Bull.Seism. Soc. Am., 68, 283-300. Bouchon, M., 1978. A dynamic source model for the San Fernando earthquake, Bull.Seism. Soc. Am., 68, 1555-1576. Bouchon, M., 1979. Discrete wave number representation of elastic wave fields in three-space dimension, J. Geophys. Res. 84, 3609-3614. 191 Bouchon, M., 1982. The rupture mechanism of the Coyote lake earthquake of August 6, 1979 inferred from near field data, Bull. Seism. Soc. Am. 72, 745-759. Brune, J. N., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes; J. Geophy. Res. 75, 4997-5009. Burridge, R., and L. Knopoff, 1964. Body force equivalent for seismic dislocation, Proc. Roy. Soc. A ., 54, 1875-1888. Cerveny, V., I. A. Molotkov, I. Psencik, 1977. Ray method in seismology, Univ. Karlova, Praha. Chapman, C. H., 1978. A new method for computing synthetic seismograms, Geophys. J.,5 4, 481-518. Chin, B. H. and K. Aki, 1991. Simultaneous study of the source, path and site effects on strong ground motion during the 1989 Loma Prieta earthquake: a preliminary result on pervasive non-linear site effects, Bull. Seism.. Soc. Am., 81, submitted. Clayton, R. W. and R. P. Comer, 1983. A tomographic analysis of mantle heterogeneities from body wave travel times, Eos Trans. AGU, 64, 776. Das, S. and K. Aki, 1977. Fault planes with barriers: a versatile earthquake model, J. Geophys. Res., 82, 5648-5670. Davison, B., 1957. Neutron Transport Theory, Oxford University Press, London. Deutsch, R., 1965. Estimation Theory, Prentice Hall, Englewood Cliffs, New Jersy. Dietz, L. D. and W. L. Ellsworth, 1990. The October 17, 1989, Loma Prieta, California, earthquake and its aftershocks: geometry of the sequence from high- resolution locations, Geophys. Res. Lett., 17, 1417-1420. Ebel, J. and D. V. Helmberger, 1982. P-wave complexity and fault asperities: The Borrego mountain, California, earthquake of 1968, Bull.Seism. Soc. Am., 72, 413- 438. 192 Engdahl, E. R., S. Billington and C. Kisslinger, 1989. Teleseismically recorded seismicity before and after the May 7, 1986, Andreanof Islands, Alaska, earthquake, J. Geophys. Res., 94, 15481-15498. Frankel, A., and R. W. Clayton, 1984. A finite-difference simulation of wave propagation through random media, Bull. Seism.. Soc. Am. 74, 267-2186. Frankel, A., and R. W. Clayton, 1986. Finite-difference simulations of seismic scattering: implications for the propagation of short-period seismic waves in the crust and models of crustal heterogeneity, J. Geophys.Res. 90, 6465-6489. Frankel, A., and L. Wennerberg, 1987. Energy-flux model of seismic coda: separation of scattering and intrinsic attenuation, Bull. Seism.. Soc. Am. 77, 1223-1251. Frankel, A., 1989. A review of numerical experiments on seismic wave scattering, PAGEOPH, 131, 639-685. Franklin, J. N., 1970. Well-posed stochastic extensions of ill-posed linear problems, Journal o f Mathematical Analysis and Applications, 31: 682-716. Fukuyama, E. and K. Irikura, 1986. Rupture process of the 1983 Japan Sea, (Akita- Oki) earthquake using a waveform inversion method, Bull.Seism. Soc. Am., 76, 1623-1640. Gao, L. S., L. C. Lee, N. N. Biswas, and K. Aki, 1983a. Comparison of the effects between single and multiple scattering on coda waves for local earthquakes, Bull. Seism.. Soc. Am. 73, 377-390. Gao, L. S., N. N. Biswas, L. C. Lee, and K. Aki, 1983b. Effects of multiple scattering on coda waves in a three-dimensional medium, PAGEOPH, 121, 3-15. Gao, L. S. and S. L. Li, 1990. Time domain solution for multiple scattering and the coda envelopes, PAGEOPH, 132, 123-150. Gordon, R., R. Bender, and G. T. Herman, 1970. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography, J. Theor. Biol., 29, 4741-481. 193 Gusev, A. A. and I. R. Abubakirov, Monte-Carlo simulation of record envelope of a near earthquake. Phys. Earth Planet. Inter., 49: 30-36, 1987. Hanks, T. C., 1979. Seismological aspects of strong motion seismology proc. 2nd U. S. National Conference on Earthquake Engineering, Earthquake Engineering Institute, 898-912. Hanks, T. C., 1982. Fmax, Bull. Seism. Soc. Am. 12, 1867-1879. Hartzell, S. H. and D. V. Helmberger, 1982. Strong-motion modeling of the Imperial Valley earthquake of 1979, Bull.Seism. Soc. Am., 72, 571-596. Hartzell, S. and M. Iida, 1990. Source complexity of the 1987 Whittier Narrows, California, earthquake from the inversion of strong motion records, J. Geophys. Res., 95, 12475-12485. Hartzell, S. and T. Heaton, 1983. Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California, earthquake, Bull.Seism. Soc. Am., 73, 1553-1583. Hartzell, S. and T. Heaton, 1986. Rupture history of the 1984 Morgan Hill. California, earthquake from the inversion of strong motion records, Bull.Seism. Soc. Am., 76, 649-674. Haskell, N. A., 1969. Elastic displacements in the near-field of a propagating fault, Bull.Seism. Soc. Am., 59, 865-908. Hauksson, E. and L. M. Jones, 1989. The 1987 Whittier Narrows earthquake sequence in Los Angeles, Southern California: seismological and tectonic analysis, J. Geophys. Res., 94, 9569-9589. Heaton, T. H. and D. V. Helmberger, 1979. Generalized ray models of the San Fernando earthquake, Bull.Seism. Soc. Am., 69, 1331-1341. Hoshiba, M., 1990. Simulation of multiple scattered coda wave excitation adopting Energy conservation law, Submitted to Phys. Earth Planet. Inter.. 194 Housner, G. W. and M. D. Trifunac, 1967. Analysis of accelerograms-Parkfield earthquake, Bull.Seism. Soc. Am., 57, 1193-1220. Humphreys, E., R. W. Clayton, and B. H. Harger, 1984. A tomographic image of mantle structure beneath southern California, Geophys. Res. Lett., 11, 625-627. Ishimaru, A., 1978. Wave Propagation and Scattering in Random Media, Vols. 1 and 2, Academic Press, New York. Ivanson, S., 1983. Remark on an earlier proposed tomographic algorithm, Geophys. J. R. Astr. Soc., 75, 855-860. Kalz, S., 1988. Global optimization methods in seismic data modeling and inversion, Submitted to J. Geophys. Res.. Kanamori, H. and G. S. Stewart, 1976. Mode of the strain release along the Gibbs fracture zone, Mid-Atlantic ridge, Phys. Earth. Planet. Inter., 11, 312-332. Kawase, H. and K. Aki, 1990. Topography effect at the critical SV wave incidence: possible explanation of damage pattern by the Whittier Narrows, California earthquake of October 1, 1987, Bull.Seism. Soc. Am. 79, 1-22. Kikuchi, M., and Y. Fukao, 1985. Iterative deconvolution of complex body waves from great earthquakes-The Tokuchi-Oki earthquake of 1968, Phys. Earth. Planet. Inter., 37, 235-248. Kopnichev, Y. F., 1977. The role of multiple scattering in the formation of a seismogram's tail, Izv. Acad. Sci., USSR, Phys. Solid Earth, 13, 394-398. Kostrov, B., 1966. Self-similar problems of propagation of shear cracks, J. Appl. Math. Mech., 28, 1077-1078. Lanczos, C., 1961. Linear Differential Operators, London: Van Nostrand. Langston, C. A., 1978. The February 9, 1971 San Fernando earthquake: a study of source finiteness in teleseismic body waves, Bull.Seism. Soc. Am., 68, 1-29. 195 Lin, J. and R. Stein, 1989. Coseismic folding, earthquake recurrence, and the source mechanism at Wittier Narrows, Los Angeles basin, California, J. Geophys. Res., 94, 9614-9632. Liu, J. C. and A. Ishimaru, 1974. Multiple scattering of waves by a uniform random distribution of discrete isotropic scatterers. J. Acoust. Soc. Am., 56, 1995-1700. Madariaga, R., 1976. Dynamics of an expanding circular fault, Bull.Seism. Soc. Am., 66, 639-666. Madariaga, R., 1977. High-frequency radiation from crack (stress drop) models of earthquake faulting, Geophys. J. R. astr. Soc., 51, 625-651. Madariaga, R., 1983. High frequency radiation from dynamic earthquake fault models, Ann. Geophys., 1, 17-23. Mason, I. M., 1981. Algebraic reconstruction of a two-dimensional velocity inhomogeneity in the High Hazles seam of Thoresby colliery, Geophysics, 46, 298- 308. Matsu'ura, M. and Y. Hasegawa, 1987. A maximum likelihood approach to nonlinear inversion under constraints, Phys. Earth Planet. Inter., 47, 179-187. McMechan, G. A., 1983. Seismic tomography in boreholes, Geophys. J. R. Astron. Soc., 14, 601-612. Mendel, J. M., 1983. Lessons in Digital Estimation Theory, Prentice-Hall, Inc., Englewood Cliffs, New Jersey. Mendez, A. J. and J. G. Anderson, 1991. The temporal and spatial evolution of the 19 September 1985 Michoacan earthquake as inferred from near-source ground-motion records, Bull.Seism. Soc. Am., 81, 844-861. Mendoza, C. and S. H. Hartzell, 1988. Aftershock patterns and main shock faulting, Bull.Seism. Soc. Am., IS , 1438-1449. 196 A Menke, W., 1984. Geophysical data Analysis: Discrete Inverse Theory, Academic Press, Inc., Orlando. Olson, A. and R. Apsel, 1982. Finite faults and inverse theory with applications to the 1979 Imperial Valley earthquake, Bull.Seism. Soc. Am., 72, 1969-2001. Olson, A. H., 1982. Forward simulation and linear inversion of earthquake ground motion, Ph.D. thesis, Univ. of California, San Diego, California. Olson, A. H., 1987. A Chebyshev condition for accelerating convergence of iterative tomographic methods - solving large least squares problems, Phys. Earth Planet. Inter., 47, 333-345. Olson, A. H. and J. G. Anderson, 1988. Implications of frequency-domain inversion of earthquake ground motions for resolving the space-time dependence of slip on an extended fault, Geophysical Journal, 94, 443-455. Oppenheimer, D. H., 1990. Seismicity in the twenty years preceding the Loma Prieta, California earthquake, Geophys. Res. Lett., 17, 1199-1202. Papageorgiou, A. S. and K. Aki, 1983a. A specific barrier model for the quantitative description of inhomogeneous faulting and the prediction of strong motion, Part I: description of the model, Bull.Seism. Soc. Am., 73, 693-722. Papageorgiou, A. S. and K. Aki, 1983b. A specific barrier model for the quantitative description of inhomogeneous faulting and the prediction of strong motion, Part II: applications of the model, Bull.Seism. Soc. Am., 73, 953-978. Peng, J. Y., 1989. Spatial and temporal variation of coda Q"1 in California, Ph.D. Thesis, University o f Southern California, Los Angeles, California. Phillips, W. S. and K. Aki, 1986. Site amplification of coda waves from local earthquake in central California, Bull. Seism. Soc. Am., 76, 627-648. Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, 1986. Numerical Recipes, The art of scientific computing: Cambridge Univ. Press, Cambridge. 197 i Richards, P. G. and W. Menke, The Apparent Attenuation of a Scattering Medium, Bull. Seism.. Soc. Am. 73, 1005-1022, 1983. Rodgers, C. D., 1976. Retrieval of Atmospheric Temperature and Composition From Remote Measurements of Thermal Radiation, Rev. o f Geophys. and Space physics, 14, 609-624. Sato, H., 1977. Energy propagation including scattering effects, single isotropic scattering approximation, J. Phys. Earth, 25, 27-41. Schwartz, S. Y., J. W. Dewey and T. Lay, 1989. Influence of fault plane heterogeneity on the seismic behavior in the southern Kurile Islands Arc, J. Geophys. Res., 94, 5637-5649. Shang, T. L. and L. S. Gao, 1988. The transportation theory of multiple scattering and its application to seismic codas of impulsive source, Sinica B l, 85-94. Shakal, A. F., 1979. Analysis and modelling of the effects of the source and medium on strong motion, Ph.D. thesis, Massachusetts Institute o f Technology, Cambridge, Mass. Somerville, P. and J. Yoshimura, 1990. The influence of critical Moho reflections on strong ground motions recorded in San Francisco and Oakland during the 1989 Loma Prieta earthquake, Geophysical Research Letter, Vol. 17 No. 8, 1203-1206. Spudich, P. and E. Cranswick, 1984. Direct observation of rupture propagation during the 1979 Imperial Valley earthquake using a short baseline accelerometer array, Bull. Seism. Soc. Am., 74, 2083-2114. Spudich, P. andL. N. Frazer, 1984. Use of ray theory to calculate high-frequency radiation from earthquake sources having spatially variable rupture velocity and stress drop, Bull.Seism. Soc. Am., 74, 2061-2082. Simila, G. W., K. C. McNally, E. Nave, M. Protti-Quesada and J. Yellin, 1990. Evidence of very early aftershock activity along the northwest portion of the 18 October 1989 earthquake rupture zone, Geophys. Res. Lett., 17, 1785-1788. 198 Su, F., K. Aki, T. Teng, Y. Zeng, S. Koyanagi and K. Mayeda, 1991. The relation between site amplification factor and surficial geology in central California, Bull.Seism. Soc. Am ., submitted. Takeo, M., 1987. An inversion method to analyze the rupture processes of earthquakes using near-field seismograms, Bull.Seism. Soc. Am ., 77, 490-513. Tarantola, A., 1987. Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estimation, Elsevier, New York. Trifunac, M. D. and F. E. Udwadia, 1974. Parkfield, California, earthquake of June 27, 1966: a three-dimensional moving dislocation model, Bull.Seism. Soc. Am., 64, 511-533. Wald, D., P. G. Somerville, M. Eeri, and L. J. Burdick, 1988. The Whittier Narrows, California earthquake of October 1, 1987 - simulation of recorded accelerations, Earthquake Spectra, 4, 139-156. Wallace, T. C., D. V. Helmberger and J. E. Ebel, 1981. A broadband study of the 13 August 1978 Santa Barbara earthquake, Bull.Seism. Soc. Am., 71, 1701-1718. Wesley, J. P., 1965. Diffusion of seismic energy in the near range, J. Geophy. Res. 70, 5099-5106. Wu, R. S., 1985. Multiple scattering and energy transfer of seismic waves - separation of scattering effect from intrinsic attenuation, I. Theoretical modeling, Geophys. J. R. Astr. Soc., 82, 57-80. Wu, R. S. and K. Aki, Multiple scattering and energy transfer of seismic waves - separation of scattering effect for intrinsic attenuation, II. Application of the theory to Hindu Kush region, PAGEOPH, 128, 49-80, 1988. U. S. Geological Survey Staff, 1990. The Loma Prieta, California, earthquake: and anticipated event, Science, 247, 286-293. 199 Xu, G. M. and F. Su, 1988. Displacement field of seismic fault dislocation influenced by the surface soft layer and its inverse problem, Chinese Journal o f Geophysics, 31, 223-249. Yoshida, S., 1986. A method of waveform inversion for earthquake rupture process, J. Phys. Earth, 34, 235-255. Zeng, Y., F. Su and K. Aki, 1991. Scattered wave energy propagation in random isotropic scattering medium - Part I: Theory, J. Geophy. Res. 96, 607-619. 2 0 0

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

A surface wave study of crustal and upper mantle structures of Eurasia

PDF

Surface wave profiling for 3-D shear wave structure of the Pacific Ocean basin

PDF

The surface wave study of crustal and upper mantle structures of mainland China

PDF

Seismic stratigraphic study of the California Continental Borderland basins: Structure, stratigraphy, and sedimentation

PDF

Scattering and attenuation of seismic waves with frequencies 1 to 10 Hz in the lithosphere under southern California

PDF

Spatial and temporal variation of coda Q[-1] in California

PDF

A gravity and magnetic study of the Tehachapi Mountains, California

PDF

Stratigraphic and sedimentologic analysis of the Monterey Formation: Santa Maria and Pismo basins, California

PDF

Geometry, kinematics, and a mechanical analysis of a strip of the Lewis allochthon from Peril Peak to Bison Mountain, Glacier National Park, Montana

PDF

Late Precambrian diabase intrusions in the southern Death Valley region, California: Their petrology, geochemistry, and tectonic significance

PDF

The fabrics of deep-sea detrital muds and mudstones: A scanning electron microscope study

PDF

Seismic sequence stratigraphy and structural development of the southern outer portion of the California Continental Borderland

PDF

Structure of Santa Cruz-Catalina Ridge and adjacent areas, Southern California Continental Borderland from reflection and magnetic profiling: Implications for late Cenozoic tectonics of Southern...

PDF

Metamorphism in the Big Maria Mountains, southeastern California

PDF

The rheology of single crystal sodium chloride at high temperatures and low stresses and strains

PDF

Stable isotopes in live benthic foraminifera from the Southern California borderland

PDF

Quaternary geomorphic surfaces on the northern Perris Block, Riverside County, California: Interrelationship of soils, vegetation, climate and tectonics

PDF

Processes influencing transportation and deposition of sediment on the continental shelf, Southern California

PDF

Trends in extent and depth of bioturbation in Great Basin Precambrian-Ordovician strata, California, Nevada and Utah

PDF

Late Quaternary erosional and depositional history of Sierra del Mayor, Baja California, Mexico

Asset Metadata

Creator Zeng, Yuehua (author)

Core Title Deterministic and stochastic modelling of the high frequency seismic wave generation and propagation in the lithosphere

Degree Doctor of Philosophy

Degree Program Geological Sciences

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag geophysics,OAI-PMH Harvest

Language English

Contributor Digitized by ProQuest (provenance)

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c29-354495

Unique identifier UC11221101

Identifier DP28598.pdf (filename),usctheses-c29-354495 (legacy record id)

Legacy Identifier DP28598.pdf

Dmrecord 354495

Document Type Dissertation

Rights Zeng, Yuehua

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the au...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus, Los Angeles, California 90089, USA