University of Southern California Dissertations and Theses

Detailed properties of seismic waveforms and earthquake sources in Southern California

(USC Thesis Other)

Detailed properties of seismic waveforms and earthquake sources in Southern California

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Doctoral Dissertation Haoran Meng

1

Detailed properties of seismic waveforms and earthquake
sources in Southern California

by

Haoran Meng

A Dissertation Presented to the
Faculty of the University of Southern California
Graduate School

In Partial fulfillment of the requirements for the degree of
Doctor of Philosophy
in
Geological Sciences

August 2019

Doctoral Dissertation Haoran Meng

2
Acknowledgements

I would like to take this opportunity to express my gratitude to my Ph.D. advisor,
Professor Yehuda Ben-Zion for his supervision, advice and guidance. Yehuda’s inspiration,
intuition and passion have made this work impressive and joyful learning experience. It is a
privilege to work with and learn from him.
Sincere appreciations are given to Jeffery McGuire for generously sharing his ideas and
knowledge and warmly hosting me during my staying in Woods Hole.
I am grateful that William Frank, Ken Alexander, Julien Emile-Geay, Charles Sammis,
Aiichiro Nakano and Francis Wu took the time out of their busy schedules to help advance
my career, having served on my qualifying exam or 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. Thanks to Cindy Waite, John Yu, John McRaney
and Miguel Rincon for helping me patiently on matters related to the graduate program and
logistics.
I would also like to thank my research collaborators: Christopher Johnson and Jeffery
McGuire. I have learned a great deal from all of you.
Lastly, I would like to convey my deep gratitude to my family and friends, for all their love
and encouragement. They have made my life colorful and unforgettable.

Doctoral Dissertation Haoran Meng

3
Table of Contents

Abstract ................................................................................................................................................................. 5
Introduction ........................................................................................................................................................ 7
1. Detection of small earthquakes with dense array data: Example from the San Jacinto
fault zone, Southern California ........................................................................................................ 12
1.1 Introduction ..................................................................................................................................... 12
1.2 Method and Results ....................................................................................................................... 14
1.3 Discussion and conclusions ....................................................................................................... 24
Appendix 1.A ........................................................................................................................................... 27
Appendix 1.B ........................................................................................................................................... 36
Appendix 1.C............................................................................................................................................ 37
2. Characteristics of airplanes and helicopters recorded by a dense seismic array near
Anza California ....................................................................................................................................... 39
2.1 Introduction ..................................................................................................................................... 39
2.2 Analysis .............................................................................................................................................. 44
2.3 Discussion and conclusions ....................................................................................................... 53
Appendix 2.A ........................................................................................................................................... 57
Appendix 2.b ............................................................................................................................................ 59
3. Detection of random noise and anatomy of continuous seismic waveforms in dense
array data near Anza California ....................................................................................................... 65
3.1 Introduction ..................................................................................................................................... 65
3.2 Analysis .............................................................................................................................................. 69
3.3 Anatomy of recorded waveforms ............................................................................................ 78
3.4 Discussion ......................................................................................................................................... 80
Doctoral Dissertation Haoran Meng

4
4. Evaluation of diffuse waveforms in dense array data near Anza California ................. 88
4.1 Introduction ..................................................................................................................................... 88
4.2 Method ............................................................................................................................................... 89
4.3 Results and discussion ................................................................................................................. 93
5. Semi-automated estimates of directivity and related source properties of small to
moderate Southern California earthquakes with second seismic moments ...............105
5.1 Introduction ...................................................................................................................................105
5.2 Data ....................................................................................................................................................107
5.3 Method ..............................................................................................................................................109
5.4 Discussion and conclusions......................................................................................................125
Appendix 5.A .........................................................................................................................................128
References ......................................................................................................................................................136

Doctoral Dissertation Haoran Meng

5
Abstract
I develop a collection of methodologies detecting small earthquakes, modeling air-traffic
events, separating continuous seismic waveforms into random noise, not-random-noise,
evaluating diffuse waveform and estimating of directivity, rupture area, duration, and
centroid velocity of earthquakes using second seismic moments and stacked empirical
Green’s functions. These techniques are suitable for systematic processing of large seismic
waveform data sets to extract detailed properties of seismic waveforms and earthquake
sources.
A technique is developed to detect small earthquakes not included in standard catalogs
using data from a dense seismic array. The technique is illustrated with continuous
waveforms recorded in a test day by 1108 vertical geophones in a tight array on the San
Jacinto fault zone. The number of validated detected events in the test day is > 5 times larger
than that in the standard catalog. Using these events as templates can lead to additional
detections of many more earthquakes.
Frequent air-traffic events are observed in continuous seismic waveforms recorded for
about 30 days in the dense array on the San Jacinto fault zone. Time-frequency analysis
reveals clear Doppler effects that can be modeled with basic equations and fitted well with
parameters corresponding to airplanes and helicopters. The total event duration covers > 7%
of the day, which is likely to exceed considerably the total time covered by earthquakes.
A methodology is developed to separate continuous seismic waveforms into random noise
(RN), not-random-noise (NRN) produced by earthquakes, wind, traffic and other sources,
and the undetermined mixture of signals. The analysis is applied to continuous records of the
dense seismic array on the San Jacinto fault zone. For the days examined, the relative
fractions of RN, NRN and mixed signals in local day (night) times are about 26% (42%), 40%
(33%) and 34% (25%), respectively.
A technique is developed to evaluate whether the recorded waveforms are diffused or not
by examining the expectation of waveform spectra and the outer products with themselves
and their conjugates. The method is applied to the detected RN and earthquake coda. The
results show the detected RN are well diffuse while the recorded waveforms in the coda
window of the 2010 Mw 7.2 El Mayor Cucapah are not necessarily diffuse because of the
triggered events or instrumental effect.
Doctoral Dissertation Haoran Meng

6
A method is developed for semi-automated estimation of directivity, rupture area,
duration, and centroid velocity of earthquakes with second seismic moments. The method is
applied to 14 earthquakes with magnitude in the range 3.5-5.2 in Southern California. Most
analyzed target events in the Trifurcation area of the San Jacinto Fault system have
directivities towards the northwest, while events around Cajon Pass and San Gabriel
Mountain propagate towards the southeast. The results are generally consistent with
predictions for dynamic rupture on a bimaterial interface and the imaged velocity contrasts
in the study area. The second moment inversions are also used to explore constraints on
upper and lower bounds of rupture areas in our dataset. We note that the stress drop can be
significantly underestimated by assuming a circular crack. The second moment method with
stacked eGfs appears capable of semi-automated, routine application to moderate
earthquakes in southern California.

Doctoral Dissertation Haoran Meng

7
Introduction
In the last decade it became a common practice to record seismic waveforms continuously.
The continuous recordings, combined with the increasing number of seismometers and
array deployments, provide important new opportunities for detecting small events, thereby
increasing the ability to resolve seismic processes and structural properties. One highly
successful technique is the matched-filter processing (or template matching) that relies on
cross-correlating seismograms of events detected by standard techniques with the
continuous waveforms (e.g., Gibbons & Ringdal 2006). Wave portions with high similarity to
the templates on multiple stations are assumed to reflect smaller events located very close
to the event templates. This method enables detection of events with Signal-to-Noise-Ratio
(SNR) below 1, and leads to 10-100 times more events than those listed in typical catalogs
(e.g., Meng et al. 2013; Shelly et al. 2016; Ross et al. 2017). However, a basic shortcoming of
this technique is that it allows only detection of events with hypocenters very close to events
already existing in the catalog.
The growing quantity and quality of continuous waveforms recorded by seismic networks
and temporary arrays at various scales, coupled with improved sensors and processing
techniques, provide increasing opportunities for detection and analyses of weak sources of
ground motion. Detection and analyses of events with signal-to-noise (SNR) ratios below 1
is now commonly done in studies of earthquakes, tremor and induced seismicity around the
world (e.g., Shelly et al., 2007; Zigone et al., 2012; Ross et al., 2017; Chen et al., 2017). The
association of weak signals with SNR ≤ 1 with genuine seismic sources below the surface
(earthquakes, tremor) requires detailed understanding of the differences between such
events and sources of energy at and above the surface. These include various types of cultural
noise at the surface (e.g., cars and trains), interactions of the wind with surface structures
and sources in the atmosphere (e.g., storms and air-traffic).
In addition to developing reliable catalogs of weak sub-surface seismic events, the ability
to detect and classify different sources that produce small episodic ground motion is
important for tomographic imaging based on the ambient seismic noise. These studies,
which are now done routinely in many environments (e.g., Campillo et al., 2011; Nakata et
al., 2015), require a diffuse wavefield (e.g., Weaver, 1982; Shapiro et al., 2000) and removal
from observed waveforms of all other sources (e.g., Bensen et al., 2007). More generally,
Doctoral Dissertation Haoran Meng

8
detailed characterization of different types of common sources of ground motion contributes
to a more complete understanding of the observed seismic waveforms.
Seismic ground motion is generated by a variety of mechanisms including tectonic sources
such as earthquakes and tremor (e.g., Ross et al. 2017; Shelly et al. 2007), interactions of
ocean waves with the solid Earth (e.g., Gerstoft & Tanimoto 2007; Hillers et al. 2013), storms
and wind shaking of structures, trees and other obstacles (e.g., Tanimoto & Valovcin 2015,
Johnson et al. 2019) and various anthropogenic sources such as trains, cars, airplanes,
helicopters and wind turbines (e.g., Salvermoser et al. 2015; Eibl et al. 2015; Neuffer &
Kremers 2017; Meng & Ben-Zion 2018b; Inbal et al. 2018). It is generally assumed that
continuous waveforms processed to remove impulsive sources (e.g. Bensen et al. 2007)
consist primarily of a random wavefield, referred to as the ambient seismic noise, that is
suitable generally for imaging and monitoring of sub-surface structures (e.g., Shapiro et al.
2005; Brenguier et al. 2008; Zigone et al. 2015). Extracting accurate Green’s functions from
the ambient seismic noise requires a stronger condition that the wavefield is fully diffuse in
addition to being random (e.g., Weaver 1982; Shapiro et al. 2000; Sá nchez-Sesma et al., 2008).
Seismic records are convolutions of source, propagation path and instrument response.
To estimate the source properties such directivity and rupture size, the path and instrument
effects should be removed from the data. This is often done by deconvolving seismograms of
the target earthquake with its empirical Green’s function (eGf) given by seismograms of a
suitable small-collocated event (e.g. Mueller, 1985; Hutchings & I, 1990; Hough & Dreger,
1995; Ross & Ben-Zion, 2016). An eGf is often required to be collocated with the target event
and to have a Focal Mechanism (FM) very close to the event targeted for time domain
analysis. In this paper, we develop a semi-automated technique using a weighted stack of
eGfs to get a more accurate eGf for deriving source properties of small to moderate
earthquakes by analysis of the second-degree seismic moments.
Earthquake catalogs (e.g. Hauksson et al. 2012) typically represent events as point sources
with location, time, magnitude and focal mechanism. The “cut and paste” technique was
developed to estimate the centroid location, time, focal mechanism using broadband
regional seismograms (e.g. Zhao & Helmberger, 1993; Zhu & Helmberger, 1996). The centroid
location and time are the spatial and temporal mean (first degree moment) of the earthquake
moment release distribution and routinely cataloged from both regional and global seismic
Doctoral Dissertation Haoran Meng

9
networks. The simplest general representation of the finite source information about
rupture extent and directivity is the second-degree moments (variances) of the moment-
release distribution (e.g. Backus & Mulcahy, 1976a, b; Backus, 1977a, b; McGuire, 2004;
McGuire, 2017). This important information, however, is not yet included in the standard
catalogs due to a lack of automated estimation techniques for the second moments.
For simplicity, many studies assume that earthquake kinematic properties, such as
rupture velocity and stress drop, are generally similar everywhere. However, this
assumption can suppress recognition of persistent variations of source processes related to
different properties of fault zones. For instance, systematic directivity of earthquakes
ruptures on bimaterial faults has been proposed by various dynamic simulation and
observations (e.g., Weertman, 1980; Ben-Zion & Andrews, 1998; Ampuero & Ben-Zion, 2008;
Brietzke et al., 2009). Directivities of large earthquakes can significantly amplify the ground
motion and therefore cause serious damage at sites in the forward direction.
Calderoni et al. (2012, 2015) showed that if eGf events are large enough to contain
directivity effects in the frequency band of interest, that can result in apparent directivity of
the target event. Therefore, using a single eGf can bias the results if the eGf has non-generic
source effects such as directivity in the band of interest. Although eGfs can effectively account
for the path and site effects, it is hard to find a perfect eGf with collocated hypocenter and
exact identical FM to a target earthquake. The difference of the FMs between the target event
and its eGf can also result in bad station coverage and biased estimation of the source
proprieties.
The dissertation is composed of five chapters: two of which are published in peer-
reviewed journals, one is currently in revision and the remaining two are in preparation.
Chapters 1 to 4 I discuss methodologies for systematic analyses waveforms recorded
mainly by a dense seismic array with 1108 vertical ZLand geophones (10 Hz) centered on
the Clark branch of the San Jacinto fault zone (SJFZ) at the Sage Brush Flat (SGB) site
southeast of Anza, California. The array recorded continuously at 500 Hz from 7 May 2014
to 13 June 2014, and covered an area of about 600 m × 600 m with nominal sensor spacing
of 10 m normal to the fault and 30 m along-strike (Ben-Zion et al., 2015). Addition surface
and borehole broadband stations in Southern California are analyzed for validations. The
methodologies include separations continuous seismic waveforms into random noise (RN),
Doctoral Dissertation Haoran Meng

10
not-random-noise (NRN) produced by earthquakes, wind, traffic and other sources, and
undetermined mixture of signals. For the days examined, the relative fractions of RN, NRN
and mixed signals in local times are about on 33%, 42% and 25% average (Chapter 3). As
part of an effort to understand the main contributions to NRN observed at the SGB site, Meng
& Ben-Zion (2018b, Chapter 2) estimate the duration of earthquake likely cover <1% of a day
in Southern California. They also analyze waveforms generated by air-traffic events. The
time-frequency analysis reveals clear Doppler effects that can be modeled with basic
equations and fitted well with parameters corresponding to airplanes and helicopters. The
flying traces can be inverted by fitting the parameters at each station across the entire array.
They find that air-traffic events occupy >7% of the recorded wavefield. Johnson et al. (2019)
showed that wind interactions with obstacles above the surface produce ground velocities
larger than those expected to be generated by earthquakes with magnitudes M ≤ 1.5 about
6% of the day. Together with the detected RN, the developed methodologies allow us to
interpret the origins of ~49% of the recorded waveforms in SGB site. Motivated by the cross
frequency correlation study by Liu and Ben-Zion (2016), I derive three conditions which are
both sufficient and necessary by examining the expectation of waveform spectra and the
outer products with themselves and their conjugates (Chapter 4). As an application, we
perform this analysis to the detected RN, NRN and the Mixture. Both RN and the Mixture are
diffuse while NRN is not diffuse. This indicate ~68% of the waveform are diffuse.
In terms of the properties of earthquakes sources, I develop technique for detecting small
events not included in a seismic catalog using data of the dense seismic array (Meng & Ben-
Zion, 2018a and Chapter 1). Waveforms are first stacked without time shift in 9 non-
overlapping sub-arrays to increase the signal to noise ratio. The 9 envelope functions of the
stacked records are then multiplied with each other to suppress signals associated with
sources affecting only some of the 9 sub-arrays. Running an STA/LTA detection algorithm on
the product leads to 723 triggers in the test day. Using a local P-wave velocity model derived
for the surface layer from Betsy gunshot data, 5 s long waveforms of all sensors around each
STA/LTA trigger are beamformed for various incident directions. Of the 723 triggers, 220 are
found to have localized energy sources and 103 of these are confirmed as earthquakes by
verifying their observation at 4 or more stations of the regional seismic network. This
demonstrates the general validity of the method and allows processing further the validated
Doctoral Dissertation Haoran Meng

11
events using standard techniques. The number of validated events in the test day is > 5 times
larger than that in the standard catalog. Using these events as templates can lead to
additional detections of many more earthquakes. This indicates that a future network
configuration consisting of several spatially dense supersites can increase significantly the
detectibly of small events. Data from such a network will allow tracking active seismicity
structures with much greater detail and increasing considerably the resolution of
tomographic and other earthquake-based studies. To estimate of directivity, rupture area,
duration, and centroid velocity of small to moderate earthquakes, I develop semi-automated
estimation using second seismic moment (Chapter 5). I further use weighted stacked eGf to
improve the waveform fits to get better station coverage. The stress drop estimates are
reasonable and results of directivities are generally consistent with predictions for dynamic
rupture on a bimaterial interface and the imaged velocity contrasts in San Jacinto fault zone.
These developed methods are suitable for systematic processing of large seismic
waveform data sets to extract detailed properties of seismic waveforms and earthquake
sources. With these developed tools, it is possible for an geoscientist to understand a larger
fraction of the recorded seismic waveforms and add detailed information (i.e. rupture size,
directivity and stress drop) to earthquake catalogs.

Doctoral Dissertation Haoran Meng

12
1. Detection of small earthquakes with dense array data: Example from the San
Jacinto fault zone, Southern California
(Meng & Ben-Zion, 2018a)

1.1 Introduction
In the last decade it became a common practice to record seismic waveforms continuously.
The continuous recordings, combined with the increasing number of seismometers and
array deployments, provide important new opportunities for detecting small events, thereby
increasing the ability to resolve seismic processes and structural properties. One highly
successful technique is the matched-filter processing (or template matching) that relies on
cross-correlating seismograms of events detected by standard techniques with the
continuous waveforms (e.g., Gibbons & Ringdal 2006). Wave portions with high similarity to
the templates on multiple stations are assumed to reflect smaller events located very close
to the event templates. This method enables detection of events with Signal-to-Noise-Ratio
(SNR) below 1, and leads to 10-100 times more events than those listed in typical catalogs
(e.g., Meng et al. 2013; Shelly et al. 2016; Ross et al. 2017). However, a basic shortcoming of
this technique is that it allows only detection of events with hypocenters very close to events
already existing in the catalog.
Array techniques may be used to detect sources of radiation without pre-existing
templates by processing moving time windows of continuous waveforms. One example is the
matched-field processing technique (not to be confused with the matched-filter processing),
which combines beamforming and back-projection of moving time windows for detection
and location of energy sources. This technique was developed in underwater acoustics
(Turek & Kuperman, 1997) and was recently adapted to geophysical applications (Cros et al.
2011; Corciulo et al. 2012; Ben-Zion et al. 2015). A somewhat related method adopted from
medical imaging (Fink 2006) is reverse time migration of waveforms and identification of
focusing spots as locations of seismic sources (e.g., Larmat et al. 2008; Nakata & Beroza, 2015;
Inbal et al. 2015). An important shortcoming of these techniques is sensitivity to the used
velocity model. Also, in most applications the detected events are not subjected to validation
and may correspond to scatterers or other non genuine seismic sources.

Doctoral Dissertation Haoran Meng

13

Figure 1.1. A location map for the dense array providing the data used in this study (Ben-
Zion et al. 2015). Red circles in the main plot denote vertical geophones and yellow circles
mark Betsy gunshots. The background grey colors represent topography. The inset at the top
left provides a regional view with the dense array indicated by a red square and regional
seismic stations denoted by black triangles.

In this paper we present a simple new technique for detecting small events not included
in a seismic catalog using data of a dense seismic array. The technique is applied to seismic
waveforms recorded by a dense array with 1108 vertical geophones (10 Hz) centered on the
Clark branch of the San Jacinto fault zone (SJFZ) at the Sage Brush Flat site in the complex
trifurcation area southeast of Anza, California (Fig. 1.1). The array covered an area of about
600 m × 600 m with a core grid consisting of 20 rows perpendicular to and centered on the
Clark fault trace. Each row had 50 sensors at a nominal 10 m interstation spacing and the
nominal separation between rows was 30 m. The remaining 108 sensors were deployed as
extensions to some rows. The array recorded earthquake and noise data continuously at 500
Hz from 7 May 2014 to 13 June 2014. Additional controlled source signals were generated
Doctoral Dissertation Haoran Meng

14
by 33 surface Betsy gunshots (yellow circles in Fig. 1.1) during June 2 and 3, 2014 (Ben-Zion
et al. 2015).
The basic goal of the present paper is methodology development, rather than creating a
detailed catalog of events from all the recorded data. The developed methodology is
summarized in a flowchart diagram (Fig. 1.2) and is illustrated with data of one test day,
Julian day 146 (2014 May 26), not associated with Betsy gunshots or field work at the site.
The analysis includes a validation step involving searching for seismic signals in data
recorded by standard regional stations at times predicted by the detections. The various
steps of the analysis procedure are described and demonstrated in the next section. The
methodology includes slant-stacking and beamforming of waveforms using a detailed local
velocity model for the subsurface material derived in an Appendix. The detected/validated
events may be used as templates that can lead to further detections.

1.2 Method and Results
Stacking the waveforms recorded by the spatially dense array can reduce the random
noise. However, strong surface sources, large occasional instrumental glitches and other
incoherent sources cannot be suppressed by stacking the data across the entire array. The
presented methodology has several steps to resolve this problem and focus on signals likely
generated by sources below the surface (Fig. 1.2). We first divide the array of 1108 sensors
into 9 non-overlapping sub-arrays covering each an area of about 200 m × 200 m (Fig. 1.1).
The choice of sub-arrays is motivated by the fact that amplitudes of waveforms generated
(in a different day) by the Betsy gunshots decay by 2 orders of magnitude of more during a
propagation distance of 100 m (Fig. 1.S1). The raw waveforms within each sub-array are
stacked without time-shift to increase the SNR of coherent signals in the sub-array. The
stacking without time-shift suppresses signals from sources at the surface or very shallow
depth. Figs 1.3a and 1.3b provide examples of stacked waveforms within the 9 sub-arrays
for two representative 18 s time windows and the entire array. Black dashed boxes bracket
time intervals where stacked coherent signals are seen in all 9 sub-arrays, while dashed grey
boxes mark intervals with large amplitude signals in 4 or less sub-arrays. The 10th traces in
Figs 1.3a and 1.3b labeled “all” show that some of these latter signals are also present in a
direct stack of data across the entire array.
Doctoral Dissertation Haoran Meng

15

Figure 1.2. A flow diagram summarizing the main steps in the detection method.

Doctoral Dissertation Haoran Meng

16

Figure 1.3. (a) Example associated with a time window on the test day indicated on the
horizontal axes. The first 9 time-series displays plain stacked waveforms within each sub
array while the 10th trace is the stacked waveform across the entire array. The 11th trace is
the Product of the envelopes of the top 9 traces. The bottom trace displays the result of
running an STA/LTA on the Product function with a threshold for detection indicated by the
horizontal grey dashed line. The vertical dashed black box highlights a detection while the
vertical dashed grey box shows strong signals in only some sub arrays. (b) Same as (a) for
another example time window. See text for additional details.

Doctoral Dissertation Haoran Meng

17
Next we calculate the envelopes of the stacked waveforms within each of the 9 sub-arrays
using the square root of 𝑣 ( 𝑡 )
2
+ 𝐻 [𝑣 ( 𝑡 ) ]
2
, where 𝑣 ( 𝑡 ) is the original signal and 𝐻 [𝑣 ( 𝑡 ) ] is its
Hilbert transform. The resulting non-negative envelope functions are normalized for each
sub-array to the range 0 to 1 during the entire day. Multiplying the 9 envelopes functions
with each other produces a signal referred to as Product function in the range 0 to 1 that
suppresses significantly large amplitude signals present only in the stacks of some of the 9
sub-arrays. As seen in the 11th traces of Fig. 1.3, the Product functions have clear signals for
the time windows associated with the black dashed boxes and have very low amplitude in
all other time intervals. The bottom traces in Fig. 1.3 display results of running a standard
detector on the Product function consisting of the ratio (e.g., Allen, 1978) of a short term
moving average (STA) to a long term moving average (LTA). Since we target events that are
very small and relatively deep, we use 1 s for the short time window and 10 s for the long
one. The STA/LTA signals have significant peaks only in the time intervals associated with
the black dashed boxes.
Fig. 1.4a displays the Product function for the entire test day using a log scale. There are
clear overall daily variations of the signals, with lower level at night time (from 21:00 to
05:00 local time), and sharp spikes that indicate possible detections. Figs 1.4b and 1.4c
present zoomed-in views of one night and one day 10-minute time windows. Fig. 1.4d gives
the frequency-size statistics of the amplitudes of the Product function in Fig. 1.4a. The results
span 25 orders of magnitudes and resemble overall typical frequency-size statistics of
earthquakes, with approximate power law distribution over about 13 orders of magnitudes.
To decide on potential detections, we run an STA/LTA detector on the Product function with
the time intervals mentioned above. The results are shown in Fig. 1.5 along with a threshold
(horizontal dashed red line) selected arbitrarily at a level of 5 times the median of the
calculated STA/LTA values. This threshold choice leads to 723 potential detections (triggers)
of small events tested further below using a beamforming approach.

Doctoral Dissertation Haoran Meng

18

Figure 1.4. (a) The Product function in the range 0 to 1 for the test day with amplitudes
in a log scale. The dashed black boxes are night and day 10 min intervals, shown with zoomed
in view in (b) and (c), respectively. (d) Histogram of Product function amplitudes. Note that
the horizontal scale covers 25 orders of amplitudes.
Doctoral Dissertation Haoran Meng

19

Figure 1.5. Results of running an STA/LTA on the Product function. The red dashed line
denotes 5 times the daily median used as a threshold to accept STA/LTA triggers.

Figure 1.6. (a) A beamforming diagram showing the energy (Equations 3 and 4) of slant-
stacked waveform within 5 s around the STA/LTA trigger of Fig. 3a. The radii on the polar
plot denote incident angles and the numbers on the outer circle indicate back azimuth angles.
The white cross marks the maximum energy. (b) Slant-stack waveforms corresponding to
the incident angle and back azimuth of the white cross. (c) and (d) Same as (a) and (b) for
Doctoral Dissertation Haoran Meng

20
the STA/LTA trigger of Fig. 3b. Note that the incident angle is defined on the top of bottom
medium (Fig. 1.A1b). The 1-D model used for beamforming is displayed in Fig. 1.A3b.

Figure 1.7. (a) A histogram of all 723 STA/LTA triggers for the test day. The ratio value
on the horizontal axis indicates the quality of the detection in terms of the degree of
localization in the beamforming diagram (Fig. 1.6). Blue bars denote high quality accepted
STA/LTA triggers. (b) A rose diagram of the back azimuths of the 220 high quality triggers.
(c) The incident angle and back azimuth of high quality triggers with ratio values denoted by
color.

The beamforming is performed together with slant stacking where, for each of the 723
potential detections, 5 s waveforms recorded by the entire array are shifted, aligned and
stacked for different possible incident directions. The relatively short duration of the
waveforms is appropriate for small target events, and each waveform starts 1 s before and
ends 4 s after a trigger. The analysis employs a detailed 1D-piecewise P-wave velocity model
for the shallow crust consisting of a surface layer of variable width over a half space. The
velocity model is derived from data associated with the Betsy gunshots and presented in
Appendix 1.A. The beamforming assumes a plane P-wave with incident angle i and back
azimuth baz at the top of the bottom medium (Fig. 1.A1b and 1.A3a), and uses second root
slant stacking to suppress noise and enhance the coherent signals (Rost & Thomas, 2002).
Specifically, we calculate
Doctoral Dissertation Haoran Meng

21

𝑢 ′
( 𝑡 ; 𝑡 , 𝑏𝑎𝑧 ) =
1
𝑁 ∑ sgn( 𝑣 𝑗 )∙ √|𝑣 𝑗 ( 𝑡 − 𝑡 𝑗 ( 𝑖 , 𝑏𝑎𝑧 ) ) |
𝑁 𝑗 =1
(1.1)

