University of Southern California Dissertations and Theses

Applying automated techniques to large seismic datasets for systematic analyses of phases, source, and structure

(USC Thesis Other)

Applying automated techniques to large seismic datasets for systematic analyses of phases, source, and structure

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Applying automated techniques to large seismic datasets for systematic analyses of phases, source, and structure Zachary Elias Ross A Dissertation Presented to the Faculty of the University of Southern California Graduate School In partial fulfillment of requirements for the degree of Doctor of Philosophy in Geological Sciences August 2016 1 Acknowledgements I would like to take this opportunity to thank all of my family, friends, and colleagues who have helped me get to where I am today. My parents have been incredibly supportive throughout all of my years of schooling and helped in guiding me towards science from an early age. My wife, Melissa, is always a source of positive energy in my life and has been very understanding of my career. I am extremely fortunate to have such a wonderful family that always encouraged me to pursue this dream without regrets. My Ph.D. advisor, Professor Yehuda Ben-Zion, always found time for me and my research in the middle of his countless other obligations. Because of his generosity, I was able to attend 15 different conferences during my time at USC, which allowed me to establish important relationships with scientists from around the world. I cannot begin to list the many things I have learned from him, but that he has taught me how to be a scientist. I am grateful that Tom Jordan, Charlie Sammis, and Ken Alexander took the time out of their busy schedules to help advance my career, having served on both my qualifying exam and dissertation defense committees. I would like to express my gratitude to the University of Southern California for affording me the privilege of earning my Ph.D. here. USC gave me everything I could possibly need to succeed. Lastly, I would also like to thank my research collaborators: Frank Vernon, Francis Wu, David Okaya, Malcolm White, Dimitri Zigone, Shiqing Xu, and Lupei Zhu. I have learned a great deal from all of you. 2 Table of Contents Abstract…………………………………………………………………………………….. 3 Introduction………………………………………………………………………………… 5 1. Automatic picking of direct P, S seismic phases and fault zone head waves…………… 9 2. An improved algorithm for real-time S-wave picking with application to the (augmented) ANZA network in southern California…………………………………... 36 3. An earthquake detection algorithm with pseudo probabilities of multiple indicators……………………………………………………………………….............. 56 4. An algorithm for automated identification of fault zone trapped waves………………... 67 5. Towards reliable automated estimates of earthquake source properties from body wave spectra……………………………………………………………………… 88 6. Analysis of earthquake body wave spectra for potency and magnitude values: Implications for magnitude scaling relations ………………………………………….112 7. Spatio-temporal variations of double-couple aftershock mechanisms and possible volumetric earthquake strain………………………………………….............123 8. Isotropic source terms of San Jacinto fault zone earthquakes based on waveform inversions with a generalized CAP method………………………………...140 Discussion………………………………………………………………………….............161 References………………………………………………………………………….............163 3 Abstract I develop a collection of automated techniques suitable for systematic processing of large seismic waveform data sets. The techniques allow for extracting event, phase, and source information from the data sets with the goal of constructing databases for downstream analyses. I propose algorithms for real-time automated picking of P- and S-wave arrivals that are capable of detecting earthquakes at the same time. The methods apply polarization filters to the 3-component data to enhance signal to noise ratios and use moving averages in tandem with a kurtosis-based detector to lock in on the phase arrival. The techniques are used to construct earthquake catalogs for Southern California from scratch by re-processing the original waveform archives. A method is developed for detecting earthquake signals using pseudo-probabilities, by combining different indicator variables together in a non-parametric manner. The method is tested on data recorded by the Southern California Seismic Network and is found to detect more than three times as many earthquakes than the network’s catalog for the examined period of time. Methods are proposed for automated identification of fault zone head and trapped waves. These are special waves that are observed exclusively in fault zones and carry high-resolution information on the nature of the geologic structures which generated them. The techniques are applied to thousands of events recorded at the Parkfield section of the San Andreas fault, the Landers rupture area, the San Jacinto fault zone, and the Hayward fault. The identified head waves are used to estimate the extent of the velocity contrast across the faults, and the identified trapped waves are used to estimate various properties of the damage zones. A general-purpose methodology is developed for calculating numerous source quantities from body-wave spectral ratios including stress/strain drop, radiated seismic energy, and the extent of directivity. The technique is demonstrated on a set of five earthquakes with M > 3.4 in Southern California on hundreds of stations. Roughly 30,000 seismograms are examined in the process with no user involvement. Four of the five earthquakes show evidence of unilateral directivity. A catalog of more than 13,000 earthquakes is constructed for the San Jacinto fault zone and used to study the scaling relationship between seismic potency and local magnitude. Moment magnitudes are calculated for more than 11,000 local earthquakes using P- and S-wave spectra. It is found that moment magnitudes are statistically larger than local magnitudes for the same set of events. Further, the b-value measured from moment magnitudes is 25% larger than the b-value measured from local magnitudes. 4 Full moment tensors are determined for a set of 7 earthquakes in the San Jacinto fault zone. Small, explosive type isotropic components are measured in all of the events. A set of tests are performed in which the station distribution and velocity model are perturbed to assess the robustness of these values. It is found that the CLVD components change sign frequently in this process and are likely the result of artifacts. The the isotropic components persist at being explosive in type for all of the events in more than 90 percent of the perturbation scenarios. Double couple focal mechanisms are studied in different regions of the Landers rupture zone. I find that the short term aftershocks in the edges of the Landers rupture are less aligned with the mainshock focal mechanism than in the central portion of the rupture. Over time, the statistical orientation of the average focal mechanism approaches the long-term background mechanism. It is shown that a possible explanation for these observations is small isotropic components present in the early aftershocks of the edge regions that occasionally flips polarities measured on nodal planes. 5 Introduction Observational seismology is a fundamental component of earthquake science. It has a rich history that now spans more than a century, and revolves around the extraction of information from recordings of earthquake ground motions. These records, which are commonly referred to as seismograms, are the accumulation of earthquake source, propagation, site, and instrument effects. As a result, seismograms contain a wealth of information that is of interest to researchers. Making a measurement from a seismogram is often complicated by the fact that many of these cumulative effects are unknown, or understood within only a narrow range of conditions. In the early days of seismology, considerable fundamental progress was achieved through empirical analyses. These works include the development of earthquake magnitude scales (e.g. Richter, 1935), observations on the frequency-magnitude statistics of earthquakes (e.g. Gutenberg and Richter, 1944), and the productivity-decay of aftershocks (e.g. Omori 1894). These types of analyses identified general features of seismicity with relatively few assumptions. As a result, it is perhaps no surprise that these features are still commonly studied at the present day. By the mid 1980’s, considerable progress had been made on the theory of seismic sources, which allowed for details about the source processes of individual earthquakes to studied in ways not previously possible. Aki (1966) made the first estimate of the seismic moment of an earthquake, which is a fundamental metric of the size of an earthquake. Gilbert (1971) proposed the concept of a seismic moment tensor characterized by the stress drop. Backus and Mulcahy (1976) and Backus (1977) further developed the theory of moment tensors and introduced the quantity stress glut. Brune (1970) and Aki (1967) introduced ideas of scale invariant rupture processes and key expectations for earthquake spectra under these conditions. Within a few years, these theories were being utilized to obtain more detailed estimates of the source processes of large earthquakes (e.g. Dziewonski et al. 1981; Silver and Jordan, 1983). Analyses of these types often had labor intensive components, such as manually identifying suitable data windows or the best filter band. It was around this time that the first digital seismometers were developed, which opened the doors for automatic processing of seismic data. Allen (1978; 1982) introduced the STA/LTA algorithm, which is the ratio of a short-term moving average to a long-term moving average. The STA/LTA algorithm allowed for earthquake phases to be picked automatically in real time so that events could be detected and located without user involvement. The term “picking” describes the 6 process by which the phase onset time is measured from the time series data. It became the de facto procedure employed by regional seismic networks around the world. Before this technique was developed, events generally needed to be detected, picked, and located by hand. STA/LTA dramatically enhanced the detection capabilities of regional seismic networks and defined a standard for earthquake detection research that exists to the present day. Additional techniques have since been studied for the purposes of detecting earthquakes and picking phase arrivals. These include auto-regressive methods (Sleeman and van Eck, 1999), neural networks (e.g. Gentili and Michelini, 2006; Wang and Teng, 1997), spectrogram pattern recognition (e.g. Joswig, 1995), damped predominant period (Hildyard et al., 2008) and higher order statistics (Saragiotis et al., 2002). Automated techniques are objective since they are defined by a sequence of algorithmic steps. They can often nearly duplicate what trained experts would do, but many still require human verification at some level. While robust techniques for picking P-waves have been around for decades, there has been less success with methods for picking S-wave arrivals. This is because the S-wave arrives after the P-wave and is generally mixed with P-wave energy; this often makes the arrival time difficult to identify even for an expert seismologist. As such, earthquake locations are determined primarily from P-wave information because the P to S pick ratio is usually much larger than 1 (e.g. Kurzon et al., 2014). For many years, the main value of these tools was to the permanent seismic networks, as the data storage, memory requirements, hardware limitations, and lack of software infrastructure often precluded researchers from fully utilizing their potential. Over the last thirty years, the availability of computing resources and access to quality instrumentation have increased considerably. Regional seismic networks, such as the Southern California Seismic Network (SCEDC, 2013), have shifted from only preserving small segments of ground motion records (as recently as 2008) to archiving continuous data at up to hundreds of samples per second, twenty-four hours a day. Additionally, in recent years there has been a trend by researchers toward developing software that is highly-portable (e.g. Beyreuther et al., 2010). General purpose seismological software, which often was written by hand during previous generations, is widely available today and is even provided by some commercial vendors. Relational databases have been created to facilitate working with large volumes of waveform data and the associated meta-data (e.g. Pavlis et al., 2004). Advances in telemetry have made it feasible for temporary instrument deployments to have real-time delivery of waveform data. These technological achievements have led to an explosion of of seismic data being recorded worldwide 7 at the present day. With all of the gains made on the hardware side, it may be surprising that a considerable number of earthquakes around the world are recorded each day but never detected or identified. For an event to be detected, some minimum number of stations must each find evidence that an event has likely occurred (e.g. Withers et al., 1998). Part of this is the result of simply not having perfect instrumental coverage. Another critical issue is that the current best tools have fundamental limitations in their ability to identify an earthquake signal in the continuous time series data. These limitations are often due to signal quality (e.g. Nippress et al., 2010) or the fact that earthquake seismograms change with event size, source-receiver distance, geologic structure, and many other factors. Transient signals often appear in the continuous data that look like earthquakes, but are generated from other sources. Making a detection typically involves a threshold to be met at each station or for an array at-large. Choosing thresholds is often a compromise between detecting more events and having fewer false detections. Analyses of large earthquake data sets often involve systematic averaging over large spatio- temporal domains to overcome data limitations. This occurs in many areas of seismology, including seismic hazard analyses and ground-motion prediction equations (e.g. Abrahamson and Silva, 2008), derivation of earthquake source properties (e.g. Shearer et al., 2006), and statistical seismology (e.g. Hauksson et al., 2012). Considering the state of instrumentation, theory, and computing resources available today, tools for automatic processing of seismic data may in fact be the last remaining piece necessary for systematic data mining of earthquakes to become a reality. In this dissertation, I introduce a set of automated techniques for earthquake detection, phase recognition and onset picking, and source property estimation. These techniques are applied to numerous seismic data sets collected in California for comprehensive analyses. The dissertation is composed of eight chapters: seven of which are published in peer-reviewed journals, while the remaining one is currently under review. In Chapter 1, I introduce methods for automated picking of P- and S-waves as well as fault zone head waves, which are used throughout the dissertation (Ross and Ben-Zion, 2014a). In Chapter 2, I propose a key modification to the S- wave picking algorithm and extensively test the method’s picking accuracy against thousands of handmade picks (Ross et al., 2016). In Chapter 3, I propose an earthquake detection algorithm that combines a set of time series indicator variables into a single “pseudo-probability” series in a non- parametric manner (Ross and Ben-Zion, 2014b). Chapter 4 develops an algorithm for automated 8 identification of fault zone trapped waves, which is applied to two dense linear seismometer arrays across the Landers rupture zone and the San Jacinto fault zone (Ross and Ben-Zion, 2015). In Chapter 5, I discuss a technique for estimating numerous earthquake source properties from spectral ratios—including strain drop, extent of directivity, and seismic energy—which is suitable for application to earthquakes with M > 3.5 with no user involvement (Ross and Ben-Zion, 2016). Chapter 6 develops a simple procedure for estimating the low-frequency asymptote of body wave spectra, and is used to calculate M w for more than 11,000 earthquakes in the San Jacinto fault zone (Ross et al., 2016). In Chapter 7, I analyze the spatio-temporal behavior of focal mechanisms of the Landers aftershocks, and find that aftershocks in the edge regions have different temporal behavior than in the center which may reflect non-double couple components (Ross and Ben-Zion, 2013). Lastly, in Chapter 8, I determine full moment tensor solutions for a set of seven earthquakes with M > 4.2 which occurred in the complex San Jacinto fault zone trifurcation area and test the non-double couple components for robustness (Ross et al., 2015). The developed algorithms provide the collective ability to form databases of seismological measurements from virtually unlimited amounts of seismic data. Some of these quantities were only able to be measured previously through labor intensive visual examination of waveform data. These types of databases are suitable for large-scale inversions and statistical analyses. With sufficient tools to automatically process large volumes of data, it not only becomes possible to perform detailed fault-specific analyses, but this may ultimately lead to experimental earthquake seismology. Given a clear understanding of the number of phase picks or earthquakes necessary to perform a particular type of study—such as tomography—it becomes possible to design a temporary instrument deployment such that the amount of desired information is obtained. 9 1. Automatic picking of direct P, S seismic phases and fault zone head waves (Ross and Ben-Zion, 2014a) 1.1 Background Systematic analyses of seismic data sets require algorithms for automatic event detection and phase picking. This is especially important with the vast amount of data that are now continuously recorded by various seismic networks. Analysts' inspections and corrections can refine the automatic results, but are unfeasible with the large modern data sets and also less objective if done alone. In the present study we attempt to develop reliable algorithms for automatic detection of seismic phases in the early portions of P and S waveforms recorded in close proximity to large faults and earthquake rupture zones. The algorithms process earthquake waveform data sets that have already been selected by event detection algorithms (e.g. Withers et al., 1998; Nippress et al., 2010; Ross and Ben-Zion, 2014). Much effort has been devoted to algorithms for automatic picking of direct P-waves which are often the first arrivals on seismograms. Proposed techniques and signals include the short term average/long term average (STA/LTA; Allen, 1978), envelope functions (Baer and Kradolfer, 1987), neural networks (e.g. Gentili and Michelini, 2006) and auto-regressive methods (e.g. Sleeman and van Eck, 1999; Rastin et al., 2013). Higher order statistics including kurtosis and skewness have also been used for automatic phase picking (Saragiotis et al., 2002). The kurtosis was found to be more accurate than the traditional methods by a number of studies (e.g. Langet et al., 2014; Nippress et al., 2010; Kuperkoch et al., 2010) since it is primarily sensitive to abrupt changes in the character of a time series. A similar statement holds for the skewness. Nippress et al. (2010) utilized two picking techniques in tandem to obtain more reliable P-wave picks. They used an STA/LTA function or the damped predominant period T pd function (Hildyard et al., 2008) to find an initial picking region, and then refined the pick using a kurtosis characteristic function (CF). Automatic picking of S-wave arrivals is more challenging due to contamination of the direct S-wave arrival by the P coda and converted phases. A number of studies have incorporated polarization analysis into their picking algorithms to distinguish between P- and S-waves. Baillard et al. (2014) used for this purpose the kurtosis function with a “dip-rectilinearity” function derived from the covariance matrix. Cichowicz (1993) calculated a set of polarization quantities from the covariance matrix, and constructed from them a CF used to pick phases. Kurzon et al. (2014) used 10 a singular value decomposition based on the algorithm of Rosenberger (2010) to filter seismograms into polarized P and S channels, and picked direct P and S phases on characteristic functions derived from them. In various situations in seismology, the direct P-wave is not the first arrival on the seismogram. For example, events recorded at regional and teleseismic distances often have first arriving head waves (Pn) that spend most of the travel path propagating (refracting) along the Moho with the faster velocity of the mantle (e.g. Aki and Richards, 2002). Head waves are associated with emergent arrivals on the slower side of the lithology interface, in contrast to the more impulsive body waves. Similar waves that refract along horizontal interfaces are encountered in exploration seismology and exploited for imaging purposes. Head waves may also refract along dipping or vertical velocity contrast (bimaterial) interfaces of major fault structures. Such waves generated by earthquakes located close to the bimaterial interface are termed fault zone head waves (FZHW), and they have opposite first motion polarity from the following direct P arrivals on the slower side of the fault (e.g., Ben-Zion, 1989; Ben-Zion and Aki, 1990). In the last few decades numerous dense seismic networks have been deployed around large faults and earthquake rupture zones, in an effort to obtain high resolution results on structure and source properties (e.g. Bakun and Lindh, 1985; Fletcher et al., 1987; Seeber et al., 2000). Seismic data recorded near faults allow the detection and use of FZHW in diverse tectonic settings including subduction zones (e.g., Fukao et al., 1983; Hori et al., 1985; Hong and Kennett, 2003) and large strike-slip faults (e.g., Ben-Zion and Malin, 1991; Hough et al., 1994; McGuire and Ben- Zion, 2005; Bulut et al., 2012). Although FZHW exist only in relatively narrow regions on the slower velocity media (see Eq. (6) in section 4), they are recorded routinely by near-fault seismic networks (e.g., Ben-Zion and Malin, 1991; Allam et al., 2014). Misidentification of FZHW as direct P arrivals can lead to systematic errors in derived velocity structures, event locations and focal mechanisms (e.g., McNally and McEvilly, 1977; Oppenheimer et al., 1988), while proper use of FZHW can increase the resolution of the derived results (e.g., Ben-Zion et al., 1992; Bennington et al., 2013). Ruptures along bimaterial interfaces have fundamentally different properties than classical ruptures in a homogenous solid (e.g., Weertman, 1980; Andrews and Ben-Zion, 1997; Ranjith and Rice, 2001; Ben-Zion, 2001). Significantly, bimaterial ruptures tend to evolve for various conditions to pulses with a preferred propagation direction and dynamic reduction of normal stress 11 leading to small frictional heat and possible tensile component of faulting (e.g., Ben-Zion and Huang, 2002; Shi and Ben-Zion, 2006; Ampuero and Ben-Zion, 2008; Brietzke et al., 2009). These and related observational studies (e.g., Rubin and Gillard, 2000; Lewis et al., 2005; Dor et al., 2008; Zaliapin and Ben-Zion, 2011; Lengline and Got, 2011) highlight the importance of detecting and imaging properties of fault bimaterial interfaces. FZHW provide the clearest evidence for the existence, and highest imaging resolution, of fault bimaterial interfaces. In the following sections we describe and illustrate a set of algorithms for automatic picking of P and S body waves and FZHW generated by local earthquakes and recorded near faults. The algorithms can be used as a single unit to pick all these phases or as a set of separate phase picking techniques. The picking of direct P and S phases utilizes polarization analysis to distinguish between these waves, while the picking of direct P and FZHW utilizes expected properties of these waves to make picks with a high degree of accuracy. We apply the techniques to waveforms generated by several thousand events at the Parkfield section of the San Andreas fault, as well as synthetic waveforms and data recorded near the Hayward fault. The algorithms separate P and S waveforms for a large portion of the data, and identify early FZHW preceding the direct P body waves for a significant fraction of the waveforms recorded on the slow sides of the San Andreas and Hayward faults. 1.2 Data and pre-processing The algorithm for picking the various phases of interest (direct P, S and fault zone head waves) is primarily illustrated and tested using seismograms recorded at 5 stations of the High Resolution Seismic Network (HRSN) at the Parkfield section of the San Andreas fault (Figure 1). The data are recorded by short-period borehole seismometers operating at depths ranging from 100-300 m. In addition, we use synthetic seismograms in a model consisting of a vertical low velocity zone between two quarter spaces (Ben-Zion and Aki, 1990; Ben-Zion, 1998) to verify the performance of the direct P and head wave pickers, and we test the algorithm further on data of the Northern California Earthquake Data Center recorded around the Hayward fault. The Parkfield section of the San Andreas fault separates generally faster granite rocks to the SW from slower Franciscan rock to the NE (e.g., Eberhart-Phillips and Michael, 1993; Thurber et al., 2006). The 5 HRSN stations utilized in this work were chosen because they were used in several previous studies of FZHW near Parkfield (Ben-Zion and Malin, 1991; Ben-Zion et al., 12 1992; Zhao et al., 2010). These studies demonstrated that the three stations on the NE side of the fault (MMNB, GHIB, and EADB) frequently record FZHW, and the two stations located on the SW side (VCAB and FROB) do not. In sections 3-5 we analyze seismograms generated by 2489 events that occurred in 2004 (grey dots in Fig. 1.1). Figure 1.1. The Parkfield section of the San Andreas Fault with five HRSN stations (blue triangles) and epicenters of ~2500 earthquakes (grey circles and colored stars) used in this study from the NCEDC catalog. The events marked by stars are utilized to illustrate various aspects of the algorithm. The town of Parkfield, CA, and surface traces of faults are indicated by black square and lines. Prior to the phase picking we perform basic pre-processing involving removing the mean from each seismic record and band-pass filtering the data to remove long period noise and high frequency glitches that are sometimes present in the data. For the HRSN Parkfield records we use a (single-pass) causal band-pass filter between 0.5-30 Hz; these values may require modifications in other applications based on data quality. A causal filter is important because 2-pass acausal filters can introduce ringing to the front of the seismogram which can reduce the accuracy of 13 picking the arrival of FZHW (Allam et al., 2014). We also note that since FZHW are emergent with lower frequencies than direct P-waves, they can be significantly affected by filters above 0.5 Hz. The instrument response is not removed from the seismograms. 1.3 Picking P and S Phases The algorithm for picking P and S body waves relies on polarization analysis to facilitate their separation within recorded seismograms. Kurzon et al. (2014) recently used a singular value decomposition technique to separate P from S waves. Here we use a covariance matrix for the polarization analysis, because we find that the phase separation is more effective for some recordings with this method. However, in most cases both techniques of polarization analysis produce similar separation of P and S phases. The covariance matrix for a finite sample of 3-component data can be computed as, ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ = ) , Cov ) , Cov ) , Cov ) , Cov ) , Cov ) , Cov ) , Cov ) , Cov ) , Cov Z (Z E (Z N (Z Z (E E (E N (E Z (N E (N N (N σ , (1) where Z is the vertical component, N is north-south component, E is the east-west component, and the covariance between any two components X and Y is given by ∑ = = M i i i y x M Y (X 1 1 ) , Cov , (2) with M being the number of samples. For separating the P and S phases, we compute the covariance matrix of the pre-processed seismograms for a 3 s sliding window. This duration has been found to be stable for use with covariance matrices by Baillard et al. (2014). The corresponding eigenvalues ( ) 3 2 1 λ λ λ ≥ ≥ and eigenvector matrix ( ) 3 2 1 , , u u u u = of σ can be used to calculate various aspects of particle motion (e.g. Jurkevics, 1988). The rectilinearity, r, or degree of linear polarization, is calculated as, ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + − = 1 3 2 2 1 λ λ λ r , (3) and has a range of [0,1]. For body waves, r should be approximately equal to 1. The apparent vertical incidence angle, φ, for the particle motion is defined as 11 ) cos( u = ϕ . (4) 14 For local earthquakes, ) cos(ϕ should be close to 1 for P-waves and close to 0 for S-waves. From these quantities, we construct two polarization filters, ) cos(ϕ r p= , (5a) )) cos( 1 ( ϕ − = r s . (5b) These filters are used to modulate the 3-component data (e.g. Rosenberger, 2010) by simple multiplication at each time step. Specifically, p is used to modulate the vertical component and s is used to modulate the two horizontal components. We refer to these modulated signals as polarized traces since they have suppressed energy with a particular polarization. This part of the procedure is similar to that used by Cichowicz (1993), who also employed a third (transverse energy) quantity involving the ratio of the horizontal vector magnitude to the total vector magnitude. An example of our P and S separation process is illustrated in Figure 1.2 for an event recorded at station MMNB (red star, Fig. 1.1). The polarization filters p and s calculated in this case are shown in Figures 1.2d and Figure 1.2e, respectively. It is clear that p is large during the P-wave and small during the S-wave and vice versa for s. The polarized versions of the traces are plotted in Figure 1.2f with the polarized Z component in red and the polarized N component in blue. The P-wave signal originally present in the N component has been almost entirely removed, and the S- wave signal is mostly removed from the Z component. This leads to a clear visual separation of the phases. Next we run detectors on the polarized traces to pick various phases automatically. To obtain an initial pick of the P-wave arrival, we calculate a STA/LTA characteristic function on the polarized vertical trace with STA and LTA windows of 1 s and 30 s, respectively. These parameters were chosen so that the STA/LTA function would have general applicability to source-receiver distances and event magnitudes common to local earthquakes, without being overly sensitive to short transient signals not arising from earthquakes. The initial P arrival pick is made when the STA/LTA reaches an ad hoc level of 5. This value was found in exploratory analyses of several data sets to provide a good balance between the numbers of correct and false detections. The initial pick is used to approximately locate the P-wave arrival, but must be refined further as there is some amount of delay time associated with the STA/LTA recognizing the P-arrival. This refinement stage is described in detail in the subsequent section where we also account for the possibility of FZHW. An example of the P-wave picking process discussed so far is illustrated in 15 Figure 1.3 for the same event used in Figure 1.2, along with results for the S-wave pick described next. Figures 1.3a-c show the unpolarized 3-component traces with P and S wave picks denoted by red and blue vertical lines, respectively. Figure 1.3d displays the STA/LTA function calculated on the P-polarized vertical trace, and Figures 1.3e-i are associated with the S-wave pick. Figure 1.2. Application of the polarization analysis to separate a seismogram into P and S phases for the event denoted by red star in Fig. 1.1 and recorded at station MMNB. a) North component velocity seismogram. b) East component velocity seismogram. c) Vertical component velocity seismogram. d) A p polarization filter (Eq. 5a) determined from the covariance matrix with high values during the P-wave and low values during the S-wave. e) A corresponding s polarization filter determined from the covariance matrix (Eq. 5b). f) A P-polarized vertical trace (red) plotted together with the S-polarized east-west trace (blue). To pick the S-wave arrival, we first calculate an STA/LTA characteristic function on each of the S-polarized horizontal traces (Figures 1.3e-f). The STA and LTA windows are 1 s and 30 s in duration, respectively. These values were chosen for general applicability to local earthquakes, and changing the durations somewhat around these values does not have an appreciable effect on the results. If the polarization filters are working well, the S-polarized horizontal traces will have minor P-wave energy remaining (as in Figure 1.2). In such cases, the STA/LTA detectors will 16 primarily recognize the S-wave and can be used to make an initial pick. This is done on the S- polarized N and E components independently by making a pick at the location of the maximum of the CF. As seen in Figures 1.3e-f, the S arrival is not clearly defined on the STA/LTA detectors, and we found that using the peak values works best at locating the initial picking region. Figure 1.3. Demonstration of the P and S picking algorithm for the event denoted by red star in Fig. 1.1 and recorded a MMNB. The red and blue bars denote the automatic P- and S-wave picks, respectively. a) North component velocity seismogram. b) East component velocity seismogram. c) Vertical component velocity seismogram. d) STA/LTA function calculated on P-polarized vertical component. e) STA/LTA function calculated on S-polarized N component. f) STA/LTA function calculated on S-polarized E component. g) Kurtosis function calculated on S-polarized N component. h) Kurtosis function calculated on S-polarized E component. i) Time derivative of kurtosis function calculated on S-polarized E component. The STA/LTA detector is very effective at recognizing changes in amplitude, but less so at accurately identifying phase arrivals. For this purpose, we refine the initial S picks with a moving kurtosis function 17 3 ) 1 ( ) ( 4 1 4 − − − = ∑ = x M i i s M x x K , (6) where M is the number of samples, } ,..., , { 2 1 M x x x x= is the finite length P-wave sequence, x is the sample mean, and s x is the standard deviation. The kurtosis function measures peakedness and is useful for identifying abrupt signal changes such as those produced in seismic records by earthquake phase arrivals. We calculate the kurtosis for a sliding window of 5 s duration on each of the S-polarized horizontal traces (Figures 1.3g-h), and then calculate the derivative of each (Figure 1.3i). We tested sliding window durations in the range of 1-10 s and found that 5 s typically was the most accurate, but changing the window length somewhat around this value does not have a significant effect on the results. For each component, we search for the peak of the derivative in a 0.5 s window around the initial STA/LTA S pick and tentatively move the pick to this location. The derivative is used over the kurtosis because it helps to localize the arrival better. Since the kurtosis function is slightly delayed in time relative to the motion onset, we look for a local minimum in the kurtosis function preceding the tentative pick location by no more than 0.25 s. If multiple minima are present in this time window, the one with the smallest value is taken. The 0.25 s value was chosen to allow a small amount of flexibility in refining the pick, without moving it too far. If one is found, the pick is updated to the minima’s location (blue vertical line), and if not the pick remains unchanged. In some cases, the P and S phases do not separate well and considerable energy for each phase remains in the polarized traces. This can lead to the P-wave arrival being recognized strongly by the STA/LTA detectors on the N and E components, which can cause picking problems. To suppress such cases we check that the P and S picks are more than 0.3 s apart and discard any S picks which do not satisfy this requirement. The separation of P and S phases and quality of the S picks are generally better at stations on the faster side of the fault, without the additional complexity of fault zone waves that are present at near-fault stations on the slow side. To demonstrate the accuracy of the technique in picking phase arrivals, the automatic picks should be compared quantitatively with manual ones. Toward this end, S-wave arrivals were picked manually for 196 events and the S-picking algorithm was run on the same set. These events were chosen because they have P, S, and FZHW present in each (the P pick accuracy is discussed separately in the next section where the topic is treated in detail). The S arrivals were picked on each of the horizontal components independently. The median absolute time difference between 18 automatic and manual picks of the S-waves was 0.17 s, with a standard deviation of 0.72 s. Approximately 80% of the automatic picks were within 0.5 s of the manual pick. When FZHW are present they are likely to be misidentified by the STA/LTA and kurtosis pickers as direct P arrivals. Indeed, a closer examination of the first arrival picked in Figure 1.3 indicates that it is a FZHW followed by a direct P arrival that is delayed by about 0.1 s (Figure 1.4). This problem was noted by the early observational papers on FZHW (e.g. Ben-Zion and Malin, 1991) and was demonstrated clearly by Allam et al. (2014) in the context of the Hayward fault; the P-wave picks of the Northern California Earthquake Data Center at near-fault stations on the slow side of the Hayward fault are commonly FZHW rather than direct P waves (Figs. 2 and 6 of Allam et al. 2014). Developing an automatic picking algorithm capable of distinguishing between these two phases can be very useful for processing seismic data recorded near large faults. This is described in the next section. Figure 1.4. Close-up of vertical component velocity and displacement seismograms for the data used in Figure 1.3. Dashed line indicates automatic P-wave pick, while solid line indicates analyst pick of the direct P-wave. The automatic picker locks on a FZHW instead of the direct P-wave that is delayed by ~0.1 s. 1.4 Picking Fault Zone Head Waves Fault zone head waves are emergent phases that propagate along fault bimaterial interfaces with first motion polarity of the faster velocity medium, and are radiated from the interface to the slow side of the fault (Ben-Zion, 1989, 1990). They are the first arrivals at stations on the slower 19 velocity medium located within a critical distance x c from the fault given by (Ben-Zion and Malin, 1991), ( ) ( ) f s c r x α α / cos tan 1 − ⋅ = , (7) where r is the distance the FZHW travels along the fault and f α , s α are the average P-wave velocities of the faster and slower media, respectively. Several characteristics distinguish direct P body waves from FZHW: they are more impulsive, have opposite first motion polarity, and are usually larger in amplitude. However, since the amplitude of the FZHW may be significantly larger than the ambient noise (e.g., McGuire and Ben-Zion, 2005; results below), they are routinely picked erroneously by automatic algorithms for detection of P body waves. Our algorithm for detection and picking of FZHW and direct P arrivals exploits the different characteristics of these phases, and the procedure is applied with several stages of quality control to minimize false detections. The algorithm for separating between and picking FZHW and direct P phases only needs one component seismograms and below we apply it to the vertical component of the unpolarized seismograms. This is because the polarization filters of the previous section can potentially change the character of the waveforms, leading to less accurate picks of FZHW and direct P arrivals. The initial stage of P-wave picking is used in our algorithm only to find the relevant portion of the seismogram for detailed analysis of direct and possible FZHW arrivals. The first stage of the FZHW picking algorithm involves making a pick on the vertical component for the earliest onset of the seismic motion over the noise. This pick should lock on the head wave arrival if one is present or the P-wave arrival otherwise. It is performed using an STA/LTA detector with an STA length of 0.1 s, LTA length of 10 s, and trigger level of 4. These values were found to produce good results with the Parkfield data (see next section). The STA/LTA detector is used at this stage because it is more sensitive to the emergent onset of the FZHW than the kurtosis function. An example of this process applied to the data of Figure 1.3 is shown in Figure 1.5. To demonstrate the process we use raw data from station MMNB on the slow side with no band-pass filtering. Figure 1.5a and 1.5b show the raw velocity and displacement seismograms for the used event, respectively. The STA/LTA function calculated on the velocity seismogram is shown in Figure 1.5c and the first-motion pick is indicated by a red vertical line. 20 Figure 1.5. Demonstration of the FZHW and direct P picking algorithm for the event denoted by red star in Fig. 1.1 and recorded at MMNB. The red and blue lines indicate the automatic FZHW and direct P picks, respectively. No filtering is used in the algorithm for this example. The time difference between the FZHW pick and P-wave pick is 0.068 s. a) Vertical component velocity seismogram. b) Vertical component displacement seismogram with 0.5 Hz high-pass filter for plotting purposes only. c) STA/LTA function calculated from a). d) Kurtosis function calculated from a). e) Skewness function calculated from a). Note the reversed polarities between the head and P-waves, indicated most clearly by the skewness. f) Time derivative of kurtosis function calculated from d). g) Time derivative of skewness function calculated from e). We now attempt to specifically pick a direct P-wave arrival. This is achieved using two additional pickers in unison to find the sharpest arrival in the early portion of the seismogram. The first of these pickers uses the kurtosis function (6), and the second picker uses a moving skewness function (Saragiotis et al., 2002), 3 1 3 ) 1 ( ) ( x M i i s M x x S − − = ∑ = , (8) where the various symbols are the same as in (6). The skewness measures statistical asymmetry and in particular indicates whether the sample is left-skewed or right-skewed. Both the kurtosis and skewness are quite sensitive to abrupt changes in the character of a time series, which should occur more strongly at the P-wave arrival than for an emergent FZHW. We calculate both of these 21 statistical functions with a 5 s sliding window (Figures 1.5d-e), and find the peak of each in the vicinity of the examined portion of the seismogram. The value of 5 s is used because it is half the duration of the LTA in the first-motion pick (using windows of 3-7 s lead to similar results). We then calculate the derivative of each function within one second to the left of the peak’s location, and take the position of the maximum derivative as the pick of the direct P-wave (Figures 1.5f-g). The derivative is used rather than the kurtosis/skewness functions themselves to indicate where these functions change most rapidly. The one second window restricts the pick to the vicinity of the seismogram onset; the precise value of the window does not have an appreciable effect on the results. There may be various circumstances not associated with FZHW leading to several distinct phases in the early portion of seismograms (e.g., free-surface reflections at borehole instruments, arrivals from different earthquakes, etc.). To distinguish such early phases from fault zone head waves we use the property that a FZHW has opposite first motion polarity from the direct P-wave (Ben-Zion, 1989, 1990). We include this constraint in the algorithm by examining the polarity of the calculated skewness function at both the first motion pick and the later skewness pick. The skewness function provides a reliable indicator of motion polarity that is insensitive to small fluctuations. In many cases it is more stable than measuring the polarities from the velocity or displacement seismograms, because it only recognizes polarity changes when the moving window fully changes its statistical asymmetry. The algorithm requires further that the skewness does not change polarity other than near the P-wave arrival. This is because a genuine FZHW should maintain the same polarity until the direct P-wave arrives. As an example, the skewness in Figure 1.5e has only a negative polarity during the portion associated with FZHW, while it has only a positive polarity during the P-wave. In some cases, the direct P pick may be a few time samples earlier than the polarity reversal in the skewness function, and we account for this by looking for a nearby zero crossing within half the rise time around the pick (here about 0.02 s). The half-rise time is approximated as the time between the pick (where the derivative is greatest) and the peak of the CF. If a polarity reversal is present, we measure the polarity immediately after the zero crossing. Otherwise, the polarity is measured at the direct P pick. We then require that the polarities of the tentative FZHW and direct P arrival are opposite. Next we calculate the time difference between each tentative pick of the direct P-wave and the first motion pick, and compare these values to basic theoretical expectations for FZHW. If no head 22 wave is present and there are no unusual early arrivals (e.g. prominent free-surface reflections), the differences between the picks should be near zero. For the example in Figure 1.5, the obtained time difference is around 0.07 s for both the kurtosis and skewness signals. In a simple model consisting of two different homogeneous quarter spaces (Ben-Zion, 1989), the FZHW travel time is x r t f s f H 2 2 − − − + = α α α , (9) where x is the perpendicular distance from the fault to the station and 2 2 x R r − = is the component of the source-receiver separation along the fault with R being the hypocentral distance. The travel time for the direct P-wave is s D R t α = . (10) From (9) and (10), we can calculate the time difference between FZHW and direct P arrivals as D H t t t − = Δ . (11) We use (11) in conjunction with assumed P wave velocities of the media across the fault to determine maximum allowed time differences between FZHW and direct P arrivals at different stations. For the Parkfield area we assume α f = 5.5 km/s and α s = 4.95 km/s, associated with a velocity contrast of 10%, which is at the upper end of imaging results for the region (e.g., Ben- Zion et al., 1992; Zhao et al., 2010; Bennington et al., 2013). A minimum allowed time difference is defined from the approximate length for two swings of a direct P-wave, estimated here to be 0.065 s. This value is an effective resolution limit of the algorithm where a FZHW is no longer distinguishable from two swings of a direct P-wave. Reducing or eliminating the minimum time would increase the number of detections, while increasing also the number of false FZHW detections. If source locations are unknown for data of interest, it is possible to simply replace the hypocentral distance with a maximum value expected for sources in the region. This will act to increase the values given by Eq. 11, leading generally to the same number of correct picks and also more erroneous picks. Similarly, when analyzing data of fault zone arrays or near-fault stations without precise information on the distance from the fault, it is best to set a fixed (small) distance from the fault for all stations (e.g. 0.25 km), which will increase the values estimated by Eq. 11. Applications to other regions, or efforts to maximize the possible detections, may require changing the assumed velocity values and minimum allowed time separation. 23 To continue with the next stage of the algorithm, we require the time difference for the direct P picks based on the kurtosis and skewness to be within the range associated with the predicted difference for the FZHW. To ensure that the direct P picks based on the kurtosis and skewness are in good agreement, we further require that both are locked on the same portion of the seismogram by allowing them to be separated in time by no more than 0.03 s. The value was chosen to correspond roughly to one swing of a P-wave. This restriction decreases the likelihood that the direct P picks are in the allowable range by chance, and additionally gives confidence to the picks as recognizing a genuine phase arrival. If the picks are separated in time by more than 0.03 s, it essentially indicates that the arrival is not well defined. The time value of 0.03 s was chosen based on inspection of the rise time in the kurtosis/skewness functions for numerous local earthquakes, and testing indicates that a majority of events will fall within this range. This is primarily used to suppress outliers at this stage. If all requirements have been satisfied, the first-motion pick is identified as a FZHW and the direct P pick is taken as the average of the kurtosis and skewness picks (blue vertical line, Fig. 1.5). Otherwise, the first motion pick is identified as a direct P pick, and the later P picks are discarded. The algorithm up to this point focused on robust identification of phases using finite-size time windows, leading to picks that are slightly later in time than the arrival onsets. In the final stage of the algorithm we attempt to refine the picks to lock on more sharply defined time points. For the first-motion pick, we start at its current value and move backward in time until an STA/LTA value of 2 is reached, at which point the pick is moved to that location. This is an attempt to get the pick as close to its departure from the noise level as possible, as STA/LTA is effectively a measure of signal-to-noise ratio. If the first-motion phase was flagged as a FZHW candidate, we also refine the later direct P pick by returning to the original picks made on the kurtosis and skewness. These picks were made at the peak of the derivative of the respective CF, which is roughly half-way between when the CF begins to increase and the maximum of the CF. The time between where the derivative is maximum and the CF is maximum is approximated as half the rise time. We search for a local minimum on the kurtosis function that precedes the skewness polarity reversal by no more than the rise time. If one is found, the kurtosis pick is revised to this value and otherwise it remains unchanged. The skewness may have a local minimum or maximum preceding the polarity reversal as the polarity of the arrival is variable. Therefore a local extremum is searched for on the skewness function that precedes the polarity reversal by no more than the rise time. If one is found, 24 the pick is updated to the location of the extremum; otherwise the pick remains unchanged. For the kurtosis function, if multiple minima are present then the one with the smallest value is taken. The skewness function has two possible cases because the polarity can be both positive and negative. If the polarity is positive, multiple minima are checked for and if present the one with the smallest value is taken. Similarly, if the polarity is negative, multiple maxima are checked for and if present the largest is taken. The final direct P pick is set to the average of the new kurtosis and skewness picks, and at this point the algorithm comes to an end. The various stages of the algorithm are summarized by a flow chart in Figure 1.6. Figure 1.6. A flow chart describing the algorithm for identification and picking of P, S, and FZHW phases. To verify that the algorithm correctly identifies and correctly picks both FZHW and direct P arrivals, we test it on synthetic seismograms. The seismograms are calculated using the analytical solution of Ben-Zion and Aki (1990) for the scalar wave equation with a line dislocation source in a structure consisting of a 100 m wide vertical low velocity zone between two different quarter spaces (Figure 1.7). The results correspond to acoustic P or SH waves depending on the assumed 25 material properties. Here we calculate P waveforms using for the 3 media, denoted by subscripts 1, 2, 3 from left to right, a common mass density of 2500 kg/m 3 and the following wave velocities and attenuation factors: α 1 = 5500 km/s, α 2 = 3700 km/s, α 3 = 5000 km/s, Q 1 = Q 3 = 1000, and Q 2 = 30. The source is located at a depth of 10 km on the interface between the left quarter space and fault zone layer (Figure 1.7). In the examples below, velocity and displacement seismograms are calculated at two receivers located on the opposite sides of the fault zone, at epicentral distances of 200 m from the overall (left) velocity contrast interface. A small amount of Gaussian noise is added to the synthetic seismograms to produce more realistic records with small fluctuations mimicking ambient noise between the main arrivals. Figure 1.7. A vertical low velocity layer between two different quarter spaces with a line dislocation source on the fault (circle). Synthetic seismograms are calculated at the two receivers (triangles) on the opposite sides of the fault with the solution of Ben-Zion and Aki (1990). See text for media parameters. Figures 1.8a and 1.8b show synthetic velocity and displacement seismograms calculated at a receiver on the faster quarter-space 200 m from the left interface. In this case the basic phases present in the seismograms are a direct P-wave and trapped-like waves resulting from the low velocity layer. The algorithm correctly makes only a direct P pick (dashed line) within 0.01 s of the true arrival time. The three CF of the algorithm based on STA/LTA, kurtosis and skewness are given in Figures 1.8c-e, respectively. Figure 1.9 shows corresponding results associated with seismograms at a receiver on slower quarter-space 200 m from the left interface. In this case the first ballistic phase is a FZHW, followed by a direct P-wave, fault zone reflected phases and 26 trapped-like waves. The algorithm correctly identifies both the FZHW and direct P-wave and makes picks (red and blue vertical bars) that are within 0.01 s of the true arrivals. We tested the algorithm on numerous additional synthetic seismograms, using cases that have an additional transition layer between the two quarter spaces (Ben-Zion, 1998) and obtained similar reliable results as those shown in Figures 1.8 and 1.9. Figure 1.8. Testing the FZHW and direct P pickers on synthetic seismograms on the faster side of the fault, where only direct P and trapped waves are present. All picks made by the algorithm lock on the direct P-wave (dashed line) and no FZHW is detected. a) Synthetic velocity seismogram. b) Displacement seismogram. c) STA/LTA function calculated on a). d) Kurtosis function calculated on a). e) Skewness function calculated on a). 1.5 Systematic Tests on Parkfield Data To perform stronger tests on the performance of the algorithm we apply it to seismograms of ~2500 earthquakes recorded by 5 HRSN stations at the Parkfield section of the San Andreas fault (Figure 1.1). The San Andreas fault near Parkfield provides a good testing ground, because it has both a prominent velocity contrast and local complexities that may include a reversal in the sense of the contrast (e.g., Eberhart-Phillips and Michael, 1993; Thurber et al., 2006). It is known from Ben-Zion and Malin (1991) and Zhao et al. (2010) that stations MMNB and EADB on the generally slower NE block record commonly FZHW over a wide range of distances, while at stations FROB and VCAB on the generally faster SW block FZHW have not been observed. Station GHIB on the NE block is located above the hypocenter location of the 2004 M6 event associated with a local reversal of the velocity contrast (Bennington et al., 2013), and does not record head waves from events in that portion of the seismogenic zone (Zhao et al. 2010). 27 Figure 1.9. Testing the FZHW and direct P pickers on synthetic seismograms on the slower side of the fault, where a FZHW is present in addition to direct P and trapped waves. The red and blue lines indicate the picked FZHW and direct P-wave, respectively. a) Synthetic velocity seismogram. b) Displacement seismogram. c) STA/LTA function calculated on a). d) Kurtosis function calculated on a). e) Skewness function calculated on a). We find that considerably more FZHW are identified and picked by the algorithm on the generally slow NE side of the fault than the opposite SW side. In particular, 10.8% of the examined events recorded at MMNB are flagged as having candidates for FZHW, while the detection level at VCAB is only 0.2% (54 times fewer). Figure 1.5 illustrated the automatic detection and picking of FZHW at MMNB using data generated by an event close to and below the station (red star, Fig, 1.1). Figure 1.10 provides another example of automatic detection and picking of FZHW at MMNB using data generated by an event to the NW of the station (yellow star, Fig. 1.1). The seismograms recorded at VCAB and FROB for the events used in Figures 1.5 and 1.10, are shown in Figures 1.11 and 1.12, respectively. As seen in these examples, the recordings at FROB and VCAB typically have a single immediate peak in the kurtosis and skewness functions, while the recordings at MMNB often have two peaks that are about 0.1 s apart or more. We also note that these events have right lateral strike-slip focal mechanisms (Thurber et al., 2006), but they produce the same first motion polarity on both sides of the fault as expected for FZHW. 28 Figure 1.10. Illustration of the FZHW and direct P picking algorithm for the event denoted by yellow star in Fig. 1.1 and recorded at MMNB. The red and blue lines indicate the automatic FZHW and direct P picks, respectively. a) Vertical component velocity seismogram. b) Vertical component displacement seismogram c) STA/LTA function calculated from a). d) Kurtosis function calculated from a). e) Skewness function calculated from a). f) Time derivative of kurtosis function calculated from d). g) Time derivative of skewness function calculated from e). Figure 1.13 summarizes the detection of FZHW at the 5 employed HRSN stations, using corresponding colors for stations and FZHW detections, in a map view (Fig. 1.13a) and vertical cross section (Fig. 1.13b). For station FROB on the generally fast SW side of the fault, the algorithm flags candidate FZHW in 1.6% of the events. For stations GHIB and EADB on the nominally slow NE side, the FZHW detection rates are 9.1% and 4.6%, respectively. We recall that the algorithm parameters (e.g., minimum and maximum time separation between FZHW and direct P arrivals) are chosen to reduce the number of false detections rather than maximize the detection rate. As seen in Figure 1.13, the stations NE of the fault detect candidate FZHW (red, orange, and yellow circles) along most of the fault, with notable exceptions in the region below and around station GHIB having a local reversal of the velocity contrast (Thurber et al., 2006; Bennington et al., 2013) and a region of structural complexity to the NW of station EADB. The locations of automatic FZHW detections within the seismicity patches are very consistent with the manual results of Zhao et al. (2010). The blue and green symbols, signifying FZHW detections by 29 stations VCAB and FROB on the SW of the fault, are few in number and are not found to have FZHW upon closer inspection. Figure 1.14 provides examples of false detections at these stations as well as at the 3 stations NE of the fault. As mentioned, the false detections are associated with spurious phases such as reflections from the free surface or other early arrivals that are reminiscent of FZHW. They may be identified by checking that the first motion polarities at different stations are consistent with expectations for direct body waves (i.e., opposite polarities across the fault) and lack of other characteristic features of FZHW. Figure 1.11. Corresponding results to Figure 1.5 generated by the event denoted by red star in Fig. 1.1 and recorded at station VCAB. Only a direct P-wave is picked by the algorithm (blue line). a) Vertical component velocity seismogram. b) Vertical component displacement seismogram c) STA/LTA function calculated from a). d) Kurtosis function calculated from a). Note the single immediate peak after the motion onset. e) Skewness function calculated from a). As with the S picks, we test the accuracy of the technique in picking P and FZHW arrivals by comparing the automatic picks with manual ones. The same data set of 196 events was used, and each event was manually selected to ensure it had a FZHW as the first arrival. P and FZHW arrivals were picked manually on the vertical component. The median absolute time difference between automatic and manual FZHW picks was found to be 0.016 s, with a standard deviation of 0.035 s. For the P picks, the median absolute time difference between automatic and manual picks was found to be 0.004 s with a standard deviation of 0.023 s. 30 Figure 1.12. Corresponding results to Figure 1.9 generated by the event denoted by yellow star in Fig. 1.1 and recorded at station FROB. Only a direct P-wave is picked by the algorithm (blue line).b) Vertical component displacement seismogram c) STA/LTA function calculated from a). d) Kurtosis function calculated from a). Note the single immediate peak after the motion onset. e) Skewness function calculated from a). 1.6 Discussion Reliable techniques for automatic identification and picking of seismic phases are becoming increasingly important for analyzing large seismic data sets. We have developed a set of algorithms in this study which have the ability to accurately identify and pick arrivals of direct P and S body waves, as well as early arriving fault zone head waves (Figure 1.6). The algorithms utilize characteristic functions based on STA/LTA, kurtosis, and skewness of waveforms in moving time windows. The S-wave picking algorithm uses polarization analysis of 3-component waveforms to considerably reduce P-wave energy from the seismograms, so that the employed detectors lock on the S arrival (Figures 1.2, 1.3). The pick of the S phase itself is made in two stages so that the arrival is targeted more accurately. 31 Figure 1.13. Summary of automatic detections for the 5 HRSN stations (colored triangles). Events producing detections of FZHW at the various stations are marked by corresponding colors, while the other events are denoted by grey. (a) Results in a map view. (b) Results projected on a vertical plane. The stars denote events producing false detections shown in Figure 1.14. 32 Figure 1.14. Examples of false detections at each of the stations generated by the events marked with stars in Figure 1.13. The recordings are vertical component velocity seismograms. The red and blue lines indicate candidate FZHW and P-wave picks, respectively. The direct P-wave and FZHW picking algorithm analyzes the early portions of vertical component seismic waveforms. It is easy to analyze all 3 components, but this was not found to improve the results in the test cases examined here. The algorithm makes use of basic characteristics of the two phases (relative motion polarities, expected time difference, relative sharpness and amplitudes) to distinguish between them in a rigorous series of tests. In doing so, the first arrival is assumed to be a P-wave unless evidence for a FZHW is provided. This gives more reliability to the FZHW picks, which are not assumed to be present in all recordings as direct P-waves are. Testing the algorithm on synthetic waveforms indicates that it correctly identifies whether FZHW are present, and accurately picks the arrival times of the FZHW and direct P waves (Figures 1.8, 1.9). Systematic analysis of waveforms recorded by the 5 employed HRSN Parkfield stations yields results that are consistent with previous studies on the velocity structure and FZHW in the Parkfield region (e.g., Ben-Zion and Malin, 1991; Thurber et al., 2006; Zhao et al., 2010; Bennington et al., 2013). The automatic algorithm identifies FZHW on the overall slow NE side of the fault generated by earthquakes along most of the study area (Figure 1.13). Important 33 exceptions are the section associated with the hypocenter and largest slip of the 2004 M6 Parkfield earthquake and a small additional area of complexity. Over 50 times as many events are identified as producing FZHW at station MMNB close to the fault on the slow side compared to station VCAB on the fast side. The percentage of all events that are identified as having FZHW on the SW side is less than 2% for FROB and VCAB, which is reasonable for automatic processing. The percentage of near-fault events producing FZHW at the stations on the NE side is likely to be higher than the 10.8% observed at MMNB, because the assumed velocity values were designed to provide a reliable set of parameters for automatic picking with relatively small number of false detections. The sensitivity of the algorithm can be increased by increasing the assumed velocity contrast and reducing the minimum allowed time separation between FZHW and direct P waves, leading to a greater number of false and true picks. To test further the performance of the FZHW picking algorithm, we examine data recorded near the Hayward fault. Specifically, we run the algorithm on waveforms generated by a set of 116 events containing FZHW picks based on manual picks made by Allam et al. (2014). For stations BKS and CSP on the nominally slow side of the fault, we find that FZHW are detected in 19.5% and 11.9% of the recordings, respectively. For stations CPM and RFSB that are directly across the fault on the fast side, the algorithm (falsely) detects FZHW in 2.3% and 3.4% of the same events, respectively. In this analysis we use during the pre-processing stage a 1 Hz high-pass filter instead of 0.5 Hz, since the data quality at the used stations is not as good as in the HRSN Parkfield stations. From this we conclude that the algorithm can be applied successfully to other data sets with minimal tuning necessary. The automatic identification and picking of FZHW should be followed by analysts confirmations and possible adjustments. In addition to the basic differences between FZHW and direct P waves assumed by the algorithm, there should be indicative patterns associated with analysis of many events. These include systematic relations between minimum propagation distance along the fault and the fault-station distance for head wave detections (e.g., Ben-Zion and Malin, 1991 McGuire and Ben-Zion, 2005) and systematic moveout between the FZHW and direct P arrivals (e.g. Lewis et al., 2007; Zhao et al., 2010). In addition, since the FZHW are radiated from the fault they should have horizontal particle motion with significant fault-normal component, while the particle motion of the direct P waves should point to the epicenter direction 34 (Bulut et al., 2012; Allam et al., 2014). Similarly, an analyst should confirm and improve if needed the automatic results associated with the picking of the direct S waves. Earthquake locations are generally calculated by comparing observed travel times with predicted ones using a given velocity model. The more S picks that are available, the better the location accuracy, and this is especially true for events that are near or outside the network. A high percent of S picks is also essential for tomographic inversions for V P /V S ratios. Misidentification of FZHW as direct P waves may produce a bias in arrival times and errors in first motion polarities. If several stations have first arriving FZHW that are erroneously picked as direct P waves, there will be a location bias toward the fast side of the fault (e.g., Ben-Zion and Malin, 1991). Similarly, using stations with FZHW first arrivals to obtain first-motion polarities will produce errors in focal mechanisms (e.g. McNally and McEvilly, 1977). All the proposed picking algorithms were designed and calibrated specifically for local earthquakes recorded in fault zone environments, but they may have utility in other circumstances. We found that the S picking algorithm can work for regional earthquakes in some cases, but more work is needed to ensure that it performs consistently on these and other types of data. It is likely that some of the parameters used in this work will need to be changed to account for the different waveform character in larger scale source-receiver settings. The head wave picker was designed with the intention of identifying FZHW, and has not been tested on more classical head waves generated by horizontal interfaces such as Pn phase. A significant difference, in addition to first motion polarity, is that Pn and other such phases are expected in all events recorded beyond a certain distance. On the other hand, the basic operating principle of the FZHW picker is that no head wave is present unless a series of rigorous tests suggest otherwise, so it is not directly translatable to Pn type phases. Adaptations of the algorithms for identification of regional phases will be the subject of future work. The accuracy of the S-wave picking algorithm is strongly tied to how well the polarization filter can remove P-wave energy. The filter used here was constructed using only the polarization quantities rectilinearity and vertical incidence angle. Additional polarization quantities, such as ratio of transverse energy (e.g. Cichowicz, 1993) may help to more effectively separate the phases. In particular, this would be effective for earthquakes with vertical incidence angles that differ from the expected values for P and S. Including information on the expected radiation patterns of P and S waves and noise level (e.g. Rastin et al., 2013) can also improve the performance of the automatic 35 pickers. Typical event detection scenarios (e.g. Joswig, 1995; Gentili and Michelini, 2006; Ross and Ben-Zion, 2014) use "association" to combine detection information from different stations and/or components before making a decision. Similar types of analysis could be used for the output of the FZHW picker as well as the direct P and S pickers. For example, this could involve requiring a minimum number of near-by stations to flag the same event before making a given pick. Ultimately, the number and types of additional analyses, as well as the set of algorithm parameters used, will depend on a balance between detection rate, reliability of having correct phase identification, and required accuracy of the picked arrival times 36 2. An improved algorithm for real-time S-wave picking with application to the (augmented) ANZA network in southern California (Ross et al., 2016) 2.1 Introduction Phase picking algorithms are becoming increasingly valuable in seismology for automated processing of large seismic data sets. Historically, phase arrivals at seismic networks were picked manually by trained analysts, but the present rate of data accumulation worldwide makes this now unfeasible in most circumstances. Automatic picks are critical for many routine calculations made by regional seismic networks including earthquake detection, hypocenter location, magnitude calculation, and derivation of source mechanisms. High quality phase picks are necessary for numerous other purposes outside the normal scope of seismic network operations, such as in seismic tomography and analysis of shear wave anisotropy. A variety of reliable methods for automatic P-wave picking have been available for decades (e.g., Allen, 1982; Sleeman and van Eck, 1999; Saragiotis et al., 2002), but S-wave picking algorithms have generally lagged behind in reliability. The majority of these methods operate in the time domain, using characteristic functions such as STA/LTA (Allen, 1982) giving the ratio of a short-term moving average to a long-term moving average. A number of recent studies provided additional developments for S-wave picking algorithms. Baillard et al. (2014) used moving window kurtosis functions to identify transient signal onsets combined with polarization analysis to identify and distinguish P-waves from S-waves. Gentili and Michelini (2006) and Wang and Teng (1997) used machine learning algorithms with time series detectors to identify S-wave arrivals. Hildyard et al. (2008) and Nippress et al. (2010) employed predominant period and/or kurtosis characteristic functions to pick S-waves. Rawles and Thurber (2015) used a nearest- neighbor and similarity approach for automatic picking of S arrivals. Kurzon et al. (2014) employed an iterative SVD algorithm to construct polarization filters for separating P and S energy in seismograms, and an SVD-derived function to make S picks from the filtered seismograms in real-time. Ross and Ben-Zion (2014, referred to hereinafter as RBZ14) utilized a sliding covariance matrix to construct polarization filters for separating P and S energy, and STA/LTA combined with kurtosis detectors in tandem to lock in on a well-defined S-arrival. The method, denoted as the DBSHEAR algorithm, was shown to perform well on various data sets recorded by regional 37 and local seismic networks (RBZ14; Wu et al., 2015; Qiu et al., 2015; Share et al., 2015) and is developed further in this work. Many picking algorithms have been designed to work on event-style traces, which are short time segments cut around earthquake signals from the continuous waveform data. This requires that the earthquake was detected previously, such as by a regional network. These methods cannot be used to detect new earthquakes themselves and are further unsuitable for use in real-time. A real-time picking and detection algorithm must be capable of scanning through continuous waveforms and identifying whether a given wavelet likely contains an earthquake with no prior knowledge. Since the majority of recorded seismograms consist of background noise or transient signals not related to earthquakes, a real-time picking algorithm must determine when a seismic phase has likely arrived at the instrument. In this paper we describe an improved DBSHEAR algorithm for automatic S-wave arrival picking in data recorded by local seismic networks. The method is capable of operating in streaming real-time on packets of continuous waveform data without prior knowledge about whether an earthquake has occurred. The algorithm is tested on data from 123 stations around the San Jacinto fault zone (SJFZ) that are archived by the ANZA network in southern California. The instruments (Figure 2.1) are a mixture of broadband and short-period seismometers, and accelerometers. The network includes the core AZ network stations, a temporary PASSCAL deployment (YN) around the SJFZ, several Plate Boundary Observatory (PB) borehole stations, adjacent Southern California Seismic Network (CI) stations, and the BVDA and GVDA (SB) borehole arrays. This collective set of stations is denoted below as the ANZA+ network for simplicity. The employed data were recorded during two separate time periods. The first data set consists of all the continuous waveform data recorded from May 1, 2014 to June 6, 2014, chosen because it coincided with a 1108 instrument nodal array that was temporarily deployed in the SJFZ (Ben- Zion et al. 2015). More than 31,000 S-wave picks were handmade by analysts during this period, yielding an extensive high quality pick and event database, and we used this set to systematically test the quality of our method’s picks. The second data set consists of the entire year 2013 of continuous waveform data recorded by the ANZA+ network. This data set is used to systematically study how the method improves upon the existing/previous real-time detection capabilities and results of the ANZA+ network. 38 Figure 2.1. Map of the San Jacinto fault zone in southern California. A total of 123 stations (triangles) in the region were used by this study, which consist of AZ, YN, PBO, SB, and CI networks. The star indicates the hypocentral location of the example event used in Figures 2.2 and 2.3. Black lines denote the surface expressions of faults in the region. 2.2 Methodology Our technique for automatic S-wave picking is based on the DBSHEAR algorithm of RBZ14. The method has three stages, which can be summarized as follows: 1) use polarization analysis on a 3-component time series to remove as much of the P-wave energy from the seismogram as possible, 2) use STA/LTA detectors on the resulting polarized horizontal traces to make a trial S- pick in the vicinity of the S-wave, and 3) use kurtosis detectors in conjunction with the trial pick to lock on the S-arrival. One problem with the method is that if the polarization analysis fails to remove enough P-wave energy, the trial S pick (in step 2) can mistakenly lock on the P-wave 39 arrival. In this study, we introduce an alternative method for making the trial S-pick more reliable, leading to significantly more S-picks that are also more accurate. Figure 2.2. A step by step demonstration of the picking method. A) vertical component velocity seismogram, b) N-S component velocity seismogram, c) E-W component velocity seismogram, d) S-polarization filter We describe the three stages of the algorithm in detail so that the entire method, including the implemented modifications, is explicitly clear. The first stage of the algorithm is designed to identify and remove as much of the P-wave energy from the horizontal seismograms as possible. This part of the procedure is unchanged from RBZ14 and uses a polarization filter constructed from particle motion analysis. Figure 2.2 contains a step-by-step visual demonstration of the algorithm, which is referred to throughout the remainder of the methodology section. In Figures 2.2a-c, there is an example 3-component trace from a local earthquake (red star, Fig. 2.1) recorded at 40 Figure 2.2, continued. e) S-polarized N component velocity seismogram f) smoothed STA/LTA calculated on e), g) kurtosis function calculated on e), h) kurtosis rate function. station JFS1 (red triangle, Fig. 2.1). To obtain the polarization filter, a covariance matrix is first calculated from a finite sample of 3-component data as ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ = ) , Cov ) , Cov ) , Cov ) , Cov ) , Cov ) , Cov ) , Cov ) , Cov ) , Cov Z (Z E (Z N (Z Z (E E (E N (E Z (N E (N N (N σ , (1) where Z, N and E are the vertical, north-south and east-west components, respectively. The covariance between any two components X and Y is given by ∑ = = M i i i y x M Y (X 1 1 ) , Cov , (2) where M is the number of samples. RBZ14 used a 3 s sliding window for instruments recording with a sampling rate of 40-200 Hz. This parameter may need to be adjusted somewhat for significantly different data sets. From the covariance matrix for a particular window, the eigenvalues ( ) 3 2 1 λ λ λ ≥ ≥ and eigenvector matrix ( ) 3 2 1 , , u u u u = can be used to calculate 41 various aspects of the particle motion (e.g. Jurkevics, 1988; Vidale, 1986). Specifically we use the rectilinearity, ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + − = 1 3 2 2 1 λ λ λ r , (3) and the apparent vertical incidence angle, ) ( cos 13 1 u − = ϕ . (4) These two quantities are then used to construct an S-wave polarization filter, )) cos( 1 ( ϕ − = r s . (5) The polarization filter has a range of [0,1] and is large when the particle motion is primarily in a horizontal plane, and small otherwise. Applying this process to the 3-component data in Figure 2.2 results in the polarization filter shown in Fig. 2.2d. Since this filter is calculated every sample as the window slides, it is actually a time series with the same length as the data traces. The polarization filter is used to modulate the horizontal traces by simple multiplication at each time step, attempting to remove as much of the P-wave motion as possible (Fig. 2.2e). Indeed, when comparing Fig. 2.2e with Figs. 2.2a-c, the P-wave energy is considerably weaker. For these reasons, RBZ14 called these S-polarized traces. For more discussion of this stage, see RBZ14. With the polarization filter applied to the data, it generally removes enough of the P-wave that its amplitude is small relative to the S-wave. RBZ14 calculated the ratio of a short-term moving average to a long-term moving average (STA/LTA; Allen, 1978, 1982) on the S-polarized horizontal traces (e.g. Fig. 2.2e). The STA/LTA effectively measures signal-to-noise ratio, or in some sense, the total energy in the window as a function of time. RBZ14 used an STA window length of 1 s and an LTA length of 10 s. Mulder et al. (2015, submitted) found that adding at least 1 count white noise to the polarized traces helps to stabilize the STA/LTA function, which minimized spurious detections. An STA/LTA function with these parameters is shown for the S- polarized trace in Fig. 2.2f. After calculating the STA/LTA function, they applied a trigger mechanism that turned on when an STA/LTA threshold of 5 was met, and turned off when the STA/LTA function fell to a threshold of 1. The trigger duration was required to be more than 1 s. This ensured the trace contained a transient signal with sufficiently large amplitude. Then, they made a trial S-pick when the STA/LTA function was maximized, which generally occurred in the vicinity of the S-wave if the polarization analysis was successful. This is also the case in Figure 2.2f. If the polarization analysis is unsuccessful at removing the P-wave, or the P-wave arrival is 42 simply too strong relative to the S-wave arrival, the STA/LTA function’s maximum potentially occurs during the P-wave instead, leading to a mispick. Figure 2.3. The modified procedure for selecting the trial S-pick to aid in cases where the P-wave has the largest amplitude. A threshold is set at 60% of the peak amplitude of the STA/LTA function on the polarized horizontal trace (dashed line). All local maxima above the dashed line are identified, and the one that is the latest in time is used as the trial pick (diamond). The trial pick is used to define the center of a 1 s duration window for the kurtosis rate function. The pick is ultimately refined to the largest value of the kurtosis rate within the 1 s window. This allows the S-wave arrival to be correctly picked when the polarization analysis does not work and the peak of the STA/LTA function is not near the S-wave. We modify this stage of the DBSHEAR algorithm to account for this shortcoming. In our tests, we found that there are generally two large broad peaks in the STA/LTA function (Fig. 2.2f). In some cases the first peak is larger, while in others the second is. This is often an indication of how well the polarization filter works, so we develop a procedure for identifying each of the largest peaks in the window to use the most appropriate peak for the trial S-wave pick. This is demonstrated in Figure 2.3 for a different station in the network where the polarization analysis did not work well. First we smooth the STA/LTA function using a simple moving average with a 43 window length equal to the STA length (1 s). This removes all of the high frequency oscillations from the STA/LTA function and helps identify only the major peaks in the time series. Next, we define a threshold equal to 60% of the maximum STA/LTA amplitude (dashed line, Fig. 2.3b). Then we identify all local maxima with a value above the threshold; there are commonly 1-2 maxima above the threshold and rarely 3 or more. Instead of making the trial pick on the single largest maximum like in RBZ14, we make the trial pick on the peak that is latest in time. For the example in Figure 2.3, this maximum (and trial pick) is indicated by a red diamond. The motivation for this is that the latest peak is most likely to be closer to the S-wave than the P-wave because STA/LTA measures moving energy. Setting a threshold at 60% of the maximum imposes the requirement that the only maxima considered had relatively high energy. Two high-energy maxima in the STA/LTA function are generally the result of the two body wave phases and a drop in energy between the arrivals. As a result, a phase with energy comparable to a P- or S-wave typically needs to be present to create an additional peak above the 60% level. The most important aspect here is that this procedure is very effective at bypassing the P-wave even in the rare cases with more than two maxima. The details on overall effectiveness of the method are discussed extensively in the next section, but first the remainder of the method is described. With the trial pick obtained from stage 2, the last stage of the algorithm is now possible. As the trial pick tends to be located somewhere around the S-wave, the pick must be refined to the location of the actual phase arrival. For this stage, we follow the exact procedure of RBZ14, who used a moving kurtosis function (Saragiotis et al., 2002) to lock in on the arrival. Kurtosis is the 4 th order statistical moment and measures peakedness; it is quite effective at identifying abrupt changes in the amplitude of a time series because it is extremely sensitive to outliers. Mathematically it is defined as, 3 ) 1 ( ) ( 4 1 4 − − − = ∑ = x M i i s M x x K , (6) where M is the number of samples, } ,..., , { 2 1 M x x x x= is the finite length data sequence, x is the sample mean, and s x is the standard deviation. RBZ14 used equation (6) with a sliding window length of 5 s on the S-polarized horizontal traces, which is shown for the example traces in Figure 2.2g. They found that the 5 s value was appropriate for sampling rates 50-200 Hz. From here, the kurtosis rate was calculated (Fig. 2.2h) and used to identify the phase arrival. This is accomplished by using the trial S-pick (diamond, Fig. 2.3b) to define the center of a 1 s duration window on the 44 kurtosis rate (Fig. 2.3c). The pick is then refined to where the kurtosis rate is maximized within the 1 s window, and adjusted slightly to the closest zero preceding the maximum. At this point, further P-wave discrimination could potentially be achieved by examining the apparent incidence angle, φ, in the vicinity of the final pick (e.g. Baillard et al., 2014; Kurzon et al., 2014). A value of φ that is greater than 45 degrees is indicative of non-S-wave energy and corresponding picks could be rejected, however such a test has not been implemented in this work. 2.3 Real-time implementation and picking accuracy RBZ14 demonstrated their S-picking algorithm only on waveforms of earthquakes that were detected and located previously. This was done with the intent of re-picking phase arrivals for an event at different stations of interest, relying on data traces to be cut around the seismogram beforehand. To use the method for detecting and locating events, however, it must be capable of working on continuous data volume with no information a priori. Here we demonstrate that the RBZ14 algorithm is suitable for real-time processing of continuous data, and furthermore capable of building a high quality seismicity catalog from scratch when combined with the Antelope Environmental Monitoring Software package (BRTT Inc.). For these purposes, we use all the continuous waveform data recorded by the ANZA+ network (Fig. 2.1) from May 1, 2014 to June 6, 2014. To aid in the task of detection and location we use the Antelope platform, which is presently employed by a number of permanent seismic networks around the world, including the ANZA+ network. We first use Antelope’s DBDETECT software to run on all the raw continuous data and make P-wave arrival detections. The parameters used for this are the default settings provided by the manufacturer, and are the same as those used internally by the AZ network to detect earthquakes. For each P-wave detection made on the continuous data, the 3-component data is cut starting 10 s before and 14 s after the P-wave detection. Then, the resulting trace is passed to the DBSHEAR algorithm including our modifications to stage 2. The DBSHEAR algorithm attempts to make S-wave detections on each of the horizontal components, and if both components are successful, the one with the larger signal-to-noise ratio is chosen. We note that a phase detection, as described here, is not the same as a phase arrival pick. This is because we are detecting new events from scratch on the continuous data, and only a subset of the phase detections will be later converted into actual phase arrival picks. The distinction between P-wave and S-wave detections 45 is therefore a naming convention for whether the phase detection is made by DBDETECT or DBSHEAR, respectively. Using this process on all the continuous data for the analyzed period results in 1,851,108 P-detections and 422,900 S-detections. It is unknown at this time whether any of these phase detections are actually associated with a seismic event. For example, a phase detection could turn out to be an instrument glitch or arbitrary transient signal not related to an earthquake. The steps described are capable of working on a completed seismic database or on packets of data streaming in real-time. In both situations, a piece of continuous data is passed to the phase detection algorithms to see if it contains something that looks like an earthquake, however most of the time it will just contain noise. We have successfully implemented the DBSHEAR algorithm in real-time on the AZ network, and configured DBSHEAR to run only when a phase detection is made by DBDETECT on a given station somewhere. This reduces the overall amount of computation considerably, since the polarization filters and kurtosis functions only need to be calculated for a small percentage of the total data. The DBSHEAR method works exactly the same regardless of whether it's working on archived data, or packets of data streaming in in real time. After a sufficient number of phase detections were made, they were processed by Antelope’s grid associator, DBGRASSOC, to check whether any are associated with a coherent earthquake- type source. The grid association is performed by back-projecting all the phase detections to different nodes in a grid around the SJFZ by subtracting their predicted travel times, calculated using a 1D velocity model, from the observed detection time and formally making an event detection when a minimum number of back-projected detections cluster within a 1.5 s window. The clustered phase detections are then promoted to phase arrival picks at this point in the process. Whenever an event is detected, the phase arrivals associated with the coherent source location are inverted using the GENLOC algorithm of Pavlis et al. (2004) to get a formal hypocenter and origin time. The described process resulted in 1,034 earthquakes detected and located, along with 15,055 P-picks and 10,163 S-picks. We then hand-picked more than 31,000 S-wave arrivals for every single earthquake detected by the ANZA+ network during this period. For every automatic pick in common with this manual set of picks, which totaled 8,246 S picks, we calculated the error relative to the analyst pick. Figure 2.4 shows a histogram of the picking errors for this data set. The mean of this distribution is very close to zero and the picks are not overly biased towards one side or the 46 other. Various additional statistics are listed in Figure 2.4 as well. Figure 2.5 shows the empirical cumulative distribution function (ECDF) for the values in Figure 2.4. From this ECDF, it is clear that 64% of the 8246 picks are within 0.10 s of the analyst pick. Further, 90% of all picks are within 0.25 s of the analyst pick, and 97% are within 0.5 s. Analyst picks at the ANZA+ network typically have an error of around 0.1 s, and almost two-thirds of the picks made by our algorithm are as accurate. Figure 2.4. Histogram of picking errors compared to analyst picks for 2014 dataset. A total of 8,246 manual S-picks were used in the comparison, and the average automatic pick had an accuracy of -0.02 +/- 0.17 s. The picking errors can further be studied as a function of signal-to-noise ratio (SNR) and epicentral distance. Figure 2.6 shows that the median pick error, indicated by the blue line, is about 0.1 s or less up to about 60 km. Beyond this distance, the median pick error is slightly higher but still remains less than 0.2 s up to about 100 km where the data becomes scarce. Figure 2.7 presents the same picks as a function of SNR and the median pick error is seen to be essentially flat with 47 regards to SNR. Taken together with the results in Figure 2.6, this indicates that the method is more affected by the changes produced in the waveforms due to longer propagation rather than magnitude itself. However these changes are not very pronounced, and the method works well over the full range of these different examined scales. The SJFZ has a wide diversity of focal mechanism types (e.g. Bailey et al., 2010) leading to greater variability in the character of observed waveforms. We therefore have no reason to expect the quality of the picks to worsen on other data sets, and indeed our preliminary testing elsewhere supports this expectation. Figure 2.5. Empirical cumulative distribution function (ECDF) describing the same results in Figure 2.4. The 8,246 picks being compared here are within 0.1 s 64% of the time, 0.25 s 90% of the time, and 0.5 s 97% of the time. 2.4 Systematic application to the ANZA+ network With the accuracy of the method established in the previous section, we now proceed to discuss the application of the algorithm to a data set over a longer time period. The data set used for this part consists of the entire year 2013 of continuous waveform data recorded by the ANZA+ 48 network. We applied the same detection and location procedure from the previous section to this larger data set to detect and locate every single earthquake from scratch. In this manner, side-by- side analysis was possible to study the effect that introducing DBSHEAR to the workflow had on the seismicity catalog. Two separate seismicity catalogs were produced for the entire year. The first catalog was produced using just the DBDETECT program alone. The resulting catalog is identical to the catalog produced by the AZ network’s real-time system for 2013, void of any analyst interaction. It is therefore a fully automatic seismicity catalog. The second seismicity catalog used both DBDETECT and DBSHEAR together, and as a result, the only differences are due to the DBSHEAR algorithm. In this manner, we can directly study how our method alters the locations and picks in a seismicity catalog. Figure 2.6. Picking method accuracy as a function of epicentral distance. The median absolute pick error relative to the analyst pick is given by the line. From 0-60 km epicentral distance, the median pick error is about 0.05 s, and then slightly increases to about 0.15 s beyond this range. Overall the errors are quite stable over the full distance range tested. The results of this process for the two catalogs are shown in Figure 2.8. In both catalogs, the number of events detected is 11,197. For the first catalog using DBDETECT alone, 151,839 P- 49 picks and 14,526 S-picks are made, resulting in a P to S pick ratio of 10.45. In the second catalog the number of P-picks decreases slightly to 151,771, due to various phases now being re-associated as S-waves as the result of having more information. The number of S-picks, however, increases to 110,430 (660%). The resulting P to S ratio drops from 10.45 to 1.37. Another way of looking at this is that for each P-pick that is successful on the continuous data, an S-pick is successful 73% of the time. Figure 2.7. Picking method accuracy as a function of signal to noise ratio (SNR). The median pick accuracy is nearly constant at 0.07 s over the entire range of SNR values observed in the dataset. With such a considerable increase in the number of S-picks between the two catalogs, the most important question is how the locations and origin times changed as a result. Figure 2.9a contains information on how various location resolution quantities changed between the two catalogs. The four quantities shown are all from the GENLOC algorithm (Pavlis et al., 2004), and are the major and minor axis of the covariance matrix, depth, and origin time resolution. For each parameter, negative values indicate that the uncertainty decreased, and that the degree of fit increased. For all 50 four of the parameters listed in Figure 2.9a, about 99% of the events have the same or better location accuracy. Conversely, addition of S-wave picks from DBSHEAR resulted in degradation of location accuracy for ~1% or less of the events. Since the location fit increased for the vast majority of cases, another question of interest is how these resolution changes map into 4D hypocentral shifts. Figure 2.9b shows the difference in latitude, longitude, depth, and origin time between the two catalogs. These histograms indicate that the average epicentral shift between the two catalogs is on the order of 0.5 km, and the average depth change is 1.1 km. This is measured over the 11,197 earthquakes tested. Further, the average origin time shift is 0.2 s. The distributions of these quantities, however, have very long tails. A number of events, for example, have depth changes as large as 3-4 km and origin time shifts of 0.75 s. Since the only difference between the two data sets is the inclusion of the DBSHEAR algorithm, it is entirely responsible for the improvement to the locations. Figure 2.8. Comparison of number of phase picks made when using the stock Antelope DBDETECT program alone, as compared with DBDETECT and the DBSHEAR algorithm together. Nearly the same amount of P-wave arrivals were obtained with both methods. Roughly 12,000 S-picks were made by DBDETECT alone, however when adding in DBSHEAR, this number increased to 110,400. The resulting P to S ratio has decreased from 10.45 to 1.37 (660%). Lastly we compare the hypocentral locations between the two catalogs visually around the SJFZ. Figure 2.10 shows a map of the SJFZ with events from the catalog produced with DBDETECT alone, alongside events from the catalog produced with DBDETECT and 51 DBSHEAR. A projection of the seismicity along the strike of the SJFZ is included as well. Compared to locations produced by a relocation algorithm such as HypoDD (Waldhauser and Ellsworth, 2000), where entire populations of events tend to shift systematically in particular directions, no clear patterns are observed among groups of events. Rather, the events appear to shift in directions that are uncorrelated in space. Figure 2.9. Comparison of location uncertainty statistics between the two 2013 data sets, which contained 11,197 events. (a) histograms of the change in hypocenter resolution as estimated by the GENLOC algorithm (Pavlis et al., 2004). For these histograms, a negative number indicates that the uncertainty in the fitted location decreased by including the DBSHEAR algorithm, compared to DBDETECT alone. In a very small number of events the uncertainty increased, reflecting S- picks that were of poor quality. (b) histograms indicating how the actual hypocenter quantities changed before and after including DBSHEAR. The average depth changed by 1.1 km, the average epicenter changed by roughly 500 m, and the average origin time changed by 0.2 s. 2.5 Discussion The method described herein is capable of making high-quality S-wave picks on continuous data with no knowledge of whether earthquakes have occurred. It is thus a real-time picking algorithm and is appropriate for use by an operational seismic network as well as for processing existing archives of continuous waveform data to detect earthquakes. The method comes with 52 several caveats, which should be mentioned for completeness. It was designed for picking phases on local earthquakes and requires 3-component instruments. The method has been successful at picking S-waves up to about 120 km. Beyond this distance, the development of regional phases complicates the signal and makes automatic detection more challenging. It is possible that re- parameterization of the method may produce useful results at regional scales, however no testing at regional scales has been performed. Figure 2.10. Location comparisons for the two 2013 data sets. For both data sets, the number of events is 11,197. The black set of dots were detected using the DBDETECT program only, without DBSHEAR, while the other set of dots used DBDETECT and DBSHEAR to detect 53 events. In general, no systematic shifting of groups of events is observed. The events that included DBSHEAR have less scatter overall In this study we propose a modification to the second stage of the algorithm to overcome the main flaw of the original method described in RBZ14. The original algorithm is dependent on the polarization analysis removing enough of the P-wave energy that the STA/LTA function peaks during the S-wave. In many cases, however, this does not happen and results in an erroneous pick or no pick at all. Since the method requires 3-component motion at each time step, this could also result from one of the components being out at the time of interest. As a result, the original version of the DBSHEAR method generally requires the polarization filters to work in order to make a successful pick. With our recommended modifications, this is no longer the case. To show this, we reprocess all of the continuous data during May 2014 in the manner of Section 4, except we do not use any polarization filters. All of the STA/LTA and kurtosis detectors are therefore calculated on raw traces rather than polarized. After comparing the resulting picks to the database of analyst picks as done in Section 4, the total number of S picks decreases by only 79 after processing with the grid associator. The mean absolute error of the picks increases by only 4.5%, with a 3% increase in the standard deviation of the picks. We can then conclude that the modified version of the algorithm is no longer dependent on the polarization filters working well to make a pick of good quality. We still recommend using the polarization filters because they improve the results slightly, but this is no longer a strict requirement. The original DBSHEAR algorithm typically led to a P to S pick ratio between 2 and 3. Wu et al. (2015) used the algorithm on the TAIGER arrays deployed around Taiwan in 2009 and obtained a P to S pick ratio of 2.56 over more than 8000 earthquakes. In this study we showed that the proposed modifications decreased this ratio to 1.37 for the 11,197 earthquakes detected for the year 2013. Further, the average pick’s accuracy is dramatically improved, and 90% of the picks tested are within 0.25 s of the corresponding analyst pick. In this study we make the formal distinction between a real-time algorithm and a “re-picking” algorithm. Many of the picking algorithms discussed in the literature (e.g. Gentili and Michelini, 2006; Wang and Teng, 1997; Baillard et al., (2014)) are re-picking methods and are unsuitable for real-time use in their present form. A re-picking algorithm expects an event to have been previously detected through some alternative means. The data are then cut around the earthquake, and the method is applied to pick the phases. Consequently, these types of methods are generally 54 not capable of being applied to continuous data, where the goal is to detect earthquakes for the first time. The continuous data frequently contains noise bursts and other transient signals, such as from vehicles, animals, or cultural sources. A real-time picking algorithm must be capable of distinguishing between all of these types of signals and deciding when to make a pick. If false picks are made too commonly, the association process can produce spurious event detections. Real-time picking algorithms are used by permanent networks to locate earthquakes as part of the detection process, and are responsible for making the first picks available to a network. Further development on this subject can therefore have a significant impact on the locations provided in seismicity catalogs. In this study it was shown that by including our S-picking algorithm in the ANZA network’s workflow, the average hypocenter was improved by more than 1 km. These improvements will propagate through to downstream applications, such as tomography studies. An important aspect of picking phase arrivals before any events are detected or located is that it allows for a fully self-consistent database. If a re-picking algorithm is applied to events that were detected previously by a permanent network, the picks will be different from what the network produced internally. All of the information derived from the network’s original picks -- which includes locations, magnitudes, and source properties -- will be inconsistent with the new picks produced. The fact that the results are not self-consistent across the entire database introduces a form of uncertainty that is generally not discussed. By further developing algorithms that can operate on the continuous data directly, this uncertainty can be eliminated while maintaining self- consistency throughout the database. 2.6 Conclusions We developed a modification to the DBSHEAR algorithm introduced originally by Ross and Ben-Zion (2014) for picking S-wave arrivals on seismic waveform data. The method is a three- stage procedure designed to account for a wide variety of expected waveform characteristics. It first uses polarization analysis to remove as much P-wave energy from a seismogram as possible. Next, STA/LTA detectors are run on S-polarized traces to identify the neighborhood of the S- wave. Lastly, a moving kurtosis function is used to refine the final pick to a well-defined region of the seismogram. In this study we demonstrated for the first time that the method works on packets of data streaming in real-time, with no-prior knowledge of whether an event has occurred. We then applied the method to more than a year of continuous waveform data recorded at over 55 100 stations around the San Jacinto fault zone, and built a seismicity catalog from scratch. The resulting phase picks associated with events in this catalog were compared with more than 8,000 handmade phase picks to rigorously examine the accuracy of the method. In doing so, 90% of all S picks were found to be within 0.25 s of the corresponding analyst pick. We demonstrated that by using the DBSHEAR method together with the Antelope platform, the average event hypocenter improved by more than 1.3 km, and that the number of S-picks made increased by more than 600%. 2.7 Data and Resources Seismic data used in this study are gathered and managed by the following networks: the Anza network (AZ) and the SJFZ (YN) network are operated by IGPP, University of California, San Diego <http://eqinfo.ucsd.edu/> (last accessed August 2015); the SCSN network (CI) is the Southern California Seismic Network operated by Caltech and USGS <http://www.scsn.org/> (last accessed August 2015); the PBO network (PB) is the Plate Boundary Observatory operated by UNAVCO <http://www.earthscope.org/science/observatories/pbo> (last accessed August 2015); the SB network is operated by University of California, Santa Barbara <http://nees.ucsb.edu/publications> (last accessed August 2015); Gaps in data were filled by the IRIS Data Management Center (DMC) <http://www.iris.edu/dms/dmc/> (last accessed August 2015); instrumentation for seismic stations was provided by IRIS PASSCAL <http://www.passcal.nmt.edu/> (last accessed August 2015). 56 3. An earthquake detection algorithm with pseudo probabilities of multiple indicators (Ross and Ben-Zion, 2014b) 3.1 Introduction Automated algorithms for earthquake detection require the use of indicator variables that tend to correlate in amplitude with the arrival of earthquake waveforms. These indicators can operate in either the time or frequency domain. One common example is the ratio (e.g. Allen, 1978) of a short term moving average (STA) to a long term moving average (LTA). Polarization analysis using singular value decompositions (e.g. Rosenberger, 2010) or covariance matrices (e.g., Jurkevics, 1988) can provide additional indicator variables for detections. Kurzon et al. (2014) have used polarization analysis to separate real time data into channels associated with P and S waves, and analyze various indicator variables to achieve significantly improved P and S wave picks. Withers et al. (1998) compared many different indicator variables and find that the STA/LTA is often the most reliable. At present seismometers commonly record in a continuous mode, providing new opportunities for detecting many small events. An ideal detection framework will identify all events recorded by a network down to the noise level. However, the current methods fall considerably short of this goal, and specialized techniques such as using templates (e.g. Yang et al., 2009; Peng and Zhao, 2009) often detect many more small events than the standard automated methods. To improve the performance of automated algorithms, multiple STA/LTA detectors can be used with different STA and LTA windows and/or different frequency bands. Each STA/LTA indicator can be evaluated independently with different trigger thresholds. If these thresholds are met simultaneously by a sufficient number of indicators across enough stations (in the association phase of the algorithm), a detection can be made. The primary issue associated with this approach is that any one indicator may often trigger erroneously when an earthquake is not present. The cutoff threshold for triggering is typically adjusted so that a suitable number of false detections are allowed (e.g. Nippress et al., 2010), and the smallest detected earthquakes are therefore limited by the accuracy of a single indicator. As such, the collective amount of information available from running multiple STA/LTA detectors is not used in a way that takes full advantage of the redundancy. In this study we describe a method that circumvents this problem by combining the information from each indicator together before any thresholds are applied to the data. This is done in a non- 57 parametric way that does not involve weighting schemes and offers great flexibility in choosing which indicator variables to use. The primary advantage is that the joint detection threshold used on the combined set of information can be lowered significantly compared to any single indicator. This results in more detections at a more reliable rate. We demonstrate the potential of the algorithm in a small region around the San Jacinto fault zone, and compare the results to those obtained by the Southern California Seismic Network (SCSN). 3.2 Methodology Our methodology for automated detection consists of several simple steps. We first choose the type and number of indicators to be used in the detection process. We focus on STA/LTA indicators with 10 different short- and long-term windows and different frequency bands (Table 3.1). The specific employed indicators are chosen to detect small earthquakes by a local network. In this section we describe the approach used after testing the methodology; the details of the tests are documented in the subsequent section. Continuous waveforms are separated into one hour segments that overlap slightly to ensure that no events are missed at the edges. We find that one hour of data is long enough to characterize the wavefield appropriately relative to the time scales of M < 4 earthquakes, and is computationally efficient. Also, one hour is short enough to retain the unique temporal characteristics of the noise background for that interval. We assume that the majority of the signals in one hour of data are associated with ambient noise rather than earthquakes. Fig. 3.1a shows an example of a 1 hour vertical component trace on 03-27-13 containing 21 earthquakes. Nearly all of these earthquakes have M < 2.0 and only several are visible with the scale of the figure. The ten different STA/LTA detectors are used on the vector magnitude (i.e., 2 2 2 Z E N + + ) of each hour window (Fig. 3.1b). Similar results are obtained if the individual components are used and the requirement of N stations in the association phase is replaced with 3N traces. The result of a single STA/LTA detector on the vector magnitude for this time window is shown in Fig. 3.1c, using an STA window of 3 s, LTA window of 10 s and a high pass filter at 3 Hz. Large values of the STA/LTA ratio indicate significant changes in the amplitudes of the wavefield over time, while accounting for variability in the noise level. We now consider only the values of STA/LTA during this time window while temporarily discarding the time associated 58 with each value. This leads to a sampling distribution (Fig. 3.1d) that characterizes the range and likelihood of different STA/LTA ratios that were observed during that hour. Figure 3.1. A visual illustration of the detection algorithm. a) Example 1-hour vertical seismogram from the orange station in Fig. 3.2 with 21 earthquakes. Two example events detected are shown with an amplified scale. b) Vector magnitude of the 3 components for the same 1-hour data. c) STA/LTA detector trace with an STA window of 3 s and LTA window of 10 s. d) Histogram of the obtained STA/LTA values. e) Empirical CDF of the STA/LTA values, giving the percentile of each value relative to all others in that hour. f) A pseudo-probability time series for the examined hour of data. g) A joint PPTS obtained by multiplying 10 individual PPTS at each time step. A probability threshold of 0.8 (dashed line) is generally found to distinguish a signal from noise at a given station. There are 17 events above the threshold with the missed 4 buried in the coda of the largest ones. Circles indicate the peak values for the two events magnified in a). 59 Since we employ here a total of 10 different STA/LTA detectors, we obtain 10 different sampling distributions. For each detector, we use its distribution as a reference and then compare each value of the STA/LTA time series to its relative ordering within the sample. To do this, each STA/LTA distribution is used to compute an empirical cumulative distribution function (ECDF) (e.g. D’Agostino and Stephens, 1986). The ECDF describes the cumulative fraction of values less than or equal to a given number, } { 1 ) ( 1 s x n s F i n i ≤ = ∑ = 1 , (1a) where F is the value of the ECDF at a given value of STA/LTA (denoted by s), x i is the ith value of the particular sample with n total values, and ⎩ ⎨ ⎧ > ≤ = ≤ s x s x s x i i i if 0 if 1 } { 1 . (1b) Fig. 3.1e shows the ECDF corresponding to the distribution in Fig. 3.1d. The horizontal axis indicates the STA/LTA size and the vertical axis indicates the fraction of values in the one hour window that are less than or equal to the corresponding STA/LTA amount. The lowest STA/LTA ratios on this plot have ECDF values close to zero, whereas the largest ones have ECDF values close to 1. We use the ECDF to define a mapping, through which we transform each STA/LTA value, s, in the time series (e.g., Fig. 3.1c). with the corresponding ) (s F . This converts the entire time series into a new one having a range of [0, 1], with each value exactly equal to the original value’s percentile relative to all others. We call this a pseudo-probability time series (PPTS), because at each time step we take the corresponding probability value from the ECDF. The "pseudo" distinction signifies that we use the PPTS solely to describe relative ordering within a sample. The PPTS corresponding to the example time series (Fig. 3.1b) is shown in Fig. 3.1f. It is clearly more erratic than the STA/LTA series in Fig. 3.1c and regularly reaches values greater than 0.9 when no earthquakes are present. Since by definition no more than 10% of the values are allowed to reach 0.9, the mapping becomes very useful when used with multiple indicators. In the present application we have ten PPTS for each one hour window with values in the range [0, 1], which can be combined by simple multiplication at each time step. While a single PPTS frequently reaches high values over the course of the time window (Fig. 3.1f), the likelihood of all ten PPTS having large values is extremely low. Furthermore, low to medium values tend to cancel 60 each other, leading to clear robust peaks in the product PPTS (Fig. 3.1g). Instead of a set of ten STA/LTA detectors that are used with ten ad hoc detection thresholds, we have one joint PPTS detector with a single detection threshold. We can compute this for each station just like a standard STA/LTA detector, but the advantage is that with our method, the redundancy built into the product PPTS allows the single cutoff threshold to be lowered more than any one STA/LTA detector. It may seem that the product PPTS in Fig. 3.1g is much noisier than the basic STA/LTA series in Fig. 3.1c, but there are distinct advantages to the former. First, the PPTS is bounded between [0, 1] whereas the STA/LTA is only required to be greater than or equal to zero. The bounds imposed by the PPTS, along with the fact that it is derived from multiplication, make it progressively more difficult to attain high values. As an example, when all 10 PPTS have a value of 0.99 or higher the product is ~0.90, while if all 10 have a value of 0.93 or lower the product is ≤ 0.49. If any one PPTS has a value of 0.5 (the median) at a given moment, the product will automatically be less than or equal to this. We find that when the product PPTS value reaches ~0.8, there is essentially always an earthquake or strong departure from the preceding noise background that should be recognized for that station. The existence of this feature at other stations can be used to refine the detection further. However, we find that the product PPTS even at a single station performs well in flagging out unusual portions of the data. In the following we adopt a very simple association scheme using the product PPTS detectors for testing the methodology. For each product PPTS we initiate a trigger window when a value of 0.3 in probability units is reached, and the window is turned off when this falls to 0.1. We require that the window be a minimum of 2 s in duration and reach a peak value of at least 0.82. If these criteria are not met the window is discarded. From the set of all remaining trigger windows for each station, we define detections when a minimum of six stations have overlapping trigger windows (no minimal overlap time is required). These detection criteria are found to have a false detection rate of ~2% based on the examples presented in section 4. 3.3 Indicator Variable Testing Here we discuss the tests leading to the best parameter combination for the detectors used in this work, and the choice of the number of indicators. We performed an exploration of the STA, LTA, and high-pass filter parameter space for ten detectors with a 24 hour window of time, starting 61 on 2013-01-01:00:00:00. The employed recordings are from 13 ANZA stations around the San Jacinto fault zone (green triangles in Fig. 3.2). We initially focused on finding parameters that detected all the 19 local earthquakes listed by the SCSN (available at http://www.scecdc.scec.org/) in the considered region of space and time. This was achieved using first a coarse exploration of the parameter space, followed by more detailed exploration in the neighborhoods of well performing parameters, using many combinations including those listed in Table 3.S1. The best performance was obtained by the 10 parameters of Table 3.1, leading to detection of 82 earthquakes in this 24 hour window. Of these, 81 events were visually confirmed at six or more stations, with only one false detection for that combination. For comparison, using only one of the ten detectors with STA and LTA windows of 3 s and 15 s, and the same association scheme with on and off trigger thresholds of 3.5 and 1, detected just 21 earthquakes which is close to the SCSN results. Figure 3.2. Events detected by the developed method (blue dots) during 7 days of data recorded by the 13 ANZA stations (green triangles) that form a subset of the SCSN stations (purple triangles). The algorithm detects 3.13 times the number of earthquakes listed in the SCSN catalog for the examined area and time (red dots). The red boxes mark the locations of two events shown in Fig. 3.3. 62 Using the parameters of Table 3.1, we tested the number of indicators necessary to achieve stability in the product PPTS for the first hour of the data previously examined. The full 24 hours was not used because a large number of false detections are present when using only a few indicators. We started with two indicators, and ran the algorithm on each station. We recorded the number of earthquakes detected and each was visually confirmed. It was also noted how many earthquakes were missed. The results of repeating this process while progressively increasing the number of indicators are summarized in Table 3.2. The ordering of the indicators was found to be relatively unimportant as each indicator is an STA/LTA detector. It can be seen that it takes eight indicators of the form used to achieve a stable result. We have, however, used ten consistently throughout this study because we find that the extra stability does not affect the performance of the algorithm. Table 3.1. Best performing set of 10 different STA/LTA detectors in the 24 hour test period. STA (s) LTA (s) High Pass (Hz) 3 10 3 3 15 3 3 20 3 3 25 3 3 30 3 2 5 5 2 7 5 2 9 5 2 11 5 2 13 5 3.4 Illustration from the San Jacinto Fault Zone Region With the parameters and number of indicator variables optimized, we now demonstrate the algorithm in comparison to the SCSN. The region of study is the Southern California plate boundary area depicted in Fig. 3.2. As before, we use only the 13 ANZA stations shown. We focus on a 7 day window of time beginning on 2013-01-01:00:00:00, in which the SCSN detected 91 earthquakes (red dots in Fig. 3.2). The stations used by the SCSN include the green triangles as well as all the purple ones. The SCSN uses more complex event association than we do, but their detection criteria requires essentially 5 or more high quality phase picks (E. Hauksson, personal comm., 2013), which is similar to our approach. We ran our algorithm on the vector magnitude of the three component traces for each station using the detection criteria described in section 2 and 63 the 10 STA/LTA detectors listed in Table 3.1. Each detection was visually confirmed, and P and S wave arrival times were picked by hand. Figure 3.3. Seismograms of two example detected earthquakes and corresponding product PPTS. a) Three-component traces and joint PPTS (at station CRY) for the event in the SE red box of Fig. 3.2. The red and blue dotted lines indicate, respectively, the P and S arrivals of the earthquake detected by a minimum of six stations. The black horizontal dashed line indicates the detection threshold. b) Corresponding results for the event in the NW red box of Fig. 3.2 at station KNW. The detection algorithm identified 456 possible earthquakes, of which 20 did not have at least six clear P and S picks. Of the remaining earthquakes, 419 were located successfully using the HYPOINVERSE code (Klein, 2003) with the velocity model of Hadley and Kanamori (1977), and a total of 285 were found to be in the region of Fig. 3.2 (blue dots). Three component seismograms for two of these earthquakes are shown in Fig. 3.3 along with the example PPTS at individual stations used to make the detection. The detected earthquakes are indicated by red and blue lines and are associated with a clear peak of the PPTS. In Fig. 3.3a, the algorithm identified only the S wave, while in Fig. 3.3b with higher signal to noise ratio both P and S waves are detected. The S wave picks in Fig. 3.3 were obtained using the peak PPTS value as a starting position, while the P 64 wave picks were obtained using a lower PPTS trigger threshold of 0.3. Both picks were refined using nearby peaks of the time derivative of PPTS values (Fig. 3.S1) derived from the individual trace components (vertical for P and horizontal for S). These results were found to be in good agreement with phase picks based on the Rosenberger (2010) algorithm. There are additional transient signals visible in Fig. 3.3 that were not detected at six or more stations. Some of these are likely to be earthquakes and may be detectable by reducing the number of stations with corresponding signals to a smaller (and spatially more compact) subset and/or using other indicators. However, our present purpose is to demonstrate the utility of the method and we leave efforts to detect additional smaller events to future work. Most of the newly-detected events are concentrated in areas where there was already some seismicity detected by the SCSN, suggesting earthquake sequences with considerably more events than known before. This significant increase in the number of events was obtained using only 13 of the 41 SCSN stations in this region. Table 3.2: Results of testing the number of indicators needed for PPTS stabilization. # Indicators # Detected % True % False # Missed 2 23 4 96 8 3 53 9 91 4 4 35 17 83 3 5 22 27 73 3 6 16 50 50 1 7 14 57 43 1 8 9 100 0 0 9 9 100 0 0 10 9 100 0 0 3.5 Discussion Our detection algorithm combines statistically normalized information from multiple indicators to achieve a more robust indicator variable. It offers considerable flexibility in terms of which indicators to employ and what parameters to use them with. In this study we used exclusively STA/LTA detectors, but other indicator variables that produce a time series can also be used. These include polarization information such as pseudo-incidence angle (e.g. Rosenberger, 2010), or the various types of amplitude-based detectors summarized by Withers et al. (1998). While we focused on local earthquake detection in the San Jacinto fault zone area, other users may select indicator variables that are more relevant to different regions and scales, or signals. Overlapping events present a challenge for the PPTS and may result in a single continuous signal, 65 but some earthquakes buried in coda are identifiable. The method can work with larger earthquakes but requires appropriately tuned parameters. The algorithm is based on exploiting the redundancy in the peaks of a set of indicator variables at common time points, which should occur only for unusual portions of the seismograms. With enough indicator variables, the likelihood of them all producing large values is low. This can be used to produce a time series that has much more robust peaks than those from any one indicator variable. The core of the method uses multiplication of a set of PPTS at each time step, each varying in the range [0, 1]. This common normalized range facilitates separation of the background noise from earthquake arrivals by multiplication of the individual PPTS. We also tested adding the PPTS, but this does not yield clear results even with various weighting schemes (which we try to avoid). The idea of multiplication was motivated by the set theory intersection operation, which describes the likelihood of a set of outcomes being true simultaneously. We use the term pseudo- probability, because there is an obvious correlation neglected between these indicators as they all measure a quantity derived from the same time series. This is not a problem here because we are not interested in actual probabilities or forecasting, but rather a quantity that measures the simultaneous coherency of a set of indicators. The choice of running the algorithm on one hour intervals is due to several reasons. The first is that the ambient noise changes on many time scales, and we are interested in a window of time that is representative of the noise around the time of earthquake arrivals. The algorithm employs one hour intervals to obtain a sample of STA/LTA values, which are used to characterize the state of the wavefield over that hour via the ECDF. Each value of the original STA/LTA time series is then transformed into its percentile relative to all other values. By using the window of time surrounding each value, we get a relevant reference distribution that is representative of each examined time. If the one hour window is extended significantly, it may ultimately lead to increased sensitivity in the PPTS because some of the values generating the ECDF may be too different than those observed at a given moment. A time window that is too narrow could lead to less sensitive results because the ECDF may be biased towards the most recent values. The choice of one hour seems to be a good compromise, but changing the duration around that value has little effects on the results. As demonstrated in section 4, using the technique with a simple detection scheme on a small network near the SJFZ detects over 3 times as many earthquakes as listed in the SCSN catalog. 66 The method may perform even better by combining different types of indicators, instead of just STA/LTA detectors with varying parameters. Phase picking is easier if an earthquake has already been detected near a particular point in time. The PPTS of individual components and their time derivatives can be useful for automatic picking of P and S phases because the phases generally arrive at the same relative location within the PPTS peak for each event. Improving earthquake detection algorithms is essential for utilizing the full potential of recorded data. Reliable detections of smaller events will lower the completeness magnitude of seismicity and focal mechanism catalogs, and may reveal key new information on the physics of faulting. Adaptation of the algorithm for real time detection may be possible with recursive updating techniques 67 4. An algorithm for automated identification of fault zone trapped waves (Ross and Ben-Zion, 2015) 4.1 Introduction Fault zone structures can produce, in addition to the usual P and S body waves, head waves that propagate along bimaterial interfaces and trapped waves that are associated with resonance modes in a low velocity fault zone waveguide (e.g., Ben-Zion and Aki, 1990). Such waves have been observed along subduction zones (e.g., Fukao et al., 1983; Shapiro et al., 2000; Hong and Kennett, 2003), large strike-slip faults (e.g., Li et al., 1990; Ben-Zion et al., 2003; McGuire and Ben-Zion, 2005) and normal faults (e.g. Rovelli et al., 2002; Calderoni et al., 2012; Avallone et al., 2014). The observation and modeling of head and trapped waves provide high-resolution information on the internal structure of fault zones, motivating deployments of dense seismic arrays across faults (e.g., Li et al., 1994; Mamada et al., 2004; Ozakin et al., 2012). A growing number of dense deployments near faults are recording vast amounts of data (e.g., Bulut et al., 2009; Kurzon et al., 2014; Ben-Zion et al., 2015). Systematic analyses of such data require automatic algorithms for detecting various phases. These types of algorithms can sift through large data sets and add objectivity to the results. Ross and Ben-Zion (2014) developed an automatic algorithm for detecting and phase picking of P and S body waves and fault zone head waves. In the present paper we introduce an algorithm for automatic identification of fault zone trapped waves (FZTW). The method searches for a set of features characteristic of FZTW in time windows identified by the S-wave picking algorithm of Ross and Ben-Zion (2014). These features are examined in seismograms generated by many events and recorded at near-fault stations. Frequent occurrence of these features at given stations are identified as reflecting candidate FZTW. The technique is independent of array geometry and station position, and it assumes no prior knowledge about which stations (if any) are located inside low velocity fault zone layers. The algorithm is verified and demonstrated using data recorded by two dense linear arrays that cross large fault zones in southern California. 4.2 Data and Pre-processing The data used in this study were recorded by two linear arrays deployed across the surface rupture of the 1992 Landers earthquake in the Eastern California Shear Zone and across the San Jacinto Fault Zone (SJFZ). The examined data from the Landers array were recorded during 68 Figure 4.1. Location map for data (circles and stars) recorded during October 14-17, 1992, by a 22-station linear array (blue triangle and inset) across the rupture zone of the 1992 Landers M 7.3 earthquake. The bottom plot shows projection on a vertical cross section along the line A-A’. The locations of the 207 earthquakes shown are from Peng et al. (2003). Data generated by these events were analyzed for trapped waves. Events with at least one detection at the array are colored red. A total of 70 events were flagged by the algorithm. The purple, green, yellow and orange stars indicate location of events used in Figures 4.3, 4.4, 4.6 and 4.9, respectively. The shading indicates the regional topography, and black lines denote faults. 69 Figure 4.2. Location map for data (circles and stars) recorded during 2013 by a 9-station linear array (blue triangle and inset) across the San Jacinto fault zone at Jackass Flats. The bottom plot shows projection on a vertical cross section along the line A-A’. The 5203 earthquakes shown were detected and located by the ANZA network. Data generated by these events are analyzed for trapped waves. Events with at least one detection at the array are colored red. A total of 540 events were flagged by the algorithm. The yellow and orange stars indicate location of events used in Figures 4.7 and 4.10, respectively. The shading indicates the regional topography, and black lines denote faults. 70 October 14-17, 1992, on 22 three-component short-period seismometers (Figure 4.1). The instruments were spaced 25 m apart within 200 m of the surface rupture and 50-100 m further away. A total of 207 events recorded by the Landers array in that period and located by Peng et al. (2003) are considered in the automatic detection analysis performed below. Their analysis indicated that ~6-7 of the stations were located inside the fault damage zone (stations C00-E06). The examined data from the SJFZ were recorded at Jackass Flats (JF) during January 1, 2013 to December 31, 2013 by 9 three-component broad-band seismometers across the surface expression of the Clark fault (Figure 4.2). This array has instrument spacing of 25-50 m and the events were detected and located by the ANZA network (eqinfo.ucsd.edu). Results from this study, which are described in detail subsequently, indicate that stations JF00-JFS3 are likely to be inside or on the edge of the fault damage zone. A total of 5203 earthquakes located within 100 km of the JF array are considered in the analysis below. Basic pre-processing was performed on the data for both arrays. The mean and trend were removed and the horizontal components of motion were rotated to a fault-parallel, fault-normal coordinate system. This is because typical Love-type trapped waves energy exists preferentially in a component of motion oriented parallel to the fault plane (e.g. Ben-Zion and Aki, 1990; Peng et al., 2003; Lewis and Ben-Zion, 2010). For the Landers array, the data are used in a raw unfiltered form to demonstrate the general applicability of the method. For the JF array dataset, which is over 20x larger, the instrument response was removed and a Butterworth band-pass filter was applied between 2-20 Hz to eliminate instrument noise. 4.3 Methodology Figure 4.3a shows an example set of waveforms recorded by the Landers array containing FZTW at stations C00 through E04. The seismograms containing FZTW (black ellipse) have longer period, higher amplitude phases arriving shortly after the S-wave that are about half a second in duration. These features are common among FZTW, but the waveform shape, arrival time, and other characteristics can vary depending on the properties of the trapping structure (e.g., Ben-Zion, 1998; Jahnke et al., 2002). Our methodology for automatic detection of FZTW can be generally stated as a procedure for detecting outliers of various waveform features across a fault zone array that indicate common aspects of FZTW. 71 Figure 4.3. Example waveforms with FZTW (enclosed by ellipse) in the Landers array generated by the event denoted by the purple star in Fig. 4.1. a) Fault-parallel component velocity seismograms across the array. FZTW are observed at stations C00-E04 shortly after the S-wave arrival (vertical black line; automatic pick). A 1.0 sec window (solid horizontal line) starting at the S pick is used to calculate predominant period, wave energy, arrival delay, and relative peak strength at each station in the array. A 6.0 sec window (dashed horizontal line) centered on the median S pick is used to calculate energy and identify records with possible site amplification. b) Plot of the Y statistic values for each of the features. Stations E01-E05 have Y statistics for all of the 1 sec window measurements that are 1 deviation above the median or more. As shown by Ben-Zion and Aki (1990) and Lewis and Ben-Zion (2010), FZTW start from the end of the S body wave, or slightly behind in structures with overall velocity contrast across the fault, so a reliable estimate of the S arrival is crucial. We therefore first run the S-wave picking 72 algorithm of Ross and Ben-Zion (2014) on seismograms at each station of a given array for each event. The picking method first employs polarization analysis to remove P-wave energy from the seismograms. Then STA/LTA and kurtosis detectors are run in tandem to lock on a well-defined S-arrival. Picks are made on both horizontal components (N and E), if possible, and the median value is calculated over all picks and components to get a single, robust estimate of the S-arrival at the array. This provides a good approximation for relatively short arrays (length < 1 km) of the type analyzed here. For longer arrays, it may be desirable to use individual picks at each station instead due to a potential move-out across the array. For the Landers and JF arrays, the method has been tested and works both ways, but using individual picks is in general less reliable since the picks vary in accuracy. After the picks are made, the traces are restored to the original raw form for detecting trapped waves, since the polarization process used during picking is unnecessary for FZTW detection. The algorithm of Ross and Ben-Zion (2014) requires a minimum signal to noise ratio (SNR) of 5 to make a valid S pick. We require further that more than half of the stations have S picks in order to calculate the median; otherwise, the event is skipped. The use of the median is to prevent outliers from dominating the result. Subsequent references to S picks in the paper denote the median S pick across the array. The S pick for a given event is used to define the start of a 1.0 sec long window (Fig. 4.3a). The duration of the trapped waves windows is expected to increase with the width of the trapping structure, propagation distance within the fault zone layer, and velocity contrast with the bounding rocks (e.g. Ben-Zion, 1998). Previous studies indicate that observed FZTW in large strike-slip faults associated with low velocity zones with widths on the order of 100 m and velocity reduction of 30-50% generally have a duration of about 1 sec or less (e.g. Li et al., 1994; Ben-Zion et al., 2003; Lewis et al., 2005; Mizuno et al., 2008). We tested window lengths of 1.0, 1.5, and 2.0 sec and found that, for the Landers rupture zone and JF site at the SJFZ, a 1.0 sec window identified the most trapped waves at the most reliable rate. This window is referred to below as the trapped wave window. Additional information about the window length is given in the discussion section. The trapped wave window is used to look for four features common to FZTW at each station. These characteristics are: 1) longer period energy than S body waves expected for a resonance mode in a waveguide, 2) typically larger amplitudes following the S-wave arrival (e.g. Fig. 4.3a), 3) arrival time delay relative to S-wave, and 4) peak amplitude of trapped waves window relative to the average amplitude (termed “relative peak”). Other characteristics such as dispersion are 73 weaker second-order features for shallow trapping structures (e.g. Ben-Zion et al., 2003), and as such are not considered here. The features we use are analyzed in the time domain, which is equivalent to but more computationally efficient than Fourier domain techniques (e.g. Kanamori, 2005). To measure the approximate period of the trapped waves window, we use the recursive predominant period algorithm of Nakamura (1988), i i P i D X T π 2 = , (1a) 2 1 i i i x X X + = − α , (1b) 2 1 i i i t x D D ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ + = − α , (1c) where P i T is the predominant period, i x is the ground velocity at time index i, and α is a damping constant equal to 0.999. These equations are recursively calculated for each sample over the trapped waves window to obtain the final measurement. Equations (1a)-(1c) are used to perform a single measurement of the predominant period for the trapped waves window in each examined seismogram. The second feature is an estimate of the seismic energy in the trapped waves window. This is calculated simply as the sum of the squares of each data point in the time window. The third feature (relative peak) measures how strong the peak of the trapped waves window is relative to the entire window. It is calculated by taking the ratio between the absolute peak amplitude of the trapped waves window to the average amplitude of the window. This feature is small numerically when the entire trapped waves window has similar amplitude, and large when the peak value is considerably greater than the average value. If a FZTW is indeed present in the window, its amplitude should be generally larger than the surrounding S-wave, leading to a large relative peak. The fourth feature measured is the time delay between the peak of the trapped waves window and the S-wave arrival. This property of FZTW can change due to a velocity contrast across the fault (Lewis and Ben-Zion, 2010) and thus may not be present in all cases. As a result, it is used in our algorithm but given relatively low weight. More details are given about this feature in a subsequent section. 74 We found that using jointly these four characteristics is useful for identifying FZTW candidate waveforms, but similar features can be produced in some cases by other propagation and site effects. For example, highly localized site amplification can occur within an array, regularly yielding large energy measurements for a given site relative to the rest of the array. To identify situations where this occurs, a fifth measurement is made by calculating the energy in a longer time interval around the S arrival. Here we use a window duration of 6 sec, centered on the S arrival (see Fig. 4.3a). The long duration of this window, combined with the fact that it is not restricted to just the S-wave, makes this a suitable quantity for identifying site amplification. This duration is also considerably longer than the expected duration of FZTW, so that if one is present it will not dominate the energy measurements in this window. The precise length of the window is not too important as long as the total duration is a significant portion of the entire seismogram and much longer than the duration of a trapped waves group. While the first four features described tend to indicate the likelihood of FZTW, the fifth one helps the algorithm to suppress erroneous cases with likely site amplification. With the five features measured at each station for a given event, the problem becomes one of outlier detection to identify events with candidate FZTW. A common statistical method of outlier identification is to test whether a measurement is sufficiently far from the mean of a set of measurements, after normalizing by the sample’s dispersion. This method is ineffective when the sample size is small and some expected values are statistical outliers themselves. Since the sample size here (the number of stations in a given array) is generally small, we must include all the stations available for each array in the calculations. Having possible unknown outliers decreases the detection sensitivity by biasing the mean and standard deviation. To address this limitation, we use a quantile based variation of the outlier identification scheme that uses instead the median and median-absolute-deviation (MAD). The MAD for the jth feature is defined as, ( ) ( ) j j i j X X MAD median median − = for 5 ,..., 2 , 1 ; ..., , 2 , 1 = = j N i (2) where X is the feature type and N is the total number of stations in the array. The MAD is intuitively a measure of dispersion that incorporates how far each measurement is from the median (rather than the mean) in an absolute value sense (e.g. Hampel, 1974). It is thus insensitive to outliers that may be present, just like the median or interquartile range. For a sample with a large number of measurements, the median and MAD would only need to be calculated once per event for each feature. However, since the number of stations may generally 75 be small, and can still be affected if multiple stations are outliers, the median and MAD are calculated independently at each station. This is done by removing the station in question from these calculations, resulting for each feature in a slightly different MAD and median per station. Using these parameters the normalized statistic Y is calculated as, ( ) j i j i j i j i MAD X X Y median − = (3) where X i j is the jth feature at the ith station. Y i j represents the number of deviations from the median for the jth feature at the ith station. This provides a simple convenient method for normalizing each measurement so that the same scale can be used for all features and events. The statistic Y is calculated for each of the five features at each station in the array, yielding a total of 5N Y-values. Figure 4.3b illustrates the Y statistic values for each feature, at each station, for the waveform set shown in Figure 4.3a. For stations C00-E04, a clear bump is visible for all features indicating that the group of stations are collectively producing features that are above the median. With the measurements transformed into a normalized form, detection thresholds can be easily set. We systematically tested many different combinations of thresholds for each of the five features, and found the set given in Table 4.1 to produce the most FZTW detections at the most reliable rate. Specific details about how these parameters were chosen are given in the next section. To make a detection, the four features measured on the 1 sec window are all required to be above the listed threshold, whereas the energy measurement over the 6 sec window is required to be below the listed threshold. The reason for the two different windows is now explicitly clear: the 1 sec window is designed to target FZTW characteristics in a localized window after the S-arrival, while the 6 sec window enables the method to avoid detection of site amplification that may be present in a particular trace. All records that satisfy the above criteria are flagged for further investigation. Table 4.1. Y-statistic thresholds used for trapped wave detection. Feature type Threshold Energy (1 sec) 0.75 Dominant period 1.25 Relative peak 0.0 Arrival delay 0.0 Energy (6 sec) 2.75 76 The thresholds given in Table 4.1 may need to be changed to fit different datasets or applications. For example, since arrival delay for the trapped waves group is not always present, this parameter can be varied based on whether a velocity contrast is expected across the fault. For the example in Figure 4.3, a detection is made on station E04. As several more stations would have been flagged if the 6 sec energy requirement was relaxed or not imposed, it is clear that there is a tradeoff between detection rate and ensuring that false detections by site amplification are kept to a minimum. We note that trapping structures and other low velocity fault zone layers generally produce themselves some site amplifications (e.g. Ben-Zion et al., 2003; Kurzon et al., 2014). This can be observed in the higher amplitude P-waves for stations C00-E05 in Figure 4.3a, which explains why the energy measurements in the 6 sec window are higher for these stations in Figure 4.3b. Nevertheless, the 6 sec energy requirement prevents false detections due to other effects (e.g. topography and sedimentary basin) that may locally amplify the ground motion more strongly. One example of such non FZTW amplification is shown in Figure 4.4. Figure 4.4. Example waveforms (fault-parallel component velocity seismograms) with large site amplification at station W07 generated by the event denoted by the green star in Fig. 4.1. The vertical black line indicates the median S pick across the array. Station W07 was found to have a Y statistic of 3.53 for the energy measurement in the 6 second window. 77 4.4 Application to data sets at the Landers and JF arrays We now demonstrate the method on the raw unfiltered dataset recorded by the array across the Landers rupture zone. In total, 190 events met all of the quality criteria (minimum signal-to-noise ratio and minimum number of S-picks) described in the previous section and were passed to the detector for investigation. First, we used the observations of Peng et al. (2003), who found that FZTW only were seen on stations W01-E06, to test all possible combinations of Y-statistic thresholds from 0 to 3.0 deviations above the median, in 0.25 deviation increments. The optimal parameter set, given in Table 4.1, was chosen by requiring the success rate to be greater than 90% and then maximizing the number of detections on stations W01-E06. To give a sense for how sensitive the method is to the actual thresholds used, nearly 200,000 different combinations were tested during this process and yet the mean success rate over all of them was 88.3%. As a result, we conclude that the method is quite robust with regards to the exact thresholds used. Figure 4.5. Results of the automatic detection algorithm applied to the Landers array data set. A total of 108 records were flagged for by the detector as likely containing trapped waves. Of these detections, 8.3% were found to be false based on visual inspection. Nearly all detections occur between stations C00-E06. Using the optimized set of thresholds provided in Table 4.1, we examined the corresponding detection results in greater detail. Figure 4.5 displays a histogram with the number of detections at each station deemed to be true, as well as the number of false detections. We visually inspected 78 each detection in both the frequency and time domains to thoroughly examine which ones were false. In total, 70 events out of 190, which corresponds to 36.8%, were flagged as FZTW candidates by the detector. The detected events are shown in map view and vertical cross section as red circles in Figure 4.1 on the background of other examined events (black circles). From the set of 108 seismograms flagged, 9 were deemed to be false detections (~8%). Of the 99 detections we consider to be true, ~60% were also identified by Peng et al. (2003) as having a large energy ratio between stations inside and outside the damage zone. It can be seen that the events with FZTW are broadly distributed in space rather than being highly localized near the rupture zone (Fig. 4.1). A deep low velocity fault zone layer that reaches the bottom of the seismogenic zone produces trapped waves only from events within the fault zone, since waves impinging on the fault zone from the outside are largely reflected away (e.g., Ben-Zion and Aki, 1990; Igel et al., 1997; Jahnke et al., 2002). In contrast, a shallow low velocity layer can generate trapped waves from deeper regional events that inject wave energy into the fault zone from below (Ben-Zion et al., 2003; Peng et al., 2003; Fohrmann et al., 2004). The detection results in Fig. 4.1 therefore indicate that the Landers trapping structure extends primarily over the top few km of the crust. One example set of FZTW generated by a given event and detected by the algorithm was given in Figure 4.3. Another example is shown in Figure 4.6; the stations with FZTW detections are denoted by a black star and the event location is marked by the yellow star in Figure 4.1. These detections, as well as nearly all of the ones summarized in Figure 4.5, are in close agreement with the general observations of FZTW detected by Peng et al. (2003). They found that FZTW in the Landers array were primarily observable at stations C00-E06. Our detection algorithm, like most other detection methods in seismology (e.g. Nippress et al., 2010; Ross and Ben-Zion, 2014; Langet et al., 2014) is based on thresholds being met, and thus will miss FZTW with weaker attributes, especially for receivers near the edges of the trapping zone. This is clearly demonstrated in Figure 4.6 by the fact that FZTW are visible at stations C00 through E05, but were only identified by the algorithm at stations E04 and E05. The method is designed to detect for each examined event the clearest FZTW at a few stations, rather than detect weak candidate FZTW on many array stations, since the latter can lead to many false detections. Next, the method is applied to the JF array dataset across the SJFZ recorded during 2013. This dataset consists of several thousand events, and with 9 stations in the array (Fig. 4.2) is a very large volume of waveforms to sort through manually in a search for FZTW. Due to the large number of 79 source-receiver combinations it is expected that there will also be a significant number of records with poor overall quality. To limit the number of seismograms with large noise and signal outside Figure 4.6. Fault-parallel component velocity seismograms with trapped waves detected at the Landers array generated by the event denoted by the yellow star in Fig. 4.1. The stations with detections are indicated by black stars and the vertical black line indicates the median S pick across the array. Visual inspection indicates candidate FZTW at stations C00-E05 (enclosed by ellipse). 80 the frequency band of interest, we first apply to this data set a 2-20 Hz band-pass filter. Figure 4.7 provides an example set of waveforms with FZTW detected at stations of the JF array. The location of the event generating this data set is marked with a yellow star in Figure 4.2. It can be seen that long period, high amplitude FZTW about 0.5 sec in duration start about half a second after the S arrival. The phases are only present clearly on JF00, JFS1, JFS2, and JFS3. These FZTW look visually overall similar to the ones in the Landers rupture zone (Figures 4.3 and 4.6), and are also overall similar to FZTW observed across a different location on the SJFZ (Lewis et al., 2005) as well as trapped waves at the Karadere-Duzce section of the North Anatolian fault (Ben-Zion et al., 2003) and the Parkfield section of the San Andrea fault (Li et al., 1990; Lewis and Ben-Zion, 2010). Figure 4.8. Results of applying the automatic detection algorithm to the JF array data set. A total of 582 records were flagged by the method with roughly 90% of them concentrated between JF00- JFS3. FZTW are generally observed most strongly in the waveforms at station JFS1, explaining the large number of detections made by the algorithm. The detections at other stations are likely associated with various occasional site effects. The ability of the algorithm to identify trapped waves for fault zone structures with different properties (here the Landers rupture zone and the JF site of the SJFZ) is important, as it allows the algorithm to be applied to different datasets with limited changing of parameters. Figure 4.8 summarizes the number of detections in each station of the JF array for the examined dataset. A total of 2255 events met all of the quality criteria (minimum signal-to-noise ratio and minimum 81 number of picks) needed for processing. In contrast to the Landers dataset, we have no previous information on FZTW in this region to compare with. In total, 526 of the 582 detections (90.4%) are concentrated at stations JF00-JFS3 as in the example shown in Figure 4.7. These stations are also believed to span the extent of the damage zone. Station JFS1 is found to have the most detections by far; this is related to the fact that this station frequently records the most pronounced trapped waves across the array. The locations of all events generating detected FZTW at the JF array are shown in Figure 4.2 (red dots). As with the Landers data set, the locations of events producing candidate FZTW at the JF array are distributed broadly, pointing again to a trapping structure that exists primarily over the top few km of the crust. While it appears from Figure 4.2 that there are a few large clusters of detections, this is visually misleading as a substantial portion of all events are concentrated in these clusters. We have inspected visually about 300 randomly selected detections at the JF array and estimate our false detection rate on this dataset to be around 10%, which is in agreement with the results from the Landers array. None of the detections on the north side were found to contain clear trapped waves and are thus considered false detections. Visual inspection of the entire dataset is the subject of future work; in particular the events with weakly generated FZTW require detailed examination to properly confirm. We further tested the method on raw JF data and found the general properties of the results to be unchanged. However, with raw data there are a larger number of additional detections for stations on the north side of the fault. An important aspect of detection algorithms is to provide information on false detections. For the Landers data set, we examined visually all the automatic detections and also compared our detections to the events identified by Peng et al. (2003) as having a large energy ratio between on- and off-fault stations. Figure 4.9 shows an example of a detection at station W07 of the Landers array determined by visual inspection to be false. For this particular event, the 5 features just barely met the thresholds necessary for flagging the record at station W07 and no detections were made at other stations. An example of a false detection on the JF array is shown in Figure 4.10. Here, the false detection occurred for station JFN2, and in this case the set of 5 features also just barely met all of the required thresholds. However, as evidenced by the histogram in Figure 4.8, detections on the north side of the JF array are uncommon. For these types of isolated cases, highly- localized site effects are the most probable cause of longer period S-waves. 82 Figure 4.9. Fault-parallel component velocity seismograms with falsely detected trapped waves at the Landers array generated by the event denoted by the orange star in Fig. 4.1. The station with the false detection is denoted by a star and the vertical black line indicates the median S pick across the array. The four features measured on the 1 sec window are just barely above the required thresholds and are not observed on nearby stations like in Figures 4.3 or 4.6. 83 Figure 4.10. Fault-parallel component velocity seismograms with falsely detected trapped waves at the JF array generated by the event denoted by the orange star in Fig. 4.2. The station with the false detection is denoted by a star and the vertical black line indicates the median S pick across the array. The four features measured on the 1 sec window are barely above the required thresholds at station JFN2 but visual inspection of the waveform at this and neighboring stations leads to rejecting the detection. 84 4.5 Discussion FZTW provide high-resolution information on internal components of fault zone structures, but identification of records containing FZTW has traditionally been a tedious task. While several automatic techniques have been used previously to aid in their identification, they were limited by the time needed to manually pick S phases and assumptions made. Peng et al. (2003) calculated energy ratios between different stations and Lewis et al. (2005) used a variant technique involving energy ratios over a particular bandwidth between different stations. The methods in these studies did not formally classify events or recordings as to whether or not they contained trapped waves, which is a key ingredient of a detection algorithm. With classification also comes success rates that provide an estimate of how reliable the method is. In this work we generalize and combine elements from these past techniques with the automatic S picking algorithm of Ross and Ben-Zion (2014), and use the method to classify each individual recording. The developed algorithm examines each seismogram in comparison to records at other stations of the array, and determines whether expected features of FZTW in the seismogram are statistical outliers of the array. Our detection algorithm adds to the energy estimates of previous techniques explicit measurements of the predominant period, relative peak strength, arrival delay of the considered phase, and possible site amplification. These extra features provide important additional metrics for comparing the wavefields at different stations. By using the five features together, lower thresholds can be used, leading to more detections at a similar degree of reliability. Another challenge for identification of FZTW is that other wave propagation phenomena can occasionally produce (e.g. for some source mechanisms and source-receiver geometries) highly localized site effects that resemble FZTW. As seen in the histograms of Figures 4.5 and 4.8, we identify multiple occurrences of such signals (e.g. station JFN2 in Fig. 4.7). However, the systematic analysis performed by our algorithm recognizes that they are not generated by many events, whereas resonance modes associated with waveguides below the stations are expected to be frequently produced. One of the most serious limitations of past studies was the need for a trained analyst to pick S wave arrivals to define the window for the FZTW search. The data volume collected now around the world is vast and rapidly increasing, making it unfeasible to sort through the data manually. The incorporation of automatic S picks in our algorithm allows for rapid processing of years of data, as demonstrated for the JF array. 85 The application of the automated algorithm should be followed by careful additional analyses to confirm that the identified phases are FZTW. The primary goal of the algorithm is to systematically and objectively analyze large datasets and reduce them to a dramatically smaller set of records that are likely to contain FZTW at certain stations. The method helps to identify which stations are worth examining first, which events are likely to be strong FZTW candidates, and which stations may have occasional signals with features similar to those of trapped waves (Figures 4.5 and 4.8). Additional follow up analyses of candidate FZTW include examining the moveout between the S-wave and trapped waves group, analyzing the dispersion and spectral properties of candidate FZTW, and performing waveform inversions for fault zone properties (e.g., Ben-Zion et al., 2003; Peng et al., 2003; Lewis et al., 2005). Identification of the highest-quality FZTW candidates within a set of detections is often an important first step for modeling waveforms to invert for structural properties of fault zones. We found the number of detections made per event to be generally a reliable indicator of FZTW quality. By sorting the events based on the number of stations with detections, one can start with the events having most detections and work backwards. For both the Landers and JF arrays, events with 2-3 detections typically had FZTW clearly visible at many stations, but this number will vary depending on the number of stations within the damage zone. Nearly all detection methods, including the one presented here, involve the use of some parameters. In the current method there are 5 detection thresholds and two window durations that may be adjusted for different applications. The values used here for these parameters (Table 4.1) were selected to provide many detections while keeping the false detection rate relatively low. We used the same parameter values for both the Landers and JF arrays and found these to work well in both cases. If the thresholds for the 1 second window features are lowered, additional weak candidate trapped waves will be recognized at the cost of more false detections. The method has several limitations that should be noted. First, the number of detections that can be made per event depends on the number of stations both inside and outside the fault zone trapping structure. For future deployments of linear arrays with fixed number of stations, our results imply that fewer arrays with more stations may be better than more arrays with fewer stations for the purpose of detecting FZTW. Our testing suggests that less than ~30% of the stations should be located inside the damage zone to facilitate identifying stations with trapped waves as outliers. Another limitation of the method is the accuracy of the automatic S picks. The algorithm 86 of Ross and Ben-Zion (2014) is able to make S-picks to within 0.25 sec of the analyst pick roughly 90% of the time. If the used trapped wave window length is 1 sec, as in this work, this allows some room for error in the S pick, and should not be an issue unless the FZTW are arriving quite late (~1 sec or more). Further, the S-wave picks tend to be late, rather than early, which helps lessen this problem. Using the median S pick and requiring a majority of stations to have valid picks improves the stability of the detections. Most of the false detections in the examined data at the Landers and JF arrays were in fact due to mispicks of S arrivals. However, the overall false detection rate is reasonably low enough to produce overall robust results. The mean S pick across the array could alternatively be used instead of the median if outliers are discarded beforehand. In cases with strong overall velocity contrast across the fault, the variability of S picks at different stations as well as the delay between the S arrival and FZTW increase (e.g. Lewis and Ben-Zion, 2010). In such cases it may be better to use S picks at individual stations and to lengthen the trapped wave window. The detection algorithm is not very sensitive to the precise number of array stations. We have tested the method on arrays with the number of stations varying from 7 to 21 and find similar success rates. The decision to rotate seismograms to fault-parallel component of motion for performing the analysis may be viewed as a pre-processing step. Love-type trapped waves involve particle motion parallel to the fault zone layer (Ben-Zion and Aki, 1990; Ben-Zion, 1998), so FZTW are expected to be strongest on the fault-parallel (and for vertical fault zones to some extent also vertical) component seismograms. For fault zones with known geometry that is approximately linear, it is straight forward to rotate seismograms to the fault-parallel direction. For fault zones with curvature, such as the Landers rupture zone (Fig. 4.1), the choice of a trend is more difficult. We found, however, that using the north component of motion for all of the Landers data yielded results that were similar to an average fault parallel component by drawing a line through the rupture end points. Alternatively, the vertical component seismograms can be used (e.g. Li et al., 1994). In any case, the choice of rotation is not critical and should be seen more as a technique for enhancing the contrast between trapped waves and S-waves rather than a strict requirement. The developed automated method for identification of fault zone trapped waves was tested in this work in the context of S-type FZTW, but it could in principle be also applied to P-type FZTW (Ellsworth and Malin, 2011) with appropriate changes of parameters. Our algorithm was tested 87 partially on several other linear arrays deployed across the SJFZ, and seems to produce satisfactory results on detection of S-type FZTW without changing the parameters (Qiu et al., 2015; Share et al., 2015). Work on comprehensive detection of FZTW in data of multiple dense arrays across the SJFZ and modeling the data for high-resolution information on the internal fault zone structure is in progress. 4.6 Conclusions We developed a method for automated identification of fault zone trapped waves produced by constructive interference in low-velocity fault zone layers. The method identifies the S-wave group using an automatic S-wave arrival picking algorithm (Ross and Ben-Zion, 2014), and computes a set of features that are representative of trapped waves on each record for a given event. The features are compared among all stations across a fault zone array, and are used to identify statistical outliers. Records that satisfy a set of criteria are flagged for further visual examination. A set of detection thresholds was optimized over thousands of different combinations to obtain the best results in the examined data sets. Applying the method to datasets in the Eastern California Shear Zone and San Jacinto Fault Zone yielded clear detections at subsets of array stations consistent with manual analyses. The method may be applied to data sets at other locations and other wave types with some changes of parameters. 88 5. Towards reliable automated estimates of earthquake source properties from body wave spectra (Ross and Ben-Zion, 2016) 5.1 Introduction Recorded seismograms are associated with convolutions of source, propagation-path and instrument effects. Estimating source properties from observations requires removal of the other effects that influence the data. This is often done by deconvolving seismograms of the event targeted for analysis with an empirical Green’s function (EGF) given by seismograms of a suitable small event located close to the hypocenter of the target event (e.g. Berckhemer, 1962; Mueller, 1985; Hutchings and Wu, 1990; Hough and Dreger, 1995). The deconvolution and subsequent analyses can be done in either the time domain or frequency domain. In the present paper we discuss techniques for deriving earthquake source properties by analysis of earthquake spectra. The goal is to reliably estimate a set of source properties with an automatic procedure from the spectral ratios of target over EGF events. These include the scalar potency/moment associated with the zero frequency asymptote, strain/stress drop involving the corner frequency between the spectral level and high frequency decay, radiated seismic energy estimated by integrating the source spectrum, directivity associated with azimuthal variations of source spectra, and apparent stress given by the ratio of the radiated energy over the potency (e.g., Aki, 1966, Wyss and Brune, 1968; Brune, 1970, Ben-Menahem, 1961; Madariaga, 1976; Ben-Zion, 2003). Deriving reliable estimates of earthquake source properties from spectra has long been a difficult problem, and observational results of stress drops and other spectral-based properties are associated with large scatter (e.g., Shearer et al., 2006; Goebel et al., 2015). A key question is how much of the scatter is real and how much results from errors in techniques and lack of knowledge such as mixing different event populations. There are partially conflicting statistical and physical requirements for resolving genuine properties of earthquakes (and other fault processes). Robust statistical estimates require that the data set is sufficiently large to suppress statistical fluctuations. This usually leads to analysis of results associated with large spatial domains. However, using overly large regions can mix different classes of events and increase the scatter (e.g., Ben-Zion, 2008). Many studies assume that earthquake properties are essentially the same everywhere, apart from statistical fluctuations, and average data from large domains. However, this approach can 89 suppress recognition of persistent variations of source processes related to different properties of fault zones and the crust, and it may increase the scatter. One example of a fault-related source behavior is systematic directivity of earthquake ruptures on bimaterial faults that separate different crustal blocks (e.g., Weertman, 1980; Ben-Zion and Andrews, 1998; Ampuero and Ben-Zion, 2008; Brietzke et al., 2009). Recent observational studies demonstrated the existence of systematic directivity of earthquakes on various fault sections even for small events (e.g., Lengliné and Got, 2011; Kane et al., 2013a; Kurzon et al., 2014; Calderoni et al., 2015). Unrecognized directivity can lead to scatter and biases in stress drops (and other properties) inferred from corner frequencies at stations affected by the directivity (e.g. Calderoni et al., 2013). Directivity of moderate and large events can significantly increase the amplitude of ground motion at stations in the forward direction, especially when coupled with structural effects (e.g. Olsen et al., 2006; Avallone et al., 2014), so it is important to clarify the possible existence of persistent earthquake directivity on given structures. Another fault-specific process that can change the amount and frequency content of seismic radiation is co-seismic rock damage in source volumes (Ben-Zion and Ampuero, 2009). This can be important for earthquakes in geometrically complex regions that break rocks (Castro and Ben-Zion, 2013; Kwiatek and Ben-Zion, 2013; Ross et al., 2015), but is less relevant for events on well developed faults that fail essentially by frictional sliding. Additional examples are possible correlations of stress drops with inter- vs. intra-plate environment, faulting style and hypocentral depths (e.g., Kanamori and Anderson, 1975; Scholz et al., 1986; Mori and Frankel, 1990; Kanamori et al., 1993; Shearer et al., 2006] and dependency of source processes on the heat flow in a region (e.g. Ben-Zion and Lyakhovsky, 2006; Enescu et al., 2009; Zaliapin and Ben-Zion, 2013a,b). Another important unresolved question is the scaling of source properties associated with earthquakes of different size. Aki (1967) suggested that earthquakes are scale invariant and that the stress drop, apparent stress and other source quantities are independent of the earthquake size (but may still depend on some properties of faults and the crust as discussed above). Earthquake scaling relations have important implications for many aspects of earthquake physics and are directly related to seismic hazard analysis. For example, if the physical process and amount of radiated energy per unit potency (i.e., unit fault area and slip) during large earthquakes are essentially the same as those of small events, one can study the source properties and effects of the abundant small earthquakes and scale them up to estimate what will happen during the rare large 90 events. The scale-invariance of earthquakes has been the subject of considerable debate (e.g. Abercrombie, 1995; Walter and Mayeda, 1996; Ide and Beroza, 2001; Shearer et al., 2006; Kwiatek et al., 2011; Oth, 2013). Several recent studies extensively analyzed how incomplete recordings and other factors can affect the estimates of derived earthquake source properties. Kaneko and Shearer (2015) performed numerical simulations of earthquakes with different rupture velocities, geometries and degrees of asymmetry, along with different take-off angles and network geometries. They showed that these factors can affect derivations of the mean corner frequency and estimates of stress drops. Kwiatek and Ben-Zion (2013, 2016) analyzed the dependency of source time functions, frequency content and radiated energy on shear vs. tensile faulting, source-receiver geometry, attenuation coefficients and properties of the recording system. Abercrombie (2015) and Kane et al. (2013b) examined different criteria for selecting EGF, including spatial separation, cross-correlation coefficients and verification of a clear source-time function after deconvolution. Abercrombie (2015) advocated a single EGF approach over regional event-station stacking approaches (e.g. Prieto et al. 2004; Shearer et al. 2006) due to the possibility of bias. However, using a single EGF can also produce biases if the EGF event has non-generic source effects such as directivity or small tensile component of faulting. Calderoni et al. (2013, 2015) showed that EGF events having directivity will result in apparent directivity in spectral ratios and erroneous inferences on stress drops of the target events. Additional problems can arise in the course of fitting the obtained spectral results. Real seismic spectra are susceptible to random heterogeneities, directivity, bandwidth limitations, signal to noise variations within the analyzed bandwidth, overlapping of seismic phases and other possible effects. These issues can lead to erroneous estimates of seismic potency, corner frequency and exponent of the spectral falloff, which are obtained by fitting a parameterized model to an observed source spectrum. They can also result in serious parameter tradeoff if the applied fitting process does not include constraints. To address this tradeoff, an omega-square model (Brune, 1970) is often used to help stabilize the fitting of source spectra. However, this constraint may not be correct for stations at all azimuths and for all reasonable source cases (e.g., Kaneko and Shearer, 2015; Kwiatek and Ben-Zion, 2016). Previous studies deriving earthquake source properties from spectra via EGF analysis have generally been limited to relatively small datasets because of user involvement in the various 91 stages. This includes the picking of phase arrivals, adjusting window durations, frequency bands, and possibly confirming a deconvolved pulse in the time domain visually (e.g. Abercrombie, 2015). As a result, these factors can add up to significantly impede the ability to study large seismic datasets in an automated way. One notable exception to this are the types of studies (e.g. Shearer et al., 2006; Prieto et al., 2004; Yang et al., 2009) which perform a large iterative stacking procedure to separate source, propagation and site effects. A shortcoming of this class of techniques is that they average results associated with groups of events, so they may suppress genuine differences of source properties at nearby locations (e.g. on/off-fault events or earthquakes on neighboring faults). In the following sections we describe a general-purpose technique designed for automated analyses of source properties of individual target events from spectra of P and S waves. The technique uses stacked EGFs to average over and minimize possible source effects present in individual EGF spectra. We demonstrate the method’s capability by obtaining spectral ratios at hundreds of stations for a set of small to medium size earthquakes southern California. The spectral ratios are used to estimate the scalar potency/moment, strain and stress drops, radiated energy, apparent stress, ratio of P over S corner frequencies and extent of rupture directivity. The simultaneous derivation of potency/moment and strain/stress drops from P and S waves provides a consistency check on the results. 5.2 Methods and Results Overview The combined source, s, propagation-path, g, and instrument, i, effects on seismograms can be written in the frequency domain as ) ( ) ( ) ( ) ( f i f g f s f u ⋅ ⋅ = . (1) Following Brune (1970), the displacement spectrum of a seismic source is typically parameterized by a corner frequency, f c , low-frequency asymptote, Ω, and spectral falloff, n n c f f f s ) / ( 1 ) ( + Ω = . (2a) A source model with n = 2 is commonly referred to as a Brune model. An alternative version of the displacement spectrum that is better suited for cases with a sharp corner frequency (Boatwright, 1980) is given by, 𝑠 𝑓 = $ %& '/' ) *+ ,/* (2b) 92 As mentioned, to study s we must first remove g(f) and i(f) using the EGF approach or another technique. The EGF method assumes that the corner frequency of the smaller earthquake is much higher than that of the target event, resulting in a delta function response below the higher corner frequency. However, since the EGF f c is finite, it is more appropriate to use instead of (2) a spectral ratio model, r [e.g. Abercrombie, 2015], γ γ γ 1 1 2 2 1 ) / ( 1 ) / ( 1 ) ( ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + Ω Ω = n c n c f f f f f r , (3) where superscripts 1 and 2 refer to the target and EGF events, respectively. The main difference between (2) and (3) is that the latter flattens off as it approaches f c 2 , which is generally observed in spectral ratios. While this adds an additional model parameter, the inclusion of f c 2 helps determining the range of useable frequencies for analyzing source properties of the target event and is generally not very important otherwise. Table 5.1. Target events analyzed in this study Name Date M w Latitude Longitude Depth SCSNID Big Bear 2015-09-16 4.0 34.137 -116.858 9.6 37243591 SJFZ Trifurcation 2013-03-11 4.7 33.507 -116.458 10.2 15296291 Cajon Pass 2015-12-30 4.4 34.191 -117.413 7.0 37507576 Elsinore 2012-08-08 4.05 33.9 -117.791 9.4 15189073 SJFZ Hot Springs 2013-06-28 3.41 33.634 -116.692 14.3 11327386 Data We present the method using data provided by the Southern California Seismic Network (SCSN) (SCEDC, 2013). The method is illustrated using a set of five events (Table 5.1) shown in Figure 5.1 along with the employed seismic stations. The target events are the 2013 M w 4.7 earthquake in the complex trifurcation area of the San Jacinto fault zone (TRIF), the 2015 M w 4.0 Big Bear earthquake (BB), the 2015 M w 4.4 Cajon Pass earthquake (CP), the 2012 M w 4.05 Elsinore fault earthquake (EL) and the 2013 M L 3.41 Hot Springs earthquake on the San Jacinto fault zone (HS). The first 4 events have similar right-lateral strike-slip focal mechanisms, while the HS earthquake has a significant oblique component (Yang et al., 2012). The events were well recorded by hundreds of stations in the SCSN. We chose these events to demonstrate the method, out of many that we have tested in detail, because they exhibit a variety of different source properties. This allows us to show how the method works under a range of circumstances. For 93 example, the TRIF event exhibits strong directivity, while the BB event appears to be essentially a circular rupture. The TRIF event has an average strain/stress drop (compared with other studies), while the BB event has a rather high strain/stress drop. For each target earthquake, EGF events are selected from the relocated seismicity catalog of Hauksson et al. (2012) and are listed in Table 5.S1. The details of how these EGF events are selected and processed are documented in the next section. Figure 5.1. Map of the Southern California plate boundary region. SCSN stations used in this study are denoted by blue triangles. The target events are indicated by stars (Table 5.1). The Trifurcation, Hot Springs, Cajon Pass, Elsinore, and Big Bear earthquakes are colored yellow, pink, orange, brown, and red, respectively. The green triangle is the location of the station used in Figure 5.2. In analysis of each target event, we initially consider all available stations archived by the SCSN in southern California (blue triangles in Figure 5.1). After running automatic phase pickers 94 on all stations (Ross and Ben-Zion, 2014; Ross et al., 2016), quality control procedures are used to automatically determine which records and spectra to retain for subsequence analysis (rather than making arbitrarily choices such as a distance cutoff). For the target and EGF events in the southern California data used here, we only work with high-sample rate broadband and strong- motion records (HH and HN channels), with the exception of HS. All three components are used throughout this study for all events. If both HH and HN channels have data, broadband (or strong- motion) records are used for hypocentral distances larger (or smaller) than 30 km due to the possibility of clipping. Short period data (EH channels) are only considered for the HS event due to its smaller magnitude. The seismograms at nearly all stations within 30 km from the TRIF event were clipped on the broadband and short period channels. For smaller events, there is less likelihood of clipping and thus HH or EH channels may be desirable at closer distances. Stacked EGF Method Choosing EGFs for obtaining spectral ratios is not a simple task. A number of recent studies have examined how different EGFs can affect the derived source properties. Commonly, EGF events are 1-2 magnitude units below that of the target event and they are ideally located as close as possible to the target event to have similar propagation paths (e.g. Kane et al., 2013b). Several studies used cross correlation to select suitable events (e.g. Abercrombie, 2015), with the idea that a good EGF event should have similar character to the target earthquake over a particular bandwidth. For time-domain studies, the EGF event should have very similar focal mechanism to that of the target event because source/structure interactions can significantly affect the phase information (e.g. Helmberger, 1983). For spectral analyses, the spatial proximity and focal mechanism similarity of the EGF and target events are less important (e.g., Calderoni et al., 2015). As mentioned, small events can have significant directivity (e.g., Boatwright, 2007; Lengliné and Got, 2011; Kurzon et al., 2014; Calderoni et al., 2015) and other non-generic source properties such as anomalous P/S radiation and isotropic source terms (e.g. Castro and Ben-Zion, 2013; Stierle et al., 2014a, b; Ross et al., 2015; Boettcher et al., 2015; Yang and Ben-Zion, 2016). Using an EGF generated by such events could bias the estimates of source properties of the target event, as demonstrated by Calderoni et al. (2013, 2015) for directivity and stress drops. Stacking EGFs produced by numerous small events can average out non-generic source effects. The need to 95 average out such features implies that the EGFs should not be selected with cross correlation, since this involves using similar events that may aggravate the problem. We select EGF events with magnitudes 1-2 units smaller than the target earthquake and hypocentral separation of no more than 5 km. If fewer than 5 such EGF events exist, we expand the separation distance to 7 km. Kane et al. (2013b) showed that the effect of separation distance on the average corner frequency is essentially constant between 2-14 km. If fewer than 5 events are available at distance up to 7 km, we skip analyzing the target earthquake. The maximum number of EGF events used is a computational choice rather than a scientific one. We do not include requirements on the focal mechanisms since we are strictly working in the frequency domain, where this is less of an issue. We tested the validity of this by comparing the spectra for different EGFs with a range of focal mechanisms, and found that each spectrum had the same general shape with different absolute amplitudes (Figure 5.S1). For the target earthquakes analyzed in this study, we use (Table 5.S1) 11 EGFs for the BB event, 14 EGFs for the TRIF event, 27 EGFs for the HS event, 24 EGFs for the CP event, and 11 EGFs for the EL event. Figure 5.S2 illustrates the reduction of scatter in the source spectral ratios by comparing results for the BB events associated with one EGF (left panel) with those derived with a stack of 1-11 EGFs at the different stations (right panel). Next, we run the P- and S-wave picking algorithms of Ross and Ben-Zion (2014) and Ross et al. (2016) on all available seismograms. These methods combine polarization filters with STA/LTA and kurtosis detectors in tandem. The P-wave picks are statistically accurate to ~0.1 sec, while the S-wave picks are statistically accurate to 0.25 sec (Ross et al. 2016). If a P-pick is successful for a given station/event, but an S pick is not, a V P /V S ratio of 1.7 is assumed to get an estimate of the S-wave arrival time. The same is applied if an S pick is successful but a P pick is not. All picks are checked to see if they fall within 1.5 sec of a predicted arrival time from a 1D velocity model (Hadley and Kanamori, 1977), and are discarded otherwise. This is done not to validate pick quality, but rather to discard serious outliers. The seismograms for the target event are windowed for P- and S-waves beginning 0.15 sec before each respective arrival pick. For P-waves, we use only the vertical components, while for S-waves, we use only the horizontal components. Any station without a pick for the target event is discarded. The window length used for calculating spectra is then determined automatically for each event and phase type. This is critical for an automated method applied to a large data set. In 96 general, the window length should be tied to f c of the target event; however at this point f c is unknown, so our procedure has two stages. The first is designed to get an initial estimate of the corner frequency. In the second stage, this estimate is used to refine the window length and frequency ranges, so the signal to noise (SNR) bandwidth is optimized and the quality and number of spectra is maximized. To obtain the initial window length, the magnitude of the target event is converted to seismic moment M 0 or potency P 0 . If the event magnitude is given in the moment magnitude scale M w , the relation of Hanks and Kanamori (1979) is used. Otherwise, if the magnitude scale is the local magnitude M L , it is converted first to seismic potency using the scaling relation of Ben-Zion and Zhu (2002) and then multiplied by a rigidity of 30 GPa. The scalar moment is used, together with an assumed stress drop, σ, to determine an estimate of f c by combining the equations of Madariaga (1976) and Brune (1970), 3 1 0 7 16 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = M k f c σ β , (4) where k is 0.34 for P-waves and 0.26 for S-waves (Kaneko and Shearer, 2015), and β is the shear wave velocity at the hypocenter depth. We calculate a lower bound for f c from (4) using a stress drop of 0.1 MPa. Then a minimum allowed frequency is calculated by subtracting 0.2 log(f) units from the minimum f c to ensure a sufficient number of points are below the corner frequency. The window length for this particular phase and event is then set equal to the inverse of the minimum allowed frequency. As examples, for the TRIF and BB events, the initial P-wave windows are found to be 5.1 sec and 2.3 sec in duration, respectively. The corresponding S-wave windows are 6.6 sec and 3.0 sec in duration, respectively. As mentioned previously, these windows are adjusted later after estimates of the f c values are obtained. Stations are only included when the S-P time is larger than the window length. This requirement is imposed for both phases, rather than just P, to ensure that there is negligible overlap between them. A pre-event noise window is then obtained for the purposes of checking SNR values. This window is chosen to end 2.0 sec before the P-wave pick to ensure no overlap and it has the same duration as used for windowing the particular phase. (The automatic pickers are used to verify that the noise windows are not likely to contain small events). From the signal and noise windows, we calculate spectra with a multitaper algorithm (Thompson, 1982; Prieto et al., 2009). For S-wave 97 data, the average of the two horizontal spectra is used. Example S-wave spectra for the TRIF earthquake are shown in Figure 5.2 for the station indicated by a green triangle in Figure 5.1. We require the SNR from the lowest frequency up to 30 Hz to be at least 3. For events with magnitude below 3, a maximum frequency of 30 Hz may be too low and need adjusting. If the SNR criteria are not satisfied at a given station, the data from that station are omitted. Figure 5.2. Examples of the different S-wave spectra derived for the Trifurcation event at station BLA2 (green star, Fig. 5.1). (a) Target event velocity spectrum (solid line) and noise spectrum (dashed line). (b) EGF velocity spectra (black) and stacked EGF (red). (c) Spectral ratio (in potency units) after deconvolving the stacked EGF from the target spectrum. Notice the sigmoidal shape resulting from a finite corner frequency in the EGF. We then apply this process to each of the EGF events and all stations. As the EGF events are lower in magnitude than the target earthquakes, the SNR is generally much lower and typically not all EGF records for a given station satisfy all the criteria. For some stations, no EGF records survive so the station is skipped. Since the SNR for S-waves is generally higher than for P-waves, this typically results in more spectra available for S-waves. The EGF spectra are then normalized individually by the seismic potency calculated from the potency-magnitude scaling relation of Ben-Zion and Zhu (2002). This converts the scale of each spectrum to one of unit potency, so that 98 all spectra are on equal footing. Examples of normalized EGF spectra are shown in Figure 5.2b. We stack all of the EGFs to get a smoother EGF (red line, Fig. 5.2b) that is less susceptible to directivity or other source effects than any single EGF is a priori. The target spectrum is then divided by the stacked EGF, yielding a spectral ratio for the given station (Fig. 5.2c). By following the steps outlined previously, the units of the spectral ratios are naturally transformed into those of the scalar seismic potency (cm·km 2 ) used in the scaling relation of Ben-Zion and Zhu (2002). This is noteworthy considering that no knowledge of the instruments gain, absolute original amplitude scales at different stations or material properties is necessary. It is simply the result of converting the EGF spectra to unit potency and taking the ratio between target and EGF spectra. One could further transform the scale into moment units by assuming a nominal rigidity of 30 GPa. A comparison between the zero frequency asymptotes for the P and S source spectra provides a check on the consistency of the derived results. Preliminary model fitting to estimate the corner frequency To extract information from a source spectrum, a model like (3) is commonly fit to a spectral ratio using a grid search, non-linear least squares or another method (e.g. Abercrombie, 1995). Often the high-frequency falloff rate, n, is fixed at 2 based on the model of Brune (1970), to avoid trade-offs between n and fc. On the other hand, arbitrarily fixing n to be a constant value is not fully justified, since it depends on the possible existence of tensile component of faulting and may vary azimuthally due to directivity effects (Kaneko and Shearer, 2015; Kwiatek and Ben-Zion, 2016). At the first stage of the analysis, we have sets of spectral ratios for both P- and S-waves based on window lengths and frequency range that are generally not optimized. The top panels in Figures 5.3-5 show spectral ratios for both phases for the TRIF, BB and EL earthquakes, respectively. These spectral ratios can be improved with a more appropriate window length for each specific event. To prepare for the fitting process, spectral ratios are re-sampled in the log domain with a constant log(f) spacing equal to 0.05 units (Ide et al., 2003). This gives more weight to the lower frequencies, which otherwise are overwhelmed by the high frequencies when fitting the model. After the spectral ratios are resampled, they are stacked over all of the stations to get a single spectral ratio that is representative of an average station. The stacking is performed by calculating the median value at each respective frequency, rather than the mean, which is designed to be robust 99 against outliers. For processing large seismic data sets, a scheme that is robust against outliers is necessary. The stacked spectra are indicated by the solid red lines in Figures 5.3-5.5. We then fit (3) to the stack spectrum using a grid search over the frequency range for which the required SNR is satisfied. The lower bound for f c 1 is the inverse of the window length, while the upper bound for f c 1 is determined from (4) by using a maximum stress drop of 100 MPa (equivalent with the assumed rigidity to a strain drop of 3.3·10 -3 ). This frequency range is divided into 100 constant log(f) increments. The lower and upper bounds for f c 2 are set to be half of the upper bound of f c 1 and 30 Hz, respectively. This range is also divided into 100 constant log(f) increments. The range for n is [1.5, 3.0] with steps of 0.05 and the range for Ω 0 is [max(r), 1.25*max(r)]. A value of γ=2 is used in (3) because it was found to estimate the corner frequency more reliably than a model with γ=1. These search ranges for the four parameters are used to minimize the sum of squared residuals. The best fitting models for the 3 example events are plotted in Figures 5.3-5.5 as red dashed lines. Figure 5.3. Spectral ratios (black lines) in potency units of both P and S phases for the TRIF event at all surviving stations. The red solid lines indicate the median spectrum, while the dashed red lines indicate the best-fitting spectrum to the stack using Eq. (3). (a) P-wave spectral ratios (51 total) obtained from the first analysis stage. (b) S-wave spectral ratios (53 total) obtained from first analysis stage. (c) P-wave spectral ratios (53 total) obtained from the second analysis stage after refining the window length. (d) S-wave spectral ratios (55 total) obtained from second analysis stage after refining the window length. The best-fitting estimates of the corner frequencies are indicated by the short vertical green marks and are listed in Table 5.2. 100 Figure 5.4. Spectral ratios (black lines) for both phases for the BB event (all surviving stations). The red solid lines indicate the median spectrum, while the dashed red lines indicate the best- fitting spectrum to the stack. (a) P-wave spectral ratios (52 total) obtained from first stage. (b) S- wave spectral ratios (65 total) obtained from first stage. (c) P-wave spectral ratios (74 total) obtained from second stage after refining window length. (d) S-wave spectral ratios (80 total) obtained from second stage after refining window length. The best-fitting estimates of the corner frequencies are indicated by the short vertical green marks and are listed in Table 5.2. Estimating strain/stress drops and radiated energy Since the window lengths used in the previous section are designed only to provide a reasonable initial estimate of the corner frequency, they rarely are optimal. In particular, the employed window lengths tend to be too long (by design), and overly restrictive on which spectra are used. This has the effect of making the lowest observed frequency too small, and as a result fewer stations pass the SNR criteria. In the second stage of the spectral analysis discussed here, we adjust the window length and minimum frequency to values that are tied to the initial estimate of the corner frequency. Extensive testing on different events suggests that window lengths equal to 4/f c 1 provide good balance between having sufficient frequency points below f c 1 and having enough stations/spectra with sufficient SNR. This sets the window length to be roughly one-half log unit below the source duration. We have tested values of 3-6 in the numerator and did not find a dramatic difference in 101 the derived source property values, but 4 seems to be the most robust. For the TRIF event, the new durations are 4.2 sec for P-waves and 4.1 sec for S-waves. For the BB event, the new durations are 1 sec for both phases, but the minimum frequency used was 1.3 Hz for P-waves and 1.0 Hz for S- waves. For the EL event, they are 1.6 sec for P-waves and 1.3 sec for S-waves. Requiring the window length to be much longer than the source duration is particularly important for stations with directivity to have a robust spectrum calculated. As example, for the TRIF event, which has strong directivity, using a window length shorter than this caused the spectra at some stations to never flatten at the low frequencies. A longer window at the same stations, however, was flat. Figure 5.5. Spectral ratios (black lines) for both phases for the EL event (all surviving stations). The red solid lines indicate the median spectrum, while the dashed red lines indicate the best- fitting spectrum to the stack. (a) P-wave spectral ratios (62 total) obtained from first stage. (b) S- wave spectral ratios (69 total) obtained from first stage. (c) P-wave spectral ratios (71 total) obtained from second stage after refining window length. (d) S-wave spectral ratios (82 total) obtained from second stage after refining window length. The best-fitting estimates of the corner frequencies are indicated by the short vertical green marks and are listed in Table 5.2. Using these new window lengths for P- and S-waves, we recalculate spectral ratios at all stations following the procedure of section 2.3. For cases where the adjusted window length is less than 1.0 sec, we fix the duration at this value. The resulting spectral ratios are shown in Figures 5.3-5.5c,d. For the TRIF, BB and EL events, there are now roughly 20% more spectral ratios available for use. As in the previous section, a stacked spectrum is calculated from all available 102 stations (red lines, Figs. 5.3-5.5). We then proceed to fit (3) to the stack using the same fitting procedure. The best fitting corner frequency values for P- and S-waves are shown in Table 5.2 for all five events, along with the corner frequency ratios and the number of spectral ratios per phase that were used in the stacking process. The obtained corner frequencies are used in conjunction with (4) to derive stress drops for the events based on the P and S source spectra. (Replacing in (4) the moment and stress drop with potency and strain drop, or dividing the obtained stress drops by the assumed nominal rigidity of 30 GPa, gives source strain drops for the events.) Table 5.2. Derived source properties for the set of analyzed earthquakes. For the Hot Springs event, the potency value is calculated using M L and the scaling relation of Ben-Zion and Zhu (2002), while for the other events they are calculated by converting M W values to moments using the relation of Hanks and Kanamori (1979) and dividing by a rigidity of 30 GPa. Energy values are calculated with Eq. (6). Directivity indices, D, are calculated with Eq. (9) by measuring the statistical splitting of spectra in opposing directions. Four of the five earthquakes are found to have directivity (from analysis of both P and S phases). An additional source quantity of interest is the radiated seismic energy. The J integral (e.g. Snoke, 1987) is proportional to the radiated energy and defined as, ∫ ∞ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + Ω = 0 2 0 ) / ( 1 2 df f f f J n c π . (5) Following Prieto et al. (2004) we calculate the J integral from the best fitting model parameters rather than the data, which allows extrapolating the integral to high frequencies. The best-fitting value of f c is from fitting a Boatwright model (Eq. 2b) with γ=2. We found that while the Quantity Big Bear Trifurcation Hot Springs Cajon Pass Elsinore # P-waves 74 53 104 100 71 # S-waves 80 55 114 111 83 f c P (Hz) 5.55 1.02 7.28 2.05 3.22 f c S (Hz) 4.63 0.80 6.65 1.73 3.18 f c P / f c S 1.20 1.28 1.10 1.18 1.01 σ P (MPa) 49.92 3.50 4.31 9.96 11.63 σ S (MPa) 64.69 3.75 7.32 13.43 24.91 σ avg (MPa) 57.30 3.63 5.82 11.69 18.27 E tot (J) 9.76E+10 6.48E+10 4.14E+08 8.07E+10 4.34E+10 E P /E S 12.36 10.21 16.20 12.85 20.43 E tot /P 0 2.30 0.14 0.25 0.48 0.86 D P 1.29 2.87 2.27 2.79 2.12 D S 1.31 4.30 4.97 3.35 4.34 103 Boatwright model was better at estimating f c 1 , the value of n was often not robust because, as seen in Figures 5.3-5.5, the upper corner frequency leads to a narrow frequency band over which the high frequency decay can be fit. For the purposes of estimating the seismic energy only, where a reliable value of n is critical, we use a standard reference value of n=2 (Brune, 1970). For a double- couple source in a homogenous whole space, the radiated seismic energy is defined as (e.g. Boatwright and Fletcher, 1984; Prieto et al., 2004), 2 0 0 5 2 4 Ω = J M c R E c C πρ , (6) where <R c 2 > is the average squared radiation pattern for the P- or S-waves (4/15 and 2/5, respectively). The total seismic energy is strongly affected by directivity (e.g. Venkataraman and Kanamori, 2004; Ma and Archuleta, 2006), and large deviations from a uniform station distribution. Moreover, as shown by Venkataraman and Kanamori (2004), correcting the seismic energy for directivity is difficult without additional knowledge such as a finite fault model, so the individual estimates of E P or E S may have appreciable errors. The ratio of these two quantities, however, is much less susceptible to these problems since the directivity effects are common to both phases. This results in similar energy flux distributions on the focal sphere. As such, we focus on E S /E P rather than E S or E P . The seismic energy ratios obtained for the five example events are shown in Table 5.2. The ratio of the total seismic energy to the scalar potency, often called the apparent stress (e.g., Aki, 1966; Brune, 1970), is also listed in Table 5.2. Quantifying rupture directivity Directivity alters the source duration at stations with different azimuths with respect to the rupture direction while preserving the area underneath the pulse. In the frequency domain this is equivalent to shifting the corner frequency up or down but keeping the low frequency asymptote the same. This process causes spectral splitting above the corner frequency, in which the spectra no longer overlap. We found that fitting (3) to spectral ratios at individual stations, and attempting to quantify directivity from azimuthal variations of f c 1 , is generally not a reliable procedure. Instead we focus on a technique that does not require fitting a model and uses statistics to directly quantify whether spectral ratios between different azimuthal directions are distinctly different. A similar approach was used by Calderoni et al. (2015), who analyzed aftershocks of the 2009 L'Aquila 104 earthquake in Italy and defined a rupture directivity index from an overlap percentage of spectral ratios between opposite directions along the fault. Figure 5.6. Analysis of spectral splitting for directivity in P-wave spectra of the TRIF event. (a) P-wave spectral ratios (black lines) along with median spectra in the NW (red line) and SE (blue line) directions. The grey vertical lines denote the range between the corner frequencies of the target and EGF events over which the spectra are analyzed for splitting. The median spectra are formed from stations within 45 degrees of the strike in the NW and SE directions. More than a factor of two difference in the corner frequencies is observed, indicating rupture directivity to the NW. (b) t-test values at each frequency measuring the statistical separation (Eqs. 7-8). A value of 2.0 (red dashed line) indicates statistical significance. The blue line shows the mean t-statistic (D) over the indicated frequency range (Eq. 9). Our method for estimating the degree of spectral ratio splitting first subsets the spectra that fall within 45 degrees of the strike of the event. For the five target earthquakes analyzed here, the strike is determined from the focal mechanism plane close the strike of the main nearby fault. For events with no prior knowledge of the fault trace, the analysis below can be done for both fault planes. Figures 5.6-5.11 show the subset of spectral ratios that satisfy this criterion for the TRIF, BB and EL events. The thicker solid lines indicate the median of these spectra in northwest (red) and southeast (blue). It is clear that while the spectra for the BB event overlap closely and have similar corner frequencies, the TRIF and EL events have different corner frequencies between the NW and SE directions. For the TRIF event there is more than a factor of two difference while for the EL event there is approximately a 50% difference. These features hold for both the P- and S-waves. The directivity effect is manifested visually most clearly by splitting in the spectral range above f c 2 , with lower amplitude in the forward rupture direction over this frequency range. This feature 105 is seen clearly for the TRIF and EL earthquakes. It is likely that the TRIF event had a unilateral rupture to the northwest, while the EL event had significant directivity to the SE. As expected, the high-frequency falloff rate for the TRIF and EL events with directivity are significantly different in the opposite directions. Figure 5.7. Same as Figure 5.6 for S-waves of the TRIF earthquake. An even larger value of D is obtained for S-waves than P-waves. Figure 5.8. Analysis of spectral splitting for directivity in P-wave spectra of the BB event. (a) P- wave spectral ratios with notations as in Figure 5.6. There is no observable spectral splitting suggesting a circular or weakly bilateral rupture (b) t-test values at each frequency measuring the statistical separation. The blue line indicating the mean t-statistic (D) over the used frequency range is below the value of 2.0 (red dashed line) used to measure statistical significance. 106 Figure 5.9. Same as Figure 5.8 for S-waves of the BB earthquake. No spectral splitting is observed, which is reflected in the D-statistic (blue line) below the dashed red line. Figure 5.10. Analysis of spectral splitting for directivity in P-wave spectra of the EL event. (a) P- wave spectral ratios with notations as in Figure 5.6. A modest spectral splitting in the opposite directions is observed suggesting rupture directivity to the southeast. (b) t-test values at each frequency measuring the statistical separation. A value of 2.0 (red dashed line) is used to measure statistical significance. The mean t-statistic (D) over the used frequency range (blue line) is above the dashed red line marking statistical significance. 107 Figure 5.11. Same as Figure 5.10 for S-waves of the EL earthquake. To quantify the state of directivity for both P- and S-waves, we try to identify if there is a statistically significant difference in the mean spectral amplitudes between the opposite along- strike directions. We calculate the t-statistic at each frequency value using, 2 1 2 1 X X s X X t − − = , (7) 2 2 2 1 2 1 2 1 n s n s s X X + = − , (8) where X is the mean spectral ratio value at a particular frequency. Eq. 7 indicates the number of deviations between the mean values on each side. Figures 5.6-5.11 show the t-test values at each frequency, over the range [0.5f c 1 , f c 2 ]. A value of ~2.0 (dashed line) is roughly the threshold that one expects for the means to be statistically different. Our directivity index, D, is calculated as the mean over this frequency range, ∑ = i i t N D 1 , (9) where t i is the t-statistic at the ith frequency. The value of D is indicated by the solid blue line. The TRIF earthquake has D > 2 for both P and S phases, while the BB earthquake rightly has D < 2. The D statistic only measures the statistical separation of spectra without indicating the direction 108 of rupture propagation. We identify the rupture direction from the plateau region above f c 2 by checking which direction has the lower average value. It is clear that the TRIF earthquake had a roughly unilateral rupture to the northwest. The BB event, on the other hand, does not exhibit any significant azimuthal variation of the corner frequency. It appears that this rupture was essentially circular in nature. The EL earthquake also has D values for both phases that are well above the threshold of 2.0, indicating there is a significant directivity component to the southeast. The D values for the remaining 2 events indicate that they are also likely to have directivity to the southeast, which we have confirmed visually from the spectral splitting. The values of the D index for the different events are given in Table 5.2. 5.3 Discussion We describe and implement a general-purpose automated methodology for estimating earthquake source properties from body wave spectra. An overview of the procedure is given by the flowchart in Figure 5.12. The primary goal of the method is to measure a set of earthquake source properties with no user involvement. Such a technique is necessary for systematic analyses of large seismic datasets. Our method works by starting with only the raw seismic waveform data and an earthquake catalog. For the southern California data used in this paper, waveforms of both P and S waves recorded by over 100 stations were analyzed for each target event. The 5 target earthquakes were analyzed with a total of 86 EGFs for calculating spectral source ratios. In deriving source properties for these five events, more than 30,000 spectra were examined automatically. The method employed the automated phase picking algorithms of Ross and Ben-Zion (2014) and Ross et al. (2016) rather than predicting arrivals from a model. The picking algorithms were found to work well and produced a sufficient number of picks for the desired source analyses. The algorithms produce roughly 1.3 P-wave picks for each S-wave pick, allowing for systematic analyses with both phases rather than just P. While the picking methods have errors, they are likely statistically smaller than those involved in predicted arrival times because of errors in locations and the velocity model. Furthermore, phase picks provide information on signal quality, since the picking methods operate on abrupt changes in the amplitude distribution of a time series. This is in contrast to predicting arrivals from a model and only then being able to check the SNR or other quality measures. 109 Figure 5.12. A flow chart summarizing the different components of the analysis method. We chose to use stacked EGFs produced by small events required only to be near the hypocenter of the target earthquake, over selecting highly-similar EGFs via cross-correlation (e.g. 110 Abercrombie, 2015), for several reasons. The first, and perhaps most important, is that small earthquakes can be affected by non-generic source effects just like large earthquakes. Using EGFs with such effects can bias the spectral ratios of the target events. Calderoni et al. (2015) paid considerable attention to this problem in the context of directivity, but other source mechanisms could also bias the EGF spectrum. One example is brittle rock damage in source volumes expected to produce (Ben-Zion and Ampuero, 2009) enhanced high-frequency P-wave radiation. Castro and Ben-Zion (2013) compared P/S spectral ratios between co-located event pairs and found that some had considerably more P-wave energy at high frequencies (> 10 Hz). Possibly related phenomena were observed by Kwiatek and Ben-Zion (2013), Boettcher et al. (2015), Ross et al. (2015) and Yang and Ben-Zion (2016). Similar source effects might exist in fluid-rich environments or bimaterial interfaces producing tensile components of faulting (e.g. Andrews and Ben-Zion, 1997; Ben-Zion and Huang, 2002; Minson et al., 2008; Nayak et al., 2014). Stacking non-similar EGFs reduces the influence of such source effects without the need for detailed examination of EGF spectra beforehand. Although there may be some sacrifice of quality in the stacked EGF compared to a single EGF with a large average cross-correlation coefficient, the risk of potential biases from non-generic EGF source features can have more adverse effects on the results. The minimum and maximum number of EGFs used in each stack is a choice between higher quality estimates per station and more stations overall. Not all stations will have the same number of EGF events satisfying the SNR criteria or have phase picks available. After extensive testing, we decided it was better to have no required minimum number of EGFs rather than omit certain stations altogether. The stacking of results from different stations helps to minimize the importance of spectral ratio quality, compared with fitting a model to each station independently. Additionally, having more stations included in the stack helps to ensure the focal sphere is better sampled. Our outlier detection scheme (median stack) appears to be sufficiently excluding erroneous spectral values. Although the analyzed five events are a very small data set, it is worthwhile discussing briefly the derived results (Table 5.2). The obtained strain/stress drops for the TRIF and HS events on the San Jacinto fault zone are close to median observed values (e.g. Abercrombie, 1995), somewhat higher for the CP earthquake on the San Jacinto and Elsinore event, and relatively high for the Big Bear earthquake consistent with other results in that area (Goebel et al., 2015). The derived P/S ratios of corner frequencies and radiated seismic energy are consistent overall with theoretical 111 expectations for shear-dominated ruptures (e.g., Molnar et al., 1973; Sato and Hirasawa, 1973; Madariaga, 1976) and observations (e.g., Venkataraman and Kanamori, 2004; Kwiatek and Ben- Zion, 2013). The directivity to the NW of the TRIF event, and the opposite directivity of the CP and HS earthquakes, are consistent with expectations for bimaterial ruptures (e.g., Weertman, 1980; Ben-Zion and Andrews, 1998) and the observed opposite polarity of the velocity contrast across the San Jacinto fault zone at these locations (Allam and Ben-Zion, 2012; Zigone et al., 2015), along with time-domain results of Kurzon et al. (2014). The directivity to the SE of the Elsinore event is also consistent with expectations for bimaterial ruptures and the observed velocity contrast across the fault (Zigone et al., 2015). It is important to test in future work whether these results characterize large data sets at the different locations. We have focused in this paper on developing a robust automated technique for estimating a variety of different source properties with no user involvement (Figure 5.12). While user inspection and adjustment of the frequency ranges, window lengths, EGFs, and many other parameters can improve the estimates of source properties of individual events, this reduces from the objectivity of results and can not be done on large data sets. The developed automated approach can be used to perform systematic analysis of large data sets. This improves the ability to address fundamental long-standing questions such as whether source properties vary between different regions, how much of the currently observed scatter is genuine and scaling of source properties with event size. These issues will be addressed in follow up studies. 112 6. Analysis of earthquake body wave spectra for potency and magnitude values: Implications for magnitude scaling relations (Ross et al., 2016) 6.1 Introduction The size of an earthquake is quantified with different parameters. One basic measure of the earthquake size is the scalar seismic potency P 0 given by the integral of the slip over the failure area or more generally the integral of inelastic strain over the source volume (e.g., Ben-Menahem and Singh, 1981; Ben-Zion, 2003). A related more commonly used parameter is the scalar seismic moment, M 0 , given by the product of the potency and effective rigidity in the source volume. The seismic potency and moment are proportional to the low-frequency asymptote of the displacement source spectrum (e.g. Aki, 1966; Ben-Zion, 2003). In addition to quantifying the size of earthquakes, the seismic potency/moment is also a critical ingredient in deriving many other source properties including stress/strain drop and radiated seismic energy (e.g. Aki and Richards, 2002; Prieto, et al. 2004). Estimating P 0 for small earthquakes is difficult because the wavefield recorded by a seismometer is affected by additional factors beyond the source radiation. It requires knowledge of how propagation and site conditions have altered a given spectrum so that they can be corrected for (e.g. Edwards et al., 2010; Abercrombie, 1995). This, in turn, requires a robust estimate of an earthquake’s location. The final step of extracting the long-period asymptote generally requires fitting a parameterized model to the corrected spectra, which has its own uncertainties. In contrast, the local magnitude, M L , requires only an estimate of the source receiver distance and a measurement of the peak amplitude of a seismogram filtered with the instrument response of a Wood-Anderson seismograph (Richter, 1935). It is a far more stable measurement for small to moderate events as it only has two dependent variables. This likely contributed to its longevity; it is still regularly used by seismic networks, even at the present day. However, the magnitude scale is not tied to a physical source quantity, so is only useful for quantifying relative differences between events, and it saturates for large events. Various studies have tried to develop scaling relationships between M L and M 0 or P 0 . The best known relation was developed by Hanks and Kanamori (1979) for moderate to large earthquakes and led to the moment magnitude scale M w . Several later studies (e.g. Bakun, 1984; Hanks and Boore, 1984; Abercrombie, 1996; Ben-Zion and Zhu, 2002; Shearer et al., 2006; Edwards et al., 2010) noted the apparent deviation of the M L 113 and M w scales for magnitudes below ~3.5. These studies analyzed events having M L above about zero. Earthquakes are well-known to be exponentially distributed with respect to magnitude (Gutenberg and Richter, 1944). The slope of the logarithm of this distribution is commonly referred to as the b-value. The b-value quantifies the relative frequency of events of different sizes and has been used widely in many applications (e.g., Utsu, 1999; Ben-Zion, 2008). Abercrombie (1996) found that the b-value of data recorded in a deep borehole is constant down to M L = 0. Given that the M w and M L scales diverge for small events, it is important to compare the b-values and completeness magnitudes of earthquakes quantified by both the different scales. In this study, we calculate seismic potency/moment and M w values for >11,000 small to moderate earthquakes which occurred in the Southern California plate-boundary around the San Jacinto fault zone (SJFZ) in 2013. The phase picks, event detections, and spectral fitting process are entirely automated, providing an objective and systematic set of results. In agreement with previous studies, the M L and M w scales are shown to diverge for M < 3.5. The results indicate systematic differences in the b-values and completeness magnitudes of the respective scales for same data, making these quantities scale dependent. A scaling relationship is developed between M L and M w that is suitable for converting between magnitudes of events in Southern California. 6.2 Data We analyze more than 13,570 earthquakes that occurred in 2013 and listed in a recently produced seismicity catalog for the region around the SJFZ (Ross et al., 2016). The catalog was built from scratch using the raw waveform archives. In addition to locations, the catalog contains 223,938 P-wave arrival picks and 199,647 S-wave arrival picks. The P-wave picks are made using the ratio of a short term moving average to a long term moving average (e.g. Allen, 1978; 1982), while the S-wave picks are made using the algorithm of Ross and Ben-Zion (2014) updated by Ross et al. (2016). These picks and the hypoDD method (Waldhauser and Ellsworth, 2000) are used to obtain high quality locations and origin times. The resulting relocated events (totaling 11,175) are compared with the SCSN catalog (SCEDC, 2013) to find events in common (8,082). The local magnitudes are taken from the SCSN catalog. In the following analysis, each event is required to have at least 5 P-wave and 5 S-wave picks to ensure sufficient quality of the derived 114 potency values. The waveform data are comprised of HH, EH, and HN channels at 87 stations. The 2013 earthquakes satisfying these criteria are shown in Fig. 6.1. Figure 6.1. Map of the San Jacinto fault zone region and >13,000 earthquakes which occurred during 2013, detected and located by Ross et al. (2016). Blue triangles indicate stations used. The yellow star is the location of the example event in Figure 6.2. The bottom panel shows a seismicity profile along the strike of the San Jacinto fault zone. 6.3 Methods A commonly used model for the spectrum of a seismic source is parameterized by a corner frequency, f c , low-frequency asymptote, Ω 0 , and spectral falloff, n (e.g. Abercrombie, 1995), 𝐴 𝑓 = $ . ∙012 (56'7 ∗ ) %& '/' : + (1) 115 where A is the displacement amplitude spectrum, t * is the whole-path attenuation operator, and f is frequency. The seismic potency is proportional to Ω 0 . Estimating Ω 0 for small earthquakes is difficult because of many issues. For example, there is often complex attenuation along the path and especially near the surface (e.g. Kilb et al., 2012; Liu et al., 2015), and the spectra are susceptible to noise issues (Edwards et al., 2010) and random heterogeneities (Shearer et al., 2006). Furthermore, the usable bandwidth is always limited, which makes determining the appropriate frequency range non-trivial. For each available event and station, the raw waveform data are demeaned and corrected for instrument response. We then obtain signal and noise windows that are 1.25 sec in duration. The signal window starts 0.25 sec before the pick, while the noise window ends 2 seconds before the P-wave pick. The choice of window length for earthquakes with M < 4 does not have an appreciable effect on the results. A multitaper algorithm is used (Park et al., 1987; Thompson, 1982) to calculate amplitude spectra for each of the three components for both P- and S-waves. The spectra are re-sampled in the log domain with Δlog(f) = 0.05 (e.g. Ide et al., 2003) and the 3- component spectra are combined via, 2 2 2 Z E N M + + = , (2) where M represents the vector magnitude spectrum and N, E, Z are the north, east, and vertical components, respectively. From here on, references to spectra are in the form of (2). Signal-to- noise ratio (SNR) spectra are obtained for each station-event pair by squaring the signal and noise spectra. For subsequent analysis we use only spectra that have at least 75% of the SNR values above 5.0. If fewer than 5 spectra were left for either P- or S-waves, the event is skipped. For the remaining events, spectra are corrected for radiation pattern using an average value 0.55 for P- waves and 0.62 for S-waves (Boore and Boatwright, 1984). At this point, all remaining spectra are integrated to displacement and then corrected for propagation and site effects. These correction factors include geometric spreading, path-dependent attenuation (t * ), as well as a free-surface correction of 2.0. We calculate t * for each source-receiver combination by tracing rays using a 1D version of the Allam and Ben-Zion (2012) velocity model (Table 6.S1). For these calculations, we use Q values for P- and S-waves from Abercrombie (1997) for depths less than 2.5 km and Q = 1000 at depths larger than 2.5 km (Table 6.S1). 116 Figure 6.2. Demonstration of the spectral fitting process for event 15269481 (SCSNID; yellow star, Fig. 6.1). The upper panel shows stacked spectra for P- (red) and S- (black) waves. The dashed lines are the best-fitting source models. The frequency ranges over which the spectra are plotted correspond to the ranges in which the SNR criteria are satisfied (bottom panel, dashed line). The lower panel shows the stacked SNR spectra for P- and S-waves. The determined moment magnitude for this earthquake is 1.35, while the SCSN M L is 0.90. The analysis employed 10 P- waves and 10 S-waves. These steps result in a set of P- and S-wave displacement source spectra for each earthquake. Next, we stack each event’s spectra for each phase type separately (Figure 6.2). The stacking is done by calculating the median logarithmic spectral value at each frequency bin rather than the mean. This keeps the stack from being biased by outlier spectra, without the need to discard them. We also stack the individual SNR spectra together to get an estimate of the average SNR at each frequency. This stacking is done in the log(SNR) domain to keep the closest stations from dominating the result (Fig. 6.2). The stacked SNR spectrum is used to determine the usable 117 frequency range for each phase stack automatically. While this range varies between different stations for a given event, the stacked SNR spectrum was found to be an excellent indicator of where artifacts began to appear in the stack. The artifacts at the low frequencies tend to cause the displacement spectra to tail upward rather than become asymptotic. A lower frequency bound, f L , is chosen to be the frequency value at the left end of the spectrum at which SNR > 5.0 was first satisfied. An upper frequency bound, f H, is chosen in the same manner while searching from the right, but is limited to a maximum value of 40 Hz due to peaks in the response of broadband seismometers above this value. Examples of the SNR stack and frequency range are shown in Figure 6.2. This imposes a minimum event size with a corner frequency of P-wave less than 40 Hz, which translates roughly into a moment magnitude of ~0.7. Next, the source model (1) is fit separately to each stacked phase spectrum for a given event. This is performed using a grid search minimizing the sum of the squared logarithm of the residuals. Using the logarithm of the residuals keeps the low frequencies from dominating the sum of squares. The search range used for Ω 0 is [0.75p, 1.25p] where p is the peak of the stacked displacement spectrum. The search range for the corner frequency [f L , f H ] and the search range for n is [1.5, 3.0]. The seismic potency is proportional to Ω 0 (e.g. Aki and Richards, 2002; Ben-Zion, 2003), 𝑃 < =4𝜋𝛺 < 𝑣 A , (3a) 𝑃 < =4𝜋𝛺 < 𝑣 2 B /𝑣 A C , (3b) where v S and v P are based on a 1D version of the velocity model of Allam and Ben-Zion (2012), and Eqs 3a and 3b are used for the S and P spectra, respectively. Examples of best-fitting spectral results are given in Figure 6.2 for the event indicated by a yellow star in Figure 6.1. The final potency value calculated for a given event is derived from a weighted average of the values obtained for the P and S waves, with the weights inversely related to the number of phases used. The obtained seismic potency value is then converted to moment by assuming a nominal rigidity of 30 GPa, and used to calculate M w with the scaling relation of Hanks and Kanamori (1979). 6.4 Results The discussed spectral fitting algorithm was applied to 11,175 earthquakes from the Ross et al. (2016) catalog. The method was able to calculate potency values for 11,091 events. Of these events, 8,078 were also detected by the SCSN, which form the final subset we use for analysis. Figure 6.3 shows frequency-size statistics of the obtained results using for event size both M w 118 derived from our results (top panel) and M L values taken from the SCSN catalog (bottom panel). For both magnitude scales, the results generally follow the Gutenberg-Richter statistics, bM a N − = ) log( , (4) where a and b are coefficients to be estimated and M is magnitude. To estimate the b-values of the results in Figure 6.3, we manually estimate the magnitude of completeness and find values of 1.0 and 1.6 (vertical blue lines) for the M L and Mw scales, respectively. Using these values, we estimate the b-value for the data with the maximum likelihood method (Aki, 1965) and obtain values of 0.93 and 1.16 for the M L and Mw scales, respectively. The magnitude values in the M w scale are generally shifted to the right (increased), but the shift is not constant for all magnitudes leading to a b-value that is about 25% higher. In contrast to many studies that compare b-values between different samples of events, the change of b-value in Fig. 6.3 results solely from using the two different magnitude scales for the same events. Figure 6.3. Frequency-magnitude distributions (FMD) for M w and M L values of the same data. In total, 8,078 of the 11,091 events with moment magnitudes were also detected by the SCSN. Only these common events are used. Cumulative and non-cumulative FMD are shown as black and red lines, respectively. The bin spacing for both plots is 0.1 units. The M w results have a b-value of 1.16, while the M L results have a b-value of 0.93. The completeness magnitude for M w is estimated at 1.6, while for M L it is 1.0 (blue lines). Note that since the same events are in the top and bottom panels, the larger magnitudes produce a b-value increase rather than a decrease. 119 Figure 6.4 shows more clearly the relationship between the M L and M w scales for our dataset. There are few events with M L > 4.0 and few values with M L < 0, so we focus on this magnitude range. A scatterplot of log(P 0 ) against M L is presented in Figure 6.4a along with a best-fit line (red), log(P 0 ) = c + dM L, (5) The slope of the best fitting line through these data is 1.13 (Table 6.1). This is similar to the value obtained by Ben-Zion and Zhu (2002) and other studies analyzing small events, and significantly different from the 1.5 value in the scaling relation of Hanks and Kanamori (1979) characterizing moderate and large events. Figure 6.4b shows the same information in a different form by plotting (M w - M L ) on the y-axis. It is clear that the magnitude scales are only approximately equal at M 3.5, and that for the majority of the range shown the M w values are larger statistically than the M L values. The difference between magnitude scales appears to increase linearly with decreasing magnitude, and at M L 0 the average earthquake has an M w value that is 0.88 units larger. From this data set, we also calculate a best-fitting line (Table 6.1). The standard deviation of the residuals is 0.19 units. The best-fitting values in Table 6.1 can be used to convert local magnitudes into moment magnitudes (or potency) for earthquakes in Southern California. However due to regional differences in attenuation, these values are not likely applicable elsewhere. 6.5 Discussion We derived potency estimates for >11,000 earthquakes with 0 < M L < 4 around the SJFZ from spectral analysis of both P- and S-waves. The potency calculations were entirely automated, involving as many as 85 stations in estimating each value of P 0 , which leads to a more objective and robust set of results. While the estimates of potency may be improved by manual adjustment of frequency ranges, window lengths, and other various parameters, this is labor intensive and difficult to perform for more than a few hundred events. The effect of attenuation on the low- frequency asymptotes is generally not appreciable. We performed the same analysis using Q = 1000 at all depths and found that the results were largely unchanged, with only slightly lower potency values (leading to less than 0.1M w difference on average). However, using more realistic estimates of Q in the upper 2.5 km from Abercrombie (1997) made the results between P- and S- 120 waves agree better. This demonstrates the need for including near-surface attenuation in the calculations. Figure 6.4. Comparison of derived potency, M w and M L values for more than 6400 earthquakes which occurred during 2013. The M L values were determined by SCSN, while M w values were determined with the developed procedure. The slope of the best-fitting line (red) in the upper panel is 1.13, significantly different from the 1.5 slope of Hanks and Kanamori (1979). The M w and M L values are roughly the same around M L 3.5, but below this magnitude M w progressively increases relative to M L . For M L = 0, the average M w is 0.88 units larger. The results confirm previous findings that the scaling relation between the potency (or moment) and magnitude of events with M L < 3.5 strongly deviates from the relation of Hanks and Kanamori (1979) for moderate and large events. The magnitude of completeness and b-value of events with size measured by M L differ from the values of the same events with size measured by M w . This is important to note because numerous studies in the literature have discussed variations of b-value with space and time with no attention to the scale used to measure the magnitudes. 121 Similarly, some studies have compared b-values across different magnitude ranges ignoring the different scales used. The derived scaling relations between log(P 0 ), M L and M w (Table 6.1) can be used to convert event size from one scale to the others. We note that this is strictly appropriate for the region around the SJFZ where the analyzed events occurred, since regional changes in attenuation can affect the M L values. However, the overall similar slopes for small events between log(P 0 ) or log(M 0 ) and M L found in other studies (e.g. Bakun, 1984; Abercrombie, 1996; Ben-Zion and Zhu, 2002; Shearer et al., 2006; Edwards et al., 2010) indicate that the regional variations are considerably smaller than changes in the scaling of events with M L less than about 3 and moderate to large events (e.g., M L > 5). Table 6.1. Magnitude scaling relation parameters. Quantity c d σ P 0 1.13 -4.06 0.29 M w -M L -0.246 0.884 0.19 The significant scatter in Figure 6.4 reflects the difficulties in calculating potency and moment values for small events. Measuring P 0 requires a high quality location, phase picks for both phases at many stations, precise estimates of the path and site contributions to a spectrum, and a robust technique for fitting a source model. Poor focal sphere coverage can further contribute to uncertainty. Here, we have focused on the statistical trend of these values for a large data set, rather than the values for any given event. The most meaningful conclusions are drawn from the trend of the curves shown in Figure 6.4, since the uncertainty in this trend is much less than the uncertainty in the potency of any single event. The standard deviation of the residuals about the curve is about 0.19 units. This should be kept in mind when converting the size of events from M L into P 0 , M 0 or M w , especially when done for individual events. However, converting an entire earthquake catalog with M L values into P 0 , M 0 or M w , and performing a statistical analysis on the results, should be more robust since much of this variability will average out. Local magnitude has been a valuable metric for quantifying the size of earthquakes for decades because of its ease of calculation and general stability. However, the local magnitudes saturate for large events and can underestimate the size of small events because of attenuation and recording issues. The moment magnitude scale was introduced several decades ago, but it is still used primarily only for moderate to large earthquakes. The scaling relation of Hanks and Kanamori (1979) does not characterize the relation between moments and magnitudes of small events. All 122 these issues have the potential of producing errors and biases in analyses of b-values or other parameters of earthquake catalogs. As another example, we note that in calculating spectral ratios for extracting earthquake source properties, a commonly employed rule of thumb is to use an empirical Green’s function from events (e.g. Hough and Dreger, 1995) that are 1-2 units smaller than the main earthquake to be studied. This rule of thumb was probably meant for moment magnitude, not local magnitude, and the results obtained in this study indicate using M L for this determination may lead to EGFs with corner frequencies that are considerably lower than expected. The source of the discrepancy between local magnitudes and moment magnitudes has been discussed by several studies in the past. Bakun (1984) and Hanks and Boore (1984) attributed the deviation to the Wood-Anderson filter that is applied to seismograms in determining local magnitude. The Wood-Anderson filter decreases in amplitude as the frequency increases and artificially lowers the amplitude of smaller events more than larger events due to the higher event corner frequency. Edwards et al. (2010) pointed out that a similar effect can be produced by seismic attenuation which strongly affects high frequencies of small events more than lower frequencies of larger events. See also Figures 6.11-6.12 of Kwiatek and Ben-Zion (2016) and related text. Ben- Zion and Zhu (2002) suggested that the different slopes of log(P 0 ) or log(M 0 ) for small and large events may reflect evolving stress heterogeneities leading to crack-like ruptures of moderate to large events and patchy fractal areas of small events. This can contribute a physical origin to the differences between the M L and M w scales of small events. In any case, the differences between the local and moment magnitudes should be kept in mind when analyzing earthquake data sensitive to the scale. 123 7. Spatio-temporal variations of double-couple aftershock mechanisms and possible volumetric earthquake strain (Ross and Ben-Zion, 2013) 7.1 Introduction Most studies on the physics and properties of crustal earthquakes focus on deviatoric/shear stress-strain components. This holds for routine derivations of earthquake source parameters, laboratory experiments, and simulations of earthquake and fault dynamics over both short and long time scales. While the changes of stress and strain during crustal earthquakes are predominantly deviatoric (e.g. Reid, 1910; Riedesel and Jordan, 1989; Dufumier and Rivera, 1997), even small dynamic variations of volumetric and normal stress-strain components can have fundamental implications for many aspects of earthquake and fault mechanics. For example, the evolution of fault resistance to motion with slip, and more generally the energy partitioning during faulting, can be altered dramatically by changes of normal stress and volumetric deformation. Transient changes of normal stress may be produced within earthquake rupture zones by mechanisms ranging from sliding on rough surfaces, ruptures on bimaterial interfaces and various fluid-thermal effects (e.g. Brune et al., 1993; Ben-Zion, 2001; Rice, 2006). Dynamic changes of elastic moduli in earthquake source volumes can generate damage-related radiation associated with a significant isotropic component (Ben-Zion and Ampuero, 2009). This is expected especially near rupture ends, geometrically-complex fault sections and regions without large preexisting faults. Recent observations of enhanced high-frequency P waves from aftershocks of the 2010 El Mayor- Cucapah earthquake may reflect damage-related isotropic radiation from earthquake source volumes (Castro and Ben-Zion, 2013). Individual earthquakes may be viewed as sensors of inelastic strain associated with slip patches of various sizes distributed throughout the seismogenic zone. The potency tensors P ij of earthquakes, computed from catalogs of focal mechanisms, can be used to describe 4D patterns of seismic strain fields. The degree of heterogeneity of a seismic strain field in a region can be quantified by calculating the distribution of rotations among sets of potency tensors. The stability of results can be increased using summed potency tensors instead of individual ones (e.g., Bailey et al., 2010). Typical focal mechanism inversions constrain the solutions to be deviatoric (e.g., Hardebeck and Shearer, 2002; Yang et al., 2012), so the isotropic components of individual source tensors are by default zero. Nevertheless, populations of derived deviatoric mechanisms can be 124 used to track collectively the 4D variations of the seismic strain field, including volumetric components, around earthquake rupture zones. In the present paper we attempt to find signatures of volumetric components of faulting by analyzing rotations of double-couple-constrained focal mechanisms. In order to have a high density of events in a geometrically-complex area, we examine aftershocks of the 1992 M w 7.3 Landers earthquake (Figure 7.1). In regions surrounding the north and south ends of the Landers mainshock rupture we find clear deviations in the distribution of rotations of early aftershock mechanism from the long term pattern. The deviations decay gradually within about a year from the time of the Landers mainshock to the common regional pattern. The results may be explained in terms of small isotropic source components in the early aftershocks near the rupture ends, where they are expected to exist, which were neglected in the derivation of the double-couple mechanisms. 7.2 Methodology and Results We examine spatio-temporal variations in the orientations of aftershock focal mechanisms, constrained to be double-couples, for evidence of isotropic source components. We anticipate the size of the isotropic source component of regular tectonic earthquakes to be generally small, but to increase in areas (e.g. near rupture ends) where the faulting process is expected to produce rock fracturing. To have statistically significant tests we need a sufficiently large number of events. In this work we analyze aftershocks of the 1992 M w 7.3 Landers earthquake that was associated with complex faulting and produced a high density of aftershocks in a narrow space-time window. We assume that early aftershocks respond to and largely reflect stress-strain changes produced by the mainshock, and that with increasing time later aftershocks return progressively to a response that reflects the regional background field. The source of each earthquake is quantified with the seismic potency tensor given by (Ben- Zion, 2003) ∫ = V T ij ij dV P ε , [1) where T ij ε is the transformational strain (e.g. Eshelby, 1957) and is the source volume. The potency tensor is a strain-based analog of the seismic moment tensor ij M (Backus and Mulcahy, € V 125 Figure 7.1. Aftershocks of the 1992 M 7.3 Landers earthquake (red star and focal mechanism) in the eastern California shear zone. The study area is divided into three approximately equal regions surrounding the two rupture ends and the central section. The epicenter locations (color-coded by depth), and associated double-couple-constrained mechanisms used in the analysis, are from the regional 1980-2010 focal mechanism catalog of Yang et al. [2012]. 1976; Ampuero and Dahlen, 2005) given by the product of the potency tensor with the tensor of elastic moduli at the source (Ben-Zion, 2008; Chapman and Leaney, 2012). We prefer to use potency tensors as they represent the seismic source without making assumptions on values of elastic moduli at the source volumes. The analysis employs the relocated Southern California focal 126 mechanism catalog of Yang et al. (2012), constrained to be double-couples and derived with the program HASH using first motion polarities and amplitude ratios (Hardebeck and Shearer, 2002). The seismic potency of each event is obtained from the magnitude using the empirical magnitude- potency scaling relation of Ben-Zion and Zhu (2002). The magnitude of completeness of earthquakes in the Landers area is estimated to be about M 2.2. There are 4139 recorded aftershocks around the Landers rupture with M ≥ 2.2 (small circles in Figure 7.1) from the time of the mainshock (06/28/1992) until the end of the employed catalog (12/31/2010). The number of M ≥ 2.2 aftershocks within 20 days of the mainshock is 1031. To examine aftershock focal mechanism orientations in regions expected to be associated with different amounts of rock fracturing, we divide the study area into three portions. Regions A and C are around the ends of the Landers rupture while region B is around the central section of the rupture zone (Figure 7.1). An earthquake source tensor can be defined by the orientations of its principal axes, P, T, and B, and the set of corresponding eigenvalues. To examine signatures of possible evolution or heterogeneity in the orientations of these axes, we need a measure of the difference between two orientation tensors A and B. This is provided by the angle necessary to rotate one set of principal axes into another, which is defined by (Kagan, 1991) ) arccos( 2 0 q = Ω , (2) where is the scalar component of a normalized quaternion k q j q i q q q 3 2 1 0 ' + + + = that describes the rotation. The angle of rotation is defined for ° ≤ Ω ≤ ° 120 0 , which means that no two sets of principal axes can be more than 120º apart. There are four different ways of rotating a source tensor to the same final orientation. Here we use the minimum of the four angles, which can be obtained by (Kuipers, 2002) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − = Ω 2 1 )) ( ), ( ), ( ), ( max( arccos 4 3 2 1 min R R R R Tr Tr Tr Tr , (3) where T B B B A A A ) , , )( , , ( 3 2 1 3 2 1 1 e e e e e e R = T B B B A A A ) , , )( , , ( 3 2 1 3 2 1 2 e e e e e e R − − = € q 0 127 T B B B A A A ) , , )( , , ( 3 2 1 3 2 1 3 e e e e e e R − − = T B B B A A A ) , , )( , , ( 3 2 1 3 2 1 4 e e e e e e R − − = , with and being the ith eigenvectors of tensors A and B, respectively. Below we denote for simplicity the minimum rotation angle as Ω. As the rotation angle does not take into account the eigenvalues of the source tensors, it solely represents their angular misalignment. In the following we use the mainshock potency tensor to define a reference orientation, and compute the set of minimum rotation angles between the reference and each of the aftershock tensors within given spatio-temporal bounds. The sampling distribution of these angles is an indicator of the heterogeneity of the aftershock focal mechanism orientations and traditional statistics can be used to describe the aspects of the sample. Using the Landers focal mechanism as a reference allows us to interpret the distribution with respect to the geometry of the mainshock. If the reference focal mechanism is changed, the axes of rotation also change leading to a different distribution of angles. We tested alternative reference tensors derived from weighted and normalized Kostrov summations of the aftershocks (Kostrov, 1974; Bailey et al., 2010) and find similar results. Figure 7.2a shows the histogram of rotation angles for the double-couple mechanisms of events within 20 days of the Landers mainshock in region A, along with the probability density function (dashed black line) for randomly oriented double couples (Kagan, 1992). Figure 7.2b presents corresponding results for all events in Region A with M ≥ 2.2, which provide a reference expectation for the background seismicity in the examined area. There are 412 events in the 20 day window and 1074 events that occurred in this region between the mainshock and the end of the catalog. The mean rotation angle, , is 65° for the 20 day window and 61° for the long term data. Additionally, the histogram for the 20-day event window shows a greater percentage of rotation angles above 70° than the histogram for the long term window. This suggests that early aftershocks in this region have focal mechanisms less aligned with the mainshock than the background seismicity in the region. To further quantify the extent of the differences between these two histograms, we use a two-sample Mann-Whitney U test (e.g. Hollander and Wolfe, 1973) to assess whether one of two independent samples has larger values than the other. When conducting this test, we first remove from the long term set the events of the 20-day set. The p-value for the U test is 10 − 5 and nearly any confidence level chosen will lead to a rejection of the null hypothesis € Ω 128 that the two samples are from the same population. Some of the aftershocks in Region A are in a zone with east-west orientation (Figure 7.1) that has a high angle to most of the other seismicity. To test that the scatter in the rotation angles for the short term data relative to the long term events is not strongly influenced by the seismicity in this east-west branch we remove those events and repeat the above statistical tests. The results indicate that the same behavior is observed with these events removed. Figure 7.2. Histograms of rotation angles of the double-couple-constrained mechanisms of aftershocks in the northern region A of Figure 7.1. The black dashed lines show the analytic probability density function for rotation angles between randomly oriented double couples (Kagan, 1992). (a) Results for events within 20 days of the Landers mainshock. (b) Results for all events in the 1980-2010 catalog. The results show considerable scatter from the orientation of the Landers mainshock mechanism and overall migration of angles from high to low values. The bars in the inset show the net changes of probability density for each bin between (a) and (b). In (a) there are approximately 9% more the rotation angles above 70° than in (b) and nearly all bins below 70° show a net increase of probability density. Using the same procedure for the central portion of the Landers area (Region B of Figure 7.1) we have 1384 events in the long term set and 349 events in the short term set (Figure 7.3). The histograms appear visually similar, suggesting that the focal mechanisms of early aftershocks here are produced in no different proportions than the background events. The mean rotation angle for events in Region B is 50.6° over the 20-day window and 49.2° over the long term. The U test for these data sets yields a p-value of 0.37, from which we conclude that there is no statistical evidence that one sample has larger values than the other. Performing similar analyses for events in the southern portion (Region C) of the Landers rupture, we find that the mean rotation angle in the 129 short term is 51.3° and 48.7° over the long term (Figure 7.4). The p-value for the U test in this region is 0.03, which suggests that the early aftershocks have statistically larger rotation angles than the background seismicity, as with the other edge region. Figure 7.3. Similar to Figure 7.2 for events in the central region B of Figure 7.1. The results exhibit considerably less scatter and temporal evolution than in Region A, with lack of statistically significant difference between the two populations. The inset in (b) shows that the net changes between (a) and (b) are relatively small compared to those in Figure 7.2. Figure 7.4. Similar to Figure 7.2 for events in the southern region C of Figure 7.1. The results exhibit intermediate scatter and temporal evolution compared to the results for Regions A and B. The inset in (b) shows the angles mostly decrease from (a) to (b) above 55°. It is useful to examine the magnitude range over which the results remain statistically significant. For Region A, the p-values are statistically significant at the 95% confidence level or greater for minimum magnitude M min in the range M min = 2.2-2.7. Above M min = 2.7, the number 130 of events in each sample is relatively small and, as a result, the statistical tests are no longer sensitive enough to resolve these effects. Region C is statistically significant at the 95% confidence level or greater for M min = 2.2-2.6, while Region B is found not to be statistically significant at the 95% confidence level for any minimum magnitude. The results of Figures 7.2-7.4 exhibit spatio-temporal variations in focal mechanism orientations. To obtain additional information on temporal evolution of these distributions in each region, we take first focal mechanisms within 5 days of the mainshock and compute the corresponding sets of rotation angles. This provides roughly 150-200 focal mechanisms for each region. Then we iteratively increase the time window duration while keeping the starting time (mainshock occurrence) fixed. For each window, we compute and continue this process until four years after the mainshock (6/28/1996). This procedure employs progressively larger number of events that provide increasingly more stable results that approach gradually (Figure 7.5) the long term values in each of the three regions. In each region, decreases with time and stabilizes after a certain amount of time at the long term values of Figures 7.2b, 7.3b and 7.4b. The long term mean values provide a measurement of how much the average focal mechanism deviates from the mainshock source configuration in a given space region. Since these values are obtained using time windows of four years, they should include some events that are not aftershocks. The mean rotation angles generally decrease monotonically in time, other than for a single data point in Region C at 10 days after the mainshock. This suggests that as time increases new events are produced with focal mechanisms that relate more closely to the Landers mainshock than the earlier aftershocks. Region A has the longest temporal decay and requires 1-2 years to stabilize. The central Region B has the fastest decay and takes only about six months to stabilize at the long term value. The mean rotation angle for Region B appears to decay with time like the other two areas, although this change is not detectable with the statistics used in the context of Figure 7.3. As with the northern edge, Region C has a longer decay than the central section and takes over one year to stabilize at the long term value. An alternative method of analyzing the rotation angles over time involves using a fixed number of events in the averaging window instead of an increasing amount of time. Using this type of sampling with various numbers of events per window produces the overall features of Figure 7.5, superposed with temporal oscillations that depend on the number of used events. A detailed € Ω Ω 131 analysis of the temporal oscillations requires careful analysis of various potential artifacts and is left for future work. Figure 7.5. The mean rotation angles (solid lines) and 1-standard deviation bounds (dashed lines) as a function of time for the double-couple-constrained mechanisms of aftershocks in the three regions of Figure 7.1. All rotation angles are computed with respect to the Landers mainshock focal mechanism. The horizontal axis denotes the length of the employed cumulative time window after the mainshock. The time scale for the mean rotation angles to stabilize in the two edges regions A and C is over one year, whereas for the central region B it is about six months. 7.3 Fault Plane Solution Uncertainties There is a difference of approximately 4° between the mean rotation angle measured 5 days after the mainshock and the mean angle for the long-term seismicity in Region A (Figure 7.2). The elevated mean rotation angles in the 20 day window relative to the background seismicity, and longer decay time to the background values near the main rupture ends (Figure 7.5) are statistically robust results. The set of faults associated with the Landers aftershocks have highly complex geometry and are prime settings for the production of rock damage during failure. Since the focal 132 mechanism catalog of Yang et al. (2012) was derived with the constraint that each event is a pure double-couple, an isotropic component in the true source mechanism will lead to an error in the fault plane solution. Because isotropic radiation adds a constant value to every point on the focal sphere, some P-wave radiation pattern values near the nodal planes may end up having their sign flipped relative to what is expected for a pure double-couple (Figure 7.6, inset). This can produce a large error in the double-couple--constrained solution. Below we test how significant such errors can be in the HASH algorithm by using synthetic first motion polarities. We represent the geometry of an earthquake source mechanism with a unit potency tensor (e.g. Ben-Zion, 2008). The far-field P-wave radiation pattern of a general can be written (e.g. Pujol and Herrmann, 1990; Aki and Richards, 2002) as ( ) j ij i P P R γ γ ˆ 2 = , (4) where are elements of the direction cosine vector. A general earthquake source tensor can be decomposed into isotropic and deviatoric components, with the deviatoric part having a double- couple term and a remainder term (e.g. CLVD). Given our focus on isotropic radiation, we assume that the deviatoric component is a pure double-couple (DC) and write as (Zhu and Ben-Zion, 2013) I ij DC ij ij P P P ˆ ˆ 1 ˆ 2 ζ ζ + − = . (5) Here ij I ij P δ 3 1 ˆ = with ij δ being the identity matrix, ζ is the fraction of the isotropic part of the source tensor taking values from −1 to 1 to account for both implosive and explosive cases, and DC ij P ˆ is the potency tensor for an arbitrarily oriented double couple. Using (4), the P-wave radiation pattern for the double-couple component of ij P ˆ can be written as ( ) P DC j DC ij i R P 2 2 1 ˆ 2 1 ζ γ γ ζ − = − , (6) where is the P-wave radiation pattern for an arbitrarily oriented double-couple (Aki and Richards, 2002) varying from −1 to 1. Similarly, the P-wave radiation pattern for the isotropic component of can be written as ij P ˆ ij P ˆ i γ ij P ˆ P DC R ij P ˆ 133 ζ γ δ ζγ γ γ ζ 3 2 3 2 ˆ 2 = = j ij i j I ij i P . (7) We define an effective radiation pattern for a source mechanism that includes a fraction ζ of isotropic component as ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + − ≡ ζ ζ 3 2 1 2 P DC P R R . (8) Along the double-couple nodal planes, where values are close to zero, a non-zero value of ζ can produce a shift of the actual nodal planes in the observed waveforms. As grid search algorithms that use polarities expect each quadrant in the focal sphere to have the same size, any flipped polarities introduced by an isotropic component will lead to error in the derived fault plane solutions and cause a rotation in the obtained eigenvectors. To obtain quantitative estimates, we write P DC R as a function of takeoff angle and azimuth as in Eq. (4.89) of Aki and Richards (2002), and use Monte Carlo simulations in two test cases to generate sets of radiation patterns for given ζ values. In the first case, we select 30 takeoff angles drawn from a uniform distribution over [0,90°] and 30 azimuths from a uniform distribution on [0,360°]. These are used to determine the radiation pattern for a specified strike, dip, and rake of the double-couple component; the particular sample size of 30 is an attempt to account for finite sampling of the focal sphere during any event. We then use the sign of the radiation pattern, along with the takeoff angles and azimuths, to determine the fault plane solution with HASH. This process is repeated 2000 times to get a set of fault plane solutions indicative of the range possible by the grid search algorithm. The assumed fault plane geometry of each event is used as a reference for computation of rotation angles with the solutions obtained from HASH. If this is done initially with 0 = ζ , the mean rotation angle, , represents the total uncertainty in the algorithm itself with 30 stations used. The value can be used to estimate the added uncertainty, ζ ΔΩ , to a fault plane solution for a given value of ζ: 0 Ω − Ω = ΔΩ ζ ζ , (9) where ζ Ω is the mean rotation angle for a given ζ. We investigated three basic focal mechanism cases to document how the HASH algorithm behaves. Figure 7.6 shows the obtained ζ ΔΩ values € Ω 0 € Ω 0 134 for each focal mechanism case with assumed ζ ranging from 0 to 0.15. For the Landers mainshock, the difference between the short (5 day) and long term mean rotation angles is approximately 4 degrees. If this difference is purely due to isotropic component of radiation, the early aftershocks have on the average ζ ~ 0.03 assuming strike slip mechanisms and ζ ~ 0.14-0.15 for the reverse and normal cases. As the distribution of focal mechanisms in the Landers area (Figures 7.2-7.4) is very diverse, the average ζ value for this region as estimated by the above technique would likely fall somewhere in the range of 0.05-0.10. The focal mechanisms retrieved in the synthetic tests with ζ = 0.03 have an increase of the mean inversion misfit by 3.14% compared to the focal mechanisms with ζ = 0. Examining the DC-constrained mechanisms of Yang et al. (2012) for the Landers aftershocks, we find that the mean inversion misfit for the 20 day window is 3% higher than the misfit for the long term window. This is in very good agreement with our synthetic results. The synthetic results leading to Figure 7.6 correspond to an ideal situation with a uniform coverage of the focal sphere. To estimate effects of isotropic radiation in a more typical observational setting, we consider a second test case with takeoff angles sampled from a distribution that is concentrated near the edge of the focal sphere. In this case we draw takeoff angles from a normal distribution with a mean of 90 degrees and a standard deviation of 30 degrees. This corresponds to a situation where the majority of stations are near the edge of the focal sphere (Figure 7.7, inset). As before, we perform sets of 2000 inversions for realizations associated with 30 stations each, and calculate the mean rotation angle for varying fractions of isotropic component. The inversion errors for different isotropic components are shown in Figure 7.7. As the station configuration is significantly different from that of the first case, the estimates of isotropic components for a given rotation error in Figure 7.7 differ from those of Figure 7.6, but both test cases produce similar trends. The reverse and normal focal mechanisms are more sensitive to the edge of the focal sphere than strike slip mechanisms. Based on the observed difference in the mean rotation angles of the short and long term results for region A (Figures 7.2a, 7.2b), we estimate average isotropic components of 0.11, 0.07, and 0.04 for events with strike-slip, normal, and reverse mechanisms, respectively. These numbers are in the same range as those obtained with uniform focal sphere sampling. 135 Figure 7.6. The uncertainty produced with the HASH algorithm by forcing a source with double- couple and isotropic components to be purely double-couple. The beach balls on top illustrate schematically how a small isotropic component can flip the polarities at stations near the nodal planes, leading to a rotation (error) in estimating a double-couple mechanism. The values shown correspond to an input range for the isotropic parameter ζ (equation 5) in the range 0-0.15. The thick lines show how the inversion errors (in degrees) increase with ζ for the three basic double- couple mechanisms. The results indicate that the mean error of 4º observed in the focal mechanism catalog (Figures 7.2-7.4) may reflect neglected isotropic components of 0.03-0.15. 7.4 Discussion A general faulting process, not limited to a pre-existing planar surface, produces brittle rock damage associated with changes of elastic moduli within and around the failure zone (e.g. Lockner et al., 1977; Hamiel et al., 2004; Stanchits et al., 2006). Seismic source tensors in regions sustaining coseismic changes of elastic moduli can have non-negligible isotropic components (Ben-Zion and Ampuero, 2009). This is expected especially in regions without pre-existing through-going fault traces, near the rupture ends, and around fault sections with large geometrical complexity (Castro and Ben-Zion, 2013). In section 2 and 3 we examine rotation angles of double-couple-constrained focal mechanisms of the 1992 Landers aftershock sequence (Yang et al., 2012). The analysis aims to clarify spatio-temporal variations of earthquake source tensors and possible existence of 136 volumetric strain changes in the source volumes. The examined mechanisms are derived with the HASH algorithm (Hardebeck and Shearer, 2002) as double-couples, so the isotropic components of individual source tensors are by default zero. Nevertheless, the population of the double-couple mechanisms can be used to assess likely variations of the coseismic strain fields, including volumetric components, at different locations. Figure 7.7. The uncertainty produced with the HASH algorithm for various levels of isotropic radiation using a ring-like station configuration with a sample realization illustrated by the focal sphere on top. The normal and reverse mechanisms are more sensitive than strike slip to the station configuration as more of their nodal regions are sampled with these source orientations. The estimated isotropic components for the three mechanisms are 0.04 for reverse, 0.07 for normal, and 0.11 for strike slip. The Landers mainshock propagated primarily in the north-northwest direction through a complex set of faults in the eastern California shear zone (Wald and Heaton, 1994). The curved geometry of the Landers rupture, which is especially prominent in the northern section, should lead to co-seismic generation of rock damage. The rupture and arrest processes near the north and south mainshock edges should include significant volumetric components that are especially prominent in the main propagation direction (e.g. Ben-Zion et al., 2012). We thus expect significant generation of rock damage around the northern end of the Landers rupture (Figure 7.1, 137 region A), relatively minor damage in the central section (Figure 7.1, region B) and intermediate damage near the southern rupture end (Figure 7.1, region C). These expectations are consistent with the results of Figures 7.2-7.5. The rotation angles of the double-couple mechanisms of the early aftershocks in region A are closer to the random distribution than in regions B and C; the difference between the mean rotation angle of the early aftershocks and the long term value is largest in the northern region A, intermediate in the southern region C and smallest in the central region B; the same order characterizes the time scales for the mean rotation angle to stabilize at the long term value. Figure 7.8. The mean number of polarities used for derivation of focal mechanisms in regions A- C of Figure 7.1 for each month following the Landers mainshock. Numerous studies have pointed out that faulting on non-planar geometries can produce deviatoric source tensors with CLVD components (e.g., Julian et al., 1998; Bailey et al., 2010). In the present work we focus on isotropic source terms and show that deriving source tensors of earthquakes with small isotropic components (equations 5-8) using an inversion algorithm that assumes purely double-couple mechanisms can lead to the observed variations in the rotation angles of the double-coupled-constrained mechanisms (Figures 7.6 and 7.7). We note that the analysis of section 3 intends to demonstrate effects associated with neglecting isotropic source 138 term. These results do not imply that the observed variations in rotation angles are produced exclusively or predominantly by isotropic source components. A more complete analysis should examine the likely contributions to scatter in double-couple-constrained mechanisms of neglecting isotropic and CLVD components (along with other possible mechanisms). Such analysis requires additional methods not used in the present paper and is deferred to a future work. Our results indicate that aftershock focal mechanisms become more like the mainshock as time after the parent event increases. To verify that the observed pattern is not produced by a changing stations configuration, we examine (Figure 7.8) the mean number of polarities used in the grid search per event for each month following the mainshock. The results show very similar evolution for regions A and B where we observe the largest and smallest temporal variations. Moreover, for the first seven months the station coverage does not increase significantly and for two of the regions it even decreases slightly. As a significant portion of the temporal decay present in Figure 7.5 occurs during this time window, we conclude that the observed phenomena are not likely the result of an increase in station coverage. It could be possible that spatio-temporal migration of the examined events led to progressive failures on different structures that contribute significantly to the observed changes in the rotation angle distributions. However, a plot of the aftershocks in the region color-coded by the occurrence day after the mainshock (supplementary Figure 7.S1) does not show any systematic migration of earthquake epicenters. Kagan (2000) investigated temporal changes of rotation angle distributions between pairs of earthquakes in the Harvard CMT catalog with lapse times between event pairs from 0 – 120 days. In contrast to our results, Kagan (2000) found that the mean rotation angle increases slowly with time and interpreted this as indication that earthquakes occurring closer in time have more similar source processes than earthquakes occurring further apart. The different observed patterns may reflect different processes occurring at the different involved scales. Our study was designed specifically to analyze low magnitude aftershocks, covering small spatio-temporal scales that may reflect local processes within and around the mainshock rupture zone. The global study of Kagan (2000) with only M > 5.5 events likely average over the phenomena observed in this study. Furthermore, if the results of Kagan (2000) were to hold at the same scales of our study, the observed increase of the mean rotation angle immediately after the Landers mainshock may provide stronger evidence of a neglected isotropic component in the source region. 139 When conducting the synthetic tests, we estimated the grid search uncertainty produced by neglecting ζ by shifting the total uncertainty for each value of ζ down with an amount equal to the mean rotation angle for ζ = 0. This value is a measurement of the accuracy of the HASH grid search algorithm given the number of used stations, assumed station sampling over the focal sphere, and assumed focal mechanism. Subtracting this value is equivalent to linearizing the uncertainty as a function of ζ to estimate the added uncertainty for each value of ζ. For ζ values between 0 and 0.2, a linear approximation seems reasonable for all three basic focal mechanisms (Figures 7.6 and 7.7); a simple linear regression to these data yields R 2 values larger than 0.90. While the linear trend may not continue beyond 0.20, we have only considered values below this threshold so the results provide adequate estimates of the uncertainty produced due to neglecting ζ in the grid search. The results of this study show that early aftershocks of the 1992 Landers earthquake have source mechanisms that are less aligned with the mainshock than the background seismicity occurring in the same region over multiple years. This could reflect, at least partly, a true phenomenon resulting from complex co-seismic and early post-seismic strain fields near the rupture ends. However, the variations in rotation angles of the double-couple-constrained mechanisms may also be produced, at least partially, as artifacts of the inversion process neglecting small isotropic source terms. These are expected to be generated by damage-related radiation (Ben- Zion and Ampuero, 2009) and other sources of volumetric deformation in fault sections with large geometrical complexity and near rupture ends (e.g., Sibson, 1986; Ben-Zion et al., 2012). A more direct estimate of isotropic earthquake source terms, and their contribution to scatter in double- couple-constrained mechanisms, can be done by combining full derivation of earthquake source tensors with analysis of the type done here. This will be done in a future work. 140 8. Isotropic source terms of San Jacinto fault zone earthquakes based on waveform inversions with a generalized CAP method (Ross et al., 2015) 8.1 Introduction Typical inversions of seismic waveforms for earthquake source properties assume pure deviatoric deformation with no volume change in the source region. The main reason for this is the dominance of shear deformation during earthquake failures associated with elastic rebound around the source volume (e.g., Reid, 1910). In addition, incorporation of small isotropic components in derivations of source tensors can lead to instabilities in the inversion process. A number of studies (Kawakatsu 1991, 1996; Hara et al., 1995; Dufumier and Rivera, 1997) performed detailed analysis of the correlation structure of different source components and inherent resolution of the isotropic components. The results from these studies suggest that in many cases, including the isotropic component in the inversion makes the problem intractable. Therefore, routine derivations of earthquake mechanisms are often done with inversions that are constrained to produce pure double-couple (or sometimes pure deviatoric) tensors (e.g. Hardebeck and Shearer, 2002; Yang et al., 2012). Inversions that include isotropic components are typically limited to analyses of deep events (Kawakatsu, 1991), explosions (Patton and Taylor, 2011; Xu et al., 2012), and earthquakes in geothermal, volcanic and other areas with elevated fluid content (Dreger et al., 2000; Minson et al., 2007; Ma et al., 2012). Standard inversions of seismic data for moment tensors (e.g., Backus and Mulcahy, 1976a,b; Tape and Tape, 2013) assume that the elastic moduli in the source volume are unchanged during the failure process. While this may hold for some cases such as deep seismic events, rock failure in the brittle lithosphere is accompanied by significant (transient) changes of elastic moduli (e.g. Lockner et al., 1977; Stanchits et al., 2006; Hamiel et al., 2009). The repeating occurrences of brittle failures in earthquake fault zones produce belts of damaged rocks (e.g. Dor et al., 2008; Wechsler et al., 2009; Mitchell et al., 2011) with seismic velocities that can be reduced by 30% or more (e.g., Lewis et al., 2005; Yang et al., 2011; Allam and Ben-Zion, 2012). The spatio-temporal variations of elastic moduli within and across earthquake rupture zones (and other regions sustaining brittle failure) motivates quantifying the sources of radiation with potency tensors that are based solely on the inelastic (transformational) strain in the source volumes (e.g. Ben-Zion 2003, 2008; Ampuero and Dahlen, 2005; Chapman and Leaney, 2012). 141 Ben-Zion and Ampuero (2009) derived a seismic representation theorem that includes, in addition to the classical inelastic strain (and corresponding stress glut), coseismic changes of elastic moduli in the source volume. The dynamic changes of elastic moduli are predicted to produce generically “damage-related-radiation” that can have an appreciable isotropic component. While the overall source may still be dominated by deviatoric deformation, a small isotropic component can lead to dynamic changes of normal stresses at the source region that can have significant effects on many aspects of earthquake physics (e.g. Brune et al., 1993; Ben-Zion, 2001). It is therefore important to resolve the possible existence of (small) isotropic earthquake source terms. Zhu and Ben-Zion (2013) provided decompositions of general (deviatoric plus isotropic) potency and moment source tensors that can be used efficiently in waveform inversions based on the Cut-and-Paste (CAP) method (Zhu and Helmberger, 1996). A number of other source tensor decompositions have been proposed and used by various authors (e.g., Julian et al., 1998; Chapman and Leaney, 2012; Tape and Tape, 2013; Vavryčuk, 2014). Different decompositions may be more suitable for different purposes and inversion schemes. In the present paper we use a generalized CAP (gCAP) algorithm based on the decomposition of Zhu and Ben-Zion (2013) to derive source tensor properties of seven M w > 4.2 earthquakes. The analyzed events occur in a geometrically complex region of the San Jacinto fault zone (SJFZ), and are therefore expected to produce damage-related-radiation with some isotropic source terms. Each solution is rigorously tested using a set of statistical techniques to quantify the uncertainty present in the inverse problem. The results indicate that small, but statistically robust, explosive isotropic components are likely present in at least six of the earthquakes examined. 8.2 Methodology and Results Seismic sources can be described in the point-source approximation by a symmetric second- order tensor that can be decomposed into constant-volume (deviatoric) and volumetric (isotropic) terms (e.g. Aki and Richards, 2002). This representation has many practical applications since the source tensor elements are linearly related to the generated ground motions. There are numerous algorithms for determining seismic source tensor properties with appropriate Green’s functions, using various decompositions and inversion schemes. The most common decomposition parameterizes the source as the sum of double-couples (DC), compensated linear vector dipoles 142 (CLVD), and isotropic (ISO) components (e.g. Julian et al., 1998). The choice of parameterization is non-unique, and is often done to provide bounds to individual source components that do not exist for the elements of the source tensor themselves. These bounds can be useful for both computational purposes as well as conceptual ones. In this study we use the source decomposition of Zhu and Ben-Zion (2013): ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + = ij ij ij D P P 2 0 1 3 2 ζ δ ζ , (1a) where ij P is the seismic potency tensor (e.g. Ben-Zion 2003), 0 P P < is the scalar potency, ij δ is the Kronecker delta, and ζ is the isotropic parameter. The latter is defined by 0 ) ( tr 3 2 P P = ζ , (1b) and it varies from -1 for implosive sources to 1 for explosive sources. ij D is a traceless tensor quantifying the deviatoric component of the source: CLVD ij DC ij ij D D D χ χ + − = 2 1 , (1c) with χ representing the strength of the CLVD source and having a range of [-1/2,1/2]. In (1c), DC ij D and CLVD ij D are normalized potency tensors for the double-couple and CVLD sources, respectively. The relative strengths of the individual components can be quantified (Zhu and Ben-Zion, 2013) by 1 = Λ + Λ + Λ DC CLVD ISO , (2a) where 2 ) sgn( ζ ζ = Λ ISO , (2b) ) 1 )( 1 ( 2 2 χ ζ − − = Λ DC , (2c) 2 2 ) 1 )( sgn( χ ζ χ − = Λ CLVD . (2d) The above decomposition is suitable for source tensor determination by grid search as the number of parameters is not too large. The CAP inversion technique breaks each waveform into a Pnl window and a surface wave window (Zhao and Helmberger, 1994; Zhu and Helmberger, 1996). One advantage of doing so is that the two windows can be filtered in different frequency bands to optimize the degree of fit during the inversion. For each parameter combination, a set of synthetic waveforms is generated and windowed into the same Pnl and surface wave groups. In this study the synthetics are formed using Green’s functions computed with the method of Zhu and Rivera 143 (2002) based on the 1D southern California velocity model of Hadley and Kanamori (1977) (referred to below as HK77). The data and synthetic waveforms are cross-correlated and the synthetics are shifted to maximize the cross-correlation coefficients to help overcome inaccuracies in the assumed velocity model and event location. Table 8.1. Ranges and step sizes for each of the parameters used in the inversions Parameter Range Step Size M w [1, 10] 0.1 Strike [0, 360] 1 Dip [0, 90] 1 Rake [-180, 180] 1 ζ [-1.0, 1.0] 0.05, 0.075, 0.1, 0.125, 0.15 χ [-0.5, 0.5] 0.05, 0.075, 0.1, 0.125, 0.15 Depth [3 km, 29 km] 2 km In addition to the traditional parameters of strike, dip, rake, M w , and hypocentral depth, the gCAP inversions performed here also include ζ and χ to account for possible isotropic and CLVD source terms. The parameters are listed along with their respective grid-search ranges and step sizes in Table 8.1. The strike, dip, and rake are obtained by a brute-force grid search through the whole parameter space volume. The grid search for ζ, χ, and M w starts from zero initial values for ζ and χ and the catalog magnitude for M w , and continues with given step sizes until a minimum is found. The process is repeated at all possible source depths to find the best set of parameters, and is iterated at several step sizes to reduce the possibility of locking on local minima (Table 8.1). At the end of each grid search, a quadratic interpolation is performed in the neighborhood of the grid point with minimum misfit to find the best parameter value and calculate the local curvature of the misfit surface for estimating the uncertainty of the value. As in Zhu and Helmberger (1996), the error function that is minimized is given by, ( ) ( ) ∑ = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = N i sw sw p i pnl pnl pnl p i s u r r s u w r r e sw pnl 1 2 0 2 0 , (3) where N is the number of stations, r i is the source-receiver distance, u and s denote the observed and synthetic seismograms, respectively, and parameters p pnl , p sw , w pnl , and r 0 have values of 1, 0.5, 2, and 100 km, respectively. Minimal pre-processing is performed on the seismograms used for the inversions. The Pnl window length is 35 sec long and phases are filtered between 0.05 – 0.3 Hz, while the surface 144 wave window is 70 sec long and phases are filtered from 0.02 – 0.1 Hz. The Green’s functions are computed up to 10 Hz, and are used to identify P and S-wave arrivals for windowing the seismograms into the Pnl and surface wave groups. Each component for every station is visually checked to confirm the quality of fit, and poorly fitting phases (e.g., those recorded on instruments located in basins) are often turned off. Additional inversion details including component weighting are documented in Zhu and Helmberger (1996). Figure 8.1. Map of the study area around the San Jacinto Fault Zone (SJFZ) in southern California (insert). Seven earthquakes with Mw > 4.2 in the complex trifurcation region (red box) of the SJFZ, denoted by a letter designation, are analyzed (Table 8.2). The shown focal mechanisms are double-couple-constrained from the catalog of Yang et al. (2012). The text near the mechanisms gives the date and magnitude of each event. The blue triangles denote the locations of stations used in the inversion of event G. We have chosen the region of study to be the geometrically-complex trifurcation area of the San Jacinto fault zone (Figure 8.1), where the primary fault structure splits into three subparallel strands. The trifurcation area also has one of the highest seismicity rates in California, and in particular produced 10 earthquakes with M w > 4 in the last 15 years. Nearly all of these earthquakes 145 were recorded by over a hundred stations at local and regional scales, and the large amount of data helps to constrain the inversion results. The CAP method has been found to work best for stations up to approximately 300 km from the source, and this cutoff distance is used in our data analysis. The closest stations used are 40 km away from the epicenter to be consistent with the point-source and far-field approximations. The focal mechanisms shown in Figure 8.1 are double-couple- constrained, derived by Yang et al. (2012) from first-motion polarities and amplitude ratios. Figure 8.2. Inversion results for the analyzed earthquakes. The blue focal mechanisms are unconstrained full source tensor solutions determined from gCAP waveform inversions. Each earthquake is found to have a small explosive isotropic component (see Table 8.2 for details). The yellow mechanisms are double-couple-constrained and shown with red “+” symbols that indicate the locations of stations used in the derivations. The blue dots are epicenters of background seismicity taken from the 1981-2011 catalog of Hauksson et al. (2012). Using the gCAP procedure and parameters described above, we perform source tensor inversions for the seven M w > 4 events shown in Figure 8.1 without imposing deviatoric constraints. The results for the set of inversions are visually displayed in Figure 8.2, with two focal mechanisms plotted for each earthquake. The unconstrained blue focal mechanisms include isotropic and CLVD components and can be compared to the yellow double-couple-constrained mechanisms obtained by the gCAP inversions. The red marks on each yellow focal mechanism 146 show the station locations on the focal sphere used in the inversion for that event. An example map-view of the station configuration used for the 2013-03-13 M 4.7 is given by the blue triangles in Figure 8.1. A sample of typical waveform fits is shown in Figure 8.3 for a set of stations at two distance ranges. In general, the fitness is very good for each of the inversions performed. Figure 8.3. Sample of data (black) and calculated (red) waveforms for 10 stations used for the 2013-03-13 earthquake (Event G). Phase components from left to right are vertical Pnl, radial Pnl, vertical surface waves, radial surface waves and tangential surface waves. The shown focal mechanism is unconstrained with station locations used in the derivations marked by red “+” symbols. The positive directions of motion for the different windows are: vertical (up), radial (outward), transverse (clockwise). Table 8.2 summarizes quantitative information on each event including retrieved isotropic and CLVD components, 95% confidence limits (explained in section 3), and results for statistical significance. The isotropic components generally have small ζ values (Eq. 1a) except for event D with ζ of 0.28. The corresponding percent of the total potency (Eqs. 2a-2b) taken up by the isotropic components ranges from 0.1-7.8%. We note that all the isotropic terms are explosive in type as expected for brittle damage production during faulting (Ben-Zion and Ampuero, 2009). The isotropic components remain explosive and similar in values when using a 1D version of the double-difference velocity model of Allam and Ben-Zion (2012) for the plate-boundary region around the SJFZ. The estimated CLVD components are also small, but frequently reverse sign. 147 Table 8.2. Results of the source tensor gCAP inversions for the analyzed earthquakes (Fig. 2). Confidence intervals for ζ (equations 1a-1b) are derived from bootstrapping,. The entries ζ (χ=0) correspond to best-fitting isotropic terms when the CLVD is constrained to be zero. All other inversion parameters are derived from full (unconstrained) tensor inversions. VR denotes variance reduction and NS denotes number of stations used for each inversion. Event Date Time Lat. Lon. Depth (km) M w Strike Dip Rake ζ ζ: 95% CI χ ζ (χ=0) % ISO VR NS A 10/31/2001 07:56:16.630 33.51 - 116.51 17 5 218 82 -47 0.10 [0.05, 0.14] 0.00 0.10 1.0 70 104 B 1/2/2002 12:11:28.680 33.37 - 116.43 15 4.2 223 85 -18 0.11 [0.4, 0.24] 0.04 0.11 1.2 76 71 C 6/12/2005 15:41:46.540 33.53 - 116.57 13 5.2 216 90 -28 0.06 [-0.01, 0.16] 0.01 0.08 0.4 80 102 D 5/1/2008 03:55:35.970 33.44 - 116.45 11 4.2 49 77 29 0.28 [0.18, 0.35] - 0.03 0.28 7.8 64 85 E 6/13/2010 03:08:57.090 33.38 - 116.42 11 4.9 234 77 -9 0.06 [0.0, 0.10] - 0.04 0.04 0.4 75 85 F 7/7/2010 23:57:23.040 33.42 - 116.48 13 5.4 49 90 5 0.03 [-0.02, 0.08] - 0.09 0.00 0.1 73 96 G 3/13/2013 16:56:06.040 33.50 - 116.46 13 4.7 216 81 -16 0.16 [0.07, 0.20] - 0.04 0.14 2.6 78 126 148 We show in the subsequent section that these CLVD components are generally not statistically significant. While a significant number of stations have been used in each inversion, this alone does not give confidence for the small derived isotropic components. In the next section we describe a number of statistical techniques used to quantify the uncertainty involved in the inverse problem, and show how the uncertainties are compared to the values summarized in Table 8.2. Table 8.3. Parameters of the HK77 velocity model. Layer depth (km) V P (km) V S (km) Q P Q S 5.5 5.5 3.18 1200 600 16.0 6.3 3.64 1200 600 32.0 6.7 3.87 1200 600 9999 7.8 4.50 1800 900 8.3 Tests for Statistical Robustness Bootstrap Analysis It is well known that using classical statistics to estimate the uncertainty for a single inversion is quite limited. However, the bootstrap theorem (Efron, 1979) allows one to resample the stations used in an inversion to quantify the variability present. We applied the bootstrap theorem to each of the events investigated in this study by randomly choosing stations with replacement from the set of stations used in the inversion. The number of stations resampled is equal to the number used for the results described in Figure 8.2; this leads to some being repeated multiple times while others are omitted. For example, the inversion of event G involved 126 stations, so for each bootstrap iteration 126 stations were randomly chosen from this original set of 126, with replacement, to create a resampled station configuration. A new inversion is performed using the resampled data set to obtain a new estimate for each of the inversion parameters. This process is repeated 1000 times, giving a statistical description of the variability possible. The inversion values for the strike, dip, rake, depth, and moment are well resolved and consistent with results reported by the double- couple-constrained focal mechanism catalogs of Yang et al. (2012). Further, the strike, dip, rake, and M w are almost identical for different values of ζ and χ. In contrast, the isotropic and CLVD components are small and variable, and thus are the focus of the bootstrap analysis. 149 Figure 8.4 presents a histogram of isotropic components resulting from the bootstrap procedure for each of the events (Fig. 8.4a), along with an empirical cumulative distribution function (ECDF) Figure 8.4. (a) Bootstrap distributions of the isotropic components of each event (different panels) obtained from resampling stations with replacement from the set used to perform the inversions shown in Fig. 8.2. For each event 1000 inversions are performed with a resampled station configuration. (b) Empirical cumulative distribution functions (ECDF) corresponding to each bootstrap distribution in (a). Each ECDF is used to obtain confidence limits and perform hypothesis tests for an explosive isotropic component. Six out of the seven earthquakes are found to have statistically significant values at the 95% confidence level, with the last one being significant at the 90% confidence level. For each event (Fig. 8.4b). The bootstrap histograms give a sense for how frequently the inversions are converging to a particular value when the station configuration is varied. If the standard 150 deviation of the histogram is small, it suggests uniformity among the results for all stations and adds greater confidence to the best-fitting parameter value. The ECDF gives the percentage of times that a resampled inversion produces an isotropic or CLVD component in a given range, which allows calculating confidence limits for the parameters. The vertical axis of the ECDF indicates percentile (in probability units) and the horizontal axis indicates the values of a given parameter (ζ or χ). Confidence intervals and hypothesis tests can be computed by reading off the value of the horizontal axis of the ECDF that corresponds to a specific percentile. We first perform a statistical hypothesis test that the isotropic component is explosive in type. Formally stated, we test a null hypothesis that the isotropic component is less than or equal to zero, and an alternative hypothesis that the isotropic component is greater than zero (explosive). This is based on the expectation from laboratory results and theory that brittle deformation (in the earthquake source volume) is associated with rock dilation (Lockner et al., 1977; Hamiel et al., 2009; Ben-Zion and Ampuero, 2009). If a 95% confidence level is desired for a 1-tailed hypothesis test, we use the ECDF to determine whether 95% of the resampled inversions produced an isotropic component greater than zero. If this is true, then the isotropic component is statistically significant with regards to the bootstrap procedure at the 95% confidence level and the null hypothesis is rejected. Such a procedure is non-parametric and avoids the stringent requirement that the data are drawn from a normally distributed population. Examination of the ECDF for the 7 analyzed events suggests that the isotropic components are statistically significant at the 95% confidence level, except for event F (Figs. 8.1-8.2) that has a statistically significant isotropic term at a 90% confidence level. One primary reason we can demonstrate significance for these events is that the amount of data used for each inversion is unusually large for regional inversions. We note, however, that it is the explosive nature of the isotropic components that is found to be robust rather than the best-fitting numerical values. Further, while statistical significance via the bootstrap theorem is encouraging, these tests implicitly assume that the Green’s functions are correct. Perturbations to the Green’s functions are addressed in section 3.4. Figure 8.5 shows corresponding results associated with applying bootstrap resampling to the CLVD components for each event. In contrast to the calculations for the isotropic components, the ECDF panels (Fig. 8.5b) indicate that 5 of the 7 events do not have statistically significant CLVD components at the 95% confidence level. It is clear from the histograms (Fig. 8.5a) that the sign 151 of the CLVD component regularly flips depending on the stations used, which is in direct contrast to the isotropic component. The bootstrap resampling calculations indicates that the non-zero CLVD for these earthquakes may be an artifact of the inversion process. Figure 8.5. Same as Fig. 8.4, but for CLVD components. Five of the seven events have CVLD components that are not statistically significant, which indicates they are likely artifacts of the inversion process. ISO to CLVD Partitioning A number of studies (Kawakatsu, 1991, 1996; Hara et al., 1995, 1996) examined the possibility of isotropic components in deep earthquakes primarily as evidence for polymorphic phase transitions. These studies have paid careful attention to parameter tradeoff in the inversion process, most notably between an isotropic component and vertically oriented CLVD sources. This tradeoff can occur when there is limited focal sphere coverage that is constrained primarily to the sphere’s 152 edges. This is because the radiation pattern for a vertically oriented CLVD is reminiscent of an explosion pattern in these areas, and in particular is constant in an azimuthal direction. However, the results from these studies of deep earthquakes at teleseismic distances may not necessarily translate at a regional scale with shallow crustal earthquakes. To examine the tradeoff present in the gCAP inversion process for the San Jacinto fault zone earthquakes, the inversions are first repeated while constraining the CLVD component to be zero. This allows us to understand how much the isotropic component for these specific results is influenced by the CLVD component. As indicated in Table 8.2, the isotropic components generally do not change in any significant way by the constraints χ = 0, and each remains explosive in type (except for event F, which is less statistically significant). To more rigorously study whether the isotropic component can map into the CLVD term, we conduct a series of synthetic inversions using a realistic station configuration. We assume the exact configuration used for the inversion of the 2013 M w 4.7 event G, and generate synthetic waveforms with pure deviatoric (DC + CLVD) sources. Three different types of DC orientations are used corresponding to the basic source types of strike-slip, reverse, and normal with 45 degree dips. For each, we generate synthetics with different values of χ and then perform an inversion while constraining the source to have only DC and isotropic terms. The range of χ values is from -0.5 to 0.5 in increments of 0.01. Non-zero values of ζ retrieved by these inversions are artifacts of CLVD components mapping into ζ. Figure 8.6 summarizes the results of the synthetic tests, with retrieved ζ values on the vertical axis and used χ values on the horizontal. For CLVD components smaller than 0.15, which is a generous upper limit on our derived χ values (Table 8.2), the isotropic components retrieved are no larger than 0.08. The retrieved ζ values corresponding to χ = 0.1 are around 0.05, which is smaller than all of the statistically significant isotropic components found for the SJFZ earthquakes (Table 8.2). Furthermore, this test assumes that all of the CLVD is mapping into the DC and isotropic components, whereas we allow χ to be a free parameter in our data inversions. These findings are in agreement with the fact that the isotropic components in the SJFZ earthquakes remain largely unchanged when constraining the CLVD component to be zero. From the CVLD bootstrap tests (Figure 8.5), it is clear that the CLVD term regularly takes on both positive and negative values for all of the events. For strike-slip focal mechanisms, however, only a positive CLVD component can map into an explosive isotropic component. This provides additional 153 evidence that CLVD components are unlikely to be responsible for derived isotropic components that are all explosive. Given the lack of significant trade off between the CLVD and ISO components in our case, and the fact that the CLVD terms of the analyzed events are not statistically significant (Figure 8.5), we perform additional studies of significance in the next sub- sections only on the isotropic components. Figure 8.6. Summary of synthetic tests for isotropic-CLVD tradeoff. Synthetic deviatoric waveforms are created for the station configuration used for event G (Fig. 8.2) for a given CLVD strength (horizontal axis). The synthetic waveforms are inverted with the constraint of zero CLVD term. The obtained isotropic values shown on the vertical axis describe parameter mapping and isotropic component error. For the range of CLVD values typical of tectonic earthquakes, the artificially retrieved isotropic components are generally smaller than isotropic components obtained in this study (Table 8.2). F-test Analysis One popular method of assessing the statistical significance of including an additional parameter in a regression model is the F-test. Using the F-test in source tensor inversions for isotropic components may be helpful (e.g. Dreger et al., 2000), but should not be thought of in the same way as a traditional regression model for several reasons. First, seismic data are known to not have normally distributed errors, which is a requirement of the F-test for true statistical significance. Second, the isotropic component has a very weak effect on the variance reduction. For example, if a true source has an isotropic component of ζ = 0.3, inverting the waveforms while 154 constraining ζ = 0 will only lower the variance reduction by a couple percent as compared with an unconstrained inversion. The F-test is directly based on large changes in the variance reduction, so it a poor choice from a theoretical standpoint for testing small isotropic components, regardless of how accurate the Green’s functions are. We conducted synthetic inversions with varying degrees of ζ, and found that even synthetic waveforms are not statistically significant with the F- test for ζ < 0.4. As a result, it is not informative to use the F-test for events with ζ < 0.1 as those analyzed in this study. Velocity Model Perturbations It is possible that the derived isotropic components result from inaccuracies in the velocity model that propagate through the inversion. This has been shown to be a valid explanation for some CLVD components (Kuge and Lay, 1994), in addition to their artificial generation by non- planar fault geometries, in derivations using the point source approximation (e.g. Julian et al., 1998). To examine whether reasonable velocity perturbations can produce uncertainties larger than the derived isotropic components, we perform two types of tests. The first test involves generating 500 different velocity models that are randomly perturbed from the original HK77 model. The perturbations are done to the thickness, depth, shear wave velocity, and the ratio of P over shear velocity of each layer. The perturbations to each parameter are drawn from a uniform distribution between [-5%, 5%] of the original value for that quantity. The value of 5% is generally the upper bound on uncertainty in the HK77 model over the analyzed frequency range compared with a number of recent 3D tomography studies in the study region (Allam and Ben-Zion, 2012; Allam et al., 2014; Zigone et al., 2014). For each event, we generate Green’s functions for each of the perturbed 500 velocity models and create 500 sets of synthetic double-couple seismograms. The DC orientation is taken from the focal mechanism catalog of Yang et al., (2012), and we invert each set of synthetics using the non-perturbed HK77 model. The synthetic seismograms are created for the exact station configuration of each earthquake as used previously (Fig. 8.2). If the events are really due to pure double-couple sources, and the true velocity model corresponds to one of the perturbed cases, the error introduced by using the non- perturbed HK77 model can be quantified. This is done by comparing the isotropic components obtained from these synthetic tests with a test distribution, such as the bootstrap results, in the form of a hypothesis test to see whether they are distinctly different. 155 Figure 8.7. Results of the first type of velocity model perturbation tests. The DC-constrained focal mechanism of each earthquake (Fig. 8.1) is used to create synthetic DC waveforms from 500 different perturbed velocity models at the same station configuration as the original inversion (Fig. 8.2). Each set of synthetics for a perturbed model is inverted using the HK77 velocity model while allowing isotropic and CLVD components. The obtained values of isotropic terms (blue histogram) are plotted alongside the values obtained from the bootstrap distribution (red) of Fig. 8.4. Comparisons via the Mann-Whitney U-test indicate that all 7 of the examined earthquakes have bootstrap distributions that are statistically larger than the distributions obtained from the perturbed velocity models at the 99% confidence level. The resulting distribution of ζ values (blue) is plotted in Figure 8.7 alongside the previously obtained bootstrap histogram (red) of each event. We use a Mann-Whitney U test (Hollander and Wolfe, 1973) to examine whether the values in one sample are statistically larger than the other. The tests indicate that for each event, the bootstrap distribution is statistically larger than the velocity perturbation distribution at the 99% confidence level. This suggests that the values obtained from the data inversions (Table 8.2) are unlikely to be the result of velocity model errors. 156 These results are visually seen by the clear separation (and relatively little overlap) between the red and blue histograms for each event in Figure 8.7. It is notable that the velocity model perturbations create both positive and negative isotropic components from the errors. If we approximate the distribution as being 50% explosive with the sign chosen at random for an artificial isotropic component produced by velocity model errors, the probability of all events having an explosive type as the result of velocity model errors is 0.5 7 which is 7.81E-03. An alternative method of testing the sensitivity of the isotropic component to velocity model errors involves performing a new inversion on the observed waveforms with each perturbed model. This is the opposite of the previous test, which inverted synthetic waveforms calculated for double- couple sources and the perturbed models. The resulting isotropic components in this second test can also represent the uncertainty limits due to errors in the velocity model on the values indicated in Table 8.2. Figure 8.8 shows the obtained histograms for each event (Fig. 8.8a) and the corresponding ECDF (Fig. 8.8b). Four of the seven events have isotropic components that are explosive (ζ > 0) for more than 95% of the perturbed models, while the remaining three have isotropic components that are explosive over 70% of the cases. In two events (D and G) the isotropic components are explosive over 99% of the time, so clearly several of these earthquakes have statistically significant explosive isotropic components with respect to the velocity model perturbations. We also performed both of the model perturbation tests using a different randomly drawn model at each station. For both test scenarios the statistical conclusions remain the same as discussed above. The conducted velocity model tests lead us jointly to conclude that the basic aspects of the derived results on isotropic source terms are unlikely due to errors in the HK77 model. 8.4 Discussion Resolving isotropic source terms of regular tectonic earthquakes is very challenging since the events are dominated by shear deformation. Trade-offs between parameters of the CLVD and isotropic components and uncertainties in the used velocity model can produce additional 157 Figure 8.8. Results of the second type of velocity model perturbation tests. Each earthquake is inverted 500 times using Green’s functions from a randomly perturbed velocity model to quantify the uncertainty in the isotropic component due to model choice. (a) Histograms of isotropic components for each event (different panels). (b) Empirical cumulative distribution functions corresponding to each sample in (a). The ECDF are used to perform hypothesis testing for explosive isotropic components. Four of the seven events have statistically significant isotropic components at the 95% level, one at the 90% level and two at about the 70% level. 158 difficulties. For these reasons routine inversions of seismic waveforms for source mechanisms constrain the inversions to be purely deviatoric (and often purely double-couple). While this is appropriate for studies interested in the overall kinematics of the source, brittle fracturing is expected to be associated generally with an isotropic source term (Ben-Zion and Ampuero, 2009). The existence of small isotropic radiation in earthquake rupture zones can change dramatically the physics of the local failure process, with significant implications for many topics ranging from the heat flow paradox to the crack vs. slip mode of earthquake rupture (Brune et al., 1993; Ben-Zion, 2001, and references therein). Clarifying the extent of isotropic earthquake source terms is also important for discrimination of small explosions from earthquakes, which rely largely on P/S ratios of recorded seismograms (e.g. Walter et al., 1995). A number of recent studies found evidence in support of isotropic components of regular tectonic earthquakes in areas likely associated with brittle fracturing. Ross and Ben-Zion (2013) examined double-couple constrained focal mechanism catalogs of the 1992 Landers, CA earthquake sequence. They found systematic spatio-temporal variations in rotations of aftershock focal mechanisms with strong variations near the edges of the rupture zone (especially around the NW end). They showed with synthetic calculations that the observed spatial heterogeneity of rotations can result from neglecting isotropic components with ζ = 0.03-0.15 when constraining the focal mechanisms to be double-couples. This estimated range of ζ is in agreement with the values obtained in this study (Table 8.2). Castro and Ben-Zion (2013) analyzed spectral ratios of seismograms generated by pairs of co-located aftershocks of the 2010 Mw 7.2 El Mayor-Cucapah earthquake. They observed amplification of P radiation of about 5-10 in the frequency range 1-10 Hz that is consistent with isotropic source components produced by rock damage. Kwiatek and Ben-Zion (2013) derived ratios of S-to-P radiated energy (E S /E P ) of >500 well-recorded small earthquakes in the Mponeng deep gold mine, South Africa. Many E S /E P ratios were relatively low (median value <5), suggesting tensile source components and enhanced P radiation for many events consistent with rock damage in the source volumes. The examined frequency ranges by Castro and Ben-Zion (2013) and Kwiatek and Ben-Zion (2013) are well beyond the 0.3 Hz maximum used in this work, and may be where the brittle damage is manifest more strongly. This is because the networks of microfractures created during the damage process may be significantly smaller than the overall size scales of the shearing sources. 159 Our results indicate small explosive isotropic source terms of crustal earthquakes in a complex region of the San Jacinto fault zone in southern California. The analysis is based on waveform inversions of data recorded by many 3-component stations at distances no more than 300 km from the sources. Previous studies calling attention to parameter correlation between isotropic and CLVD sources (Kawakatsu, 1991, 1996; Hara et al., 1995, 1996) primarily focused on earthquakes with focal depths of hundreds of kilometers, using recordings obtained at teleseismic distances. These studies typically used centroid moment tensor inversions with a greater number of parameters than used in the present work, and considered data in a considerably lower frequency range (20-1000s) than the ones used in our work (3-20s for body waves and 10-50s for surface waves). Also, in several of the aforementioned studies the retrieved isotropic components regularly flip sign between implosive and explosive sources, which is typical for inversion artifacts, in contrast to our results where all isotropic components are found to be explosive (Table 8.2). Dufumier and Rivera (1997) extensively analyzed the capability of extracting information related to the isotropic component from seismograms in terms of resolution and correlation matrices, with emphasis on inversions at teleseismic distance scales in different frequency bands than the ones used here. They discussed the instability of the isotropic component under certain conditions, but such instabilities are not observed in this work. Also, the joint inversion performed here weights the Pnl and surface wave groups differently to keep the larger amplitude arrivals from dominating the inversion (Zhu and Helmberger, 1996), an issue raised by Pearce and Rogers (1989) and Dufumier and Rivera (1997). The earthquakes analyzed in this work, as well as in the studies of Ross and Ben-Zion (2013), Castro and Ben-Zion (2013), and Kwiatek and Ben-Zion (2013), did not occur on geometrically- simple faults that may fail by frictional sliding alone. Some rock fracturing, and hence damage- related-radiation, is expected to accompany the earthquake source process in such locations. However, we can not rule out that additional mechanisms such as seismic anisotropy, faulting on non-planar surfaces, and fluid effects contribute to the observed isotropic radiation (e.g., Julian et al., 1998). Some of these mechanisms are likely to coexist; for example, faulting on non-planar faults and near rupture ends should lead to rock damage (e.g., Ben-Zion 2008, section 6) and brittle rock damage can produce seismic anisotropy in the source region (e.g. Hamiel et al., 2009). The damage-related-radiation is expected to be relatively large in area without well-defined preexisting faults and relatively small on geometrically-simple fault sections. We also expect shallower events 160 to produce more damage-related-radiation than deeper ones. These expectations should be tested in future studies involving larger areas and larger number of events. To check whether the isotropic source terms derived for the seven analyzed events (Figs. 8.1- 8.2) could be considered robust, we conducted bootstrap analysis, examined isotropic to CLVD partitioning, and explored the effects of velocity model perturbations on the retrieved isotropic components. Collectively, these tests suggest that the isotropic components derived for most or all of the events are statistically robust. Each of the examined earthquakes is found to have a small explosive isotropic component that may reflect damage-related-radiation stemming from brittle cracking in the source volumes. To increase the confidence that regular tectonic earthquakes (in geometrically complex regions) have isotropic source terms that may reflect damage-related- radiation, it would be useful to incorporate in the analysis 3D Green’s functions, include higher frequencies in the inversion process and examine data of additional events. Work on these topics will be done in future studies. 161 Discussion I have introduced a suite of algorithms for systematic processing of large seismic waveform data sets. The techniques include methods for earthquake detection and phase picking, identification of fault zone generated phases, and estimation of earthquake source quantities. When used together, the techniques are capable of constructing not just seismicity catalogs from scratch, but entire databases of seismological measurements that are suitable for large-scale downstream analyses. These methods have been applied to a number of other data sets not previously mentioned. They have been used on the TAIGER arrays in Taiwan (Wu et al., 2016), the SAHKE array in New Zealand (e.g. Henrys et al., 2013), and the MESONET array in Tokyo, Japan (e.g. Kasahara et al., 2010). The earthquakes and phase picks obtained by Wu et al. (2016) were used in a crustal shear wave splitting study (Okaya et al., 2016), where it was found that the shear wave fast directions correlated with tectonic terranes produced by plate convergence. The TAIGER and SAHKE data sets were collected as part of temporary deployments by researchers, and no formal catalog was previously produced for them. In the case of the TAIGER array, the 3-month data set was collected not for studying earthquakes but to record airguns as part of an onshore-offshore experiment. Wu et al. (2016) were able to detect nearly 9,000 earthquakes with the techniques developed herein, two-thirds of which were unknown beforehand. Thus, automatic techniques for systematic processing of large seismic waveform data sets are invaluable for fully utilizing the data sets that are collected around the world. Applying automatic tools to a data set requires a deep understanding of their limitations, however. There is often some amount of visual inspection necessary to ensure that the results are as expected for a new region or data set examined. This may be performed randomly on a small percentage of cases to ensure that the methods are working properly. All of the techniques developed in this paper have thresholds that are used at various stages to ensure that the quality of the results is sufficient. These thresholds are usually chosen to maximize the number of successes, while operating at a suitable success rate. In many cases the used threshold values can be data or region dependent, and may require significant time investment on the part of the user to optimize. The threshold values suggested in my studies were chosen to work well under a wide range of conditions, but may still require adjustment. Errors are also a key part of applying techniques like these, and will occur. As a result, automated techniques are most suitable for performing 162 systematic studies on large data sets, so that statistics can minimize the influence of errors or outliers. I believe that seismology is on the cusp of a fundamental transition to an era marked by data mining. This change will ultimately lead to detailed statistical analyses of individual fault zones for spatio-temporal variations in phases, sources, and structure. At the present day, the theory, data availability, instrumentation, and computing power are all more than sufficient to accomplish these goals. However, the automatic techniques for extracting information are often still not up to the task. It is my hope that the methods developed in this dissertation can bring the field a step closer to accomplishing this. 163 References Abercrombie, R. E., 1995, Earthquake source scaling relationships from− 1 to 5 ML using seismograms recorded at 2.5-km depth: Journal of Geophysical Research - Solid Earth, v. 100(B12), p. 24015-24036. Abercrombie, R. E., 2015, Investigating uncertainties in empirical Green's function analysis of earthquake source parameters: Journal of Geophysical Research - Solid Earth, v. 120(6), p. 4263-4277. Abercrombie, R.E., 1996, The magnitude-frequency distribution of earthquakes recorded with deep seismometers at Cajon Pass, southern California: Tectonophysics, v. 261, p. 1-7. Abercrombie, R.E., 1997, Near-surface attenuation and site effects from comparison of surface and deep borehole recordings: Bulletin of the Seismological Society of America, v. 87, p. 731-744. Abrahamson, N.A., and Silva, W., 2008, Summary of the Abrahamson and Silva NGA Ground- Motion Relations: Earthquake Spectra, v. 24, p. 67-97. Aki, K., 1965. Maximum likelihood estimate of b in the formula logN= a-bM and its confidence limits, Bulletin of the Earthquake Research Institute, Tokyo Univ., 43, 237-239. Aki, K., 1966, Generation and propagation of G waves from the Niigata Earthquake of June 16, 1964: Part 2. Estimation of earthquake moment, released energy, and stress-strain drop from the G wave spectrum: Bulletin of the Earthquake Research Institute, v. 44, p. 73-88. Aki, K., 1967, Scaling Law of Seismic Spectrum: Journal of Geophysical Research, v. 72 (4), p. 1217-1231. Aki, K., and Richards, P.G., 2002, Quantitative Seismology: University Science Books, 2nd edition, Sausalito, CA, USA. Allam, A. A., and Ben-Zion, Y., 2012, Seismic velocity structures in the Southern California plate-boundary environment from double-difference tomography: Geophysical Journal International., v. 190, p. 1181–1196, doi:10.1111/j.1365-246X.2012.05544.x. Allam, A. A., Ben-Zion, Y., Kurzon, I., and Vernon, F.L., 2014, Seismic velocity structure in the Hot Springs and Trifurcation Areas of the San Jacinto Fault Zone, California, from double-difference tomography: Geophysical Journal International, v. 198, p. 978–999, doi: 10.1093/gji/ggu176. 164 Allam, A.A., Ben-Zion, Y. and Peng, Z., 2014, Seismic Imaging of a Bimaterial Interface Along the Hayward Fault, CA, with Fault Zone Head Waves and Direct P Arrivals: Pure and Applied Geophysics, v. 171, doi:10.1007/s00024-014-0784-0. Allen, R., 1978, Automatic Earthquake Recognition and Timing from Single Traces: Bulletin of the Seismological Society of America, v. 68, p. 1521-1532. Allen, R., 1982, Automatic phase pickers: their present use and future prospects: Bulletin of the Seismological Society of America, v. 72, S225-S242. Ampuero, J.-P. and Y. Ben-Zion, 2008, Cracks, pulses and macroscopic asymmetry of dynamic rupture on a bimaterial interface with velocity-weakening friction: Geophysical Journal International., v. 173, p. 674–692, doi: 10.1111/j.1365-246X.2008.03736. Ampuero, J.P., and Dahlen, F.A., 2005, Ambiguity of the moment tensor: Bulletin of the Seismological Society of America, v. 95, p. 390–400. Andrews, D. J. and Y. Ben-Zion, 1997, Wrinkle-like Slip Pulse on a Fault Between Different Materials: Journal of Geophysical Research, v. 102, p. 553-571. Avallone, A., Rovelli, A., Di Giulio, G., Improta, L., Ben-Zion, Y., Milana, G. and Cara, F., 2014, Waveguide effects in very high rate GPS record of the 6 April 2009, Mw 6.1 L'Aquila, central Italy earthquake: Journal of Geophysical Research: Solid Earth, v. 119, p. 490–501, doi:10.1002/2013JB010475. Backus, G. and Mulcahy, M., 1976b, Moment Tensors and other Phenomenological Descriptions of Seismic Sources—II. Discontinuous Displacements: Geophysical Journal of the Royal Astronomical Society, v. 47(2), p. 301-329. Backus, G. and. Mulcahy, M., 1976a, Moment Tensors and other Phenomenological Descriptions of Seismic Sources—I. Continuous Displacements: Geophysical Journal of the Royal Astronomical Society, v. 46(2), p. 341-361. Backus, G., 1977, Interpreting the seismic glut moments of total degree two or less: Geophysical Journal of the Royal Astronomical Society, v. 51(1). Baer, M. and Kradolfer, U., 1987,. An Automatic Phase Picker for Local and Teleseismic Events: Bulletin of the Seismological Society of America, v. 77, p. 1437-1445. Bailey, I. W., Ben-Zion, Y., Becker, T.W., and Holschneider, M., 2010, Quantifying focal mechanism heterogeneity for fault zones in Southern and Central California: Geophysical Journal International., v. 183, p. 433–450, doi: 10.1111/j.1365-246X.2010.04745.x. 165 Baillard, C., Crawford, W.C., Ballu, V., Hibert, C. and Mangeney, A., 2014, An Automatic Kurtosis-Based P- and S-Phase Picker Designed for Local Seismic Networks: Bulletin of the Seismological Society of America, v. 104, doi:10.1785/0120120347 Bakun, W.H. and Lindh, A.G., 1985, The Parkfield, California earthquake prediction experiment: Science, v. 229, 619-624. Bakun, W.H., 1984, Seismic moments, local magnitudes, and coda-duration magnitudes for earthquakes in central California: Bulletin of the Seismological Society of America, v. 74, p. 439-458. Ben-Menahem, A., 1961, Radiation of seismic surface-waves from finite moving sources: Bulletin of the Seismological Society of America, v. 51(3), p. 401-435. Ben-Menahem, A., and Singh, S.J., 2012, Seismic waves and sources: Springer Science and Business Media, New York. Ben-Zion Y. and Huang, Y., 2002, Dynamic Rupture on an Interface Between a Compliant Fault Zone Layer and a Stiffer Surrounding Solid: Journal of Geophysical Research, v. 107 (B2). Ben-Zion Y. and Zhu, L., 2002, Potency-magnitude Scaling Relations for Southern California Earthquakes with 1.0 < ML < 7.0: Geophysical Journal International, v. 148, p. F1-F5. Ben-Zion, Y. and Aki, K., 1990, Seismic Radiation from an SH Line Source in a Laterally Heterogeneous Planar Fault Zone: Bulletin of the Seismological Society of America, v. 80, p. 971-994. Ben-Zion, Y. and Ampuero, J., 2009, Seismic radiation from regions sustaining material damage: Geophysical Journal International, v. 178(3), p. 1351-1356. Ben-Zion, Y. and Andrews, D., 1998, Properties and implications of dynamic rupture along a material interface: Bulletin of the Seismological Society of America, v. 88(4), p. 1085- 1094. Ben-Zion, Y. and Lyakhovsky, V., 2006, Analysis of aftershocks in a lithospheric model with seismogenic zone governed by damage rheology: Geophysical Journal International, v. 165(1), p. 197-210. Ben-Zion, Y. and Malin, P., 1991, San-Andreas Fault Zone Head Waves Near Parkfield, California: Science, v. 251, p. 1592-1594. 166 Ben-Zion, Y., 1989, The response of two joined quarter spaces to SH line sources located at the material discontinuity interface: Geophysical Journal International, v. 98, p. 213-222. Ben-Zion, Y., 1990, The response of two half spaces to point dislocations at the material interface: Geophysical Journal International, v. 101, p. 507-528. Ben-Zion, Y., 1998, Properties of seismic fault zone waves and their utility for imaging low- velocity structures: Journal of Geophysical Research-Solid Earth, v. 103, p. 12567- 12585. Ben-Zion, Y., 2001, Dynamic Rupture in Recent Models of Earthquake Faults: J. Mech. Phys. Solids, v. 49, p. 2209-2244. Ben-Zion, Y., 2003, Appendix 2: Key Formulas in Earthquake Seismology: in International Handbook of Earthquake and Engineering Seismology, eds. W. HK Lee, H. Kanamori, P.C. Jennings, and C. Kisslinger, Part B, 1857-1875, Academic Press. Ben-Zion, Y., 2008, Collective Behavior of Earthquakes and Faults: Continuum-Discrete Transitions, Progressive Evolutionary Changes and Different Dynamic Regimes: Reviews of Geophysics, v. 46, RG4006, doi:10.1029/2008RG000260. Ben-Zion, Y., Katz, S., and Leary, P., 1992, Joint Inversion of Fault Zone Head Waves and Direct-P Arrivals for Crustal Structure Near Major Faults: Journal of Geophysical Research-Solid Earth., v. 97, p. 1943-1951. Ben-Zion, Y., Peng, Z., Okaya, D., Seeber, L., Armbruster, J. G., Ozer, N., Michael, A. J., Baris, S. and Aktar, M., 2003, A shallow fault zone structure illuminated by trapped waves in the Karadere-Duzce branch of the North Anatolian Fault, western Turkey: Geophysical Journal International, v. 152, p. 699-717. Ben-Zion, Y., Rockwell, T., Shi, Z., and Xu, S., 2012, Reversed-polarity secondary deformation structures near fault stepovers: Journal of Applied Mechanics, v. 79, 031025, doi:10.1115/1.4006154. Ben-Zion, Y., Vernon, F.L., Ozakin, F.L., Zigone, D., Ross, Z.E., Meng, H., White, M., Reyes, J., Hollis, D., and Barklage, M., 2015, Basic data features and results from a spatially- dense seismic array on the San Jacinto fault zone: Geophysical Journal International., v. 202, p. 370–380, doi:10.1093/gji/ggv14. 167 Bennington, N.L., Thurber, C., Peng, Z., Zhang, H. and Zhao, P., 2013, Incorporating fault zone head wave and direct wave secondary arrival times into seismic tomography: application at Parkfield, California: Journal of Geophysical Research., v. 118, p. 1008-1014. Berckhemer, H., 1962, Die Ausdehnung der Bruchflӓche im Erdbebenherd und ihr Einfluss auf das seismische Wellenspektrum: Gerlands Beitr. Geophys., v. 71, p. 5–26. Beyreuther, M., Barsch, R., Krischer, L., Megies, T., Behr, Y., and Wassermann, J., 2010 ObsPy: A Python Toolbox for Seismology: Seismological Research Letters, v. 81(3), p. 530-533. Boatwright, J., 1980, A spectral theory for circular seismic sources; simple estimates of source dimension, dynamic stress drop, and radiated seismic energy: Bulletin of the Seismological Society of America, v. 70(1), p. 1-27. Boatwright, J., 2007, The persistence of directivity in small earthquakes: Bulletin of the Seismological Society of America, v. 97(6), p. 1850-1861. Boatwright, J., and Fletcher J.B., 1984, The partition of radiated energy between P and S waves: Bulletin of the Seismological Society of America, v. 74(2), p. 361-376. Boettcher, M. S., Kane, D.L., McGarr, A., Johnston, M., and Reches, Z., 2015, Moment Tensors and Other Source Parameters of Mining-Induced Earthquakes in TauTona Mine, South Africa: Bulletin of the Seismological Society of America, v. 105(3), p. 1576-1593. Boore, D. M., and Boatwright, J., 1984, Average body-wave radiation coefficients: Bulletin of the Seismological Society of America, v. 74(5), p. 1615-1621. Brietzke, G.B., Cochard, A., and Igel, H., 2009, Importance of biomaterial interfaces for earthquake dynamics and strong ground motion: Geophysical Journal International, v. 178(2), p. 921–938. Brune, J. N., 1970, Tectonic stress and the spectra of seismic shear waves from earthquakes: Journal of Geophysical Research, v. 75(26), p. 4997-5009. Brune, J., Brown, S., and Johnson, P., 1993, Rupture mechanism and interface separation in foam rubber models of earthquakes: a possible solution to the heat flow paradox and the paradox of large overthrusts: Tectonophysics, v. 218, p. 59-67. Bulut, F., Ben-Zion, Y., and Bohnhoff, M., 2012, Evidence for a bimaterial interface along the Mudurnu segment of the North Anatolian Fault Zone from polarization analysis of P waves: Earth and Planetary Science Letters, v. 327, p. 17-22. 168 Bulut, F., Bohnhoff, M., Ellsworth, W.L., Aktar, M., and Dresen, G., 2009, Microseismicity at the North Anatolian Fault in the Sea of Marmara offshore Istanbul, NW Turkey: Journal of Geophysical Research., v. 114, B09302, doi:10.1029/2008JB006244. Calderoni, G., Di Giovambattista, R., Vannoli, P., Pucillo, S., and Rovelli, A., 2012, Fault- trapped waves depict continuity of the fault system responsible for the 6 April 2009 MW 6.3 L'Aquila earthquake, central Italy: Earth and Planetary Science Letters, v. 323-324, p. 1--8, doi:10.1016/j.epsl.2012.01.003. Calderoni, G., Rovelli, A., and Singh, S.K., 2013, Stress drop and source scaling of the 2009 April L’Aquila earthquakes: Geophysical Journal International, v. 192(1), p. 260-274. Calderoni, G., Rovelli, A., Ben-Zion, Y., and Di Giovambattista, R., 2015, Along-strike rupture directivity of earthquakes of the 2009 L'Aquila, central Italy, seismic sequence: Geophysical Journal International, v. 203(1), p. 399-415. Castro, R. R. and Ben-Zion, Y., 2013, Potential Signatures of Damage-Related Radiation from Aftershocks of the 4 April 2010 (Mw 7.2) El Mayor-Cucapah Earthquake, Baja California, México: Bulletin of the Seismological Society of America, v. 103, p. 1130– 1140, doi: 10.1785/0120120163. Chapman, C., and Leaney, W., 2012, A new moment-tensor decomposition for seismic events in anisotropic media: Geophysical Journal International, v. 188(1), p. 343-370. Cichowicz, A., 1993, An Automatic S-Phase Picker: Bulletin of the Seismological Society of America, v. 83, p. 180-189. D'Agostino, R.B., and Stephens, M.A., 1986, Goodness-of-Fit Techniques, Vol. 68, CRC press. Dor, O., Yildirim, C., Rockwell, T.K., Ben-Zion, Y., Emre, O., Sisk, M., and Duman, T.Y., 2008, Geologic and geomorphologic asymmetry across the rupture zones of the 1943 and 1944 earthquakes on the North Anatolian Fault: possible signals for preferred earthquake propagation direction: Geophysical Journal International, v. 173, p. 483-504, doi: 10.1111/j.1365-246X.2008.03709.x. Dreger, D. S., Tkalcic, H., and Johnston, M., 2000, Dilatational Processes Accompanying Earthquakes in the Long Valley Caldera: Science, v. 288, p. 122-125. Dufumier, H. and Rivera, L., 1997, On the resolution of the isotropic component in moment tensor inversion: Geophysical Journal International, v. 131(3), p. 595-606. 169 Dziewonski, A., Chou, T., and Woodhouse, J., 1981, Determination of earthquake source parameters from waveform data for studies of global and regional seismicity: Journal of Geophysical Research, v. 86(B4). Eberhart-Phillips, D. and Michael, A., 1993, Three-Dimensional Velocity Structure, Seismicity, and Fault Structure in the Parkfield Region, Central California: Journal of Geophysical Research.-Solid Earth, v. 98, p. 15737-15758. Edwards, B., Allmann, B., Fäh, D., and Clinton J., 2010, Automatic computation of moment magnitudes for small earthquakes and the scaling of local to moment magnitude: Geophysical Journal International, v. 183(1), p. 407-420. Efron, B., 1979, Bootstrap methods: Another look at the jackknife: Ann. Statist., v. 7, p. 1–26. Ellsworth, W.L., and Malin, P, 2011, Deep rock damage in the San Andreas Fault revealed by P- and S-type fault-zone-guided waves: Geological Society of America, special publication, v. 359 (1), p. 39. Enescu, B., Hainzl, S., and Ben-Zion, Y., 2009, Correlations of seismicity patterns in Southern California with surface heat flow data: Bulletin of the Seismological Society of America, v. 99(6), p. 3114-3123. Eshelby, J., 1957, The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems: Proceedings of the Royal Society A, v. 241(1226), p. 376-396. Fletcher, J., Haar, L., Vernon, F., Baker, L., Hanks, T., and Brune, J., 1987, The digital array at Anza, California:Processing and initial interpretation of source parameters: Journal of Geophysical Research., v. 92, p. 369-382, 1987. Fohrmann, M., Igel, H., Jahnke, G., and Ben-Zion, Y., 2004, Guided waves from sources outside faults: an indication for shallow fault zone structure?: Pure and Applied Geophysics, v. 161, p. 2125-2137. Fukao, Y., Hori, S., and Ukawa, M., 1983, A seismological constraint on the depth of basalt– eclogite transition in a subducting oceanic crust: Science, v. 303, p. 413-415, doi:10.1038/303413a0 Gentili, S., and Michelini, A., 2006, Automatic picking of P and S phases using a neural tree: Journal of Seismology, v. 10, p. 39-63. Gilbert, F., 1971, Excitation of the normal modes of the Earth by earthquake sources: Geophysical Journal of the Royal Astronomical Society, v. 22(2), 223-226. 170 Goebel, T., Hauksson, E., Shearer, P., and Ampuero, J., 2015, Stress-drop heterogeneity within tectonically complex regions: a case study of San Gorgonio Pass, southern California: Geophysical Journal International, v. 202(1), p. 514-528. Gutenberg, B., and Richter, C.F., 1944, Frequency of earthquakes in California: Bulletin of the Seismological Society of America, v. 34, p. 185-188. Hadley, D., and Kanamori, H., 1977, Seismic structure of the Transverse Ranges, California: Geological Society of America Bulletin, v. 88, p. 1469-1478. Hamiel, Y, Lyakhovsky, V., Stanchits, S., Dresen, G., and Ben-Zion, Y., 2009, Brittle Deformation and Damage-Induced Seismic Wave Anisotropy in Rocks: Geophysical Journal International, v. 178, p. 901–909, doi: 10.1111/j.1365-246X.2009.04200.x. Hamiel, Y., Liu, Y., Lyakhovsky, V., Ben-Zion, Y., and Lockner, D., 2004, A Visco-Elastic Damage Model with Applications to Stable and Unstable fracturing: Geophysical Journal International, v. 159, p. 1155-1165, doi: 10.1111/j.1365-246X.2004.02452.x. Hampel, F. R., 1974, The Influence Curve and its Role in Robust Estimation: Journal of the American Statistical Association, v. 69, p. 383-393. Hanks, T. C. and Boore, D. M., 1984, Moment-magnitude relations in theory and practice: Journal of Geophysical Research.: Solid Earth, v. 89, p. 6229-6235. Hanks, T. C. and Kanamori, H., 1979, A Moment Magnitude Scale: Journal of Geophysical Research., v. 84(B5), p. 2348-2350. Hara, T., Kuge, K., and Kawakatsu, H., 1995, Determination of the isotropic component of the 1994 Bolivia Deep Earthquake: Geophysical Research Letters, v. 22(16), p. 2265-2268. Hara, T., Kuge, K., and Kawakatsu, H., 1996, Determination of the isotropic component of deep focus earthquakes by inversion of normal-mode data: Geophysical Journal International., v. 127(2), p. 515-528. Hardebeck, J., and Shearer, P., 2002, A new method for determining first-motion focal mechanisms: Bulletin of the Seismological Society of America, v. 92(6), p. 2264—2276. Hauksson, E., Yang, W., and Shearer, P.M., 2012, Waveform Relocated Earthquake Catalog for Southern California (1981 to June 2011): Bulletin of the Seismological Society of America, v. 102(5), p. 2239-2244. 171 Helmberger, D., 1983, Theory and application of synthetic seismograms: in Earthquakes: observation, theory and interpretation, proceedings of the International School of Physics ‘Enrico Fermi’, v. 37, p. 173-222. Henrys, S., Wech, A., Sutherland, R., Stern, T., Savage, M., Sato, H., Mochizuki, K., Iwasaki, T., Okaya, D., Seward, A., Tozer, B., Townend, J., Kurashimo, E., Iidaka, T., and Ishiyama, T., 2013, SAHKE geophysical transect reveals crustal and subduction zone structure at the southern Hikurangi margin, New Zealand: Geochemistry, Geophysics, Geosystems, v. 14(7), p. 2063-2083. Hildyard, M.W., Nippress, S.E., and Rietbrock, A., 2008, Event detection and phase picking using a time-domain estimate of predominate period Tpd: Bulletin of the Seismological Society of America, v. 98, p. 3025-3032. Hollander, M. and Wolfe, D., 1973, Nonparametric Statistical Methods: John Wiley and Sons, Inc. New York, New York, USA. Hong, T.-K., and Kennett, B.L.N., 2003, Modelling of seismic waves in heterogeneous media using a wavelet-based method: application to fault and subduction zones: Geophysical Journal International, v. 154, p. 483-498. Hori, S., Inoue, H., Fukao, Y. and Ukawa, M., 1985, Seismic detection of the untransformed ‘basaltic’oceanic crust subducting into the mantle: Geophysical Journal of the Royal Astronomical Society, v. 83, p. 169-197. Hough, S. E., Anderson, J. G., Brune, J., Vernon, F., Berger, J., Fletcher, J., Haar, L, Hanks, T., and Baker, L., 1988, Attenuation near Anza, California: Bulletin of the Seismological Society of America, v. 78, p. 672-691. Hough, S. E., Ben-Zion, Y., and Leary, P., 1994, Fault zone waves observed at the southern Joshua Tree earthquake rupture zone: Bulletin of the Seismological Society of America, v. 84, p. 761-767. Hough, S., and Dreger, D., 1995, Source parameters of the 23 April 1992 M 6.1 Joshua Tree, California, earthquake and its aftershocks: Empirical Green's function analysis of GEOS and TERRAscope data: Bulletin of the Seismological Society of America, v. 85(6), p. 1576-1590. 172 Hutchings, L., and Wu, F.T., 1990, Empirical Green's functions from small earthquakes: a waveform study of locally recorded aftershocks of the 1971 San Fernando earthquake: Journal of Geophysical Research - Solid Earth, 95(B2), 1187-1214. Ide, S., and Beroza, G.C., 2001, Does apparent stress vary with earthquake size: Geophysical Research Letters, v. 28(17), p. 3349-3352. Ide, S., Beroza, G.C., Prejean, S.G., and Ellsworth, W., 2003, Apparent break in earthquake scaling due to path and site effects on deep borehole recordings: Journal of Geophysical Research - Solid Earth, v. 108(B5). Igel, H., Ben-Zion, Y. and Leary, P., 1997, Simulation of SH and P-SV Wave Propagation in Fault Zones: Geophysical Journal International, v. 128, p. 533-546. Jahnke, G., Igel, H. and Ben-Zion, Y., 2002, Three-dimensional calculations of fault-zone- guided waves in various irregular structures: Geophysical Journal International, v. 151(2), p. 416-426, doi: 10.1046/j.1365-246X.2002.01784.x. Joswig, M., 1995, Automated classification of local earthquake data in the BUG small array: Geophysical Journal International, v. 120, p. 262-286. Julian, B., Miller, A.D., and Foulger, G.R., 1998, Non-double-couple earthquakes 1. Theory: Reviews of Geophysics, v. 36(4), p. 525-549. Jurkevics, A., 1988, Polarization analysis of three-component array data: Bulletin of the Seismological Society of America, v. 78, p. 1725-1743. Kagan, Y., 1991, 3-D rotation of double-couple earthquake sources: Geophysical Journal International, v. 106(3), p. 709-716. Kagan, Y., 1992, Correlations of earthquake focal mechanisms: Geophysical Journal International, v. 110(2), p. 305-320. Kagan, Y., 2000, Temporal correlations of earthquake focal mechanisms: Geophysical Journal International, v. 143, p. 881-897. Kanamori, H. and Anderson, D.L., 1975, Theoretical basis of some empirical relations in seismology: Bulletin of the Seismological Society of America, v. 65(5), p. 1073-1095. Kanamori, H., 2005, Real-time seismology and earthquake damage mitigation: Annual Reviews of Earth and Planetary Science, v. 33, p. 195-214. 173 Kanamori, H., Mori, J., Hauksson, E., Heaton, T.H., Hutton, L.K., and Jones, L.M., 1993, Determination of earthquake energy release and ML using TERRAscope: Bulletin of the Seismological Society of America, v. 83(2), p. 330-346. Kane, D. L., Kilb, D.L., and Vernon, F.L., 2013a, Selecting Empirical Green’s Functions in Regions of Fault Complexity: A Study of Data from the San Jacinto Fault Zone, Southern California: Bulletin of the Seismological Society of America, v. 103(2A), p. 641-650. Kane, D. L., Shearer, P.M., Goertz-Allmann, B.P., and Vernon, F.L., 2013b, Rupture directivity of small earthquakes at Parkfield: Journal of Geophysical Research., v. 118, doi: 10.1029/2012JB009675 Kaneko, Y., and Shearer, P.M., Variability of seismic source spectra, estimated stress drop and radiated energy, derived from cohesive-zone models of symmetrical and asymmetrical circular and elliptical ruptures: Journal of Geophysical Research., v. 120, doi: 10.1002/2014JB011642, 2015. Kasahara, K., Sakai, S., Morita, Y., Hirata, N., Tsuruoka, H., Nakagawa, S., Nanjo, K., and Obara, K., 2010, Development of the Metropolitan Seismic Observation network (MeSO-net) for Detection of Mega-thrust beneath Tokyo Metropolitan Area: Bulletin of the Earthquake Research Institute, University of Tokyo, v. 84, p. 71-88. Kawakatsu, H., 1991, Insignificant isotropic component in the moment tensor of deep earthquakes: Nature, v. 351, p. 50-53. Kawakatsu, H., 1996. Observability of the isotropic component of a moment tensor: Geophysical Journal International, v. 126(2), p. 525-544. Kilb, D., Biasi, G., Anderson, J. G., Brune, J., Peng, Z., and Vernon, F. L., 2012, A comparison of spectral parameter kappa from small and moderate earthquakes using southern California ANZA seismic network data: Bulletin of the Seismological Society of America, v. 102, p. 284-300. Klein, F.W., 2003, The HYPOINVERSE2000 earthquake location program: International Geophysics, v. 81, p. 1619-1620. Kostrov, V.V., 1974, Seismic moment and energy of earthquakes, and seismic flow of rock: Earth Physics, v. 1, p. 23-40. 174 Kuge, K., and Lay, T., 1994, Data-dependent non-double-couple components of shallow earthquake source mechanisms: Effects of waveform inversion instability: Geophysical Research Letters, v. 21(1), p. 9-12. Kuipers, J., 2002, Quaternions and Rotations: a Primer with Applications to Orbits, Aerospace, and Virtual Reality: Princeton University Press, Princeton, NJ, USA. Kuperkoch, L., Meier, T., Lee, J., Friederich, W. and EGELADOS Working Group, 2010, Automated determination of P-phase arrival times at regional and local distances using higher order statistics: Geophysical Journal International, v. 181, p. 1159-1170. Kurzon, I., Vernon, F., Ben-Zion, Y., and Atkinson, G., 2014, Ground motion prediction equations in the San Jacinto fault zone: significant effects of rupture directivity and fault zone amplification: Pure and Applied Geophysics, v. 171(11), p. 3045-3081. Kurzon, I., Vernon, F.L., Rosenberger, A. and Ben-Zion, Y., 2014, Real-time Automatic Detectors of P and S Waves Using Singular Value Decomposition: Bulletin of the Seismological Society of America, v. 104, p. 1696–1708, doi: 10.1785/0120130295. Kwiatek, G. and Ben-Zion, Y., 2013, Assessment of P and S wave energy radiated from very small shear-tensile seismic events in a deep South Africa mine: Journal of Geophysical Research., v. 118, p. 3630–3641, doi: 10.1002/jgrb.50274. Kwiatek, G. and Ben-Zion, Y., 2016, Theoretical limits on detectability of small earthquakes: Journal of Geophysical Research., in review. Kwiatek, G., Plenkers, K., Dresen, G., and JAGUARS Research Group, 2011, Source parameters of picoseismicity recorded at Mponeng deep gold mine, South Africa: implications for scaling relations: Bulletin of the Seismological Society of America, v. 101(6), p. 2592- 2608. Langet, N., Maggi, A., Michelini, A. and Brenguier F., 2014, Continuous Kurtosis-Based Migration for Seismic Event Detection and Location, with Application to Piton de la Fournaise Volcano, La Reunion: Bulletin of the Seismological Society of America, v. 104, p. 229-246, doi:10.1785/0120130107 Lengline, O., and Got, J.L., 2011, Rupture directivity of microearthquake sequences near Parkfield, California: Geophysical Research Letters, v. 38, L08310. Lewis, M.A., and Ben-Zion, Y., 2010, Diversity of fault zone damage and trapping structures in the Parkfield section of the San Andreas Fault from comprehensive analysis of near fault 175 seismograms: Geophysical Journal International, v. 183, p. doi: 10.1111/j.1365- 246X.2010.04816.x. Lewis, M.A., Ben-Zion, Y. and McGuire, J.J., 2007, Imaging the deep structure of the San Andreas Fault south of Hollister with joint analysis of fault zone head and direct P arrivals: Geophysical Journal International, v. 169, p. 1028-1042. Lewis, M.A., Peng, Z., Ben-Zion, Y. and Vernon, F.L., 2005, Shallow Seismic Trapping Structure in the San Jacinto fault zone near Anza, California: Geophysical Journal International, v. 162, p. 867-881, doi:10.1111/j.1365-246X.2005.02684.x. Li Y.G., Aki, K., Adams, D., Hasemi, A., Lee, W.H.K., 1994, Seismic guided waves trapped in the fault zone of the Landers, California, earthquake of 1992: Journal of Geophysical Research, v. 99, p. 11705-11722. Li, Y.G., Leary, P., Aki, K. and Malin, P., 1990, Seismic trapped modes in the Oroville and San Andreas fault zones: Science, v. 249, p. 763-766. Liu, X., Ben-Zion, Y., and Zigone, D., 2015, Extracting seismic attenuation coefficients from cross-correlations of ambient noise at linear triplets of stations: Geophysical Journal International, 203, 1149–1163, doi: 10.1093/gji/ggv357. Lockner, D., Walsh, J., and Byerlee, J., 1977, Changes in Seismic Velocity and Attenuation During Deformation of Granite: Journal of Geophysical Research., v. 82, p. 5374-5378. Ma, K., Lin, Y., Lee, S., Mori, J., and Brodsky, E.E., 2012, Isotropic Events Observed with a Borehole Array in the Chelungpu Fault Zone, Taiwan: Science, v. 337, p. 459-463. Ma, S. and Archuleta, R., 2006, Radiated seismic energy based on dynamic rupture models of faulting: Journal of Geophysical Research - Solid Earth, v. 111(B5). Madariaga, R., 1976, Dynamics of an expanding circular fault: Bulletin of the Seismological Society of America, v. 66(3), p. 639-666. Mamada, Y., Kuwahara, Y., Ito, H. and Takenaka, H., 2004, Discontinuity of the Mozumi– Sukenobu fault low-velocity zone, central Japan, inferred from 3-D finite-difference simulation of fault zone waves excited by explosive sources: Tectonophysics, v. 378, p. 209-222. McGuire, J. and Ben-Zion, Y., 2005, High-resolution imaging of the Bear Valley section of the San Andreas fault at seismogenic depths with fault-zone head waves and relocated seismicity: Geophysical Journal International, v. 163, p. 152-164. 176 McNally, K.C. and McEvilly, T.V., 1977, Velocity contrast across the San Andreas fault in central California: small-scale variations from P-wave nodal plane distortion: Bulletin of the Seismological Society of America, v. 67, p. 1565-1576. Minson, S. E., Dreger, D.S., Bürgmann, R., Kanamori, H., and Larson, K.M., 2007, Seismically and geodetically determined nondouble-couple source mechanisms from the 2000 Miyakejima volcanic earthquake swarm: Journal of Geophysical Research., v. 112, B10308, doi:10.1029/2006JB004847. Mitchell, T. M., Ben-Zion, Y., and Shimamoto, T., 2011, Pulverized fault rocks and damage asymmetry along the Arima-Takatsuki Tectonic Line, Japan: Earth and Planetary Science Letters, v. 308, p. 284-297. Mizuno, T., Kuwahara, Y., Ito, H. and Nishigami, K., 2008, Spatial Variations in Fault-Zone Structure along the Nojima Fault, Central Japan, as Inferred from Borehole Observations of Fault-Zone Trapped Waves: Bulletin of the Seismological Society of America, v. 98, p. 558-570, doi:10.1785/0120060247. Molnar, P., Tucker, B., and Brune, J.N., 1973, Corner frequencies of P and S waves and models of earthquake sources: Bulletin of the Seismological Society of America, v. 63(6-1), p. 2091-2104. Mori, J., and Frankel., A., 1990, Source parameters for small events associated with the 1986 North Palm Springs, California, earthquake determined using empirical Green functions: Bulletin of the Seismological Society of America, v. 80(2), p. 278-295. Mueller, C. S., 1985, Source pulse enhancement by deconvolution of an empirical Green's function: Geophysical Research Letters, v. 12(1), p. 33-36. Mulder, T.L., Kilb, D., White, M., Vernon, F.L., 2015, Analysis of the 28 October 2012 Haida Gwaii Aftershock Sequence Using Automated Earthquake Locations from Singular Value Decomposition Detections: Submitted to Bulletin of the Seismological Society of America Nakamura, Y., 1988, On the urgent earthquake detection and alarm system (UrEDAS): Proceedings of the 9th World Conference on Earthquake Engineering, Japan, Vol.VII, p. 673-678. Nayak, A. and Dreger, D.S., 2014, Moment tensor inversion of seismic events associated with the Sinkhole at Napoleonville Salt Dome, Louisiana: Bulletin of the Seismological Society of America, v. 104(4), p. 1763-1776. 177 Nippress, S.E.J., Rietbrock, A. and Heath, A.E., 2010, Optimized automatic pickers: application to the ANCORP data set: Geophysical Journal International, v. 181, p. 911-925. Olsen, K., Day, S., Minster, J., Cui, Y., Chourasia, A., Faerman, M., Moore, R., Maechling, P., and Jordan, T.H., 2006, Strong shaking in Los Angeles expected from southern San Andreas earthquake: Geophysical Research Letters, v. 33(7). Omori, F., 1894, On after-shocks of earthquakes: J. Coll. Sci. Imp. Univ. Tokyo, v. 7, p. 111- 200. Oppenheimer, D.H., Reasenberg, P.A., and Simpson, R.W., 1988, Fault plane solutions for the 1984 Morgan Hill, California, earthquake sequence: Evidence for the state of stress on the Calaveras fault: Journal of Geophysical Research., v. 93, p. 9007-9026. Oth, A., 2013, On the characteristics of earthquake stress release variations in Japan: Earth and Planetary Science Letters, v. 377, p. 132-141. Ozakin, Y., Ben-Zion, Y., Aktar, M., Karabulut, H., and Peng, Z., 2012, Velocity contrast across the 1944 rupture zone of the North Anatolian fault east of Ismetpasa from analysis of teleseismic arrivals: Geophysical Research Letters, v. 39, L08307, doi:10.1029/2012GL051426. Park, J., Lindberg, C.R., and Vernon, F.L., 1987, Multitaper Spectral Analysis of High-Frequency Seismograms: Journal of Geophysical Research.: v. 92(B12), p. 12,675-12,684. Patton, H. J., and Taylor, S.R., 2011, The apparent explosion moment: Inferences of volumetric moment due to source medium damage by underground nuclear explosions: Journal of Geophysical Research., v. 116, B03310, doi:10.1029/2010JB007937. Pavlis, G. L., Vernon, F.L., Harvey, D., and Quinlan, D., 2004, The generalized earthquake- location (GENLOC) package: an earthquake-location library: Computers and Geosciences, v. 30 p. 1079-1091, doi:10.1016/j.cageo.2004.06.010. Pearce, R. G., and Rogers, R.M., 1989, Determination of earthquake moment tensors from teleseismic relative amplitude observations: Journal of Geophysical Research., v. 94(B1), 775–786, doi:10.1029/JB094iB01p00775. Peng, Z., and Zhao, P., 2009, Migration of early aftershocks following the 2004 Parkfield earthquake: Nature Geoscience, v. 2, p. 877-881. 178 Peng, Z., Ben-Zion, Y., Michael, A.J. and Zhu, L., 2003, Quantitative analysis of seismic trapped waves in the rupture zone of the 1992 Landers, California earthquake: Evidence for a shallow trapping structure: Geophysical Journal International, v. 155, p. 1021-1041. Prieto, G., Parker, R., and Vernon, F., 2009, A Fortran 90 library for multitaper spectrum analysis: Computers and Geosciences, v. 35(8), p. 1701-1710. Prieto, G., Shearer, P., Vernon, F., and Kilb, D., 2004, Earthquake source scaling and self- similarity estimation from stacking P and S spectra: Journal of Geophysical Research - Solid Earth, v. 109(B8), B08310, doi: 10.1029/2004JB003084. Pujol, J., and Herrmann, R., 1990, A student's guide to point sources in homogeneous media: Seismological Research Letters, v. 61(3-4), p. 209-221. Qiu, H., Ben-Zion, Y., Ross, Z. E., Share, P.-E. and Vernon, F, 2015, Internal structure of the San Jacinto fault zone at Jackass Flat from data recorded by a dense linear array: Abstract of the annual SSA meeting. Okaya, D., Christensen, N., Ross, Z. E., and F. T. Wu, 2015, Terrane-controlled crustal shear wave splitting in Taiwan: Geophysical Research Letters, 43(2), doi: 10.1002/2015GL066446 Ranjith, K. and Rice, J., 2001, Slip dynamics at an interface between dissimilar materials: Journal of Mechanics of Physical Solids, v. 49, p. 341-361. Rastin, S.J., Unsworth, C.P., Benites, R., and Gledhill, K.R., 2013, Using Real and Synthetic Waveforms of the Matata Swarm to Assess the Performance of New Zealand GeoNet Phase Pickers: Bulletin of the Seismological Society of America, v. 103, p. 2173-2187. Rawles, C., and Thurber, C., 2015, A nonparametric method for automatic determination of P- wave and S-wave arrival times: Application to local microearthquakes: Geophysical Journal International, v. 202 (2), p. 1164-1179, doi: 10.1093/gji/ggv218. Reid, H.F., 1910, The mechanics of the earthquake, in The California Earthquake of April 18, 1906: Vol. 2, Carnegie Institute, Washington, DC, USA. Rice, J.R., 2006, Heating and weakening of faults during earthquake slip: Journal of Geophysical Research., v. 111, B05311, doi:10.1029/2005JB004006 Richter, C.F., 1935, An instrumental earthquake magnitude scale: Bulletin of the Seismological Society of America, v. 25, p. 1-32. 179 Riedesel, M. and Jordan, T., 1989: Display and assessment of seismic moment tensors: Bulletin of the Seismological Society of America, v. 79(1), p. 85-100. Rosenberger, A., 2010, Real-Time Ground-Motion Analysis: Distinguishing P and S Arrivals in a Noisy Environment: Bulletin of the Seismological Society of America, v. 100, p. 1252- 1262. Ross, Z. E., and Ben-Zion, Y., 2013, Spatio-temporal variations of double-couple aftershock mechanisms and possible volumetric earthquake strain: Journal of Geophysical Research., v. 118, p. 2347–2355, doi:10.1002/jgrb.50202 Ross, Z. E. and Ben-Zion, Y., 2014a, Automatic picking of direct P, S seismic phases and fault zone head waves: Geophysical Journal International, v. 199, p. 368–381, doi: 10.1093/gji/ggu2672014. Ross, Z.E., and Ben-Zion, Y., 2014b, An automatic detection algorithm with pseudo- probabilities of multiple indicators: Geophysical Journal International, doi: 10.1093/gji/ggt516 Ross, Z. E., Ben-Zion, Y., and Zhu, L., 2015, Isotropic source terms of San Jacinto fault zone earthquakes based on waveform inversions with a generalized CAP method: Geophysical Journal International, v. 200(2), p. 1269-1280. Ross, Z. E., and Ben-Zion, Y., 2015, An algorithm for automated identification of fault zone trapped waves: Geophysical Journal International, v. 202 (2), p. 933–942, doi:10.1093/gji/ggv197. Ross, Z. E. and Ben-Zion, Y., 2016, Towards reliable automated estimates of earthquake source properties from body wave spectra: Journal of Geophysical Research: Solid Earth, in review. Ross, Z. E., White, M.C, Vernon, F.L., and Ben-Zion, Y., 2016, An improved algorithm for real- time S-wave picking with application to the (augmented) ANZA network in southern California: Bulletin of the Seismological Society of America, in review. Ross, Z. E., Ben-Zion, Y., White, M. C., and Vernon, F.L., 2016, Analysis of earthquake body wave spectra for potency and magnitude values: Implications for magnitude scaling relations, Geophysical Journal International, in review. 180 Rovelli, A., Caserta, A., Marra, F., and Ruggiero, V., 2002, Can Seismic Waves Be Trapped inside an Inactive Fault Zone? The Case Study of Nocera Umbra, Central Italy: Bulletin of the Seismological Society of America, v. 92, p. 2217-2232. Rubin, A. and Gillard, D., 2000, Aftershock asymmetry/rupture directivity along central San Andreas fault microearthquakes: Journal of Geophysical Research., v. 105, p. 19,095- 19,109. Saragiotis, C., Hadjileontiadis, L. and Panas, S., 2002, PAI-S/K: A robust automatic seismic P phase arrival identification scheme, IEEE Transactions on Geoscience and Remote Sensing, v. 40, p. 1395-1404. Sato, T. and Hirasawa, T., 1973, Body wave spectra from propagating shear cracks: Journal of the Physics of the Earth, v. 21, p. 415-431 SCEDC, 2013, Southern California Earthquake Center. California Institute of Technology. Dataset. doi: 10.7909/C3WD3xH1 Scholz, C. H., Aviles, C., and Wesnousky, S.G., 1986, Scaling differences between large interplate and intraplate earthquakes: Bulletin of the Seismological Society of America, v. 76(1), p. 65-70. Seeber L., Armbruster J.G., Ozer N., Aktar M., Baris S., Okaya D., Ben-Zion, Y. and Field, E., 2000, The 1999 Earthquake Sequence along the North Anatolia Transform at the Juncture between the Two Main Ruptures, in The 1999 Izmit and Duzce Earthquakes: preliminary results, edit. Barka A., O. Kazaci, S. Akyuz and E. Altunel, Istanbul technical university, 209-223. Shapiro, N. M., Olsen, K. B. and Singh, S. K., 2000, Wave-guide effects in subduction zones: Evidence from three-dimensional modeling: Geophysical Research Letters, v. 27(3), p. 433-436. Share, P.-E., Ben-Zion, Y., Ross, Z.E., Qiu, H. and Vernon, F., 2015, Characterization of the San Jacinto Fault Zone northwest of the trifurcation area from earthquake data recorded by a dense linear array: Abstract of the annual SSA meeting. Shearer, P. M., Prieto, G.A., and Hauksson, E., 2006, Comprehensive analysis of earthquake source spectra in southern California: Journal of Geophysical Research - Solid Earth, v. 111(B6). 181 Shi, Z. and Ben-Zion, Y., 2006, Dynamic rupture on a bimaterial interface governed by slip- weakening friction: Geophysical Journal International, v. 165, p. 469-484. Sibson, R. H., 1986, Rupture interactions with fault jogs: Earthquake Source Mechanics, Geophys. Monogr. Ser., S. Das, J. Boatwright, and C. H. Scholz, eds., AGU, Washington, DC, 37, pp. 157–167. Silver, P.G., and Jordan, T.H., 1983, Total-Moment Spectra of Fourteen Large Earthquakes: Journal of Geophysical Research, v. 88, p. 3273-3293. Sleeman, R. and van Eck, T., 1999, Robust automatic P-phase picking: an on-line implementation in the analysis of broadband seismogram recordings: Physics of the Earth and Planetary Interiors, v. 113, p. 265-275. Snoke, J., 1987, Stable determination of (Brune) stress drops: Bulletin of the Seismological Society of America, v. 77(2), p. 530-538. Stanchits, S., Vinciguerra, S., and Dresen, G., 2006, Ultrasonic velocities, acoustic emission characteristics and crack damage of basalt and granite: Pure and Applied Geophysics, v. 163, p. 975-994. Stierle, E., Bohnhoff, M., and Vavryčuk, V., 2014a, Resolution of non-double-couple components in the seismic moment tensor using regional networks—II: application to aftershocks of the 1999 Mw 7.4 Izmit earthquake: Geophysical Journal International, v. 196(3), p. 1878-1888. Stierle, E., Vavryčuk, V., Šílený, J., and Bohnhoff, M., 2014b, Resolution of non-double-couple components in the seismic moment tensor using regional networks—I: a synthetic case study: Geophysical Journal International, v. 196(3), p. 1869-1877. Tape, W., and Tape, C., 2013, The classical model for moment tensors: Geophysical Journal International., doi: 10.1093/gji/ggt302. Thomson, D. J., 1982, Spectrum estimation and harmonic analysis: Proceedings of the IEEE, v. 70(9), p. 1055-1096. Thurber, C., Zhang, H., Waldhauser, F., Hardebeck, J., Michael, A. and Eberhart-Phillips, D., 2006, Three-dimensional compressional wavespeed model, earthquake relocations, and focal mechanisms for the Parkfield, California, region: Bulletin of the Seismological Society of America, v. 96, S38-S49. 182 Utsu, T., 1999, Representation and analysis of the earthquake size distribution: a historical review and some new approaches: in Seismicity Patterns, their Statistical Significance and Physical Meaning, p. 509-535, Springer. Vavryčuk, V., 2014, Moment tensor decompositions revisited: Journal of Seismology, v. 18(4), doi: 10.1007/s10950-014-9463-y Venkataraman, A., and Kanamori, H., 2004, Observational constraints on the fracture energy of subduction zone earthquakes: Journal of Geophysical Research - Solid Earth, v. 109(B5). Vidale, J. E., 1986, Complex Polarization Analysis of Particle Motion: Bulletin of the Seismological Society of America, v. 76 (5) p. 1393-1405. Wald, D., and Heaton, T., 1994, Spatial and temporal distribution of slip for the 1992 Landers, California, earthquake: Bulletin of the Seismological Society of America, v. 84, p. 668- 691. Waldhauser, F., and Ellsworth, W.L., 2000, A double-difference earthquake location algorithm: Method and application to the northern Hayward fault, California: Bulletin of the Seismological Society of America, v. 90, p. 1353-1368. Walter, W.R., Mayeda, K., and Patton, H., 1995, Phase and spectral ratio discrimination between NTS earthquakes and explosions, Part I: Empirical Observations: Bulletin of the Seismological Society of America, v. 85, p. 1050-1067. Wang, J., and Teng, T., 1997, Identification and picking of S phase using an artificial neural network: Bulletin of the Seismological Society of America, v. 87, p. 1140-1149. Wechsler, N., Rockwell, T.K., and Ben-Zion, Y., 2009, Application of high resolution DEM data to detect rock damage from geomorphic signals along the central San Jacinto Fault: Geomorphology, v. 113, p. 82-96, doi:10.1016/j.geomorph.2009.06.007. Weertman, J., 1980, Unstable Slippage Across a Fault that Separates Elastic Media of Different Elastic-Constants: Journal of Geophysical Research., v. 85, p. 1455-1461. Withers, M., Aster, R., Young, C., Beiriger, J., Harris, M., Moore, S. and Trujillo, J., 1998, A comparison of select trigger algorithms for automated global seismic phase and event detection: Bulletin of the Seismological Society of America, v. 88, p. 95-106. Wu, F.T., Ross, Z.E., Okaya, D., Ben-Zion, Y., Wang, C.-Y., Kuo-Chen, H., and Liang, W., 2015, Dense Network, Intense Seismicity and Tectonics of Taiwan, Tectonophysics, in review. 183 Wyss, M. and Brune, J.N., 1968, Seismic moment, stress, and source dimensions for earthquakes in the California-Nevada region: Journal of Geophysical Research., v. 73(14), p. 4681- 4694. Xu, H., Rodgers, A.J., Lomov, I.N., and Vorobiev, O.Y., 2012, Seismic source characteristics of nuclear and chemical explosions in granite from hydrodynamic simulations: Pure and Applied Geophysics, doi: 10.1007/s00024-012-0623-0. Yang, H., Zhu, L. and Chu, R., 2009, Fault-plane determination of the 18 April 2008 Mount Carmel, Illinois, earthquake by detecting and relocating aftershocks: Bulletin of the Seismological Society of America, v. 99, p. 3413-3420. Yang, W. and Ben-Zion, Y., 2016, Corner frequency ratios of P and S waves and strain drops of earthquakes recorded by a tight network around the Karadere segment of the North Anatolian Fault Zone: Evidence for non-classical source processes: Geophysical Journal International., v. 205, p. 220–235, doi: 10.1093/gji/ggv560. Yang, W., Hauksson, E., and Shearer, P., 2012, Computing a Large Refined Catalog of Focal Mechanisms for Southern California (1981–2010): Temporal Stability of the Style of Faulting: Bulletin of the Seismological Society of America, v. 102(3), p. 1179-1194. Yang, W., Peng, Z. and Ben-Zion, Y., 2009, Variations of strain drops in aftershocks of the 1999 İzmit and Düzce earthquakes along the Karadere-Düzce branch of the North Anatolian fault: Geophysical Journal International, v. 177, p. 235-246. Zaliapin, I., and Ben-Zion, Y., 2011, Asymmetric distribution of aftershocks on large faults in California: Geophysical Journal International, v. 185, p. 1288-1304. Zaliapin, I., and Ben-Zion, Y., 2013a, Earthquake clusters in southern California I: Identification and stability, Journal of Geophysical Research - Solid Earth, 118(6), 2847-2864. Zaliapin, I., and Ben-Zion, Y., 2013b, Earthquake clusters in southern California II: Classification and relation to physical properties of the crust: Journal of Geophysical Research - Solid Earth, v. 118(6), 2865-2877. Zhao, L., and Helmberger, D.V., 1994, Source estimation from broadband regional seismograms: Bulletin of the Seismological Society of America, v. 84(1), p. 91-104. Zhao, P., Peng, Z., Shi, Z., Lewis, M.A., and Ben-Zion, Y., 2010, Variations of the velocity contrast and rupture properties of M6 earthquakes along the Parkfield section of the San Andreas fault: Geophysical Journal International, v. 180, p. 765-780. 184 Zhu, L., and Ben-Zion, Y., 2013, Parameterization of general seismic potency and moment tensors for source inversion of seismic waveform data: Geophysical Journal International, v. 194, p. 839-843, doi:10.1093/gji/ggt137. Zhu, L., and Helmberger, D.V., 1996, Advancement in source estimation techniques using broadband regional seismograms: Bulletin of the Seismological Society of America, v. 86(5), p. 1634-1641. Zhu, L., and Rivera, L., 2002, A note on the dynamic and static displacements from a point source in multi-layered media: Geophysical Journal International, v. 148, p. 619-627, doi:10.1046/j.1365-246X.2002.01610.x. Zigone, D., Ben-Zion, Y., Campillo, M., and Roux, P., 2015, Seismic tomography of the Southern California plate boundary region from noise-based Rayleigh and Love waves: Pure and Applied Geophysics, v. 172(5), p. 1007-1032.

