University of Southern California Dissertations and Theses

A unified methodology for seismic waveform analysis and inversion

(USC Thesis Other)

A unified methodology for seismic waveform analysis and inversion

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content NOTE TO USERS This reproduction is the best copy available. ® UMI Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. A UNIFIED METHODOLOGY FOR SEISMIC WAVEFORM ANALYSIS AND INVERSION by Po Chen A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (GEOLOGICAL SCIENCES) August 2005 Copyright 2005 Po Chen Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. UMI Number: 3196790 INFORMATION TO USERS The quality of this reproduction is dependent upon the quality of the copy submitted. Broken or indistinct print, colored or poor quality illustrations and photographs, print bleed-through, substandard margins, and improper alignment can adversely affect reproduction. In the unlikely event that the author did not send a complete manuscript and there are missing pages, these will be noted. Also, if unauthorized copyright material had to be removed, a note will indicate the deletion. ® UMI UMI Microform 3196790 Copyright 2006 by ProQuest Information and Learning Company. All rights reserved. This microform edition is protected against unauthorized copying under Title 17, United States Code. ProQuest Information and Learning Company 300 North Zeeb Road P.O. Box 1346 Ann Arbor, Ml 48106-1346 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Acknowledgements First, I would like to express my deep gratitude to my advisor Professor Tom Jordan for providing me guidance in the last five years of my study at USC, for his willingness to share his ideas with me, for the time and energy he has put into developing my research, presentation and writing skills, and for all those research opportunities he has provided me in Southern California Earthquake Center. Second, I would like to thank Dr. Li Zhao for kindly providing me advices on my research and software development, and for sharing with me his research and life experiences in USA. I would also like to thank all my friends here at USC, their constant support and encouragement have made my days at USC pleasant and unforgettable. Special thanks to Cindy Waite for kindly helping me with so many registration problems, degree process and OPT paperwork. Last, I would like to thank my parents in China. I hope my efforts and achievements here in USA could make them proud. Numerical examples of Frechet kernels described in this dissertation were provided by Li Zhao (see Bibliography and References in Chapter 1, 2, and 6). Numerical computations were supported in part by University of Southern Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. California Center for High Performance Computing and Communications (http://www.usc.edu/hpcc'). All observed data seismograms, Figure 1.1, 1.4, 4.9 and 6.2 were provided by Southern California Earthquake Center. The LSQR program for tomgraphic inversion described in Chapter 6 was developed by Sebastien Chevrot at CNES, France (sebastien.chevrot@cnes.fr). Other data analysis and inversion tools were developed using Python (http://www.python.org) and SciPy (http://www.scipv.org). Geographic figures were made using the Generic Mapping Tools of Wessel & Smith (1998). Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Content Acknowledgements ii List of Tables viii List of Figures ix Abstract xii Chapter 1. Introduction 1 1.1 Predictive Forward Earthquake Simulations........................................... 1 1.2 Seismological Inverse Problems...............................................................3 1.3 Formulations...............................................................................................9 1.3.1 Finite Source Properties..........................................................9 1.3.2 Full 3D Tomography............................................................. 12 1.4 References................................................................................................. 15 Chapter 2. Generalized Seismological Data Functionals and Time-frequency-dependent Sensitivity Kernels 18 2.1 Abstract......................................................................................................18 2.2 Introduction............................................................................................... 19 2.3 An Example Application..........................................................................24 2.3.1 The Measurements Process..................................................... 25 2.3.2 The Goodness of Waveform Fitting......................................31 2.3.3 Interpretation of Frequency-dependent Phase-delay............34 2.4 Theoretical Development.........................................................................36 2.4.1 Effects of Windowing and Narrow-band Filtering on Target Waveform............................................................................... 37 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. V 2.4.2 The Summation of Gaussian Wavelets.................................. 43 2.4.3 Sensitivity Kernels for Tomography.......................................50 2.5 Discussion.................................................................................................57 2.6 References.................................................................................................59 Chapter 3. A Quantitative Evaluation of SCEC Community Velocity Model Version 3.0 Using Generalized Seismological Data Functionals 61 3.1 Abstract.....................................................................................................61 3.2 Introduction...............................................................................................62 3.3 Methodology.............................................................................................67 3.3.1 GSDF Processing..................................................................... 68 3.3.2 Sensitivity Kernel..................................................................... 74 3.4 GSDF Data Analysis................................................................................78 3.4.1 Improvement of CVM3.0 Synthetics over AID or SoCaL Synthetics...............................................................................82 3.4.2 Spatial Variations of Phase-delay Anomaly Stp.................... 86 3.4.3 Searching for Anisotropy.........................................................87 3.4.4 Attenuation Effect.................................................................... 91 3.4.5 Comparison between SCEC CVM3.0 and Harvard3D 92 3.5 Conclusions...............................................................................................95 3.6 References................................................................................................. 97 Chapter 4. Finite Moment Tensors of the 3 September 2002 Yorba Linda Earthquake 99 4.1 Abstract.................................................................................................. 99 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. vi 4.2 Introduction............................................................................................100 4.3 The Finite Moment Tensor...................................................................105 4.4 GSDF Analysis......................................................................................108 4.5 Single Event Inversion......................................................................... 112 4.6 Multiple Event Inversion......................................................................115 4.7 Results....................................................................................................118 4.8 Conclusions........................................................................................... 124 4.9 Acknowledgements.............................................................................. 125 4.10 References............................................................................................ 126 Chapter 5. Resolving Fault-plane Ambiguity for Small Earthquakes 130 5.1 Abstract..................................................................................................130 5.2 Introduction...........................................................................................131 5.3 Receiver Green Tensors and Reciprocity............................................135 5.4 CMT Inversion using 3D Synthetics.................................................. 139 5.5 Removing Fault-plane-ambiguity from CMT Solutions.................. 142 5.5.1 Finite Moment Tensors........................................................142 5.5.2 Sensitivity Derivatives for FMT Parameters..................... 144 5.6 Discussion.............................................................................................149 5.7 References.............................................................................................151 Chapter 6. Fully Three-dimensional Tomography of Los Angeles Basin Area 153 6.1 Abstract................................................................................................. 153 6.2 Introduction...........................................................................................154 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. vii 6.3 Reference Model.................................................................................... 157 6.4 Data......................................................................................................... 159 6.5 Frechet Kernels.......................................................................................161 6.6 Three-dimensional Seismic Structure of Los Angeles Basin............ 169 6.6.1 Parameterization and Inversion Procedure.......................... 169 6.6.2 Model LABF3D......................................................................175 6.7 Conclusion..............................................................................................176 6.8 References............................................................................................... 177 Chapter 7. Summary 180 7.1 Theoretical Development.......................................................................180 7.2 Model Assessment................................................................................. 182 7.3 Invert Data for Finite Source Parameters............................................. 184 7.4 Full 3D Tomography for LA Basin......................................................186 Bibliography 188 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. List of Tables 3.1 LI- and L2-norm based statistics made on GSDF measurements... 4.1 Inverted FMT parameters for 09/03/2002 Yorba Linda earthquake Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ix List of Figures 1.1 Terashake simulation at 250 seconds after initiation..................................... 1 1.2 Generalized Seismological Data Functionals (GSDF) technique.................4 1.3 Los Angeles basin and surrounding areas....................................................... 7 1.4 SCEC Community Fault Model....................................................................... 8 1.5 Four source representations commonly used in seismology....................... 10 1.6 Numerical example of Frechet kernel for P wave in uniform halfspace... 13 2.1 Philosophy of direct waveform fitting technique......................................... 20 2.2 Philosophy of GSDF technique.................................................................... 21 2.3 The source-station path from Yorba Linda earthquake of station DLA.. .26 2.4 GSDF processing for the P wave................................................................... 27 2.5 GSDF processing for the X phase.............................................................. 28 2.6 Waveform fitting by minimizing GSDF measurements.......................... 33 2.7 The sensitivity kernel for the P wave measurements................................... 34 2.8 The sensitivity kernel for the X phase measurements.................................. 35 2.9 The amplitude spectrum of C ff and Cg........................................................ 49 2.10 The time domain waveform for Cff and .................................................50 2.11 GSDF measurements...................................................................................51 2.12 The kernel C for the GSDF measurements................................................... 52 2.13 The I{t) function for the P wave isolation filter............................................56 3.1 The 72 small events and 44 Trinet stations using in this study..................63 3.2 Examples of the three velocity models..........................................................64 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. X 3.3 Transverse components of observed and synthetic waveforms.................. 65 3.4 An example of GSDF analysis....................................................................... 71 3.5 An example of using GSDF analysis.............................................................73 3.6 Frechet kernel for phase-delay time at station LAF..................................... 78 3.7(a) Distribution of GSDF measurements.............................................................83 3.7(b) Distribution of phase-delay time measurements...........................................84 3.8 Box-plots for GSDF measurements............................................................... 85 3.9 The correlation plots....................................................................................... 86 3.10 Phase-delay time as a function of basin thickness........................................87 3.11 Sensitivity kernel for the P wave recorded at DLA......................................88 3.12 Correlation between phase-delay time of SV-SH pairs................................ 89 3.13 Distribution of differential phase-delay time between SV and SH.............90 3.14 The differential phase-delay as a function of basin thickness.................... 90 3.15 Amplitude-reduction time as a function of hypocenter distance................ 91 3.16 Histograms for Havard3D measurements and other models.......................93 3.17 Histograms for phase-delay measurements with Havard3D........................94 3.18 Correlation between Havard3D and CVM3.0 measurements......................95 4.1 Four source representations discussed in the text.......................................101 4.2 Source-station configuration for the Yorba Linda event........................... 102 4.3 Some examples of recorded and synthetic seismograms........................... 103 4.4 Plots of sensitivity coefficients for Stq as functions of azimuth.................110 4.5 Examples of GSDF processing.....................................................................I l l 4.6 Measurements of S tp for P waves................................................................118 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 4.7 GSDF measurements and the fit of the inverted model.............................119 4.8 The inverted FMT parameters for the Yorba Linda event........................ 120 4.9 Projection of a 3D rendering of the Yorba Linda mainshock and its aftershocks relative to Whittier Fault..........................................................120 4.10 The solution space of the characteristic dimension and directivity 121 5.1 Los Angeles basin and surrounding areas...................................................133 5.2 The spatial source-centered grid for computing source-coordinate gradients.........................................................................................................135 5.3 Synthetics computed using reciprocity and RGT...................................... 136 5.4 CMT solutions determined using 3D synthetics........................................ 139 5.5 Plots of sensitivity coefficients for 8tq as functions of azimuth................143 5.6 Preferred fault planes for 40 small events...................................................146 6.1 Path coverage of Los Angeles basin and surrounding area...................... 154 6.2 The fence diagram of SCEC CVM3.0 S velocity model...........................156 6.3 The GSDF processing for P wave...............................................................160 6.4 The GSDF processing for X phase...............................................................161 6.5 The sensitivity kernel for the vertical component P wave........................ 164 6.6 The sensitivity kernel for the vertical component X phase....................... 165 6.7 Kernel gallery................................................................................................. 166 6.8 LABF3D, inverted 3D P velocity perturbation...........................................170 6.9 Checker board test with cell size 12 km by 12 km.................................... 171 6.10 Checkerboard test with cell size 10 km by 10 km.................................... 171 6.11 Checker board test with cell size 8 km by 8 km........................................ 172 6.12 Spike test will cell size 12 km by 12 km.................................................... 172 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. xii Abstract A central problem of seismology is the inversion of regional waveform data for models of earthquake sources and earth structure. In regions such as Southern California, preliminary 3D earth models are already available, and efficient numerical methods have been developed for solving the point-source forward problem. We describe a unified inversion procedure that utilizes these capabilities to improve 3D earth models and derive centroid moment tensor (CMT) or finite moment tensor (FMT) representations of earthquake ruptures. Our data are time- and frequency-localized measurements of the phase and amplitude anomalies relative to synthetic seismograms computed from reference seismic source and structure models. Our analysis on these phase and amplitude measurements shows that these preliminary 3D models provide substantially better fit to observed data than either laterally homogeneous or path-averaged ID structure models that are commonly used in previous seismic studies for Southern California. And we found a small but statistically significant polarization anisotropy in the upper crust that might be associated with basin layering effect. Using the same type of phase and amplitude measurements, we resolved finite source properties for about 40 earthquakes in the Los Angeles basin area. Our results on a cluster of events in the Yorba Linda area show left-lateral faulting conjugate to the nearby right-lateral Whittier fault and are consistent with the “escaping-block” hypothesis about regional tectonics around Los Angeles basin. Our analysis on 16 events in a seismicity trend that extends southwest from Fotana to Puente Hills show right- Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. lateral mechanism that is conjugate to the trend of the hypocenter distribution, suggesting a developing weak-zone that might be related to such “escaping” deformation. To set up the structural inverse problem, we computed 3D sensitivity kernels for our phase and amplitude measurements using the 3D SCEC CVM as the reference model and derived a preliminary 3D perturbation. Our results show that CVM is too slow on average and event slower in the basin area. To our knowledge, this is the first application of this “fully 3D” technique on a regional dataset. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 1 1 Introduction 1.1 Predictive Forward Earthquake Simulations Like in many other branches of predictive sciences, one important task of seismology is to predict ground motion caused by earthquakes using physics- based quantitative models. For example, what will happen if a magnitude 7.7 earthquake strikes Southern California? Recent advances both in seismology and in cyber-infrastructure have made possible large-scale, high-resolution, physics- Figure 1.1 The Terashake simulation at 250 seconds after initiation. Color shows peak ground velocity. W hite lines show major faults in this area. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 2 based predictive earthquake simulations. The largest and most detailed earthquake simulation so far, which is known as “TeraShake”, tries to answer questions such as which areas will be hit hardest if a 230 kilometer section of the San Andreas Fault ruptures from north to south, beginning near Wrightwood, California and producing a magnitude 7.7 earthquake and what are the ground velocities that can be expected to shake buildings and infrastructure in the Los Angels Basin metropolitan area (Figure 1.1). Simulations as large as TeraShake require not only large computation capability but also massive storage and fast visualization. To simulate the big earthquake, TeraShake scientists used the Anelastic Wave Model (AWM), which is a fourth-order staggered grid finite-difference code developed by K. Olson (1994), and a regular mesh with 1.8 billion grid points. The segment of the San Andreas Fault that slipped during the simulation was modeled using three right-lateral strike-slip fault planes. And the Southern California Earthquake Center (SCEC) Community Velocity Model Version 3.0 (CVM3.0), which is a 3D seismic velocity model used for regional waveform modeling in Southern California, was employed to propagate the seismic waves. Advances in cyber- infrastructure in the past two decades have made such simulations possible. TeraShake was carried out on the new 10 teraflops DataStar supercomputer and large-scale data resources of the San Diego Supercomputer Center (SDSC) at UC San Diego. In about four days, the simulation generated more than 47 TB of data that describes the details of the generation and propagation of seismic waves in a very complicated geological environment. As articulated by Tom Jordan, leader of the TeraShake collaboration and director of SCEC, TeraShake simulation not Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3 only “enriches our understanding of the basic science of earthquakes,” but also “contributes to estimating seismic risk, planning emergency preparation, and designing the next generation of earthquake-resistant structures, potentially saving lives and property.” 1.2 Seismological Inverse Problems Achieving predictive forward-modeling simulations such as TeraShake, which combines kinematic source models with dynamic models of wave propagation, requires the ability to solve two joint seismological inverse problems: from observations of ground motions infer the structure of (a) the source that generates the seismic waves and (b) the geological medium through which the waves propagate. Although these two problems are coupled, seismologists usually attempt to separate variables through a judicious choice of ground-motion observations and appropriate averaging rules. For example, the travel-times of wave groups are used to constrain focal positions and seismic velocities, while their amplitudes (and polarities) are used to constrain the focal mechanisms and attenuation parameters. These procedures work well when applied to first motion data, but they must be generalized to capture the rich information now available from the complex waveforms recorded by modern broadband instruments. In this dissertation, we develop a unified methodology for the analysis and interpretation of regional broadband waveform information in terms of seismic Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 4 source parametefs that describe the kinematic rupture process and 3D subsurface seismic velocity structures. The philosophy of this methodology is depicted in Figure 1.2. We first derive time- and frequency-dependent generalized seismological data functionals (GSDF) (Gee & Jordan, 1992; Chapter 2) that effectively capture the differences between observed and synthetic waveforms. In Chapter 2, we will discuss GSDF Seismograms M E A S U R E M E N Data Functionals 1 N V E R S w I 0 N \J Structure Model Source Model f W aveform ' V Synthesis Synthetics Figure 1.2 Generalized Seismological Data Functionals (GSDF) technique. In GSDF method, the synthetic seismograms constructed from the reference structure and source models are used to extract “generalized data functionals” from a selected wave-group at a particular frequency. And these GSDF measurements of dtV A can be used to quantitatively assess the quality of the reference models and subsequently inverted for model parameter perturbations. Any wave effects incorporated into the synthetics such as dispersion, interferences among different arrivals are all accounted for in the measurement process and the sensitivity kernels for these data functionals. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 5 processing of seismograms as well as the theoretical background behind this waveform analysis procedure. The two types of GSDF measurements, the phase-delay time (< 5 tp) and amplitude- reduction time (dtq), resemble the travel-time and the t in classic methodology, except that GSDF measurements are frequency-dependent and they are measured relative to synthetics, which incorporate our prior knowledge about seismic source and earth structure. The quality of our prior knowledge can be quantitatively evaluated using our GSDF measurements. In Chapter 3, we give a quantitative assessment of a 3D reference velocity model (CVM3.0) for Los Angeles basin using our GSDF measurements. If the prior models are not good enough according to our model evaluation, we then compute sensitivity (Frechet) kernels that linearly relate the same GSDF measurements used for model evaluation to perturbations in model parameters using our prior source and structure models. Chapter 4 and 5 give our inversion results on finite source properties. Chapter 6 discusses our regional tomography results for Los Angeles basin area. The inverted model parameter perturbations refine our prior knowledge and new synthetics can be computed using improved models and the whole process can be iterated until the differences between synthetic and observed waveforms are small enough according to our model evaluation. The methodology developed in this dissertation can be understood under the very general framework of data assimilation. By definition, data assimilation is an Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. estimation problem for the state variables, model parameters, or both. It involves combining observational data with underlying dynamical principles governing the system. Schemes for solving the estimation problem can be related to either estimation theory, usually with a Bayesian statistical approach, or control theory that solves constrained optimization problems. In estimation theory, both observations and model dynamics are treated as a priori knowledge, while in control theory, model dynamics are usually treated as strong constraints. Important estimation theory methods include Kalman Filter and Optimal Interpolation. Important control theory schemes include Generalized Inverse and Adjoint Methods. Tarantola (1986) demonstrated how to use adjoint method to solve seismic tomography problems. More recently, Tromp et al. (2004) systematically adopted adjoint theory to solve both seismic source and structure inversions. A key aspect of Bayesian analysis is the ease with which sequential analysis can be performed. As new observations become available, the entire calculation does not need to be redone; rather we use the previous posterior information as the new prior information to obtain new posterior information conditioned by the additional observations. In general, our GSDF methodology adopts a Bayesian approach. It improves our prior knowledge about seismic source and structure, which is encoded into the synthetic seismograms, by iteratively assimilating observations (GSDF measurements). The likelihood associated with the observations is computed using the dynamic wave equation that accounts for all elastic interaction effects. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 7 The posterior information of the previous inversion is then used as the prior information in the next inversion. By repeating such a data assimilation process, we expect to bootstrap our starting models, which might only be a very sketchy representation of reality, to a certain optima that can best fit the observations according to certain criteria. To demonstrate all elements of our unified methodology, we choose for our study the Los Angeles Basin area, which situates between the Transverse and Peninsular Ranges of Southern California and right beneath metropolitan Los Angeles (Figure 1.3). Formed in early Neogene as a product of regional crustal extension ■4000 240” 2 4 1 ” 242° 243° Figure 1.3 Los Angeles basin and surrounding areas. Blue triangles show the locations of California Integrated Seismic Network (CISN) three-component broadband stations. Black solid lines show the surface traces o f active faults in this area. Background color shows the thickness of sedimentary basin. The inner box is the area where we have conducted seismic source and structure inversions. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 8 associated with the opening of California Borderland and rotation of the Transverse Ranges (Luyendyk and Homafius, 1987; Wright, 1991; Crouch and Suppe, 1993; McCulloh et al., 2000), this area has undergone a complex tectonic history including extension and strike-slip faulting in the Oligocene/Miocene and oblique contraction, through thrusting and strike-slip faulting, in Pliocene and Quartemary (Davis et al., 1989; Hauksson, 1990; Wright, 1991; Shaw and Shearer, 1999), resulting in a very complex 3D seismic velocity structure underground. It is still tectonically very active today, deforming under the interactions among complex 3D fault systems. Figure 1.4 shows the Community Fault Model (CFM) of SCEC. Armed with modem broadband instruments, this area has been under intense seismological studies for decades, providing us Figure 1.4 SCEC Community Fault Model. This map shows the 3-dimensional structure of major faults beneath Southern California. Vertical faults such as the San Andreas (yellow band from top left to bottom right) are shown as a thin strip. Faults that are at an angle to the surface are shown as wider ribbons of color. The nearest fault to you might be a few miles beneath your home. Areas that seem to have few faults can still experience strong shaking from earthquakes on unmapped faults or from large earthquakes on distant faults. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 9 abundant information on basin structure (Magistrate et al, 2000; Suss and Shaw, 2003) and well-documented earthquake catalogs. 1.3 Formulations Generally speaking, the source and structure effects are coupled in our GSDF measurements. As a result, in our unified methodology, we formulate the linearized forward problem in terms of both source and velocity structure variations, + fK rf(r,m 0)-5 m (r)JF (r). (1.1) v Here, d is our GSDF measurement, p is comprised of seismic source parameters we are inverting for, K,/ is the Frechet kernel for measurement d with respect to velocity model parameter m, r is the spatial coordinate, mo is the prior velocity structure model, 5m is the perturbation to model parameter m and the integration is over V, the volume of the earth under study. Our task in this dissertation is to determine both the seismic source parameter p and seismic velocity perturbation 5m from the GSDF measurements d in a mutually consistent manner. 1.3.1 Finite Source Properties In previous seismic studies, three types of source representations are commonly used (Figure 1.5): isotropic point source (IPS), specified by an origin time to, a Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 10 hypocenter ro, and local magnitude Ml; Centroid moment tensor (CMT), specified by a symmetric seismic moment tensor M, a centroid time t\, and a centroid location ri (Dziewonski et al., 1981); fault slip distribution (FSD), specified by a displacement discontinuity Au(r,t) across a predetermined fault plane. Isotropic Point Centroid Moment Tensor (CMT) Finite Moment Tensor (FMT) Fault Slip Distribution (FSD) 1 --• --- 1 + < Number of -» (5) parameters ^ V (8-10) 1 (14-20) 1 1 © — O s 1 2 3 4 5 M agnitude 6 7 Figure 1.5 Four source representations discussed in the text, arrayed near the lower end of the magnitude range of their applicability in Southern California. We have developed the methodology for recovering finite moment tensor (FMT) representations of earthquakes in the magnitude range 4-6 from regional broadband data. The number of FM T parameters varies from 14 to 20, depending on the constraints on the source model. In Southern California, IPS parameters are routinely determined for earthquakes as small as M l ~ 1.0-1.5 (Mori et al., 1998). CMT parameters are regularly determined from first motions down to M l ~ 2.0-2.5 (Hauksson, 2000). Complete CMT solutions can be obtained from regional broadband waveform data for earthquakes with local magnitude larger than 3.0-3.5 (Zhu and Helmberger, 1996; Pasyanos et al., 1996; Liu et al., 2004). Due to the large number of parameters, Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 11 complete FSD imaging is usually applied to large earthquakes (M l > 6) with abundant high-quality observations, such as 1992 Landers earthquake (Wald and Heaton, 1994; Cohee and Beroza, 1994), 1994 Northridge earthquake (Wald et al., 1996), and 1999 Hector Mine earthquake (Dreger and Kaverina, 2000). In this dissertation, we discuss new procedures for determining source parameters using a “finite moment tensor” (FMT) — the lowest order representation of a finite fault. The FMT resolves fault plane ambiguity of the CMT and provides several additional parameters of seismological interest, including the characteristic length Lc, width Wc, and duration Tc of the faulting, as well as the directivity vector V d of the fault slip. We expect the FMT representation to be most useful for earthquakes in the magnitude range M l ~ 4-6, where finite-source effects can be resolved from regional broadband data but are too small to warrant a full FSD analysis (Figure 1.5). We have successfully applied our method to resolve fault-plane-ambiguity for 40 small to medium-sized earthquakes (3 < Ml < 5) using low-frequency waveform data and we have recovered the full FMT parameters for 3 September 2002 Yorba Linda earthquake (M l 4.8) using broadband high-quality waveform data from the California Integrated Seismic Network (CISN). The details about the inversion will be discussed in Chapter 4 and 5. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 12 1.3.2 Full 3D Tomography The formulation of the tomography problem can be characterized by the dimensionalities of K* mo and 8m. A common method for determining subsurface seismic velocity structure is by analyzing travel-time residuals, which is the difference between observed and theoretical travel-times. Even though the inverted model perturbation 5m could be 3D, this type of technique usually assumes a ID reference velocity model and calculates theoretical travel-times using ray theory, which is the high-frequency approximation to wave equation. This approach has been applied both to regional (Aki et al, 1977) and to global (Dziewonski et al., 1977) tomography studies. Because the Fermat approximation implied in this type of technique presumes no ray-path change due to model perturbation, the sensitivity kernel is nonzero only on the geometric ray-path predetermined by the reference model and source-receiver configuration. Advances in computing technology in the past two decades have made 3D ray- tracing much more affordable, Bijwaard and Spakman (2000) have extended the traditional linearized approach to a nonlinear fashion by taking into account seismic ray bending due to inferred 3D velocity heterogeneities and iteratively updating ray paths and travel-times during the inversion. If the scale of the heterogeneities is comparable with the wavelength, tomography techniques based on ray theory and Fermat principle might not be adequate. In such a case, we need to consider finite-frequency wave propagation and associated constructive and destructive interference effects in our inversion. By Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 13 combining ray theory with the single-scattering Bom approximation (Marquering et al, 1999), Dahlen et al. (2000) derived 3D sensitivity kernels for finite- frequency body-wave travel-times. The resulting 3D kernels have the appearances of hollow bananas: appearing as bananas if visualized in the planes of wave X-component pS pP H r - 5 10 15 time (sec) free surface 5tp 8tp H l^ S jiiik ■ T il vJ free surface Figure 1.6 Numerical results for the P wave from an explosive source in a uniform half space. Plotted are the X-component finite-difference synthetic and the Frechet kernels of the P wave phase-delay tim e Stp and amplitude-reduction tim e 5tq. The source and receiver, denoted as S and R, respectively, have a distance of 32.2 km and are both at the same depth of 24 km. X is the direction from the source to receiver and Z is upward. The banana-doughnut phenomena can be seen in the two plots for Stp in the source-receiver plane vertical plane and the transverse plane midway between the source and the receiver indicated by the dash line. W arm colors represent negative sensitivity indicating a velocity increase leads to an advance in arrival time and an increase in am plitude and cool colors represent positive sensitivity indicating a velocity increase leads to a delay in arrival time and a reduction in amplitude. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 14 propagation, but as doughnuts on cross-sections perpendicular to ray-paths (Figure 1.6). Somehow counter-intuitively, this suggests body-wave travel-times have zero sensitivity on the ray-paths. This type of 3D banana-doughnut kernels were generalized by Zhao et al. (2000) to arbitrary phases using normal-mode theory in a ID reference model. So far, tomography studies for crustal structures in Los Angeles Basin area have been largely based on inverting P-wave onset time using ray theory and Fermat approximation (Hauksson, 2000). The rich information encoded in waveforms is not well utilized. This problem is partly due to the complexity of regional waveform data. At regional distances, near-surface 3D heterogeneity plays a more important role on wave propagation, which makes it difficult to separate direct, scattered, converted and surface wave phases. To utilize such waveforms, we need to consider all elastic interactions that generate these phases. So we adopted Kim Olsen’s 3D finite-difference method as our wave propagation method and calculated 3D sensitivity kernels for our GSDF measurements using 3D starting models. As a result, our “full 3D” methodology can utilize waveform information for an arbitrary phase on a seismogram since our sensitivity kernel accounts for all elastic interaction effects that generate this phase. Using this “full 3D” technique, we have obtained 3D perturbations to the SCEC CVM3.0, which we use as our 3D reference model for Los Angeles Basin. In Chapter 6, we will discuss the details about our tomography experiment. To our knowledge, this is the first application of the “full 3D” seismic tomography to a regional dataset. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 15 1.4 References Aki, K., Christoffersson, A. and Husebye, E.S., 1977. Determination of the three- dimensional seismic structure of the lithosphere, J. Geophys. Res., 82, 277-296. Backus, G., and M. Mulcahy, Moment tensors and other phenomenological descriptions of seismic sources, I, Continuous displacements. Geophys. J. Roy. Astron. Soc. 46, 341-361, 1976. Bijwaard H. and W. Spakman, 2000. Nonlinear global P-wave tomography by iterated linearized inversion, Geophys. J. Int., 141, 71-82. Cohee, B., and G. Beroza, Slip distribution of the 1992 Landers earthquake and its implications for earthquake source mechanics. Bull. Seismol. Soc. Am. 84, 692- 712, 1994. Crouch, J. K., and J. Suppe, 1993, Late Cenozoic tectonic evolution of the Los Angeles basin and inner California borderland: A model for core complex-like crustal extension: Geological Society of America Bulletin, v. 105, p. 1415-1434. Dahlen, F. A., S.-H. Hung and G. Nolet, 2000. Frechet kernels for finite- frequency traveltimes - I. Theory. Geophys. J. Int. 141, 157-174. Davis, T. LJ. Namson and R. F. Yerkes, 1989. A cross-section of the Los Angeles area: seismically active fold-and-thrust belt, the 1987 Whittier Narrows earthquake and earthquake hazard. J. Geophys. Res., 94, 9644-9664. Dreger, D., and A. Kaverina, Seismic remote sensing for the earthquake source process and near-source strong shaking: A case study of the October 16, 1999 Hector Mine earthquake. Geophys. Res. Lett. 27, 1941-1944, 2000. Dziewonski, A. M., B. H. Hager, and R.J. O'Connell, 1977. Large-scale heterogeneities in the lower mantle, J. Geophys. Res., 82 (B2), 239-255. Dziewonski, A. M., T. A. Chou, and J. H. Woodhouse, Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. J. Geophys. Res. 86, 2825-2852, 1981. Gee, L. S. and T. H. Jordan, 1992. Generalized Seismological data functionals, Geophys. J. Int., Ill, 363-390. Hauksson, E., 1990. Earthquakes, faulting and stress in the Los Angeles basin. J. Geophys. Res., 95, 15365-15394. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 16 Hauksson, E., 2000. Crustal structure and seismicity distribution adjacent to the Pacific and North America plate boundary in southern California. J. Geophys. Res. 105, 13875-13903. Liu, Q.-Y., J. Polet, D. Komatitcsh and J. Tromp, Spectral-element moment tensor inversions for earthquakes in Southern California. Bull. Seismol. Soc. Am., 94, 1748-1761, 2004. Luyendyk, B. P. & J. S. Homafius, 1987. Neogene crustal rotations, fault slip and basin development in southern California, in Cenozoic basin development of coastal California, edited by R. V. Ingersoll & W. G. Ernst, pp 259-283, Prentice- Hall, Old Tappan, New Jersey. Magistrale, H., S. Day, R. W. Clayton and R. Graves, 2000. The SCEC Southern California Reference Three-dimensional Seismic Velocity Model Version 2. Bull. Seismol. Soc. Am. 90, S65-S76. Marquering, H., F. A. Dahlen and G. Nolet, 1999. Three-dimensional sensitivity kernels for finite-frequency traveltimes: the banana-doughnut paradox, Geophys. J. Int., 137, 805-815. McGuire, J. J., L. Zhao, and T. H. Jordan, Teleseismic inversion for the second- degree moments of earthquakes. Geophys. J. Int. 145, 661-678, 2001. Mori, J., H. Kanamori, J. Davis, E. Hauksson, R. Clayton, T. Heaton, L. Jones, A. Shakal, R. Porcella, Major improvements in progress for southern California earthquake monitoring, Eos 79, 217-221, 1998. Olsen, K. B., 1994. Simulation of three-dimensional wave propagation in the Salt Lake Basin, Ph.D. Thesis, University of Utah, Salt Lake City, Utah, 157p. Pasyanos, M. E., D. S. Dreger, and B. Romanowicz, Towards real-time determination of regional moment tensors. Bull. Seismol. Soc. Am. 86, 1255-1269, 1996. Shaw, J. H. and P. Shearer, 1999. An elusive blind-thrust fault beneath metropolitan Los Angeles, Science, 283, 1516-1518. Suss, M. P. and J. H. Shaw, 2003, P-wave seismic velocity structure derived from sonic logs and industry reflection data in the Los Angeles basin, California, J. Geophys. Res., 108, NO. B3, 2170, ESE 13. Tarantola, A., A strategy for nonlinear elastic inversion of seismic reflection data, Geophysics (1986), 51, 1893-1903. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 17 Tromp, J., C. Tape and Q. Liu, Seismic tomography, adjoint methods, time reversal, and banana-doughnut kernels. Geophys. J. Int. (2005) 160, 195-216. Wald, D. J., and T. H. Heaton, Spatial and temporal distribution of slip for the 1992 Landers, California earthquake. Bull. Seismol. Soc. Am. 84, 668-691, 1994. Wright, T. L., 1991. Structural geology and tectonic evolution of the Los Angeles basin, California. In Active Margin Basins, edite by K. T. Biddle, vol. 52, pp. 35- 134. Am. Assoc. Pet. Geol. Memoir. Zhao, L., T. H. Jordan and C. H. Chapman, 2000. Three-dimensional Frechet differential kernels for seismic delay times. Geophys. J. Int., 141, 558-576. Zhao, L., T. H. Jordan, K. B. Olsen and P. Chen, 2005. Frechet Kernels for Imaging regional earth structure based on three-dimensional reference models. Bull. Seismol. Soc. Am, submitted. Zhu, L. and D. Helmberger, Advancement in source estimation techniques using broadband regional seismograms, Bull. Seismol. Soc. Am., 86, 1634-1641, 1996. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 2 18 Generalized Seismological Data Functionals and Time- frequency-dependent Sensitivity Kernels 2.1 Abstract We formulate a waveform-analysis procedure to recover time-frequency- dependent phase and amplitude information from individual seismograms that makes use of the ability to compute complete synthetic seismograms from realistic seismic source and structure models. The basic tool is the “isolation- filter”, a synthetic waveform constructed to select data from a desirable portion of the seismogram. When the cross-correlation between the isolation-filter and the observed seismogram is localized in time domain through windowing and in frequency domain through narrowband filtering, the resulting cross-correlagram can be well approximated using a Gaussian wavelet parameterized by the apparent differential phase-delay time dtp and amplitude-reduction time Stq. We have developed procedures to measure these two quantities by fitting a Gaussian wavelet to the windowed and narrowband filtered cross-correlagram. These two time-frequency-dependent quantities are functionals of seismic source and structure models and we have derived linearized equations to relate these data functionals to perturbations in model parameters. The effects of windowing and Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 19 narrowband filtering on the target waveform are completely accounted for in the corresponding sensitivity kernels for these data functionals. Compared with the “banana-doughnut” kernel for a cross-correlation travel-time (Dahlen et al. 2000) and the 3D Frechet kernel for a finite frequency amplitude measurement (Dahlen & Baig 2002), the sensitivity kernels derived for Stp and Stq depend not only on the time window selected but also on the frequency-domain localization at which we make our StP A measurements. Using this type of kernel, we can better utilize the waveform information than the traditional cross-correlation technique that derives a single measurement over a relatively broad frequency band. 2.2 Introduction Earthquake hazard analysis depends on seismology to quantify both the seismic source and the subsurface velocity structure. A major focus of seismological research today concerns how to construct earthquake source and earth structure models from the rich information carried by the complex waveforms recorded by modem broadband instruments. Techniques that deal with waveform inversion problems can be loosely separated into two categories: (a) direct waveform fitting techniques (Figure 2.1) that map waveform mismatches directly to model parameter perturbations (Tarantola, 1986; Nolet, 1990; Ak^elik et al., 2003; Tromp et al., 2005) and (b) data functional techniques (Figure 2.2) that isolate the waveform-fitting step, which is nonlinear, from the inversion step, which is linear or can be well linearized, by estimating data functionals from seismograms that on one hand can quantify waveform misfit and on the other hand can be used as Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 20 linear constraints on model parameters (Sipkin & Jordan, 1976; Gee & Jordan 1992; Gaherty et al., 1996; Katzman et al., 1998). Structure Model Seismograms Waveform Inversion i i Source Model $ Waveform Synthesis Synthetics Figure 2.1. Direct waveform fitting technique. Advances in computing and numerical techniques in the past two decades have made accurate waveform synthesis affordable. Synthetic seismograms can be constructed from prior knowledge about earth structure and seismic source. In direct waveform fitting algorithms, the differences between synthetic and observed seismograms are inverted directly for model perturbations. The whole process can be iterated until the misfits between synthetic and observed waveform are small enough according to a certain criterion or no further improvement to the seismic models can be gained from the data. Due to the harmonic nature of a band-limited signal, the manifold of the models consistent with the direct waveform data may have many local minima. Those direct waveform-fitting processes are highly nonlinear with respect to model parameters (Luo & Schuster, 1991). To avoid converging to local minima, direct waveform fitting techniques usually require strong low-pass filtering to the whole Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 21 seismogram in early iterations. One disadvantage of such low-pass filtering process is that it sacrifices resolution in time domain and may mix signals with noise that are beyond the explanation ability of model parameterization. Because of the existence of such low-frequency noise, models consistent with one seismogram might not be compatible with models consistent with other seismograms and a simple averaging among incompatible models might render the final inversion results meaningless. Seismograms A M E A S U R E M E N W Data Functionals S u Sta \r A 1 N V E R S w I o N Structure Model Source Model W aveform Synthesis Synthetics Figure 2.2. Generalized Seismological Data Functionals (GSDF) technique. In GSDF method, the synthetic seismograms constructed from the reference structure and source models are used to extract “generalized data functionals” from a selected wave-group at a particular frequency. And these GSDF measurements of < )fpq can be used to quantitatively assess the quality of the reference models and subsequently inverted for model parameter perturbations. Any wave effects incorporated into the synthetics such as dispersion, interferences among different arrivals are all accounted for in the measurement process and the sensitivity kernels for these data functionals. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 22 For a transient source such as an earthquake, energy is unevenly distributed on the seismogram. To enhance signal-to-noise ratio, a robust waveform inversion technique should be able to exploit resolution both in time domain and in frequency domain. The operations that we usually adopt to isolate signals from surrounding noise are correlation, localization in time domain by windowing and localization in frequency domain by narrowband filtering. These operations can isolate the desirable signals from noise but can also modify the signals themselves, so we need to account for their effects on the target wave-groups in our analysis. The Generalized Seismological Data Functionals (GSDF) technique (Figure 2.2) developed by Gee & Jordan (1992) systematically takes into account any effects of windowing, filtering and interferences from other arrivals on the seismogram. It derives from seismograms frequency-dependent phase-delay time dtp, which is a generalized differential travel-time, and amplitude-reduction time Stq, which is a generalized differential t*. Just like traditional amplitude information or travel-time measured from cross-correlagrams, these two types of generalized data functionals are much more linear with respect to model parameters than direct waveform data. Whereas traditional cross-correlation techniques consider only a single travel-time delay defined by the cross- correlagram peak (Dahlen et al. 2000), the GSDF method extracts from each cross-correlagram a more complete set of information that encodes all the waveform information. In other words, the GSDF method reorganizes information Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 23 unevenly distributed on the whole seismogram into a more condensed data set that can be more directly interpreted in terms of model parameterization. The original GSDF analysis (Gee & Jordan, 1992) was developed for signals with relatively narrow bandwidth, such as low-pass filtered GDSN recordings. If the frequency band is too wide, at the boundaries of the band, the corrections for the effects of windowing and narrowband filtering on the target wave-group become unstable. To account for this disadvantage, we used to prefilter a broadband signal using a series of band-pass filters with relatively narrow bandwidth. On each prefiltered signal, GSDF analysis can be successfully applied. The disadvantage of such a prefiltering process is that it sacrifices the resolution in time domain, spreads out the originally concentrated energy and mixes signals with noise even further. In this study, we show that by substituting the Gram-Charlier expansion [Equation (15), Gee & Jordan, 1992] with a discrete Fourier transform and taking advantage of certain properties of Gaussian wavelets, GSDF analysis can easily be adapted to broadband signals without prefiltering. Using this new formulation of GSDF, we also derived sensitivity kernels for waveform tomography and seismic source parameter inversion. Numerical experiments conducted to test this new theory will also be presented. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 24 2.3 An Example Application GSDF analysis has been successfully applied to the analysis and inversions of teleseismic data. By inverting about 1500 frequency-dependent phase-delay times measured on 233 waveforms including both body waves and surface waves, Gaherty et al. (1996) derived a ID structure model (PA5) for Tonga-Hawaii corridor. Using PA5 as the reference model, Katzman et al. (1998) obtained a 2D high-resolution tomogram by inverting 1122 frequency-dependent phase-delay times of three-component turning and surface waves as well as travel-times of ScS reverberations. The resulting image revealed fine-scale patterns of shear-speed highs and lows within upper mantle, suggesting various kinds of sublithospheric convection beneath a sufficiently long-lived plate. GSDF methodology was also adopted to study seismic source properties. Ihmle and Jordan (1994) used GSDF method to estimate source-time functions for “slow earthquakes”. McGuire et al. (2001) inverted both amplitude-reduction times and phase-delay times measured on teleseismic data to resolve finite source properties of large earthquakes. The first full-fledged application of GSDF method to regional and local seismic inverse problems was conducted in Chen et al. (2005). In that paper and Chapter 4-5 of this dissertation, GSDF measurements made on local broadband waveforms were used to estimate finite source properties of small to medium sized earthquakes. In Chapter 6, we use GSDF method to obtain a 3D tomogram of Los Angeles Basin area. To our knowledge, this is the first full 3D tomography ever done using a local broadband dataset. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 25 The power of GSDF technique lies in its generality: the same procedure can be applied to any phases the synthetic can model, and in its precision: it automatically accounts for the frequency-dependent nature of waves and the various interferences among different wave constituents. In the following, we illustrate various aspects of this methodology with an example. 2.3.1. The Measurement Process The seismogram is from the Ml4.8 (Mw4.3) Yorba Linda earthquake happened on 3 September 2002 in the Los Angeles Basin area recorded at the California Integrated Seismic Network (CISN) station DLA (Figure 2.3). The depth of this earthquake is about 8.5 km and the epicenter distance is about 30 km. Figure 2.4(a) plots the vertical component BHZ of the observed data s(t). The synthetic seismogram [Figure 2.4(a), s(t) ] was computed using the fourth-order staggered- grid finite-difference method (Olsen, 1994) and the Southern California Earthquake Center (SCEC) Community Velocity Model Version 3.0 (CVM3.0), which is a 3D reference model used for regional waveform modeling (Magistrale et al., 2000). The agreement between the observed waveform and the synthetic waveform is reasonably good for the P wave and the “X” phase arriving at aboutlO sec except that the observed waveforms arrive earlier than the synthetics. Both phases are potentially useful for constraining basin structure. We will carry out a step-by-step application of the GSDF method to both phases to derive the phase-delay times and amplitude-reduction times. And we will interpret our Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 26 GSDF phase-delay measurements in terms of characters of our reference structure model CVM3.0. 241° 30' 34° 00' 242° 00' M^/C 242° 30’ DJJ GP A Rjp B[S P^U m C G ^ C PfiR ' RUS C£P " H i-T CljIN MJ.S P 1 S RPV 34° 00' 241° 30' 242° 00' 242° 30' Figure 2.3. The epicenter-station path from the 09/03/02 M L 4.8 Yorba Linda earthquake to CISN station DLA. The background color shows the thickness of the sedimentary basin. CISN stations are shown as black triangles. Focal mechanism as determined by Hauksson (2000) is plotted at the epicenter. Red line shows the epicenter-station path. The measurement process starts by constructing isolation filters f(t) [Figure 2.4(a)]. In previous studies of teleseismic data, isolation filters were computed by summing normal modes weighted according to their group velocities. In our analysis of local broadband data, we simply apply a time window on the complete synthetic seismogram and isolate the part we are interested in. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 27 s(i). m f(t) TO 20 S A A A A A A A A A A A A > i r > O O Q O o ° o o C o o O (f) 1 .2 0.4 0.6 O .f Figure 2.4. The GSDF processing for the P wave, (a) the observed seismogram, the complete synthetic and the isolation filter, (b) Solid line: the autocorrelation of the isolation filter; dash line: the windowed autocorrelation, (c) Solid line: the cross correlation between the isolation filter and the observed data; dash line: the windowed cross-correlation, (d) Solid line: the narrowband filtered windowed autocorrelation; dash line: the best-fit Gaussian wavelet, (e) Solid line: the narrowband filtered windowed cross-correlation, central frequency 0.2 Hz, bandwidth 0.02 Hz; dash line: the best-fit Gaussian wavelet, (f) Circles: 8tp\ triangles: < 5 tq. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 28 0.04 0.02 s (t) 0 -0.02 s(t) -0.04 -0.06 -0.08 f i t ) -0.1 -0 .1 2 1 0 5 1 5 2 0 25 \ i I I (d) i i i I \ W , , , j f v ^ U I i i ! \ I ! V V Figure 2.5. The GSDF processing for the X phase, (a) the observed seismogram, the com plete synthetic and the isolation filter, (b) Solid line: the autocorrelation of the isolation filter; dash line: the windowed autocorrelation, (c) Solid line: the cross-correlation between the isolation filter and the observed data; dash line: the windowed cross-correlation, (d) Solid line: the narrowband filtedash line windowed autocorrelation; dash line: the best-fit Gaussian wavelet, (e) Solid line: the narrowband filtedash line windowed cross-correlation, central frequency 0.2 Hz, bandwidth 0.02 Hz; dash line: the best-fit Gaussian wavelet, (f) Circles: r»p; triangles: Stq. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 29 The discrepancy between the isolation filter and the observed waveform can be concentrated by computing the cross-correlation Cfs(t) [Figure 2.4(c)]. The peak of the cross-correlagram is displaced from zero-lag by about - 0.8 second and the shape of the cross-correlagram is quite different from the autocorrelation of the isolation filter, CJf (t) [Figure 2.4(b)], which is an even function centered at zero- lag. The phase delay measured from the broadband cross-correlagram peak is the basis for a class of signal processing techniques that are commonly employed to estimate differential travel-times of finite-frequency waves (Su & Dziewonski, 1992; Masters et al., 1996; Dahlen et al., 2000). The major problem with these techniques is that the shape of the waveform is not well utilized. In GSDF analysis, this problem is overcome by a series of localization operations both in time domain and in frequency domain. To reduce the interferences from other arrivals on the seismogram, we multiply the cross-correlation with a time window centered at the peak of the cross- correlagram. The windowed cross-correlation WCfS (t) = W(t)CfS (t) [Figure 2.4(c)] and the auto-correlation WCff(t) = W(t)Cff(t) [Figure 2.4(b)] are then narrowband filtered at a set of frequencies 0), that sample the bandwidth of the windowed cross-correlagrams. Typically, we use zero-phase Gaussian filters of the form Fi(co) ~ exp[— (| to \ -ry, )2 / 2cr21 , where the ratio o i la>i « \ . The windowed and narrowband filtered cross-correlation FtWC fs(t) = Fi(t)*[W(t)Cfs(t)],i = 1,2,...I [Figure 2.4(e)] and autocorrelation Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 30 FjWCff [Figure 2.4(d)] can be well approximated using five-parameter Gaussian wavelets of the form g(t) = A Ga[ < r s(t - tg)] cos[cos(t-tp)]. (2.1) Here, Ga(x) = exp(-x2 /2). Two of the five parameters, the amplitude A and the phase tp are of our interest. The five parameters of the Gaussian wavelet model are estimated by fitting the expression (2.1) to the windowed and narrowband filtered cross-correlagrams using a weighted lease-squares method. We seek the wavelet g(t) that minimizes the quantity , f[F]WC fs (f) - g(t)]2 exp [ - f i t - ti)2]dt X = J TT------------- ;---------;--------------• (2-2) j g 2(t)exp[-r2( t - t i)2]dt The weighting is controlled by which is usually chosen close to the peak of FiWCfS (t), and y, which usually approximates o,. The fitting is usually excellent where the amplitude is large [Figure 2.4(de)] and a well-defined global minimum can be obtained. If the amplitude and phase delay for FiWCfS {t), as determined from such an algorithm, are A and tp, and for Ft WCff , are A and lp, the GSDF measurements of phase-delay time and amplitude-reduction time are defined as Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 31 (2.3) dtg = - In (2.4) The frequency-dependent GSDF measurements for the P wave are shown on Figure 2.4(f), and the measurements for the X phase are shown on Figure 2.5(f). 2.3.2 The Goodness of Waveform Fitting Minimizing GSDF measurements of Stm is consistent with matching waveforms point by point given enough frequency bandwidth and time span. Generally speaking, minimizing certain measures such as the L2 norm of GSDF measurements is a much more efficient way to reduce waveform misfit than direct waveform fitting techniques based on minimizing certain norms of the differential waveform Sf(t) = f ( t ) - f ( t ) . We illustrate this point using the P wave as an example. The observed and synthetic waveforms for this P wave are shown on Figure 2.6(a). The L2 norm of the differential waveform ,":l 2 is about 0.13 [Figure 2.6(b)], which is even larger than the L2 norm of the observed waveform ||/[ ? ) ||l 2 = 0.1. We decompose the frequency-dependent phase-delay time St'p [Figure 2.4(f)] measurements into a frequency-independent mean 8 t = < 8t' > and a Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 23 32 frequency-dependent residual Si5t' = St' - S t . We perturb the synthetic waveform f (t) to match the observed waveform f(t) in four steps. First, we correct S t by multiplying the Fourier transform of the synthetic waveform f(a )) with a phase shift exp\ico8tp]. As shown on Figure 2.6(d), reduced dramatically to about 0.06. Second, we apply the same procedure to amplitude measurements and correct for the frequency-independent mean of amplitude anomaly. The result is shown on Figure 2.6(ef), ]"L2 reduced to about 0.037. Third, we correct the frequency-dependent residual SS tp by multiplying f(aj)exp[ia)St] with a phase shift exp[ia)SSt (co)], where SS t (co) is a smooth function that interpolates SS tp at those frequencies where we made our GSDF measurements. The result is shown on Figure 2.6(gh), w o n l 2 reduced to about 0.02 and the perturbed synthetic waveform can model the shape of the observed waveform much better. And finally, we correct for the frequency-dependent residual of amplitude anomaly, ||< 5/(0||l2 reduced to about 0.001. We obtain an almost perfect fit to the observed waveform [Figure 2.6(ij)]. One important aspect of GSDF analysis is that, as pointed out in Gee & Jordan (1992), large travel-time shifts and large amplitude perturbations are accommodated by the linearized equations, provided that the differential variations with frequency is not too large. In most applications involving compact waveforms, the requirement for GSDF analysis is much less restrictive than Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 23 23 33 linearized waveform inversion based on minimizing the amplitude of differential seismograms. 0.04 0.02 - 0.02 -0.04 9 10 5 6 7 8 0.04 0.02 -0.02 ■0.04 8 9 10 5 6 7 0.04 0.02 - 0.02 -0.04 8 6 7 9 10 5 0.04 0.02 - 0.02 -0.04 6 7 8 9 10 5 0.04 0.02 - 0.02 -0.04 5 6 7 8 9 10 0.04 0.02 0 - 0.02 -0.04 5 6 7 8 9 10 0.04 0.02 0 - 0.02 -0.04 6 5 7 8 9 10 Time (s) Time (s) Figure 2.6. W aveform fitting by minimizing GSDF measurements. Left column: the observed (solid line) and synthetic (dash line) waveforms. Right column: the differential waveform e.g. the observed waveform minus the synthetic waveform. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 34 2.3.3 Interpretation of Frequency-dependent Phase-delay Times GSDF measurements can be related to perturbations of model parameters using linearized equations. In the second part of this chapter, we will show how to derive such equations using the Bom approximation (Marquering et al., 1999). The relation between Stp and perturbations of P wave velocity a(r) can be expressed as (2.5) Figure 2.7. The sensitivity kernel for the P wave measurements at 0.2 Hz (left) and 1.0 Hz (right). Top row: map view at depth 8.0 km. Bottom row: looking perpendicular to the source-receiver plane. Blue: negative sensitivity; red: positive sensitivity. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. C:+54C 35 Here, Kp(r, mo) is the sensitivity kernel for P wave velocity perturbation Sa(r) at position r, which is computed using the reference model mo- The integration is over the whole modeling volume V. The sensitivity kernels for the P wave phase-delay time at 0.2 Hz and 1.0 Hz are shown on Figure 2.7. First of all, we see strong sensitivity in the basin underneath the station, which suggests the measured phase-delay times of about negative 0.9 second at 0.2 Hz and 0.5 second at 1.0 Hz are likely acquired from the basin beneath the station. Second of all, a close examination of these two kernels shows that the kernel for the 0.2 Hz measurement is wider than the kernel for the 1.0 Hz measurement. As predicted by Dahlen et al. (2000), scattering, wave front healing and diffraction effects render the travel-time of a finite-frequency wave sensitive to 3D structures off the ray path. At lower frequency, the wavelength is longer, Figure 2.8. The sensitivity kernel for the X phase measurements at 0.2 Hz (left) and 1.0 Hz (right). Top row: map view at depth 8.0 km. Bottom row: looking perpendicular to the source-receiver plane. Blue: negative sensitivity; red: positive sensitivity. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 36 thus the seismic wave is sensitive to a larger area off the ray path. The frequency- dependence of the sensitivity kernel is more obvious for the X phase (Figure 2.8). At different frequency, the kernel is sensitive to quite different areas. We can also gain some insight about how the wave is propagated within complicated 3D structure models by examining this type of kernels. For example, the X phase is probably generated from scattering caused by heterogeneities around the basin near the surface. Measurements made on such phases are potentially useful for constraining shallow structures within the basin. 2.4 Theoretical Development The two localization operations, windowing in time domain and narrowband filtering in frequency domain, can isolate the desirable signal from surrounding noise, but can also distort the signal itself, so we need to account for their effects on the target wave-group in our analysis. The original GSDF analysis as formulated in Gee & Jordan (1992) corrects windowing and filtering effects by solving a set of equations [Equation (56)-(59), Gee & Jordan, (1992)]. The corrections for dtP A [Equation (57), (59), Gee & Jordan (1992)] are proportional to the fractional shift in the centre frequency due to filtering, (a> f - coc)/cbf , where 0)f is the narrowband filtering frequency and Q )( is the centre frequency of the isolation filter. Although this factor vanishes at the centre of the band, its magnitude can approach unity near the low- and high- Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 37 frequency ends for broadband signals, and the correction, which is expected to be small, can actually be significant. The problem lies in the Gram-Charlier expansion of the autocorrelation Cff (co) , which is a local approximation and only accurate in the vicinity of centre frequency cdc. In the following, we re-formulate GSDF analysis in terms of discrete Fourier expansion. The resulting theory is more elegant mathematically and more stable for local broadband signals such as those recorded by CISN stations. 2.4.1 Effects o f W indowing and Narrowband Filtering on Target W aveform We denote the target wave-group on the observed seismogram as f(t). It could be and phase, such as the P wave, the S wave, surface waves or just an arbitrary segment that constitutes a number of phases. The isolation filter, / (t) , is constructed to model the target wave-group. The autocorrelation of the isolation filter Cff (t) and the cross-correlation between the isolation filter and the target wave-group C$(f) can both be expanded as a summation of cosine functions of different amplitudes, frequencies and phase shifts: Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Cff (0 = 7 (0 ® 7 (0 = X Aj cos[Q )j (t + r / )], j 38 (2.6) cff (0 = 7(0 ® /(0 = X cos[«.(t + TJ P)\. (2.7) All functions and variables here are real and we assume coj > 0 and t e (-« ,+ » ). For each digital signal, the values of £ O j can be chosen according to the sampling theorem of discrete Fourier transform. The functionals that can be related to perturbations of model parameters are the differential phase shifts f r J p = T J p - ? j , (2.8) and the amplitude anomalies 8Tj q = -\n (A j/A j )/0Jj. (2.9) We first consider the windowing operation in time domain. Following Gee & Jordan (1992), to isolate the signal, we apply Gaussian windows on both correlagrams. The windowed correlagrams are then: WCff(t) = exp — a 2( t - t )2 ^ XV X W ' Cff(t), (2.10) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 39 < 7 ( t - t )2 W V W J Cff (t). (2.11) Here, the centers of the time windows are t and tw and the half-widths are a ' W ' w and (7W . Bring Equation (2.1)-(2.2) into Equation (2.3)-(2.4) and transform to frequency domain, we obtain: WClf(0J) = 2 ^ - Aj {exp 2<r W J 2 * {CQ-Uj) 2d,,2 exp[ia)jTl - iTjco - (o.)] + exp (a)+Q)j) 2 d ,2 (2.12) exp[-ieojfj - itjo j + cO j)]} W J 2 ' 2a,.2 exp [iQ}.r] - i t j o j - + exp 0co + cO j) 2a 2 exp[- (2.13) i(OjT] p - i t w(a) + (Oj) ]} Here, ffl6 ( - “ d 00). and we have WC a (~0J) = conjugate[lEC/? (0))\ and WCff{-co) = conjugate[WC# (ry)]. The factor exp[-(ry- )2 H 2 a 2 )\ actually measures the spreading of the spectrum caused by the windowing in time domain. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 40 The spectrum at a particular frequency, say m , is mixed with the spectrum at surrounding frequencies and the larger < J w is, the narrower the time window is, the more serious such a mixing effect is. We now apply the Gaussian narrowband filter centered at a set of frequencies of interest a > i > 0 , with half-bandwidths cr > 0 : Ft (of) = exp (ry-ry,.)2 2 + exp (co+coy . 2cr< 2cr. After some algebraic manipulation, we obtain the windowed and narrowband filtered cross-correlagrams in frequency domain: (ffl, - 0 1 j ) 2 (0 1 - 0 1 j ) 2 (Ol+Ol'j)2 ■ia)j t i p - i t w (oi+ojj ) (ejj+ W j)2 (o i-tQ j)2 ( o i+ a ] f . 2 ( a 2 +/J:2) r 2(7: + e " ' ’\e ■iCOjT} p - i t w (CO+COj ) (2.15) F . W C / oj) (W j-W j)2 (O l- O J j)2 (a i+ c o ’ j )2 e w (a> i+a> j) 2 (0 1 -0 1 j ) 2 ( 0 1 + 0 1 j )2 2a]" - i a i j t 1 p - i t w (ai+ai] ) (2.16) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The bandwidth of the windowed and narrowband filtered waveform is cr/, the contribution from each frequency (O j to the final waveform is controlled by the factor exp[-(a>j - o^)2/ 2((jJ + cr2)], thus, the farther c O j is from a)i , the less the contribution from (O s is. Transform Equation (2.15)-(2.16) back to time domain and after some algebraic manipulations we obtain the series for the windowed and narrowband filtered cross-correlagrams: Fftcff(t) = e 2 Y j cosi<S'( (t + f ’ J ) 1 + A J cosf<5)'(r +r / )1, (2.20) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 42 Fl WCff(t) = e 2 ^A)cos[ffl((f + r,; )] + A"cos[ft>J(f-fr,')] • (2.21) i Here, A . = ^ A te 2 ( f l F » 2 + < r '2 ) < 7 4 = — v (a ij-iO j)2 2(aA + cr2) (2 .22 ) (W i+ a > j)2 & , ~ 2 (ff„ 2 + o - ,2 ) Aj = ^ A je cr 4 = fV 2(0-lv 2 + C T ,.2 ) (2.23) ^ = L7 ‘ ' p VV ’ <27 . < 2 > ; - OJ j T J = —- T 1 — P P C O ; C O : (2.24) „. < 2 7 r — — — 7 J - P 0 0: P 00: - 00: 00: T-‘ = ^T > " °> ] ' ^2 ~ ^2 C O : (2.25) The cross-correlagrams after windowing and narrowband filtering can be expressed as summations of Gaussian wavelets parameterized by the quantities defined in Equation (2.17)-(2.19) and (2.23)-(2.25). The first half of the series modulated by A} and A; accounts for the contribution of the positive side of the correlagram spectrum to the positive peak of the narrowband filter and the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 43 negative side of the spectrum to the negative peak of the narrowband filter. The second half of the series modulated by Aj and A; accounts for the cross contributions, that is the contribution of the negative side of the spectrum to the positive peak of the filter and the positive side of the spectrum to the negative peak of the filter. The sufficient condition under which we can ignore the second half is that A . / A . « 1 and A . / A . « 1, which gives exp « 1, exp 2o jio j] 2 " a + a w 1 « 1. (2.26) This condition is not necessary though, because when A. and A] are small enough to be ignored, Aj and A j, which are certainly not larger than A . and A ., can also be ignored. Intuitively, a smaller cr means a narrower bandwidth in frequency domain, thus the second half of the summation is smaller; a smaller a w or (7W means a wider time window in time domain [Equation (2.10)-(2.11)], thus a narrower bandwidth in frequency domain and a smaller contribution from the second half of the summation. 2.4.2 The Sum m ation of Gaussian W avelets The right-hand sides of Equation (2.20) and (2.21), which are summations of Gaussian wavelets weighted according to their frequencies, can be approximated Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 44 by a single Gaussian wavelet. How good the approximation is can be adjusted by changing cr , aw or aw . In practice, cr. should be less than about one-tenth of Q ) i to enforce a good approximation. We express the composite Gaussian wavelets as: Fl WCjf(t)~g(t) = Aex p 2a. cos (2.27) 2a.. cos k (*-*«,)]• (2.28) The parameters that define these composite Gaussian wavelets can be found by minimizing the quantity defined in Equation (2.2). And the GSDF measurements of phase-delay time and amplitude-reduction time are defined in Equation (2.3)- (2.4). If we consider FjWCff (t) as a perturbation from F.WCff (t), the differences between their Gaussian wavelet parameters can be related back to perturbations of the parameters of each constituent Gaussian wavelet. We first define some parameters here: Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 45 B t = -^-ex p A G (t ~ t ) I v vv g / < p 'j = (dj'jf'J - ajstp) + (8)v - 5) )tg, (2.29) S. = ^ -e x p A ^ = (a>]f "j - a)s tp) + (0)s - )tg. (2.30) Let C and S be two row vectors with their elements defined by C = [B] cos (pj, B] cos q)j ], (2.31) S = [Bj sin < p 'j, Bj sin q >\]. (2.32) Then we have: Stq = C -Stq + S -d tp, (2.33) Stp = C -Stp - S - S t q. (2.34) Here, 8tq and 5tp are column vectors with their elements defined as: ln(A, /A ;) ln(A / A ) 8t=[-------L - ,----- (2.35) (O , \0)} Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 46 Stp = \r'l - f j ;t"J - f pJ]. (2.36) Simplifications can be achieved if we make some assumptions about our GSDF measuring process. If we assume a w = a w, Equation (2.10)-(2.12) give us: <7, ( T . a)j= 0)j. ajj - (0r (2.37) Considering Equation (2.8)-(2.9), we have Stq = & S r i ^ 8 r J q], (2.38) ( O j \ Q J j | a), oo-a). _ a): (O:-co, Stp = [ ^ S r J p ; (,tw — tw ); vSTp - t j ] . (2.39) (O j a)j (o! a), If we further assume 7 w=tw=0, e.g., we center both time windows at zero lag of the correlagrams, we have ft) , . f t ) . Stp ( O j ( O j (2.40) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 47 For the autocorrelation Cff(t), we have f ! = f j = T i -0, thus we know tp = tg = 0, so we have b = ^ (2.4i) A A and <Pj=<p"j= 0, (2.42) the two row vectors thus become C = [B),B]], S = 0. (2.43) In this case, we achieve full decoupling between the phase delay time and the amplitude reduction time, Equation (2.33)-(2.34) become & ,= C - 5 tq, (2.44) * „ = C - 8 t p. (2.45) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 48 Equation (2.44)-(2.45) is the perturbation form of the Wavelet Summation Theorem [Appendix E, Gee & Jordan (1992)]. These two equations are excellent /2 2 approximations when the effective bandwidth ^jai + criv is small enough. In practice, we use a cosine time window, which has a flat part in the middle, rather than a Gaussian time window. The effect is that the target wave-group is disturbed as little as possible while the extraneous phases are excluded from the time window. The effective bandwidth is actually controlled by c r. To test the applicability of Equation (2.44)-(2.45), we conducted a numerical C0:-a)0 experiment. In Equation (2.6), we specify A = Ga[--------- ] and f p = 0. We assume a linear dependence of SrJ pq on frequency: Srq - acOj l(ln ) , Srp = P a> i K in) • We chose a)0 = 0.6n , cr0 = 1.0?r, a - 2, [3 = 1. The amplitude spectrum of Cff and Cff are shown on Figure 2.9, their time domain waveforms are shown on Figure 2.10. For the windowing and narrowband filtering parameters, we chose < 7 W = 0.2, cr. = 0.027T, and Q )i ranging from 0.2n to 2.0n . The GSDF measurements of StP A and the predicted measurements made using the specified S rJ pq and the right-hand sides of Equation (2.44)-(2.45) are plotted on Figure 2.11. The vectors C for the 10 measurements on Figure 2.11 are shown on Figure 2.12 as functions of frequency. Significant deviation of the predicted value from the measured value starts at Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 49 about 0.9 Hz, where Stp is more than 70% of the period of the wave we are measuring. E 5 o < D Q _ “ 0.5 © 73 3 Q. E < a 0.5 frequency (Hz) Figure 2.9. The amplitude spectrum of C g (solid line) and C g (dash line). The GSDF amplitude-reduction anomaly dtq measures amplitude difference between two waveforms. Another important aspect of Equation (2.44)-(2.45) is that, as pointed out in Gee & Jordan (1992), large travel-time shifts and large amplitude perturbations are accommodated by the linearized equations, provided that the differential variations with frequency are not too large. In most applications involving compact waveforms, the requirement for GSDF analysis is much less restrictive than linearized waveform inversion based on minimizing differential seismograms. One example is the P wave on Figure 2.4 and Figure 2.6. The amplitude of the differential waveform is even larger than the observed waveform [Figure 2.6(ab)]. It is not a small perturbation from the synthetic waveform. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 50 Direct waveform inversion based on perturbations of waveform cannot be easily applied. In GSDF analysis, we shift and scale the synthetic waveform into the right position to minimize the misfit rather than simply minimize the amplitudes of the differential waveform. 1 0.5 0 0 5 -5 time (s) Figure 2.10. The time-domain waveform for C g (solid line) and C g (dash line). C g is dispersed relative to C g . 2.4.3 Sensitivity Kernels for Tomography In this part, we are going to relate our frequency-dependent GSDF measurements of Stv,q back to perturbations of seismic velocity. The sensitivity kernel is dependent upon the parameters we choose for our windowing and narrowband filtering process. These time-frequency-dependent sensitivity kernels will tell us what we can see at different time windows and frequencies. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 51 The autocorrelation of the isolation filter can be expressed in frequency domain as: Cff(co) = f\c o )f(o j). (2.46) The observed waveform can be express as a perturbation from the isolation filter as f(co) = f(co) + df(co). (2.47) 1 A cr d . % 0.8 0.4 0.4 0.2 0.8 1 Figure 2.11 GSDF measurements made by fitting five-param eter Gaussian wavelets to the windowed and narrow-band filtered cross-correlagrams (solid lines) and the prediction made using Equation (2.44)-(2.45) (dash lines). Circles: 6tp ; triangles: 8tq. Significant deviation of the predicted value from the measured value starts at 9 Hz, where it is difficult to make accurate measurements due the small signal to noise ratio for C g (Figure 2.9). Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 52 Here, Sf can be related to perturbations of seismic velocity structures using the standard Bom approximation (Marquering et al., 1999). Thus, the cross correlation between the isolation filter and the observed waveform can be express as: = f\(o)f{co) = Cff{co) + f\cQ)8f{co) . (2.48) On the other hand, we have 0.075 d 0.025 0.2 0.3 0.4 0.5 0.6 0.7 0.8 frequency (Hz) Figure 2.12. The kernel C for the 10 5tpq measurements on Figure 2.11. They are Gaussian functions of frequency. Their centoids and widths are controlled by windowing and narrow band filtering parameters. At each narrow-band filtering frequency c o s, only 5Xpq at surrounding frequencies have significant contributions to the measurement 5tpq made at frequency o),. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Cff (co) = C , (00) exp[i co8xp (co) - coSrq (cd)\, 53 (2.49) where Srp(co) and Srq((D ) are the same as the phase-delay time and the amplitude-reduction time as defined in Equation (2.8)-(2.9) except that here we express them as continuous functions of frequency. Under the assumption that these perturbations are small, we can linearize the exponential as: exp[i<w5rp (co) - coStq (oj)] = 1 + i coSrp (co) - coSrq (co) (2.50) and consider Equation (2.48) we can arrive at: Sra(co) = — ^-Re q co / » C ff (co) Sf(co) (2.51) (co) = — Im co f \ c o ) Cs (C D ) Sf(CD) (2.52) Based on the results of the last section, we can express the GSDF measurements of Stm using 5tva(c d ) : Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 54 + ^cr'icoCJf(co) Stq = j exp cr,.,Q ) A (co — ( J t )i ) 2 (a w2 +cr 2) St (co)dco + P < j , . t y C # ( o j ) + ----- „ , -ex p o I 0) | A (CO+COt ) W + r r , 2) (2.53) Srq (co)dco ^ a tcoCff(co) Stp = J exp „ (7w (oA + ° } < 7 '.(O C ff (C O ) + „ ~ exp n CT,.,CO A (a)-Mi? 2(ow2+ e 2) 5 rp (co)dco (co+cot) 2 ( c r j + c r 2) (2.54) S rp (co)dco Here, analogous to Equation (2.11)-(2.12), we have c o = ■ o 2 c o + c t 2 o j i 2 , _ 2 cr,., + cr. (2.55) co = a 2 co - crjcpj cr,2 + cr,.2 (2.56) Bring Equation (2.49)-(2.50) into Equation (2.51)-(2.52) and consider Cff(co) is real, we obtain: Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 55 St = - I"— Re] exp n 1 (co-oo^ 2 2 itrJ+v,2) f \ c o ) 00 Sf(co) \doo ~ J : t t r T xp (O O + CO^ 2 (^w 2 +CT,2) f { ^ S f(o o )\d c o 0 0 + 7 & St = — ^ Im < e x p o * w A { C T + J— ^Im jexp (|(0-00,.) 2 « T w2 + ° i 2\ f {^ -S f(oo)\dco C O 0<V* (oo + oo^2 2(0-w 2 + ^ ,2). f ( P -S f(o o )\d o j 00 We can express Sf(co) using its time domain expression through ; transform: -t- oo Sf(oo) = jSf(t) exp(icot)dt. Bring Equation (2.59) into Equation (2.57)-(2.58), we obtain: (2.57) (2.58) Fourier (2.59) (2.60) (2.61) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1 . 0 H z 0.6 H z 0.2 H z 56 2 0 ■ 2 5 10 15 2 1 0 ■ 2 5 10 15 2 0 ■ 2 10 5 15 2 1 0 •2 2 0 •2 5 10 15 2 0 ■ 2 5 10 15 Time (sec) Time (sec) Figure 2.13. The I(t) function for the P wave isolation filter shown on Figure 2.4a. The narrowband filtering frequency is 0.2 Hz (top), 0.6 Hz (middle) and 1.0 Hz (bottom). The real part is shown on the left column, the imaginary part is shown on the right column. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 57 Here, the weighting function 7(f) is defined as: exp(iax)da> exp(iax)dco (2.62) H(a>) is the Heaviside function. When the effective bandwidth narrow enough the second integral can be ignored. I(t) is a compact, oscillatory function with its width controlled by the effective bandwidth. An example of lit) is shown on Figure 2.13. The width of the corresponding sensitivity kernel is also controlled by this effective bandwidth, the oscillatory character of the kernel depends on the narrowband filtering frequency (O i as shown on Figure 2.7-2.8. 2.5 Discussion In this chapter, we have reformulated and extended the original GSDF analysis (Gee & Jordan, 1992) to better utilize regional and local broadband waveform information. One important difference from the original formulation is that instead of correcting raw GSDF measurements of r)7 p ,q for windowing and filtering effects to get unbiased 6 tva and computing sensitivity kernels for < 5 rP ;q , we incorporate the metadata generated from the measurement process into the sensitivity kernel through the weighting function I(t), which is derived from the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 58 isolation filter. Because the weighting function has compact support (localized) both in time and in frequency, we can achieve good resolution in one domain without sacrificing too much resolution in the other. The tradeoff between resolutions in two domains should conform to the Heisenberg Uncertainty Principle. In the following chapters, we will apply this new theory to various observational problems in seismology. In Chapter 4-5, we will derive sensitivity kernels for seismic source parameters and invert GSDF measurements for finite source properties. In Chapter 6, we will use the sensitivity kernels derived in this chapter for tomography problems to invert GSDF measurements for a 3D perturbation to the 3D reference model CVM3.0. As we will see in those applications, GSDF methodology is a very powerful tool to utilize broadband seismic waveform data. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 59 2.6 R eferen ces Akgelik, V., J. Bielak, G. Biros, I. Epanomeritakis, A. Fernandez, O. Ghattas, E. J. Kim, J. Lopez, D. O'Hallaron, T. Tu, and J. Urbanic, High Resolution Forward and Inverse Earthquake Modeling on Terasacale Computers. Proceedings of the ACM/IEEE SC2003 Conference November 15 -21, 2003 Phoenix, Arizona. Chen, P., T. H. Jordan and L. Zhao, Finite moment tensors of 3 September 2002 Yorba Linda earthquake, Bull. Seismol. Soc. Am. in press. Dahlen, F. A., S.-H. Hung and G. Nolet, Frechet kernels for finite-frequency traveltimes - 1. Theory. Geophys. J. Int. (2000) 141, 157-174. Dahlen, F. A. and A. M. Baig, Frechet kernels for body wave amplitude. Geophys. J. Int. (2002) 150, 440-466. Gaherty, J.B., T. H. Jordan and L. S. Gee, Seismic structure of the upper mantle in a central Pacific corridor. J. Geophys. Res. (1996), vol 101, No. B10, 22291- 22309. Gee, L. S., and T. H. Jordan, Generalized seismological data functionals. Geophys. J. Int. (1992) 111, 363-390. Ihmle, P. F., and T. H. Jordan, Teleseismic Search for Slow Precursors to Large Earthquakes, Science (1994), 266: (5190) 1547-1551, Dec 2. Katzman, R., L. Zhao, T. H. Jordan, High-resolution, two-dimensional vertical tomography of the central Pacific mantle using ScS reverberations and frequency- dependent travel times. J. Geophys. Res. (1998), vol 103, No. B8, 17933-17971. Luo, Y. and G. T. Schuster, Wave equation travel-time inversion. Geophysics (1991), 56, 645-653. Magistrate, H., S. Day, R. W. Clayton, and R. Graves, The SCEC southern California reference 3D seismic velocity model Version 2, Bull. Seismol. Soc. Am. (2000) 90, no. 6B, S65-S76. Marquering, H., F. A. Dahlen and G. Nolet, Three-dimensional sensitivity kernels for finite-frequency traveltimes: the banana-doughnut paradox, Geophys. J. Int. (1999), 137, 805-815. Masters, G., S. Johnson, G. Laske and H. Bolton, A shear-velocity model of the mantle, Phil. Trans. R. Soc. Lond. (1996), A354, 1385-1411. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 60 McGuire, J. J., L. Zhao, and T. H. Jordan, Teleseismic inversion for the second- degree moments of earthquakes. Geophys. J. Int. 145, 661-678, 2001. Nolet, G., Partitioned waveform inversion and two-dimensional structure under the network of autonomously recording seismographs. J. Geophys. Res. (1990) vol. 95. No. B6. 8499-8512. Olsen, K.B., Simulation of three-dimensional wave propagation in the Salt Lake Basin, Ph.D. Thesis (1994), University of Utah, Salt Lake City, Utah, 157 p. Sipkin, S. A. and T. H. Jordan, Lateral heterogeneity of the upper mantle determined from the travel times of multiple ScS. J. Geophys. Res., (1976) 81, 6307-6320. Su, W.-J. and A. M. Dziewonski, On the scale of mantle heterogeneity, Phys. Earth Planet. Inter. (1992), 74, 29-54. Tarantola, A., A strategy for nonlinear elastic inversion of seismic reflection data, Geophysics (1986), 51, 1893-1903. Tromp, J., C. Tape and Q. Liu, Seismic tomography, adjoint methods, time reversal, and banana-doughnut kernels. Geophys. J. Int. (2005) 160, 195-216. Zhao, L., T. H. Jordan, K. B. Olsen and P. Chen, 2005. Frechet Kernels for Imaging regional earth structure based on three-dimensional reference models. Bull. Seismol. Soc. Am, submitted. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 3 61 A Quantitative Evaluation of SCEC Community Velocity Model Version 3.0 Using Generalized Seismological Data Functionals 3.1 Abstract We present a systematic methodology for evaluating and improving 3D seismic velocity models using broadband waveform data from regional earthquakes. The operator that maps a synthetic waveform into an observed one is expressed in the Rytov form D(co) = exp[im<5rp(< y)-co(S T q(<y)]. We measure the phase delay time Stp(cO i) and the amplitude reduction time Stq(oJi) at discrete frequencies of interest cot using the GSDF method (Gee & Jordan, 1992; Chen et al., 2005b). Each GSDF measurement of c)rp ,q can be expressed as a linear integration of dTm (co) weighted by a Gaussian-like function of frequency, whose parameters are determined by the GSDF measuring process (Chen et al. 2005b). We have applied this procedure to SCSN recordings of 72 small events (3.0 < M l < 4.8) in the Los Angeles region. Synthetic seismograms were computed using 3 types of velocity models: the 3D SCEC Community Velocity Model version 3.0(CVM3.0) (Magistrale, et al., 2000), the ID Standard Southern California Crustal Model (SoCaL) (Dreger and Helmberger, 1993), and a set of path-averaged ID models (AID). We measured Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 62 more than 1000 body waves from 1073 source-station paths, at five frequencies ranging from 0.2 Hz to 1 Hz, yielding a data set comprising -20,000 dtP A (co) observations. Overall, the CVM3.0 provided a substantially better fit to the waveform data than either laterally homogeneous or path-averaged ID models. Relative to SoCaL, CVM3.0 reduced the variance in Stp by 61% to 0.375 s2, and reduced the variance in Stq by 64% to 0.162 s2. The variance reductions of CVM3.0 relative to AID were 55% and 53%, respectively. The variance reductions for S wave measurements are as large as for P wave measurements suggesting the scaling relations relating S wave velocity to P wave velocity in CVM3.0 are good approximations. Our analysis on 90 SV-SH pairs shows that the polarization anisotropy is very small, only 4 ms/km in the upper crust. Our measurements indicate that CVM3.0 wave speeds might be too low in Los Angeles Basin area. 3.2 Introduction For areas with complex subsurface structures such as the Los Angeles basin area (Figure 3.1), 3D seismic features can have large effects on ground motion. The SCEC Community Velocity Model Version 3.0 (CVM3.0) (Magistrale et al., 2000) serves as a unified 3D reference model for regional waveform modeling in Southern California. It is a combination of existing knowledge about geological structure, regional tomography and sonic logs in Los Angeles, San Fernando, San Bernardino, Ventura Basins and surrounding areas. The quality of CVM3.0 directly affects strong motion prediction and seismic hazard assessments as well Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 63 as seismic source studies in these areas. The purpose of this paper is to quantify the ability of CVM3.0 in predicting recorded broadband waveform data by comparing synthetic seismograms calculated using CVM3.0 with observed ones. To demonstrate the complexity of 3D heterogeneity in Los Angeles basin area, we 241“ 00' 241” 30' 242“ 00' 242” 30' 243” 00' Hollywood /Y o rb a Linda 241 ”00' 241” 30' 242” 00' 242“ 30' 243" 00' Figure 3.1 The 72 small events and the 44 Trinet broadband 3-component stations used in this study. The focal mechanisms were determined by Hauksson (2000). The stations are shown as black triangles. The background color shows the sediment depth in and around Los Angeles basin area. Red lines show the source-station paths used in this study. show a profile of CVM3.0 P wave velocity across the basin from the epicenter of the 09/03/2002 Yorba Linda event to station LAF (Figure 3.1). The wave speed varies from 339 m/s at the top to 7.95 km/s at the Moho (Figure 3.2a). The basin is clearly indicated by the low velocity area in the middle. To demonstrate the importance of 3D structural heterogeneity on the propagation of seismic waves, we compared synthetic seismograms computed using CVM3.0 with synthetics Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 64 computed using the ID Standard Southern California Crustal Model (SoCaL) (Dreger and Helmberger, 1993) (Figure 3.2b) and the path-averaged ID models (AID) (Figure 3.2b), which were extracted from CVM3.0 by averaging wave slowness along the source-receiver path. 1040 2000 3000 4000 5000 6000 7000 V ( k m /s ) p Figure 3.2 Examples of the three velocity models employed in our study, (a) The 2D profile of CVM3.0 P wave velocity along the source-station path from the epicenter of 09/03/2002 Yorba Linda event to station LAF. (b) The ID P wave velocity models along the same source-receiver path, solid line shows the A ID model, dash line shows the SoCaL model. Some examples of observed and synthetic seismograms for the 09/09/01 Hollywood earthquake (Figure 3.1) are shown on Figure 3.3. We compare the ground motion recorded at two different sites: station PAS, a rock site outside the basin, and station WTT, a basin site. The transverse component at PAS is dominated by a simple SH wave and all three structure models were able to model both the shape and the arrival time of this SH pulse at frequencies up to 1.0 Hz. The transverse component at station WTT is much more complicated, suggesting Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 3.3 The transverse components of the observed and synthetic seismograms for the 09/09/2001 Hollywood earthquake as recorded at station PAS and WTT. First trace from top of each panel shows the observed seismogram; second, CVM 3.0 synthetic seismogram; third, A ID synthetic; fourth, SoCaL synthetic, (ab): seismograms for rock site PAS. (cd) seismograms for basin site WTT. (ac): low- pass filtered to 1.0 Hz using a 4th-order Butterworth filter, (bd): low-pass filtered to 0.2 Hz using a 4th-order Butterworth filter. strong interactions between the waves and structural heterogeneities in or around the basin. At 0.2 Hz, the 3D synthetic was still able to model both the SH phase and the arrival at about 13 second, which did not show up on the two ID synthetics. At 1.0 Hz, all three models did a poor job modeling the recorded waveform. The SH arrival-time predicted by SoCaL is more than two seconds earlier than recorded arrival-time, suggesting SoCaL S wave velocity is too fast along this basin path. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 66 So far, regional tomography in Southern California has primarily based on the inversion of P wave onset time using geometric ray theory (Hauksson, 2000). Large amount of information carried by the waveform data is not utilized. The inversions for seismic source characters are largely based on direct waveform- fitting techniques such as “cut and paste” (Zhu and Helmberger, 1996; Graves and Wald, 2000). For a few model parameters, some searching algorithms can be used to cover the entire model space to find a global optimum. If the gradient method is employed, direct waveform fitting tends to lock the inversion into false local minima if the starting model is not close enough to the true model. The problems associated with seismic inversion are three-fold: first, how to parameterize the model so that the information in the observation is neither over interpreted nor wasted; second, how to quantify the misfit between observed data and model prediction, which is termed data functionals in the paper, so that the relation between the data functionals and model parameters can be well linearized; third, how to construct sensitivity kernels to relate misfit information to perturbations on model parameters and minimize the cost in computation and storage. To quantify the differences between waveforms, we adopt Gee and Jordan’s generalized seismological data functional (GSDF) technique (Gee and Jordan, 1992; Chapter 2), which provides a unified methodology for the analysis and inversion of regional waveform data. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 67 3.3 Methodology Regional seismic networks such as California Integrated Seismic Network (CISN) have been providing us numerous seismograms containing rich information about underground seismic structures and earthquake rupture processes. Progress of regional seismology will be depending on our ability to interpret these seismograms. At regional distances, near-surface 3D heterogeneity plays a more important role on wave propagation, which makes it difficult to separate direct, scattered, converted and surface wave phases. For moderate to large earthquakes, regional seismograms are further complicated by phases generated from finite faulting processes, such as the stopping phases. GSDF method is particularly suitable for the analysis of these complicated seismograms. Its robustness lies in its ability to enhance signal to noise ratio and it is flexible enough to utilize phases of any kind. It is adopted in our study to fulfill three goals: First, provide a unified methodology for the analysis and inversion of regional waveform data; Second, make measurements relative to synthetic seismograms that incorporate existing knowledge about seismic source and structure; Third, invert the measurements for source and earth structure using full 3D methodology. GSDF methodology separates the process of waveform fitting from the inversion and provides a quantitative assessment of data quality before the inversion. Instead of fitting the wiggles, we extract generalized functionals such as frequency-dependent phase and amplitude anomalies from waveforms and derive Frechet kernels that linearly relate these functionals to the model parameters. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 68 Small errors in the Frechet kernels won’t cause instability in the inversion. Moreover, since it is based on the gradient method, it is much more efficient than certain searching algorithms when the number of parameters are large. The GSDF methodology has been successfully applied for inverting seismic source parameters using regional broadband seismograms recorded by CISN (Chen et al., 2005a). Its full potential will be demonstrated in the evaluation and refinement of the seismic structure model of Southern California described in Chapter 6. 3.3.1 GSDF Processing GSDF analysis is a type of windowed-Fourier analysis. It consists of 6 steps: (1) synthesize the target wave-group into an isolation-filter f i t ) , which could be direct P, S, surface waves or any phases on the full synthetic seismogram s(t). In previous studies of global tomography problems (Gaherty and Jordan, 1996; Katzman, et al, 1998) and global studies of seismic sources (McGuire et al., 2001, 2002), isolation-filters were computed by summing normal-modes weighted according to the group and phase velocities of the target phases. In this study, we obtain the isolation-filters by windowing the target wave-groups on the full synthetic seismograms. (2) Cross-correlate the isolation-filter f i t ) both with the full synthetic seismogram s(t) and with the observed seismogram sit) to obtain the auto-correlagrams Cfs (t) and the cross-correlagram Cfs it). Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 69 (3) Multiply the crosscorrelagrams with a time window centered at the peaks of the correlagrams to obtain WCfs (f) and WCfs (t). This operation reduces interferences from other arrivals on the seismograms since near peak the correlagrams are dominated by the phase within the isolation filter. (4) Narrow-band filter the windowed correlagrams using a set of zero-phase Gaussian filters F ] with centre frequencies (t)i and bandwidths c i. (5) When certain conditions about windowing and narrow-band filtering are enforced, the resulting waveforms FiWCfs (t) and FiWCfs (?) can always be well matched by five-parameter Gaussian wavelets in the form FiWCfs(t) = A Ga[<rs(t - tg)]cos[5)s ( t - t p)] (3.1) for the synthetics and F.WCfs (t) = A Ga[crs (* - tg )]cos[tys (t - tp)] (3.2) for the observed ones. Here, Ga[v] = exp(-v2 / 2). In practice, we find the five parameters of the Gaussian-wavelet model by fitting the time-domain expression (right-hand side of Equation (3.1) and (3.2)) to the windowed and narrow-band filtered correlagrams (left-hand side of Equation (3.1) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 70 and (3.2)) using a weighted least-squares method. (6) The differences between observed and synthetic waveforms can be parameterized by two frequency-dependent time-like quantities: the phase delay anomaly, StP = tp - t p ’ (3-3) and the amplitude reduction anomaly, 6tq = -(In A - In A) / a> s. (3.4) These two frequency-dependent data functionals measure the differential dispersion of the observed waveform relative to the synthetic one. The phase-delay time and amplitude-reduction time resemble travel-time and t* in traditional data analysis except that GSDF measurements of dtm are frequency-dependent and they are measured relative to synthetic waveforms. Since this is the first time we discuss its application to regional seismology, we illustrate GSDF processing with two examples. Figure 3.4a shows the observed seismogram from 09/09/01 Hollywood earthquake as recorded on the radial component of station GR2 (Figure 3.1) and the complete synthetic seismogram computed using CVM3.0. Both seismograms have been low-pass filtered using a Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 71 fourth-order Butterworth filter with comer frequency at 0.2 Hz. At such a low frequency, there is no separation between direct P, SV and surface waves. We thus employ the whole wave-group as our isolation filter and cross-correlate it both with the observed record and with the complete synthetic seismogram. The cross- correlagrams are then windowed using a Gaussian taper with half-width of 10 seconds (Figure 3.4b) and narrow-band filtered at frequencies ranging from 0.1 to 0.3 Hz. The resulting waveforms can be well approximated by five-parameter (a) (b) (c) x 1Q'a Time (second) -20 10 0 1 Time (second) 10 20 0.2 0.1 0.3 frequency (Hz) Figure 3.4 An example o f GSDF analysis, (a) shows the radial component of the observed and synthetic seismograms for the 09/09/2001 Hollywood event as recorded at station GR2 as well as the isolation filter obtained from the synthetic. All seismograms are low-pass filtered to 0.2 Hz using a 4th -order Butterworth filter, (b) The cross-correlagram between the isolation filter and observed seismogram is shown on top, the cross-correlagrams between the isolation filter and the synthetic seismogram is shown at bottom. Dash lines show the windowed cross-correlagrams. (c) shows the cross-correlagrams narrow-band filtered at 0.2 Hz. (d) shows the GSDF measurements, circles: phase-delay time; triangles: amplitude-reduction time. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 72 Gaussian wavelets (Figure 3.4c). Our measurements of differential phase-delay and amplitude-reduction anomaly are then computed using Equation (3.3)-(3.4). Figure 3.4d shows differential dispersion of the observed waveform relative to the synthetic one. At 0.1 Hz, the wave arrives about 1.1 seconds earlier than predicted by CVM3.0; at 0.3 Hz, it arrives about 0.4 second earlier. GSDF allows us to make measurements on any wave-group without identifying its type. The complicated process that generates this wave-group is accounted for in our sensitivity kernel (Zhao et al, 2005), which is computed by solving the wave equation using certain numerical schemes such as finite-difference or spectral-element method. This is particularly advantageous in analyzing regional and local seismograms, which are characterized by small separations between different arrivals. Using GSDF, we do not have to separate these arrivals and invert them separately; we can invert the whole wave-group that contains many different arrivals, because the interferences between these arrivals have been correctly accounted for in our sensitivity kernel (Zhao et al, 2005). To show the differences between different synthetic waveforms, we compare the observed record, CVM3.0 synthetic, AID synthetic with SoCaL synthetic. Figure 3.5ab displays the observed seismogram from 09/03/02 Yorba Linda earthquake as recorded on the vertical component of station MWC (Figure 3.1) and its synthetics computed using CVM3.0, AID and SoCaL. All seismograms have been low-pass filtered using a fourth-order Butterworth filter with comer Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 73 frequency at 1 Hz. We intend to compare the P waves on the observed and synthetic seismograms. We use SoCaL synthetic as our reference and obtain the isolation filter by multiplying the full SoCaL synthetic with a time window A ID CVM3.0 Data -0" ' SoCal x 10 h 2 o - & ■ — 2 3 4 5 0 . S 0.6 0.7 0.8 05 1.1 02 0.3 0.4 frequency (Hz) Figure 3.5 An example of using GSDF analysis, (a) shows the observed waveform, CVM3.0 and A ID synthetic waveforms, (b) shows the SoCaL synthetic and the isolation filter, (c) shows the GSDF measurements, upper panel: phase-delay time, lower panel: amplitude-reduction time. centered at the P wave (Figure 3.5b). This isolation filter is then cross-correlated with the observed record, the CVM synthetic, the AID synthetic as well as the SoCaL synthetic itself. The cross-correlagrams are then windowed and narrow band filtered at frequencies ranging from 0.2 Hz to 1.2 Hz. On Figure 3.5c we can see CVM3.0 synthetic has similar frequency-dependence as the observed waveform, while the two ID synthetics do not match the observed curves very well. This simply indicates that CVM3.0 synthetic waveform match the observed Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 74 waveform better than the two ID synthetics, which can be seen from the seismograms on Figure 3.5ab. 3.3.2 Sensitivity Kernel GSDF measurements of StP A can not only parameterize waveform differences but also be related back to perturbations on model parameters. In GSDF analysis, we can represent the mapping of a synthetic waveform u into an observed waveform u using a differential operator u(a))=D(co)u(a) (3.5) The differential operator can be expressed in Rytov form using two frequency- dependent, time-like quantities < 5 rp,q(co) D(co) = exp [ia>Srp (co) - (o&rq (<«)]. (3.6) As shown in Zhao et al. (2005), Chen et al. (2005b) and Chapter 2, JTp ,q (ft> ) can be expressed as linear combinations or integrations of perturbations of model parameters. The coefficients in front of model parameter perturbations are usually termed “Frechet (sensitivity) kernels” or “Frechet derivatives”. For tomography problems, the linearized relation is in the following form: Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 75 Stx{(o) = ^Kx° (r,a))Sm(r)d3 r , (3.7) where, subscript x represents p or q, dm is velocity perturbation we are trying to invert for and K is the sensitivity kernel, which is a function of space r and frequency and can be computed given a reference model mo. As shown in Chen et al. (2005b) and in Chapter 2, each GSDF measurements of r)7p,q(cO i) can be expressed as an integration of SrP A (co) weighted by a Gaussian-like function: Here, Ix(co) is a bell-shaped function whose width and center are controlled both by the spectrum of the isolation filter and the parameters we choose during making the GSDF measurement. The center of Ix(a> ) usually lies close to narrow band filtering frequency & > , and its width is usually controlled by the narrow-band filter width The relation between Stx and model parameter perturbations is thus: dtx = J / x{(0)Stx(0i)d(O. (3.8) (3.9) The sensitivity kernel for Stx, which is Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Lx (r) = j l x(co)K”°(r,co)dco, 76 (3.10) can be calculated either in frequency domain or in time domain. Some examples of such time-frequency-dependent sensitivity kernels can be found in Chen et al. (2005b), Zhao et al. (2005), Chapter 2. For Stp measurements, these kernels may have a “banana-donut” shape, whose decay and oscillatory characters depend on Uco)- GSDF measurements of dtm can be directly inverted for model parameter perturbations using the sensitivity kernels described above. As shown in Chen et al. (2005b), Chapter 2 and Gee & Jordan (1992), equation (3.8) is applicable even to large travel-time or amplitude perturbations as long as the differential variation with frequency is not too large. In practice, this means we can pick any phase on the seismogram for our GSDF analysis no matter how far apart the synthetic and the observed waveforms are or how big the constant amplitude difference between them is as long as the differential variations in the shapes of the two waveforms are not too large. In linearized waveform inversion (Tromp et al., 2004), the synthetic waveform is simply subtracted from the observed waveform and the resulting differential waveform is back-projected to model parameter perturbations. Even small shifts in travel-time can cause large amplitudes of the differential waveform and a linearized gradient optimization such as Conjugate- Gradient might not be effective in solving the inverse problem. In GSDF analysis, Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 77 the inversion procedure tries to shift and scale the synthetic waveform into the right position to minimize the misfit rather than simply minimizes the amplitudes of the differential waveform. In most inverse problems involving compact waveforms, GSDF analysis is much less restrictive than linearized waveform inversion using differential waveforms. In traditional cross-correlation technique for travel-time tomography, arrival-time difference between the synthetic and observed waveforms is determined by the location of the broadband cross- correlagram peak. The differences between the shapes of the waveforms (e.g. differential dispersion between two waveforms) are ignored. In GSDF analysis, we take into account the frequency-dependent nature of waves, the waveform information is well utilized. Because GSDF measurements can quantify waveform differences and because they can be related to perturbations in model parameters, we can use GSDF measurements to examine the quality of the reference model mo in predicting observed ground motion. One example is shown on Figure 3.6. The seismic waves generated by the 09 Sept. 2002 Yorba Linda earthquake propagate across the Los Angeles basin and arrive at station LAF. The sensitivity kernels, which are calculated using CVM3.0 and K Olsen’s finite difference code (Zhao et al., 2005), show complicated and intensive interactions between the seismic waves and the heterogeneities in and around the basin, especially for later arrivals such as the kernel for S wave. The GSDF measurements associated with this kind of kernels encode a large amount of information about the basin structure, when Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 78 properly interpreted, some preliminary data analysis on < 5 fp > q (coi) can give us some idea about the general quality of the reference model. Z-component Z-component Pj\ 20 5 10 15 _AF time (sec) free surface j \ F time (sec) free surface ww*r."<-n» * — 6 X n > ’ s 5tp Z-component Y-component I I ! I 1 ------ 1 ------ 1 ------1 ------ 1 ------1 ------1 ------ 1 ------ 1 ------ L- 10 15 20 time (sec) free surface 10 15 20 time (sec) free surface ^ t 5x,j Figure 3.6 Frechet kernels for phase-delay tim e at station LAF from the 09/03/2002 Yorba Linda earthquake. The first three kernels are for the P-wave velocity for the time windows highlighted in red in the vertical component waveforms. The last plot is the kernel for the shear-wave velocity for the S wave arriving around 22 sec on the north- south com ponent waveform. 3.4 GSDF Data Analysis GSDF analysis can take full advantage of a limited data set. On a single seismogram, we can make measurements on different arrivals; for a single arrival, we can make measurements at different frequencies. Zhao et al., (2005) and Chen et al. (2005b) and Chapter 2 show that sensitivity kernels for GSDF Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 79 measurements are both time dependent and frequency dependent. Using GSDF analysis, even a few seismograms can provide us a large amount of information about underground seismic structure. We collected about 4000 CISN 3-component broadband velocity seismograms from 72 small earthquakes with local magnitude ranging from 3.0 to 4.8 in the Los Angeles region. The sampling interval of these seismograms is 0.05 s. The north and east components were rotated into the radial and transverse components. We kept only the records with a good signal-to-noise ratio. The locations and focal mechanisms of these events as determined by Hauksson (2000) and the 44 SCSN stations used in this study are shown on Figure 3.1. The focal depth of these events varies from 3 km to 21 km. The Los Angeles basin is well covered by our source-station paths. The epicentral distance ranges from 10 km to 60 km. We computed synthetic seismograms using CVM3.0, AID and SoCaL structure models. Figure 3.2 displays the P wave velocity on a cross-section through CVM3.0, the path-averaged AID model along this profile and the laterally homogeneous SoCaL model. The 3D synthetics were computed using K Olsen’s fourth-order staggered-grid finite-difference code (Olsen, 1994). The ID synthetics were computed using frequency-wavenumber integration (Zhu and Rivera, 2003). In general, 3D Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 80 synthetics calculated using CVM3.0 provide better fit to the observed seismograms than ID synthetics computed using either laterally homogeneous or path-averaged ID models. We lowpass-filtered all seismograms using a 4th-order Butterworth filter with corner frequency at 1.0 Hz. Then we compared observed seismograms with their corresponding synthetics. If we found a certain well-isolated arrival on the complete synthetic seismogram whose shape is in good agreement with the corresponding arrival on the observed seismogram, we constructed an isolation filter by windowing the complete synthetic using a cosine taper. Arrivals that were not well isolated or in poor agreement with their observed opponents were not included in the isolation filters. We obtained 494 isolation filters including 165 P waves, 110 SV, 170 SH and 48 other phases (Table 2.1) in the final data set. This process of selecting isolation filters was done manually. The isolation filters were cross-correlated both with the observed seismograms and the completed synthetics and the resulting cross-correlagrams were then windowed and narrow-band filtered in 0.02-Hz bands at 0.2, 0.4, 0.6, 0.8 and 1.0 Hz. The phase-delay and amplitude-reduction times were obtained by fitting Gaussian wavelets to the narrow-band filtered cross-correlagrams. We corrected our phase-delay time measurements for cycle-skipping using the broadband cross correlation travel-time anomalies as calibrators. The final GSDF data set comprised 9955 dtv and 9955 ^m easurem ents. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. An observational error is assigned to each GSDF measurement based on its narrow-band filtering frequency, the signal to noise ratio of the seismograms, the separation of the target arrival from other arrivals and the agreement between synthetic and observed waveforms. For most of our GSDF measurements, the observational error is within 10% of the measurements. CVM3.0 AID SoCaL Mean (Stp) Total -0.3136 -0.1413 0.6844 P -0.2934 0.0829 0.6035 SV -0.2623 -0.3992 0.6803 SH -0.3408 -0.2241 0.7916 Median ( & p ) Total -0.3347 -0.2125 0.6316 P -0.3127 -0.0445 0.5153 SV -0.3355 -0.5126 0.6473 SH -0.3465 -0.3231 0.7358 V ariance (* P ) Total 0.37499 0.71727 0.94167 P 0.28950 0.51549 0.58469 SV 0.52026 1.03760 1.42430 SH 0.39309 0.83763 0.93968 Variance (&q) Total 0.16019 0.22608 0.22799 P 0.0942 0.2151 0.2627 SV 0.1636 0.3566 0.4327 SH 0.1162 0.2371 0.3554 IQR (Stp) Total 0.4835 0.9419 1.0329 P 0.3453 0.6790 0.8135 SV 0.9762 1.0778 1.4438 SH 0.5918 1.0655 1.0687 IQR (&q) Total 0.2370 0.3364 0.3913 P 0.2299 0.3266 0.3513 SV 0.3364 0.4717 0.5236 SH 0.2373 0.2843 0.3800 Table 3.1. L I- and L2-norm based statistics made on GSDF measurements. The L I- (median and IQR) and L2-norm (mean and variance) based statistics of GSDF measurements were made on observed waveforms using CVM3.0, A ID and SoCaL synthetics as references. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 82 3.4.1 Improvement of CVM 3.0 synthetics over A ID or SoCaL synthetics As summarized in Table 2.1 and Figure 3.7, the dtm measurements made using CVM3.0 synthetics have the smallest variances both for P and for S waves, suggesting that in general CVM3.0 synthetic waveforms provide a better fit to the observed waveforms than AID or SoCaL synthetics. The longer tail in SoCaL phase-delay measurements is probably due to the basin that is not accounted for in SoCaL. The variance reduction from SoCaL to CVM3.0 is 61% in dtp and 64% in dtq, the variance reduction from AID to CVM3.0 is 55% and 53%, respectively. And the variance reductions for measurements of S waves are as large as those for P waves, even though S velocity below 300 m depth in Los Angeles Basin area was inferred from P velocity using empirical rules based on Poisson’s ratio rather than S wave tomography (Magistrale et al., 2000). This result suggests that those empirical rules are good approximations. S wave 3tp measurements do have larger variances than P wave measurements, suggesting the S velocity model of CVM3.0 is not as good as its P velocity model. Statistics based on L2-norm such as “mean” and “variance” is sensitive to outliers. To give a more robust evaluation about the structure models, we also computed LI-norm based statistics, the median and the Interquartile Range (IQR), which is the range that contains the middle 50% of a distribution (Table 2.1 and Figure 3.8). The IQR for CVM3.0 Stv measurements is 46% less than SoCaL and 43% less than AID. For dtq measurements, the reduction in IQR is 35% relative to Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 83 SoCaL and 25% relative to AID. Since changes in the upper and lower 25% of the data do not affect it, the IQR is a more robust estimate of the spread of the data. 600 o rd > u 400 200 Var=0.37499 Var=0.16019 i 1000 1 1 500 0 1 400 Q < 200 ro U O o n 600 400 200 -5 St (s) Var=0.71727 1000 Var=0.22608 I JL . 500 0 1 - 5 Var=0.94167 1000 Var=0.22799 500 0 L St (s) Figure 3.7 (a) Distribution of GSDF measurements made on observed waveforms using CVM 3.0 (first row), A ID (second row) and SoCaL (third row) synthetics as references, (a) shows the phase-delay (left column) and amplitude-reduction times (right column) for the three types of synthetics. CVM3.0 has not only the smallest variances in < 5 f P jq , but also the smallest shift in Stp (about -0.3 s) relative to observed data (Table 2.1 and Figure 3.7). Both CVM3.0 and AID are too slow, while SoCaL is too fast. If we use SoCaL synthetic as the reference, we see a strong correlation between CVM3.0 synthetics Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 84 and observed waveform (Figure 3.9). The correlation coefficient is 0.87 for Stp measurements and 0.80 for Stq measurements. The correlation coefficients provide another way to measure the improvement of CVM3.0 synthetics over SoCaL synthetics. 200 Q < 1 0 0 Var=0.2895 Var=0.51549 Var=0.58469 Var=0.39309 5 -5 0 Var=0.83763 200 1 0 0 5 -5 200 Var=0.93168 1 0 0 A 5 -5 SH Var=0.52026 60 40 20 0 5 - 5 5 - 5 5 - 5 V ar=l .4243 Var=1.0376 40 20 0 Figure 3.7 (b) Distribution of GSDF measurements made on observed waveforms using CVM 3.0 (first row), A ID (second row) and SoCaL (third row) synthetics as references, (b) shows the phase-delay time for P, SH and SV waves. The numbers on top of each histogram are the variance of this distribution. The variance for Stq measurements is about two times smaller than the variance for Stp measurements, which is not the case for previous teleseismic studies of seismic structures (Gaherty and Jordan, 1996; Katzman, et al, 1998) and Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 85 earthquake sources (McGuire et al., 2001, 2002). As shown in Zhao et al. (2005) and Chen et al. (2005b), the amplitude-reduction time is less sensitive to velocity perturbations than phase-delay time in our local seismic studies on Los Angelese area. The major factors that have large effects on dtq measurements are the attenuation model and focal mechanisms. Because of short wave propagation distances and relatively high Q values in this study, the variation in dtq caused by variation in attenuation is small. And an examination of first-motion polarities of synthetic and observed seismograms shows the focal mechanisms determined by Hauksson (2000) are excellent. c r C L c £ > CVM AID SoCaL CVM AID SoCaL F igure 3.8 Box-plots for GSDF measurements. A standard box plot generally shows the median (red line in the middle of the box), upper hinge (75th percentile, upper boundary of the box), lower hinge (25th percentile, lower boundary o f the box), higher adjacent value (black line above the box, the largest value in the data below upper hinge plus 1.5 times IQR), lower adjacent value (black line below the box, the smallest value in the data above lower hinge minus 1.5 times IQR) and outliers (red crosses). The height of the box gives the IQR. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 86 as - 0.5 Q . - 1.5 - 5 -2 5t'( DAT-SoCaL) - 5 §t (DAT-SoCaL) Figure 3.9 The correlation plots between GSDF measurements m ade on observed waveforms using SoCaL synthetics as references (horizontal axes) and GSDF measurements made on CVM3.0 synthetics using SoCaL synthetics as references (vertical axes). Left: phase-delay times; right: amplitude-reduction times. 3.4.2 Spatial variations of phase-delay anom aly Stp The mean of Stp measurements for CVM3.0 synthetics is about -0.3 seconds (Table 2.1 and Figure 3.7). CVM3.0 might be too slow in general, or it might be due to errors in earthquake origin times. To resolve the phase-delay anomaly Stp as spatial features, we used the differential phase-delay anomaly <S*p C V M SoC aL = r)Zp(SoCaL) - d/p(CVM3.0) as a measure of sedimentary basin thickness and we ploted dtv measurements for CVM3.0 synthetics as a function of & 5fp C V M ~ S o C aL (Figure 3.10). We observed a negative correlation between c ) 7 p(CVM3.0) and < 5 ( 5 r p C V M ~ S o C aL , the slope of a least- Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. square fit to the measurements is about -0.2, suggesting the thicker the basin, the slower CVM3.0 is. p Figure 3.10 Phase-delay times (circles) measured on observed waveforms using CVM3.0 synthetics as references (vertical axis) plotted as a function of the differential phase-delay times between SoCaL and CVM 3.0 measurements (horizontal axis), which are used as a proxy for sedimentary basin thickness. The straight line shows a least- square fit to the measurements. An example is given on Figure 3.11. The seismograms shown are from the 09 Sept. 2002 Yorba Linda earthquake at station DLA. The sensitivity kernel (Zhao et al., 2005) shows that the contribution to the phase-delay anomaly mostly comes from the basin underneath the station. The measured phase-delay anomaly 3tp of about-0.55 s suggests CVM3.0 might be too slow in the basin area. 3.4.3 Searching for anisotropy CVM3.0 is an isotropic model. To look for possible anisotropy in the crust, we compared our 8tv measurements made on 189 SV-SH pairs. Figure 3.12 shows the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 88 correlation between the two Stp measurements made at 0.2, 0.4, 0.6, 0.8 and 1.0 Hz. In spite of strong scattering, the measurements are slightly biased towards Stp(SV) > r)'/p(SH) side. To demonstrate the significance of this bias, we carried out a Binomial test using the number of data points lying in the c)/p(SV) > (5 /p(SH) region as test statistic. The null hypothesis that there is no bias towards Stp(SV) > 4fp(SH) side can be rejected at 50% or larger confidence level (Figure 3.12). a kernel for 5rp , source-receiver vertical plane • |: hs = 8.48 km W-E component, m odified SCEC C V M v3.0 Figure 3.11 Left: comparison of W -E component record (top) with finite-difference synthetics in modified SCEC CVM 3.0 with minimum P velocity of 1.25 km/s (middle) and 2.5 km/s (bottom). Right: the sensitivity kernel (Zhao et. al., 2004) shows that the contribution to phase-delay anomaly comes mostly from the basin underneath the station. The measured phase-delay anomaly of about -0.55 s suggests that (modified) CVM3.0 might be too slow in the basin (red color indicates and advance in arrival time for velocity increase). To give a quantitative estimate of the amount of anisotropy, we computed the differential phase-delay anomaly c )c )/p S V ~ S H = [d/p(SV) - rVp(SH) ] / epicenter distance. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 89 The mean and median of this differential phase-delay anomaly (Figure 3.13) give an estimate of about 0.5-4 ms/km. A more detailed analysis of Sdt^ SV-SH shows CVM-SoCaL that it has a positive correlation with the differential phase-delay SStp (Figure 3.14), which we used as a measure of basin thickness. This slight polarization anisotropy as detected in our GSDF data analysis might be due to basin layering effect. Comparing with phase-delay anomaly caused by 3D seismic velocity variations, it is not necessary to include polarization anisotropy in our future tomography inversions. 0.2 Hz 0.4 Hz 0.6 Hz OUJD O c - 4 - 2 0 2 St(SH) (s) P 0.8 Hz -2 -3 -3 - 2 -1 0 1.0 Hz SV slow & -2 SH slow -3 -3 -2 1 0 1 po -2 -2 1 0 0.2 0.4 0.6 0.8 frequency (Hz) Figure 3.12 Correlation between phase-delay anomalies of SV waves < S tp(SV) and corresponding SH waves c 5 fp(SH) at 0.2, 0.4, 0.6 and 0.8 Hz. There are more points lying in the Vsv < VS H region (above the straight line) than in the region Vsv > VS H (below the straight line). The lower-right plot shows a binomial test of the significance o f this bias. The four horizontal straight lines show (from bottom to top) 50%, 80%, 90% and 95% confidence level at which we accept the hypothesis that this bias is significant. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 90 x 10 ’ 4.5 — m edian 3.5 E in < U U c T O 2.5 X o n o . 4— < C O > u " > o . c S 0.5 0.2 0.3 0.4 0.5 0.6 0.7 0.9 0.8 frequency (Hz) F ig u re 3.13 The mean and median of ddtvsw'su as a function of frequency. The amount o f polarization anisotropy is on the order of 0.5~4ms/km. - 2 —3 - 4 - 5 - 4 — 3 - 2 - 1 0 St (CVM )- 5 t (SoCaL) £s) 1 2 3 4 5 (SoCaL) Figure 3.14 The differential phase-delay tim e between SV and corresponding SH waves, which we use to measure the amount of anisotropy, plotted as a function of the differential phase-delay tim e between CVM3.0 and SoCaL synthetics, which are used as proxy for sedimentary basin thickness. The straight line shows a least-square fit to all the measurements. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 91 3.4.4 Attenuation Effect So far, we have not included attenuation in our 3D synthetics. Our measurements of amplitude-reduction time can help to evaluate distance-attenuation scaling relationships. 4 3 2 1 0 -1 -2 -3 0 20 40 60 80 100 120 140 distance (km) F igure 3.15. The amplitude-reduction time measurements as a function of hypocenter distance. M easurements at all five frequencies are included. Figure 3.15 shows our amplitude-reduction time measurements as a function of hypocenter distance. A linear least-square regression to all observations shows a positive slope of about 0.0010, which means the amplitude ratio between observed waveform and synthetic waveform decreases with hypocenter distance at a very slow rate within our frequency band. We did not consider attenuation when we computed our 3D synthetics. This analysis on amplitude-reduction time Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 92 shows that in general attenuation has only a very small effect for the source- station paths under study. The slope of the least-squares fit of about 0.001 implies an attenuate factor Q of about 300 for S waves. 3.4.5 Comparison between SCEC CVM 3.0 and Harvard3D Harvard3D is a 3D structure model of Southern California developed by Suss and Shaw (2003). Its P velocity in the Los Angeles basin area is constrained by sonic logs and reflection profiles from oil industry. Its S velocity is scaled from P velocity using empirical relations. To compare SCEC CVM3.0 with Harvard3D, we measured the same set of observed waveforms, which were used to make SCEC CVM3.0 measurements in previous sections, relative to synthetics computed using Harvard3D model. Figure 3.16 shows the histograms of Harvard3D measurements together with measurements for SCEC CVM3.0 and the two types of ID structure models. The measurements for Harvard3D have less scattering than the measurements for SCEC CVM3.0, suggesting Harvard3D provides a slightly better fit to the observed waveforms than SCEC CVM3.0. A more detailed analysis on phase-delay time measurements (Figure 3.17) shows that Harvard3D provides better fit to P wave measurements than SCEC CVM3.0, while SCEC CVM3.0 probably has a slightly better S velocity model. We also found a high correlation between SCEC CVM3.0 measurements and Harvard3D measurements. The correlation coefficient is about 0.6 for phase-delay Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 93 measurements and about 0.8 for amplitude-reduction time measurements (Figure 3.18). Such high correlation between GSDF measurements might suggest high similarity between the two 3D models. mean=-0.28012 std=0.61405 o o "P i I m ean= -032332 std=0.6615 > 1 0 0 U mean=-0.14132 std=G.94635 so 60 O < 40 mean=0.68441 std=1.0139 1 0 0 80 — 1 60 20 -5 mean=— 0.0043617 std=0.34345 250 200 150 1 0 0 50 m ean=-0.0071119 std=0.34952 200 100 50 - 5 m ean=-0.05924 std=0.51231 2 0 0 150 too 50 m ean=-0.098856 std=0.58182 200 150 100 50 -s Stp (s) Stn (S) F igure 3.16. Histograms showing the fit of the models to the phase-delay measurements (left column) and amplitude-reduction time measurements (right column). Each histogram comprises of about 5000 measurements on P, SV and SH waves at 0.2, 0.4, 0.6, 0.8 and 1.0 Hz for source station paths shown in Figure 3.1. The mean and the standard deviation for each histogram are shown on the top of each plot. The GSDF residuals were measured using the synthetic waveforms for each model as isolation filters. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 94 P 200 1.50 ro 100 so -5 0.3718: 0.20879 0.3578: 150 100 Q < .6604J -s 0 5 ISO 100 -5 5t (s) 100 40 1 2 0 100 80 60 40 20 8t (s) Figure 3.17. Histograms showing the phase-delay measurements for P (left column) and SH (right column) waves. The numbers shown are variances for each histogram. For P wave measurements, Harvard3D is slightly better than CVM3.0. For SH wave measurements, CVM3.0 is slightly better than Harvard3D. For both P and SH waves, the two 3D models have much smaller variances than either path-averaged ID model (AID) or laterally homogeneous ID model (SoCaL). Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 95 Corrcoef - 0.6 Corrcoef = 0.8 Havard3D - 2 - 3 -2 CVM3.0 1.5 0.5 -05 -1.5 -2 0 1 2 CVM3.0 Figure 3.18. The correlation between CVM3.0 and Havard3D phase-delay measurements (left) and amplitude-reduction tim e measurements (right). The correlation coefficient is 0.6 for phase-delay measurem ent and 0.8 for amplitude-reduction time measurements. 3.5 Conclusions GSDF measurements of travel time and amplitude anomalies show that in general SCEC CVM3.0 and Harvard3D provide much better results than AID or SoCaL in predicting observed waveforms. Both 3D models can be used as reference models in linearized tomography inversions. These observations of phase and amplitude anomalies StP A also facilitate us to discover possible problems within our structure model. To resolve the velocity discrepancy as a spatial feature and the existence of anisotropy in the crust, we intend to invert these GSDF measurements for structural perturbations. In particular, Zhao, et al. (2005) has developed codes that can efficiently compute Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 96 the Frechet kernels for the GSDF measurements using CVM3.0 as a starting model. In fact, we have shown how the same measurements can also be employed to invert for seismic source properties, including the finite moment tensor (FMT) (Chen et al., 2005a). The GSDF methodology thus provides a unified framework for regional studies of seismic sources and Earth structure in Southern California and elsewhere. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 97 3.6 References Chen, P., L. Zhao and T. H. Jordan, Finite Moment Tensor of the 3 September 2002 Yorba Linda Earthquake. Bull. Seismol. Soc. Am. 2005a, in press. Chen, P., L. Zhao and T. H. Jordan, Generalized Seismological Data Functionals for Broadband Seismic Waveform Analysis and Inversion. Geophys. J. Int. 2005b, in preparation. Dreger, D. S. and D. V. Helmberger, Determination of Source Parameters at Regional Distances With Three-component Sparse Network Data. J. Geophys. Res. 98, 8107-8125, 1993. Gaherty, J. B., T. H. Jordan and L. S. Gee, Seismic Structure of the Upper Mantle in a Central Pacific Corridor, J. Geophys. Res. 101, 22291-22309. Gee, L. S. and T. H. Jordan, Generalized Seismological Data Functionals. Geophys. J. Int. Ill, 363-390, 1992. Graves, R. W. and D. J. Wald, Resolution Analysis of Finite Fault Source Inversion Using ID and 3D Green’s Functions, Part I: Strong Motions. J. Geophys. Res. 106, 8745-8766. Hauksson, E., Crustal structure and seismicity distribution adjacent to the Pacific and North America plate boundary in southern California. J. Geophys. Res. 105, 13875-13903, 2000. Katzman, R., L. Zhao and T. H. Jordan, High-resolution Two-dimensional Vertical Tomography of the Central Pacific Mantle Using ScS Reverberations and Frequency-dependent Travel Times. J. Geophys. Res. 103, 17933-17971. Magistrate, H., S. Day, R. W. Clayton and R. Graves, The SCEC Southern California Reference Three-dimensional Seismic Velocity Model Version 2. Bull. Seismol. Soc. Am. 90, S65-S76. McGuire, J. J., L. Zhao, and T. H. Jordan, Teleseismic inversion for the second- degree moments of earthquakes. Geophys. J. Int. 145, 661-678, 2001. McGuire, J. J., L. Zhao, and T. H. Jordan, Predominance of unilateral rupture for a global catalog of earthquakes. Bull. Seismol. Soc. Am. 92, 3309-3317, 2002. Olsen, K. B., Simulation of Three-dimensional Wave Propagation in the Salt Lake Basin. Ph.D. Thesis, University of Utah, Salt Lake City, Utah, 157p, 1994. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 98 Tromp, J., C. Tape and Q. Liu, Seismic Tomography, Adjoint Methods, Time Reversal and Banana-Doughnut Kernels. Geophys. J. Int. 2004, accepted. Zhao, L., T. H. Jordan, K. B. Olsen and P. Chen, 2005. Frechet Kernels for Imaging regional earth structure based on three-dimensional reference models. Bull. Seismol. Soc. Am, submitted. Zhu, L. and D. Helmberger, Advancement in source estimation techniques using broadband regional seismograms, Bull. Seismol. Soc. Am., 86, 1634-1641, 1996. Zhu, L., and L. A. Rivera, A note on the dynamic and static displacements from a point source in multilayered media Geophys. J. Int. 148, 619-627, 2003. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 99 Chapter 4 Finite Moment Tensor of the 3 September 2002 Yorba Linda Earthquake 4.1 Abstract We have developed procedures for inverting broadband waveforms for the finite moment tensors (FMTs) of regional earthquakes. The FMT is defined in terms of second-order polynomial moments of the source space-time function; it removes the fault-plane ambiguity of the centroid moment tensor (CMT) and yields several additional parameters of seismological interest: the characteristic length Lc, width W c, and duration Tc of the faulting, as well as the directivity vector V d of the fault slip. Following McGuire et. al. (2001), we represent the observed waveform relative to the synthetic in terms of two frequency-dependent differential times, a phase delay Sz^(co) and an amplitude reduction time STq(co), which we measure using the GSDF method (Gee and Jordan, 1992, Chapter 2). We numerically calculate the FMT partial derivatives, which allows us to use synthetics computed using any forward modeling tools. We have tested our methodology on Southern California Seismic Network (SCSN) recordings of the 03 Sept 02 Yorba Linda earthquake (Mw = 4.3). Using ID synthetic Green functions, we determined the CMT and resolved fault plane ambiguity. To resolve the details of source Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 100 finiteness, we employed a joint inversion technique that recovers the CMT parameters of the aftershocks, as well as the CMT and FMT parameters of the mainshock. The joint system of equations relating the data to the source parameters of the mainshock-aftershock cluster is denuisanced for path anomalies in both observables; this projection operation effectively corrects the mainshock data for path-related anomalies in a way similar to, but more flexible than, empirical Green function (EGF) techniques. Our results indicate that the Yorba Linda rupture occurred as left-lateral slip on a fault patch conjugate to the nearby, right-lateral Whittier fault. We obtained source dimensions of Lc = 0.7 ±0.1 km, W c = 0.4 ±0.1 km, and Tc - 0.2 ± 0.05 s, implying a stress drop of about 3.2 MPa, and we found a directivity of 0.8 ± 0.2, oriented up and to the northeast. The inferred fault plane is consistent with the mainshock-aftershock distribution relocated by Hauksson (2002). 4.2 Introduction Seismic hazard analysis depends on seismology to quantify earthquake sources. Although any indigenous seismic source can be exactly represented by a space time distribution of stress glut T(r, t) (Backus and Mulcahy, 1976), this general representation is rarely used. Instead, seismologists approximate T(r, t) in various ways, depending on the observational circumstances. Three representations are commonly employed in Southern California, which is our primary area of interest (Figure 4.1): Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 101 1. Isotropic point source (IPS), specified by an origin time to, a hypocenter ro, and local magnitude ML . The Southern California Seismic Network (SCSN) routinely recovers and catalogs the IPS representation for events as small as M l ~ 1.0-1.5 (Mori et al., 1998). 2. Centroid moment tensor (CMT), specified by a symmetric seismic moment tensor M, a centroid time t\, and a centroid location ri (Dziewonski et al., 1981). The 10-parameter CMT representation is the most general indigenous point source, although supplementary restrictions are often applied to reduce the parameter set; e.g., only 8 parameters are needed to describe a pure double-couple. Source mechanisms are regularly determined from SCSN first motions down to M l ~ 2.0-2.5 (Hauksson, 2000), and complete CMT solutions can be automatically recovered from broadband waveform data for M l above 3.0-3.5 (Zhu and Helmberger, 1996; Pasyanos et al., 1996). 3. Fault slip distribution (FSD), specified by a displacement discontinuity Au(r, t) across a predetermined fault plane r e S. Depending on the parameterization and kinematic assumptions, the number of unknown parameters can easily exceed 100. Full FSD imaging is usually applied only to large (ML > 6) earthquakes, such as 1992 Landers earthquake (Wald and Heaton, 1994; Cohee and Beroza, 1994), 1994 Northridge earthquake (Wald et al., 1996), and 1999 Hector Mine earthquake (Dreger and Kaverina, 2000). Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 102 In this paper, we discuss new procedures for determining source parameters using a “finite moment tensor” (FMT)—the lowest order representation of a finite fault. The FMT resolves fault plane ambiguity of the CMT and provides several additional parameters of seismological interest, including the characteristic length Lc, width Wc, and duration Tc of the faulting, as well as the directivity vector V d of the fault slip. We expect the FMT representation to be most useful for earthquakes in the magnitude range M l ~ 4-6, where finite-source effects can be resolved from regional broadband data but are too small to warrant a full FSD analysis (Figure 4.1). Number of parameters 1 M agnitude Isotropic Centroid Finite Fault Point Moment Moment Slip Source Tensor Tensor Distribution (IPS) (CMT) (FMT) (FSD) 1 --• --- 1 +« 1 -> (5) 1 V (8-10) 1 (14-20) 1 1 (>100) 1 1 2 3 4 5 6 7 Figure 4.1. Four source representations discussed in the text, arrayed near the lower end of the magnitude range of their applicability in Southern California. W e have developed the methodology for recovering finite moment tensor (FMT) representations of earthquakes in the magnitude range 4-6 from regional broadband data. The number o f FM T parameters varies from 14 to 20, depending on the constraints on the source model. We have successfully applied our method to recover the FMT for 3 September 2002 Yorba Linda earthquake (Ml 4.8) using high-quality waveform data from the TriNet seismic network (Hauksson et al., 2001). The 24 broadband stations Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 103 used in this study are shown on Figure 4.2, and some examples of recorded and synthetic seismograms are shown on Figure 4.3. 241 30' 24200' 24230' 24300' 243 30' 3430' 34 30' TA2 D E C J J P X PDU 3400' 34 00' [N M L S LGB O LI B B S SRN STG 33 30' 33 30’ 243 30' 24300' 24200' 242 30' Figure 4.2. The beachball shows the location and focal mechanism of the Yorba Linda earthquake determined by us using ID synthetic Green functions. The 24 TriNet stations used in this study are shown as triangles. Major faults in this region are shown as solid lines. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 104 Vertical(EGF1) North(EGF1) SRN A ^ a a x a/V wv^a /\/\/v /\/ 0 . 0 0 4 5 8 5 7 D . 0Q1317- CHN 0 . 0 0 1 6 8 2 5 CHN 0 . 0 0 7 2 1 9 5 CPP 0 . 0 0 2 2 3 5 2 . 5Z 3e -05 0 ■ 0 0 9 2 7 4 8 OLI WLT seconds seconds E a s t (EGFl) V e r t i c a l (EGF2) 0 . 0 0 0 4 3 4 3 4 R .3 9 0 8 ^ 0 5 a |V \A /'A /A aJ v a ^ ~ - — 0 . 2 2 2 5 1 0 ...00746.8.1 CHN 0 . 3 4 0 7 1 y ^ /w \A A A \^ A /y ^ i V V v V V u ^ /'A ''^“ V \/~ '" 0 . 0 0 0 2 3 8 8 0 . 1 9 4 8 3 0 . 0 0 9 2 2 2 1 CPP 0 . 3 2 6 2 2 CPP n . n n f i i R 0 7 0 . 0 0 0 4 4 7 1 8 4 . 9 8 3 8 e 0 5 0 . 3 0 1 8 1 ■ /yV w \/vv-w \y\i |^-v— aW v 'V W A 0 . 0 0 1 8 3 5 8 WLT 0 . 0 0 3 1 5 8 8 OLI 0 • 5 54 6 5 —'/ v ^ V V * 0 . 0 0 0 9 5 6 7 2 5 . 7 7 9 2 e 0 5 0 . 2 2 0 1 5 P..Q Q '2AS51 seconds seconds F ig u re 4.3 Some examples of recorded and synthetic seismograms show P waves on vertical components, S waves on east and north components. Seismograms on “Vertical (E G F l)”, “North (E G F l)” and “East (E G F l)” are band-pass filtered from 0.5 Hz to 2.5 Hz; seismograms on “Vertical (EGF2)” are bandpassed from 1.5 Hz to 5 Hz. Upper trace: mainshock; middle trace: one of the aftershocks used in our study; lower trace: synthetic Green function. Numbers show the peak amplitudes in cm/s. In our FMT inversions, we assume the aftershocks have the same focal mechanism as the mainshock, consistent with the data shown here. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 105 4.3 The Finite Moment Tensor Our formulation of the FMT representation follows McGuire et al. (2001) in assuming the stress glut is everywhere proportional to a constant seismic moment tensor, r(r,f) = M ?f(T,t')dt'. (4.1) ‘0 The excitation of seismic waves is thus proportional to the source space-time function f(r,t), which we assume is zero before an origin time to- We expand f(r,t) in terms of its space-time polynomial moments \i< m J > (Backus and Mulcahy, 1976) and approximate it by the low-order terms. The zeroth-order moment is unity by definition, //(0 ,0 ) = \\ f(r,t)dVdt = 1, and the first-order moments provide the spatial and temporal centroids: ^ (1 0 ) = ^f{Y ,t)rd V d t= rx, (4.2a) y 0-1 ) = \\f(v,t)tdVdt = tx. (4.2b) The tensor M, vector ri, and scalar t\ form the CMT (point-source) representation. The lowest order description of source finiteness is given by the second moments. The second central moment in time specifies the characteristic duration of the event (e.g., Silver and Jordan, 1983), Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 106 ju(0 '2 ) = \ \ f ( r , m - t x)2d,Vdt=(Tc IT)2 (4.3) The second central moment in space is a symmetric second-order tensor, £ (2 ,o > = ^ f t j j X r - r . X r - r ^ d V d t = UAUT . (4.4) Here A is a diagonal matrix of eigenvalues, and U is an orthogonal matrix of eigenvectors. For the special case of a planar rectangular dislocation of length L and width W, the eigenvalue matrix reduces to A=diag[(Lc / 2)2, (Wc / 2)2,0], (4.5) where Lc = L /V 3 is the characteristic length and Wc = W / a/3 is the characteristic width of the source. The mixed space-time moment gives the directivity vector (Ben-Menahem, 1961), p (U ) = j ] / ( r ,0 ( r -! ))(; - t , ) d V d t ^ 2 ) v d. (4.6) The directivity velocity vector V d will lie in the plane of a simple dislocation. The magnitude of this vector is zero for a symmetric bilateral rupture and LJTC for a unilateral rupture along the length of the fault. Hence, the ratio d = ||vd|| TJLC measures the magnitude of the directivity as a dimensionless number between zero and unity (McGuire et al., 2002). Under the strictures of equation (4.1), the FMT representation comprises 20 parameters, 10 for the CMT plus 10 for the second central moments of f(r,t). For Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 107 most small and many moderate-sized earthquakes, a planar dislocation may provide an adequate source model, in which case the total number of FMT parameters reduces to 14: eight CMT parameters plus Lc , Wc, Tc, d, one principal- axis orientation parameter in U, and one vector orientation parameter in va. To recover FMTs for moderate-sized earthquakes using regional networks, we must correct the waveforms for propagation effects, which in Southern California can be strongly three-dimensional. Most previous researchers have accounted for the inadequacy of available Earth models by using empirical Green functions (EGFs). One EGF method is to deconvolve the propagation effects using the waveforms from proximate events of negligible spatial and temporal extent, such as small aftershocks, and then invert the pulses resulting from deconvolution for a slip distribution (e.g., Mori and Hartzell, 1990; Mori, 1996; Courboulex et. al., 1998). McGuire’s (2004) recent study of small earthquakes shows that this type of EGF analysis can also give robust estimates of FMT parameters. A more direct EGF method is to match the mainshock waveforms by superposing EGFs constrained by kinematic rupture models, thus avoiding deconvolution (Hellweg and Boatwright, 1999; Okada et. al., 2001). Our objective is to develop a methodology for FMT inversion that does not rely on deconvolution and can take full advantage of model-based Green functions— synthetic seismograms—as well as EGFs. In this preliminary study, we have used synthetics from ID structural models, but our formulation can be easily adapted to Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 108 3D structures now available for Southern California, such as the SCEC Community Velocity Model (Kohler et al., 2003; Magistrale et al, 2000; Hauksson, 2000) and the Harvard Model (Suss & Shaw, 2003; Komatitsch et. al., 2004). Our procedure requires the computation of synthetic seismograms and partial derivatives to set up the inversion, which can be avoided in traditional EGF analysis, but it pays off by allowing us to incorporate more information about Earth structure and the distribution of seismic sources in the determination of the FMT. Moreover, it provides a uniform methodology to invert waveforms for both Earth and source structure. 4.4 GSDF Analysis In the frequency domain, an observed waveform u(tr, a S ) at a receiver location Tr can be related to a point-source synthetic waveform u (rR, (0) by a differential operator expressed in the exponential (Rytov) form: /• \ ~ / \ \( o S T n { c o ) - c o S r tA c o ) _ N n(rR , C O ) = u (rR , cd) e (4.7) where 8 tv( o )) is a frequency-dependent phase delay time, and 8Tq{co) is an “amplitude reduction” time, so named because a positive value of 8 r q reduces the amplitude of u relative to u . McGuire et al. (2001) have linearized the forward problem of relating the functionals 8TpA to the FMT parameters to obtain: Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 109 (4.8b) Here Vs is the gradient operator with respect to the source coordinates, and the second-order spatial derivatives are We computed the first and second-order spatial derivatives by applying a symmetrized finite-difference formula (Press et. al., 2001) to point-source synthetic seismograms, which were calculated for one-dimensional (ID), multi layered velocity models using Zhu and Rivera's (2003) frequency-wavenumber method. Examples of the partial derivatives of with respect to the source parameters (sensitivity coefficients) are shown for P waves from the Yorba Linda earthquake in Figure 4.4. Some intuitive aspects can be seen; e.g., the sensitivity coefficients for directivity parameters increase with frequency and are negative in the direction of rupture propagation, because a negative Srq implies an amplitude increase. The kernel for second-order spatial moment has a cos 2 9 dependence on azimuth because of constructive and destructive interference from waves generated by a finite source. For more examples of sensitivity coefficients, we refer to McGuire et. al. (2001). (4.9a) (4.9b) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 110 The phase delay times and amplitude reduction times were measured using Gee and Jordan's (1992) method of generalized seismological data functionals (GSDF), as illustrated in Figure 4.5. In the GSDF analysis, the targeted wave group— a body wave, surface wave or any other phase—is computed for a point- source model, and this “isolation filter” is cross-correlated both with the observed seismogram and with the full synthetic seismogram. The resulting cross- correlagrams are windowed in time domain in the vicinity of their peaks, and narrow-band filtered at a target frequency 0 % . The resulting waveforms can be approximated by Gaussian wavelets, and the parameters of these wavelets can be used to recover estimates of SzpiC D i) and < 3 Y q( & > ,■ ) corrected for filtering bias and interference effects (Gee and Jordan, 1992; Chen et al., 2003). The GSDF methodology provides a unified approach to the analysis and interpretation of the information on seismograms. It can be applied to any waveform on any component and, through the use of synthetic seismograms, can incorporate existing knowledge of source and Earth structure. The processing algorithm can use synthetic waveforms from ID or 3D models as isolation filters and accounts for waveform interference, P-SV coupling, and other complications modeled by the synthetics. Moreover, the algorithm automatically encodes the metadata needed for computing the Frechet kernels that specify the sensitivity of the phase and amplitude data to lateral as well as radial Earth structure (Zhao & Jordan, 1998; Zhao, Jordan & Chapman, 2000). The regional waveform data Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. I l l measured by this technique can therefore be used in joint inversions for source and Earth structure, including the evaluation and improvement of 3D regional Earth models (Chen et al., 2003). A ( 1,0) i a sx /sna o ) q *z * asya^™ a s a y a ^ 'i ) + astW 1 '1 * q 0.3 0.2 05 A* A a A A ^ a-A- 0.2 Hz 0 * A - 0.1 -0.5 - 0.2 180 360 360 180 3 2 1 -1 -2 -3 * a a s t ^ + 3 & y 9 n : am x y am * f * • * + A -* • & * A A + A A A i - - A - A - - - ty + a a A A ! -1 - T 180 360 0.06 0.3 0.04 0.2 0.02 0.1 2.0 Hz - 0.02 - 0.1 -0.04 0.2 -0.06 180 360 180 360 180 360 * A A A A A A 4 azimuth azimuth azimuth Figure 4.4 Plots of the sensitivity coefficients for Stq as a function of azimuth for P waves from the Yorba Linda earthquake. First column: coefficients for a centroid shift (triangles: derivative relative to a perturbation towards north; crosses: a perturbation upward). Second column: coefficients for the directivity moment p (1 ,1 ) (triangles: rupture propagates towards north; stars: rupture propagates towards east; crosses: rupture propagates upward). Third column: coefficients for second-order spatial mom ent p,< 2 ,0 ) (stars: (a*/2,0 ) triangles: ( txy< 2 '0); crosses: pK(2 '0)). Upper row: 0.2 Hz, lower row: 2.0 Hz. Dash lines indicate zero. The sensitivity of Stq with respect to directivity increases with frequency, while the sensitivity to centroid shift and second-order spatial mom ent decreases with frequency. The apparent scatter in the first column is caused by singularities at 30°, 120°, 210° and 300° azimuth, which are the nodal planes of P wave radiation pattern. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 112 main 0 5 shock ° data o s 3 2 synthetic , Green function 0 1 2 2 3 4 (b) 3 2 1 0 1 2 "2 3 4 6 t =4.597 t = 0.1 4 2 0 0 5 2 '1 0 1 x 10 6 t =4.600 t = 0.04 0 5 empirical2 Green 1 function o seconds t =4.583 t = 0.12 seconds seconds F igure 4.5. Examples o f GSDF processing. Observed and synthetic seismograms are shown on column (a); the isolation filter shown on column (b) was obtained by windowing the full synthetic. The cross-correlagrams shown on column (c) were windowed (dash lines) and narrow-band filtered. The resulting waveforms as shown on column (d) can be approximated by Gaussian wavelets parameterized by tq and tp. The differential measurements of StpA are then recovered from these parameters. 4.5 Single Event Inversion The linearized forward problem described by equation (8ab) can be written in the form A om o + Aimi + A2H I 2 = d. (4.10) Here d is a data vector of GSDF measurements of and Srq relative to the point-source synthetic seismograms, and the integer subscripts on the model Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 113 vectors mk denote the polynomial order of the source moment; i.e., the vector mo contains the perturbations to the zeroth moment (the moment tensor), mi specifies the centroid space and time shifts [p,(10), / / 0 ,1 )], m2 comprises the second-order moments [p(2’ 0), pn'lj, //° ’ 2 )], and the matrices Ak contain the partial derivatives(sensitivity coefficients) of the data with respect to the model parameters. We note that the relative magnitudes of the sensitivity coefficients for 8tv and 8rq alternate with the order number k; the Srq coefficients are generally larger for k even, whereas the Stp coefficients are larger for k odd. The geologic structure of the upper crust in Southern California is known to be very heterogeneous laterally as well as vertically. To account for the grosser aspects of the geographic heterogeneity, we extracted path-specific ID velocity structures from Version 3.0 of the 3D SCEC Community Velocity Model (Kohler et al., 2003; Magistrale et al, 2000; Hauksson, 2000) by horizontally averaging wave slownesses along source-receiver paths. We refer to these path-dependent ID structures collectively as the AID model. The data were inverted using a hierarchical, iterative method that began with observations at lower frequencies. We low-pass filtered both the recorded seismograms and the point-source synthetics from the AID model below 0.2 Hz and made GSDF measurements on 48 waveforms of P, S and surface waves from TriNet records with good azimuthal distribution (Figure 4.2). The process was Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 114 initiated by inverting the low-frequency measurements of the amplitude reduction times Srq for an improvement to the moment tensor M. We obtained a static moment of 3.5 x 101 5 Nm (Mw = 4.3) and the double-couple solution shown in Figure 4.2. The latter was in good agreement with the SCSN focal mechanism determined from first-motion data (Hauksson, 2002). Our best estimate of M had a small non-double-couple component, but a bootstrap test showed that the deviation from a pure double-couple was statistically insignificant. In the next step, we inverted the phase delays Stp for improvements to the source centroid. The spatial centroid shifted by only 0.2 km from the original SCSN location, which we again found to be statistically insignificant. Finally, we inverted the low-frequency amplitude data for the second-order moments to obtain the full FMT. The solution yielded a northeast-striking fault plane with a characteristic duration of about 0.1-0.2 s and a characteristic dimension of 1-2 km. The bootstrap test identified the northeast nodal plane as the fault plane with reasonably high (81%) confidence, but we found that the uncertainties in Lc were comparable to its magnitude and that the low-frequency data placed no useful constraints on the directivity vector V d. Using higher-frequency data did not improve the solution. Our results show that at low frequencies, where directivity caused by fi(1 ,1 ) has only a small effect, it is still possible to resolve fault-plane-ambiguity by careful analysis of the constructive and destructive interference effects caused by source Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 115 finiteness n(2 ,0 ), whose relative contribution to amplitude-reduction time Srq is larger at lower frequency (Figure 4.4). For a pure bilateral rupture with a strike of 30°, the amplitude of finite-source waveform relative to the point-source synthetic waveform should be maximized at azimuths of 120° and 300°, owing to constructive interference, while slip on the conjugate fault plane would generate a maximum at 30° and 210°. At frequencies below 0.2 Hz, the synthetic Green functions computed using path-specific ID structure models are good enough to detect such interference effects. In most previous studies, the fault planes of smaller earthquakes have been identified from directivity effects in the amplitudes and the widths of EGF- deconvolved source pulses (e.g., Mori and Hartzell, 1990; Mori, 1996). Our results are similar to McGuire (2004), who showed that the fault plane ambiguity can be resolved using EGF-derived estimates of /rs(0 ,2 ) even if the directivity is small. Our approach has the advantage of being applicable to events lacking suitable EGFs. 4.6 Multiple Event Inversion FMT estimation uncertainties for the single-event inversion are dominated by signal-generated noise. The finite-source perturbations to the GSDF observations for an event as small as Yorba Linda are less than about 0.3 s. Chen et al. (2003) have shown that, relative to AID synthetics, lateral heterogeneity in Southern Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 116 California introduces fluctuations in 5rq which have a root-mean-square (RMS) magnitude of about 0.45 s for P and S phases in the frequency band 0.2-1.2 Hz. This scatter can be reduced using SCEC CVM3.0 to generate the reference seismograms, but the RMS residuals are still on the order of 0.35 s (Chen et al., 2003), larger than the finite-source signal. To improve the signal-to-noise ratio, we formulated a joint inversion procedure that recovers the CMT parameters of the aftershocks as well as the FMT parameters of the mainshock. For each event in the cluster (j = 1 for the mainshock; j = 2, 3, ... for aftershocks), we write the linearized perturbation equation Ao m o' + A /m / + A /m / + By n = d7 . (4.11) The vector n is a set of path anomalies, which we treat as nuisance parameters. We assumed the dimensions of the cluster are small compared to the path lengths and therefore n is approximately the same for all events in a cluster. Since the aftershocks were small, we further assumed the effects of source finiteness were negligible; i.e., m / = 0 for ; = 2, 3, ... With these assumptions, the path anomalies can be effectively removed by applying a projection operator that annihilates the terms B; n in (11). Let A ^, B*, Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 117 and d be the concatenations of {A<7}, {B*7}, and {d; } for k = 0, 1,2 and j = 0, 1, 2, 3, .... Then, the appropriate annihilator is Qb = I - B (BT B)-1 BT, and the reduced system of equations can be written in the form Q b Aoiiio + Q b Aimi + Qb A2 m2 = Qb d. (4.12) This type of “denuisancing” procedure has been commonly used in relative- location techniques (e.g., Jordan and Sverdrup, 1981). When observations are available for all events from all stations in the network, the B; matrices are simply identity matrices; then, under our assumptions, Qb A2 = A2, and the projected data for the mainshock are equal to the original data, d°, corrected for path anomalies by subtracting the average value of the aftershock data. When some stations do not provide data for all events, the results are optimal in a least- squares sense (Jordan and Sverdrup, 1981). The multiple-event method is closely related to traditional EGF analysis, but it offers several advantages. No waveform deconvolution is involved, and it can utilize data from more than one aftershocks, including aftershocks recorded by only a subset of stations, allowing better averaging in the removal of the path anomalies. Moreover, because the perturbations are computed relative to synthetic seismograms, our methods can utilize the information in 3D seismic velocity models. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 118 4.7 Results More than 50 Yorba Linda aftershocks with local magnitudes ranging from 0.9 to 3.0 were located by the SCSN. For this preliminary multiple-event inversion, we selected a M l 2.9 aftershock that occurred 7 minutes after the mainshock (07:15:51.26 UT) and a M l 2.0 aftershock 4 minutes after the mainshock (07:12:55.20 UT); these were designated EGF1 and EGF2, respectively. We measured the amplitude reduction times Srq on 28 P waves and 28 S waves at frequencies ranging from 0.5 Hz to 2.5 Hz. Based on relocations by E. Hauksson (personal communication, 2002), the distances between mainshock and aftershock hypocenters are 0.2 km for EGF1 and 0.4 km for EGF2. A comparison of the three events across their common band of observation showed that Hauksson's relative location was good and that the moment tensors of both the mainshock and the two aftershocks were well- described by the same focal mechanism. We therefore approximated the aftershocks as point sources with our low-frequency mechanism (Figure 4.2). The differences in scalar moment between the mainshock and the two aftershocks can be denuisanced by normalizing mainshock and aftershock data to the same scalar moment. The amplitude-reduction times used in this study were insensitive to the relative distances between mainshock and aftershock centroid location, so we droped the first two terms in equation (4.12). The effects of denuisancing on the frequency-dependent data from two of the stations are shown in Figure 4.6. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 119 0 .4 5 0.4 0.35 0.3 0.25 0.2 O " 0.15 0.1 0.05 - 0.05 0.5 2.5 frequency (Hz) F ig u re 4.6 M easurements of Srq for P waves on the vertical component of station SRN (crosses) and MW C (circles). Solid line: raw observations; dash line: denuisanced for path effects using the aftershock data. To solve the inverse problem, we employed a semi-definite programming method (Vandenberghe and Boyd, 1996), which fits the data by least-squares while incorporating various nonlinear constraints, such as the positivity of the space time source ellipsoid (McGuire et al., 2001). The inversion results are given in Table 1 and Figure 4.8, and the fit to the observations is presented in Figure 4.7. With the source ellipsoid pi(2 '0 ) unconstrained in orientation, the inversion resolved the fault-plane ambiguity at a very high confidence level (> 99%). In the optimized solution, the source directivity points to the northeast and slightly upward, consistent with the source ellipsoid principal axes. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 120 0.5 Hz 1.0 Hz 1.5 Hz 2.0 Hz 25 Hz 0.3 0.3 0.2 0.2 0.2 0.2 L r , 0.1 O i -0.1 % 0 r 0.1 ! 0! -0.1 0.1 a : -0.1 0.1 O -0.1 m - - o -0.2 -0.2 -0.2 -0.2 200 2.0 Hz 200 0.5 Hz 200 1.0 Hz 2.5 Hz 1.5 Hz 200 1.8 Hz 2.0 Hz 2.2 Hz 0.3 0.3 0.2 02 a 2 0.1 $ 0 0.1 § O ai H o 0 0 l o 0.1 -0.1 -0.1 0 2 -0,2 - 0 2 Q measurements i -^•prediction j Azimuth (degree) F ig u re 4.7 GSDF measurements of amplitude reduction time as functions of azimuth and the fit of the inverted model. For EGF1, 19 P waves and 28 S waves were measured at frequencies ranging from 0.5 FIz to 2.5 Hz. For EGF2, 9 P waves were measured at frequency 1.8, 2.0 and 2.2 Hz. Circles: GSDF measurements; stars: predictions made using the inverted model and the sensitivity derivatives. Hence, the multiple-event inversion confirms that the Yorba Linda earthquake involved a left-lateral rupture on a plane conjugate to the nearby right-lateral Whittier fault (Figure 4.8). This conclusion is consistent with the aftershock distribution relocated by E. Hauksson (personal communication, 2002), which shows an extension perpendicular to the Whittier Fault (Figure 4.9). It is also consistent with recent studies showing that a high proportion of the north-south shortening in the Los Angeles metropolitan region is accommodated by conjugate Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 121 strike-slip faulting and east-west “escape” deformation (Walls et al, 1998; Yeats, 2003). Horizontal F igure 4.8. The inverted FM T parameters for the Yorba Linda event. Left: the map view, arrow shows the orientation of the fault plane and directivity; right: the characteristic dimensions o f the rupture, arrow shows the orientation of direcitivity. Whittier Fault (CF F igure 4.9. Projection o f a 3D rendering of the Yorba Linda main shock (blue sphere) and its aftershocks (gold sphere) relative to the W hittier Fault (white triangulated surfaces). The relocated hypocenters were provided by E. Hauksson, the CFM -A model by A. Plesch and the 3D rendering using the LA3D visualization system developed by the SCEC/ITR interns. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 122 To investigate the tradeoffs among various source parameters, we sampled the solution space using a bootstrap method. The characteristic duration Tc is well constrained at 0.1-0.2 s by the data, but we observe a strong tradeoff between the characteristic length Lc and the directivity d (Figure 4.10), which implies a large source (high Lc) has a low directivity, while a small source has a high directivity. Here we assume the rupture velocity, vR = L J \(2 - d)Tc\, which ranges from LC /2TC for a bilateral rupture (d = 0) to LJTC for a unilateral rupture {d = 1). Another parameter of concern is the stress drop. For Lc larger than 1.0 km, we obtain stress-drop below 10 bars (first column on Table 1). If we require rupture V =6.5 k m /s V =3.7 km/s 2 2 . 5 3 3 . 5 0 . 5 1 1 . 5 0 L (km) c Figure 4.10 The solution space o f the characteristic dimension Lc and directivity d as sampled using a bootstrap method. The triangle shows the unconstrained solution (first column of Table 1), the cross shows the solution obtained under the sub-shear constraint (second column of Table 1). The two straight lines show constant rupture velocity at P and S wave velocity in the seismic source region. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 123 velocity to be sub-shear (below shear wave velocity) in our inversion, we obtain a more compact source with a higher stress-drop of 32 bars and a directivity of 0.8 (second column on Table 1) and the fit to the observed data is equally good. Without constraints With constraints Lc (km) 1.47 0.74 Wc (km) 0.64 0.40 Tc (s) 0.15 0.19 d 0.39 0.80 Strike (°) 39.8 33.2 Dip (°) 83.2 88.3 Stress drop (bars) 4 32 Rupture velocity (km/s) 6.1 3.2 Table 4.1. Inverted FM T parameters for 09/03/02 Yorba Linda earthquake. First column: results obtained without any constraints except the positivity of source space time volumn. Second column: results obtained by constraining the rupture velocity below shear wave velocity in the inversion. Many, if not most, large earthquakes display high directivity (d > 0.5), presumably owing to the heterogeneity of major faults (McGuire et al., 2002). The results for the Yorba Linda earthquake contribute another datum in support of the notion that these symmetry-breaking processes also operate on a much smaller scale. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 124 4.8 Conclusions We have demonstrated that the FMT for a small earthquake can be recovered from regional waveform data. Our method is based on GSDF measurements of the differences between the observed and synthetic waveforms. Owing to the effects of 3D structure, we found it necessary to incorporate EGF-type constraints from aftershock data in the inversion of the differential amplitudes to resolve details of source finiteness. However, we note that even without using aftershock information, we were still able to determine the CMT solution and resolve fault- plane ambiguity using synthetic Green functions computed from path-specific ID structure models. We have already adapted our method to employ synthetics calculated from 3D models, such as the SCEC Community Velocity Model. FMT inversions using 3D synthetics for a set of earthquakes in the Los Angeles region will be presented in a subsequent paper. Our research suggests that improved versions of 3D models will eventually allow the automated recovery of FMTs for nearly all moderate-sized earthquakes in well-instrumented regions like Southern California. When applied to a large number of events, this new technique should improve our understanding of regional tectonics and earthquake scaling laws, as well as provide new data for seismic hazard analysis. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 125 4.9 Acknowledgments We thank Jeff McGuire for useful discussions and Egill Hauksson for his relocated aftershocks. Suggestions by Qinya Liu and an anonymous reviewer improved the paper. This research was supported by the Southern California Earthquake Center. SCEC is funded by NSF Cooperative Agreement EAR-0106924 and USGS Cooperative Agreement 02HQAG0008. The SCEC contribution number for this paper is 810. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 126 4.10 R eferen ces Backus, G., and M. Mulcahy, Moment tensors and other phenomenological descriptions of seismic sources, I, Continuous displacements. Geophys. J. Roy. Astron. Soc. 46, 341-361, 1976. Ben-Menahem, A., Radiation of seismic surface-waves from moving sources. Bull. Seismol. Soc. Am. 51, 403-435, 1961. Chen, P., L. Zhao and T. H. Jordan, A Quantitative Evaluation of SCEC Community Velocity Model Version 3.0, Eos, Trans. AGU, 84(46), Fall Meet. Suppl., Abstract S42H-02 (FI 108), 2003. Cohee, B., and G. Beroza, Slip distribution of the 1992 Landers earthquake and its implications for earthquake source mechanics. Bull. Seismol. Soc. Am. 84, 692- 712, 1994. Courboulex F., A. Deschamps, M. Cattaneo, F. Costi, J. Deverchere, J. Virieux, P. Augliera, V. Lanza, D. Spallarossa, Source study and tectonic implications of the 1995 Ventimiglia (border of Italy and France) earthquake (ML = 4.7). Tectonophysics 290, 245—257, 1998. Dreger, D., and A. Kaverina, Seismic remote sensing for the earthquake source process and near-source strong shaking: A case study of the October 16, 1999 Hector Mine earthquake. Geophys. Res. Lett. 27, 1941-1944, 2000. Dziewonski, A. M., T. A. Chou, and J. H. Woodhouse, Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. J. Geophys. Res. 86, 2825-2852, 1981. Gee, L. S., and T. H. Jordan, Generalized seismological data functionals. Geophys. J. Int. Ill, 363-390, 1992. Hadley, D., and H. Kanamori, Seismic structure of the Transverse Ranges, California. Geol. Soc. Am. Bull. 88, 1469 - 1478, 1977. Hauksson, E., Crustal structure and seismicity distribution adjacent to the Pacific and North America plate boundary in southern California. J. Geophys. Res. 105, 13875-13903, 2000. Hauksson, E., K. Hutton, L. Jones and D. Given, The M4.8 Yorba Linda, Orange County earthquake of 3 September 2002: California Integrated Seismic Network Preliminary Report, 2p, 2002. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 127 Hauksson, E., P. Small, K. Hafner, R. Busby, R. Clayton, J. Goltz, T. Heaton, K. Hutton, H. Kanamori, J. Polet, D. Given, L. M. Jones, and D. Wald. Caltech/USGS element of TriNet 1997-2001. Seismol. Res. Lett. 72, 690-704, 2001. Hellweg, M. and J. Boatwright, Mapping the rupture process of moderate earthquakes by inverting accelerograms, J. Geophys. Res. 104, 7319-7328, 1999. Jordan, T.H. and K.A. Sverdrup. Teleseismic location techniques and their application to earthquake clusters in the south-central Pacific, Bull. Seismol. Soc. Am. 71, 1105-1130, 1981. Komatitsch, D., Q. Liu, J. Tromp, P. Suss, C. Stidham and J. H. Shaw, Simulations of ground motion in the Los Angeles basin based upon spectral- element method, Bull. Seismol. Soc. Am. 94, 187-206, 2004 Kohler, M. D., H. Magistrale, R. W. Clayton, Mantle heterogeneities and the SCEC reference three-dimensional seismic velocity model version 3, Bull. Seismol. Soc. Am. 93, 757-774, 2003. Magistrale, H., S. Day, R. W. Clayton, and R. Graves, The SCEC southern California reference 3D seismic velocity model Version 2, Bull. Seismol. Soc. Am. 90, no. 6B, S65-S76, 2000. McGuire, J., Estimating finite source properties of small earthquake ruptures, Bull. Seismol. Soc. Am. 94, no 2, 377-393, 2004. McGuire, J. J., L. Zhao, and T. H. Jordan, Teleseismic inversion for the second- degree moments of earthquakes. Geophys. J. Int. 145, 661-678, 2001. McGuire, J. J., L. Zhao, and T. H. Jordan, Predominance of unilateral rupture for a global catalog of earthquakes. Bull. Seismol. Soc. Am. 92, 3309-3317, 2002. Mori, J., Rupture directivity and slip distribution of the M 4.3 foreshock to the 1992 Joshua Tree earthquake, Southern California. Bull. Seismol. Soc. Am. 86, 805-810, 1996. Mori, J. and S. Hartzell, Source inversion of the 1988 upland, California, earthquake: determination of a fault plane for a small event, Bull. Seismol. Soc. Am. 80, 507-518, 1990. Mori, J., H. Kanamori, J. Davis, E. Hauksson, R. Clayton, T. Heaton, L. Jones, A. Shakal, R. Porcella, Major improvements in progress for southern California earthquake monitoring, Eos 79, 217-221, 1998. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 128 Okada, T., N. Umino, Y. Ito, T. Matsuzawa, A. Hasegawa and M. Kamiyama, Source Processes of 15 September 1998 M5.0 Sendai, Northeastern Japan, Earthquake and Its M3.8 Foreshock by Waveform Inversion, Bull. Seismol. Soc. Am. 91, 1607-1618, 2001. Pasyanos, M. E., D. S. Dreger, and B. Romanowicz, Towards real-time determination of regional moment tensors. Bull. Seismol. Soc. Am. 86, 1255-1269, 1996. Press, W. H., S. A. Teukolsky, W. T. Vetterling and B. R. Flannery, Numerical Recipes in Fortran 77 - the art of scientific computing. Cambridge University Press, 2001. Silver, P., and T. H. Jordan, Total-moment spectra of fourteen large earthquakes. J. Geophys. Res. 88, 3273-3293, 1983. Suss, M. P. and J. H. Shaw, P-wave seismic velocity structure derived from sonic logs and industry reflection data in the Los Angeles basin, California, J. Geophys. Res., 108, NO. B3, 2170, ESE 13, 2003. Tukey, J. W., The collected works of John W. Tukey. Wadsworth, Belmont, CA., 1984. Vandenberghe, L. and S. Boyd, Semidefinite programming. SIAM Rev. 38, 49-95, 1996. Wald, D. J., and T. H. Heaton, Spatial and temporal distribution of slip for the 1992 Landers, California earthquake. Bull. Seismol. Soc. Am. 84, 668-691, 1994. Wald, D. J., T. H. Heaton and K. W. Hudnut, The slip history of the 1994 Northridge, California, earthquake determined from strong ground motion, teleseismic, GPS, and leveling data Bull. Seismol. Soc. Am. 86, S49-S70, 1996. Walls, C., T. Rockwell, K. Mueller, Y. Bock, S. Williams, J. Pfanner, J. Dolan and P. Fang, Escape tectonics in the Los Angeles metropolitan region and implications for seismic risk. Nature 394, 356—360, 1998. Yeats, R. S., Tectonics of the San Gabriel Basin and surroundings, southern California, GSA Bulletin 116(9/10): 1158-1182, 2004. Zhao, L. and T. H. Jordan, Sensitivity of frequency-dependent traveltimes to laterally heterogeneous, anisotropic Earth structure, Geophys. J. Int., 133, 683- 704, 1998. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 129 Zhao, L., T. H. Jordan and C. H. Chapman, Three-dimensional Frechet differential kernels for seismic delay times, Geophys. J. Int., 141, 558-576, 2000. Zhu, L. and D. Helmberger, Advancement in source estimation techniques using broadband regional seismograms, Bull. Seismol. Soc. Am., 86, 1634-1641, 1996. Zhu, L., and L. A. Rivera, A note on the dynamic and static displacements from a point source in multilayered media Geophys. J. Int. 148, 619-627, 2003. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 5 130 Resolving Fault-Plane-Ambiguity for Small Earthquakes 5.1 Abstract A central problem of seismology is the inversion of regional waveform data for models of earthquake sources. In regions such as Southern California, preliminary 3D earth structure models are already available, and efficient numerical methods have been developed for solving the point-source forward problem. We describe an automated procedure that utilizes these capabilities to derive centroid moment tensor (CMT) and to resolve the fault-plane-ambiguity associated with CMT solutions. The procedure relies on the use of receiver-side Green tensors (RGTs), which comprise the spatial-temporal displacements produced by the three orthogonal unit impulsive point forces acting at the receiver. We have constructed a RGT database for 64 broadband stations in the Los Angeles region using the SCEC Community Velocity Model (CVM) and K. Olsen’s finite-difference code. Three-dimensional synthetic seismograms for any earthquake can be simply calculated by extracting a small, source-centered volume from the RTG database and applying the reciprocity principle. The partial derivatives needed for the CMT inversion and for resolving fault-plane-ambiguity can be generated in the same way. Using this new technique, we have successfully determined CMT solutions for 149 small to moderate-sized earthquakes and resolved fault-plane-ambiguity Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 131 for 40 of these earthquakes. Our CMT solutions for earthquakes with local magnitude larger than 3.0 in and around Los Angeles basin area show 2 distinct earthquake clusters with similar double-couple focal mechanisms, one is in the Yorba Linda area, the other extends from Fontana area to Puente Hills. After removing fault-plane-ambiguity from our CMT solutions, the events in Yorba Linda area show left-lateral mechanism conjugate to the nearby right-lateral Whittier fault. Most of the events in the Fontana cluster show right-lateral mechanism that is parallel to the nearby San Jacinto fault but conjugate to the southwest trend of hypocenter distribution. Our results might indicate a developing weak zone associated with the “escaping” deformation in and around Los Angeles basin area. 5.2 Introduction The study of seismic sources has many important applications. Earthquake hypocenters have been used to delineate faults; earthquake focal mechanisms have been used to indicate relative plate motions and to study the orientation and state of the stress on subsurface faults; a better understanding of seismic sources can help us mitigate seismic hazard. The most general point source representation is the Centroid Moment Tensor (CMT), specified by a symmetric seismic moment tensor, a centroid time and a centroid location [Dziewonski et al., 1981]. For most tectonic earthquakes that involve only shear dislocation, the seismic moment tensor represents a double Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 132 couple source mechanism with two mathematically equivalent orthogonal nodal planes. The inability to distinguish the real fault plane on which slip occurred from the auxiliary plane, which is inherent in a double-couple source mechanism, is referred to as the “fault-plane-ambiguity” [Aki & Richards, 2002], Recent advancements both in seismic networks and in computing technology have made possible routine determination of focal mechanisms or complete CMT solutions. In Southern California, source mechanisms are regularly determined from first-motion data for earthquakes down to Ml2.0~2.5 [Hauksson, 2000], and complete CMT solutions can be automatically recovered from broadband waveform data (typical period range 5 to 100 seconds) for ML above 3 .0 -3 .5 [Zhu & Helmberger, 1996; Pasyanos et al., 1996; Liu et al., 2004], To provide better constraints on regional tectonics and to improve the assessment of stress-transfer models, which require knowledge of the fault plane to compute Coulomb stress increments, it is highly desirable to remove the fault-plane-ambiguity from our routinely determined source mechanisms. Resolving fault-plane-ambiguity will help to delineate complicated fault structures, especially when the faults are subsurface and the seismicity is diffusive. For example, Webb and Kanamori [1985] and Seeber and Armbruster [1995] proposed the existence of a regional decollement in the middle crust in the Transverse Ranges area based on their study of low-angle faulting or shallow- dipping focal mechanisms in that region. One nodal plane of some earthquakes at Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 133 the depth of 10 to 15 kilometers can dip at less than 25 degree. They thus suggested a structure model that a detachment fault floors the imbricate thrust system associated with the Big Bend of the San Andreas Fault. The seismicity in that area is diffusive, the floor of the seismicity could mark a detachment structure but could also be just the brittle-ductile transition zone, the problem is whether the shallow-dipping nodal plane is the fault plane. Middle-crust detachment structures were also suggested elsewhere, such as the Brevard zone of southern Appalachians [Edelman et al., 1987], the Vienna basin and Czechoslovakia [Royden, 1984], based on reflection profiles and other evidences. By removing the fault-plane ambiguity, we could help resolve whether such detachments really exist, find out the relative motion between the upper crust and the lower crust, and provide better constrains on regional tectonics. A generalized earthquake interaction model is the stress-transfer model, which explains not only the interactions between a mainshock and its aftershocks or foreshocks but also interactions between earthquakes separated by much larger spatial and temporal scales. The long-range correlated stress field capable of producing a large earthquake is established by stress-transfer of numerous small to intermediate earthquakes in a region scaling with the magnitude of the oncoming large earthquake [Bowman and King, 2001], Usually, we evaluate the stress-transfer model by computing the Coulomb stress increments [Harris, 1998]. For finite sources, the calculation of the Coulomb stress increments requires the knowledge of the fault plane. By routinely resolving fault-plane-ambiguity and Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 134 computing Coulomb stress increments, it is possible to detect and monitor the establishment of long-range correlation of the stress field in the crust and provide important clues for earthquake prediction. For large earthquakes, we can resolve the fault-plane-ambiguity by observing the spatial distribution of its aftershocks’ hypocenters. For small to moderate-sized earthquakes with few aftershocks, resolving fault-plane-ambiguity in a routine manner is a much more challenging task. In most previous studies, the fault planes of small to moderate-sized earthquakes have been identified from directivity effects using empirical Green function (EGF) techniques [e.g., Mori and Hartzell, 1990; Mori, 1996]. But for earthquakes lacking suitable EGFs or having no significant horizontal directivity, this type of techniques is less effective. In this study, we introduce a new automated procedure to determine CMT solutions and to resolve fault-plane-ambiguity for small earthquakes using 3D synthetic seismograms. Our procedure relies on the use of Receiver Green Tensors (RGTs) and numerically computed sensitivity derivatives for second- order Finite Moment Tensor (FMT), which is the lowest order representation of a finite fault [Chen et al., 2005, Chapter 4], This new technique does not depend on detecting directivity effects and is much more general than classical EGF methods. Using regional broadband waveform data collected from 64 stations of the California Integrated Seismic Network (CISN) (Figure 5.1), we have Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 135 successfully determined CMT solutions for 149 small earthquakes in Los Angeles basin area and resolved fault-plane-ambiguity for 40 of them. 5.3 Receiver Green Tensors and Reciprocity The Receiver-side Green Tensor (RGT) G;*(r, Tr; t) is the 1 th component wave- field at position r and time t for a kth component impulsive force at the receiver position T r . In regions such as Southern California, preliminary 3D earth structure models are already available, and efficient numerical methods have been developed for solving the point-source forward problem. Using Southern California Earthquake Center Community Velocity Model version 3.0 (CVM3.0) 35° 34° 33° 240° 241° 242° 243° Figure 5.1 Los Angeles basin and surrounding areas. Blue triangles show the locations of California Integrated Seismic Network (CISN) three-component broadband stations. Black solid lines show the surface traces of active faults in this area. Background color shows the thickness of sedimentary basin. The inner box is the area where we have computed Receiver Green Tensors for CISN stations. iP lp- EQW U jflR F tg ^ ' s | c - \ 'B B C ALP ( £ > 1 J-iL W G B A A SES t 0 -2000 -4000 -6000 s^c L£U - MS' Agsjy cM F L ! A G . M W C SBPX B I R T < 2 V dp® - |E ' !k s sjr, BB S S JD D£ RR gr Y A H£C £a R X jg w p^Sa C g P N G £R p y w LVg\2 B J R X Ai DPP Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 136 [Magistrale et al., 2000] and K. Olsen’s fourth-order staggered-grid finite- difference code [Olsen, 1994], we have constructed a RGT database for 64 broadband stations in the Los Angeles region (Figure 5.1). The synthetic seismogram excited by a point source with moment tensor My at position rs at time ts and recorded by a receiver located at rR is Sk it) - My dS jGki (fr, rs; t-ts), (5.1) where ds) is the gradient with respect to source coordinate rs. Seismic reciprocity implies: Gk i (rR , rs; t-ts) = Gik (rs, rR ; t-ts). (5.2) Therefore, using RGTs and the principle of reciprocity, we can obtain complete 3D synthetic seismograms at the 64 stations for any event located within our modeling box. c The source-coordinate gradient d 7 is computed by extracting Glk on a small source-centered volume and using a five-point-formula [Press, 2001] fi = (fi-2 -m -i+ m + i-fi+ 2yi2h + o(h \ (5.3) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 137 where / represents a certain component of a RGT, I is the grid index in each direction and h is the grid-spacing (Figure 5.2). The gradient computed using the five-point-formula is 4th -order. As shown in Figure 5.3, the synthetic seismograms computed using RGTs and reciprocity have almost no difference from the synthetics computed from a direct forward simulation using the same seismic source and structure models and the same finite-difference code. (/, m, f t + 1 ^ (/, m-1, n) (/, m, n) (/-l, m, n) ♦ (/+ 1 , m, n) ( I , m + 1 , n ) (ly m, 72-1) Figure 5.2 The spatial source-centered grid for computing source-coordinate gradients. The red dot in the center is the point-source location. The indices (I, m, n) represent values of spatial coordinates (x, y, z), respectively. The grid spacing h is defined as the length between the centers of two adjacent grid points. The RGTs for the 64 stations were computed on a regular mesh of about 36 million grid points with a grid spacing of 200 meters and a sampling rate of 0.01 second for 1 minute. On a Pentium 4 cluster with 64 nodes, it takes about 12,500 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 138 node-hours of total CPU time to compute all the RGTs. The computed RGTs were then sampled on the same mesh with a sampling rate of 0.1 second and archived on the Storage Resource Broke (SRB) system at the San Diego Supercomputer Center, where they occupy a total data volume of about 24 TB. The SRB provides high-level digital library functionality, including a maintained association of data and metadata and tools for rapid data queries and subset retrieval, which proves to be essential both for our seismic source inversions and for our tomography studies in the Los Angeles basin area [Zhao et al., 2005, Chapter 6], 0.02 0.01 «3 4-* 0 < D > -0.01 -0.02 < 0.05 2 " O < T J C C 9818433 BRE - A ... / \ A A / \ A ^ A ' 1 r 1 / \ / V / y ' v / y V I / V J \J \ r x A v -------- V 10 15 20 30 .......1 .....................................1 ......................................!.... ^ A A A ,, .... A / X / X .A V \j v \ / A 10 15 - 0.01 15 20 Time (s) F igure 5.3 Synthetic seismograms from the 03 September 2002 M L 4.8 Yorba Linda earthquake at station BRE (Figure 5.1). Blue solid line: synthetics computed directly from forward simulations by propagating seismic waves from the source location to the station at surface. Red dash line: synthetics computed using Receiver Green Tensors (RGTs), which were computed beforehand and stored at San Diego Supercomputer Center’s Storage Resource Broker (SRB) and by the application of reciprocity. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 139 5.4 CMT Inversion Using 3D Synthetics In general, the moment tensor My has 6 independent elements. For a pure- deviatoric source, we require the trace of My to vanish. For a pure double-couple source, we further require the determinant of My to be zero. Following Kikuchi & Kanamori [1991], we decompose a general moment tensor My into 6 elementary basis tensors: '0 1 o' '1 0 o' " 0 0 o' 1 0 0 ; M 2 : 0 -1 0 ; m 3 : 0 0 1 0 0 0 0 0 0 0 1 0 '0 0 1' '-1 0 o ' '1 0 o' 0 0 0 ; M 5 : 0 0 0 ; M 6 : 0 1 0 1 0 0 0 0 -1 0 0 1 And the synthetics for a general moment tensor My can be expressed as a linear combination of synthetics for each basis moment tensor: s(0 = f iang n(t). (5.5) n=1 Here gn(t), the synthetic seismogram for basis moment tensor Mn , is computed using equation (5.1)-(5.3) and the method described in the last section, and a„ is the coefficient we are trying to determine using observed seismograms. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 140 To parameterize the misfit between observed and synthetic waveforms, we adopt the Generalized Seismological Data Functionals (GSDF) technique [Gee & Jordan, 1992; Chapter 2], The two types of data functionals: the phase-delay time dtv and amplitude-reduction time Stq are measured on a set of body waves with good signal-to-noise ratio. The sensitivity derivatives of Stm with respect to source parameters such as moment tensor coefficients a„, shifts in centroid time t\ and centroid location n are computed by numerically differentiating synthetic seismogram s(t) with respect to these source parameters. This numerical differentiation is implemented by applying the reciprocity principle and equations analogous to (5.3). The linearized inverse problem can be expressed in the following form: Aomo + Aimi = d. (5.6) Here d is the data vector comprising GSDF measurements of dtV A , mo contains perturbations to the moment tensor coefficients an, mi specifies the centroid space and time shifts, Ao and Ai contain partial derivatives (sensitivity derivatives) of the data with respect to corresponding source parameters. A linear least-squares inversion can be carried out to solve this linear system for mo.i, then new synthetics can be computed using updated source parameters and the whole process can be iterated until the norm of the data vector d is small enough or no further improvement can be derived from the data. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Reproduced w ith permission o f th e copyright owner. Further reproduction prohibited without permission. 3 < J Q O 3 P- P O C D P . P 3 “ p p o 3 P . P O ' P P 3* "3 O P P 3 c ^ - C D C L 2 cd 3 cd p- cr n h - H C / D 2! S 3 cd r - f 3 o p 3 P - r ^ - =r cd o o £ L 3 cd o 3 “ p 3 Q- C D C D 3 C D a H o p cd =r c d cd P < C D 3 C D 3 < C D C/3 o ’ 3 O o C D C/3 cr C D O c b' S' p 3 P - 3 “ C D 34° w F igure 5.4 Using 3D synthetics computed from Receiver Green Tensors and reciprocity, we have determined CMT solutions for 149 local small to moderate-sized earthquakes (M L 3.0-4.9) in Los Angeles basin area. Focal mechanisms determined using first-motion data were used as initial solutions in an automatic gradient optimization procedure. Beach balls are plotted at the epicenters o f the earthquakes. 142 Shearer’s (2002) algorithm (HASH) was adopted to carry out the grid-search for suitable dip, rake and strike that best fit the P wave polarity information. The focal mechanism determined from first-motion data using HASH is then used as the initial solution in a gradient-based optimization procedure that minimizes GSDF measurements of amplitude-reduction time. In practice, the optimization procedure converges after about 10 iterations. We have successfully determined CMT solutions for 149 small to moderate-sized earthquakes (M l3.0 ~ 4.9) in Los Angeles basin area. Our solutions are shown in Figure 5.4. For those events that are also in Egill Hauksson’s focal mechanism catalog (Hauksson, 2000), our solutions are in good agreement with Hauksson’s solutions. 5.5 Removing Fault-plane-ambiguity from CMT Solutions 5.5.1 Finite M om ent Tensors To resolve fault-plane-ambiguity, we need to go beyond the point source representation of a seismic source. The lowest-order representation of a finite earthquake source is the second-order polynomial moment of the source-space time function termed Finite Moment Tensor (FMT) as in Chen et al., [2005] and Chapter 4. The Finite Moment Tensor extends the CMT to include the characteristic space time dimensions and orientation of the source. Our formulation of the FMT representation follows McGuire et al. [2001] and Chen et al., [2005]. If we Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 143 assume the stress glut is everywhere proportional to a constant seismic moment tensor, ( T(r,t) = M J / (r,t')dt’, (5.7) then the first-order polynomial moments of the source-space-time function provide the spatial and temporal centroids: H( 1 ,0 ) = ||/(r,O rrfV c?f=rj, (5.8a) = ^ f(r ,t)td V d t= tl . (5.8b) The lowest-order description of source finiteness is given by the second moments. The second central moment in time specifies the characteristic duration Tc of the event [Silver & Jordan, 1983] f i (0 ’2 ) = \ j f ( r , t ) ( t - tl)2dVdm (Tc/2 )2. (5.9) The second-order central moment in space is a symmetric second-order tensor, ^ < 2 0 ) = |J / ( r ,0 ( r -T jX r - r j ) TdV7fr = UAUT . (5.10) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 144 Here A is a diagonal matrix of eigenvalues, which gives us the characteristic dimensions of the seismic source, and U is an orthogonal matrix of eigenvectors, which gives us the orientation of the fault plane in space. For a simple rectangular dislocation with length L and width W, the eigenvalue matrix reduces to A=diag[(Lc / 2)2, (W c / 2)2,0]. (5.11) Here, Lc = L /V 3 is the characteristic length and Wc= W / a/3 is the characteristic width. The mixed space-time moment gives the directivity vector [Ben-Menahem, 1961] £ (U ) - J J y (r,t)(r-rj)(t - t j )dVdt =fi(0 '2 ) v d. (5.12) The directivity vector V d will lie in the plane of a simple dislocation. The magnitude of this vector is zero for a symmetric bilateral rupture and L fT c for a unilateral rectangular rupture along the length. 5.5.2 Sensitivity Derivatives for FM T Parameters To resolve fault-plane-ambiguity in our CMT solutions, we calculated partial derivatives of GSDF measurements with respect to second-order central moment in space ji(2’ 0 ) by applying reciprocity and numerically differentiating synthetic Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 145 seismograms with respect to source coordinates up to second order. The exact expressions for these sensitivity derivatives are given in McGuire et al., [2001], Chen et al., [2005] and Chapter 4. A 35t 73^' ( 1,0 ) + 35r /3p q ( 1,0 ) * asxq/a^1 '1 ’ A 3 5 ^ /3 ^ 3) + 38t /3|i( 1 ,1 ) q z aAa =HAA++++++-H-++-I * 3 8 t /3p q m- . (2,0 ) 180 360 A 3ox /3p q “ xy 3 8 t / 3 ( i q 180 360 0.06 A 3 0.3 * A * * 0.04 2 a Aa A A A 0.2 + + A * A * 0.02 < 1 < 1 < 1 1 A * A * 0.1 A A + A A Q i A a ^ aA - ---------Av 0*& — A -----$------ A ------ 0« * * * A & — A - --A- - - A ------- * A * A A * A 0.02 A A A -1 -0.1 A A + A A A A 0.04 -0 .2 + + A + + + 0 06 - 3 0 180 360 0 180 360 0 180 360 F igure 5.5 Plots of the sensitivity coefficients for 5tq as a function of azimuth for P waves from the Yorba Linda earthquake. First column: coefficients for a centroid shift (triangles: derivative relative to a perturbation towards north; crosses: a perturbation upward). Second column: coefficients for the directivity moment p,(11) (triangles: rupture propagates towards north; stars: rupture propagates towards east; crosses: rupture propagates upward). Third column: coefficients for second-order spatial moment |i,< 2 ,0 ) (stars: pxt(2 ’ 0); triangles: px / 2,0); crosses: pK< 20)). Upper row: 0.2 Hz, lower row: 2.0 Hz. Dash lines indicate zero. The sensitivity of Stq with respect to directivity increases with frequency, while the sensitivity to centroid shift and second-order spatial moment decreases with frequency. The apparent scatter in the first column is caused by singularities at 30°, 120°, 210° and 300° azimuth, which are the nodal planes of P wave radiation pattern. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 146 Examples of partial derivatives of amplitude-reduction time Stq with respect to first- and second-order moments are shown for the P waves generated by the 3 September 2002 ML 4.8 Yorba Linda earthquake in Figure 5.5. This earthquake has a strike-slip focal mechanism. The strike is 30° NE and the rake angle is 0°. The sensitivity for second-order spatial moment p (2 ,0 j has a cos(20) dependence on station azimuth 0 because of constructive and destructive interference effects caused by waves generated from a finite source. And this example also shows that at lower frequency (0.2 Hz), where directivity effect caused by p (1 1 ) is small (second column of Figure 5.5), it is still possible to resolve fault-plane-ambiguity by careful analysis of the interference effect caused by source finiteness p (2 0), which has a relatively larger contribution to amplitude-reduction time at lower frequency. After correcting GSDF measurements for errors in 0th order moment (the moment tensor) and 1st order moments (centroid shifts in space and time), we may attribute Stp,q residues to contributions from finite moment tensors A2 m2 = d, (5.13) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 147 where d is the data vector of dllxq, m2 contain finite moment tensor parameters and A2 is the corresponding sensitivity derivative matrix. At lower frequency (0.1 ~ 0.3 Hz), we may ignore contributions from the directivity moment p (U). The sensitivity for second-order temporal moment p (0 '2 ) has no azimuthal dependence [equation 8a, Chen et al., 2005; Chapter 4], Therefore, the azimuthal variations of 8tm are largely due to second-order spatial moment p (2,0). The two nodal planes given by the CMT solution have mutually orthogonal second-order spatial moments and will predict different azimuthal dependences of GSDF measurements. For example, for a pure bilateral rupture with a strike of 30°, the amplitude of finite-source waveform relative to the point-source synthetic waveform should be maximized at azimuths of 120° and 300° (minimum of amplitude-reduction time Stq), owing to constructive interference, while slip on the conjugate fault plane would generate a maximum at 30° and 210°. We therefore implemented an automated procedure to systematically test the two nodal planes of our CMT solution using the sensitivity kernel for second-order spatial moment p (2,0). The nodal plane that provides a better azimuthal correlation with Stm observations is identified as the preferred fault plane. Our results for 40 small to moderate-sized earthquakes in Los Angeles basin area are shown on Figure 5.6. The probability of our solution is computed using a bootstrap method. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Reproduced w ith permission o f th e copyright owner. Further reproduction prohibited without permission. 241 ° 30' 242° 00' 242° 30' 243° 00' Figure 5.6 We resolved fault-plane-ambiguity for 40 earthquakes in Figure 5.4. The preferred fault planes are highlighted with red-lines on the beachballs. The numbers besides each beach ball show the probability for our preferred fault plane. They were computed using a boot strap method. £ oo 149 5.6 Discussion The uncertainty of our solution highly depends upon the azimuthal coverage of the seismic network. And our results are more robust for strike-slip events than for dip-slip ones due to the fact that most of our seismic stations are on the surface of the earth. The errors in our 3D structure model might contaminate our solution for source parameters. To improve the signal-to-noise ratio, we used only low-frequency data measured at 0.1 ~ 0.3 Hz to resolve fault-plane-ambiguity. As we improve regional structure model, we may invert GSDF measurements made at higher frequencies ranges and we will be able to routinely recover complete FMT solutions for small to moderate-sized earthquakes in well-instrumented areas such as Southern California. When applied to a large number of earthquakes, this new technique will provide additional constraints on regional tectonics. One example is the incipient faulting on the Puente Hills-Fontana seismicity trend located to the southwest of the right- lateral San Jacinto fault and the west-striking and north-dipping Cucamonga fault as shown on Figure 5.4 and 5.6. This seismicity trend is not associated with any mapped fault traces. We resolved fault-plane-ambiguity for 18 small events in this swarm of earthquakes and they all show right lateral source mechanisms conjugate to the southwest trend of the epicenter distribution (Figure 5.6) and parallel to San Jacinto fault. This Fontana seismicity trend stops quite abruptly at Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 150 the Puente hills area. About 20 miles further to the southwest, in the Yorba Linda area, our study of the 03 September 2002 ML 4.8 Yorba Linda earthquake and several surrounding smaller earthquakes shows left-lateral faulting conjugate to the nearby right-lateral Whittier fault. Whether there is any connection between these two seimicity clusters is interesting from a tectonic point of view. Recent studies show that a high proportion of north-south shortening in the Los Angeles metropolitan region is accommodated by conjugate strike-slip faulting and east- west escaping of crustal blocks [Walls, et al., 1998; Chapter 3], The Puente Hills- Fontana seismicity trend and the left-lateral faulting in Yorba Linda area might indicate a developing weak zone associated with such “escaping” deformation. By recovering finite source properties for small to moderate-sized earthquakes, our technique provides new information for constraining regional tectonics and fault structure, which are essential for seismic hazard analysis and mitigation. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 151 5.7 References Aki, K and P. G. Richards, Quantitative Seismology, 2n d edition. University Science Books, Sausalito, CA, 2002. Bowman, D. D. and G. C. P. King, Accelerating seismicity and stress accumulation before large earthquakes. Geophys. Res. Lett., 28, 4039-4042, 2001. Chen, P., T. H. Jordan and L. Zhao, 2005. Finite Moment Tensor of 3 September 2002 Yorba Linda Earthquake. Bull. Seismol. Soc. Am., 2005, in press. Dziewonski, A. M., T. A. Chou, and J. H. Woodhouse, Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. J. Geophys. Res. 86, 2825-2852, 1981. Edelman, S. H., A. Liu, and R. D. Hatcher, The Brevard in South Carolins and adjacent areas: An Alleghanian orogen-scal dextral shear zone reactivated by a thrust fault. J. Geology 95, 793-806, 1987. Gee, L. S., and T. H. Jordan, Generalized seismological data functionals. Geophys. J. Int. Ill, 363-390, 1992. Hardebeck, J. L., P. M. Shearer, A new method for determining first-motion focal mechanisms. Bull. Seismol. Soc. Am. 92, 2264-2276, 2002. Harris, R. A., Introduction to special section: Stress triggers, stress shadows, and implications for seismic hazard. J. Geophys. Res., 103, 24347-24358, 1998. Hauksson, E., Crustal structure and seismicity distribution adjacent to the Pacific and North America plate boundary in southern California. J. Geophys. Res. 105, 13875-13903, 2000. Kikuchi, M. and H. Kanamori, Inversion of complex body waves -III, Bull. Seismol. Soc. Am., 81, 2335-2350, 1991. Liu, Q.-Y., J. Polet, D. Komatitcsh and J. Tromp, Spectral-element moment tensor inversions for earthquakes in Southern California. Bull. Seismol. Soc. Am., 94, 1748-1761, 2004. Magistrale, H., S. Day, R. W. Clayton, and R. Graves, The SCEC southern California reference 3D seismic velocity model Version 2, Bull. Seismol. Soc. Am. 90, no. 6B, S65-S76, 2000. McGuire, J. J., L. Zhao, and T. H. Jordan, Teleseismic inversion for the second- degree moments of earthquakes. Geophys. J. Int. 145, 661-678, 2001. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 152 Mori, J., Rupture directivity and slip distribution of the M 4.3 foreshock to the 1992 Joshua Tree earthquake, Southern California. Bull. Seismol. Soc. Am. 86, 805-810, 1996. Mori, J. and S. Hartzell, Source inversion of the 1988 upland, California, earthquake: determination of a fault plane for a small event, Bull. Seismol. Soc. Am. 80, 507-518, 1990. Olsen, K. B., Simulation of three-dimensional wave propagation in the Salt Lake Basin, Ph.D. Thesis, University of Utah, Salt Lake City, Utah, 157p, 1994. Pasyanos, M. E., D. S. Dreger, and B. Romanowicz, Towards real-time determination of regional moment tensors. Bull. Seismol. Soc. Am. 86, 1255-1269, 1996. Press, W. H., S. A. Teukolsky, W. T. Vetterling and B. R. Flannery, Numerical Recipes in Fortran 77 - the art of scientific computing. Cambridge University Press, 2001. Royden. L., 1984. Thin-skinned extension of Vienna basin [abs], Amer. Assoc. Petroleum Geologists, Bull. 68, 523. Seeber, L., and J. G. Armbruster, The San Andreas fault system through the Transverse Ranges as illuminated by earthquakes. J. Geophys. Res. 100, 8285- 8310, 1995. Silver, P., and T. H. Jordan, Total-moment spectra of fourteen large earthquakes. J. Geophys. Res. 88, 3273-3293, 1983. Webb, T. H., and H. Kanamori, Earthquake focal mechanisms in the eastern Transverse Ranges and San Emigdio Mountains, Southern California, and evidence for a regional decollement. Bull. Seismol. Soc. Am. 75, 737-757, 1985. Zhao, L., T. H. Jordan, K. B. Olsen and P. Chen, 2005. Frechet Kernels for Imaging regional earth structure based on three-dimensional reference models. Bull. Seismol. Soc. Am, submitted. Zhu, L. and D. Helmberger, Advancement in source estimation techniques using broadband regional seismograms, Bull. Seismol. Soc. Am., 86, 1634-1641, 1996. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 6 153 Fully Three-Dimensional Tomography of Los Angeles Basin Area 6.1 A b stract Reliable prediction of strong motion caused by damaging earthquakes requires seismologists to provide realistic seismic structural models through which seismic waves propagate. In regions such as Southern California, preliminary 3D earth models are already available, and efficient numerical methods have been developed for solving the point-source forward problem. We describe an inversion procedure that utilizes these capabilities to improve 3D earth models. Our data are time- and frequency-localized measurements of the phase and amplitude anomalies relative to synthetic seismograms computed from the 3D elastic starting model. And we have computed 3D sensitivity kernels for our time- frequency-dependent measurements by convolving source-side wave fields with receiver-side Green tensors, which comprise the spatial-temporal displacement fields produced by the three orthogonal unit impulsive point forces acting at the receiver positions. We have obtained a preliminary 3D perturbation to the 3D reference model SCEC CVM3.0 for Los Angeles Basin area. Our result shows Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 154 that in general SCEC CVM3.0 is too slow and might even be lower in the basin area. To our knowledge, this is the first “fully 3D” inversion of waveform data for regional earth structure. And the method developed in this study should be widely application to other areas of seismological interest. 6.2 Introduction Earthquake hazard analysis depends on seismology to quantify ground motions. Over the last two decades, engineering techniques that enhance protection against earthquake damages have advanced considerably and engineers are expected to design buildings that can sustain a certain level of functionality under various scenarios of earthquake shaking. This type of “performance-based-engineering” requires earthquake scientists to provide a more detailed description of ground motion than classical intensity measures at frequencies important for building performance (0.3 to 10 Hz). With the advancement of high-performance computing and numerical techniques, such as the finite-difference method (Frankel & Vidale, 1992; Olsen, 1994; Graves, 1996), the finite-element method (Bao et al., 1998; Aagaard et al., 2001) and the spectral-element method (Komatitsch & Vilotte, 1998; Komatitsch et al., 2003), physics-based simulations of ground motions for relatively well-characterized regions, such as the Los Angeles metropolitan area, are now feasible. The recent “Terashake” simulation tries to model the complete time history of the ground velocities that can be Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 155 expected to shake buildings and infrastructure in the Los Angels basin area if a magnitude 7.7 earthquake strikes Southern California. The accuracy of such ground-motion simulations depends heavily upon the seismic velocity model employed to propagate seismic waves within. It is well known that sedimentary basins can amplify seismic waves (Kawase, 1996; Graves et al., 1998; Pitarka et al., 1998; Davis et al., 2000). Earthquakes can excite resonance in deep basins and shake the sediments like jelly. Seismic waves can also be amplified through focusing effects by some lens-like 3D heterogeneity (Gao et al., 1996) and fault-bounded edges of sedimentary basin can generate destructive surface waves known as the basin-edge effects (Graves et al., 1998). Therefore, for areas with complex 3D subsurface structures and high earthquake risk such as the Los Angeles basin area (Figure 6.1), it is very important to map out small-scale (on the order of a few kilometers or less) 3D heterogeneities in both P and 5 wave speeds in order to reliably predict strong ground-motion at frequencies important for engineering purposes. Up to now, three-dimensional studies of Los Angeles basin area are primarily focused on the P wave speed and available 3D velocity models largely come from three sources: (a) P-wave and/or 5-wave travel-time tomography (Hauksson, 2000); (b) information about the age and depth of sediments from borehole measurements, oil and water exploration and geological studies (Magistrale et al., 1996; Magiatrale et al., 2000); (c) linear-array active-source surveys (Fuis et al., Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 156 2001; Godfrey et al., 2002), direct measurements from sonic logs and stacking velocities from petroleum industry reflection profiles (Suss & Shaw, 2003). Thus, current available 3D velocity models for Los Angeles basin area either do not have enough resolution [the tomogram from Hauksson (2000) has a horizontal resolution of 15 km] or do not cover the entire area. Furthermore, the S velocity models in all those studies cited in (b) and (c) are not directly constrained by data, but derived from P velocity through certain empirical scaling laws that relate S 241" 00' 241" 30' 242" 00' 242" 30' 243*00' 241‘ 00' 241° 30’ 242" 00' 242° 30' 243" 00' F igure 6.1 Los Angeles basin and surrounding area. Black triangles show the location of California Integrated Seismic Network broadband stations. The focal mechanisms o f 72 earthquakes used in this study are plotted at their epicenters. Red lines show the 1073 source-station paths that cover the Los Angeles basin area. The inner box is the region of our tomography model. velocity to P velocity. As a result, current 3D models are insufficient for reliable prediction of ground-motion that can lead to effective seismic hazard analysis. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 157 The purpose of this study is to conduct a detailed and comprehensive investigation for 3D seismic velocity structure in Los Angeles basin area using a full 3D approach: imaging the 3D structure using true 3D Frechet kernels computed in a 3D reference model. This new full 3D approach eliminates the high frequency and structural averaging approximations that are ubiquitous in current tomography studies and provides a suitable tool for high-resolution structural imaging at regional scale. 6.3 Reference Model The SCEC Community Velocity Model Version 3.0 (CVM3.0) (Magistrale et al., 2000) serves as a unified 3D reference model for regional waveform modeling in Southern California. Figure 6.2 shows a fence diagram of P-wave velocity in and around the Los Angeles basin area. Within major sedimentary basins (Los Angeles basin, Ventura basin, San Gabriel Valley, Chino basin, San Bernardino Valley and the Salton Trough), the model is constructed using an empirical relation between sediment age, depth and P-wave seismic velocity. The sediment age and depth information were obtained from oil and water exploration activities and other geologic studies. Outside the basins, the model is based on three- dimensional P-wave travel-time tomography results of Hauksson (2000), which has a horizontal resolution of 15 km and a vertical resolution of 2-9 km. The Moho is represented by a three-dimensional surface with variable depth, which is determined using the receiver-function technique (Zhu & Kanamori, 2000). The Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 158 density is derived from P velocity by assuming the empirical relation of Nafe & Drake (1960) and Poisson’s ratio is determined from density using the empirical relation of Ludwig et al. (1970). The S velocity is then computed using the P velocity and Poisson’s ratio. F igure 6.2 The fence diagram o f SCEC CVM 3.0 S velocity model (Magistrale et al., 2000). Cross section locations are shown as red lines in lower left panel. S velocity is scaled from P velocity using empirical relations. P velocity in the basins are controlled by rule-based models and by tomographic results outside the basins (Hauksson, 2000). This hybrid three-dimensional velocity model CVM3.0 appears to be an acceptable reference model. We can derive three-dimensional perturbations using this reference model and refine it through inversions. In Chapter 2, we carried out a detailed quantitative evaluation of CVM3.0 by comparing synthetic seismograms computed using CVM3.0 and Olsen’s (1994) finite-difference code Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 159 with observed seismograms. Our analysis on about 9955 body waves including 4070 P waves, 1235 SV and 4110 SH waves shows that CVM3.0 provides a substantially better fit to observed waveforms than either laterally homogeneous one-dimensional model SoCaL (Dreger & Helmberger, 1993) or path-averaged one-dimensional model AID and we concluded that CVM3.0 is sufficient to be used as a reference model in linearized tomography inversion. In addition, we also found that the polarization anisotropy in observed data is very small and it is not necessary to consider anisotropy in initial inversions. 6.4 D ata We selected 72 small to moderate-sized earthquakes (3.0 < Ml < 5.0) in the Los Angeles basin area and we collected 5168 seismograms with high signal-to-noise ratio from 45 California Integrated Seismic Network (CISN) three-component broadband stations (Figure 6.1). The total number of source-station paths is 1073 and the basin area is well covered. The horizontal seismograms were rotated into radial and transverse components and all recorded seismograms were low-pass filtered with a fourth-order Butterworth filter having comer frequency at 1.0 Hz. To determine the focal mechanisms for these 72 events, we implemented a two- step procedure. First, we determined an initial focal mechanism using first-motion data picked on vertical seismograms. Hardebeck & Shearer’s (2002) algorithm (HASH) was adopted to carry out the grid-search. Three-dimensional synthetic Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 160 seismograms were then computed using this initial solution. Second, the focal mechanism obtained from the first step was then refined in a multi-step gradient optimization process based on minimizing the amplitude-reduction times (Gee & Jordan, 1992; Chapter 2) measured on various body waves (Chapter 5). To account for three-dimensional structural heterogeneity, the Frechet derivatives with respect to moment tensor parameters were computed numerically using three-dimensional synthetics, which were calculated using CVM3.0 and Olsen’s finite-difference code. Generally, after about 10 iterations, we could obtain a stable moment tensor solution. Three-dimensional synthetics were then computed using the refined moment tensor solution. Horizontal synthetic seismograms were rotated into radial and transverse components and all synthetics were low-pass filtered using the same filter we used for observed seismograms. The data used for our tomographic inversion in this study are 9955 frequency- dependent phase-delay times and 9955 amplitude-reduction times (Gee & Jordan, 1992; Chapter 2) measured on 4070 P waves, 4110 SH waves, 1235 SV waves and 540 other phases. Different from previous applications of the GSDF method on teleseismic data, the isolation filters used in this study were not constructed by summing normal modes weighted according to group velocity but by windowing the complete synthetic seismograms and isolating the segments that we were interested in. The isolation filters were then cross-correlated with the data seismograms, and the resulting cross-correlagrams were windowed and narrow band filtered at 0.2, 0.4, 0.6, 0.8 and 1.0 Hz. The phase-delay times and Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 161 amplitude-reduction times were determined by fitting Gaussian wavelets to the windowed and narrow-band filtered cross-correlagrams. Details about the GSDF analyzing procedure can be found in Gee & Jordan (1992) and Chapter 2 of this dissertation. Some key features in this data set are illustrated in Figure 6.3-6.4, which shows examples of observed and synthetic seismograms as well as isolation filters for some representative paths. The implication of these GSDF observations for basin structure can be understood by examining a subset of their Frechet kernels (Figure 6.5-6.7). 6.5 Frechet Kernels The Frechet kernel provides a linear relationship between the measurement and the velocity perturbation. Generally, we can express such a linear relation as * P ,q = fK M(r,m 0)-8m(r)dV(r). (6.1) V Here, < 5 tp,q is the GSDF measurements of phase-delay time or amplitude-reduction time, KP ,q is the corresponding Frechet kernel, which is computed using the reference model mo, and 6m represents the perturbation to the reference model. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 162 o -0.5 -1.5 *- 0.2 0.8 1 0.' 1 .4 F igure 6.3 The GSDF processing for the P wave, (a) the observed seismogram, the complete synthetic and the isolation filter, (b) Solid line: the autocorrelation of the isolation filter; dash line: the windowed autocorrelation, (c) Solid line: the cross correlation between the isolation filter and the observed data; dash line: the windowed cross-correlation, (d) Solid line: the narrowband filtered windowed autocorrelation; dash line: the best-fit Gaussian wavelet, (e) Solid line: the narrowband filtered windowed cross correlation, central frequency 0.2 Hz, bandwidth 0.02 Hz; dash line: the best-fit Gaussian wavelet, (f) Circle: Stp; triangle: 5tq. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 163 0.04 0.02 s(t) o 0.02 s(t) -0.04 -0.06 -0.08 -0.1 -0.12 5 1 O 15 20 25 0 A A F igure 6.4 The GSDF processing for the X phase, (a) the observed seismogram, the complete synthetic and the isolation filter, (b) Solid line: the autocorrelation of the isolation filter; dash line: the windowed autocorrelation, (c) Solid line: the cross correlation between the isolation filter and the observed data; dash line: the windowed cross-correlation, (d) Solid line: the narrowband filtedash line windowed autocorrelation; dash line: the best-fit Gaussian wavelet, (e) Solid line: the narrowband filtedash line windowed cross-correlation, central frequency 0.2 Hz, bandwidth 0.02 Hz; dash line: the best-fit Gaussian wavelet, (f) Circle: Stp; triangle: Stq. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 164 The integration is over V, which is the volume of the earth under study, and r is the spatial coordinate. The three major components in this equation: model perturbation 8m, sensitivity kernel Kp ,q and the reference model mo, are intrinsically three-dimensional functions of space coordinate r. But depending on research interests or computational resources available, previous tomographic studies usually assume either the reference model or the sensitivity kernel or both to be one- or two- dimensional objects. For example, in Hauksson (2000), P and S wave travel-time anomalies determined by the instantaneous onset times of the signals on relatively short-period or broadband seismic records were inverted for three-dimensional structures in the central Transverse Ranges and the Los Angeles basin area. The geometric ray theory was adopted to describe signal propagation and the Fermat approximation (e.g. no ray path change due to velocity perturbation) was assumed in computing Frechet kernels for P and S travel time measurements. Geometric ray theory together with Fermat approximation produces kernels that are non-zero only on the geometric ray path by averaging structure in the plane perpendicular to the ray path. It is valid only if the length scale of the heterogeneity is much larger than the width of the first Fresnel zone. To reconstruct small-scale heterogeneities as attempted in this study, the effect of finite-frequency of wave propagation is expected to be significant. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 165 To achieve higher resolution of the three-dimensional structure, we need to account for the finite-frequency effects of wave propagation. Therefore we adopted Olsen’s (1994) fourth-order, staggered-grid, finite-difference code to compute synthetic seismograms and Green tensors in our three-dimensional reference model CVM3.0. To fully utilize finite-frequency waveform information, we computed Frechet kernels for our frequency-dependent phase-delay times and amplitude-reduction times using three-dimensional synthetics and Green tensors together with the Bom approximation (Marquering et al., 1999), which provides a linearized relation between waveform perturbation and the three-dimensional velocity perturbation. In our fully three-dimensional approach, the three major components of the tomography problem (e.g. the reference model, the model perturbation and the Frechet kernel) are all three-dimensional. The formulation and some numerical experiments for our true 3D kernel can be found in Zhao et al., (2005). To compute the kernel, we need the full wave-field generated by the earthquake source and the Green tensors from every grid point within the simulation volume to the receiver. In practice, these Green tensors from every grid-point to receiver were obtained by applying the reciprocity principle and the Green tensors computed using unit impulsive point sources located at the receiver. The kernel was calculated by convolving the gradient of the forward wave-field with the Green tensors. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 166 DLA. surface * vertical 15 10 5 < r r - 5L DLA. Time (sec) surface S Figure 6.5 The sensitivity kernel for the vertical com ponent P wave generated by 09/03/2002 Yorba Linda earthquake and recorded at CISN station DLA. Upper panel: sensitivity kernel for phase-delay time; middle panel: the synthetic seismogram with the P wave highlighted in red; lower panel: sensitivity kernel for amplitude-reduction time. W arm colors represent negative sensitivity indicating that a velocity increase leads to an advance in arrival tim e and an increase in amplitude; cool colors represent positive sensitivity indicating a velocity increase leads to a delay in arrival time and a reduction in amplitude. Some examples of our true three-dimensional Frechet kernels are displayed in Figure 6.5-6.7. In Figure 6.5 we show the kernels for the vertical-component P wave generated by the September 3, 2002 Yorba Linda earthquake and recorded at basin station DLA. These kernels are for the GSDF measurements made at 1.0 Hz. The sensitivity for phase-delay time have the famous “banana-doughnut” shape, it is distributed around the P wave ray path and spreads over a wide area. The seismic waves traverse both basin and non-basin structures and most of the phase-delay time sensitivity is concentrated in the shallow zone in the Los Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 167 DLA surface •% 5 tP • S OL/^Jj**# t fit, .... p l\ ^ /v, l\ I\ ft — H w F D L A ^ f c . W t 5 t, 5 10 15 D LA. Time (sec) surface l s 8 t? Figure 6.6 The sensitivity kernel for the vertical component X phase generated by 09/03/2002 Yorba Linda earthquake and recorded at CISN station DLA. Upper panel: sensitivity kernel for phase-delay time; middle panel: the synthetic seismogram with the X phase highlighted in red; lower panel: sensitivity kernel for amplitude-reduction time. Left column: looking perpendicular to the source-receiver plane; right column: map views Angeles basin beneath the receiver. The measured phase-delay time of about -0.8 second (Figure 6.3) suggests the slowing-down of the synthetic wave-field is likely acquired in the basin near the station. In other words, the P-wave velocity of CVM3.0 is perhaps too low in the basin. This example demonstrates that this kind of kernels is very useful for identifying origins of potential problems in the reference model. Another example of our three-dimensional kernels is shown on Figure 6.6. These kernels are for the arrival at about 10 seconds (Figure 6.4) on the same seismogram as Figure 6.5. The distribution of the sensitivity shows that most of the energy in this phase comes from the scattering caused by the heterogeneities near the surface around the Los Angeles basin. This kind of kernel is very helpful Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 168 for us to identify and analyze arrivals that cannot be easily characterized using geometric ray theory, and provides us a convenient tool for studying wave propagation in complex geological environment such as the Los Angeles basin area. BRE(Z). 113W30 118W I17WJ0 12 14.5 LLS(Z) 7.5 10.5 0U(T) AM. iff naW 3o n s w o o ii7W aa TT— ......... RPV(Z) dtn BRE(T) Stq " 3 ^ nr 5a T T T 113W30 118WG0 117W30 8.7 11.2 8 (3 5a sp DLA(Z) iisW 30 n e w tso H 7W io Figure 6.7 Kernel gallery. First column: map views of source-station path; second column: observed (black) and synthetic (red) seismograms; third column: phase-delay time sensitivity kernels for phases within the tim e window indicated in the second column; fourth column: amplitude-reduction time sensitivity kernels. The kernels for BRE(T) and OLI(T) are sensitivities to shear wave velocity perturbation 8p, the rest o f the kernels are sensitivities to P velocity perturbation 8a. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 169 More examples of our Frechet kernels both for P waves and for SH waves are shown in Figure 6.7. The sensitivity for the two SH waves on BRE(T) and OLI(T) does not have the “banana-doughnut” shape as we have seen in the sensitivity kernel for the P wave on Figure 6.3. This is because the reference model is a highly variable three-dimensional model and the SH waveforms as we identified in the time windows are no longer from the direct S wave arrivals only. Zhao et al., (2005) provided a detailed theoretical formulation for our true 3D Frechet kernels and conducted comprehensive numerical experiments to investigate the behavior of these 3D sensitivity kernels. The algorithm developed in that paper provides a very general and automatic process for regional waveform tomography and is the theoretical foundation for the 3D tomographyic inversion in this study. 6.6 Three-dimensional Seismic Structure of Los Angeles Basin Using the three-dimensional sensitivity kernels for our -20,000 GSDF dtm measurements, we have derived three-dimensional perturbations to the three- dimensional reference model CVM3.0. To our knowledge, this is the first successful application of the fully three-dimensional tomographic technique to a regional data set. 6.6.1 Param eterization and Inversion Procedure The geographic area for our inversion is shown on Figure 6.1. Our model covers the region from 118.6°W to 117.4°W and from 33.6°N to 34.3°N and extends Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 170 down to the depth of 35 km. Along longitude, latitude and depth, this 111 km x 78 km x 35 km volume is divided into 55 x 37 x 17 (34,595 in total) uniform blocks. The size of each block is about 2 km in each direction. The sensitivity kernels were computed on a much denser grid with 400-meter grid spacing. To compute the sensitivity for each block, we simply average the sensitivity on all the grids within this block. A constrained linear least-squares inversion procedure is used to obtain velocity perturbations in all the blocks. The perturbations were assumed to satisfy the linear system A 5m + B 5n = d (6.2) where 5m is the model perturbation vector of dimension M = 34,595 and d is the data vector of dimension D = 19,910 (9955 frequency-dependent phase-delay times and 9955 amplitude-reduction times). The D x M matrix A is the sensitivity matrix whose general element Ay is the sensitivity of model block j for the ilh datum. The vector 5n comprises “nuisance parameters” that account for quantities not included in model parameter 5m but could potentially contaminate the data. In this study, 5n contains the source origin time and perturbations to source centroids. The matrix B contains the sensitivity derivatives of data vector with respect to nuisance parameters. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 171 The effect of nuisance parameters can be removed by applying a projection operator that annihilates the term B 8n. The appropriate annihilator is QB = I - B (b t B r 1 b t , where I is the identity matrix. And the reduced system of equations can be written in the form A' 5m = d' (6.3) where A' - Qb a and d' - Qb d. The details about this type of “denuisancing” procedure can be found in Jordan & Sverdrup (1981). We include weighting and damping by appending M rows to the above linear system: A '^ l „ ( C'1 d’ ^ V d AS 5m = v d 0 (6.4) Here, Cd ' is the a priori data covariance matrix that weights the data and S is a M x M matrix of smoothing constraints that minimize the difference between velocity perturbations in adjacent model blocks. The damping parameter X controls the tradeoff between satisfying the smoothing constraints and fitting the data. The damped linear least-squares inversion was carried out using Paige & Saunders’s (1982) LSQR method. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 172 D E Depth 5km Depth 9km 34° 00’ 33° 40' -118° 20’ -118° 00’ -117° 40' -117° 20' Depth 13km 118° 20’ -118° 00' -117° 40’ -117° 20' Depth 17km 34° 00' - 33° 40’ - -118° 20’ -118*00’ -117” 40’ -117° 20’ -118° 20’ -118° 00’ -117° 40’ -117*20' 118° 20’ -118*00’ -117° 40’ -117° 20’ -118° 20’ -118° 00’ -117° 40’ -117*20’ -0.04 -0.02 0.00 0.02 0.04 Sot/a {%) C • I f . IT _ X 1 ______ 1............. L...... ........._L -1184 -118.2 -118 -1178 -117.6 -117.4 , 10- 15- i.S fc J E i. ■ i.„." „ i l 1184 -118.2 -118 -1178 -117.6 -117.4 5 - .1 0 - ! ...............1 ..............1 .............. 1 J I -1184 -118.2 -118 -1178 -117.6 -1174 5 E io I 33.8 33.9 34 34.1 F ig u re 6.8 LABF3D, inverted 3D P velocity perturbation relative to 3D reference model CVM3.0. W arm colors represent negative perturbation; cool colors represent positive perturbation. Upper plots: map views at 4 different depth; lower plots: cross sections along 5 profiles, the locations of these profiles are shown as green lines on the upper-left panel of the map view plot. 33.8 33.9 34 34.1 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 173 Input ----------- g- ----------------------- ,------------------------ ,------------------------ , -118° 20' -118" 00' -117° 40' -117° 201 Depth 7km Depth 5km -118° 20' -118" 00' -117" 40' -117' 20' Depth 9km F igure 6.9 Checker board test with cell size 12 km by 12 km. The input model is shown on the upper-left plot. Map views of the recovered model are shown at three different depths. Input Depth 5km I V » i" -118’ 2<r -118° 00' -117*40' -117*20' Depth 7km -118° 20' -118" 00' -117" 40' -117° 20' Depth 9km 34° 00' 33° 40’- I -118°20’ -118" 00' -117*40' -117" 20’ -118°20' -118"00' -117"40' -117° 20' F igure 6.10 Checker board test with cell size 10 km by 10 km. The input model is shown on the upper-left plot. Map views o f the recovered model are shown at three different depths. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 174 Input Depth 5km 34" 00' 33" 40' -118“ 20' -118“ 00' -117“ 40’ -117“ 20’ -118“ 20' -118“ 00' -117“ 40’ -117“ 20’ Figure 6.11 Checker board test with cell size 8 km by 8 km. The input model is shown on the upper-left plot. Map views of the recovered model are shown at three different depths. Input 34“ 00’ - 33’ 40' - 1 -118“ 20' -118“ 00' -117“ 40' -117“ 20’ Depth 7km Depth 5km -118“ 20’ -118‘ 00’ -117“ 40’ -117“ 20’ Depth 9km 34“ 00' - 33° 40' 4 -118“ 20' -118'00' -117' 40’ -117“ 20' -118“ 20' -118' 00' -117'40' -117" 20' Figure 6.12 Spike test with cell size 12 km by 12 km. The input model is shown on the upper-left plot. M ap views of the recovered model are shown at three different depths. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 175 6.6.2 M odel LABF3D The smoothness of the inverted 3D model depends on the damping parameter X we choose during the inversion. As we increase X, the fit to the data deteriorates. Our preferred model LABF3D is shown as a 3D perturbation to the 3D reference model CVM3.0 in Figure 6.8. This model gives a relatively good overall fit to the data with a total variance reduction of about 20%. This small variance reduction is partly due to the fact that our starting model is already a good 3D model. As discussed in Chapter 2, the variance reduction from laterally homogeneous ID velocity model SoCaL (Dreger & Helmberger, 1993; Chapter 2) to CVM3.0 is more than 60% for phase-delay time measurements. The primary feature of LABF3D is the strong positive anomaly in and around the basin area. This suggests CVM3.0 is too slow in the basin. As we have discussed earlier, in CVM3.0, seismic velocities in sedimentary basins were not determined using direct tomography technique but by the application of empirical rules that relate sediment age and depth to P wave velocity. Our results suggest those empirical relations might underestimate the true velocity of wave propagation. The resolving power of our data set has been investigated through a series of checkerboard tests and a spike test as shown in Figure 6.9-6.12. Patterns with horizontal wavelength of 8 km are well resolved above 13 km depth. The Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 176 resolution decreases with depth below 15 km. This is due to the fact that 82% of the earthquakes in our data set are above 15 km. 6.7 Conclusion In this study, we have derived a 3D perturbation to a 3D starting model using 3D kernels numerically computed by the method of Zhao et al. (2005) that accounts for all elastic-interaction effects. To our knowledge, this is the first application of“full 3D” seismic tomography to a regional data set. The preferred model LABF3D is represented as a 3D perturbation to the reference 3D model CVM3.0. A prominent feature of LABF3D is that the velocity in and around the basin area is higher than CVM3.0, suggesting the empirical rules that were used in CVM3.0 to relate P wave velocity to the age and depth of sediments might underestimate true wave speed. The method developed in this project should be widely applicable to other tomography studies such as the US Array component of the EarthScope project. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 177 6.8 References Aagaard, B. T., J. F. Hall and T. H. Heaton, 2001. Characterization of near-source ground motions with earthquake simulations. Earthquake Spectra, 17(2), 177- 207. Bao, H., J. Bielak, O. Ghattas, L. F. Kallivokas, D. R. O’Hallaron, J. R. Shewchuk and J. Xu, 1998. Large-scale simulation of elastic wave propagation in heterogeneous media on parallel computers, Comput. Methods Appl. Mech. Engrg., 152, 85-102. Davis, P. M., J. L. Rubinstein, K. H. Liu, S. S. Gau and L. Knopoff, 2000. Northridge earthquake damage caused by geologic focusing of seismic waves. Science, 289, 1746-1750. Dreger, D. S. and D. V. Helmberger, 1993. Determination of Source Parameters at Regional Distances With Three-component Sparse Network Data. J. Geophys. Res. 98, 8107-8125. Frankel, A. and J. Vidale, 1992. A three-dimensional simulation of seismic waves in the Santa Clara valley California, from the Loma Prieta aftershocks, Bull. Seismol. Soc. Am., 82, 2045-2074. Fuis, G. S., T. Ryberg, N. J. Godfrey, D. A. Okaya and J. M. Murphy, 2001. Crustal structure and tectonics from Los Angeles basin to the Mojave Desert, Southern California, Geology, 29, 15-18. Gao, S., H. Liu, P. M. Davis and L. Knopoff, 1996. Localized amplification of seismic waves and correlation with damage due to the Northridge earthquake: Evidence for focusing in Santa Monica. Bull. Seismol. Soc. Am., 86, S209-S230. Gee, L, S. and T. H. Jordan, 1992. Generalized Seismological data functionals, Geophys. J. Int., Ill, 363-390. Godfrey, N. J., V. Langenheim, G. S. Fuis and D. A. Okaya, 2002. Lower crustal deformation beneath the central Transverse Ranges, Southern California: Results from the Los Angeles region seismic experiment, J. Geophys. Res., 107, 10.1029/2001JB000354. Graves, R., A. Pitarka and P. Somerville, 1998. Ground motion amplification in the Santa Monica area: effects of shallow basin edge structure, Bull. Seismol. Soc. Am., 88, 337-356. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 178 Hardebeck, J. L., P. M. Shearer, 2002. A new method for determining first-motion focal mechanisms. Bull. Seismol. Soc. Am. 92, 2264-2276. Hauksson, E., 2000. Crustal structure and seismicity distribution adjacent to the Pacific and North America plate boundary in southern California. J. Geophys. Res. 105, 13875-13903. Kawase, H., 1996. The cause of the damage belt in Kobe: “The basin-edge effect,” constructive interference of the direct 5-wave with the basin-induced diffracted/Rayleigh waves, Seismol. Res. Lett., 67, 25-34. Komatitsch, D., Q.-Y. Liu, J. Tromp, M. P. Siiss, C. Stidham and J. H. Shaw, 2004. Simulations of ground motion in the Los Angeles Basin based upon the spectral-element method. Bull. Seismol. Soc. Am., 94, 187-206. Komatitsch, D., and J. P. Vilotte, 1998. The spectral-element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures, Bull. Seism. Soc. Am. 88, no. 2,368 -392. Ludwig, W. J., J. E. Nafe and C. L. Drake, 1970. Seismic refraction, The Sea, A. E. Maxwell, (Editor) Vol 4, Wiley-Interscience, New York, 53-84. Magistrale, H., K. McLaughlin and S. Day, 1996. A geology based 3-D velocity model developed using nonlinear asymptotic coupling theory. Bull. Seism. Soc. Am., 86, 1161-1166. Magistrale, H., S. Day, R. W. Clayton and R. Graves, 2000. The SCEC Southern California Reference Three-dimensional Seismic Velocity Model Version 2. Bull. Seismol. Soc. Am. 90, S65-S76. Marquering, H., F. A. Dahlen and G. Nolet, 1999. Three-dimensional sensitivity kernels for finite-frequency traveltimes: the banana-doughnut paradox, Geophys. J. Int., 137, 805-815. Nafe, J. E. and C. L. Drake, 1960. Physical properties of marine sediments. The Sea, M. N. Hill, (Editor) Vol. 3. Interscience. New York, 794-815. Olsen, K. B., 1994. Simulation of three-dimensional wave propagation in the Salt Lake Basin, Ph.D. Thesis, University of Utah, Salt Lake City, Utah, 157p. Paige C. C. and M. A. Saunders, 1982. LSQR: An algorithm for sparse linear equations and sparse least squares, TOMS 8(1), 43-71. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 179 Pitarka, A., K. Irikura, T. Iwata, and H. Sekiguchi, 1998. Three-dimensional simulation of the near-fault ground motion for the 1995 Hyogo-ken Nanbu (Kobe), Japan, earthquake, Bull. Seism. Soc. Am., 88, 428-440. Suss, M. P. and J. H. Shaw, 2003, P-wave seismic velocity structure derived from sonic logs and industry reflection data in the Los Angeles basin, California, J. Geophys. Res., 108, NO. B3, 2170, ESE 13. Zhao, L., T. H. Jordan, K. B. Olsen and P. Chen, 2005. Frechet Kernels for Imaging regional earth structure based on three-dimensional reference models. Bull. Seismol. Soc. Am, submitted. Zhu, L., and H. Kanamori, 2000. Moho depth variation in Southern California from teleseismic receiver functions. J. Geophys. Res., 105, 2969-2980. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 7 180 Summary In this dissertation, we have developed a unified methodology for the analysis and inversion of regional seismic waveform data for seismic source parameters and seismic velocity structures. This new technique is based on the capability to compute synthetic seismograms that incorporate existing knowledge about earthquake sources and earth structures. It improves our source and structure models iteratively by assimilating new information acquired from ever-increasing observations. In regions such as Southern California, preliminary three- dimensional velocity models such as CVM3.0 are already available and efficient numerical algorithms such as finite-difference and spectral-element methods have been developed to simulate dynamic wave propagation. We utilized such capabilities and developed automated procedures to recover CMT solutions and to resolve finite source properties for small to moderate-sized earthquakes, and using CVM3.0 as reference model we have conducted a fully three-dimensional tomography for Los Angeles basin area. 7.1 Theoretical Development We have refined the procedures for measuring and interpreting generalized seismological data functionals (GSDFs) from regional waveforms (Chapter 2). Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 181 The GSDFs are frequency-dependent quantities with the dimensions of time that encode the phase and amplitude differences between an observed and a synthetic waveform. For a particular waveform at frequency co, we measure two quantities: a phase-delay time Stp{co) and an amplitude-reduction time dtq(o/). We have developed processing algorithms that can (a) use synthetic waveforms from ID (using F-K integration algorithm) or 3D models (using finite-difference or spectral-element method) as isolation filters, (b) account for complex interference effects on the seismogram, as well as the filter bias, and (c) automatically encode the metadata needed for computing the Frechet kernels that specify the sensitivity of StV A (co) to regional earth structure and finite source parameters. We have reformulated GSDF theory in terms of a discrete Fourier representation of the differential operator, D(a> ) = exp[i&j£)Tp(a>) - cvr)rq (co)|, that maps the synthetic waveform into the observed waveform. Gee & Jordan’s (1992) original analysis was based on a Taylor series expansion of D(co), which is suitable for narrow-band seismograms, such as low-pass filtered GDSN recordings. We have found, however, that if the frequency band is too wide, the expressions that relate the data functionals dtV A to the model functionals r)'rp ,q become unstable. In our revised procedures, we now express the observables < 5 f p,q at a particular frequency a > i as linear combinations of r)T p ,q over a range of nearby frequencies. The new formulation is much more stable for regional broadband seismograms, and we now employ it in our data processing and in the calculation of the Frechet kernels. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 182 The 3D sensitivity kernels derived from this new formulation provide accurate descriptions between the time-frequency-dependent observations and model parameter perturbations and can well characterize essential natures of the generation and propagation of seismic waves, especially for those phases that are formed through complex interactions with 3D heterogeneities and cannot be accurately described using geometric ray theory (Figure 2.7-2.8). We have also coded our new formulation in Scientific Python (SciPy), an open- source extension built upon the open-source Python scripting language, which is widely used in high-performance scientific computing community. 7.2 Model Assessment We have applied the GSDF methodology to obtain about 20,000 measurements of 5tP tq from more than 1000 body waves sampling 1073 source-receiver paths (Figure 3.1). Synthetic seismograms were computed for the SCEC CVM3.0 using K. Olsen’s finite-difference code. Synthetics were also computed using a frequency-wavenumber integration code for two types of ID structures: the laterally homogeneous SoCaL model and a set of path-averaged structures AID obtained by averaging the slowness of CVM3.0 between the source and receiver (Chapter 3). This data set provides good coverage of the Los Angeles region, and we have used it in our preliminary comparisons among the various ID and 3D models. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 183 Overall, the 3D model CVM3.0 provided a substantially better fit to the waveform data than either laterally homogeneous (SoCaL) or path-averaged ID models (AID) (Figure 3.7). For CVM3.0, the RMS scatter is somewhat larger for S phases and for P phases, primarily owing to signal-generated noises (e.g. interference by scattered waves). But the variance reductions relative to ID models for S phases are as large as the variance reduction for P phases, which suggests the empirical scaling relations used in CVM3.0 to derive S velocity from P velocity are good approximations. We concluded that both P and S models of CVM3.0 could be used as reference models in linearized tomography inversions. By analyzing spatial distribution of our phase-delay time measurements, we observed a negative correlation between Stp and the thickness of the basin. This suggests CVM3.0 may underestimate true wave velocity in the basin area. A more detailed study of the sensitivity kernels for some representative basin paths shows that contributions to the negative phase-delay time measurements mostly come from the basin area. This further confirms our conclusion about CVM3.0’s underestimating true wave velocity in the basin area. We have used the data to evaluate the magnitude of polarization anisotropy. As shown in Figure 3.12, the data at all frequencies show a small but consistent bias indicating that SV waves propagate slightly slower than SH waves. This may indicate a small amount of polarization anisotropy in the upper crust, perhaps Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 184 associated with the layering in sedimentary basins. It will probably not need to be included in the initial tomographic modeling. 7.3 Invert Data for Finite Source Parameters We have developed source-inversion procedures for recovering the “finite moment tensor” (FMT) from observations of dtV A using isolation filters from ID or 3D models; partial derivatives with respect to FMT parameters were evaluated numerically. The FMT removes fault-plane-ambiguity of the CMT and yields several additional parameters of seismological interest: the characteristic length, width and duration of the faulting, as well as the directivity vector of the fault slip. We resolved complete FMT parameters for 09/03/2002 Ml4.8 Yorba Linda earthquake (Chapter 4). An important accomplishment has been the computation of a collection of Receiver Green Tensors (RGT) for the 44 stations in Los Angeles region (Figure 5.1) from the CVM3.0 model. Using RGT and reciprocity, we can now compute 3D synthetics at these stations for an arbitrary source in LA region within seconds. Based on this RGT library, we have developed an automatic procedure for CMT inversion and fault-plane-ambiguity resolution. We calculated the sensitivity derivatives of the GSDF measurements to the FMT parameters using our refined CMT point source. The derivatives that describe the constructive and destructive interferences caused by the spatial extent of a seismic source can then be used to Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 185 test which of the two CMT nodal planes is the actual fault plane. For small to moderate-sized earthquakes with few aftershocks, our procedure provides a convenient and effective way to resolve fault-plane-ambiguity from low- frequency (<0.5 Hz) waves. Our results on the finite source properties of regional small to moderate-sized earthquakes revealed some very interesting tectonic structures in and around Los Angeles basin area. One example is the incipient faulting on the Puente Hills- Fontana seismicity trend located to the southwest of the right-lateral San Jacinto fault and the west-striking and north-dipping Cucamonga fault as shown on Figure 5.4 and 5.6. This seismicity trend is not associated with any mapped fault traces. We resolved fault-plane-ambiguity for 18 small events in this swarm of earthquakes and they all show right lateral source mechanisms conjugate to the southwest trend of the epicenter distribution (Figure 5.6) and parallel to San Jacinto fault. This Fontana seismicity trend stops quite abruptly at the Puente hills area. About 20 miles further to the southwest, in the Yorba Linda area, our study on the 03 September 2002 M l4.8 Yorba Linda earthquake and several surrounding smaller earthquakes shows left-lateral faulting conjugate to the nearby right-lateral Whittier fault. Whether there is any connection between these two seimicity clusters is interesting from a tectonic point of view. Recent studies show that a high proportion of north-south shortening in the Los Angeles metropolitan region is accommodated by conjugate strike-slip faulting and east- west escaping of crustal blocks [Walls, et a l, 1998; Chapter 4], The Puente Hills- Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 186 Fontana seismicity trend and the left-lateral faulting in Yorba Linda area might indicate a developing weak zone associated with such “escaping” deformation. By recovering finite source properties for small to moderate-sized earthquakes, our technique provides new information for constraining regional tectonics and fault structure, which are essential for seismic hazard analysis and mitigation. 7.4 Full 3D Tomography for LA Basin Most previous tomographic studies that derive 3D perturbations usually assume a ID starting-model, and they have generally made severe approximations to the kernel, either using a ID “path-average” approximation or a ray-theoretic 3D approximation. In this study, we derived 3D perturbation to a 3D starting-model (CVM3.0) using 3D kernels that account for all elastic-interaction effects. Our true 3D kernels (Chapter 2, Chapter 6) demonstrate that within a complex geological environment such as Los Angeles basin area, seismic waves interact with the structure in very complicated ways. Even at periods as short as 1 second, the sensitivity kernels rarely display properties expected by ray theory or the banana-doughnut phenomenon except for first-arriving P waves. The signals on a seismogram after the first arrival are usually difficult to identify and analyze, and our 3D kernels are very useful for visualizing the sources of energy that contribute to the signals and in locating potential problems in the reference model. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 187 Using our true 3D sensitivity kernels for our phase-delay time and amplitude- reduction time measurements, we have derived 3D perturbations to the reference model CVM3.0. Our tomography results indicate that CVM3.0 is too slow within the basin, suggesting those empirical rules used to establish the relation between P velocity and sediment age and depth might underestimate the true velocity in Los Angeles basin area. The application of our methodology on regional data from the earthquake-prone Los Angeles basin area proves the robustness, efficacy and generality of our theory. With the development of regional and local seismic networks throughout the world, the methodology developed in this study should be widely applicable to other areas of seismological interest. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 188 Bibliography Aagaard, B. T., J. F. Hall and T. H. Heaton, 2001. Characterization of near-source ground motions with earthquake simulations. Earthquake Spectra, 17(2), 177- 207. Ak9elik, V., J. Bielak, G. Biros, I. Epanomeritakis, A. Fernandez, O. Ghattas, E. J. Kim, J. Lopez, D. O'Hallaron, T. Tu, and J. Urbanic, High Resolution Forward and Inverse Earthquake Modeling on Terasacale Computers. Proceedings o f the ACM/IEEE SC2003 Conference November 15 -21, 2003 Phoenix, Arizona. Aki, K., Christoffersson, A. and Husebye, E.S., 1977. Determination of the three- dimensional seismic structure of the lithosphere, J. Geophys. Res., 82, 277-296. Aki, K and P. G. Richards, Quantitative Seismology, 2n d edition. University Science Books, Sausalito, CA, 2002. Backus, G., and M. Mulcahy, Moment tensors and other phenomenological descriptions of seismic sources, I, Continuous displacements. Geophys. J. Roy. Astron. Soc. 46, 341-361, 1976. Bao, H., J. Bielak, O. Ghattas, L. F. Kallivokas, D. R. O’Hallaron, J. R. Shewchuk and J. Xu, 1998. Large-scale simulation of elastic wave propagation in heterogeneous media on parallel computers, Comput. Methods Appl. Mech. Engrg., 152, 85-102. Ben-Menahem, A., Radiation of seismic surface-waves from moving sources. Bull. Seismol. Soc. Am. 51, 403-435, 1961. Bijwaard H. and W. Spakman, 2000. Nonlinear global P-wave tomography by iterated linearized inversion, Geophys. J. Int., 141, 71-82. Bowman, D. D. and G. C. P. King, Accelerating seismicity and stress accumulation before large earthquakes. Geophys. Res. Lett., 28, 4039-4042, 2001. Chen, P., L. Zhao and T. H. Jordan, A Quantitative Evaluation of SCEC Community Velocity Model Version 3.0, Eos, Trans. AGU, 84(46), Fall Meet. Suppl., Abstract S42H-02 (F1108), 2003. Chen, P., T. H. Jordan and L. Zhao, Finite moment tensors of 3 September 2002 Yorba Linda earthquake, Bull. Seismol. Soc. Am. 2005, in press. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 189 Chen, P., L. Zhao and T. H. Jordan, Generalized Seismological Data Functionals for Broadband Seismic Waveform Analysis and Inversion. Geophys. J. Int. 2005, in preparation. Cohee, B., and G. Beroza, Slip distribution of the 1992 Landers earthquake and its implications for earthquake source mechanics. Bull. Seismol. Soc. Am. 84, 692- 712, 1994. Courboulex F., A. Deschamps, M. Cattaneo, F. Costi, J. Deverchere, J. Virieux, P. Augliera, V. Lanza, D. Spallarossa, Source study and tectonic implications of the 1995 Ventimiglia (border of Italy and France) earthquake (M l = 4.7). Tectonophysics 290, 245—257, 1998. Crouch, J. K., and J. Suppe, 1993, Late Cenozoic tectonic evolution of the Los Angeles basin and inner California borderland: A model for core complex-like crustal extension: Geological Society o f America Bulletin, v. 105, p. 1415-1434. Dahlen, F. A., S.-H. Hung and G. Nolet, 2000. Frechet kernels for finite- frequency traveltimes - 1 . Theory. Geophys. J. Int. 141, 157-174. Dahlen, F. A. and A. M. Baig, Frechet kernels for body wave amplitude. Geophys. J. Int. (2002) 150, 440-466. Dahlen, F. A., S.-H. Hung and G. Nolet, Frechet kernels for finite-frequency traveltimes - 1 . Theory. Geophys. J. Int. (2000) 141, 157-174. Davis, T. L.J. Namson and R. F. Yerkes, 1989. A cross-section of the Los Angeles area: seismically active fold-and-thrust belt, the 1987 Whittier Narrows earthquake and earthquake hazard. J. Geophys. Res., 94, 9644-9664. Davis, P. M., J. L. Rubinstein, K. H. Liu, S. S. Gau and L. Knopoff, 2000. Northridge earthquake damage caused by geologic focusing of seismic waves. Science, 289, 1746-1750. Dreger, D. S. and D. V. Helmberger, Determination of Source Parameters at Regional Distances With Three-component Sparse Network Data. J. Geophys. Res. 98, 8107-8125, 1993. Dreger, D., and A. Kaverina, Seismic remote sensing for the earthquake source process and near-source strong shaking: A case study of the October 16, 1999 Hector Mine earthquake. Geophys. Res. Lett. 27, 1941-1944, 2000. Dziewonski, A. M., B. H. Hager, and R.J. O'Connell, 1977. Large-scale heterogeneities in the lower mantle, J. Geophys. Res., 82 (B2), 239-255. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 190 Dziewonski, A. M., T. A. Chou, and J. H. Woodhouse, Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. J. Geophys. Res. 86, 2825-2852, 1981. Edelman, S. H., A. Liu, and R. D. Hatcher, The Brevard in South Carolins and adjacent areas: An Alleghanian orogen-scal dextral shear zone reactivated by a thrust fault. J. Geology 95, 793-806, 1987. Frankel, A. and J. Vidale, 1992. A three-dimensional simulation of seismic waves in the Santa Clara valley California, from the Loma Prieta aftershocks, Bull. Seismol. Soc. Am., 82, 2045-2074. Fuis, G. S., T. Ryberg, N. J. Godfrey, D. A. Okaya and J. M. Murphy, 2001. Crustal structure and tectonics from Los Angeles basin to the Mojave Desert, Southern California, Geology, 29, 15-18. Gaherty, J.B., T. H. Jordan and L. S. Gee, Seismic structure of the upper mantle in a central Pacific corridor. J. Geophys. Res. (1996), vol 101, No. B10, 22291- 22309. Gao, S., H. Liu, P. M. Davis and L. Knopoff, 1996. Localized amplification of seismic waves and correlation with damage due to the Northridge earthquake: Evidence for focusing in Santa Monica. Bull. Seismol. Soc. Am., 86, S209-S230. Gee, L. S. and T. H. Jordan, 1992. Generalized Seismological data functionals, Geophys. J. Int., I l l , 363-390. Godfrey, N. J., V. Langenheim, G. S. Fuis and D. A. Okaya, 2002. Lower crustal deformation beneath the central Transverse Ranges, Southern California: Results from the Los Angeles region seismic experiment, J. Geophys. Res., 107, 10.1029/2001JB000354. Graves, R., A. Pitarka and P. Somerville, 1998. Ground motion amplification in the Santa Monica area: effects of shallow basin edge structure, Bull. Seismol. Soc. Am., 88, 337-356. Graves, R. W. and D. J. Wald, Resolution Analysis of Finite Fault Source Inversion Using ID and 3D Green’s Functions, Part I: Strong Motions. J. Geophys. Res. 106, 8745-8766. Hadley, D., and H. Kanamori, Seismic structure of the Transverse Ranges, California. Geol. Soc. Am. Bull. 88, 1469 - 1478, 1977. Hardebeck, J. L., P. M. Shearer, A new method for determining first-motion focal mechanisms. Bull. Seismol. Soc. Am. 92, 2264-2276, 2002. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 191 Harris, R. A., Introduction to special section: Stress triggers, stress shadows, and implications for seismic hazard. J. Geophys. Res., 103, 24347-24358, 1998. Hauksson, E., 1990. Earthquakes, faulting and stress in the Los Angeles basin. J. Geophys. Res., 95, 15365-15394. Hauksson, E., 2000. Crustal structure and seismicity distribution adjacent to the Pacific and North America plate boundary in southern California. J. Geophys. Res. 105, 13875-13903. Hauksson, E., K. Hutton, L. Jones and D. Given, The M4.8 Yorba Linda, Orange County earthquake of 3 September 2002: California Integrated Seismic Network Preliminary Report, 2p, 2002. Hauksson, E., P. Small, K. Hafner, R. Busby, R. Clayton, J. Goltz, T. Heaton, K. Hutton, H. Kanamori, J. Polet, D. Given, L. M. Jones, and D. Wald. Caltech/USGS element of TriNet 1997-2001. Seismol. Res. Lett. 72, 690-704, 2001 . Hellweg, M. and J. Boatwright, Mapping the rupture process of moderate earthquakes by inverting accelerograms, J. Geophys. Res. 104, 7319-7328, 1999. Ihmle, P. F., and T. H. Jordan, Teleseismic Search for Slow Precursors to Large Earthquakes, Science (1994), 266: (5190) 1547-1551, Dec 2. Jordan, T.H. and K.A. Sverdrup. Teleseismic location techniques and their application to earthquake clusters in the south-central Pacific, Bull. Seismol. Soc. Am. 71, 1105-1130, 1981. Katzman, R., L. Zhao, T. H. Jordan, High-resolution, two-dimensional vertical tomography of the central Pacific mantle using ScS reverberations and frequency- dependent travel times. J. Geophys. Res. (1998), vol 103, No. B8, 17933-17971. Kawase, H., 1996. The cause of the damage belt in Kobe: “The basin-edge effect,” constructive interference of the direct 5-wave with the basin-induced diffracted/Rayleigh waves, Seismol. Res. Lett., 67, 25-34. Kikuchi, M. and H. Kanamori, Inversion of complex body waves -III, Bull. Seismol. Soc. Am., 81, 2335-2350, 1991. Kohler, M. D., H. Magistrale, R. W. Clayton, Mantle heterogeneities and the SCEC reference three-dimensional seismic velocity model version 3, Bull. Seismol. Soc. Am. 93, 757-774, 2003. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 192 Komatitsch, D., Q. Liu, J. Tromp, P. Suss, C. Stidham and J. H. Shaw, Simulations of ground motion in the Los Angeles basin based upon spectral- element method, Bull. Seismol. Soc. Am. 94, 187-206, 2004 Komatitsch, D., and J. P. Vilotte, 1998. The spectral-element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures, Bull. Seism. Soc. Am. 88, no. 2,368 -392. Liu, Q.-Y., J. Polet, D. Komatitcsh and J. Tromp, Spectral-element moment tensor inversions for earthquakes in Southern California. Bull. Seismol. Soc. Am., 94, 1748-1761, 2004. Ludwig, W. J., J. E. Nafe and C. L. Drake, 1970. Seismic refraction, The Sea, A. E. Maxwell, (Editor) Vol 4, Wiley-Interscience, New York, 53-84. Luo, Y. and G. T. Schuster, Wave equation travel-time inversion. Geophysics (1991), 56, 645-653. Luyendyk, B. P. & J. S. Homafius, 1987. Neogene crustal rotations, fault slip and basin development in southern California, in Cenozoic basin development of coastal California, edited by R. V. Ingersoll & W. G. Ernst, pp 259-283, Prentice- Hall, Old Tappan, New Jersey. Magistrale, H., K. McLaughlin and S. Day, 1996. A geology based 3-D velocity model developed using nonlinear asymptotic coupling theory. Bull. Seism. Soc. Am., 86, 1161-1166. Magistrale, H., S. Day, R. W. Clayton and R. Graves, 2000. The SCEC Southern California Reference Three-dimensional Seismic Velocity Model Version 2. Bull. Seismol. Soc. Am. 90, S65-S76. Marquering, H., F. A. Dahlen and G. Nolet, 1999. Three-dimensional sensitivity kernels for finite-frequency traveltimes: the banana-doughnut paradox, Geophys. J. Int., 137, 805-815. Masters, G., S. Johnson, G. Laske and H. Bolton, A shear-velocity model of the mantle, Phil. Trans. R. Soc. Lond. (1996), A354, 1385-1411. McGuire, J., Estimating finite source properties of small earthquake ruptures, Bull. Seismol. Soc. Am. 94, no 2, 377-393, 2004. McGuire, J. J., L. Zhao, and T. H. Jordan, Teleseismic inversion for the second- degree moments of earthquakes. Geophys. J. Int. 145, 661-678, 2001. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 193 McGuire, J. J., L. Zhao, and T. H. Jordan, Predominance of unilateral rupture for a global catalog of earthquakes. Bull. Seismol. Soc. Am. 92, 3309-3317, 2002. Mori, J., Rupture directivity and slip distribution of the M 4.3 foreshock to the 1992 Joshua Tree earthquake, Southern California. Bull. Seismol. Soc. Am. 86, 805-810, 1996. Mori, J. and S. Hartzell, Source inversion of the 1988 upland, California, earthquake: determination of a fault plane for a small event, Bull. Seismol. Soc. Am. 80, 507-518, 1990. Mori, J., H. Kanamori, J. Davis, E. Hauksson, R. Clayton, T. Heaton, L. Jones, A. Shakal, R. Porcella, Major improvements in progress for southern California earthquake monitoring, Eos 79, 217-221, 1998. Nafe, J. E. and C. L. Drake, 1960. Physical properties of marine sediments. The Sea, M. N. Hill, (Editor) Vol. 3. Interscience. New York, 794-815. Nolet, G., Partitioned waveform inversion and two-dimensional structure under the network of autonomously recording seismographs. J. Geophys. Res. (1990) vol. 95. No. B6. 8499-8512. Olsen, K. B., 1994. Simulation of three-dimensional wave propagation in the Salt Lake Basin, Ph.D. Thesis, University of Utah, Salt Lake City, Utah, 157p. Okada, T., N. Umino, Y. Ito, T. Matsuzawa, A. Hasegawa and M. Kamiyama, Source Processes of 15 September 1998 M5.0 Sendai, Northeastern Japan, Earthquake and Its M3.8 Foreshock by Waveform Inversion, Bull. Seismol. Soc. Am. 91, 1607-1618, 2001. Paige C. C. and M. A. Saunders, 1982. LSQR: An algorithm for sparse linear equations and sparse least squares, TOMS 8(1), 43-71. Pasyanos, M. E., D. S. Dreger, and B. Romanowicz, Towards real-time determination of regional moment tensors. Bull. Seismol. Soc. Am. 86, 1255-1269, 1996. Pitarka, A., K. Irikura, T. Iwata, and H. Sekiguchi, 1998. Three-dimensional simulation of the near-fault ground motion for the 1995 Hyogo-ken Nanbu (Kobe), Japan, earthquake, Bull. Seism. Soc. Am., 88, 428-440. Press, W. H., S. A. Teukolsky, W. T. Vetterling and B. R. Flannery, Numerical Recipes in Fortran 77 - the art of scientific computing. Cambridge University Press, 2001. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 194 Royden. L., 1984. Thin-skinned extension of Vienna basin [abs]. Amer. Assoc. Petroleum Geologists, Bull. 68, 523. Seeber, L., and J. G. Armbruster, The San Andreas fault system through the Transverse Ranges as illuminated by earthquakes. J. Geophys. Res. 100, 8285- 8310,1995. Shaw, J. H. and P. Shearer, 1999. An elusive blind-thrust fault beneath metropolitan Los Angeles, Science, 283, 1516-1518. Silver, P., and T. H. Jordan, Total-moment spectra of fourteen large earthquakes. J. Geophys. Res. 88, 3273-3293, 1983. Sipkin, S. A. and T. H. Jordan, Lateral heterogeneity of the upper mantle determined from the travel times of multiple ScS. J. Geophys. Res., (1976) 81, 6307-6320. Su, W.-J. and A. M. Dziewonski, On the scale of mantle heterogeneity, Phys. Earth Planet. Inter. (1992), 74, 29-54. Suss, M. P. and J. H. Shaw, 2003, P-wave seismic velocity structure derived from sonic logs and industry reflection data in the Los Angeles basin, California, J. Geophys. Res., 108, NO. B3, 2170, ESE 13. Tarantola, A., A strategy for nonlinear elastic inversion of seismic reflection data, Geophysics (1986), 51, 1893-1903. Tromp, J., C. Tape and Q. Liu, Seismic tomography, adjoint methods, time reversal, and banana-doughnut kernels. Geophys. J. Int. (2005) 160, 195-216. Tukey, J. W., The collected works of John W. Tukey. Wadsworth, Belmont, CA., 1984. Vandenberghe, L. and S. Boyd, Semidefinite programming. SIAM Rev. 38, 49-95, 1996. Wald, D. J., and T. H. Heaton, Spatial and temporal distribution of slip for the 1992 Landers, California earthquake. Bull. Seismol. Soc. Am. 84, 668-691, 1994. Wald, D. J., T. H. Heaton and K. W. Hudnut, The slip history of the 1994 Northridge, California, earthquake determined from strong ground motion, teleseismic, GPS, and leveling data Bull. Seismol. Soc. Am. 86, S49-S70, 1996. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 195 Walls, C., T. Rockwell, K. Mueller, Y. Bock, S. Williams, J. Pfanner, J. Dolan and P. Fang, Escape tectonics in the Los Angeles metropolitan region and implications for seismic risk. Nature 394, 356—360, 1998. Webb, T. H., and H. Kanamori, Earthquake focal mechanisms in the eastern Transverse Ranges and San Emigdio Mountains, Southern California, and evidence for a regional decollement. Bull. Seismol. Soc. Am. 75, 737-757, 1985. Wright, T. L., 1991. Structural geology and tectonic evolution of the Los Angeles basin, California. In Active Margin Basins, edite by K. T. Biddle, vol. 52, pp. 35- 134. Am. Assoc. Pet. Geol. Memoir. Yeats, R. S., Tectonics of the San Gabriel Basin and surroundings, southern California, GSA Bulletin 116(9/10): 1158-1182, 2004. Zhao, L. and T. H. Jordan, Sensitivity of frequency-dependent traveltimes to laterally heterogeneous, anisotropic Earth structure, Geophys. J. Int., 133, 683- 704, 1998. Zhao, L., T. H. Jordan and C. H. Chapman, 2000. Three-dimensional Frechet differential kernels for seismic delay times. Geophys. J. Int., 141, 558-576. Zhao, L., T. H. Jordan, K. B. Olsen and P. Chen, 2005. Frechet Kernels for Imaging regional earth structure based on three-dimensional reference models. Bull. Seismol. Soc. Am, submitted. Zhu, L. and D. Helmberger, Advancement in source estimation techniques using broadband regional seismograms, Bull. Seismol. Soc. Am., 86, 1634-1641, 1996. Zhu, L., and H. Kanamori, 2000. Moho depth variation in Southern California from teleseismic receiver functions. J. Geophys. Res., 105, 2969-2980. Zhu, L., and L. A. Rivera, A note on the dynamic and static displacements from a point source in multilayered media Geophys. J. Int. 148, 619-627, 2003. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Holocene sedimentation in the southern Gulf of California and its climatic implications

PDF

A tectonic model for the formation of the gridded plains on Guinevere Planitia, Venus: Implications for the thickness of the elastic lithosphere

PDF

Evolutionary paleoecology and taphonomy of the earliest animals: Evidence from the Neoproterozoic and Cambrian of southwest China

PDF

Characterization of geochemical and lithologic variations in Milankovitch cycles: Green River Formation, Wyoming

PDF

A seismic source inversion study and non-stationary nonlinear seismic signal processing

PDF

Integrated geochemical and hydrodynamic modeling of San Diego Bay, California

PDF

Cyclostratigraphy and chronology of the Albian stage (Piobbico core, Italy)

PDF

Modeling and imaging asperities on a fault plane and characterizing spatial and temporal patterns of precursory seismicity

PDF

Numerical modeling of long period events at the Piton de la Fournaise by the use of object-oriented scientific computing

PDF

Investigation into the tectonic significance of the along strike variations of the Peninsular Ranges batholith, southern and Baja California

PDF

A cluster analysis of the differences between expert and novice counselors based on real time training

PDF

Applying aggregate-level traffic control algorithms to improve network robustness

PDF

Fractionation of nitrogen isotopes during early diagenesis

PDF

An adaptive soft classification model: Content-based similarity queries and beyond

PDF

A phylogenetic analysis of oological characters: A case study of saurischian dinosaur relationships and avian evolution

PDF

Automatically and accurately conflating road vector data, street maps and orthoimagery

PDF

Fourier grain-shape analysis of quartz sand from the Santa Monica Bay Littoral Cell, Southern California

PDF

Biotic recovery from the end-Permian mass extinction: Analysis of biofabric trends in the Lower Triassic Virgin Limestone, southern Nevada

PDF

A mutational analysis on monoamine oxidase

PDF

Expression and functional studies of the Notch signaling pathway in feather development

Asset Metadata

Creator Chen, Po (author)

Core Title A unified methodology for seismic waveform analysis and inversion

School Graduate School

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)

Advisor Jordan, Tom (committee chair)

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c16-397015

Unique identifier UC11336523

Identifier 3196790.pdf (filename),usctheses-c16-397015 (legacy record id)

Legacy Identifier 3196790.pdf

Dmrecord 397015

Document Type Dissertation

Rights Chen, Po

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