𝑢 ( 𝑡 ; 𝑡 , 𝑏𝑎𝑧 ) = sgn( 𝑢 ′)∙ |𝑢 ′
( 𝑡 ; 𝑡 , 𝑏𝑎𝑧 ) |
2
(1.2)
where N is the number of the array stations and 𝑣 𝑗 ( 𝑡 )represents the seismogram at
station j in a 5 s time window with 1 s before the trigger and 4 s after. The S wave energy is
suppressed in the stacking since we use the local P-wave velocity model and residuals at each
station (Appendix) to align possible P arrivals. The total energy recorded by the array for
each 5 s time window associated with each trigger is calculated by integrating the square of
the signal obtained from the slant stack procedure (Equation (1.2))
𝐸 ( 𝑖 , 𝑏𝑎𝑧 ) = ∫ 𝑢 2
( 𝑡 ; 𝑡 , 𝑏𝑎𝑧 ) 𝑑𝑡
𝑡 𝑘 +𝑡 2
𝑡 𝑘 −𝑡 1
(1.3)
where 𝑡 𝑘 is the STA/LTA trigger time, t
1
=1 s and t
2
= 4 s. The calculations are repeated
in a grid search over the incident angle and back azimuth in the lower half space with a 2
spacing. The corresponding energy values calculated with Equation (1.3) are normalized
using
𝑒 ( 𝑖 , 𝑏𝑎𝑧 ) =
𝐸 ( 𝑖 , 𝑏𝑎𝑧 )
max ( 𝐸 ( 𝑖 , 𝑏𝑎𝑧 ) )
(1.4)
Figs 1.6a and 1.6c present beamforming diagrams for the two STA/LTA triggers in Fig. 1.3,
where each point represents the integrated power for given incidence angle and back
azimuth, and the white crosses denote the strongest energy. The associated second-root-
slant-stacked waveforms have considerably less noise than the corresponding linear slant
stacks (Figs 1.6b and 1.6d). To evaluate the quality of a possible detection, we integrate the
beamforming results for each trigger for directions with high energy on a lower hemisphere
with unit radius using
𝑟𝑎𝑡𝑖𝑜 =
1
2𝜋 ∬ 𝑑𝑏𝑎𝑧 ∙ 𝑑𝑖
𝑒 ≥ 0.8
∙ sin ( 𝑖 )𝑒 ( 𝑖 , 𝑏𝑎𝑧 ) . (1.5)
Doctoral Dissertation Haoran Meng

22
The integration is normalized by 2𝜋 so the ratios span the range 0 to 1. Low ratios
associated with concentrated power distribution are taken to represent high-quality
detections. Fig. 1.7a presents the histogram of the calculated ratios for all 723 STA/LTA
triggers. Defining for example high quality results to have ratios less than 0.1 leads to 220
detections. These detections include all 18 events during the test day found with standard
techniques using the Southern California network, ANZA network and a dense local PASSCAL
deployment (YN) in and around the SJFZ (Vernon et al., 2014). Using a less strict criterion
associated with 0.25 leads to 339 detections. Figs 1.7b and c show, respectively, the azimuth
distribution and incident directions of the detections associated with ratios ≤ 0.1.
As a validation step of the detection analysis, we check the waveforms at stations of the
regional Southern California and ANZA networks within ±40 s from the predicted arrival
time of each of the 220 high-quality detections. The predicted arrival time are calculated
using the velocity model of Fig. 1.A3 and the waveforms at regional stations are band-pass
filtered from 2 to 15 Hz. Figs 1.8a and 8b display clear seismic signals at regional stations
(red triangles at inset maps and name on left) close in time to predicted arrivals associated
with the two STA/LTA triggers shown in Fig. 1.3 (dashed vertical red line). Of the 220 high
quality detections, 103 are confirmed as seismic events by finding corresponding phases at
4 or more seismic stations. The number of confirmed events is 5 times larger than the events
in the ANZA and SCSN catalogs during the test day; many more of the 723 STA/LTA triggers
in Fig. 1.5 may also be genuine seismic events.
To locate approximately the validated events, we manually pick the P or S arrivals and
adopt the method of Geiger (1912), associated with ray tracing in a 1D model obtained by
averaging the 3D velocity model of Lee et al. (2014) for the region shown in the top left inset
of Fig. 1.1. To improve the convergence and stability of the algorithm, we grid search the
focal depth with a 0.5 km step and derive tentative epicenter and origin time for each focal
depth. The initial epicenter is set to the horizontal location of the station with the earliest P
or S arrival and the initial origin time is set to the earliest P or S arrival. The location
generating the minimum L2 norm misfit in the grid search is taken to be the source position.
Fig. 1.9 presents the obtained locations of the validated events. The location uncertainties
are estimated to be on the average ~4 km horizontally and ~4 km vertically considering the
misfit curves and reasonable velocity variations from the adopted average 1D model. Fig. 1.9
Doctoral Dissertation Haoran Meng

23
also shows estimated magnitudes of the validated detections, based on relative amplitudes
to small reference events with available catalog magnitudes and focal mechanisms
(Appendix 1.B). Most detected events are very small with magnitudes around or below zero,
but the results include also events with magnitudes up to 2.5 that are not detected (for some
reason) by the standard techniques producing the regional catalog.

Figure 1.8. (a) Seismic records at several stations of the Southern California and Anza
networks around STA/LTA trigger time (red dashed line) in Fig. 1.3a. The waveforms are
self-normalized and band-passed from 2 to 15 Hz. Station names and components are
marked on the left. The inset in (b) shows with red triangles stations with clear phases. (c)
and (d) Same as (a) and (b) for waveforms in a time window around STA/LTA trigger time
in Fig. 1.3b.
Doctoral Dissertation Haoran Meng

24

Figure 1.9. Approximate locations (color) and magnitudes (circle size) of the 103
validated events. The bold circles mark the locations of the two example events E1 and E2 in
Figs 1.3, 1.6 and 1.8. See text for additional details.

1.3 Discussion and conclusions
We developed a multi-step procedure for detecting small earthquakes with continuous
waveforms recorded by a dense seismic array, and illustrated the technique with data
recorded by 1108 sensors in a tight configuration around the Clark branch of the SJFZ (Fig.
1.1). The different analysis steps (Fig. 1.2) are designed to suppress local sources affecting
only sub-arrays (likely associated with surface and/or instrumental noise), and then
detecting, validating and locating sources affecting the entire array. The implementation
involves stacking data of 9 sub-arrays without migration, multiplying the envelopes of the
stacked sub-array data, and performing the remaining steps on the Product function. Sources
detected on the Product function have coherent energy across the entire array so are likely
to be associated with events below the surface. The validation and location steps of the
Doctoral Dissertation Haoran Meng

25
method demonstrate that at least a subset of these detections (103 during the test day) is
indeed associated with small earthquakes in the region.
The choice of using 9 sub-arrays in this study stems from the observed significant
amplitude decay of Betsy gunshots over distance of 100 m (Fig. 1.S1). Using fewer sub-arrays
would lead to reduced suppression of surface sources, many more STA/LTA triggers of local
sources and considerably more computation in the following beamforming step. On the other
hand, using more sub-arrays would result in fewer stations in each sub-array and reduced
suppression of random noise. Appendix 1.C presents a derivation of an optimal number of
sub-arrays for detecting weak signals arriving with large incident angle (corresponding to
deep events with general locations not directly below the array). The derivation assumes
stacking in sub-arrays without time shifts and subsequent analysis of the Product function
as in our procedure. The choice of 3×3 sub-arrays made in this study helps to maximize the
array sensitivity to weak signals with large incident angles and facilitates the detection of
more events. We recall that the direct stacking of data in each sub-array in the early analysis
procedure is motivated by our focus on deep events and helps to suppress surface and
shallow sources. Efforts to detect shallow small events should incorporate migration in each
sub-array before stacking to enhance the Product function.
The detection procedure includes beamforming and slant staking using a 1D-piecewise
local P-wave velocity model for the subsurface material (Fig. 1.2). The delay times for each
pair of incident angle and back azimuth are calculated only once using equation 1.A23 and
then used for slant staking and beamforming for all triggers, making the process
computationally efficient. In our case, computing the beamforming on a CPU node with 2 GB
memory requires ~3 minutes for one trigger and under 40 hours for all triggers. The model
consists of a layer with a constant velocity gradient over a half space, and is derived in
Appendix 1.A from arrivals of Betsy gunshots and analytical expressions with 4 parameters.
Although the model is simple, it accounts for topographic variations across the array and it
explains the data well (Figs 1.A3 and 1.A4). Residuals from model predictions provide site
corrections at each station that are incorporated in the beamforming and slant staking. The
velocity variations across the array are consistent with features observed in previous
imaging studies based on the array data (Ben-Zion et al. 2015; Roux et al. 2016; Hillers et al.
2016). The rapid velocity variations with depth in the very shallow crust (P velocities
Doctoral Dissertation Haoran Meng

26
increasing from ~550 ms/ to ~1250 m/s in the top ~100 m and then jumping to ~2900 m/s)
have important implications for properties of the subsurface material (e.g., Sleep 2011;
Kurzon et al. 2014) and can be useful for future imaging and detection studies in the area.
The analysis procedure reduces the data redundancy of the dense array and zooms in
progressively on small earthquakes (Figs 1.3-1.7) by calculating STL/LTA on the Product
function representing coherent signals across the array and performing beamforming on
time windows around the STL/LTA triggers. The method generates 723 triggers during the
test day, of which 220 are shown to be associated with localized energy sources below the
surface (Fig. 1.7). Of these, 103 are validated by showing that they produce earthquake
signals in 4 or more regional stations. Figs 1.3, 1.6 and 1.8 illustrate with two examples
STL/LTA triggers how detections are validated as genuine seismic events. The remaining
117 localized energy sources are also likely to be associated with seismic events, too small
to be recorded by ≥ 4 standard seismic stations. The same likely holds for many of the ~500
additional triggers in the test day. The slant-stacking performed in this study is based on the
local P-wave velocity model and reduces signals associated with S waves. Performing in
parallel beamforming and slant-stacking associated with an S wave velocity model can lead
to additional detections and validations. Using the detected/validated events as templates
will likely increase the number of detections by a factor of 20 or more (e.g., Shelly et al. 2016;
Ross et al. 2017).
The approximate locations of the validated events (Fig. 1.9) are calculated using an
average 1D velocity model based on the 3D tomographic results of Lee et al. (2014) for
Southern California. The event depths range from ~3 km to ~25 km. The back azimuths of
the events are generally consistent with the incident directions in Figs 1.7b and 1.7c with
some discrepancies. For instance, the locations of the example events in Figs 1.6 a and c are
marked in Fig. 1.9. The back azimuth of the first event is within 2° to that estimated by the
beamforming (Fig. 1.6a), while for the second event the discrepancy is 15° (Fig. 1.6 c). The
discrepancies reflect the significant variations of seismic velocities around the SJFZ (e.g.,
Allam & Ben-Zion 2012; Zigone et al. 2015). The uncertainties of event locations are relatively
large because of the employed average 1D velocity model and the fact that many events have
only a limited number of picks. The magnitudes of the validated events are estimated by
calculating relative amplitudes to reference small catalog events near the hypocenters of the
Doctoral Dissertation Haoran Meng

27
detected events, with attention to likely radiation pattern effects on amplitudes (Appendix
1.B).
The methodology and results of this paper can be improved in several ways. The locations
can be refined by using a 3D velocity model and attempting to obtain additional high quality
picks. The study can be extended by analyzing more data and using both P and S waves in
the detection stage. Another useful future research direction would be an effort to detect
very shallow events. Focusing on shallow events can benefit from using a shorter STA
interval (e.g. 0.1 s or less) and incorporating migration in the initial stacking of the sub-
arrays. These additional improvements notwithstanding, the study demonstrates that dense
arrays of the type used in this work can form “supersites” that allow the creation of many
new validated templates (> 100 per day in the study area). Combining this additional set of
small events not included in regular catalogs with standard matched-filter technique can
increase the detection of small events further. We note that the magnitudes and source-
receiver distances of the validated detections significantly outperform the theoretical limits
to detection of Kwiatek & Ben-Zion (2016) based on simulations of wave propagation from
various sources and assumed SNR > 1 at individual sensors, because the stacking enhances
considerably the SNR of the array. A future network configuration consisting of several
spatially dense supersites can increase significantly the detectibly of small events. Data from
such a network will allow tracking active seismicity structures with much greater detail and
increasing considerably the resolution of tomographic and other earthquake-based studies.

Appendix 1.A. A 1D-piecewise P-wave velocity model for the very shallow crust
1.A1. Ray tracing in a two-media model
The P-wave structure beneath the dense deployment is simplified as a surface layer over
a half space with a constant velocity-gradient in the top layer and a constant velocity in
the bottom medium (Fig. 1.A1). Down-going waves in the top layer curve back to the surface
and the constant velocity in the bottom solid is sufficient to account for the incident plane
wave assumption of the beamforming approach. The remaining two parameters are the
velocity at the bottom of the surface layer and the depth of the interface. The velocity
Doctoral Dissertation Haoran Meng

28
v
st
at a station is a function of z, so the topography (grey lines in Fig. 1.A1a) is considered in
the ray tracing.

Figure 1.A1. A framework for deriving from Betsy gunshot data a shallow velocity model
using a surface layer with topography denoted by grey lines over a half-space. (a) A diagram
for ray tracing of a regular body wave (red line) and arrival associated with a head wave at
the bottom of the layer (purple line). The star and triangles denote the source and two
example stations, respectively. (b) A diagram for arrivals at two stations from a plane wave
(bold black line) at the bottom layer. The red and purple curves denote rays at reference and
target stations, respectively. See text for additional explanations.

Doctoral Dissertation Haoran Meng

29
1.A1.1 Down-going wave in a solid with constant velocity-gradient
We consider a model in which the P-wave velocity v increases linearly with depth (Fig.
1.A1a). By referring to the epicenter distance and travel time of up-going wave derived by
Slawinski and Slawinski (1999), the epicenter distance of down-going wave in the top layer
is
𝑋 =
1
𝑏𝑝
(√1 − ( 𝑝 𝑣 src
)
2
+ √1 − ( 𝑝 𝑣 st
)
2
), (1.A1)
where 𝑝 is the constant ray parameter, 𝑏 is the constant velocity gradient in the top layer,
and, 𝑣 src
and𝑣 st
are the velocities at source and station positions. Solving equation (1.A1)
for 𝑝 and ignoring non-physical solutions gives

𝑝 =
2𝑏𝑋
√[( 𝑏𝑋 )
2
+ 𝑣 src
2
+ 𝑣 st
2
]
2
− ( 2𝑣 src
𝑣 st
)
2
.
(1.A2)
Thus, the travel time for down going wave is
𝑡 𝑝 =
1
𝑏 |ln [
1 − √1 − ( 𝑝 𝑣 src
)
2
𝑝 𝑣 src
]| +
1
𝑏 |ln [
1 − √1 − ( 𝑝 𝑣 st
)
2
𝑝 𝑣 st
]|. (1.A3)

1.A1.2 Head wave
As the epicenter distance increases, the ray can refract along the interface between the
media before turning up and reaching the surface (Fig. 1.A1a). In this case, the first arrival is
a head wave with ray parameter 𝑝 = 1/𝑤 0
, where 𝑤 0
is the constant velocity in the bottom
solid. For the down-going segment of the head wave, the horizontal distance (𝑥 1
in Fig. 1.A1a)
is
𝑥 1
=
1
𝑏𝑝
(√1 − ( 𝑝 𝑣 src
)
2
− √1 − ( 𝑝 𝑣 0
)
2
), (1.A4)
where 𝑣 0
is the velocity at the bottom of the top layer. The corresponding travel time is
𝑡 1
=
1
𝑏 |ln [
𝑣 0
𝑣 src
1 − √1 − ( 𝑝 𝑣 src
)
2
1 − √1 − ( 𝑝 𝑣 0
)
2
]|. (1.A5)
For the up-going segment, the horizontal distance (𝑥 3
in Fig. 1.A1a) and corresponding travel
time are
𝑥 3
=
1
𝑏𝑝
(√1 − ( 𝑝 𝑣 st
)
2
− √1 − ( 𝑝 𝑣 0
)
2
), (1.A6)
Doctoral Dissertation Haoran Meng

30
and
𝑡 3
=
1
𝑏 |ln [
𝑣 0
𝑣 src
1 − √1 − ( 𝑝 𝑣 st
)
2
1 − √1 − ( 𝑝 𝑣 0
)
2
]|. (1.A7)
For the ray segment at the interface between the media and epicenter distance , the
horizontal distance (x 2 in Fig. 1.A1a) and corresponding travel time are
𝑥 2
= 𝑋 − 𝑥 1
− 𝑥 3
, (1.A8)
and
𝑡 2
= 𝑝𝑥
2
. (1.A9)
If 𝑥 2
> 0, the ray is a head wave; otherwise, the ray just propagates in the top layer. The total
travel time for the head wave is
𝑡 ℎ
= 𝑡 1
+ 𝑡 2
+ 𝑡 3
. (1.A10)

Figure 1.A2. Waveforms with automatic gain control, generated by Betsy gunshot No. 1
(Fig. 1.S1) and recorded at sensors along row 7 of the array. First arrivals are marked by
short vertical line segments.
Doctoral Dissertation Haoran Meng

31

Figure 1.A3. (a) A derived P velocity model along cross section AA’ of Fig. 1. The star and
triangles represent a Betsy gunshot and geophones. (b) A depth varying P velocities, starting
from surface points with different elevations as in (a), inverted from the Betsy gunshots data.
(c) Observed and predicted travel times of all Betsy gunshots data.

1.A1.3 Travel time from an incident plane wave at the bottom of the surface layer
To apply beamforming approach to data, we need travel times for each station from an
incident plane wave. In the Cartesian coordinate system defined in Fig. 1.A1b, the location of
a reference station is
𝑟 st
= [
𝑥 𝑦 𝑧 ], (1.A11)
while the location of a target station is
𝑟 st
′
= [𝑥 ′
𝑦 ′
𝑧 ′ ]. (1.A12)
Doctoral Dissertation Haoran Meng

32
For a plane wave with back azimuth and incident angle 𝑖 , the ray parameter is 𝑝 =
sin𝑖 /𝑤 0
and the incident direction vector is
𝑘 ̂
= [
−cos 𝑏𝑎𝑧 ∙ sin𝑖 −sin𝑏𝑎𝑧 ∙ sin𝑖 −cos 𝑖 ]. (1.A13)
The horizontal offsets between the intersections on the interface and the surface stations
(𝑥 st
and 𝑥 st
’ in Fig. 1.A1b) are
𝑥 st
=
1
𝑏𝑝
(√1 − ( 𝑝 𝑣 st
)
2
− √1 − ( 𝑝 𝑣 0
)
2
), (1.A14)
and
𝑥 st
′
=
1
𝑏𝑝
(√1 − ( 𝑝 𝑣 st
′ )
2
− √1 − ( 𝑝 𝑣 0
)
2
). (1.A15)
The corresponding travel times in the top layer are
𝑡 =
1
𝑏 |ln [
𝑣 0
𝑣 st
1 − √1 − ( 𝑝 𝑣 st
)
2
1 − √1 − ( 𝑝 𝑣 0
)
2
]|, (1.A16)
and
𝑡 ′
=
1
𝑏 |ln [
𝑣 0
𝑣 st
′
1 − √1 − ( 𝑝 𝑣 st
′ )
2
1 − √1 − ( 𝑝 𝑣 0
)
2
]|. (1.A17)
For nearly vertical incident plane wave at the bottom layer, the ray parameter 𝑝 → 0. The
limits of equations (1.A16) and (1.A17) for this case are
lim
𝑝 →0
𝑡 =
1
𝑏 |ln (
𝑣 st
𝑣 0
)|, (1.A18)
and
lim
𝑝 →0
𝑡 ′ =
1
𝑏 |ln (
𝑣 st
′
𝑣 0
)|. (1.A19)
The intersections for both rays at the interface are
𝑟 itsc
= [𝑥 + 𝑥 𝑠𝑡
cos 𝑏𝑎𝑧 𝑦 + 𝑥 𝑠𝑡
sin𝑏𝑎𝑧 𝑧 0
], (1.A20)
and
𝑟 itsc
′
= [𝑥 ′
+ 𝑥 𝑠𝑡
′
cos 𝑏𝑎𝑧 𝑦 ′
+ 𝑥 𝑠𝑡
′
sin𝑏𝑎𝑧 𝑧 0
]. (1.A21)
For the target station, the travel time for the ray segment in the bottom layer is

𝑡 𝑑 =
( 𝑟 itsc
′
− 𝑟 itsc
)∙ 𝑘 ̂
𝑤 0
.
(1.A22)
Thus, the delay time between arrivals at the target and reference stations is
Doctoral Dissertation Haoran Meng

33
∆𝑡 = 𝑡 ′
+ 𝑡 𝑑 − 𝑡 . (1.A23)

1.A2 Inversion of data for a 1D model
For a given model 𝑚 0
= [𝑏 𝑣 0
𝑤 0
𝑧 0], the travel time of the first P arrivals 𝑇 = 𝑇 ( 𝑚 0
)
at all stations are predicted by equations (1.A3) and (1.A10). Starting with a reasonable
initial model 𝑚 0
and observed arrivals𝑇 0
, the final model 𝑚 is inverted by improving data fit
iteratively using
𝑚 = 𝑚 0
+ ( 𝐺 T
𝐺 )
−1
𝐺 T
( 𝑇 0
− 𝑇 ( 𝑚 0
) ) (1.A24)
𝐺 =
𝜕𝑇
𝜕𝑚
|
𝑚 =𝑚 0
. (1.A25)
The partial derivatives can be easily calculated numerically. The mean residuals at individual
stations are used as station correction to mitigate in the beamforming effects not accounted
for in the inverted 1D-picewise model.

1.A3 Imaging the structure below the array
We hand pick first P arrivals generated by the 33 Betsy gunshots (yellow circles in Fig. 1.1)
at all 1108 sensors. To make more accurate pick, we use an automatic gain controlled (AGC)
algorithm to enhance the first arrivals of each gunshot. Fig. 1.A2 displays the AGC waveforms
and the corresponding arrival picks. In total, 22406 high quality picks are collected for 28
Betsy gunshots (3 of the other gunshots are at repetitive location and 2 are bad). The
inversion is highly overdetermined because the model has only 4 parameters. The final P
velocity model and data fits are displayed in Fig. 1.A3. The obtained model can be converted
to a local S velocity model assuming a distribution of Poisson ratios. The velocity gradient in
the top layer (5.1 s
-1
) is high and the P-wave velocity jumps from 1120 m/s to 2689 m/s
across the interface at elevation of 1373 m, which is about 110 m below the local topographic
high in the study area (Fig. 1.A3b). The 1D velocity model generally predicts well the travel
times (Fig. 1.A3c). The data scattering can mostly be predicted by the topography.
Fig. 1.A4a shows the average residuals of travel time at each station. Positive values
represent arrivals later than predicted (slow local anomalies) and negative values represent
the opposite situation. There are two major low velocity anomalies. The one around the left
fault trace is associated with a local sedimentary basin, while the one on the NE side of the
Doctoral Dissertation Haoran Meng

34
right trace is associated with fault damage zone including a trapping structure (Ben-Zion et
al. 2015). Since the topography is considered in the model, 3D wave propagation effects
contribute most of the remaining residuals. The mean residuals of Fig. 1.A4a are used as
station corrections in the slant stacking done in the paper. Fig. A4b displays the azimuthal
distribution of residuals. Assuming P-wave azimuthal anisotropy, we fit the residuals with
station correction removed (Fig. 1.A4c) to cosine function res= A×cos(2az+f) . Although
station corrections are removed, we still find weak but possible P-wave azimuthal
anisotropy with fast axis parallel to the strike of SJFZ. This is consistent with noise-based
imaging results of Hillers et al. (2016) based on the dense deployment data.

Doctoral Dissertation Haoran Meng

35

Figure 1.A4. (a) Average residuals between observed and predicted first arrivals from
Betsy gunshots data at sensors of the dense array. The background grey colors represent
topography. (b) Residuals versus azimuth with red line showing the best fit of 𝑟𝑒𝑠 = A ∙
cos ( 𝑎𝑧 + ∅) with az denoting azimuth. (c) Residuals with station correction removed.

Doctoral Dissertation Haoran Meng

36

Figure 1.B1. A histogram of magnitudes for the 103 validated events.

Appendix 1.B. Magnitude estimation
The magnitude estimation is motivated by the empirical Green’s function (eGf) approach
to account for propagation and site effects. We first find reference events within a 5 km
radius sphere of the hypocenter of the target event using the focal mechanism catalog of Yang
et al. (2012) extended to later events. The magnitudes of reference events are smaller than
3 so the finite source effects are negligible. We then measure the vector sum of the max
velocity amplitude 𝐴 𝑖 at station i for the target event and 𝐴 0
𝑖 for the reference event using the
seismic velocity time series (HH channels) filtered at 2-15 Hz to increase the signal to noise
ratio. The magnitude M
l
is estimated to be
𝑀 𝑙 = median(log (
𝐴 𝑖 /𝑅 𝑖 𝐴 0
𝑖 /𝑅 0
)) + 𝑀 𝑙 0
, (1.B1)
where 𝑀 𝑙 0
is the magnitude of the reference (eGf) event and 𝑅 0
, 𝑅 𝑖 are the radiation
patterns of P or S phases of the reference and target events at station i, respectively. We use
median instead of mean in equation 1.B1 to reduce possible effects of outliers. Since it is not
possible to determine the focal mechanism of the target event with the limited available
measurements, we estimate the mechanism from a probability distribution calculated as the
joint distribution of the focal mechanisms and their uncertainties of all the reference events.
Doctoral Dissertation Haoran Meng

37
Specifically, we draw 1,000 possible mechanisms of the target event from the distribution,
calculate for each the radiation pattern R
i
at station i and estimate the correspondingM
l
for
each reference event. The final magnitude of the target and its uncertainty are set to be the
mean and standard deviation of the 1,000 estimates of M
l
. The estimated magnitudes of the
validated detected events for Julian day 146 span the range －1.5 to 2.5 (Fig. 1.B1).

Appendix 1.C. Sub-array configuration
Here we analyze the optimal selection of sub-arrays with the procedure used on this work
(Fig. 1.2), with the goal of increasing the sensitivity to weak coherent signals with large
incident angle from sources below the surface. For simplicity, we assume that N stations
are evenly installed in a dense deployment with square shape on a flat surface. The signal
from the source is approximated by a plane wave with incident angle i , back azimuth baz ,
and duration longer than the maximum moveout across the array. The average SNR of
waveforms at the array stations is r
0
and the entire array is divided into m´m sub-arrays.
Assuming random noise, the SNR of the stacked waveform for each sub-array is

𝑟 𝑠𝑢𝑏 = 𝑐 ( 𝑖 , 𝑏𝑎𝑧 )
√𝑁 𝑚 𝑟 0
,
(1.C1)
where 𝑐 ( 𝑖 , 𝑏𝑎𝑧 ) is a stacking coefficient related to the aperture and density of the array,
signal duration and incident direction. We note that 𝑚 /√𝑁 < 𝑐 ( 𝑖 , 𝑏𝑎𝑧 )≤ 1 for stacking
without time shifts and 𝑐 ( 𝑖 , 𝑏𝑎𝑧 ) =1 when 𝑖 = 0° and the waveforms are perfectly aligned. In
general, 𝑐 ( 𝑖 , 𝑏𝑎𝑧 ) decreases with increasing 𝑖 or array aperture, and 𝑐 ( 𝑖 , 𝑏𝑎𝑧 ) increases with
increasing array density or signal duration. By multiplying 𝑚 × 𝑚 envelops of the stacked
waveform, the SNR of the Product function is

𝑟 𝑃 ( 𝑚 ) ≤ [𝑐 ( 𝑖 , 𝑏𝑎𝑧 )
√𝑁 𝑚 𝑟 0
]
𝑚 2
.
(1.C2)
Using 𝐼 ( 𝑚 ) to represent the right-hand side of equation (1.C2), and a > 1 to
represent 𝑐 ( 𝑖 , 𝑏𝑎𝑧 ) √𝑁 𝑟 0
, the maximum of 𝐼 ( 𝑚 ) is found from the derivative

∂𝐼 ( 𝑚 )
∂𝑚 = 𝑚 [
𝑎 𝑚 ]
𝑚 2
[2ln (
𝑎 𝑚 ) − 1].
(1.C3)
Doctoral Dissertation Haoran Meng

38
Setting (1.C3) to zero, the optimal choice of 𝑚 is a/√ 𝑒 . For our dataset with 𝑁 = 1108 , 𝑟 0

is slightly larger than 1 with 𝑚 = 3 (9 sub-arrays) when 𝑐 ( 𝑖 , 𝑏𝑎𝑧 ) = 0.15. Therefore, the
choice 𝑚 = 3 adopted in this paper helps maximizing the sensitivity of the array to weak
signals with larger incident angles.

Figure 1.S1. Peak ground velocity (PGV) generated by a Betsy gunshot normalized across
the array and color-coded using a log scale. The background grey colors represent
topography.