Abstract (if available)

Abstract I develop a collection of automated techniques suitable for systematic processing of large seismic waveform data sets. The techniques allow for extracting event, phase, and source information from the data sets with the goal of constructing databases for downstream analyses. I propose algorithms for real-time automated picking of P- and S-wave arrivals that are capable of detecting earthquakes at the same time. The methods apply polarization filters to the 3-component data to enhance signal to noise ratios and use moving averages in tandem with a kurtosis-based detector to lock in on the phase arrival. The techniques are used to construct earthquake catalogs for Southern California from scratch by re-processing the original waveform archives. ❧ A method is developed for detecting earthquake signals using pseudo-probabilities, by combining different indicator variables together in a non-parametric manner. The method is tested on data recorded by the Southern California Seismic Network and is found to detect more than three times as many earthquakes than the network’s catalog for the examined period of time. ❧ Methods are proposed for automated identification of fault zone head and trapped waves. These are special waves that are observed exclusively in fault zones and carry high-resolution information on the nature of the geologic structures which generated them. The techniques are applied to thousands of events recorded at the Parkfield section of the San Andreas fault, the Landers rupture area, the San Jacinto fault zone, and the Hayward fault. The identified head waves are used to estimate the extent of the velocity contrast across the faults, and the identified trapped waves are used to estimate various properties of the damage zones. ❧ A general-purpose methodology is developed for calculating numerous source quantities from body-wave spectral ratios including stress/strain drop, radiated seismic energy, and the extent of directivity. The technique is demonstrated on a set of five earthquakes with M > 3.4 in Southern California on hundreds of stations. Roughly 30,000 seismograms are examined in the process with no user involvement. Four of the five earthquakes show evidence of unilateral directivity. ❧ A catalog of more than 13,000 earthquakes is constructed for the San Jacinto fault zone and used to study the scaling relationship between seismic potency and local magnitude. Moment magnitudes are calculated for more than 11,000 local earthquakes using P- and S-wave spectra. It is found that moment magnitudes are statistically larger than local magnitudes for the same set of events. Further, the b-value measured from moment magnitudes is 25% larger than the b-value measured from local magnitudes. ❧ Full moment tensors are determined for a set of 7 earthquakes in the San Jacinto fault zone. Small, explosive type isotropic components are measured in all of the events. A set of tests are performed in which the station distribution and velocity model are perturbed to assess the robustness of these values. It is found that the CLVD components change sign frequently in this process and are likely the result of artifacts. The the isotropic components persist at being explosive in type for all of the events in more than 90 percent of the perturbation scenarios. ❧ Double couple focal mechanisms are studied in different regions of the Landers rupture zone. I find that the short term aftershocks in the edges of the Landers rupture are less aligned with the mainshock focal mechanism than in the central portion of the rupture. Over time, the statistical orientation of the average focal mechanism approaches the long-term background mechanism. It is shown that a possible explanation for these observations is small isotropic components present in the early aftershocks of the edge regions that occasionally flips polarities measured on nodal planes.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Detailed properties of seismic waveforms and earthquake sources in Southern California

PDF

Multi-scale imaging and monitoring of crustal and fault zone structures in southern California

PDF

Stress-strain characterization of complex seismic sources

PDF

Determining fault zone structure and examining earthquake early warning signals using large datasets of seismograms

PDF

Volumetric interactions between major ruptures and fault zones illuminated by small earthquake properties

PDF

High-resolution imaging and monitoring of fault zones and shallow structures: case studies in southern California and on Mars

PDF

Elements of seismic structures near major faults from the surface to the Moho

PDF

Multi-scale imaging of major fault zones in Southern California

PDF

Heterogeneity of earthquake stress drops, focal mechanisms and active fault zones

PDF

Multi-scale imaging of the fault zone velocity structure: double-difference tomography, inversion of fault zone headwaves, and fault zone sensitivity kernels

PDF

Analysis of waveform and catalog data of aftershocks for properties of earthquakes and faults

PDF

Observations and modeling of dynamically triggered high frequency burst events

PDF

Statistical analyses of ambient seismic noise spectra with applications to detecting imperfectly diffuse wave field and deriving attenuation and phase velocity information

PDF

Seismic velocity structure and spatiotemporal patterns of microseismicity in active fault zones of Southern California

PDF

Symmetry properties, pulverized rocks and damage architecture as signatures of earthquake ruptures