Doctoral Dissertation Haoran Meng

39
2. Characteristics of airplanes and helicopters recorded by a dense seismic array
near Anza California
(Meng & Ben-Zion 2018b)

2.1 Introduction
The growing quantity and quality of continuous waveforms recorded by seismic networks
and temporary arrays at various scales, coupled with improved sensors and processing
techniques, provide increasing opportunities for detection and analyses of weak sources of
ground motion. Detection and analyses of events with signal-to-noise (SNR) ratios below 1
is now commonly done in studies of earthquakes, tremor and induced seismicity around the
world (e.g., Shelly et al., 2007; Zigone et al., 2012; Ross et al., 2017; Chen et al., 2017). The
association of weak signals with SNR ≤ 1 with genuine seismic sources below the surface
(earthquakes, tremor) requires detailed understanding of the differences between such
events and sources of energy at and above the surface. These include various types of cultural
noise at the surface (e.g., cars and trains), interactions of the wind with surface structures
and sources in the atmosphere (e.g., storms and air-traffic).
In addition to developing reliable catalogs of weak sub-surface seismic events, the ability
to detect and classify different sources that produce small episodic ground motion is
important for tomographic imaging based on the ambient seismic noise. These studies,
which are now done routinely in many environments (e.g., Campillo et al., 2011; Nakata et
al., 2015), require a diffuse wavefield (e.g., Weaver, 1982; Shapiro et al., 2000) and removal
from observed waveforms of all other sources (e.g., Bensen et al., 2007). More generally,
detailed characterization of different types of common sources of ground motion contributes
to a more complete understanding of the observed seismic waveforms.
Doctoral Dissertation Haoran Meng

40

Figure 2.1. (a) A regional map showing the dense array on the San Jacinto fault zone (red
square), regional seismic stations (black triangles) and station JFS4 (red triangle). (b) A
conceptual diagram of an airplane flying over a station (black triangle). (c) Zoom in view of
the dense deployment at the SGB site. The red circles denote vertical component geophones,
the yellow circle marks a geophone used to present example data, the yellow triangle
Doctoral Dissertation Haoran Meng

41
denotes strong ground motion station SGBF0 and the green triangle marks 148 m deep
borehole station B946. Data recorded by the instruments along profile A-A' are shown in Fig.
6. The background grey colors represent topography and the lines show surface traces of the
Clark fault at the site.

Figure 2.2. (a) A raw waveform recorded by the marked geophone in Fig. 1c of the dense
deployment during Julian day 146 (2014), 17:23:20 to 18:00:00 GMT time. (b) A
corresponding spectrogram with solid and dashed white boxes marking air-traffic events.

In the present work we quantify characteristics of airplane and helicopter signals
recorded by a dense seismic array with 1108 vertical ZLand geophones (10 Hz) centered on
the Clark branch of the San Jacinto fault zone (SJFZ) at the Sage Brush Flat (SGB) site
southeast of Anza, California (red box Fig. 2.1a). The array recorded continuously at 500 Hz
from 7 May 2014 to 13 June 2014, and covered an area of about 600 m × 600 m (Fig. 2.1c)
with nominal sensor spacing of 10 m normal to the fault and 30 m along-strike (Ben-Zion et
Doctoral Dissertation Haoran Meng

42
al., 2015). The array data has been used for detailed imaging of the fault zone structure based
on the ambient seismic noise (Roux et al., 2016; Hillers et al., 2016) and detection of
microearthquakes (Meng & Ben-Zion, 2018a; Gradon et al., 2019). The focus in the present
work on air-traffic events is complementary to these studies. In addition to data of the dense
ZLand array, we examine waveforms recorded by a borehole senor ~150 m deep at station
B946 (green triangle in Fig. 2.1c), and use data of a strong ground motion station SGBF0
(yellow triangle in Fig. 2.1c) to measure durations of cataloged earthquakes in year 2014. As
a control test, we also analyze data of a broadband station JFS4 at the Jackass Flat site (red
triangle in Fig. 2.1a) with known record of a helicopter flying nearby. The borehole sensor at
B946 records continuously at 100 Hz and the sensors at SGBF0 and JFS4 record continuously
at 200 Hz.
In the next section we present and analyze example observations of airplane and
helicopter events that generate wavetrains with typical duration of several minutes. The
appearance of these events in waveforms resembles in casual analysis tremor or a sequence
of small earthquakes, but they have clear Doppler effects that can be easily recognized with
time-frequency analysis and modeled with basic analytical expressions. The air-traffic
events cover in the study area > 7% of the total time in a day, so they make a significant
contribution to the observed ground motion. In large cities with highly active airports such
as Los Angles, the air-traffic events can cover > 75% of the day. Two appendixes provide
material on estimating the flying traces from seismic array data and a scaling relation
between earthquake durations and magnitudes. The results are discussed and summarized
in section 3 of this chapter.
Doctoral Dissertation Haoran Meng

43

Figure 2.3. (a) Zoom in waveform for an airplane event in the white solid box in Fig. 2.2
recorded by the example station in Fig. 2.1. (b) A corresponding spectrogram. The grey boxes
denote 20 s window of data shown by Fig. 6a. (c) Time and frequency data (black crosses)
picked in the spectrogram in (b) and line fit to the data based on equation 2.1 (red curve).
The inferred parameters for the airplane are indicated in the box.

Doctoral Dissertation Haoran Meng

44
2.2 Analysis
2.2.1 Observations
We start the analysis by examining raw waveforms recorded by the dense array and
corresponding spectrograms. Fig. 2.2a shows example results from a 2200 s time window in
Julian day 146, 2014 recorded by the sensor marked in Fig. 2.1c. The raw waveforms in this
and other examples include frequent tremor-like wavetrains with significant amplitude over
durations that are hundreds of seconds long (Figs 2.2-2.7 and 2.S1-2.S4). To obtain
spectrograms the waveforms are transformed to the frequency domain using short moving
time intervals of 1024 samples with 512 overlapping samples. The results reveal clear
characteristic signals that start with relatively high frequency and drop to low frequency
smoothly with relatively flat portions at the high and low frequency ends (e.g., yellow curve
in Fig. 2.3b). Using different short time windows (e.g., 512 or 2048 samples) and overlap
settings (e.g., 256 or 2014 samples) give essentially the same results. Converting the
waveform in Fig. 2.3a (associated with the white solid box in Fig. 2.2b) into an audio signal
produces a sound that includes a signal resembling the Doppler effect of an aircraft flying-
over (supplementary material Audio 2.S1).
We observe many such signals with clear Doppler effects and frequencies that usually
decrease with time. This has a simple interpretation in terms of radiation from moving
sources. One type of these signals has only one overtone, as the yellow curve in the
spectrogram in Fig. 2.3b, which generally spans a frequency range above 70 Hz. Another type
of signals has multiple overtones in the spectrogram and spans a frequency range above 15
Hz (Fig. 2.5b). More complex signals that appear to be distorted versions of the above two
types (e.g., the right white dashed box in Fig. 2.2b) occur much less frequently. All these
signals have significantly larger amplitudes than the background noise and they last for
minutes.
Fig. 2.6a displays 20 s raw waveforms of the event in Figure 2.3 recorded by the
geophones indicated by profile AA’ in Fig. 2.1c. Fig. 2.6b displays corresponding waveforms
associated with a local earthquake with magnitude M L = 0.5 and hypocentral distance of
about 12 km. The apparent velocity of the coherent phase in Fig. 2.6a is ~300 m/s, while the
apparent velocity of coherent phases in Fig. 2.6b is much higher. The slow apparent velocity
of the air-traffic signal and surface locations of the geophones indicate propagation that is
Doctoral Dissertation Haoran Meng

45
primarily in the air. However, the velocity of high-frequency Rayleigh waves at the site is
~300-400 m/s (Roux et al., 2016; Hillers et al., 2016; Meng & Ben-Zion, 2018a), so there are
likely also induced signals that propagate in the sub-surface material. Fig. 2.S1 shows
waveforms and spectrograms generated by the helicopter event in Fig. 2.5 and recorded by
a 148 m deep sensor in borehole station B946. These results indicate that the acoustic signal
can indeed produce coupled vibrations in the ground. However, the borehole signals are
relatively weak compared to the signals recorded at the surface because of strong
attenuation of the shallow material and possible weak coupling of the acoustic waves to the
ground.

Figure 2.4. (a) A raw waveform of of the marked geophone in Fig. 2.1c for Julian day 156,
22:25:00 to 22:36:40 GMT time. (b) A corresponding spectrogram with solid and dashed
white boxes marking air-traffic events.
Doctoral Dissertation Haoran Meng

46

Figure 2.5. (a) Zoom in waveform for a helicopter event in the white solid box in Fig. 2.4
recorded by the marked station in Fig. 2.1c. (b) A corresponding spectrogram with vertical
axis in log scale. The white arrows in (b) and (c) denote the two strong overtones picked for
analysis. (c) Time and frequency data (black crosses) for two overtones picked in the
spectrogram in (b) and line fit to the data based on equation 1 (red curves). The inferred
parameters for the helicopter are indicated in the boxes.

Doctoral Dissertation Haoran Meng

47
2.2.2 Modeling
Time-frequency analysis and basic analytical results show that the discussed signals have
clear changing-frequency signatures that can be modeled as Doppler effects associated with
airplanes or helicopters flying nearby. To analyze quantitatively the signals highlighted by
the solid white box in Fig. 2.2b, we consider a simple model with an aircraft flying with
constant velocity along a straight path (Fig. 2.1b). The aircraft is assumed to be sufficiently
far from the seismic stations to be represented by a point source radiating acoustic waves
centered on frequency 𝑓 0
. The signal frequency recorded by a seismic station at the surface
is

𝑓 ( 𝑡 ′ ) = 𝑓 0
1
1 +
𝑣 0
𝑐 𝑣 0
𝑡 √𝑙 2
+ ( 𝑣 0
𝑡 )
2
,
(2.6)
where 𝑣 0
is the aircraft velocity, 𝑙 is the shortest distance between the station and the
aircraft path, 𝑡 is time and 𝑐 is the acoustic wave propagation velocity (~343 m/s). A wave
generated by the aircraft at time 𝑡 arrives at the station at time 𝑡 ′

𝑡 ′ = 𝑡 +
√𝑙 2
+ ( 𝑣 0
𝑡 )
2
𝑐 .
(2.7)
Solving equation (2.2) and keeping the causal solution with 𝑡 ′ > 𝑡 gives

𝑡 =
𝑡 ′ −
√
𝑡 ′
2
− (1 −
𝑣 0
2
𝑐 2
)(𝑡 ′
2
−
𝑙 2
𝑐 2
)
1 −
𝑣 0
2
𝑐 2
.
(2.3a)
To model data with arbitrary reference time, we add a reference 𝑡 0
′ value to equation
(2.3a)

𝑡 =
( 𝑡 ′ − 𝑡 0
′ )−
√( 𝑡 ′ − 𝑡 0
′ )
2
− (1 −
𝑣 0
2
𝑐 2
)[( 𝑡 ′− 𝑡 0
′ )
2
−
𝑙 2
𝑐 2
]
1 −
𝑣 0
2
𝑐 2
.
(2.3b)
The above expressions have only 4 unknown parameters that can be easily determined
from the observed data.
Doctoral Dissertation Haoran Meng

48
For a given group of parameters 𝑚 0
= [𝑓 0
𝑣 0
𝑙 𝑡 0
′ ]
T
, the frequency at a receiver
𝑓 ( 𝑡 ′ ) = 𝑓 ( 𝑚 0
; 𝑡 ′ ) is predicted by equations (2.1) and (2.3b). Starting from a reasonable initial
model 𝑚 0
and observed frequency 𝑓 ob
( 𝑡 ′ ) , the final model 𝑚 is inverted by improving the
data fit iteratively using
𝑚 = 𝑚 0
+ ( 𝐺 T
𝐺 )
−1
𝐺 T
[𝑓 ob
( 𝑡 ′ )− 𝑓 ( 𝑚 0
; 𝑡 ′ ) ] , (2.4a)
𝐺 =
∂𝑓 ∂𝑚 |
𝑚 =𝑚 0
. (2.4b)
Taking as example the data shown in Figs 2.3a and 2.3b, we extract the time and frequency
signals of the yellow curve in Fig. 2.3b by finding the frequencies of peak amplitudes for all
times slices in the spectrogram. The black crosses in Fig. 2.3c denote the picked time and
frequency data. Outliers with large misfits are removed iteratively in the inversion for the
Doppler parameters. The red curve in Fig. 2.3c is obtained by inserting equation (2.1) into
(2.3b) and it fits the observed values very well. The fitting parameters (frequency 𝑓 0
~ 131
Hz, velocity 𝑣 0
~ 377 km/h, distance 𝑙 ~ 4.5 km and 𝑡 0
′ ~ 2303 s) have reasonable values for
a small airplane or a passenger airplane that is climbing to (or descending from) a cruising
altitude. Processing the waveforms at the remaining stations indicates a clear trend of 𝑙 and
𝑡 0
′ across the entire array (Fig. 2.A2a and 2.A2b). This is used to determine the flying trace
by the technique discussed in Appendix 2.A.
Another type of air-traffic event is represented by the example in the white solid box of
Fig. 2.4b. The source of this event is an aircraft with multi-frequency radiation. Zooming on
the corresponding time window 1800-2030 s, the spectrogram is displayed in Fig. 2.5b with
a log scale for the frequency. Taking the logarithm on both sides of equation (2.1) we have

log 𝑓 ( 𝑡 ′ ) = log𝑓 0
− log(1 +
𝑣 0
𝑐 𝑣 0
𝑡 √𝑙 2
+ ( 𝑣 0
𝑡 )
2
) .
(2.5)
Equation (2.5) shows that multiple overtones with different 𝑓 0
values in the same aircraft,
moving with same velocity, distance and reference time, have identical second term on the
right side (associated with time). This implies that characteristics associated with different
𝑓 0
values radiated from the same moving object are shifted up or down on spectrograms
using log scale for the frequency. The parallel yellow curves in Fig. 2.5b are therefore
generated by the same moving aircraft. The 4 model parameters can be determined by the
Doctoral Dissertation Haoran Meng

49
same technique used for the previous example event from the observed time and frequency
data of two strong sources. This is illustrated in Fig. 2.5c that shows the picked data (with
outliers removed) and fitted curve of equations (2.1) and (2.3b). In this case, the flying
velocity is 236 km/h and the shortest distance between the flying object and the seismic
station is 1260 m. Fitting separately the two sources with different central frequency 𝑓 0
gives
same results. The obtained flying velocity, height, and existence of multiple central
frequencies that may be produced by different rotor blades indicate that the source
generating this type of data is a helicopter. Similar seismic records associated with
helicopters flying over the Hekla Volcano in Iceland were discussed by Eibl et al. (2015). As
the sensors used in this study have higher Nyquist frequency (250 Hz) than the ones used by
Eibl et al. (2015) (50 Hz), our records cover a broader frequency range and contain many
more overtones from 50 Hz to 220 Hz (Fig. 2.5b).

Doctoral Dissertation Haoran Meng

50
Figure 2.6. (a) Raw 20 s waveforms of air-traffic event in Fig. 2.3 recorded by geophones
in profile AA’ of Fig. 2.1c. (b) Corresponding 20 s waveforms of a local seismic event
(bandpass filtered from 1 to 20 Hz). The red dashed lines denote apparent slowness of 300
m/s.

Figure 2.7. (a) A raw vertical component waveform for a known helicopter event flying
over seismic station JFS4. (b) A corresponding spectrogram with vertical axis in log scale.
The arrows denote the two strong overtones picked for analysis. (c) Time and frequency data
(black crosses) for two overtones picked in the spectrogram in (b) and line fit to the data
based on equation 1 (red curves). The inferred parameters for the helicopter are indicated
in the box.
Doctoral Dissertation Haoran Meng

51

Figure 2.8. (a) A histogram of all detected air-traffic events from Julian day 129 to 157,
2014, with about 31 air-traffic events per day and relatively quiet time during the local night.
(b) A histogram of high quality small candidate earthquakes detected Julian day 146 by Meng
and Ben-Zion (2018). (c) Hourly mean absolute amplitude of stacked waveforms of the entire
dense array in Julian day 146, 2014.

To confirm that the source generating data of the type shown in Figs 2.4 and 2.5 is indeed
a helicopter, we analyze seismic waveforms (Fig. 2.7a) associated with a known helicopter
Doctoral Dissertation Haoran Meng

52
event during maintenance of seismic stations at the Jackass Flat site on the SJFZ. Field notes
indicate that the helicopter flew over station JFS4 (red triangle in Fig. 2.1a), slowed down,
hovered, picked up a seismic sensor package and then flew away from the site. Fig. 2.7a
shows a 140 s long raw waveform of a vertical component sensor associated with the
helicopter flying over the site. The corresponding spectrogram is plotted with log scale for
the frequency in Fig. 2.7b. The results are very similar to those in Figs 2.4 and 2.5, showing
clearly source signatures with multi 𝑓 0
values. Fitting frequency and time data of two 𝑓 0

sources indicates a flying velocity of ~ 85 km/h and distance of ~ 50 m (Fig. 2.7c). Since the
helicopter was very close to the sensor and did not fly with constant velocity along a straight
path, the data fits in this case are not as good as in Fig. 2.5c.
To document the rate of occurrence and duration of air-traffic events in the study area,
we examine manually the entire data of the dense deployment and detect 892 airplane and
helicopter events from Julian day 129 to 157 in 2014. Most events have a single 𝑓 0
source of
the type shown in Fig. 2.3 corresponding to airplanes, and 51 events have clear multiple 𝑓 0

sources of the type shown in Figs 2.5 and 2.8 characteristics of a helicopter source. The
events identified as airplanes usually contain only one overtone (or sometime 2) with
frequencies above 70 Hz, and based on the Doppler analysis parameters they have flying
velocities higher than 300 km/h. The events identified as helicopters contain many
overtones covering frequencies above 15 Hz and have flying velocities below 250 km/h.
Fig. 2.8a displays a histogram of the number of detected air-traffic events for each hour of
the day. As expected, there are considerably more events during the daytime and very few
events from 23:00 to 04:00 local time. The air-traffic events have an average duration of ~
200 s, so together they cover > 7% of the day. The total duration of such signals at locations
with higher air-traffic can cover a larger portion of the day. The detectability of small
earthquakes is highly affected by the background noise level. Fig. 2.8b shows a histogram of
candidate small events detected in the dense deployment data during Julian day 146, 2014,
with a technique involving several steps of stacking of sub-arrays and analyses of the stacked
records with detector functions and beamforming (Meng and Ben-Zion, 2018). There are
considerably more detections from 23:00 to 04:00 local time when the rate of air-traffic
events is low. The variation of the seismic noise during Julian day 146 is estimated by the
hourly mean absolute amplitudes of stacked waveforms across the array (Fig. 2.8c). The anti-
Doctoral Dissertation Haoran Meng

53
correlation of the results with those of Fig. 2.8a indicates that the air-traffic events contribute
significantly to the daily variations of the local noise level. The air-traffic events can affect
especially the detectability of very small earthquakes with high corner frequencies (e.g., M ≤
1) that are the target of recent detection efforts (e.g., Yoon et al., 2015; Inbal et al., 2016; Li et
al., 2018). Other human activities producing noise over lower frequency bands can affect
similarly the detectability of larger seismic events.

2.3 Discussion and conclusions
We observe frequent daily airplane and helicopter events in continuous waveforms
recorded by a dense seismic array at the SGB site on the SJFZ southeast of Anza, California
(Fig. 2.1a). Each event lasts several minutes and resembles in seismograms tremor or
collection of small earthquakes (e.g., Figs 2.2a, 2.3a). However, time-frequency analysis
reveals clear Doppler effects that can be fitted well with a basic analytical solution in terms
of flying velocity, average distance and central frequency (equations (2.1) and (2.3b)).
Airplane events have a simple signature involving typically a single central frequency (Fig.
2.3b). They generate mainly frequencies higher than 70 Hz that can be removed from the
data by low-pass filtering. Helicopter events are associated with several central frequencies
and produce more complex signatures (Figs 2.5b, 2.7b) covering a broad frequency band and,
hence, harder to remove form data. The air-traffic events increase significantly the local
noise level and reduce the detectability of small earthquakes (Fig. 2.8). They should not affect
typical noise-based imaging studies that use frequencies lower than 1 Hz, but can affect
noise-based applications with dense arrays and frequencies higher than 15 Hz (e.g., Fig. 9 of
Ben-Zion et al., 2015).
The acoustic wavelengths of typical air-traffic events (< 5 m with velocity of 343 m/s and
frequency > 70 Hz) are shorter than the interstation spacing (~ 10-30 m), so beamforming
can not be used to determine the sources trajectories. However, the flying traces of airplanes
can be determined by the technique described in Appendix 2.A. The detectable distance of
air-traffic events can be found from the simple expression √𝑙 2
+ ( 𝑣 0
𝑡 𝑑 /2)
2
, where 𝑡 𝑑 is the
duration of the air-traffic signal. The detectable distance is related to the strength of the
waves radiated from the source and the background noise level. As examples, using 𝑡 𝑑 ~ 150
Doctoral Dissertation Haoran Meng

54
s for the airplane event in Fig. 2.3 and 𝑡 𝑑 ~ 220 s for the helicopter event in Fig. 2.5, the
corresponding distances are ~ 9.0 km and ~ 7.3 km, respectively. Most detectable distances
for other air-traffic events are in the range 6-10 km.
Most airplane events have a single source 𝑓 0
, but some have two overtones as shown in
the white dashed box in Fig. 2.4b. The second source with higher 𝑓 0
is truncated by the corner
frequency of sensor sensitivity at ~ 220 Hz. Both airplanes and helicopters may generate
overtones with frequencies > 250 Hz that are not recorded by the sensors. The radiation
sources from airplane events are primarily the engines, perhaps augmented by additional
weaker vibrations associated with the aircraft body and propellers (for small airplanes). The
radiation sources from helicopters are more complex and are related to the engines and
rotating propellers. More than 10 overtones of the helicopter event in Fig. 2.5b demonstrate
the possible multiple resonances produced by helicopters. Various conditions can distort the
basic Doppler characteristics incorporated in the simple equations (2.1)-( 2.3) and produce
complex features of the type shown in right white dashed box in Fig. 2.2b. These include
significant changes in the aircraft speed, flying direction, interaction of small aircraft with
strong wind gusts etc. Figs 2.S2 to 2.S4 provide additional examples of complex signals
produced by air traffic events.
In the relatively remote location of the study area there are several air-traffic events per
hour, except during local nighttime (Fig. 2.8a), producing together ground motion over 7%
of the day. In locations closer to airports and large cities (e.g., Long Beach, Los Angeles) the
air-traffic events can cover a significantly larger portion of the day. As example, the Los
Angeles International airport handled in 2016 on average 1910 aircraft landings and takeoffs
per day (http://www.lawa.org/welcome_lax.aspx?id=40). Assuming that each event lasts
about ~ 200 s, and that the events have 75% overlap, signals generated by airplanes (and
helicopters) exist in the entire portion of the day (e.g., 16-18 hours) with significant air-
traffic activity. The ground vibration generated by air-traffic events can mask and prevent
detection of signals associated with very small seismic sources (Fig. 2.8b), or in worse cases
misidentified as small earthquakes or tremor (e.g., Eibl et al. 2015). To highlight the
importance of the results for current efforts to detect small earthquakes, we estimate below
the expected duration of ground motion generated in a day by small local earthquakes.

Doctoral Dissertation Haoran Meng

55

Figure 2.9. Measured durations of earthquakes with different local magnitudes (black
crosses) and a fitted line (red) associated with the scaling relation indicated in the box.

The total duration 𝑇 of daily earthquakes in the magnitude range 𝑀 1
≤ 𝑀 ≤ 𝑀 2
is given
by the integral
𝑇 = ∫ 𝜏 ( 𝑀 )
𝑀 2
𝑀 1
∙ 𝑛 ( 𝑀 ) 𝑑𝑀 , (2.6)
where 𝜏 ( 𝑀 ) is the duration of an earthquake with magnitude 𝑀 and 𝑛 is the number of
events around magnitude 𝑀 given by the Gutenberg-Richter relation
log𝑛 ( 𝑀 ) = 𝑎 − 𝑏𝑀 . (2.7)
The b-value of the Gutenberg-Richter relation is typically around 1 while the a-value
depends strongly on the location and time (e.g. Zaliapin & Ben-Zion, 2013, 2016). To evaluate
the expected daily earthquake duration for the study area, we use the technique illustrated
in Appendix 2.B to measure the duration of 266 cataloged earthquakes (Hauksson et al., 2012)
in the magnitude range 0 ≤ 𝑀 ≤ 3 recorded by 3-component station SGBF0 (Fig. 2.9). The
results can be summarized by the scaling relation
log𝜏 ( 𝑀 ) = 0.41𝑀 + 0.89 . (2.8)
Doctoral Dissertation Haoran Meng

56
Using logt in the scaling relation avoids getting negative duration for earthquakes with
small magnitudes. Combining equations (2.6)-( 2.8) we get,
𝑇 = ∫ 10
0.41𝑀 +0.89
𝑀 2
𝑀 1
∙ 10
𝑎 −𝑏𝑀
𝑑𝑀 . (2.9)
To evaluate (9) for daily seismicity in the trifurcation area of the SJFZ, we use 𝑏 = 1 and 𝑎
= 1.5/day. Integrating with these constants equation (2.9) over the magnitude range 0 ≤
𝑀 ≤ 3 gives total daily earthquake duration of 178 s, which covers 0.2% of the day. Reducing
𝑀 1
to –1 and –2, and assuming the scaling relation of equation (2.8) can be extrapolated to
lower magnitude than those used in Fig. 2.9, increases the estimated total duration covered
by earthquakes to 0.8% and 3.2% of the day, respectively. Reducing further the lower
magnitude limit to –4.54 gives a total duration of a full day. These total earthquake duration
estimates are gross upper bounds, since many small regional events may overlap in time and
many small events not sufficiently close to the recording site will not produce observable
signals (e.g., Kwiatek and Ben-Zion, 2016). Increasing the upper magnitude limit to 5 has
negligible affect on the results because the rate of higher magnitude events is very low. We
note that the seismicity rate, and hence the 𝑎 -value, in the trifurcation area is the highest in
Southern California (Ross et al., 2017), so corresponding estimates of earthquake durations
in most other regions will be smaller than those given above.
The results of this study demonstrate that seismic arrays can be used to detect and
quantify flight parameters of aircrafts (these should include low flying and stealth airplanes
that evade radar detection). The presented results contribute to the growing recognition of
different types of surface and atmospheric events (e.g., car-traffic, trains, storms, wind
turbines, air-traffic) producing signals that are contained in observed seismograms (e.g.,
Salvermoser et al., 2015; Eibl et al. 2015; Tanimoto & Valovcin, 2015; Neuffer & Kremers, 2017;
Fuchs & Bokelmann, 2017). In the study area near Anza California, the total duration of air-
traffic signals (> 7% of the day) is likely larger than the duration covered by earthquake
signals. In many other places, and especially large urban areas with highly active airports,
the duration of air-traffic signals can exceed 50% of the day and be significantly larger that
the duration of seismic events. Continuing efforts to recognize characteristic signatures of
various anthropogenic and atmospheric sources of ground motion, and separating these
Doctoral Dissertation Haoran Meng

57
sources from weak sub-surface earthquakes and tremor, will help developing refined
catalogs of seismic sources and improved understanding of the observed seismic wavefield.

Figure 2.A1. A conceptual diagram of an aircraft flying over two example stations.

Appendix 2.A. Finding the flying traces of airplanes
To model the flying trace of an aircraft, we consider a reference point 𝑟 ⃑ = [𝑥 0 𝑧 ] on
the trace and a direction vector 𝑒 ̂= [
𝑒 𝑥 𝑒 𝑦 𝑒 𝑧 ] with 𝑒 𝑥 2
+ 𝑒 𝑦 2
+ 𝑒 𝑧 2
= 1. The reference
times 𝑡 0
′ for all stations gives good constraint for the direction vector. Fig. 2.A1 displays the
geometry of the flying trace and two example stations. At time 𝑡 ′ = 𝑡 01
′, station S 1 receives
the signal of the aircraft radiated at time 𝑡 1
at point A
1
. From equation (2.3b), when 𝑡 ′ = 𝑡 01
′,
𝑡 = −𝑙 1
/√𝑐 2
− 𝑣 0
2
. This means it takes the acoustic wave 𝑙 1
/√𝑐 2
− 𝑣 0
2
to propagate from A
1

to S
1
. Therefore, the length of line segment |𝐴 1
𝑆 1
| = 𝑐 𝑙 1
/√𝑐 2
− 𝑣 0
2
, the length of line
segment |𝐴 1
𝐵 1
| = 𝑣 0
𝑙 1
/√𝑐 2
− 𝑣 0
2
, and

𝑡 01
′ = 𝑡 1
+
|𝐴 1
𝑆 1
|
𝑐 = 𝑡 1
+
𝑙 1
√𝑐 2
− 𝑣 0
2
.
(2.A1)
By performing similar estimates for station S
2
, the length of line segment |𝐴 2
𝐵 2
| =
𝑣 0
𝑙 2
/√𝑐 2
− 𝑣 0
2
, the length of line segment |𝐴 2
𝑆 2
| = 𝑐 𝑙 2
/√𝑐 2
− 𝑣 0
2
, and
Doctoral Dissertation Haoran Meng

58

𝑡 02
′ = 𝑡 2
+
|𝐴 2
𝑆 2
|
𝑐 = 𝑡 2
+
𝑙 2
√𝑐 2
− 𝑣 0
2
.
(2.A2)
Subtracting equation (2.A1) from (2.A2),

𝑡 02
′− 𝑡 01
′ = 𝑡 2
− 𝑡 1
+
𝑙 2
− 𝑙 1
√𝑐 2
− 𝑣 0
2
.
(2.A3)
The location of station S 1 is 𝑟 ⃑
1
, the location of station S 2 is 𝑟 ⃑
2
, and the length of |𝐵 1
𝐵 2
| = 𝑒 ̂ ∙
( 𝑟 ⃑
2
− 𝑟 ⃑
1
) . Based on simple geometrical consideration,

|𝐴 1
𝐴 2
| = |𝐴 1
𝐵 1
| + |𝐵 1
𝐵 2
| − |𝐴 2
𝐵 2
|
=
𝑣 0
𝑙 1
√𝑐 2
− 𝑣 0
2
+ 𝑒 ̂ ∙ ( 𝑟 ⃑
2
− 𝑟 ⃑
1
)−
𝑣 0
𝑙 2
√𝑐 2
− 𝑣 0
2
.
(2.A4)
Substituting 𝑡 2
− 𝑡 1
in equation (2.A3) with |𝐴 1
𝑆 1
|/𝑣 0
,
𝑡 02
′− 𝑡 01
′ =
𝑒 ̂ ∙ ( 𝑟 ⃑
2
− 𝑟 ⃑
1
)
𝑣 0
, (2.A5)
where the reference times (𝑡 01
′ and 𝑡 02
′) are given by the data and equation (2.4), and 𝑣 0
is
determined by the mean of the observed flying velocity for all stations. By picking a reference
station, we can build a matrix for inversion of the direction vector. Among the three
components of the direction vector, only two are independent. Using the constraint 𝑒 𝑥 2
+
𝑒 𝑦 2
+ 𝑒 𝑧 2
= 1 and equation (2.4), the inversion for the direction vector 𝑒 ̂ converges quickly
within a couple of iterations.

Doctoral Dissertation Haoran Meng

59

Figure 2.A2. (a) Distribution of the distance l (colored circles) based on data at all stations
for the example event 1 in Fig. 2.3. (b) Corresponding t' 0 distribution. (c) A northward cross
section showing the height of the flying trace. The black arrows mark the inferred flying trace
of the airplane.

There are only 2 independent parameters for the reference point 𝑟 ⃑ of the flying trace. By
fitting the data with Doppler effect, the distance between each station and the flying trace is
displayed as the colored circles in Fig. 2.A2a. From the Pythagoras theorem, the distance is

𝑙 = √( 𝑟 ⃑ − 𝑟 ⃑
st
)
2
− [( 𝑟 ⃑ − 𝑟 ⃑
st
)∙ 𝑒 ̂ ]
2
.
(2.A6)
By using the same inversion technique equation (2.4), 𝑟 ⃑ can be determined by the direction
vector ˆ e , station location 𝑟 ⃑
st
and the observed distances. The flying trace of the example
event in Fig. 2.3 is marked by black arrows in Fig. 2.A2. This flight trace moves to the
northeast with a height of about 4.5 km.

Appendix 2.B. Estimating earthquake durations
Doctoral Dissertation Haoran Meng

60
The durations of earthquakes depend on site conditions. To estimate a basic scaling
relation 𝜏 ( 𝑀 ) at the SGB site on the San Jacinto fault zone, we measure the durations of 266
earthquakes from the catalog of Hauksson et al. (2012) recorded by the strong ground
motion station SGBF0 (Fig. 2.1c) from Julian day 129 to 157, 2014. The procedure is
illustrated in Fig. 2.B1 for an 𝑀 L
2.41 earthquake with epicentral distance of 95 km. For each
event we first band-pass filter the accelerograms at 2-15 Hz to increase the SNR. Figs 2.B1a
and 2.B1b show, respectively, the filtered 3-component waveforms and vector sum of the 3
components for this event. A 1D model obtained by averaging the 3D velocity model of Lee
et al. (2014) for the Southern California region shown in Fig. 2.1a is used to calculate
predicted P and S arrivals (vertical blue and red lines). To estimate the end time of the coda
waves, we calculate a moving mean of the vector sum to get a smooth envelop of the
waveform (Fig. 2.B1c). Since the earthquake durations can have several orders of time
intervals, we use for each event a moving window length that depends on the source
duration.
The source duration is estimated as follows. We first estimate the rupture size of each
event, assuming it is approximately a circular crack with radius 𝑟 sustaining a uniform strain
drop ∆𝘀 in a Poissonian solid. The seismic potency (moment/rigidity) in this case
(appropriate for the analyzed small events) is given by (Eshelby, 1957; Ben-Zion, 2008),
𝑃 0
=
16
7
∆𝘀 𝑟 3
. (2.B3)
The seismic potency and magnitude of earthquakes spanning a relatively small range of
magnitudes are related empirically by a relation of the type (e.g., Hanks & Kanamori, 1979;
Ben-Zion & Zhu, 2002),
log𝑃 0
= 𝑑 ∙ 𝑀 + 𝑒 . (2.B3)
For small earthquakes in Southern California, 𝑑 = 1 and 𝑒 = −4.7 with potency units in
km
2
·cm (Ben-Zion & Zhu, 2002). Combining equation (2.B1) and (2.B2) and assuming an
average strain drop ∆𝘀 = 10
−4
, the radius in unit of km of an earthquake with magnitude 𝑀
is
𝑟 ( 𝑀 ) = [
10
−1
( 16/7)
] ∙ 10
( 𝑑 ∙𝑀 +𝑒 )
3
, (2.B3)
Doctoral Dissertation Haoran Meng

61
where the 10
-1
factor stems from using ∆𝘀 = 10
−4
and including a unit conversion of 𝑃 0

from km
2
×cm to km
3
. Assuming an average rupture propagation velocity of 0.9𝑉 𝑆 , with shear
wave velocity 𝑉 𝑆 ≈ 3 km/s, the source duration for an earthquake with magnitude 𝑀 is
𝑇 sd
=
𝑟 ( 𝑀 )
0.9𝑉 𝑆 = 0.0162∙ 10
( 𝑑 ∙𝑀 +𝑒 )
3
. (2.B4)
In the performed analysis we use 500∙ 𝑇 sd
as the length of the time window to calculate
the moving mean of the vector sum (Fig. 2.B1c). Using somewhat smaller and larger window
length gives essentially the same results. The noise level is estimated as the RMS of a 10 s
moving mean window before the P arrival (red horizontal dashed line in Fig. 2.B1c). The end
of the coda waves is picked when the moving mean drops to noise level (vertical green line).
The duration of each analyzed earthquake is estimated as the time between the P arrival and
the end of the coda waves.

Figure 2.B1. An illustration of estimating duration for example event with M L 2.41 and
epicenter distance of 95 km recorded by a 3-composnent station SGBF0. (a) Z, N and E
component waveforms band-pass filtered between 2 and 15 Hz. (b) The vector sum of the
three component waveforms. (c) The moving mean of the vector sum. The red dashed
Doctoral Dissertation Haoran Meng

62
horizontal line denotes the noise level, while the blue, red and green vertical lines show the
P, S arrivals and the end of the coda.

Figure 2.S1. (a) A raw waveform of a vertical component sensor at depth of 148 m in
borehole station B946 (Fig. 2.1c) for Julian day 156, 22:25:00 to 22:36:40 GMT time. (b) A
corresponding spectrogram with solid and dashed white boxes marking air-traffic events.

Doctoral Dissertation Haoran Meng

63

Figure 2.S2. (a) A raw waveform recorded by the geophone marked in Fig. 2.1c for Julian
day 156, 20:46:40 to 21:26:40 GMT time. (b) A corresponding spectrogram with solid white
boxes marking air-traffic events.

Figure 2.S3. (a) A raw waveform recorded by the geophone marked in Fig. 1c for Julian
day 146, 00:00:00 to 01:00:00 GMT time. (b) A corresponding spectrogram with solid white
boxes marking air-traffic events.
Doctoral Dissertation Haoran Meng

64

Figure 2.S4. (a) A raw waveform recorded by the geophone marked in Fig. 2.1c for Julian
day 146, 16:00:00 to 17:00:00 GMT time. (b) A corresponding spectrogram with solid white
boxes marking air-traffic events.

jgrb52814-sup-0002-2017JB015240-s2.wav
Audio 2.S1. An audio signal converted from the waveform in Fig. 2.3a by speeding the
signal by a factor of 10. The source is identified as an airplane.

jgrb52814-sup-0003-2017JB015240-s3.wav
Audio 2.S2. An audio signal converted from the waveform in Fig. 2.5a by speeding the
signal by a factor of 10. The source is identified as a helicopter.

Doctoral Dissertation Haoran Meng

65
3. Detection of random noise and anatomy of continuous seismic waveforms in
dense array data near Anza California
(Meng et al. 2019 in revision)

3.1 Introduction
Seismic ground motion is generated by a variety of mechanisms including tectonic sources
such as earthquakes and tremor (e.g., Ross et al. 2017; Shelly et al. 2007), interactions of
ocean waves with the solid Earth (e.g., Gerstoft & Tanimoto 2007; Hillers et al. 2013), storms
and wind shaking of structures, trees and other obstacles (e.g., Tanimoto & Valovcin 2015,
Johnson et al. 2019) and various anthropogenic sources such as trains, cars, airplanes,
helicopters and wind turbines (e.g., Salvermoser et al. 2015; Eibl et al. 2015; Neuffer &
Kremers 2017; Meng & Ben-Zion 2018a; Inbal et al. 2018). It is generally assumed that
continuous waveforms processed to remove impulsive sources (e.g. Bensen et al. 2007)
consist primarily of a random wavefield, referred to as the ambient seismic noise, that is
suitable generally for imaging and monitoring of sub-surface structures (e.g., Shapiro et al.
2005; Brenguier et al. 2008; Zigone et al. 2015). Extracting accurate Green’s functions from
the ambient seismic noise requires a stronger condition that the wavefield is fully diffuse
(e.g., Weaver 1982; Shapiro et al. 2000; Sá nchez-Sesma et al., 2008). Demonstrating that a
wavefield is diffuse requires analysis with techniques of the type developed by Margerin et
al. (2009), Sens-Schonfelder et al. (2015) and Liu and Ben-Zion (2016).
In the present paper we develop techniques based on cross-correlations and amplitude
spectra of moving time windows to separate wave packets with random signals from
waveform sections having structured signals. The results help to clarify general aspects of
observed seismic waveforms, and the separated random/not-random signals can be used in
various additional analyses involving detection of small seismic events, imaging and
monitoring sub-surface structures, etc. The study employs data of a dense seismic array with
1108 vertical component 10 Hz geophones at the Sage Brush Flat (SGB) site on the Clark
branch of the San Jacinto fault zone (SJFZ) southeast of Anza, California (red box Fig. 3.1a).
The array recorded continuously at 500 Hz from 7 May 2014 to 13 June 2014 in an area of
about 600 m × 600 m (Fig. 3.1b) with nominal sensor spacing of 10 m normal to the fault and
30 m along-strike (Ben-Zion et al. 2015). The dense array data have been used for detailed
Doctoral Dissertation Haoran Meng

66
imaging of subsurface structures at different depth ranges (Hillers et al. 2016; Qin et al. 2018;
Mordret et al. 2019) and detection of microearthquakes (Meng & Ben-Zion 2018b; Gradon et
al. 2019).

Figure 3.1. (a) Regional map showing the dense array on the San Jacinto fault zone (red
square) with the regional seismic stations (black triangles) in southern California. (b)
Colored circles show the median absolute amplitude of waveforms recorded at the array
stations on Julian day 146, 2014. The geophones with analyzed data are marked by squares.
The background grey colors represent topography.

Doctoral Dissertation Haoran Meng

67
As part of an effort to understand the main contributions to ground motion observed at
the SGB site, Meng and Ben-Zion (2018b) analyzed waveforms generated by airplanes and
helicopters, and found that air-traffic events occupy >7% of the recorded wavefield. Johnson
et al. (2019) showed that wind interactions with obstacles above the surface produce ground
velocities larger than those expected to be generated by earthquakes with magnitudes M ≤
1.5 about 6% of the day. The air-traffic and wind-related events produce earthquake- and
tremor-like waveforms that stand out clearly from the seismic noise field recorded at the site.
Less frequent sources of seismic signals that are visually distinguishable from the noise are
genuine earthquakes and local car-traffic (Fig. 3.2).

Figure 3.2. Waveforms in count generated by (a) a local M 2.29 earthquake, (b) an air
traffic event, (c) a car event and (d) a tree-shaking event driven by wind gust. (e)-(h)
Corresponding spectral amplitudes of (a)-(d). (i)-(l) Corresponding spectrograms of (a)-(d).

Inspection of waveforms at the SGB site shows a variety of emergent and impulsive signals
occurring throughout the day. Quantitative analysis of the main classes of waveforms can
improve the understanding of the observed ongoing ground motion and provide useful
information for refined seismic imaging and detection studies. Toward this goal, we use the
Doctoral Dissertation Haoran Meng

68
dense SGB data to separate Random-Noise (RN) signals from wave packets containing Not-
Random-Noise (NRN) data produced by various sources such as air- and car-traffic, wind,
and earthquakes. To identify wave packets that may be labeled as RN, we use 1000 1-s time
windows for each hour and select those with the lowest Root Mean Squared (RMS) amplitude
as initial candidate RN wave packets. The selection of RN signals is refined iteratively using
cross-correlations and the distribution of amplitude spectra. The method allows us to
separate each 1-s interval of the continuous waveforms into RN signals, NRN with higher
amplitude in time domain or some structure in the frequency domain, and mixed signals
where neither the RN nor NRN features are dominant. The RN windows can be used to obtain
better estimates of noise spectra, monitor hourly and longer changes of the background
noise level, derive structural information, and measure the degree to which the noise field is
diffuse. The RN windows detected by our method can also be used as training datasets for
machine learning analyses of phase picking and signal classification. The time windows
containing NRN signals may be further analyzed to select templates with specific signals such
as those generated by air/car-traffic, wind, earthquakes or tectonic tremor for detection and
classification of such signals.
Doctoral Dissertation Haoran Meng

69

Figure 3.3. A workflow diagram summarizing the main steps of the classification method.
The workflow includes multiple iterations to update and stabilize RN, not RN and the mixture.

3.2 Analysis
Seismograms are analyzed by using cross-correlations of waveform segments and a
classification algorithm to extract parameters in time and frequency domains. Data recorded
by station R3010 (Fig. 3.1) is used to illustrate the method of separating continuous seismic
waveforms into RN, NRN, and a mixture of the two. Station R3010 is in the center of the array
Doctoral Dissertation Haoran Meng

70
and has good data quality to demonstrate the method. Several iterations are implemented to
update and stabilize the random noise templates, mean noise spectra and separation results.
The methodology is summarized in a workflow diagram (Fig. 3.3). We analyze the seismic
records in short periods, 𝑇 𝑆 = 1 hr, and cut the waveforms into 1 s non-overlapping windows,
T, with 500 samples each to stably estimate the amplitude spectra and cross-correlation
coefficients. The window length 𝑇 should be adjusted for other data with different sampling
rates. The mean and linear trend are removed and waveforms are high-pass filtered at 2/𝑇
(2 Hz for this study) to include at least two cycles in each window.
To investigate potential techniques and parameters that can separate RN and NRN signals,
we first generate synthetic RN by Fourier transforming a low amplitude 500 s waveform
recorded at night on sensor R3010, replacing the phase information with a uniform
distribution of values from 0 to 2π, and performing an inverse Fourier transform on the
amplitude spectrum and randomized phase (Figs 3.4a and 3.4b). Next, we extract a local
M2.29 earthquake and an air-traffic event recorded by sensor R3010 to use as example NRN
waveforms (Figs 3.4c and 3.4d). The waveforms are cut into 1 s non-overlapping windows
to obtain 500 synthetic RN windows, 25 earthquake windows, and 50 air-traffic event
windows (grey boxes in Figs 3.4c and 3.4d).
The Maximum of the Absolute value of Cross-correlation Coefficients (MACCs) is
computed for all window pairs. Fig. 3.4e shows the distributions of MACCs for RN – RN
windows (grey), RN – earthquake windows (blue), and RN – air traffic windows (red). The
MACCs of RN – RN windows have a well-defined interval between 0.09-0.35, while the
distributions of RN – NRN systematically shift to intervals with lower median values
demonstrating a simple statistical difference between the RN and NRN waveforms. The
median of MACCs for these target windows is the first extracted time domain parameter to
differentiate the RN and NRN signals. Fig. 3.5 shows that the distributions and Standard
Deviations (STD) for earthquake – earthquake windows (green), air traffic – air traffic
windows (purple), and earthquake – air traffic windows (brown) are spread over a broader
range than RN – NRN windows. This indicates that the NRN windows contain both very
similar and highly dissimilar waveforms. For example, the spectrogram of the air-traffic
event (Fig. 3.2j) has a changing peak frequency due to the Doppler effect (Meng & Ben-Zion
2018a) that results in high MACCs for similar window pairs. The STDs of the MACCs for NRN
Doctoral Dissertation Haoran Meng

71
– NRN windows are significantly higher than the STDs of MACCs of RN – NRN windows (Fig.
3.5b). This feature is used as the second extracted time domain parameter to differentiate
the RN and NRN waveforms.

Figure 3.4. (a) Amplitude spectrum of a 500 s waveform at night with low amplitudes. (b)
Synthetic RN waveform using the amplitude spectrum in (a). (c) Waveform of a local M 2.29
earthquake recorded by the geophone R3010 in the center of the array. (d) Waveform of an
air-traffic event recorded by the same geophone. (e) Histograms of the MACCs of RN – RN
windows, RN – earthquake windows and RN – air-traffic windows. The median values for
each distribution are shown above. The gray boxes show the NRN waveform segments used
for the cross-correlation calculation.

Doctoral Dissertation Haoran Meng

72

Figure 3.5. (a) Histograms of the MACCs of earthquake – earthquake windows, air-traffic
– air-traffic windows, and earthquake – air-traffic windows. (b) Standard deviations of
histograms for different window pairs in (a) and Fig. 3.3e. The standard deviation of the NRN
– NRN windows are significantly larger than the standard deviation of RN – NRN windows.

To adopt the method to continuous seismic records, we assume that the RN is quasi-
stationary (i.e., statistically time invariant) within a short period of one hour. This implies
that the amplitude spectrum of RN only changes slightly within a short period 𝑇 𝑆 .
Supplementary Fig. 3.S1 shows that the RN is generally quasi-stationary within time
intervals of about 1 hour, although the noise level varies from day to night. In the frequency
domain the spectra of RN and NRN windows deviate significantly from each other. The L2
norm of the deviation of an amplitude spectrum of a given target window and the mean
amplitude spectra of RN candidates provides a third extracted frequency domain parameter
to differentiate between RN and NRN windows. These three parameters in time and
frequency domains are used to classify time windows of RN.
Doctoral Dissertation Haoran Meng

73

Figure 3.6. (a) Example raw waveform recorded by sensor R3010 on Julian day 146, 2014
from 14:00 to 15:00 (local time). (b) The median of MACCs of detected RN windows and
target waveforms (𝐶 𝑖 MDN
). (c) L2 norm deviations from the mean amplitude spectra of RN
windows. (d) The STD of MACCs of NRN windows and target waveforms (𝐶 𝑖 STD
).

The corresponding three parameters extracted in time and frequency domains are shown
in Fig. 3.6 for the 1 hour waveform from sensor R3010 in Julian day 146 from 14:00 to 15:00
(local time). The data contain frequent impulsive signals throughout the hour, and the
median MACCs of target and RN windows generally show a decrease during the strong bursts.
In the following we describe how we iteratively build libraries of RN, NRN and mixture
signals to separate the continuous waveforms into different classes.
The sets of RN, NRN and the undetermined mixtures are denoted by 𝑁 , 𝑆 and 𝑀 , where
Doctoral Dissertation Haoran Meng

74

𝑁 = {𝑖 : 𝑖 ⊆ indexes of all RN windows },
𝑆 = {𝑖 : 𝑖 ⊆ indexes of all NRN windows },
𝑀 = {𝑖 : 𝑖 ⊆ indexes of undertermined mixtrues }.
(3.1)
These three sets satisfy

𝑁 ∪ 𝑆 ∪ 𝑀 = 𝑈 = {𝑖 ; 𝑖 is an integer , 1 ≤ 𝑖 ≤ 𝑇 𝑆 /𝑇 },
𝑁 ∩ 𝑆 = 𝑁 ∩ 𝑀 = 𝑆 ∩ 𝑀 = ∅,
(3.2)
where 𝑇 𝑆 /𝑇 is 3600 in this study. The element 𝑖 is linked to the waveform of 𝑖 -th window and
its MACCs with other windows. Since windows with a small RMS value are more likely to be
RN, we first build an initial RN library 𝑁 by taking the indexes of 1000 windows with the
lowest RMS. The initial NRN library 𝑆 contains indexes of the 1000 windows with larger RMS
than the others. The initial library of the undetermined mixture, which is a transition from
RN to NRN, is the complement of the 𝑁 and 𝑆 sets.
To avoid redundant computations, we calculate the MACCs for all pairs of the 3600
windows. The MACC of the 𝑖 -th and 𝑗 -th windows is defined as
𝐶 𝑖𝑗
= max (abs(
𝑓 𝑖 ( −𝑡 )∗ 𝑓 𝑗 ( 𝑡 )
STD ( 𝑓 𝑖 ( 𝑡 ) )∙ STD ( 𝑓 𝑗 ( 𝑡 ) )
)), (3.3)
where 𝑓 𝑖 ( 𝑡 ) and 𝑓 𝑗 ( 𝑡 ) are the waveforms of the 𝑖 -th and 𝑗 -th windows, 0 ≤ 𝑡 ≤ 𝑇 , STD ( 𝑓 𝑖 ( 𝑡 ) )
is the STD of all the data points of the 𝑖 -th window. For the 𝑖 -th window, we compute the
median of the MACCs of the target window and RN windows
𝐶 𝑖 MDN
= median( 𝐶 𝑖𝑛
) |
𝑛 ∈𝑁 ,𝑛 ≠𝑖 , (3.4)
and the STD of MACCs of the target window and NRN windows

𝐶 𝑖 STD
= STD ( 𝐶 𝑖𝑛
) |
𝑛 ∈𝑆 ,𝑛 ≠𝑖 .
(3.5)
The outliers are removed by excluding windows with 𝐶 𝑖 STD
≥ 1.1 ∗ median( 𝐶 𝑗 STD
) , 𝐶 𝑖 MDN
≥
1.1 ∗ meidan( 𝐶 𝑗 MDN
) and 𝐶 𝑖 MDN
≤ 0.9 ∗ median( 𝐶 𝑗 MDN
) , where 1 ≤ 𝑗 ≤ 3600 . The method is
not sensitive to the initial RN library. All libraries are updated iteratively and stabilize
quickly. The iteration stops when the total number of changed elements in libraries 𝑁 and 𝑆
from two consecutive iterations are < 0.5% ∙ 𝑇 𝑆 /𝑇 . The 𝐶 𝑖 MDN
and 𝐶 𝑖 STD
for each window of
Doctoral Dissertation Haoran Meng

75
the target waveform in Fig. 3.6a are presented in Figs 3.6b and 3.6d, respectively. The
waveform of each window is cosine-tapered with a width of 5% on each end to stably
estimate the amplitude spectrum using Fourier transform. Fig. 3.6c shows the L2 norm
deviations from the mean spectra of RN windows.

Figure 3.7. (a) A 2D histogram of the example waveform (Fig. 3.6a) in the parameter
space of 𝐶 𝑖 MDN
and spectral deviation. The colors denote the normalized data density. (b)
Same as (a) but colored by 𝐶 𝑖 STD
. (c) Same as (a) with colors denoting the normalized data
density weighted by the inverse of the 𝐶 𝑖 MDN
. The contour lines show the criteria to
determine RN, not RN and mixture signals. (d) histograms of MACCs of RN – RN windows,
NRN – NRN and mixed – mixed windows.

Doctoral Dissertation Haoran Meng

76
Results of the classification algorithm performed in the parameter space of the spectral
deviation and 𝐶 𝑖 MDN
is shown in Fig 3.7. Each point in Fig. 3.7a represents a time window and
the colors denote the normalized data points density 𝜌 𝑖 , which is computed by the total
number of data points within a rectangle area centered at the target points with length of
0.2 ∙ STD ( 𝐶 𝑖 STD
) and width of 0.2 ∙ STD ( 𝐶 𝑖 MDN
). The hourly RN is quasi-stationary so the
results should cluster tightly (with high 𝜌 𝑖 ) in the parameter space. Fig. 3.7b shows the same
diagram but colored by the 𝐶 𝑖 STD
. The windows with smaller 𝐶 𝑖 STD
are more likely to be RN.
To make use of both 𝜌 𝑖 and 𝐶 𝑖 STD
, Fig. 3.7c uses a weighted density as

𝜌 𝑖 𝑊 =
𝜌 𝑖 𝐶 𝑖 STD
max (
𝜌 𝑗 𝐶 𝑗 STD
)
.
(3.6)
This combines 𝐶 𝑖 MDN
, 𝐶 𝑖 STD
and spectral deviations extracted from the time and frequency
domains, and therefore is more sensitive to the difference between RN and NRN windows.
Time windows with parameters in the tight cluster with high weighted density 𝜌 𝑖 𝑊 (colored
by dark red in Fig. 3.7c) are identified as RN, time windows that deviate significantly from
the RN cluster are identified as NRN, and windows with values in between are identified as
a mixture. The contours in Fig. 3.7c show the criteria to determine RN, NRN and the mixture.
The corresponding sets are then updated as

𝑁 = {𝑖 : 𝜌 𝑖 𝑊 ≥ 0.45},
𝑆 = {𝑖 : 𝜌 𝑖 𝑊 ≤ 0.15},
𝑀 = {𝑖 : 0.15 < 𝜌 𝑖 𝑊 < 0.45}.
(3.7)
With the updated RN library 𝑁 , we recalculate the mean amplitude spectrum of RN and then
update the spectral deviations for each window. The 𝐶 𝑖 MDN
and 𝐶 𝑖 STD
are updated
correspondingly with the new libraries 𝑁 and 𝑆 . This process is repeated for a few iterations
until the three libraries are stable. Solid lines in Fig. 3.S2 show the mean amplitude spectra
of RN for example waveforms (Fig. 3.6a) in different iterations. As a validation step, Fig. 3.7d
presents the histograms of MACCs of RN – RN windows, NRN – NRN and mixed – mixed
windows. The MACCs of RN – RN windows span a range of 0.09 of 0.35 while the MACCs of
not NRN – NRN windows span a much broader range. These results are consistent with the
histograms of the MACCs of the synthetic RN, earthquake and, air-traffic shown in Fig. 3.4e
Doctoral Dissertation Haoran Meng

77
and 3.5a. The thresholds in equation (3.7) should be adjusted to minimize the difference of
distribution, if the histograms of the MACCs of the detected RN (e.g. Fig. 3.7d) and the
synthetic RN (e.g. Fig 3.4e) are considerably different.

Figure 3.8. (a) Example raw waveform recorded by R3010 on Julian day 146, 2014 from
14:00 to 15:00 (local time). The RN, NRN and the mixture windows are colored in red, black
and blue. (b) Corresponding spectrogram of the waveform in (a) with colors denoting power
(dB). The data marked by the red horizontal bar between (a) and (b) illustrate that windows
with lower amplitudes are not necessarily RN, as evidenced by the 3 overtones in the
spectrogram. (c) Example raw waveform recorded by sensor R3010 at day 146 from 01:00
to 02:00 (local time). (d) Corresponding spectrogram of waveform in (c). Green vertical lines
in (a) and (c) mark the small earthquakes detected by Meng & Ben-Zion (2018b).
Doctoral Dissertation Haoran Meng

78

Figure 3.9. (a) and (b) Zoom-in views of two 100 s windows highlighted by the grey
horizontal bars in Figs 8a and b. (c) and (d) Corresponding spectrograms of (a) and (b). The
RN, NRN and the mixture waveforms are colored in red, black and blue.

3.3 Anatomy of recorded waveforms
The waveforms and corresponding spectrograms of two example hours recorded by
geophone R3010 on Julian day 146 from 14:00 to 15:00 (local time) colored, red, black, and
blue for RN, NRN, and mixed signals, respectively, are shown in Fig 3.8. The earthquakes
detected by Meng and Ben-Zion (2018a) are highlighted by the green vertical lines and occur
in windows of NRN. The air traffic events from 200 to 400 s and 1000 to 1200 s also fall in
windows of NRN and have clear Doppler effect in the spectrogram (Fig. 3.8b). The RN
windows have signals with lower amplitudes. The NRN windows have either larger
amplitudes or clear signatures that stand out from the background in the spectrograms.
Doctoral Dissertation Haoran Meng

79
However, the windows with lower amplitudes are not necessarily RN as exemplified by the
time window from 2550 to 2650 s marked with red horizontal bar in Fig. 3.8a and is
identified as NRN due to 3 overtones in the spectrogram. The analyses for nighttime hours
show similar results with more time windows defined as RN. During the daytime RN and
NRN waveforms are found to occupy 25.9% and 34.1% of the time. During nighttime the RN
fraction increases to 41.6% and the NRN fraction decreases to 25.2% (Figs 3.8c and 3.8d). As
an additional demonstration of detecting NRN, Fig. 3.S3 shows the waveforms and
spectrograms of an example hour (Julian day 146 from 5:00 to 6:00 local time), and a zoom-
in view of a M 2.86 local earthquake with an epicentral distance of 176 km recorded by
geophone R3010. The entire earthquake waveform from the first arrival to the end of the
clear coda is identified as NRN in Fig. 3.S3c, substantiating the sensitivity of the developed
methodology.
Fig. 3.9 shows the discussed features in more detail using time windows of only 100 s.
The daytime data show NRN with large amplitudes and strong energy around 150 Hz (Figs
3.9a and b). The periods of RN contain lower spectral amplitudes as expected and are not
continuous for the 100 s duration presented. In general, the energy of NRN spans a large
spectral range and the transition from RN to NRN is smooth. No hard threshold is determined
to differentiate the RN from NRN, and the mixture label is used to classify the superposition
of RN and weak NRN during the transitions.
Many NRN windows show significant energy around 150 Hz (Fig. 3.8) and the signals
possibly originate from the atmosphere or very local sources since high frequency seismic
waves attenuate quickly within distances of 20-30 meters (Meng & Ben-Zion 2018a).
Aircrafts can generate high frequency acoustic waves with little intrinsic attenuation in the
air that couple to the ground and recorded by the seismometers (Meng & Ben-Zion 2018b).
However, the Doppler effect of sources moving in the air is not observed in most of the NRN
windows indicating the vibrations are possibly generated by very local sources on the
ground. The spectra for a geophone located near a tall tree have this high frequency signal
(Fig. 3.S4) corroborating that some such sources originate from objects on the surface. Fig.
3.S4 shows a clear propagation effect from geophone R2005 to its nearby seismometers. The
corresponding spectra at R2005 has a peak around 125 Hz and the high frequency energy
attenuates within 20 meters. The vibrations of the vegetation are likely driven by wind
Doctoral Dissertation Haoran Meng

80
(Johnson et al. 2019) and the peak frequency at 125 Hz may be associated with instrumental
effect or the coupling of the tree next to the sensor to the ground, with possible generation
of small failure events in the top crust. At the SGB site many of the NRN windows with
significant energy around 150 Hz are likely wind generated ground motions if air traffic is
not observed in the spectrogram.

3.4 Discussion
The developed methodology uses waveforms from a single seismic station at a time, to
examine the statistical characteristics of data segments in time and frequency domains, and
separate the records into RN, NRN, and undetermined mixture. Seismic records consist of
RN, which may be associated partly with classical diffuse noise, and NRN containing
earthquakes and ground motions produced by atmospheric and anthropogenic activities
such as air traffic, lightning, cars, trains, and wind shaking of structures, vegetation and other
obstacles above the surface. Small earthquakes can be detected by multiple techniques such
as template matching or beamforming. Unlike tectonic activities, NRN produced by
atmospheric and anthropogenic activities can depend strongly on various aspects of the site
including the distribution and types of obstacles above the ground and attenuation in the
very shallow crust. NRN is therefore not a distinct repeating signal, making it unfavorable
for similarity-based detection techniques.
The presented waveform classification method is stable and not sensitive to the initial
criteria of building and updating the libraries 3 classes of waveforms. The method jointly
examines the amplitude spectra of moving time windows and the MACCs with RN and NRN
templates from hourly waveforms. The hourly RN is quasi-stationary (Fig. 3.S1) and the
results group in the parameter space of the median MACC and L2 norm deviations from the
mean spectra of random noise candidates (Fig. 3.7). Time windows with parameters in this
tight cluster are identified as RN, windows that deviate significantly from the RN cluster are
identified as NRN, and windows with values in between RN and NRN are identified as
undetermined mixtures. Several iterations are performed on each hourly data to update the
random noise templates and mean noise spectra until the libraries stabilize. We tested an
initial RN library by including the indexes of 500 windows with the smallest RMS without
excluding outliers and the results quickly converge in 4 iterations and are almost identical
Doctoral Dissertation Haoran Meng

81
to those shown in Fig. 3.8. The results are not sensitive to the criteria in equation 3.7 and
perturbing the lower and upper thresholds by 0.1 produces changes <1% in the fraction of
time covered by RN and NRN. The results are also not sensitive to the domain size used to
compute the data points densities in Fig. 3.7, which is 0.2 ∙ STD along each axis. Variations of
this choice from 0.1 to 0.4 ∙ STD produce no more than 0.6% changes in the fraction of time
covered by RN and NRN.
By detecting RN windows for each hour we are able to estimate the hourly mean spectral
amplitude of the RN (Fig. 3.10). The amplitude spectra change gradually from one hour to
the next, supporting the notion that the RN is quasi-stationary. Since there are more
atmospheric and anthropogenic activities in daytime hours, the energy of scattered
wavefield is also higher. Peaks at 20 Hz, 80 Hz, and 150 Hz are consistently seen in the
spectrograms, indicating various local resonators including vegetation, ground structures
and sedimentary layers. The hourly mean spectral peaks at 125 Hz in Fig. 3.S5 and 110 Hz
peak in Fig. 3.S6 may be associated, respectively, with instrumental effects and/or coupling
of the tree and antenna tower near sensors to the ground. Fig. 3.11 summarizes daily records
by the hourly percentage of RN, NRN and the mixture for (a) Julian day 146 for all analyzed
geophones and (b) all available days for a single geophone. The curves in Fig. 3.11b are
smoother and the uncertainties are slightly larger than those in Fig. 3.11a, but the overall
structures are consistent. The percentage of the day occupied by NRN decreases significantly
from ~ 50% in daytime to ~30% at night. The results are consistent with the occurrence of
the air-traffic events (Meng & Ben-Zion, 2018b) and wind-driven vibration of vegetation and
local structures with diurnal fluctuations that increase at the SGB site by mid-day local time
and subside during the evening hours (Johnson et al., 2019).
It is important to note that since ~35% of the day is covered by RN on average, 7% of the
day is occupied by air-traffic events and 6% by wind-related events, we can now classify ~48%
the continuous waveforms at the SGB site. By excluding NRN and the mixture windows, the
RN windows can be used to stably estimate better and use for imaging and monitoring
studies noise spectra with windows excluding NRN signals. The method can be generalized
to 3-component records by redefining 𝐶 𝑖 MDN
, 𝐶 𝑖 STD
and spectral deviation, and using the
results of the different components for deriving refined structural information (e.g., via H/V
Doctoral Dissertation Haoran Meng

82
ratios). Estimating the extent to which wavefields in different time windows are diffuse
require additional analyses (e.g., Margerin et al. 2009; Sens-Schonfelder et al. 2015; Liu & Ben-
Zion 2016). This topic will be addressed in a follow up study in Chapter 4. The NRN covers
on average 40% of the day and analyzing these signals further will provide improved
understating of common sources producing the observed ongoing ground motion. The RN
windows detected by our method can also provide a good training dataset for studies using
machine learning for phase picking and signal classification.

Figure 3.10. Hourly mean spectral amplitudes of the detected RN of the waveform
recorded by sensor R3010 during Julian day 146, 2014, with colors denoting the hours. (a)
Results for 0:00 to 12:00 GMT (17:00 to 5:00 local time). (b) Results for 12:00 to 24:00 GMT
(5:00 to 17:00 local time).

Doctoral Dissertation Haoran Meng

83

Figure 3.11. (a) Hourly percentage of RN, NRN and mixed signals for all analyzed
geophones in Fig. 1b for Julian day 146, 2014. (b) Hourly percentage of RN, NRN and mixed
signals for all available days at geophone R3010.

Figure 3.S1. RMS amplitude of all detected RN windows at station R3010 in Julian day
146, 2014. The hourly RN is quasi-stationary.

Doctoral Dissertation Haoran Meng

84

Figure 3.S2. Mean amplitude spectra of the RN templates in different iterations for the
two example hours shown in Figs 3.8a and c. The spectra converge quickly in 4 iterations.
Doctoral Dissertation Haoran Meng

85

Figure 3.S3. (a) Example raw waveform recorded by sensor R3010 in Julian day 146,
2014 from 5:00 to 6:00 (local time). The RN, NRN and the mixture windows are colored in
red, black and blue. (b) Corresponding spectrogram of waveform in (a). The colors denote
power (dB). (c) and (d) Zoom-in views of two 200 s windows highlighted by the grey
horizontal bars in (a) and (b). Green vertical lines in (a) and (c) highlight the small
earthquakes detected by Meng & Ben-Zion (2018b).
Doctoral Dissertation Haoran Meng

86

Figure 3.S4. Waveforms generated by vibration of a tree next to station R2005. (a) The
self-normalized waveforms recorded by R2005 and 10 nearby stations in a linear profile
perpendicular to the fault. The sensor spacing is about 10 m. (b) Corresponding self-
normalized amplitude spectra.

Figure 3.S4. Hourly mean spectral amplitude of detected RN in the waveform recorded
by geophone R2005 underneath a tree in Julian day 146, 2014, with colors denoting the
hours. (a) Results for 0:00 to 12:00 GMT (17:00 to 5:00 local time). (b) Results for 12:00 to
24:00 GMT (5:00 to 17:00 local time).
Doctoral Dissertation Haoran Meng

87

Figure 3.S5. Hourly mean spectral amplitude of the detected RN of the waveform
recorded by R0506 beneath an antenna tower for each hour on Julian day 146, 2014, with
colors denoting the hours. (a) Results for 0:00 to 12:00 GMT (17:00 to 5:00 local time). (b)
Results for 12:00 to 24:00 GMT (5:00 to 17:00 local time.

Doctoral Dissertation Haoran Meng

88
4. Evaluation of diffuse waveforms in dense array data near Anza California
(Meng & Ben-Zion, 2019, in preparation)

4.1 Introduction
This work is motivated by a previous study of Liu and Ben-Zion (2016), in which they
develop a method to estimate correlation coefficients of neighboring frequencies in ambient
seismic noise. From the assumption of fully diffuse wave, they derive that different frequency
components of diffuse wave are not correlated. This is a necessary condition of diffuse wave.
To find diffuse wave and make the theory complete, we derive three conditions which are
both sufficient and necessary by examining the expectation of waveform spectra and the
outer products with themselves and their conjugates.

Figure 4.1. (a) Regional map showing the dense array on the San Jacinto fault zone (red
square) with the regional seismic stations (black triangles) in Southern California. The red
triangles show the example station SNO and borehole station GVDA05 at 150m depth. Grey
triangles denote other analyzed broadband stations. The beach ball displays the hypocenter
and focal mechanism of the 2010 Mw 7.2 El Mayor Cucapah earthquake. (b) Zoom in view of
the dense deployment. Colored circles show the median absolute amplitude of waveforms
recorded at the array stations on Julian day 146, 2014. The analyzed geophones are
displayed as the squares. The background grey colors represent topography.

Doctoral Dissertation Haoran Meng

89
In the present work we analyze seismic waveforms recorded by a dense seismic array
with 1108 vertical ZLand geophones (10 Hz) centered on the Clark branch of the San Jacinto
fault zone (SJFZ) at the Sage Brush Flat (SGB) site southeast of Anza, California (red box Fig.
4.1a). The array recorded continuously at 500 Hz from 7 May 2014 to 13 June 2014, and
covered an area of about 600 m × 600 m (Fig. 4.1c) with nominal sensor spacing of 10 m
normal to the fault and 30 m along-strike (Ben-Zion et al., 2015). In addition to data of the
dense ZLand array, we examine waveforms recorded at 200 Hz by a borehole senor ~150 m
deep at station GVDA05 (red triangle in Fig. 4.1a) to evaluate the coda wave of the 2010 Mw
7.2 El Mayor Cucapah earthquake.

4.2 Method
A wavefield can be expressed analytically as a superposition of normal modes,

𝜓 ( 𝒓 , 𝑡 ) = ∑ 𝑎 𝑛 𝑢 ( 𝑛 )
( 𝒓 ) 𝑒 −𝑖 𝜔 𝑛 𝑡 𝑛 ,
(4.1)
where 𝑢 ( 𝑛 )
( 𝒓 ) are normal mode eigenfunctions, 𝑎 𝑛 are complex modal amplitudes, 𝒓 is
location vector, 𝑡 is time and 𝜔 𝑛 is angular frequency corresponding to the nth mode. For a
fully diffuse wavefield (Weaver & Lobkis, 2004), the statistics of the modal amplitudes
require
E[𝑎 𝑛 ] = 0, (4.2a)
E[𝑎 𝑛 𝑎 𝑚 ] = 0, (4.2b)
E[𝑎 𝑛 𝑎 𝑚 ∗
] = 𝐹 ( 𝜔 𝑛 ) 𝛿 𝑚𝑛
. (4.2c)
Since the eigenfunctions usually are not available, it is not feasible to evaluate whether
the waveform is diffused or not using equations 4.2a-c. The equivalent expressions in
frequency domain, which are much easier to observe, are derived in the following of this
section.

Doctoral Dissertation Haoran Meng

90

Figure 4.2. (a) Amplitude spectrum of a 500 s waveform at night. (b) Synthetic diffuse
noise using the amplitude spectrum in (a). (c) The spectrogram of the synthetic diffuse noise.
(d-f) The vector A, matrices B and C of the synthetic diffuse noise in (b). All elements in 𝐴
and B are close to zero. The condition numbers of 𝐵 + 𝐼 and 𝐶 are close to one.

Doctoral Dissertation Haoran Meng

91
According to Liu and Ben-Zion (2016), a truncated recording in frequency domain

𝜓 𝐹 ( 𝒓 , 𝜔 ) = ℱ (rect ( 𝑡 /𝑇 )∑ 𝑎 𝑛 𝑢 ( 𝑛 )
( 𝒓 ) 𝑒 −𝑖 𝜔 𝑛 𝑡 𝑛 )
= ∑ 𝑎 𝑛 𝑢 ( 𝑛 )
( 𝒓 ) 𝑇 sinc[
( 𝜔 − 𝜔 𝑛 ) 𝑇 2𝜋 ]
𝑛 ,
(4.3)
where 𝑇 is window length. The expectation of the record in frequency domain can be
estimated by many truncated windows,

E[𝜓 𝐹 ( 𝒓 , 𝜔 ) ] = ∑ E[𝑎 𝑛 ]𝑢 ( 𝑛 )
( 𝒓 ) 𝑇 sinc[
( 𝜔 − 𝜔 𝑛 ) 𝑇 2𝜋 ]
𝑛 .
(4.4)
Obviously, if we adopt E[𝑎 𝑛 ] = 0, then E[𝜓 𝐹 ( 𝒓 , 𝜔 ) ] = 0. In addition, since the normal
mode eigenfunctions 𝑢 ( 𝑛 )
( 𝒓 ) are not always zero, if we observe E[𝜓 𝐹 ( 𝒓 , 𝜔 ) ] = 0 at all
frequencies, then E[𝑎 𝑛 ] = 0. Therefore, E[𝑎 𝑛 ] = 0 is equivalent to E[𝜓 𝐹 ( 𝒓 , 𝜔 ) ] = 0. It is
easier to evaluate a unitless vector 𝐴 ( 𝒓 , 𝜔 𝑛 ) using
𝐴 ( 𝒓 , 𝜔 𝑛 ) =
|E[𝜓 𝐹 ( 𝒓 , 𝜔 𝑛 ) ]|
2
E[|𝜓 𝐹 ( 𝒓 , 𝜔 𝑛 ) |
2
]
= 0. (4.5)
To derive an equivalent of equation 4.2b, we first compute the covariance of 𝜓 𝐹 ( 𝒓 , 𝜔 𝑚 ) and
𝜓 𝐹 ( 𝒓 , 𝜔 𝑚 )

E[𝜓 𝐹 ( 𝒓 , 𝜔 𝑚 ) 𝜓 𝐹 ( 𝒓 , 𝜔 𝑚 ) ]
= E [∑ ∑ 𝑎 𝑝 𝑎 𝑞 𝑢 ( 𝑝 )
( 𝒓 ) 𝑢 ( 𝑞 )
( 𝒓 ) 𝑇 2
sinc[
( 𝜔 𝑛 − 𝜔 𝑝 ) 𝑇 2𝜋 ] sinc[
( 𝜔 𝑚 − 𝜔 𝑞 ) 𝑇 2𝜋 ]
𝑞 𝑝 ].
(4.6)
If 𝑇 is long enough, sinc[
( 𝜔 𝑛 −𝜔 𝑝 ) 𝑇 2𝜋 ] can be well approximated by a Kronecker delta 𝛿 𝑛𝑝
. Then

E[𝜓 𝐹 ( 𝒓 , 𝜔 𝑚 ) 𝜓 𝐹 ( 𝒓 , 𝜔 𝑛 ) ]
= ∑ ∑ E[𝑎 𝑝 𝑎 𝑞 ]𝑢 ( 𝑝 )
( 𝒓 ) 𝑢 ( 𝑞 )
( 𝒓 ) 𝑇 2
𝛿 𝑛𝑝
𝛿 𝑚𝑞
𝑞 𝑝
= E[𝑎 𝑚 𝑎 𝑛 ]𝑢 ( 𝑚 )
( 𝒓 ) 𝑢 ( 𝑛 )
( 𝒓 ) 𝑇 2

(4.7)
Similarly, E[𝑎 𝑚 𝑎 𝑛 ] = 0 is equivalent to E[𝜓 𝐹 ( 𝒓 , 𝜔 𝑚 ) 𝜓 𝐹 ( 𝒓 , 𝜔 𝑛 ) ] = 0. A unitless form
B( 𝒓 , 𝜔 𝑚 , 𝜔 𝑛 ) is
Doctoral Dissertation Haoran Meng

92
𝐵 ( 𝒓 , 𝜔 𝑚 , 𝜔 𝑛 ) =
|E[𝜓 𝐹 ( 𝒓 , 𝜔 𝑚 ) 𝜓 𝐹 ( 𝒓 , 𝜔 𝑛 ) ]|
2
E[|𝜓 𝐹 ( 𝒓 , 𝜔 𝑚 ) |
2
]E[|𝜓 𝐹 ( 𝒓 , 𝜔 𝑛 ) |
2
]
= 0. (4.8)
𝐵 is a symmetric matrix. In addition, the covariance of 𝜓 𝐹 ( 𝒓 , 𝜔 𝑚 ) and 𝜓 𝐹 ∗
( 𝒓 , 𝜔 𝑛 )

E[𝜓 𝐹 ( 𝒓 , 𝜔 𝑚 ) 𝜓 𝐹 ∗
( 𝒓 , 𝜔 𝑛 ) ]
= ∑ ∑ E[𝑎 𝑝 𝑎 𝑞 ∗
]𝑢 ( 𝑝 )
( 𝒓 ) 𝑢 ( 𝑞 ) ∗
( 𝒓 ) 𝑇 2
sinc[
( 𝜔 𝑛 − 𝜔 𝑝 ) 𝑇 2𝜋 ] sinc[
( 𝜔 𝑚 − 𝜔 𝑞 ) 𝑇 2𝜋 ]
𝑞 𝑝
= ∑ ∑ E[𝑎 𝑝 𝑎 𝑞 ∗
]𝑢 ( 𝑝 )
( 𝒓 ) 𝑢 ( 𝑞 ) ∗
( 𝒓 ) 𝑇 2
𝛿 𝑛𝑝
𝛿 𝑚𝑞
𝑞 𝑝
= E[𝑎 𝑚 𝑎 𝑛 ∗
]𝑢 ( 𝑚 )
( 𝒓 ) 𝑢 ( 𝑛 ) ∗
( 𝒓 ) 𝑇 2
.
(4.9)
If 𝑝 ≠ 𝑞 , E[𝑎 𝑝 𝑎 𝑞 ∗
] = 0 is equivalent to E[𝜓 𝐹 ( 𝒓 , 𝜔 𝑝 ) 𝜓 𝐹 ∗
( 𝒓 , 𝜔 𝑞 ) ] = 0. If 𝑝 = 𝑞 ,

E[𝜓 𝐹 ( 𝒓 , 𝜔 𝑝 ) 𝜓 𝐹 ∗
( 𝒓 , 𝜔 𝑞 ) ]
= E[𝑎 𝑝 𝑎 𝑝 ∗
]𝑢 ( 𝑝 )
( 𝒓 ) 𝑢 ( 𝑝 ) ∗
( 𝒓 ) 𝑇 2
.
(4.10)
Let E[𝜓 𝐹 ( 𝒓 , 𝜔 𝑝 ) 𝜓 𝐹 ∗
( 𝒓 , 𝜔 𝑝 ) ] = 𝐹 ( 𝜔 𝑝 ) 𝑢 ( 𝑝 )
( 𝒓 ) 𝑢 ( 𝑝 ) ∗
( 𝒓 ) 𝑇 2
. Then E[𝑎 𝑛 𝑎 𝑚 ∗
] = 𝐹 ( 𝜔 𝑛 ) 𝛿 𝑚𝑛
is
equivalent to E[𝜓 𝐹 ( 𝒓 , 𝜔 𝑚 ) 𝜓 𝐹 ∗
( 𝒓 , 𝜔 𝑛 ) ] = 𝐹 ( 𝜔 𝑚 ) 𝑢 ( 𝑚 )
( 𝒓 ) 𝑢 ( 𝑚 ) ∗
( 𝒓 ) 𝑇 2
𝛿 𝑚𝑛
. A unitless form
𝐶 ( 𝒓 , 𝜔 𝑚 , 𝜔 𝑛 ) is
𝐶 ( 𝒓 , 𝜔 𝑚 , 𝜔 𝑛 ) =
|E[𝜓 𝐹 ( 𝒓 , 𝜔 𝑚 ) 𝜓 𝐹 ∗
( 𝒓 , 𝜔 𝑛 ) ]|
2
E[|𝜓 𝐹 ( 𝒓 , 𝜔 𝑚 ) |
2
]E[|𝜓 𝐹 ( 𝒓 , 𝜔 𝑚 ) |
2
]
= 𝛿 𝑚𝑛
(4.11)
where 𝐶 is a symmetric definitive matrix. To sum up, equations 4.2a-c are equivalent to
equations 4.5, 4.8, and 4.11. These are sufficient and necessary conditions for diffuse wave.
On one hand, if we observe that 𝐴 is a zero vector, 𝐵 is a zero matrix and 𝐶 is an identity
matrix, the target waveform is fully diffuse. On the other hand, if one or more of these three
conditions are not satisfied, the target waveform is not diffuse.
To demonstrate the method, we first generate synthetic diffuse noise by Fourier
transforming a low amplitude 500 s waveform recorded at night on example sensor R3010,
replacing the phase information with a uniform distribution of values from 0 to 2π, and
performing an inverse transform on the amplitude spectrum and randomized phase (Figs
4.2a and b). The spectrogram is displayed in Fig. 4.1c. The waveforms are cut into 𝑇 = 1 s
non-overlapping windows. Each window is cosine-tapered with a width of 5% on each end
to stably estimate the spectrum using Fourier transform. All spectra are complex vectors.
Doctoral Dissertation Haoran Meng

93
The spectra and the outer products with themselves and their conjugates are averaged to
compute the vector 𝐴 , matrices 𝐵 and 𝐶 in equations 5, 8 and 11. The corresponding results
are shown in Figs 4.1d-f. As each element in 𝐴 is smaller than 0.03 (Fig. 4.1d), vector 𝐴 is
almost identical to a zero vector. To evaluate how well 𝐵 and 𝐶 are similar to a zero and an
identity matrix 𝐼 , respectively, we compute the condition numbers of 𝐵 + 𝐼 and 𝐶 . In the
following of this paper, we use cond ( 𝑋 ) to represent the condition number of a definitive
and symmetric matrix 𝑋 . Matrix 𝐵 and 𝐶 of the synthetic diffuse noise are displayed in Figs
4.2e and f. The corresponding cond ( 𝐵 + 𝐼 ) and cond ( 𝐶 ) are very close to one. Consequently,
the eigenvalues of 𝐵 + 𝐼 and 𝐶 are very close to one. Therefore, 𝐵 and 𝐶 are almost identical
to a zero and an identity matrix, respectively. These indicate the synthetic waveform in Fig.
4.2a is fully diffuse.
As an example of waveform which is not diffuse, we perform the similar analysis to the
record of a local M 2.29 earthquake recorded by the example sensor R3010. Figs 4.3a-c show
the waveforms, corresponding amplitude spectrum, spectrogram, respectively. Vector 𝐴 in
Fig. 4.3d is not zero vector. In addition, the condition numbers of matrices 𝐵 + 𝐼 and 𝐶 are
very large, indicating the record of the earthquake is not diffuse. Figs 4.4, 4.5 and 4.6 show
the example of an air traffic event, a car event and wind gust. The corresponding large
condition numbers of 𝐵 + 𝐼 and 𝐶 imply these waveforms are not diffuse.

4.3 Results and discussion
This paper present a follow up study of Meng et al. (2019) in which they develop a
methodology to separate continuous seismic waveforms into random noise (RN), not-
random-noise (NRN), and undetermined mixtures of signals. The waveform anatomy
method includes analyses of cross-correlation of waveform segments and amplitude spectra
comparisons. The hourly RN is quasi-stationary and the results cluster tightly in the
parameter space of cross-correlation coefficients and L2 norm deviations from the mean
spectra of RN candidates. To estimate the extent to which wavefields in different time
windows are diffuse, we build on the waveform anatomy method and analyze the data
recorded by the dense deployment.

Doctoral Dissertation Haoran Meng

94

Figure 4.3. (a) The waveform of a local M 2.29 earthquake recorded by the geophone
R3010 in the center of the dense array. The entire waveform is 100 s long. Here we only show
the first 40 s. (b) The amplitude spectrum of the waveform in (a). (c) The spectrogram
waveform in (a). (d-f) The vector A, matrices B and C of the waveform in (a). Elements in 𝐴
are not negligible. The condition number of 𝐵 + 𝐼 and 𝐶 are very large.

Doctoral Dissertation Haoran Meng

95

Figure 4.4. (a) The waveform of an air-traffic event recorded by the geophone R3010 in
the center of the dense array. (b) The amplitude spectrum of the waveform in (a). (c) The
spectrogram waveform in (a). (d-f) The vector A, matrices B and C of the waveform in (a).
Elements in 𝐴 are very close to zero but the condition number of 𝐵 + 𝐼 and 𝐶 are
significantly greater than one.

Doctoral Dissertation Haoran Meng

96

Figure 4.5. (a) The waveform of a car event recorded by the geophone R3010 in the center
of the dense array. (b) The amplitude spectrum of the waveform in (a). (c) The spectrogram
waveform in (a). (d-f) The vector A, matrices B and C of the waveform in (a). Elements in 𝐴
are very close to zero but the condition number of 𝐵 + 𝐼 and 𝐶 are very large indicating non-
diffuse waveform.

Doctoral Dissertation Haoran Meng

97

Figure 4.6. (a) The waveform of a series of wind gust events recorded by the geophone
R2005 in the dense array. This geophone is right beneath a tree which generate ground
motion when shaken by wind. (b) The amplitude spectrum of the waveform in (a). (c) The
spectrogram waveform in (a). (d-f) The vector A, matrices B and C of the waveform in (a).
Elements in 𝐴 are very close to zero but the condition number of 𝐵 + 𝐼 and 𝐶 are very large
indicating non-diffuse waveform.

Doctoral Dissertation Haoran Meng

98

Figure 4.7. (a) Example waveform recorded by R3010 on Julian day 146, 2014 from 11:00
to 12:00 (local time). (b-d)The RN, NRN and the mixture windows are colored in red, black
and blue. (e) Corresponding spectrogram of waveform in (a). The colors denote power (dB).

Doctoral Dissertation Haoran Meng

99

Figure 4.8. (a) Example waveform recorded by R3010 on Julian day 146, 2014 from 02:00
to 03:00 (local time). (b-d)The RN, NRN and the mixture windows are colored in red, black
and blue. (e) Corresponding spectrogram of waveform in (a).

Doctoral Dissertation Haoran Meng

100

Figure 4.9. (a) Vectors 𝐴 for waveform segments in Figure 4.7b-d. (b) The condition
numbers of 𝐵 + 𝐼 and 𝐶 for waveform segments in Figure 4.7b-d. (c) Vectors 𝐴 for
waveform segments in Figure 4.8b-d. (d) The condition numbers of 𝐵 + 𝐼 and 𝐶 for
waveform segments in Figure 4.8b-d.
Doctoral Dissertation Haoran Meng

101

Figure 4.10. (a) Raw waveform of the 2010 Mw 7.2 El Mayor Cucapah earthquake
recorded by the station SNO in vertical component. (b) Same waveform but high pass filtered
at 1 Hz. The blue and red lines (in log scale) display the Ratio of L2 norm spectra of different
windows in the frequency bands of 1-50 Hz and 20-50 Hz. The Red arrow indicates the
possible aftershocks. (c) A corresponding spectrogram of the waveform. (d) Raw waveform
of the 2010 Mw 7.2 El Mayor Cucapah earthquake recorded by the borehole station GVDA05.
(e) Same waveform but high pass filtered at 1 Hz. The blue and red lines (in log scale) display
the Ratio of L2 norm spectra of different windows in the frequency bands of 1-100 Hz and
20-100 Hz. The Red arrow indicates the possible aftershocks. The grey arrows indicate the
possible instrumental glitches. (f) A corresponding spectrogram of the waveform.

Doctoral Dissertation Haoran Meng

102

Figure 4.11. (a, b) Vectors 𝐴 and condition numbers of 𝐵 + 𝐼 and 𝐶 for the waveforms of
noise, body wave, Rayleigh wave and coda in Figure 10a. (c, d) Vectors 𝐴 and condition
numbers of 𝐵 + 𝐼 and 𝐶 for the waveforms of noise, body wave, Rayleigh wave and coda in
Figure 10d. The condition numbers of coda drop significantly after cleaning the aftershocks
and other events.

Figs 4.7 and 4.8 presents the raw waveforms, corresponding detected RN, MIX, NRN and
spectrograms of two example hours recorded by geophone R3010 on Julian day 146 from
11:00 to 12:00 and 2:00 to 3:00 (local time). On the two pieces hourly waveforms we perform
classification analysis using the method developed by Meng and Ben-Zion (2019) in section
3. We then perform the same diffuse analysis to the RN, MIX and NRN of the two pieces of
continues waveforms. Vector 𝐴 , matrices 𝐵 and 𝐶 as well as cond ( 𝐵 + 𝐼 ) and cond ( 𝐶 ) are
estimated. These results are shown in Fig. 4.9. RN and MIX are well diffuse but NRN are not
diffuse Using a threshold of 5 as the condition number to differentiate diffuse and not-diffuse
waveforms. In addition, the spikes at 125 Hz in Fig. 4.9c indicate the possible electrical noise
of the 3C Zland sensor which are consistently shown for waveforms at night. These spikes
disappear when the RN level rise in the day.
Earthquake coda, assumed to be fully diffuse, are widely used to extract the empirical
Green’s function of surface wave. As an application of the method demonstrated in section
4.2, we perform the diffuse analysis on the coda of 2010 Mw 7.2 El Mayor Cucapah
Doctoral Dissertation Haoran Meng

103
earthquake at various stations. Fig. 10 shows the raw and high passed waveforms recorded
by strong motion sensors, and corresponding spectrograms in Z components at station SNO
on surface and GVDA05 at 150 m depth. We separate the waveforms in to noise, P and S
waves, Rayleigh wave using the predicted arrival based on isap91 velocity model. The onsets
of coda are 200 s later than the arrivals of Rayleigh wave. We then perform the diffuse
analysis to each segments of the high-pass filtered waveforms. The vector 𝐴 , cond ( 𝐵 + 𝐼 )
and cond ( 𝐶 ) indicate the noise before P arrival are fully diffuse but the body and Rayleigh
waves are not diffuse. Fig. 4.10f shows the spectrogram of the raw waveform of 2010 Mw 7.2
El Mayor Cucapah earthquake recorded by station GVDA05 with possible instrumental
glitches. Therefore we perform the diffuse analysis both on the cleaned coda (events
removed) and not-clean coda. The coda at SNO is diffuse while the not clean coda is not
diffuse at station GVDA05. After excluding the events, the coda is fully diffuse.
We perform the same analysis to the waveforms of 2010 Mw 7.2 El Mayor Cucapah
earthquake recorded by strong motion sensors at 56 stations in Southern California. Fig. 12
shows that 33 out of 56 stations have diffuse coda (not cleaned) and this number increases
to 48 after excluding the events in coda. The results indicate that the developed method is
powerful in evaluating whether the coda is diffuse or not. We strongly recommend
performing this analysis before using coda to extract empirical Green’s functions.
Doctoral Dissertation Haoran Meng

104

Figure 4.12. (a) Condition numbers of 𝐵 + 𝐼 and 𝐶 for the coda of the 2010 Mw 7.2 El
Mayor Cucapah earthquake recorded by all analyzed stations (grey triangles in Fig. 4.1a). (b)
Condition numbers of 𝐵 + 𝐼 and 𝐶 for the coda after cleaning the aftershocks and other
events.

Doctoral Dissertation Haoran Meng

105
5. Semi-automated estimates of directivity and related source properties of small
to moderate Southern California earthquakes with second seismic moments
(Meng et al., 2019, in preparation)

5.1 Introduction
Earthquake catalogs (e.g. Hauksson et al. 2012) typically represent events as point sources
with location, time, magnitude and focal mechanism. The “cut and paste” technique was
developed to estimate the centroid location, time, focal mechanism using broadband
regional seismograms (e.g. Zhao & Helmberger, 1993; Zhu & Helmberger, 1996). The centroid
location and time are the spatial and temporal mean (first degree moment) of the earthquake
moment release distribution and routinely cataloged from both regional and global seismic
networks. In addition to the centroid, second-degree seismic moments (variances) of the
moment-release distribution is the simplest general representation of source finiteness
including as rupture extent, duration and directivity (e.g. Backus & Mulcahy, 1976a, b; Backus,
1977a, b; McGuire, 2004; McGuire, 2017). Descriptions of source finiteness however, are not
yet included in the standard catalogs due to a lack of techniques for estimating the second
seismic moments with little user involvement.
A seismogram is a convolution of instrument response, seismic wave propagation effects
and earthquake source terms. To estimate the source properties such directivity and rupture
size, the path and instrument effects should be removed from the data. This is often done by
deconvolving seismograms of the target earthquake with its empirical Green’s function (eGf)
given by seismograms of a suitable small-collocated event (e.g. Mueller, 1985; Hutchings &
Wu, 1990; Hough & Dreger, 1995; Ross & Ben-Zion, 2016). Although eGfs can effectively
account for the path and site effects, it is difficult to find a perfect eGf with collocated
hypocenter and almost identical FM to a target earthquake. The difference of the FMs
between the target event and its eGf can result in a failure to obtain accurate ASTFs for
stations in multiple directions and therefore poor station coverage and biased estimation of
the source proprieties. In addition, Calderoni et al. (2012, 2015) showed that if eGf events
are large enough to contain directivity effects in the frequency band of interest, that can
result in apparent directivity of the target event. Therefore, using a single eGf may slightly
bias the results if the eGf has non-generic source effects such as directivity. In this paper, we
Doctoral Dissertation Haoran Meng

106
develop a semi-automated technique using a weighted stack of eGfs to get a better
represented eGf to improve the FM and suppress the non-generic source effects for deriving
source properties of small to moderate earthquakes by the analysis of the second-degree
seismic moments.

Figure 5.1. A map of the Southern California plate boundary region. SCSN and ANZA
broadband stations used in this study are denoted by black triangles. The hypocenters of
target events are indicated by red circles. The associated beach balls display the focal
mechanisms and focal depths. Example stations RVR and IDO are denoted by dark blue
triangles. The detailed information are shown in Table 5.1.

Many studies assume that earthquake kinematic properties such as rupture velocity are
generally similar everywhere, which may be true to the first order. However, this assumption
can suppress recognition of persistent variations of source processes related to different
properties of fault zones. For instance, systematic directivity of earthquakes ruptures on
bimaterial faults has been proposed by various dynamic simulation and observations (e.g.,
Weertman, 1980; Ben-Zion & Andrews, 1998; Ampuero & Ben-Zion, 2008; Brietzke et al., 2009).
Doctoral Dissertation Haoran Meng

107
Directivities of large earthquakes can significantly amplify the ground motion and therefore
cause serious damage at sites in the forward direction.
In this paper, we describe a technique of semi-automated analyses of source properties of
individual target events using the second-degree seismic moments. We use a weighted stack
to get a better represented eGf to suppress deficiencies of using single eGfs. We applied our
method using data collected from stations of Southern California Seismic Network (SCSN)
for 14 small to moderate earthquakes (Magnitude 3.5 to 5.2) to estimate the rupture
directivity in each event. The corresponding stress drop estimates including upper and lower
bounds are determined using various constraints on the rupture dimensions in the second
moment inversion.

5.2 data
We present the method using data provided by the Southern California Seismic Network
(SCSN) and Broadband Seismic Data Collection Center (ANZA) networks collected over a set
of 41 well recorded target earthquakes (Table 5.1; Fig. 5.1). Most of these events occurred in
the San Jacinto Fault Zone (SJFZ), the Elsinore Fault Zone (EFZ) or the East California Shear
Zone (ECSZ), and they have similar right lateral strike-slip FMs. For each target earthquake,
tentative eGf candidates are selected from the relocated catalog of Hauksson et al. (2012)
with magnitudes 1.5 to 2.5 units smaller than the target event and a hypocentral separation
of no more than 5 km from the target event. The details of how these eGf candidates are
qualified and processed are documented in the next section.
For the target and eGf events analyzed in this study, we only work with broadband and
strong motion records (HH and HN channels) with sampling rates of 100 Hz and higher. A
bandpass filtering with corner frequencies at 0.5 and 20 Hz is performed to suppress the
noise possibly caused by ocean tides, wind and anthropogenic sources(e.g. Meng & Ben-Zion,
2018b; Inbal et al., 2018; Johnson et al., 2019). We primarily use the transverse component
to get relatively clean SH wave data and vertical component to get P wave data. The
deconvolution results of SH wave in general have significant better performance in
waveform fitting than both P and SV wave. To improve the computational efficiency, we
exclude records with low signal to noise ratio (SNR) and clipped waveforms detected by the
Doctoral Dissertation Haoran Meng

108
technique developed by Yang and Ben-Zion (2010). Broadband records are preferred when
both HH and HN channels are available.

Figure 5.2. A flow diagram summarizing the main steps in the semi-automated method.

Doctoral Dissertation Haoran Meng

109
5.3. Method
5.3.1 Mathematical representation of the second seismic moments
Following McGuire (2004, 2017), we make the simplifying assumption that the spatial
variations in moment rate of an earthquake can be described as:

𝑀 ̇ ( 𝑟 , 𝑡 ) = 𝑀 𝑓 ̇ ( 𝑟 , 𝑡 ) ,
∫ 𝑓 ̇ ( 𝑟 , 𝑡 ) 𝑑𝑉𝑑𝑡 = 1,
(5.8)
where 𝑀 is the seismic moment tensor, 𝑀 ̇ ( 𝑟 , 𝑡 ) is the moment rate function and 𝑓 ̇ ( 𝑟 , 𝑡 ) is a
scalar function that describes the spatial and temporal distribution of moment release along
the fault. The source time function (STF) is defined as an integral of 𝑓 ̇ ( 𝑟 , 𝑡 ) over the entire
source volume.
The first seismic moments of an earthquake are the centroid location 𝑟 0
and time 𝑡 0

representing the spatiotemporal means of its moment release

𝑟 0
= ∬ 𝑓 ̇ ( 𝑟 , 𝑡 ) 𝑟 𝑑𝑉𝑑𝑡 ,
𝑡 0
= ∬ 𝑓 ̇ ( 𝑟 , 𝑡 ) 𝑡 𝑑𝑉𝑑𝑡 .
(5.2)
The second seismic moments are mathematical representation to capture the overall
kinematic properties of a rupture that are well constrained by the far-field seismic
waveforms. The second moments are defined as the second-order space and time moments
of the normalized moment-rate distribution function (McGuire et al., 2001):

𝜇 ̂
( 2,0)
= ∬ 𝑓 ̇ ( 𝑟 , 𝑡 ) ( 𝑟 − 𝑟 0
)
T
( 𝑟 − 𝑟 0
)𝑑𝑉𝑑𝑡 ,
𝜇 ̂
( 0,2)
= ∬ 𝑓 ̇ ( 𝑟 , 𝑡 ) ( 𝑡 − 𝑡 0
) ( 𝑡 − 𝑡 0
)𝑑𝑉𝑑𝑡 ,
𝜇 ̂
( 1,1)
= ∬ 𝑓 ̇ ( 𝑟 , 𝑡 ) ( 𝑟 − 𝑟 0
) ( 𝑡 − 𝑡 0
)𝑑𝑉𝑑𝑡 ,
(5.3)
where 𝑟 0
and 𝑡 0
are column vectors representing the centroid location and time defined in
equation (2). The hat signs in equation (3) indicate these second moments are central
moments taken about the centroid. The integrals are taken over the entire source volume
and earthquake duration (Backus, 1977a,b; McGuire et al., 2001). The second spatial moment
𝜇 ̂
( 2,0)
is a 3 by 3 symmetric tensor which associated with the rupture extent; the temporal
Doctoral Dissertation Haoran Meng

110
moment 𝜇 ̂
( 0,2)
is a scalar associated with the rupture duration; and the mixed moment 𝜇 ̂
( 1,1)

is a column vector associated with rupture propagation (McGuire, 2004, 2017).
The characteristic rupture duration 𝜏 c
, rupture length 𝐿 c
, and centroid rupture velocity 𝑣 0

are defined as

𝑥 c
( 𝑛 ̂)= 2
√
𝑛 ̂
T
𝜇 ̂
( 2,0)
𝑛 ̂,
𝜏 c
= 2√𝜇 ̂
( 0,2)
,
𝑣 0
=
𝜇 ̂
( 1,1)
𝜇 ̂
( 0,2)
,
(3.4)
in which 𝑥 c
( 𝑛 ̂) is the spatial extent of the rupture in the direction 𝑛 ̂. 𝐿 c
= 𝑚𝑎𝑥 ( 𝑥 c
( 𝑛 ̂) ) is the
largest eigenvalue of the special moment and the characteristic rupture width 𝑊 c
is the
second largest eigenvalue. The directivity ratio (0 ≤ 𝑑𝑖𝑟 ≤ 1) is defined as

𝑣 c
=
𝐿 c
𝜏 c
,
𝑑𝑖𝑟 =
|𝑣 0
|
𝑣 c
,
(3.5)
where 𝑣 c
is the characteristic velocity. The directivity ratio is a scaler indicating how strong
the directivity is. 𝑑𝑖𝑟 = 0 for a perfect symmetric bilateral rupture while 𝑑𝑖𝑟 = 1 for a
uniform unilateral rupture.

5.3.2 Qualifying eGf candidates
Based on the approach demonstrated in McGuire (2004, 2017) which requires heavy user
involvement, we develop a semi-automated technique which automates most of the steps.
The methodology is summarized in a flowchart diagram (Fig. 5.2) and illustrated using
seismic data of the 2014, Mw 4.58 target event (No. 5 event in Table 1).
In addition to these data, we use the 2009 to 2018 refined catalog (e.g. Hauksson et al.,
2012) to select the eGf candidates with magnitudes 1.5 to 2.5 smaller than that of the target
and hypocentral separations of no more than 5 km. This is to guarantee the selected eGf
events are well approximated by point sources and collocated with the target event’s
hypocenter to account for the path effects. We don’t limit the selection in a time-range (i.e. a
month before and after the target earthquake) to maximize the number of eGf candidates
Doctoral Dissertation Haoran Meng

111
and increase the chances to find the best eGfs. We then choose seismic stations at
hypocentral distances less than 120 km, thereby avoiding the noisier records from stations
at farther hypocentral distances.

Figure 5.3. An example of PLD for the Mw 4.5.8, 2014 at station RVR. (a) The S waveform
and its fit in T component. (b) The S waveform of the eGf in T component. (c) The black curve
shows the normalized misfit and the grey curve shows its derivative. The red circles denote
the selected duration. (d) The ASTF.

We then average CVM-S4.26 model in Southern California (Lee, et al., 2014) to a 1-D
velocity model to predict the waveforms’ P and S arrivals. These arrivals are improved by
maximizing the Signal to Noise Ratio (SNR) on each zero-crossing of waveforms within the
±1 second window (Li et al., 2016). The SNR is estimated using a ratio of two 1-second
Doctoral Dissertation Haoran Meng

112
moving averages after and before the tentative P/S arrival. The eGf waveforms with SNR ≤ 5
since are discarded since they generally have bad performances.

Figure 5.4. Weighted stack of eGfs for example stations RVR and IDO in Fig. 5.1 for the
Mw 4.58, 2016 All waveforms are S phases in transverse component. (a) Waveforms of target
events and 3 best performed eGfs recorded by station RVR, eGf 1 is the best performed eGf
in deconvolution. (b) Grid search for the coefficients for eGf 2 and eGf 3, bold white circles
denote the best pair of weights. Color displays the normalized misfit in deconvolution. (c)
target waveforms, ASTFs and waveform fits for the weighted stacked eGfs at station RVR. (d),
Doctoral Dissertation Haoran Meng

113
(e) and (f) are same as (a), (b) and (c) but for station IDO. The circles and vertical lines in (a)
and (c) indicate the predicted and refined arrivals.

In the next step, we evaluate the qualities of eGf candidates using cross-correlation and
deconvolution. We first bandpass filter the T component raw waveforms from 0.5 to 2 Hz
and cut 10 s long S phases to compute the maximal cross correlation coefficients (CCs) of the
target earthquake and its eGf candidates (Fig. 5.S1a). These waveforms are only employed
for computing CCs, not for deconvolution. We then roughly evaluate the quality of the eGf
candidates by their average CCs and the numbers of available records (Fig. 5.S1b). To
promote the computing efficiency, we only perform the time-consuming deconvolution
analysis for no more than 10 eGf candidates with higher qualities. Fig. 5.S1c shows the cross-
correlation coefficient functions of the S phases of the target and its eGf (SCSN event ID:
15521073) at stations IDO and RVR. The pulses around the maximum CCs are very similar
to the ASTFs computed from the deconvolution. This is because the phase terms of the
devolution and cross-correlation are identical. Therefore, the CCs are higher at stations with
ASTFs of a single pulse, which are more likely to be the rupture propagation direction. For
the target event 2014, Mw 4.58, the CCs are generally higher at stations in the azimuth range
16°-145° indicating a possible directivity towards the southeast. This is proved by further
analyses in section 3.3.
In the calculation of the deconvolution, the window lengths for the P and S phases are 256
and 512 samples with a sampling rate of 100 Hz respectively. The qualities of the eGf
candidates are then finalized using the Projected Land-Weber Deconvolution (PLD)
algorithm (Betero et al., 1997; Lanza et al., 1999) which performs the deconvolution with
moment release restricted to a series of increasing-length time intervals and analyzes the
misfit as a function of the interval length. Example of the PLD processing step is shown by
Fig. 5.3 with the SH phase of the 2014, Mw 4.58 event and its best-performing single eGf at
station RVR (Fig. 5.1). The deconvoluted Apparent Source Time Function (ASTF) is then
estimated using PLD with non-negative and limited duration constraints. The normalized
misfit (Fig. 5.3c), which is calculated as the L2 norm of the difference of the observed and
fitted SH waves normalized by the L2 norm of the observed SH, decreases monotonically as
longer permissible time intervals allow moment release to occur. The optimized duration of
Doctoral Dissertation Haoran Meng

114
the ASTF is then picked at the interval where the misfit is small enough and the tradeoff
curve is relatively flat. The flat stage indicates a relative stable ASTF estimate that is not
significantly affected by a small perturbation in the permitted duration. The selection of
duration is automated by selecting the smallest duration with normalized misfit ≤ 0.3 and its
time derivative at local maxima (close to zero). The optimized ASTF with respect to the
duration picked in Fig. 5.3c is shown in Fig. 5.3d.
Since the deconvolution is also sensitive to the accuracy of the onsets of the phases
especially for eGfs, we perform grid search on the zero-crossings of the target and eGf
waveforms within a short time window (i.e. ±0.5 second) around the improved P/S picks and
run PLD for each pair of zero-crossings. The P/S picks are then further refined to the zero-
crossings with the best fit of the target waveform. The performance of eGf candidates are
then determined from high to low by the number of stations with acceptable waveform fits.
This step of evaluating each eGf is relatively time-consuming and takes most of the
computing time because of the grid search and the time domain deconvolution but it has
been automated.

5.3.3 Weighted Stack of eGfs
To ensure good station coverage, and to stabilize the second seismic moment inversion
for an target earthquake, we need an eGf that produces as many measurements of apparent
duration as possible. Here we decide to use additional eGf candidates to stack with the best-
performing eGf to obtain a better represented eGf. This is motivated by the fact that a noise-
free eGf is a linear combination of five independent Green’s function components assuming
double couple source mechanism. The combination coefficients depend on the strike, dip and
rake of the rupture. Therefore, a linear combination five collocated and well-aligned eGfs
with different FMs but similar sizes can possibly recover an eGf with any double couple FM.
In this paper, we use two additional eGfs with good performance in the deconvolution to
further improve the best-performing eGf, which usually has the closest FM to the target event.
In addition, the stacking also reduces the noise level and decreases non-generic source
effects, such as directivity of an single eGf.
Doctoral Dissertation Haoran Meng

115
After we find three qualified eGfs for the target event, i.e. 2014, Mw 4.58 earthquake, we
first normalize the eGfs individually by the seismic potency calculated from the potency-
magnitude scaling relation of Ben-Zion and Zhu (2002). The waveforms of these eGfs are
different but basically have enough coherencies to be aligned using cross-correlation at most
stations. Fig. 5.4 shows the SH phases of the target, three aligned and normalized eGfs,
associated grid searches of the stacking weights and corresponding deconvolution results
(Fig. 5.4c and f) of the 2014, Mw 4.58 target event at station RVR and IDO. The eGfs are
aligned by cross-correlations and the predicted and refined arrivals of S wave are denoted
by red circles and vertical lines in red and black. The amplitude differences at various wiggles
are primarily because of the FM differences. Here we use two additional eGfs to further
improve the best-performing eGf to obtain more measurements. We perform a grid-search
to determine the weights of eGfs 2 and 3 with a 0.1 spacing in the range from -3 to 3. The
weight of eGf 1 is fixed as 1 to reduce the search complexity. The duration of ASTF is fixed as
the duration picked in the PLD for the grid search for the best-performing eGf 1 and that
value is used in the PLD to calculate the misfit for each pair of weights. The white cross at the
local minimum of Fig. 5.4b denotes the best pair of weights for the stacking coefficients. Fig.
5.4c shows the weighted stacked eGfs with the selected pair of weights, the resulting ASTFs,
and waveform fits. In this paper, we choose to limit the stack to three eGfs instead of five to
achieve improved performance but retain computation efficiency. Stacking more than three
eGfs will dramatically increase the searching complexity but could be implemented with
much long computing time. The same processes at another example station IDO are shown
by Figs 5.4d-5.4f. The best pairs of weights are different at various stations indicating the
aleatory variability of the measurements. The ASTFs in Figs 4c and 4f are significantly
different because station RVR and IDO are at different directions to the rupture propagation
direction and record the directivity differently.

Doctoral Dissertation Haoran Meng

116

Figure 5.5. Averaged misfit distribution for different pairs of weights for all accepted
stations of (a) the 2014 Mw 4.58 and (b) the 2013 Mw 4.70 target events. The bold white
crosses denote the pair of weights with the best waveform fit.

The misfits as functions of weights of eGf 2 and 3 (i.e. Figs 5.4b and 5.4e) are then stacked
and averaged iteratively to find a pair of weights works for all stations. In the first iteration,
we stack and average the misfits for all stations to obtain a pair of weights with least misfit.
Stations with misfit > 0.3 for the selected weights are then discarded. In the following
iterations, we stack and average the misfits for the remaining stations and find a new pair of
weights with least misfit. We repeat this process until no more stations are excluded. The
iterative process is very stable and usually converges in three iterations. Fig. 5.5 presents the
final averaged misfit patterns for the 2014, Mw 4.58 and the 2013, Mw 4.7 events with the
selected weights denoted by the white crosses. The average misfit is significantly smaller
using the weighted stacked eGf compared with a single eGf.

Doctoral Dissertation Haoran Meng

117

Figure 5.6. (a) Characteristic durations, (b) ASTFs and (c) normalized misfits for the best
performed and weighted stacked eGfs for the 2014, Mw 4.58 event. P phase in vertical and S
phase in transverse directions. Texts display the station names and corresponding back
azimuths for all stations.

The characteristic durations, ASTFs and normalized misfits for the 2014 Mw 4.58 event
at different stations are presented in Figs 5.6a-5.6c. The ASTFs (Fig. 5.6b) show systematic
variations with the back azimuth thus indicate a significant directivity to the SE. We also note
that there are two clear pulses in the ASTFs for the approximate azimuth range 225°-320°,
suggesting the existence of two major sub-events. Fig. 5.6c shows the normalized misfit
before and after using the weighted stack. The stacked eGf results in 36 accepted
measurements (26 by using single eGf) for this target event and thus gives better azimuthal
and take-off angle coverage for the second moment inversion. Figs 5.7a-5.7c show the
corresponding results for a 2013 Mw 4.7 earthquake on the San Jacinto fault zone. The ASTFs
for this event indicate moderate directivity to the NW and again suggest (ASTFs in the
approximate azimuth range 69°-225°) the existence of two main sub-events. These
Doctoral Dissertation Haoran Meng

118
inferences for the 2013 Mw 4.7 are consistent with McGuire (2017). Fig 5.S2 presents the
waveforms and data fit of these two example target events. The good waveform fits indicate
the high quality of weighted stacked eGfs and good ASTF estimates. Fig 5.S3 illustrates the
characteristic durations, ASTFs and normalized misfits of the 2016, Mw 5.19 earthquake in
the Trifurcation area, which shows a strong directivity towards the northwest. This is
consistent with the results of Ross and Ben-Zion (2016) and Ross et al. (2017).

Figure 5.7. (a) Characteristic durations, (b) ASTFs and (c) normalized misfits for the best
performed and weighted stacked eGfs for the 2013, Mw 4.7 event. P phase in vertical and S
phase in transverse directions. Texts display the station names and corresponding back
azimuths for all stations.

The weighted stack technique employed here, compared with previous second moment
methods using only one eGf, results in a larger number of accepted measurements and
enhances their stability of the measurements.
Doctoral Dissertation Haoran Meng

119

5.3.4 Inversion scheme for second moment
The second moments are related to the azimuthal variations in the duration of ASTFs at a
given station as

𝜇 ̂
( 0,2)
( 𝑠 ) = 𝜇 ̂
( 0,2)
− 2𝑠 ∙ 𝜇 ̂
( 1,1)
+ 𝑠 ∙ 𝜇 ̂
( 2,0)
∙ 𝑠 ,
𝜏 c
( 𝑠 ) = 2√𝜇 ̂
( 0,2)
,
(5.6)
where 𝜇 ̂
( 0,2)
( 𝑠 ) is measured at each station and 𝑠 is the slowness of P/S phase in the source
region, 𝜏 c
( 𝑠 ) is the apparent characteristic duration (McGuire 2004, 2017). The slowness
vectors can be easily obtained with the 1-D ray tracing techniques. Employing the FM of the
target earthquake, the number of unknown parameters can be reduced from 10 to 6 by
projecting equation (5.6) onto the 2D fault plane using the rotation matrices of equations
(5.A2-5.A4) in the appendix A. By ignoring the fault perpendicular direction, equation (5.6)
can then be rewritten as

𝜇 ̂
( 0,2)
( 𝑠 ) = [𝑠 1
2
2𝑠 1
𝑠 2
𝑠 2
2
−2𝑠 1
−2𝑠 2
1] ∙ 𝑥
𝑥 = [
𝜇 ̂
11
( 2,0)
𝜇 ̂
12
( 2,0)
𝜇 ̂
22
( 2,0)
𝜇 ̂
1
( 1,1)
𝜇 ̂
2
( 1,1)
𝜇 ̂
( 0,2)
]
T

(5.7)
where 𝑠 1
(along strike) and 𝑠 2
(along dip) are parallel to the fault plane. With observed
durations of ASTFs at various stations, we can build a linear system for inversion

𝑏 = 𝐴 ∙ 𝑥 ,
(5.8)
where 𝑏 is a column vector consists of apparent durations 𝜇 ̂
( 0,2)
( 𝑠 ) and 𝐴 is the matrix
associated with the slowness components in equation (5.7). Equations (5.7) and (5.8) are
subject to
[
𝜇 ̂
( 0,2)
𝜇 ̂
( 1,1) T
𝜇 ̂
( 1,1)
𝜇 ̂
( 2,0)
] ≥ 0, (5.9)
which enforces the physical constraint forbidding a negative source volume. The problem is
solved through convex optimization (Vandenberghe & Boyd, 1996; Boyd & Vandenberghe,
2004; Grant & Boyd 2008, 2014) as

min ‖𝑏 − 𝐴 ∙ 𝑥 ‖
2

(5.10)
Doctoral Dissertation Haoran Meng

120
subject to [
𝜇 ̂
( 0,2)
𝜇 ̂
( 1,1) T
𝜇 ̂
( 1,1)
𝜇 ̂
( 2,0)
] ≥ 0,
and 𝜇 ̂
( 0,2)
≤ 2 ∙ max ( 𝑏 ) .
In addition to the physical constraint, we require second temporal moment less than twice
the largest measurement of 𝜇 ̂
( 0,2)
( 𝑠 ) , which should be satisfied with good station coverage
of a target earthquake.
Fig. 5.8 show the observed and predicted apparent 𝜏 𝑐 s at all stations and the azimuthal
and take-off angle distribution of the accepted stations using the weighted stacked eGf
method for the target events 2014, Mw 4.58 and the 2013, Mw 4.7 event. The apparent 𝜏 𝑐 s
shows strong directivity for both of the target event. To resolve the fault plane ambiguity, we
perform the second moment inversion using both possible rupture planes. The plane with
smaller misfit to the observed the apparent 𝜏 𝑐 s values is determined as the actual rupture
plane. In Fig. 5.8c, there are no significant azimuthal gaps and it contains both up-going (take-
off angle >90°) and down-going (take-off angle <90°) rays. The accepted fault plane for this
event is parallel to SJFZ and SAFZ. The arrow towards the southeast (Fig. 5.8c) indicates the
rupture propagation direction. Table 2 documents the inversion results of the 2014, Mw 4.58
event using both the single eGf and the weighted stack. The directivity ratio is 0.89 indicating
a unilateral rupture. Since we use two more eGfs, here we define the residual per degree of
freedom as the L2 norm of the differences of the observed and the predicted apparent
duration divided by 𝑀 − 3, where 𝑀 is the total number of measurements used for inversion.
The residual per degree of freedom is significantly smaller after using the weighted stack in
Table 5.2 indicating an improvement of the data fit.
The same analyses have been applied to all target earthquakes. The arrows in Fig. 5.9
presents the centroid rupture velocity both in the horizontal and vertical directions. The
colors denote the directivity ratios of analyzed target events.

5.3.5 Uncertainty estimates of rupture area and stress drop
To estimate the stress drop, we follow the approach developed by McGuire and Kaneko
(2018) to evaluate the uncertainty in the rupture dimensions. The upper and lower bounds
on rupture area, permissible for a given data set, are constrained by the 95 percent
Doctoral Dissertation Haoran Meng

121
confidence interval of data misfit based on 𝜒 2
statistics implemented within the convex
optimization. The 𝜒 2
statistic is defined as
𝜒 2
= ∑
( 𝑏 𝑖 − 𝑏 ̂
𝑖 )
2
𝜎 𝑖 2
𝑁 𝑖 (5.11)
where 𝜎 𝑖 2
is the variance of the i-th measurement. In practice, we assume a constant 𝜎 𝑖 2
for
all the measurements, which is calculated from the minimum residual, e.g., the misfit for the
optimal second moment estimates 𝜎 2
. In addition, we assume the degrees of freedom as 𝑀 −
3, where 𝑀 is the number of measurements used for inversion. Correspondingly, the 𝜒 2

statistics will be 𝜒 95%,𝑀 −3
2
. With the degrees of freedom and 𝜎 2
, we can obtain the maximum
permissible misfit as 𝜎 2
∙ 𝜒 95%,𝑀 −3
2
, which is used to constrain the misfit level when
estimating the maximum/minimum rupture area.
Estimating the maximum rupture area at a particular confidence level can be formulated
as:

max det𝜇 ̂
( 2,0)

subject to [
𝜇 ̂
( 0,2)
𝜇 ̂
( 1,1) T
𝜇 ̂
( 1,1)
𝜇 ̂
( 2,0)
] ≥ 0
and [
𝜎 2
∙ 𝜒 95%,𝑀 −3
2
(𝑏 − 𝐴 ∙ 𝑥 )
T
𝑏 − 𝐴 ∙ 𝑥 𝐼 𝑀 ] ≥ 0
and 𝜇 ̂
( 0,2)
≤ 2 ∙ max ( 𝑏 ) .

(5.12)
Solutions of equation (12) set an upper bound on the rupture area (𝐿 𝑐 ⋅ 𝑊 𝑐 ) and hence a
lower bound on stress drop. The posed problem is convex and can be solved with convex
optimization.
Minimizing det𝜇 ̂
( 2,0)
, however, is not a convex problem. We approximate the problem
instead as

min 𝐿 𝑐 2
+ 𝑊 𝑐 2

subject to [
𝜇 ̂
( 0,2)
𝜇 ̂
( 1,1) T
𝜇 ̂
( 1,1)
𝜇 ̂
( 2,0)
] ≥ 0
(5.13)
Doctoral Dissertation Haoran Meng

122
and [
𝜎 2
∙ 𝜒 95%,𝑀 −3
2
(𝑏 − 𝐴 ∙ 𝑥 )
T
𝑏 − 𝐴 ∙ 𝑥 𝐼 𝑀 ] ≥ 0
and 𝜇 ̂
( 0,2)
≤ 2 ∙ max ( 𝑏 ) .
where 𝐿 𝑐 and 𝑊 𝑐 are the maximum and second largest eigenvalues of 𝜇 ̂
( 2,0)
. Solutions of the
problem can roughly estimate the minimum rupture area at a given confidence level
(McGuire & Kaneko, 2018; Fan & McGuire, 2018) therefore give an upper bound on the stress
drop.
The rupture is approximated by an ellipse with area is 𝑆 = 𝜋 𝐿 𝑐 𝑊 𝑐 (McGuire & Kaneko,
2018). Assuming the rupture is an elliptical crack in homogeneous elastic medium (Eshelby,
1957), the stress drop is estimated as
∆𝜎 = 𝐶 ( 𝐿 𝑐 , 𝑊 𝑐 , 𝜈 )
𝑀 0
𝑊 𝑐 𝑆 , (5.14)
where 𝑀 0
is seismic moment, and 𝐶 is a constant of order unity that depends on the shape of
rupture ellipse with major axis 𝐿 𝑐 and minor axis 𝑊 𝑐 and the Poisson ratio 𝜈 of the medium.
When the aspect ratio 𝐿 𝑐 /𝑊 𝑐 is close to one, i.e. < 1.2, 𝐶 ( 𝐿 𝑐 , 𝑊 𝑐 , 𝜈 ) is singular, therefore the
stress drop is estimate by
∆𝜎 = 2.44
𝑀 0
( 𝜋 𝐿 𝑐 𝑊 𝑐 )
3/2
, (5.15)
assuming circular crack. To estimate the stress drops, we first need to convert the
magnitudes to seismic moments. The magnitudes for some target events in the catalog (e.g.,
Hauksson et al., 2012), however, are local magnitude instead of moment magnitude. To avoid
overestimating the stress drops, we first convert the local magnitude to moment magnitude
using empirical relation documented in the catalog change history of SCSN
𝑀 𝑤 = 0.853 ∙ 𝑀 𝐿 + 0.40125 , (5.16)
where 𝑀 𝑤 is moment magnitude and 𝑀 𝐿 is the local magnitude. Table 3 presents the optimal
stress drop estimate of 9.98 MPa for the 2014 Mw 4.58 earthquake and 8.12 MPa for the
2013 Mw 4.7 earthquake. The corresponding lower and upper bounds are (1.12,69.99) MPa
and (4.00,11.11) MPa. The excellent station coverage of these two events produces tight
bounds on stress drop. The stress drop (78.2 MPa) estimated by Ross et al. (2017) using finite
fault model for the 2016, Mw 5.19 earthquake is higher than but in the same order of our
Doctoral Dissertation Haoran Meng

123
optimal estimate 20.9 MPa. Their estimate is in the uncertainty range of our 95% confidence
level (3.1, 177.9) MPa.

Figure 5.8. (a, b) Observed and predicted characteristic durations of the 2014, Mw 4.58
event. The colored dots denote the τc at that station. (c) Azimuthal and take-off angle
distribution of used stations. The black arrow denotes the rupture directivity. (d, e)
Observed and predicted characteristic durations of the 2014, Mw 4.58 event. The colored
dots denote the 𝜏 𝑐 at that station. (f) Azimuthal and take-off angle distribution of used
stations. The black arrow denotes the rupture directivity.

Doctoral Dissertation Haoran Meng

124

Figure 5.9. Directives in horizontal and vertical directions of analyzed target events
shown in Fig. 5.1. The arrows show the centroid rupture velocities. The colors denote the
directivity ratio.

Doctoral Dissertation Haoran Meng

125
5.4. Discussion and conclusions
We describe and implement a semi-automated method for estimating earthquake source
properties from the second seismic moments. The method includes automated steps of
selecting the three best-performing eGfs using cross-correlation and PLD, grid search of the
stacking weights and finally the second moment inversion. Additional manual corrections
are required for some stations (usually a small subset of all available stations) where the
cross-correlation may fail to align the waveforms as the 3 best-performing eGfs can lose good
coherency. The duration picks may also need further adjustment after obtaining ASTFs (i.e.
Figs 5.6, 5.7 and 5.S3) when ASTF of a station is very different from that at nearby stations.
The average computing time of the automated steps for one target event on a CPU node with
2GB memory requires ~ 2 hour. Another 15 min on average is needed for user involvement
to make corrections to the alignments and durations.
The primary goal of this method is to measure a set of earthquake source properties with
little user involvement. Such a technique is necessary for the systematic analyses of large
seismic data sets. For the southern California dataset used in this paper, we perform second
moment inversions using weighted stacked eGfs for target events with measurements of
ASTFs accepted at 15 or more stations that provide good azimuthal coverage. The horizontal
component of the rupture directivities is usually well constrained with good azimuthal
coverage for the strike slip events, i.e., Figs 5.8c and 5.8f. The along dip component of the
rupture directivity can only be well constrained when stations available right on top of the
hypocenter because the seismic rays take off almost horizontally for most other stations. For
events with less good station coverage, we relax our criteria for stations in the azimuthal
gaps or right above the hypocenter (i.e. accepting measurements with normalized misfit <
0.35 instead of 0.3).
A detailed comparison of results of the 2014, Mw 4.58 Event using single eGf and stacked
eGf is documented in Table 5.2. For this event, using the weighted stack results in 10 more
accepted measurements and smaller residuals per degree of freedom, which implies the
improvement of the station coverage and data fit. Table 3 shows the second moment
inversion results for target events of M 3.5 and above. Smaller events (M<3.5) are not
resolvable because the SNR of eGfs are too low to get stable ASTF estimates in the PLD.
Doctoral Dissertation Haoran Meng

126
The results in Table 3 show that the target events with larger magnitudes generally result
in larger estimates of 𝐿 𝑐 and 𝜏 𝑐 as expected. The estimated 𝐿 𝑐 and 𝜏 𝑐 for all target events can
be slightly underestimated because the source time functions of eGfs are not ideal Dirac delta
functions. For example, if the eGf is 2 magnitudes smaller than the target event, the rupture
size and duration of source time function of the eGf is roughly 10% of target event. By
deconvolving the eGf from the target, the duration of the estimated ASTF is about 90% of the
real one. The rupture dimensions are directly associated with the observed characteristic
𝜏 𝑐 s which will be a little underestimated (90% of the true value). Therefore, the centroid
velocity and rupture directivity won’t be significantly affected because these two parameters
are measuring ratio of the rupture dimension and characteristic duration. The directivity
ratio will also not be affected since it is the ratio of the centroid and characteristic rupture
velocities. The stress drops, however, will be slightly overestimated.
Fig. 9 presents the centroid rupture velocities and directivities of all target events. The
directivities of events are consistent with independent results of Kurzon et al. (2014) and
Ross and Ben-Zion (2016) based on ratios of ground motion and source spectra in different
directions. Most events in middle San Jacinto fault zone have directivities towards the
northwest while most events around Cajon Pass and San Gabriel Mountain are towards the
southeast. There are also a few events showing directivities parallel to the conjugate of the
major fault. The results are generally consistent with the prediction of dynamic rupture on a
bimaterial interface (e.g. Andrews & Ben-Zion, 1997; Ampuero & Ben-Zion, 2008) and the
imaged velocity contrast in the study area (e.g. Allam & Ben-Zion, 2012; Zigone et al., 2015;
Share et al., 2016; Shear et al., 2017).
The stress drops and corresponding upper and lower bounds of the target events are
estimated by the approach developed by McGuire and Kaneko (2018). Table 3 and Fig. 10
show that the results span a range from 2 to 145 MPa. The lower bounds of the stress drop,
corresponding to maximum rupture area, are always well constrained. The upper bounds of
a considerable number of events can exceed 10
4
MPa, indicating they are not constrained.
This asymmetry results from the difficulty in ruling out very small rupture widths for some
station geometries. Stress drop estimates are also sensitive to focal depth. Changing the focal
depth will affect the take-off angles connecting the source and stations. The rupture area
Doctoral Dissertation Haoran Meng

127
estimates are well constrained if the events have both up-going (take-off angle > 90°) and
down-going rays (take-off angle < 90°).
Compared with the method of estimating stress drop using corner frequencies (e.g. Brune,
1970; Madariaga, 1976), we still assume the crack model applies (equation 5.14) but
estimate both the major and minor axes of the elliptical rupture area, e.g. 𝐿 𝑐 and 𝑊 𝑐 , instead
of only the major axes. Fig. 5.10b shows a correlation between the stress drop estimates and
the aspect ratio estimates, e.g. 𝐿 𝑐 /𝑊 𝑐 as is expected owing to the term 𝐶 ( 𝐿 𝑐 , 𝑊 𝑐 , 𝜈 ) in
equation (5.14). However, the factor of ~5 variability around this trend indicates that
natural variability in the ratio of rupture length to moment remains important even after
accounting for aspect ratio. Fig. 5.10d shows the ratio of stress drop estimates of an elliptical
crack and a circular crack with the same rupture length 𝐿 𝑐 but different aspect ratios. The
stress drops of the elliptical crack are significantly higher when the aspect ratio > 3. This
indicates that the stress drops are likely to be systematically underestimated using methods
(e.g. corner frequency method) that only constrain one rupture dimension and assume a
circular crack. Our study finds many events have aspect ratios in the 2-5 range. If this pattern
holds in larger datasets, it would imply that many previous studies that have assumed
circular cracks have systematically underestimated stress drop by about a factor of 5-10.
We have focused in this paper on developing a semi-automated technique for estimating
a variety of different source properties with little user involvement (Fig. 5.2). The developed
approach stably resolves the finite source properties of small to moderate earthquakes and
can be used to perform systematic analysis of large data sets. Our results of the rupture
directivity shows consistently with the prediction of dynamic rupture on a bimaterial
interface in southern California. This is important to the seismic hazard analysis in terms of
estimating the ground motion affected the directivity pulse of a potential large earthquake
on the biomaterial interface. By adopting our method to other data set, such as Japan and
China, it can also help to address fundamental long-standing questions such as whether
source properties vary between different regions.

Doctoral Dissertation Haoran Meng

128

Figure 5.10. (a) Stress drops and uncertainties vs. event magnitudes. The vertical arrows
indicate the upper bound of stress drops are not well constrained (> 10
4
MPa). (b) Stress
drops vs. aspect ratio 𝐿 𝑐 /𝑊 𝑐 . (c) Stress drops vs. focal depths. (d)The ratio of stress drops of
elliptical/circular cracks with the same rupture length. This ratio increases significantly for
earthquakes with aspect ratios, i.e. 𝐿 𝑐 /𝑊 𝑐 , > 3.

Appendix 5.A: Coordinate rotation
With the FM of the target earthquake, the number of unknown parameters can be reduced
from 10 to 6 by constraining the rupture onto the 2D fault plane. In this appendix we derive
the rotation matrices. To start with, the geographical Cartesian coordinate system is defined
as 𝑥 positive North, 𝑦 positive East and z positive Down. The rotated coordinate is defined
as 𝑥 ’ along strike, 𝑦 ’ fault normal, 𝑧 ’ down-dip. As a result, the direction vectors of 𝑥 ’𝑦 ’𝑧 ’ in
𝑥 yz coordinate system is

𝑒 ̂
𝑥 ′
= [cos ( 𝑠𝑡𝑟𝑖𝑘𝑒 ) sin ( 𝑠𝑡𝑟𝑖𝑘𝑒 ) 0]
T

𝑒 ̂
𝑧 ′
= [−sin ( 𝑠𝑡𝑟𝑖𝑘𝑒 )∙ cos ( 𝑑𝑖𝑝 ) cos ( 𝑠𝑡𝑟𝑖𝑘𝑒 )∙ cos ( 𝑑𝑖𝑝 ) sin ( 𝑑𝑖𝑝 ) ]
T

(5.A1)
Doctoral Dissertation Haoran Meng

129
𝑒 ̂
𝑦 ′
= 𝑒 ̂
𝑧 ′
× 𝑒 ̂
𝑥 ′
.
We use a two-step rotation to rotate a vector in 𝑥 yz coordinate to 𝑥 ′y′z’ coordinate. First, we
keep z still, rotate x and y to align x to the strike direction. The rotation matrix is
𝑅 1
= [
cos ( 𝑠𝑡𝑟𝑖𝑘𝑒 ) sin ( 𝑠𝑡𝑟𝑖𝑘𝑒 ) 0
−sin ( 𝑠𝑡𝑟𝑖𝑘𝑒 ) cos ( 𝑠𝑡𝑟𝑖𝑘𝑒 ) 0
0 0 1
]. (5.A2)
We then keep the new 𝑥 ( 𝑥 ’) still, rotate 𝑧 to the down-dip direction using
𝑅 2
= [
1 0 0
0 sin ( 𝑑𝑖𝑝 ) −cos ( 𝑑𝑖𝑝 )
0 cos ( 𝑑𝑖𝑝 ) sin ( 𝑑𝑖𝑝 )
]. (5.A3)
To rotate a vector 𝑉 in 𝑥𝑦𝑧 coordinate to 𝑉 ′ in 𝑥 ’𝑦 ’𝑧 ’ coordinate,
𝑉 ′=𝑅 2
∙ (𝑅 1
∙ 𝑉 ). (5.A4)
One can verify the rotation by checking

𝑉 𝑥 ′
= 𝑉 T
∙ 𝑒 ̂
𝑥 ′

𝑉 𝑦 ′
= 𝑉 T
∙ 𝑒 ̂
𝑦 ′

𝑉 𝑧 ′
= 𝑉 T
∙ 𝑒 ̂
𝑧 ′
.
(5.A5)
To rotate the vector 𝑉 ′ in 𝑥 ’𝑦 ’𝑧 ’ coordinate back to 𝑉 in 𝑥𝑦𝑧 coordinate,
𝑉 = 𝑅 1
T
∙ (𝑅 2
T
∙ 𝑉 ′). (5.A6)
Doctoral Dissertation Haoran Meng

130

Figure 5.S1. (a) The maximal absolute value of the Cross-correlation Coefficients (CCs) of
the bandpass filtered (0.5 to 2 Hz) S phase waveforms of the target events and its eGf
candidates. The texts show the station names and SCSN event ids. (b) The mean CC and the
completeness of data. (c). The black curves show the CC functions of the waveforms of the
target event and best performed eGf at two example stations IDO and RVR. The red curves
denote the ASTFs using PLD. We note that the major pulses of the CC functions are very
similar to the ASTFs.
Doctoral Dissertation Haoran Meng

131

Figure 5.S2. (a) The waveforms and fits of the target 2014, Mw 4.58 earthquake. (b) The
waveforms and fits of the target 2013, Mw 4.70 earthquake. P phases are in vertical and S
phases are in the transverse directions.

Doctoral Dissertation Haoran Meng

132

Figure 5.S3. (a) Characteristic durations, (b) ASTFs and (c) normalized misfits for the best
performed and weighted stacked eGfs for the 2016, Mw 5.19 event. P phase in vertical and S
phase in transverse directions. Texts display the station names and corresponding back
azimuths for all stations.

Doctoral Dissertation Haoran Meng

133
Table 5.1. Target events analyzed in this study.
No. Date Event ID Latitude Longitude
Depth
(km)
M Strike Dip Rake
1 2016-06-10 37374687 33.44031 -116.4352 12.433 5.19 303 66 179
2 2014-03-29 15481673 33.90461 -117.95333 6.005 5.1 130 55 165
3 2010-06-13 10701405 33.39017 -116.39906 7.856 4.9 331 83 -178
4 2013-01-11 15296281 33.50576 -116.4586 9.954 4.7 315 76 -167
5 2014-07-05 15520985 34.27669 -117.02524 7.972 4.58 306 72 -179
6 2018-05-08 38167848 34.021 -116.78521 13.535 4.49 277 58 132
7 2012-08-08 15189073 33.91105 -117.7897 9.532 4.46 130 86 -174
8 2012-08-08 15189281 33.91018 -117.78671 9.318 4.45 121 71 180
9 2018-08-15 38245496 33.49083 -116.78945 4.114 4.43 344 90 171
10 2014-01-15 11413954 34.14743 -117.44549 3.864 4.43 303 63 178
11 2018-08-29 38038071 34.12779 -117.78223 4.717 4.41 315 74 -162
12 2015-12-30 37507576 34.19584 -117.41829 8.49 4.4 308 81 178
13 2016-01-06 37510616 33.9715 -116.88264 15.995 4.39 307 54 140
14 2010-01-16 10530013 33.9401 -117.01966 15.099 4.28 327 84 178
15 2010-01-12 14571828 33.97229 -116.87233 9.069 4.27 345 66 171
16 2008-05-01 10321561 33.446 -116.43665 8.615 4.19 316 65 168
17 2014-03-29 15483001 33.9563 -117.89722 6.615 4.14 114 63 -168
18 2011-09-14 11006189 33.96212 -117.06751 16.492 4.14 333 57 -172
19 2008-10-02 14396336 34.09154 -116.9718 14.984 4.14 98 44 96
20 2012-08-29 15207433 33.9063 -117.79178 9.016 4.13 128 71 -178
21 2008-11-17 14403732 33.50959 -116.84798 11.557 4.11 314 87 -164
22 2010-02-13 10541957 34.01192 -117.18447 8.832 4.1 174 59 -143
23 2008-05-09 14367532 33.44844 -116.4403 7.848 4.06 314 78 175
24 2015-09-16 37243591 34.13885 -116.85416 11.858 4 304 57 -157
25 2012-06-14 15164105 33.91156 -117.78737 9.694 3.99 130 72 179
26 2008-06-23 14376612 34.05278 -117.24309 15.231 3.99 339 63 -177
27 2009-04-24 10399889 33.9045 -117.79791 0.004 3.98 114 76 172
28 2018-01-25 38092312 33.7495 -117.50117 10.711 3.97 273 43 106
29 2018-04-26 37924223 33.38645 -116.29694 11.293 3.94 288 80 169
30 2016-12-28 37772688 34.15459 -116.70373 12.7 3.92 137 89 -175
31 2012-10-28 15237073 33.70513 -116.80654 16.852 3.87 297 38 138
32 2018-04-23 37920791 33.92635 -116.31881 9.57 3.87 333 75 176
33 2016-02-16 37524376 34.29963 -116.86221 5.522 3.87 303 75 -178
Doctoral Dissertation Haoran Meng

134
34 2009-04-16 14444832 33.25538 -116.4199 3.221 3.86 318 72 171
35 2014-05-19 15503377 34.25475 -116.82586 8.448 3.84 317 84 179
36 2012-04-28 15141521 34.24089 -117.43179 12.971 3.82 283 43 102
37 2015-07-25 37213455 34.09732 -117.44811 4.849 3.81 308 77 175
38 2013-09-20 11366994 33.35229 -116.39017 12.175 3.8 334 76 179
39 2008-07-29 14384052 33.94852 -117.81039 13.9 3.8 125 56 163
40 2011-05-27 14992276 34.22025 -117.05263 0.728 3.8 286 81 -166
41 2014-08-15 15538561 34.29708 -116.44877 10.351 3.8 345 83 -177
42 2016-06-10 37377079 33.46393 -116.41829 10.404 3.46 172 68 -154

Table 5.2. Second moment inversion results of the 2014 Mw 4.58 Event (ID: 15520985).
eGf type
𝐿 𝑐
(km)
𝑊 𝑐
(km)
𝜏 𝑐 (s)
𝑣 0

(km/s)
Directivity
ratio
Residual
per degree
of freedom
(s)
Number of
measurements
Single eGf 1.12 0.48 0.37 2.70 0.87 0.0278 26
Stacked eGf 1.12 0.50 0.37 2.73 0.89 0.0223 36

Table 5.3. Results for the analyzed event using weighted stacked eGf.
Event ID M
L c
(km)
W c
(km/s)
𝜏 𝑐 (s)
Directivity
ratio
Stress
drop
(MPa)
Min
stress
drop
(MPa)
Max
stress
drop
(MPa)
37374687 5.19 1.77 0.82 0.294 0.85 20.90 3.14 177.87
10701405 4.90 1.16 0.36 0.208 0.50 58.81 10.14 1811.32
15296281 4.70 1.64 0.58 0.244 0.32 8.12 4.00 11.11
15520985 4.58 1.12 0.52 0.366 0.89 9.98 1.21 69.99
38038071 4.41 0.65 0.49 0.171 0.59 4.54 1.46 1362.45
37507576 4.40 0.73 0.31 0.172 0.63 21.75 2.46 82.01
10530013 4.28 0.29 0.19 0.074 0.55 44.83 10.59 -
14571828 4.27 0.48 0.13 0.094 0.43 57.02 2.41 870.13
10541957 4.10 0.61 0.14 0.073 0.29 24.62 2.12 187.82
14367532 4.06 0.55 0.12 0.058 0.36 41.65 9.06 -
14444832 3.86 0.28 0.06 0.065 0.31 138.51 7.26 410.03
37213455 3.81 0.24 0.11 0.052 0.39 66.73 7.60 709.86
11366994 3.80 0.26 0.15 0.053 0.28 21.29 5.14 -
37377079 3.46 0.15 0.06 0.049 0.38 82.53 9.83 -
Doctoral Dissertation Haoran Meng

135
- max stress drop estimation exceeds 10
4
MPa.

Doctoral Dissertation Haoran Meng

136
References
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.
Allen, R.V., 1978, Automatic earthquake recognition and timing from single traces: Bulletin
of the Seismological Society of America, v. 68, p. 1521-1532.
Ampuero, J.-P., and Ben-Zion, Y., 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.x.
Andrews, D.J., and Ben-Zion, Y., 1997, Wrinkle-like slip pulse on a fault between different
materials: Journal of Geophysical Research: Solid Earth, v. 102, p. 553–571, doi:
10.1029/96jb02856.
Backus, G., and Mulcahy, M., 1976a, Moment Tensors and other Phenomenological
Descriptions of Seismic Sources—I. Continuous Displacements: Geophysical Journal
International, v. 46, p. 341–361, doi: 10.1111/j.1365-246x.1976.tb04162.x.
Backus, G., and Mulcahy, M., 1976b, Moment tensors and other phenomenological
descriptions of seismic sources--II. Discontinuous displacements: Geophysical Journal
International, v. 47, p. 301–329, doi: 10.1111/j.1365-246x.1976.tb01275.x.
Backus, G.E., 1977a, Interpreting the seismic glut moments of total degree two or less:
Geophysical Journal International, v. 51, p. 1–25, doi: 10.1111/j.1365-
246x.1977.tb04187.x.
Backus, G.E., 1977b, Seismic sources with observable glut moments of spatial degree two:
Geophysical Journal International, v. 51, p. 27–45, doi: 10.1111/j.1365-
246x.1977.tb04188.x.
Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., Moschetti, M.P., Shapiro,
N.M., and Yang, Y., 2007, Processing seismic ambient noise data to obtain reliable
broad-band surface wave dispersion measurements: Geophysical Journal International,
v. 169, p. 1239–1260, doi: 10.1111/j.1365-246x.2007.03374.x.
Doctoral Dissertation Haoran Meng

137
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, doi:
10.1046/j.1365-246x.2002.01637.x.
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, doi: 10.1029/2008rg000260.
Ben-Zion, Y., and Andrews, D. J., 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., Vernon, F.L., Ozakin, Y., 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/ggv142.
Bertero, M., Bindi, D., Boccacci, P., Cattaneo, M., Eva, C., and Lanza, V., 1997, Application of
the projected Landweber method to the estimation of the source time function in
seismology: Inverse Problems, v. 13, p. 465–486, doi: 10.1088/0266-5611/13/2/017.
Boyd, S.P., and Vandenberghe, L., 2018, Convex optimization: Cambridge, Cambridge
University Press.
Brenguier, F., Campillo, M., Hadziioannou, C., Shapiro, N.M., Nadeau, R.M., and Larose, E.,
2008, Postseismic Relaxation Along the San Andreas Fault at Parkfield from
Continuous Seismological Observations: Science, v. 321, p. 1478–1481, doi:
10.1126/science.1160943.
Brietzke, G.B., Cochard, A., and Igel, H., 2009, Importance of bimaterial interfaces for
earthquake dynamics and strong ground motion: Geophysical Journal International, v.
178, p. 921–938, doi: 10.1111/j.1365-246x.2009.04209.x.
Brune, J. N., 1970, Tectonic stress and the spectra of seismic shear waves from earthquakes,
Journal of Geophysical Research, v. 75, p. 4997– 5009, doi:10.1029/JB075i026p04997.
Calderoni, G., Rovelli, A., and Singh, S.K., 2012, Stress drop and source scaling of the 2009
April L’Aquila earthquakes: Geophysical Journal International, v. 192, p. 260–274, doi:
10.1093/gji/ggs011.
Doctoral Dissertation Haoran Meng

138
Calderoni, G., Rovelli, A., Ben-Zion, Y., and Giovambattista, R.D., 2015, Along-strike rupture
directivity of earthquakes of the 2009 LAquila, central Italy, seismic sequence:
Geophysical Journal International, v. 203, p. 399–415, doi: 10.1093/gji/ggv275.
Campillo, M., Roux, P. and Shapiro, N.M., 2011, Using seismic noise to image and to monitor
the Solid Earth, in Encyclopedia of Solid Earth Geophysics, ed. Gupta, Harsh K., pp.
1230–1235, Springer.
Chen, X., Nakata, N., Pennington, C., Haffener, J., Chang, J.C., He, X., Zhan, Z., Ni, S., and Walter,
J.I., 2017, The Pawnee earthquake as a result of the interplay among injection, faults
and foreshocks: Scientific Reports, v. 7, doi: 10.1038/s41598-017-04992-z.
Corciulo, M., Roux, P., Campillo, M., and Dubucq, D., 2012, Instantaneous phase variation for
seismic velocity monitoring from ambient noise at the exploration scale: Geophysics, v.
77, doi: 10.1190/geo2011-0363.1.
Cros, E., Roux, P., Vandemeulebrouck, J., and Kedar, S., 2011, Locating hydrothermal
acoustic sources at Old Faithful Geyser using Matched Field Processing: Geophysical
Journal International, v. 187, p. 385–393, doi: 10.1111/j.1365-246x.2011.05147.x.
Eibl, E.P., Lokmer, I., Bean, C.J., Akerlie, E., and Vogfjörd, K.S., 2015, Helicopter vs. volcanic
tremor: Characteristic features of seismic harmonic tremor on volcanoes: Journal of
Volcanology and Geothermal Research, v. 304, p. 108–117, doi:
10.1016/j.jvolgeores.2015.08.002.
Eshelby, J.D, 1957, The determination of the elastic field of an ellipsoidal inclusion, and
related problems: Solid Mechanics and its Applications Collected Works of J. D. Eshelby,
p. 209–229, doi: 10.1007/1-4020-4499-2_18.
Fan, W., and Mcguire, J.J., 2018, Investigating microearthquake finite source attributes with
IRIS Community Wavefield Demonstration Experiment in Oklahoma: Geophysical
Journal International, v. 214, p. 1072–1087, doi: 10.1093/gji/ggy203.
Fink, M., 2006, Time-reversal acoustics in complex environments: Geophysics, v. 71, doi:
10.1190/1.2215356.
Fuchs, F., and Bokelmann, G., 2017, Equidistant Spectral Lines in Train Vibrations:
Seismological Research Letters, v. 89, p. 56–66, doi: 10.1785/0220170092.
Geiger, L., 1912, Probability method for the determination of earthquake epicenters from
the arrival time only, Bulletin of Saint Louis University, 8(1), 56-71.
Doctoral Dissertation Haoran Meng

139
Gerstoft, P., and Tanimoto, T., 2007, A year of microseisms in southern California:
Geophysical Research Letters, v. 34, doi: 10.1029/2007gl031091.
Gibbons, S.J., and Ringdal, F., 2006, The detection of low magnitude seismic events using
array-based waveform correlation: Geophysical Journal International, v. 165, p. 149-
166, doi: 10.31223/osf.io/5z9ew.
Gradon, C., Moreau, L., Roux, P., and Ben-Zion, Y., 2019, Analysis of surface and seismic
sources in dense array data with match field processing and Markov chain Monte Carlo
sampling: Geophysical Journal International, v. 218, p. 1044–1056, doi:
10.1093/gji/ggz224.
Grant, M.C., and Boyd, S.P., 2008, Graph Implementations for Nonsmooth Convex Programs:
Lecture Notes in Control and Information Sciences Recent Advances in Learning and
Control, p. 95–110, doi: 10.1007/978-1-84800-155-8_7.
Grant, M.C., and Boyd, S.P., 2014, CVX: Matlab Software for Disciplined Convex
Programming CVX Research, Inc., http://cvxr.com/cvx/.
Hanks, T. C. and Kanamori, H., 1979. A moment magnitude scale: Geophysical Journal
International, v. 84, p. 2348-2350, doi: 10.1029/JB084iB05p02348
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, p. 2239–2244, doi: 10.1785/0120120010.
Hough, S. E., and Dreger, D. S., 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.
Hillers, G., Ben-Zion, Y., Landès, M., and Campillo, M., 2013, Interaction of microseisms with
crustal heterogeneity: A case study from the San Jacinto fault zone area: Geochemistry,
Geophysics, Geosystems, v. 14, p. 2182–2197, doi: 10.1002/ggge.20140.
Hillers, G., Roux, P., Campillo, M., and Ben ‐Zion, Y., 2016, Focal spot imaging based on zero
lag cross ‐correlation amplitude fields: Application to dense array data at the San
Jacinto fault zone: Journal of Geophysical Research: Solid Earth, v. 121, p. 8048–8067,
doi: 10.1002/2016jb013014.
Doctoral Dissertation Haoran Meng

140
Hutchings, L., and Wu, F., 1990, Empirical Greens Functions from small earthquakes: A
waveform study of locally recorded aftershocks of the 1971 San Fernando Earthquake:
Journal of Geophysical Research, v. 95, p. 1187, doi: 10.1029/jb095ib02p01187.
Inbal, A., Ampuero, J.P., and Clayton, R.W., 2016, Localized seismic deformation in the upper
mantle revealed by dense seismic arrays: Science, v. 354, p. 88–92, doi:
10.1126/science.aaf1370.
Inbal, A., Clayton, R.W., and Ampuero, J.-P., 2015, Imaging widespread seismicity at
midlower crustal depths beneath Long Beach, CA, with a dense seismic array: Evidence
for a depth-dependent earthquake size distribution: Geophysical Research Letters, v.
42, p. 6314–6323, doi: 10.1002/2015gl064942.
Inbal, A., Cristea ‐Platon, T., Ampuero, J.P., Hillers, G., Agnew, D., and Hough, S.E., 2018,
Sources of Long ‐Range Anthropogenic Noise in Southern California and Implications
for Tectonic Tremor Detection: Bulletin of the Seismological Society of America, doi:
10.1785/0120180130.
Johnson, C., Meng, H., Vernon, F., Nakata, N. and Ben-Zion, Y., 2019, Characteristics of
ground motion generated by interaction of wind gusts with trees, structures and other
obstacles above the surface: Journal of Geophysical Research: Solid Earth, in revision.
Kurzon, I., Vernon, F.L., 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, p. 3045–3081, doi:
10.1007/s00024-014-0855-2.
Kwiatek, G., and Ben-Zion, Y., 2016, Theoretical limits on detection and analysis of small
earthquakes: Journal of Geophysical Research: Solid Earth, v. 121, p. 5898–5916, doi:
10.1002/2016jb012908.
Lanza, V., Spallarossa, D., Cattaneo, M., Bindi, D., and Augliera, P., 1999, Source parameters
of small events using constrained deconvolution with empirical Greens functions:
Geophysical Journal International, v. 137, p. 651–662, doi: 10.1046/j.1365-
246x.1999.00809.x.
Larmat, C., Tromp, J., Liu, Q., and Montagner, J.-P., 2008, Time reversal location of glacial
earthquakes: Journal of Geophysical Research, v. 113, doi: 10.1029/2008jb005607.
Doctoral Dissertation Haoran Meng

141
Lee, E.J., Chen, P., Jordan, T.H., Maechling, P.B., Denolle, M.A. and Beroza, G.C., 2014. Full‐3‐D
tomography for crustal structure in Southern California based on the scattering‐
integral and the adjoint‐wavefield methods: Journal of Geophysical Research, v. 119, p.
6421-6451, doi: 10.1002/2014JB011346.
Li, Z., and Peng, Z., 2016, An Automatic Phase Picker for Local Earthquakes with
Predetermined Locations: Combining a Signal ‐to ‐Noise Ratio Detector with 1D
Velocity Model Inversion: Seismological Research Letters, v. 87, p. 1397 –1405, doi:
10.1785/0220160027.
Li, Z., Peng, Z., Hollis, D., Zhu, L., and Mcclellan, J., 2018, High-resolution seismic event
detection using local similarity for Large-N arrays: Scientific Reports, v. 8, doi:
10.1038/s41598-018-19728-w.
Liu, X., and Ben-Zion, Y., 2016, Estimating correlations of neighbouring frequencies in
ambient seismic noise: Geophysical Journal International, v. 206, p. 1065–1075, doi:
10.1093/gji/ggw196.
Madariaga, R., 1976, Dynamics of an expanding circular fault: Bulletin of the Seismological
Society of America, v. 66, p. 639–666.
Margerin, L., Campillo, M., Tiggelen, B.A.V., and Hennino, R., 2009, Energy partition of
seismic coda waves in layered media: theory and application to Pinyon Flats
Observatory: Geophysical Journal International, v. 177, p. 571–585, doi:
10.1111/j.1365-246x.2008.04068.x.
Mcguire, J.J., 2004, Estimating Finite Source Properties of Small Earthquake Ruptures:
Bulletin of the Seismological Society of America, v. 94, p. 377–393, doi:
10.1785/0120030091.
Mcguire, J.J., Zhao, L., and Jordan, T.H., 2001, Teleseismic inversion for the second-degree
moments of earthquake space-time distributions: Geophysical Journal International, v.
145, p. 661–678, doi: 10.1046/j.1365-246x.2001.01414.x.
Mcguire, J.J., 2017, A MATLAB Toolbox for Estimating the Second Moments of Earthquake
Ruptures: Seismological Research Letters, v. 88, p. 371–378, doi:
10.1785/0220160170.
Doctoral Dissertation Haoran Meng

142
Mcguire, J.J., and Kaneko, Y., 2018, Directly estimating earthquake rupture area using
second moments to reduce the uncertainty in stress drop: Geophysical Journal
International, v. 214, p. 2224–2235, doi: 10.1093/gji/ggy201.
Meng, X., Peng, Z., and Hardebeck, J.L., 2013, Seismicity around Parkfield correlates with
static shear stress changes following the 2003Mw6.5 San Simeon earthquake: Journal
of Geophysical Research: Solid Earth, v. 118, p. 3576–3591, doi: 10.1002/jgrb.50271.
Meng, H., and Ben-Zion, Y., 2018a, Detection of small earthquakes with dense array data:
example from the San Jacinto fault zone, southern California: Geophysical Journal
International, v. 212, p. 442–457, doi: 10.1093/gji/ggx404.
Meng, H., and Ben-Zion, Y., 2018b, Characteristics of Airplanes and Helicopters Recorded by
a Dense Seismic Array Near Anza California: Journal of Geophysical Research: Solid
Earth, v. 123, p. 4783–4797, doi: 10.1029/2017jb015240.
Meng, H., Ben-Zion, Y. and Johnson, C., 2019. Detection of random noise and anatomy of
continuous seismic waveforms in dense array data near Anza California: Geophysical
Journal International, in revision.
Meng, H., McGuire, J. and Ben-Zion, Y., 2019. Semi-automated estimates of directivity and
related source properties of small to moderate Southern California earthquakes with
second moments, in preparation.
Meng, H., and Ben-Zion, Y., 2019, Evaluation of diffuse waveforms in dense array data near
Anza California, in preparation.
Mordret, A., Roux, P., Boué, P., and Ben-Zion, Y., 2018, Shallow three-dimensional structure
of the San Jacinto fault zone revealed from ambient noise imaging with a dense seismic
array: Geophysical Journal International, v. 216, p. 896–905, doi: 10.1093/gji/ggy464.
Mueller, C.S., 1985, Source pulse enhancement by deconvolution of an empirical Greens
function: Geophysical Research Letters, v. 12, p. 33–36, doi:
10.1029/gl012i001p00033.
Nakata, N., Chang, J.P., and Lawrence, J.F, and Boué, P., 2015, Body wave extraction and
tomography at Long Beach, California, with ambient-noise interferometry, Journal of
Geophysical Research: Solid Earth, v. 120, p. 1159-1173., doi:10.1002/2015JB011870.
Doctoral Dissertation Haoran Meng

143
Nakata, N., and Beroza, G.C., 2015, Stochastic characterization of mesoscale seismic velocity
heterogeneity in Long Beach, California: Geophysical Journal International, v. 203, p.
2049–2054, doi: 10.1093/gji/ggv421.
Neuffer, T., and Kremers, S., 2017, How wind turbines affect the performance of seismic
monitoring stations and networks: Geophysical Journal International, v. 211, p. 1319–
1327, doi: 10.1093/gji/ggx370.
Qin, L., Ben-Zion, Y., Qiu, H., Share, P.-E., Ross, Z.E., and Vernon, F.L., 2017, Internal structure
of the San Jacinto fault zone in the trifurcation area southeast of Anza, California, from
data of dense seismic arrays: Geophysical Journal International, v. 213, p. 98–114, doi:
10.1093/gji/ggx540.
Ross, Z.E., and Ben-Zion, Y., 2016, Toward reliable automated estimates of earthquake
source properties from body wave spectra: Journal of Geophysical Research: Solid
Earth, v. 121, p. 4390–4407, doi: 10.1002/2016jb013003.
Ross, Z.E., Hauksson, E., and Ben-Zion, Y., 2017, Abundant off-fault seismicity and
orthogonal structures in the San Jacinto fault zone: Science Advances, v. 3, doi:
10.1126/sciadv.1601946.
Rost, S., 2002, Array seismology: Methods and applications: Reviews of Geophysics, v. 40,
doi: 10.1029/2000rg000100.
Roux, P., Moreau, L., Lecointre, A., Hillers, G., Campillo, M., Ben-Zion, Y., Zigone, D., and
Vernon, F., 2016, A methodological approach towards high-resolution surface wave
imaging of the San Jacinto Fault Zone using ambient-noise recordings at a spatially
dense array: Geophysical Journal International, v. 206, p. 980–992, doi:
10.1093/gji/ggw193.
Salvermoser, J., Hadziioannou, C., and Stähler, S.C., 2015, Structural monitoring of a
highway bridge using passive noise recordings from street traffic: The Journal of the
Acoustical Society of America, v. 138, p. 3864–3872, doi: 10.1121/1.4937765.
Sánchez-Sesma, F.J., Pérez-Ruiz, J.A., Luzón, F., Campillo, M., and Rodríguez-Castellanos, A.,
2008, Diffuse fields in dynamic elasticity: Wave Motion, v. 45, p. 641–654, doi:
10.1016/j.wavemoti.2007.07.005.
Doctoral Dissertation Haoran Meng

144
Sens-Schönfelder, C., Snieder, R., and Stähler, S.C., 2015, The lack of equipartitioning in
global body wave coda: Geophysical Research Letters, v. 42, p. 7483–7489, doi:
10.1002/2015gl065108.
Shapiro, N.M., Campillo, M., Stehly, L. and Ritzwoller, M.H., 2005, High-Resolution Surface-
Wave Tomography from Ambient Seismic Noise: Science, v. 307, p. 1615–1618, doi:
10.1126/science.1108339.
Shapiro, N.M., Campillo, M., Margerin, L., Singh, S.K., Kostoglodov, V. and Pacheco, J., 2000,
The Energy Partitioning and the Diffusive Character of the Seismic Coda: Bulletin of the
Seismological Society of America, v. 90, p. 655–665, doi: 10.1785/0119990021.
Share, P.-E., Ben-Zion, Y., Ross, Z.E., Qiu, H., and Vernon, F.L., 2017, Internal structure of the
San Jacinto fault zone at Blackburn Saddle from seismic data of a linear array:
Geophysical Journal International, v. 210, p. 819–832, doi: 10.1093/gji/ggx191.
Share, P.-E., and Ben-Zion, Y., 2016, Bimaterial interfaces in the south San Andreas Fault
with opposite velocity contrasts NW and SE from San Gorgonio Pass: Geophysical
Research Letters, v. 43, doi: 10.1002/2016gl070774.
Shelly, D.R., Beroza, G.C., and Ide, S., 2007, Non-volcanic tremor and low-frequency
earthquake swarms: Nature, v. 446, p. 305–307, doi: 10.1038/nature05666.
Shelly, D.R., Ellsworth, W.L., and Hill, D.P., 2016, Fluid-faulting evolution in high definition:
Connecting fault structure and frequency-magnitude variations during the 2014 Long
Valley Caldera, California, earthquake swarm: Journal of Geophysical Research: Solid
Earth, v. 121, p. 1776–1795, doi: 10.1002/2015jb012719.
Slawinski, R.A. and Slawinski, M.A., 1999, On raytracing in constant velocity-gradient
media: calculus approach: Canadian Journal of Exploration Geophysics, v. 35, p. 24-27.
Sleep, N.H., 2011, Seismically damaged regolith as self-organized fragile geological feature:
Geochemistry, Geophysics, Geosystems, v. 12, doi: 10.1029/2011gc003837.
Tanimoto, T., and Valovcin, A., 2015, Stochastic excitation of seismic waves by a hurricane:
Journal of Geophysical Research: Solid Earth, v. 120, p. 7713–7728, doi:
10.1002/2015jb012177.
Turek, G., and Kuperman, W.A., 1997, Applications of matched-field processing to structural
vibration problems: The Journal of the Acoustical Society of America, v. 101, p. 1430–
1440, doi: 10.1121/1.418168.
Doctoral Dissertation Haoran Meng

145
Vandenberghe, L., and Boyd, S., 1996, Semidefinite programming: SIAM review, v. 38(1), p.
49-95, doi: 10.1137/1038003
Vernon, F., Ben-Zion, Y., and Hollis, D., 2014, Sage Brush Flats Nodal Experiment.
International Federation of Digital Seismograph Networks: Other/Seismic Network.
doi:10.7914/SN/ZG_2014.
Weaver, R. L., 1982, On diffuse waves in solid media: The Journal of the Acoustical Society
of America, v. 71(6), p. 1608-1609, doi:10.1121/1.387816
Weaver, R., and Lobkis, O., 2004, Diffuse fields in open systems and the emergence of the
Green’s function: The Journal of the Acoustical Society of America, v. 117, p. 2394–
2394, doi: 10.1121/1.4809391.
Weertman, J., 1980, Unstable slippage across a fault that separates elastic media of
different elastic constants: Journal of Geophysical Research: Solid Earth, v. 85, p. 1455–
1461, doi: 10.1029/jb085ib03p01455.
Yang, W., and Ben-Zion, Y., 2010, An Algorithm for Detecting Clipped Waveforms and
Suggested Correction Procedures: Seismological Research Letters, v. 81, p. 53–62, doi:
10.1785/gssrl.81.1.53.
Yang, W., Hauksson, E., and Shearer, P.M., 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, p. 1179–1194, doi:
10.1785/0120110311.
Yoon, C.E., O’Reilly, O., Bergen, K.J., and Beroza, G.C., 2015, Earthquake detection through
computationally efficient similarity search: Science Advances, v. 1, doi:
10.1126/sciadv.1501057.
Zaliapin, I., and Ben-Zion, Y., 2013, Earthquake clusters in southern California I:
Identification and stability: Journal of Geophysical Research: Solid Earth, v. 118, p.
2847–2864, doi: 10.1002/jgrb.50179.
Zaliapin, I., and Ben-Zion, Y., 2016, A global classification and characterization of
earthquake clusters: Geophysical Journal International, v. 207, p. 608–634, doi:
10.1093/gji/ggw300.
Doctoral Dissertation Haoran Meng

146
Zhao, L. S., and Helmberger, D.V., 1993, Source retrieval from broadband regional
seismograms: Hindu Kush region: Physics of the Earth and Planetary Interiors, v. 78, p.
69–95, doi: 10.1016/0031-9201(93)90085-n.
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.
Zigone, D., Rivet, D., Radiguet, M., Campillo, M., Voisin, C., Cotte, N., Walpersdorf, A., Shapiro,
N.M., Cougoulat, G., Roux, P., Kostoglodov, V., Husker, A., and Payero, J.S., 2012,
Triggering of tremors and slow slip event in Guerrero, Mexico, by the 2010 Mw 8.8
Maule, Chile, earthquake: Journal of Geophysical Research: Solid Earth, v. 117, doi:
10.1029/2012jb009160.
Zigone, D., Ben-Zion, Y., Campillo, M., and Roux, P., 2014, Seismic Tomography of the
Southern California Plate Boundary Region from Noise-Based Rayleigh and Love
Waves: Pure and Applied Geophysics, v. 172, p. 1007–1032, doi: 10.1007/s00024-014-
0872-1.

Abstract (if available)

Abstract I develop a collection of methodologies detecting small earthquakes, modeling air-traffic events, separating continuous seismic waveforms into random noise, not-random-noise, evaluating diffuse waveform and estimating of directivity, rupture area, duration, and centroid velocity of earthquakes using second seismic moments and stacked empirical Green’s functions. These techniques are suitable for systematic processing of large seismic waveform data sets to extract detailed properties of seismic waveforms and earthquake sources. ❧ A technique is developed to detect small earthquakes not included in standard catalogs using data from a dense seismic array. The technique is illustrated with continuous waveforms recorded in a test day by 1108 vertical geophones in a tight array on the San Jacinto fault zone. The number of validated detected events in the test day is > 5 times larger than that in the standard catalog. Using these events as templates can lead to additional detections of many more earthquakes. ❧ Frequent air-traffic events are observed in continuous seismic waveforms recorded for about 30 days in the dense array on the San Jacinto fault zone. Time-frequency analysis reveals clear Doppler effects that can be modeled with basic equations and fitted well with parameters corresponding to airplanes and helicopters. The total event duration covers > 7% of the day, which is likely to exceed considerably the total time covered by earthquakes. ❧ A methodology is developed to separate continuous seismic waveforms into random noise (RN), not-random-noise (NRN) produced by earthquakes, wind, traffic and other sources, and the undetermined mixture of signals. The analysis is applied to continuous records of the dense seismic array on the San Jacinto fault zone. For the days examined, the relative fractions of RN, NRN and mixed signals in local day (night) times are about 26% (42%), 40% (33%) and 34% (25%), respectively. ❧ A technique is developed to evaluate whether the recorded waveforms are diffused or not by examining the expectation of waveform spectra and the outer products with themselves and their conjugates. The method is applied to the detected RN and earthquake coda. The results show the detected RN are well diffuse while the recorded waveforms in the coda window of the 2010 Mw 7.2 El Mayor Cucapah are not necessarily diffuse because of the triggered events or instrumental effect. ❧ A method is developed for semi-automated estimation of directivity, rupture area, duration, and centroid velocity of earthquakes with second seismic moments. The method is applied to 14 earthquakes with magnitude in the range 3.5-5.2 in Southern California. Most analyzed target events in the Trifurcation area of the San Jacinto Fault system have directivities towards the northwest, while events around Cajon Pass and San Gabriel Mountain propagate towards the southeast. The results are generally consistent with predictions for dynamic rupture on a bimaterial interface and the imaged velocity contrasts in the study area. The second moment inversions are also used to explore constraints on upper and lower bounds of rupture areas in our dataset. We note that the stress drop can be significantly underestimated by assuming a circular crack. The second moment method with stacked eGfs appears capable of semi-automated, routine application to moderate earthquakes in southern California.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

PDF

Stress-strain characterization of complex seismic sources

PDF

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

PDF

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

PDF

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

PDF

Dynamic rupture processes and seismic radiation in models of earthquake faults separating similar and dissimilar solids

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

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

PDF

Multi-scale imaging of major fault zones in Southern California

PDF

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

PDF

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

PDF

Integration and validation of deterministic earthquake simulations in probabilistic seismic hazard analysis

PDF

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

PDF

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

PDF

Wave method for structural system identification and health monitoring of buildings based on layered shear beam model

PDF

Detection and modeling of slow slip events as creep instabilities beneath major fault zones

PDF

Paleoseismologic and slip rate studies of three major faults in southern California: understanding the complex behavior of plate boundary fault systems over millenial timescales

Asset Metadata

Creator Meng, Haoran (author)

Core Title Detailed properties of seismic waveforms and earthquake sources in Southern California

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Geological Sciences

Publication Date 07/24/2019

Defense Date 04/10/2019

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

Tag anthropogenic sources,dense seismic array,earthquake dynamics,earthquake source observations,empirical Green's function,OAI-PMH Harvest,rupture directivity,second seismic moment,seismic noise,signal classification,site effects,stress drop,wave propagation,wave scattering and diffraction

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Ben-Zion, Yehuda (committee chair), Alexander, Kenneth (committee member), Frank, William (committee member)

Creator Email haoranme@usc.edu,seismo.menghaoran@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-187396

Unique identifier UC11663021

Identifier etd-MengHaoran-7586.pdf (filename),usctheses-c89-187396 (legacy record id)

Legacy Identifier etd-MengHaoran-7586.pdf

Dmrecord 187396

Document Type Dissertation

Format application/pdf (imt)

Rights Meng, Haoran

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