PDF

Microseismicity, fault structure, & the seismic cycle: insights from laboratory stick-slip experiments

$Seismic radiation from an explosion in rock having an anisotropic distribution of fractures$

PDF

Seismic radiation from an explosion in rock having an anisotropic distribution of fractures

PDF

Quantifying ground deformation of large magnitude earthquakes using optical imgaging systems

PDF

Seismicity distribution near strike-slip faults in California

PDF

Spatiotemporal variations of stress field in the San Jacinto Fault Zone and South Central Transverse Ranges

Asset Metadata

Creator Ross, Zachary Elias (author)

Core Title Applying automated techniques to large seismic datasets for systematic analyses of phases, source, and structure

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Geological Sciences

Publication Date 06/16/2016

Defense Date 04/07/2016

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

Tag body waves,earthquake detection,earthquake source,moment tensor,OAI-PMH Harvest,phase picking,seismicity,seismology

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Ben-Zion, Yehuda (committee chair), Alexander, Kenneth (committee member), Jordan, Thomas H. (committee member), Sammis, Charles G. (committee member)

Creator Email zross@gps.caltech.edu,zross@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-251908

Unique identifier UC11280661

Identifier etd-RossZachar-4439.pdf (filename),usctheses-c40-251908 (legacy record id)

Legacy Identifier etd-RossZachar-4439.pdf

Dmrecord 251908

Document Type Dissertation

Format application/pdf (imt)

Rights Ross, Zachary Elias

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 a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA