University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

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

Hongrui Qiu

A Dissertation

Presented to the Faculty of the University of Southern California

Graduate School

In partial fulfillment of requirements for the degree of

Doctor of Philosophy

in

Geological Sciences

August 2019

1
Acknowledgements
I would like to first thank my parents for their love and support. From a young age, they not
only taught me to be hardworking and persistent in the study but also how to be a good person. I
was so blessed to have my beloved girlfriend, Lei, and she is always a source of happiness in my
life. She and I have been together through all the ups and downs of both undergraduate and
graduate student life and without her, these last few years would have been unquestionably much
harder and boring. Without their support and encouragements, I would not have been able to
finish a Ph.D. successfully. This achievement belongs to them as much as it does to me.
To my Ph.D. advisor, Professor Yehuda Ben-Zion, I owe an enormous amount of gratitude.
Throughout my Ph.D., he always found time for me and kept me positive and motivated when I
lost my focus from work. As a senior scientist, he taught me how to do research more accurately
and concisely that one should always keep the big picture in mind. As an elder, he is kind and
always shares personal experiences and words of wisdom with me. Because of his generosity, I
was able to gain more than 15 pounds of weight and attend many conferences in the past 6 years.
Thanks to him and all the graduate students who were and are in the research group (Lei, Pieter,
Zach, Xin, Yaman, Shiqing, Amir, Haoran, Yifang, and Niloufar) and USC earth sciences
department (Gen, En-Jui, Guang-Sin etc.), it makes me feel like living in a family although
thousands of miles away from home.
The list of people whom I have collaborated with during my Ph.D. life is long and I
appreciate and thank them all for their help and feedback. The ones most directly involved in my
work were: Yehuda Ben-Zion, Fan-Chi Lin, Gregor Hillers, Frank Vernon, Zachary Ross, Pieter-
Edward Share, Amir Allam, Gen Li, Lei Qin, and Elizabeth Berg. I have gained a lot from all of
you. I am grateful that Charlie Sammis, Francis Wu, Joshua West, and Krzysztof Pilch took the
time out of their tight schedules to serve on both my qualifying exam and dissertation defense
committees. Thanks also to the USC staff for making life much easier for us students. John
McRaney, John Yu, Cynthia Waite, Vardui Ter-Simonian, Tran Huynh, Deborah Gormley,
Karen Young, Barbara Grubb, and Miguel Rincon, you guys are the best.
Lastly, to my friends and family in the US and China. None of this would have been
possible without you. Thank you, everyone.
2
Table of Contents
Acknowledgements ................................................................................................................. 1
Abstract .................................................................................................................................. 4
Introduction ............................................................................................................................ 6
1. Eikonal tomography of the southern California plate boundary region .............................. 10
1.0 Preamble .................................................................................................................................. 10
1.1 Introduction .............................................................................................................................. 10
1.2 Data .......................................................................................................................................... 14
1.3 Automated dispersion picking ................................................................................................... 16
1.3.1 Frequency Time Analysis (FTA) ................................................................................................. 16
1.3.2 Determination of phase traveltime dispersions ....................................................................... 18
1.3.3 Determination of group traveltime dispersions ....................................................................... 20
1.3.4 Quality control .......................................................................................................................... 22
1.4 Eikonal tomography .................................................................................................................. 23
1.4.1 Methodology ............................................................................................................................ 23
1.4.2 Isotropic phase and group speed maps .................................................................................... 28
1.5 Shear wave velocity inversion ................................................................................................... 30
1.5.1 Methodology ............................................................................................................................ 30
1.5.2 3-D shear wave velocity model ................................................................................................. 33
1.6 Discussion & Conclusions .......................................................................................................... 38
2. Internal structure of the San Jacinto fault zone at Jackass Flat from data recorded by a
dense linear array ................................................................................................................. 41
2.0 Preamble .................................................................................................................................. 41
2.1 Introduction .............................................................................................................................. 41
2.2 Data .......................................................................................................................................... 44
2.3 Analysis .................................................................................................................................... 45
2.3.1 Teleseismic arrival time delay analysis ..................................................................................... 46
2.3.2 Local P wave Delay Time Analysis (DTA) ................................................................................... 52
2.3.3 Fault Zone Head Wave (FZHW) ................................................................................................. 57
2.3.4 Fault Zone Trapped Wave (FZTW) ............................................................................................ 64
2.4 Discussion ................................................................................................................................. 69
3. Internal structure of the San Jacinto fault zone at Dry Wash from seismic data of a linear
array ..................................................................................................................................... 73
3.0 Preamble .................................................................................................................................. 73
3.1 Introduction .............................................................................................................................. 73
3.2 Data .......................................................................................................................................... 76
3.3 Methods and Results ................................................................................................................ 77
3.3.1 Teleseismic earthquake delay time .......................................................................................... 77
3.3.2 Local earthquake delay time ..................................................................................................... 80
3.3.3 Fault zone head waves (FZHW) ................................................................................................. 80
3.3.4 Fault zone trapped waves ......................................................................................................... 82
3.4 Discussion ................................................................................................................................. 84
3
4. Temporal changes of seismic velocities in the San Jacinto Fault zone associated with the
2016 Mw5.2 Borrego Springs earthquake .............................................................................. 86
4.0 Preamble .................................................................................................................................. 86
4.1 Introduction .............................................................................................................................. 86
4.2 Data and processing .................................................................................................................. 90
4.3 Methodology ............................................................................................................................ 93
4.3.1 Relative Velocity Change Estimates .......................................................................................... 93
4.3.2 Curve fitting .............................................................................................................................. 97
4.4 Results ...................................................................................................................................... 99
4.4.1 Results of curve fitting ............................................................................................................ 101
4.4.2 Comparative analyses ............................................................................................................. 106
4.5 Discussion and Conclusions ..................................................................................................... 112
Discussion ........................................................................................................................... 118
A1. Chapter 1 supplementary figures .................................................................................. 121
A2. Chapter 2 supplementary figures .................................................................................. 138
A3. Chapter 3 supplementary figures .................................................................................. 146
A4. Chapter 4 supplementary figures .................................................................................. 149
References .......................................................................................................................... 167

4
Abstract
I developed a set of novel seismologically-based tools to image the San Jacinto Fault Zone
(SJFZ) and other parts of the plate boundary in southern California. The thesis includes analyses
of surface and coda waves reconstructed from seismic interferometry using ambient noise
records, P and S waveforms from local and teleseismic earthquakes, and fault zone head,
reflected, and trapped waves. Joint analysis of these signals reveals high-resolution spatio-
temporal information about bimaterial fault interface (i.e. the interface between two types of
rock), fault damage zones, and the neighboring upper crust.
A modified version of Eikonal tomography for imaging large-scale S-velocity (Vs) structure
is developed. The method is applied to data recorded by an irregular seismic network in southern
California. Reliable surface wave phase and group dispersion curves for 2.5-16 s are extracted
from ambient noise cross-correlations (ANC) and then used to derive isotropic surface wave
speeds. Joint inversion of phase and group velocities is performed to obtain local 1D Vs models,
which significantly improves the velocity structure in the top 3-20 km, particularly near faults.
Investigation of the non-uniqueness of the 1D Vs inversion suggests that imaging the top 3 km
Vs structure require either shorter period (≤ 2s) surface wave dispersion measurements or other
types of datasets such as Rayleigh wave ellipticity.
Slowness analysis of P waves from teleseismic and local seismic events is proposed for
understanding the internal structure of the fault zone at a local scale (i.e. 10-100 m). The
technique is applied on data recorded by dense linear arrays in the trifurcation area of the SJFZ at
Jackass Flat (JF) and Dry Wash (DW). The two arrays have 9-12 stations each with instrument
spacing of 25-100 m. Delay times between phase arrivals associated with all available local
earthquakes and teleseismic events are used to estimate P-velocity (Vp) variations within the
arrays. The statistically derived delay time pattern is utilized to distinguish the type of local fault
zone structure, a low-velocity damage zone and/or bi-material interface.
Fault zone head waves (FZHW) are emergent first arrivals that propagate exclusively along
fault bimaterial interfaces and only radiate energy to the low-velocity side. Such phase can be
used to constrain fault geometry and sharp across-fault structural property (e.g. Vp and Vs)
contrasts at depth. Existing automatic algorithm is first used to flag preliminary FZHW
candidates and visual inspections are later applied for verification. The detected FZHW at JF and
5
DW provide an unprecedented resolution on seismic wave velocity contrasts across the studied
fault interface.
Fault zone trapped waves (FZTW) result from constructive interference of reflected phases
within a sufficiently uniform low-velocity fault zone layer (i.e. fault damage zone). S-type
FZTW are systematically detected using data recorded by linear fault zone arrays at DW and JF,
and modeled for fault damage zone properties (i.e. depth and width, Vs, and attenuation
coefficient) as a waveguide under a 2D simplification.
Temporal changes of seismic velocities in the trifurcation area of the SJFZ are studied in
response to the 2016 Mw5.2 Borrego Springs earthquake. Ambient noise cross-correlations
computed for station pairs in JF and DW arrays, crossing the fault surface trace ~3 km northwest
(DW) and southeast (JF) of the event epicenter, are analyzed. Relative velocity changes (δv/v)
are estimated from arrival time changes in the daily correlation coda waveforms compared to a
reference stack. The obtained array-average δv/v time series exhibit seasonal variations, linear
trends and changes associated with the Borrego Springs event. The earthquake signal is
characterized by a rapid co-seismic velocity drop followed by a gradual recovery. The larger
changes at lower frequencies imply that the variations are not limited to the near subsurface
material. A decreasing co-seismic velocity reduction with coda wave lapse time indicates larger
co-seismic structural perturbations within the fault zone compared to the surrounding host rock.
The observed larger changes at the DW array compared to the changes at the JF array may
reflect larger damage/perturbation produced by the northwestward rupture directivity of the
Borrego Springs earthquake.

6
Introduction
Large continental earthquakes have been threatening communities in regions controlled by
seismogenic fault systems. Understanding the occurrence and mechanism of large earthquakes
requires detailed delineation of the subsurface structures and material properties of seismogenic
fault zones. However, limited understanding of the subsurface conditions of fault systems has
been obtained, mostly due to the lack of field observation datasets with sufficiently high spatio-
temporal resolution and effective tools that can “image” the subsurface. To fill this knowledge
gap, I develop a set of advanced seismological tools to characterize and monitor fault zones in
southern California (SC) at multiple length scales (hundreds of kilometers to tens of meters) by
taking advantages of the improved data acquisition, multi-dimensional seismological data sets,
detailed high quality earthquake catalogs, and state-of-the-art automatic detection algorithms of
seismic phases (Ross and Ben-Zion, 2014, 2015; Ross et al., 2016).
The San Jacinto fault zone (SJFZ) is one of the most seismically active fault zones in SC
(Hauksson et al., 2012; Ross et al., 2017a), and is capable of producing large Mw > 7.0
earthquakes (e.g. Petersen and Wesnousky, 1994; Onderdonk et al., 2013; Rockwell et al., 2015).
Several Mw > 5 earthquakes occurred in the SJFZ during the past decade. The most recent 2016
Mw5.2 Borrego Springs (BS) earthquake had a hypocenter at ~12 km depth in the complex
trifurcation area of SJFZ with three major fault branches. Five linear arrays that consist of ~ 10
stations with an average station spacing of ~ 20-40 m were deployed across different sections of
the SJFZ (Vernon and Ben-Zion, 2010) to study the internal structure of fault zones at different
locations along the fault strike. In addition to those local linear arrays, the plate boundary region
of SC is also well instrumented at a regional scale with typical station spacing varying from 1 km
to 30 km covered by several seismic networks, including the Anza network (FDSN), the
California Integrated Seismic network (CISN), and Southern California Earthquake Data Center
(SCEDC). Based on these seismic stations, several detailed high-quality catalogs have been
constructed for the area (e.g. Hauksson et al., 2012; White et al., 2019). Seismicity patterns in the
trifurcation area of the SJFZ indicate orthogonal structures, such as the aftershock distributions
of the 2016 BS event observed by Ross et al. (2017a) and other Mw > 5 earthquakes that
occurred in the region (Cheng et al., 2018), suggesting complex fault structures and geometry
rather than a simple vertically dipping plane.
7
Fault zone structures contain important information on various aspects of earthquake and
fault mechanics ranging from long-term evolutionary processes to brittle rock rheology and
dynamic stress fields operating during the occurrence of earthquakes (e.g. Ben-Zion and
Sammis, 2003; Xu et al., 2012; Rowe and Griffith, 2015). Some elements in the core structure of
fault zones can control (and reflect past) properties of earthquake ruptures. Specifically,
bimaterial interfaces separating different crustal blocks can influence significantly the location,
mode, speed and propagation direction of earthquake ruptures (e.g. Weertman, 1980; Ben-Zion
and Andrews, 1998; Brietzke and Ben-Zion, 2006; Ampuero and Ben-Zion, 2008; Brietzke et al.,
2009; Shlomai and Fineberg, 2016), and hence generated ground motion (e.g. Olsen et al., 2006;
Kurzon et al., 2014). A low-velocity fault zone layer can produce oscillations of rupture velocity
and slip rate (e.g. Ben-Zion and Huang, 2002; Huang et al., 2014) and promote rupture
propagation (Weng et al., 2016). Asymmetric rock damage across the seismogenic fault may
indicate repeating occurrence of large earthquake ruptures with statistically preferred
propagation direction (e.g. Ben-Zion and Shi, 2005; Lewis et al., 2005; Dor et al., 2006, 2008).
The surface expressions of fault zones often have high geometrical complexity with hierarchical
damage zones and multiple surface traces. This surface complexity and the diffuse character of
low magnitude seismicity make it difficult to identify the main seismogenic fault.
Several local and regional tomographic models of SC with a variety of spatial scales (0.5 km
to the entire SC) and resolutions (tens of meters to >5 km) have been developed in previous
studies. Methods using fault zone phases (e.g. Qiu et al., 2017; Share et al., 2017; Qin et al.,
2018), body wave travel times (e.g. Allam et al., 2012, 2014a; Share et al., 2019), surface waves
(e.g. Zigone et al., 2015; Barak et al., 2015; Berg et al., 2018), joint inversion of body and
surface waves (e.g. Fang et al., 2016) and full inversion of low-pass-filtered waveforms (e.g.
Tape et al., 2010; Lee et al., 2014) have been applied to data recorded in this region. These
studies imaged with nominal resolutions of 1–5 km several key structural features including
velocity contrasts across the SJFZ and damage zones at different locations. Each of these models
resolves different components of the crustal structures due to variations in data sensitivity and
quality (e.g. uncertainty), non-uniqueness and parameterization of the inversion process (e.g.
regularization and smoothing). A combination of tomographic models that incorporates different
types of data and seismic network configurations not only provides a clear illustration of the
8
structural complexity in the SC region (e.g., Shaw et al., 2015) but also plays an important role in
comprehensive physics-based seismic hazard evaluations (e.g. CyberShake from SCEC).
Ambient noise correlation (ANC) techniques have been used extensively to approximate the
Green’s function of the subsurface between two stations (e.g. Campillo and Paul, 2003; Shapiro
and Campillo, 2004; Shapiro et al., 2005; Zigone et al., 2019). Noise-based surface wave
tomography has been shown to be effective in resolving 3D Vs crustal structure either for the
entire SC region (e.g. Barak et al., 2015) or more focused local areas (e.g. Zigone et al., 2015;
Fang et al., 2016). Although the noise sources are not isotropically distributed in the SC plate
boundary region, the biases in surface wave dispersions measured from ambient noise cross-
correlations (ANC) have been shown to be minor (e.g. Hillers et al., 2013). However, for studies
monitoring seismic velocities (e.g. Kedar and Webb, 2005; Stehly et al., 2006), biases in the
reconstructed waves due to variations in the noise source distribution are non-negligible. Code
Wave Interferometry (CWI) uses coherent scattered arrivals in the coda part of the cross-
correlation (e.g. Snieder et al., 2002; Grêt et al., 2006; Sens-Schönfelder and Wegler, 2006) and
is less sensitive to biases associated with changing source properties (Colombi et al., 2014). The
recent advances in deployments of dense seismic arrays around fault zones and analysis
techniques (ANC, fault zone phases, CWI) provide not only sharper images of the fault zone
structure with a higher resolution but also improved opportunities to track the temporal evolution
of seismic velocities within active fault zones.
In Chapter 1, we use the ambient noise cross-correlations computed from stations of the
regional SC seismic network and apply Eikonal tomography to resolve phase and group velocity
maps. The resulting isotropic phase and group speeds are used to infer a new Vs model for the
SC plate boundary region. The final surface wave phase speed and Vs models are compared to
tomographic models obtained from previous studies, and the prominent geological structures that
observed in our models are discussed in section 1.5.2. The technical improvements and updated
geophysical knowledge achieved in our final models are summarized in Section 1.6. In addition
to resolving better some crustal components, the results complement the existing knowledge on
large-scale fault structures (velocity contrasts, fault dipping) in the SC plate boundary region.
Although the resulting Vs inversion is often not capable of resolving structures in the top few
kilometers (e.g. 3 km; Chapter 1), we showed that these shallow structures can be well
constrained by the Rayleigh wave ellipticity (Berg et al., 2018).
9
While the large-scale tomography provides a smoothed picture of fault structures, sharper
images of local fault zone structures can be achieved by analyzing data from fault zone linear
dense arrays. In Chapters 2 and 3, we image internal structures of the SJFZ at Jackass Flat (JF)
and Dry Wash (DW) based on delay times of P waves at stations across the arrays, and analyses
of fault zone head and trapped waves. Teleseismic P waves with low frequencies (e.g. < 1 Hz)
provide larger lateral and horizontal scales in structure information than the high-frequency (> 5
Hz) P waves from local earthquakes. Whereas fault zone phases that propagate exclusively along
fault bimaterial interfaces, fault zone head waves (FZHW), and in low-velocity damage zones,
fault zone trapped waves (FZTW), provide an unprecedented resolution on seismic wave
velocity contrasts across the studied fault interface, and properties within the fault damage zone,
respectively. The results indicate that the main seismogenic fault is close to the SW end of the JF
array but crossing the center of the DW array. The existence of a significant shallow local low-
velocity zone NE of the fault at JF array, with Q of ~20 and ~35% reduction in shear wave
velocity relative to the host rock, produces a local reversal to the overall large-scale velocity
contrast (Chapter 2). In contrast, delay time analysis at the DW array and lack of clear trapped
waves suggest a more localized fault zone structure (Chapter 3).
Chapter 4 documents temporal evolutions of seismic velocities in the trifurcation area of the
SJFZ in response to the 2016 Mw5.2 Borrego Springs (BS) earthquake. Two dense linear arrays
crossing the fault surface trace ~3 km northwest (DW) and southeast (JF) of the event epicenter
are analyzed. Co- and post-seismic relative velocity changes (δv/v) associated with the 2016 BS
event are observed consistently using both time- and frequency-domain analysis methods on data
of different components in various frequency bands at both arrays. The results suggest a much
larger co-seismic perturbation of seismic velocities within the fault zone compared to the
surrounding host rock. The observed frequency dependence of δv/v may indicate seismic velocity
changes that are associated with different mechanisms and depth sections. The resulting co-
seismic velocity reductions show systematically larger values at DW array for both low (0.1-0.4
Hz) and high (1-4 Hz) frequency bands, which may be associated with the northwestward
rupture directivity of the 2016 BS earthquake.
Following these chapters is a discussion on the integration of results from different chapters,
which gives a more comprehensive understanding of structural properties of the SJFZ and the
surrounding crust. I conclude with potential topics for future research in the study region.
10
1. Eikonal tomography of the southern California plate boundary
region

1.0 Preamble
This chapter is currently in review as: Qiu, H., Ben-Zion, Y., and Lin, F.-C., 2019a, Eikonal
tomography of the Southern California plate boundary region: Journal of Geophysical Research:
Solid Earth.
I was the main author, computed the ambient noise cross-correlations using all available
continuous waveforms of the year 2014 recorded by seismic stations in southern California,
analyzed the reconstructed surface wave signals, performed the shear wave velocity inversion,
and wrote the manuscript. Fan-Chi Lin provided the main supervision on the methodology.
Specifically, he advised on the technical details in developing the improved procedures of
Eikonal tomography. Yehuda Ben-Zion advised on the outline of the manuscript and provided
suggestions on the interpretation of the resulting velocity models. Both co-authors helped
improve the submitted manuscript.
The seismic data used in this work were obtained from the Southern California Earthquake
Data Center [SCEDC, 2013]. This study was supported by the Earthquake Hazards Program of
the USGS (grant G16AP00105) and the National Science Foundation (grant EAR-1841315).

1.1 Introduction
The boundary between the North American and Pacific plates in Southern California (SC)
has several major faults, including the San Andreas Fault (SAF), San Jacinto Fault (SJF),
Garlock Fault (GF), and Elsinore Fault (EF). These and other faults separate the SC crust into
several distinctive geologic provinces including the Southern Central Valley, Sierra Nevada and
Mojave Desert in the north; the Western, South-Central and Eastern Transverse Ranges plus
Eastern California Shear Zone (ECSZ) in the center; and the Ventura Basin, Los Angeles (LA)
Basin, Peninsular Ranges, and Salton Trough in the south (Fig. 1.1). Having a good 3D
tomographic model is crucial for understanding structural properties such as continuity and
dipping of the major faults, and providing an accurate framework for inversions of earthquake
source properties, calculations of seismic ground motion and other topics.
11
Several local and regional tomographic models of SC with a variety of spatial scales (0.5 km
to the entire SC) and resolutions (tens of meters to >5 km) have been developed in previous
studies. Methods using fault zone phases (e.g. Qiu et al., 2017; Share et al., 2017; Qin et al.,
2018), body wave travel times (e.g. Allam et al., 2012, 2014a; Share et al., 2019), surface waves
(e.g. Zigone et al., 2015; Barak et al., 2015; Berg et al., 2018), joint inversion of body and
surface waves (e.g. Fang et al., 2016) and full inversion of low-pass-filtered waveforms (e.g.
Tape et al., 2010; Lee et al., 2014) have been applied in this region. Each of these models
resolves different components of the crustal structures due to variations in data sensitivity and
quality (e.g. uncertainty), non-uniqueness and parameterization of the inversion process (e.g.
regularization and smoothing). A combination of tomographic models that incorporates different
types of data and seismic network configurations provides a clear illustration of the structural
complexity in the SC region (e.g., Shaw et al., 2015).

Figure 1.1. Location map of 346 (299 three-component sensor in red) seismic stations (triangles)
used for imaging the Southern California (SC) plate boundary region. Ambient noise cross-
correlations (ANC) computed at two example station pairs (green lines) are shown in Fig. 1.3.
The green square shows the location of the grid point used in Fig. 1.12. Surface traces of large
faults together with the state and national boundaries are shown as black lines. Localities of the
major faults and geologic provinces in SC are labeled. Cross sections of the final inverted shear
wave velocities (Vs) are shown for the blue lines crossing San Andreas Fault at various locations
in Fig. 1.17.
12
Noise based surface wave tomography has been shown to be effective in resolving 3D Vs
crustal structure either for the entire SC region (e.g. Barak et al., 2015) or more focused local
areas (e.g. Zigone et al., 2015; Fang et al., 2016). Although the noise sources are not
isotropically distributed in the SC plate boundary region, the biases in surface wave dispersions
measured from ambient noise cross-correlations (ANC) have been shown to be minor (e.g.
Hillers et al., 2013). Surface waves are typically assumed to propagate on the great circle path
connecting the virtual source and receiver, and the corresponding velocity structures are resolved
using all available source-receiver pairs (Barmin et al., 2001).
In contrast to the conventional surface wave tomography method, Eikonal tomography
accounts for ray bending and determines phase velocities by solving the Eikonal equation across
phase travel time maps (Lin et al., 2009). Through statistical analysis of velocity measurements
obtained from different sources, isotropic phase speeds together with azimuthal anisotropy and
corresponding uncertainty estimates can be determined. Eikonal tomography method was first
applied across the USArray (Lin et al., 2009), and the derived isotropic phase results were shown
to be slightly slower, on average, compared to those from straight-ray tomography, particularly
in SC. The differences suggest that it is necessary to take the ray bending effect into
consideration in order to obtain better phase velocity estimates.
In the present paper, we use the ANC computed from stations of the regional SC seismic
network and apply Eikonal tomography to resolve phase and group velocity maps. The results
have finer grid size (i.e. 0.05°×0.05°), broader period range toward short periods (i.e. 2.5s-16s),
and better data coverage compared to previous studies (e.g. Lin et al., 2009; Roux and Ben-Zion,
2017). In section 1.2, we describe the data used in the study and the necessary processing steps to
calculate reliable ANC for each station pair. Considering the vast number of station pairs, we
adopt the modified automatic Frequency Time Analysis (FTA) developed by Bensen et al.
(2007) and Lin et al. (2009) to obtain both phase and group dispersion measurements. Following
the flow chart shown in Figure 1.2, we describe the procedures to measure the surface wave
travel times as a function of period for every station pair in section 1.3.
In section 1.4.1, we review the methodology of Eikonal tomography and its underlying
assumptions. Different from previous Eikonal tomography studies (e.g. Lin et al., 2009, 2013;
Ritzwoller et al., 2011; Gouédard et al., 2012; Xu et al., 2016) that use evenly spacing seismic
arrays, the station spacing for the SC seismic network is rather irregular, with ~1-5 km near SJF
13
and SAF, ~5 km in LA and Ventura basins, and ~10-30 km in Mojave desert and ECSZ (Fig.
1.1). We thus also discuss the inclusion of several additional quality control criteria to ensure the
reliability of the resulting phase velocities. Furthermore, we derive group velocities from a
modified Eikonal tomography procedure, which is also discussed in section 1.4.1. The resulting
isotropic phase and group speeds (section 1.4.2) are used to infer a new Vs model for the SC
plate boundary region.

Figure 1.2. Flow chart of the procedures to obtain shear wave velocity model from the one-year
stacked ANC. NT, FA and PT are short for Negative Time lag, Fold and Average, and Positive
Time lag of the stacked ANC.

The Vs inversion is first performed at each grid cell and then assembled together to
construct the final 3-D Vs model in section 1.5. The final surface wave phase speed and Vs
models are compared to tomographic models obtained from previous studies, and the prominent
geological structures that observed in our models are discussed in section 1.5.2. The technical
14
improvements and updated geophysical knowledge achieved in our final models are summarized
in Section 1.6. In addition to resolving better some crustal components, the results complement
the existing knowledge on large-scale fault structures (velocity contrasts, fault dipping) in the SC
plate boundary region.

Figure 1.3. Daily ANC for the entire year 2014 computed at Vertical – Vertical (ZZ) component
of the (a) coast-parallel pair DJJ–GOR and (b) coast perpendicular pair GSC–SDD. Red, green
and blue represent positive, zero and negative amplitude values, respectively. Rayleigh wave
signals (bounded by dashed lines) at both positive and negative time lags are zoomed in and
shown in the insets. (c) One year stacked cross-correlation at components of ZZ, TT, RR, ZR,
ZT, and RT computed at station pair DJJ-GOR. (d) Same as (c) for pair GSC-SDD. Noise source
directionality is clearly observed in both pairs and for all components.

1.2 Data
We download all available continuous waveforms recorded by 346 stations (with 299 being
three-component; Fig. 1.1) in the SC plate boundary region during the entire year of 2014.
Seismic stations from several SC seismic networks, including the Anza network (FDSN), the
California Integrated Seismic network (CISN), Southern California Earthquake Data Center
15
(SCEDC) are used. This combined seismic network consists of 5 types of seismic instruments,
with 238 being broadbands and 108 being short period sensors, covering the ~600 km aperture
study region with typical station spacing varying from 1 km to 30 km.
Noise preprocessing steps are essential to increase the quality and accuracy of surface wave
signals extracted from the noise cross-correlation method (e.g. Shapiro and Campillo 2004;
Bensen et al., 2007; Poli et al., 2012; Boue et al., 2013; Zigone et al., 2015). In this study, we
follow closely to the method described by Zigone et al. (2015) to compute daily ANC between
all available station pairs and components. The computed multi-component ANC are then rotated
to the coordinate system of vertical (Z), radial (R), and transverse (T) directions by viewing one
station as the source and the other as the receiver (Lin et al., 2008). Figs. 1.3a&b shows the
resulting daily ANC for ZZ component at two example station pairs – the coast parallel pair DJJ-
GOR and the coast perpendicular pair GSC-SDD. For both station pairs, coherent surface wave
arrivals are observed on both positive and negative correlation time lags in the daily cross-
correlograms throughout the year. These daily correlation functions are stacked over the entire
year to further enhance the signal to noise ratio (SNR; Figs. 1.3c&d). In this paper, we use the
stacked ZZ and TT component ANC to measure the Rayleigh and Love wave travel time
dispersions, respectively. Here, as only the fundamental mode surface waves are reliably
observed in ANC, all results and discussions are referring to the fundamental mode surface
waves.
In SC, ambient noise seismic waves are mostly excited from oceanic waves in the
southwestern direction (e.g. Kedar and Webb, 2005; Hillers et al., 2013), which results in the
asymmetry of ANC particularly for coastline normal station pairs (e.g. Figs. 1.3b&d). Despite
the apparent noise directionality, earlier studies suggest that surface wave dispersion can still be
reliability extracted from ANC in this area (e.g. Hillers et al., 2013). To further enhance the
signal and effectively homogenize the noise wavefield, we calculate the “symmetric signal” by
folding and averaging the waveforms on both the positive and negative time lags (e.g. Lin et al.,
2007). We note that the Eikonal tomography approach used in this study determines local surface
phase velocities based on relative travel time measurement and is also less sensitive to
inhomogeneous noise source distribution as discussed in Lin et al. (2013).
16

Figure 1.4. The top black trace shows the one-year stacked ZZ component correlation function
recorded at station pair GSC-SDD. The corresponding symmetric signal, by folding and
averaging (FA) the positive and negative time lags, is displayed in red. The symmetric signal is
then filtered at periods 2s, 3s, 5s, 7s, 10s, 15s, and 20s, and the filtered waveform and
corresponding envelope are shown in blue and black, respectively. The surface wave window is
defined as an average velocity range of 1.5 km/s to 4.5 km/s, whereas an average velocity less
than 1.5 km/s outlines the noise window. Signal to noise ratio (SNR) is calculated for each
envelope global peak (red star). Reference phase traveltime dispersion calculated using the
CVM-S4.26 is illustrated as the red dashed curve. The blue star shows the location of a local
maximum of the envelope filtered at 2s.

1.3 Automated dispersion picking
Figure 1.4 shows the one-year stacked ZZ component ANC and the symmetric narrow
bandpass filtered signals for the example coastline normal station pair GSC-SDD. Clear period
dependent travel time and SNR can be observed. Considering the vast number of ANC (~ 40,000
for ZZ component and ~ 30,000 for TT component), we adopt a modified automated dispersion
picking algorithm of Bensen et al. (2007) to extract surface wave travel times for periods
between 2s and 20s. The procedures are described in detail below.’

1.3.1 Frequency Time Analysis (FTA)
17
Figure 1.5 illustrates the standard procedures of FTA applied on the example station pair
GSC-SDD. First, we taper the time series using a window (dashed lines in Fig. 1.5a) that outlines
the surface wave signal (i.e. between the assumed 4.0 km/s and 1.5 km/s maximum and
minimum velocities) and define the waveforms in the window with an assumed velocity lower
than 1.5 km/s as the noise. Then, a series of Gaussian narrow bandpass filters centered on
different angular frequencies, ω
k
, G ω,ω
k
( )
=exp −α ω −ω
k
( )
ω
k
⎡
⎣
⎤
⎦
2
{ }
, are applied to the
tapered waveform. Here α is a unitless parameter that controls the width of the Gaussian filter,
which we set to 20 based on trial and error. Then the amplitude and phase components of the
filtered signal, S
f
t,ω
k
( ) , can be written as:

S
f
t,ω
k
( )
= A t,ω
k
( )
⋅e
iϕ t,ω
k
( )
, (1)
where A t,ω
k
( )
, ϕ t,ω
k
( ) are the corresponding envelope and phase functions in the time
domain, and t is the lapse time. The envelope and phase functions at 7s are illustrated in Fig.
1.5b. Figure 1.5c shows a 2-D amplitude diagram that aligns the envelope functions with respect
to the corresponding central periods T
k
=2π ω
k
. This 2-D amplitude diagram is later used to
determine travel time dispersion in section 1.3.2
For different station pairs, the period range in which we can extract good quality surface
wave signals can vary significantly. To determine the proper period range for FTA, we first set a
maximum cut-off period as T
max
=Δ 2c≈Δ 6 to satisfy the far field approximation (Bensen et
al., 2007). Here Δ is the interstation distance in kilometers, and c is the assumed reference phase
velocity and is set to be 3 km/s (Fig. A1.1). In the case of T
max
>20s , we set T
max
=20s . Then,
we calculate the preliminary period dependent SNR as the ratio between the maximum amplitude
of the envelope function and the root mean square amplitude within the noise window. For each
ANC, we only perform FTA for the period range in which the SNR is larger than 5.
18

Figure 1.5. Example of frequency-time analysis performed on the symmetric correlation
function shown in Fig. 1.4. (a). The symmetric correlation (black) is first tapered using a window
bounded by the moveout range of 4 km/s and 1.5 km/s. (b) The waveform after tapering and
filtering using a Gaussian narrow bandpass filter centered at period 7s is denoted by the blue
signal. Phase and envelope functions are calculated and shown in red and black, respectively.
The white star indicates the envelope peak with the corresponding traveltime showing as green
dashed lines. (c) Frequency-time diagram. After applying a series of Gaussian narrow bandpass
filters centered on periods from 2s to 20s on the tapered signal shown in (a), envelope functions
are arranged by the corresponding center periods. The amplitudes are illustrated as colors from
blue to red indicating values from 0 to the maximum. The envelope shown in (b) is depicted at
the white dashed lines. The red dashed curve denotes the reference phase traveltime dispersion
curve calculated using model CVM-S4.26. Local and global maximums of all the envelope
functions are shown as symbols of black plus & green circles, and red circles, respectively. Here
we discard any envelope maximums (black plus) that are below the black dashed lines.

1.3.2 Determination of phase traveltime dispersions
The surface wave phase travel time dispersion can be obtained from the phase function
ϕ t,ω
k
( )

using equation derived in Lin et al. (2008) by assuming the source phase ambiguity
19
term λ

equal to 0:

t
ph
ω
( )
= ϕ t,ω
k
( )
+ωt −
π
4
−N ⋅2π
⎡
⎣
⎢
⎤
⎦
⎥
ω . (2)
Here, N is the cycle skipping ambiguity term, which can be resolved using a reference dispersion
relation. Note that the instantaneous angular frequency ω =∂ϕ t,ω
k
( )
∂t can be slightly
different from the central angular frequency ω
k
of the applied filter. To resolve cycle skipping
ambiguity, N, we take advantage of the existing high-resolution Community Velocity Models for
SC (i.e. Shaw et al., 2015 – CVM-H15.1; Lee et al., 2014 – CVM-S4.26). First, we compute
synthetic phase travel time dispersion for each station pair based on the model CVM-S4.26. At
each grid point, we extract the 1-D P- and S-wave velocity depth profiles and calculate the
corresponding phase velocity dispersion curves (Herrmann, 2013). By compiling these phase
velocity dispersion curves at all locations, we construct a series of 2-D phase velocity maps at
different periods. For each period, we calculate the synthetic surface wave phase travel time for
every station pair based on the 2-D phase velocity map using the fast marching method (Sethian,
1996). The model predicted pairwise phase travel time dispersion curve is then used as a
reference to constrain the travel times measured through the ANC.
Ideally, phase travel time can be measured at any selection of lapse time t within the surface
wave window following equation 2. But, in practice, phase travel time is often evaluated at the
time of the envelope peak, t = t
g
ω
( ) , to guarantee a maximum SNR (e.g. Aki and Richards,
2002; Lin et al., 2008). Figure 1.4 shows the envelope functions and corresponding t
g
ω
( ) for
various periods (red star), and we find the global maximum at 2s period is abnormally fast (i.e.
faster than the signals observed at longer periods). This indicates that either the noise source
distribution is not sufficiently homogeneous or the signal is dominated by strong body waves or
higher mode surface waves at this period. In this case, instead of evaluating the phase travel time
at the global peak, the lapse time of another local maximum should be used (i.e. the blue star in
Fig. 1.4). Considering the vast number of station pairs used in this study, we develop an
automatic procedure to filter out these abnormal envelope peaks. First, we set a “reasonable
travel time range” for each station pair and rule out envelope peaks (both local and global ones)
that are outside the range as well as those peaks with SNR less than 5. Considering the envelope
20
should propagate slower than the phase, we use the phase travel time predicted from the
reference CVM-S4.26 model (red dashed lines in Fig. 1.4 and Fig. 1.5c) to determine the lower
bound for group travel time t
g
ω
( ) . In addition, considering the CVM-S4.26 model is derived
using data with frequencies lower than 0.2 Hz (Lee et al., 2014) and it may yield poor predictions
at short periods, we further reduce the lower bound based on a linearly varying scale of 5% and
0% between 2s and 20s period (black dashed curve in Figure 1.5c).
After filtering out erroneous envelope peaks, we further require the candidate t
g
ω
( ) (e.g. red
stars in Fig. 1.5c) to be continuous as a function of periods. To ensure the continuity, we use the
second order derivative edge detection algorithm to find possible jumps in the t
g
ω
( ) dispersion
curve, and fix the discontinuity, if detected, following the procedure of Bensen et al. (2007).

1.3.3 Determination of group traveltime dispersions
Theoretically, the surface wave group travel time is given by the lapse time, t
g
ω
( ) , where
the envelope function A t,ω
k
( )
reaches the maximum amplitude (Aki and Richards, 2002). In
Figure 1.5, we use the surface wave from the symmetric signal of the one-year stacked ANC to
demonstrate the picking process (Fig. A1.2 for results using causal and acausal sides). Figure 1.6
shows the resulting Rayleigh wave phase travel time along with the continuous t
g
ω
( ) dispersion
curves measured at the casual side, acausal side, and the symmetric signal of the correlation
function for station pair GSC-SDD. While consistent phase travel time dispersion curves are
obtained from all three cases, significant discrepancies are observed between the t
g
ω
( )
dispersion curves, especially at longer periods (> 15s), suggesting that the peak of envelope
function is sensitive to noise and often associated with large uncertainties (Fig. 1.6d and Fig.
1.7a). Therefore, in this study, instead of determining group travel time based on t
g
ω
( ) (e.g.
Barak et al., 2015; Zigone et al., 2015), we simply derive the group travel time from the
smoothed phase dispersion using the following relation:

υ
gp
=υ
ph
−λ∂υ
ph
∂λ
. (3)
21
where υ
ph
and υ
gp
represent the smoothed phase and resulting group dispersions in the form of
average velocities (i.e. υ
ph
=Δ t
ph
, υ
gp
=Δ t
gp
) and λ is the wavelength given by 2πυ
ph
ω .
Here, we use the observed phase dispersion to first invert for a 1D Vs model and then predict the
smoothed phase velocity dispersion to stabilize the first derivative in equation 3 (Fig. 1.7).

Figure 1.6. Rayleigh wave group and phase travel time dispersion results for example station
pair GSC-SDD. (a) Black solid curve represents the group travel time dispersions measured
using waveform at the symmetric signal. The corresponding phase travel time dispersion is
shown as the red solid curve. The blue dashed curve represents the model predicted phase travel
time dispersion using CVM-S4.26. Phase travel time dispersions with one cycle skipped
(N=N
0
±1 in eq. 2; red dashed curves) are shown for visual comparison. (b) Same as (a) measured
at the negative time lag. (c) Same as (a) using the positive time lag. (d) Comparison of all the
group (black dashed) and phase (red dashed) dispersion results. The blue and green solid curves
represent the final phase and group dispersion measurements.

Although the phase derived group travel times do not provide additional independent
constraints to the earth structure, using both phase and group dispersion curves stabilizes the 1-D
Vs inversion and yields better results of Vs structure than those using phase dispersions alone.
This is mostly due to the differences in depth sensitivity kernels of the group and phase velocities
22
(Fig. A1.3) that helps reducing non-uniqueness inherent in the inversion (Moschetti et al., 2010a;
Li et al., 2012; Herrmann, 2013).

1.3.4 Quality control
As Eikonal tomography determines velocities based on first order spatial derivative of travel
time maps, the result is very sensitive to erroneous travel time measurements (e.g. Lin et al.,
2009). Because of that, it is crucial to identify and remove as much erroneous travel time
measurements as possible. In this section, we introduce the following three selection criteria to
control the quality of the travel times measured following procedures in section 1.3.2 and 1.3.3:
(1) Consistent between symmetric, causal, and acausal signals. Here we reject dispersion
measurements with inconsistent phase dispersion between causal, a-causal, and the
symmetric component signals. Here we calculate the phase travel time differences between
measurements of the symmetric and causal components at each period, and select the period
range that yields discrepancies less than 5%. Such period range can also be obtained by
comparing phase dispersions of the symmetric and acausal signals. The phase dispersion
curves within the union of these two period ranges are considered to be robust.
(2) Consistent with the reference model. As no major discrepancy is expected between the
observed and reference predicted dispersion curves in particularly at the long periods, two
additional selection criteria are introduced for phase dispersion measurements. First, the
predicted and observed travel time difference at the longest measurable period (T
1
) should be
smaller than 0.3T
1
. Second, the average predict and observed travel time difference in the top
one-third of the measurable period range should be smaller than 0.4 of the average period in
that range. Dispersion curves that do not satisfy the above criteria are discarded.
(3) Minimum period range requirement. Since the predicted reference dispersion curve is
only considered reliable at periods larger than 5s (Lee et al., 2014), we reject all dispersion
curves with either the longest measurable period T
1
smaller than 6s or the measurable period
range shorter than 2.5s.

Figure 1.8 shows the histograms of the phase and group travel times that pass the above
quality selection criteria for 3s, 7s, and 11s Rayleigh and Love waves. Since there are less three-
component stations in the seismic network and the SNR of Love wave is generally smaller than
23
that of Rayleigh wave, the total number of travel time measurement is significantly higher (i.e.
~40%-50% more) for Rayleigh wave than Love wave. Based on the distributions of the
measurements, we find that the width of the histogram (i.e. standard deviation) decreases as
period increases, and the histograms of Love wave are wider than those of Rayleigh wave. These
observations are likely due to higher degree of laterally heterogeneity at shallower depth. Since
SNR is lower at shorter period, the poor quality of short period dispersion measurements also
contributes to the large width of the corresponding histogram. The median average speed (i.e.
peak of the histogram) increases with period, which is consistent

Figure 1.7. Derivation of group traveltime dispersion curve. Panel on the left shows the
measured phase traveltime dispersion curve (solid blue curve in Fig. 1.6d) in term of average
velocity as red dots. A 1-D Vs inversion is performed to fit the phase dispersion starting with the
1-D Vs profile (black curve in the right panel) averaged over the entire CVM-S4.26 as the initial
model. The phase dispersion curve (red curve in left panel) of the best fitting 1-D Vs profile (red
curve in right panel) gives the smoothed phase dispersion curve, and the corresponding group
dispersion (black curve in the left panel) is calculated following equation 3.

1.4 Eikonal tomography
1.4.1 Methodology
Different from the traditional straight-ray tomography (e.g. Barmin et al., 2001), Eikonal
tomography accounts for ray bending in the surface wave propagation and is based on the
Eikonal equation
24

!
s
ij
=
!
e
ij
υ
ij
=∇τ
i
!
x
j
( )
, (4a)
which is derived from the 2-D Helmholtz wave equation (e.g. Wielandt, 1993; Lin and
Ritzwoller, 2011) by neglecting the term associated with the amplitude Laplacian:

1
υ
ij
2
= ∇τ
i
!
x
j
( )
2
−
ΔA
i
!
x
j
( )
A
i
!
x
j
( )
ω
2
, (4b)
where i and j indicate the virtual source and grid point indexes,
!
s and
!
e are local slowness and
the unit vectors orienting towards the wave propagation direction, τ
i
is the phase travel time,
!
x
j

is location of the j-th grid cell, A is the wave amplitude, and ω is the angular frequency. In
Eikonal tomography, phase velocity structure can be simply inferred locally by applying the
inverse operator – the spatial gradient to the phase traveltime field without constructing the
forward operator. It is straightforward and computational less intensive compared to the straight-
ray tomography. Lin and Ritzwoller (2011) refers the phase velocity derived from Eikonal
equation as the “apparent” phase velocity and that calculated via Helmholtz equation as
“corrected” phase velocity. These two velocities are approximately equal when 1) the angular
frequency, ω , is sufficiently high or 2) the amplitude field is smooth enough so that the
amplitude Laplacian is negligible. Although the group travel time does not obey the Eikonal
equation 4a, based on the assumption that the propagation of the surface wave envelope is the
same as indicated by the phase front, we can apply a modified version of Eikonal equation:

!
s
ij
g
=
!
e
ij
υ
ij
g
=∇τ
i
g
!
x
j
( )
, (5)
to infer the local group velocity, where
!
s
ij
g
, υ
ij
g
, and τ
i
g
represent the local group slowness
vector, group speed, and group travel time field, respectively.
Following the procedure developed in Lin et al. (2009), for each common station, all
available phase or group travel times associated with the central station (Fig. A1.4a) are used to
construct a travel time map on a 0.05°×0.05° grid (Fig. A1.4b). The minimum curvature
interpolation method (Smith and Wessel, 1990) is used. Despite all the quality selection criteria
we developed in section 1.3.4, to obtain smooth travel time map, we impose additional quality
control criteria to remove outlier travel times that are not consistent with their nearby
measurements. Specifically, we reject any travel time measurement that meets any of the two
following conditions: 1. The amplitude of travel time gradient at the station location is less than
25
0.25 s/km or larger than 2 s/km (green circles in Fig. A1.4c). 2. The curvature of travel time
field at the station location is identified as an outlier when larger than 2 times the standard
deviation computed over the entire map (red circles in Fig. A1.4c). After removing the outliers,
new travel time maps are then regenerated and phase propagation direction and local phase and
group slowness can be evaluated through equation 4a (Fig. A1.4d) and equation 5 at each grid
point.

Figure 1.8. Histograms of phase (blue) and group (orange) traveltime measurements for
Rayleigh (top panels) and Love (bottom panels) waves at 3s (left panels), 7s (middle panels), and
11s (right panels). Total number of the traveltime measurements for each histogram is indicated
as well.

Unlike the Eikonal tomography performed on USArray by Lin et al. (2009), the southern
California seismic network used in this study is highly irregular (Fig. 1.1) with regions that has
various distinctive station spacing configurations: ~1-5 km near the major faults, and ~5 km
within basin areas, and ~10-30 km for other regions (e.g. ECSZ). Because of this uneven station
distribution, it is essential to identify regions with robust and reliable travel time interpolation,
and only estimate the travel time gradient within these areas. We first adopt the criteria used in
Lin et al. (2009), including truncation of regions that are within two wavelengths of the virtual
source location, “three- out of four-quadrant selection criterion” with a searching radius of 50
km, and comparison of phase travel time interpolation with two different surface tension (Figs.
26
A1.5a&b). In order to further tackle with the irregular station spacing configuration, we
introduce the concept of “station configuration error” and further remove regions that are
highlighted by large error values.

Figure 1.9. Eikonal phase velocity maps computed at period 7s by using stations (a) GOR, (b)
GSC, (c) IRM, and (d) OLI as the virtual source. Azimuths of the gradient are illustrated with
arrows.

Similar to the idea of a synthetic checkerboard test (Lévěque and Wittlinger, 1993), which
provides an estimation of the spatial resolvability for specific data coverage, we perform
synthetic tests to evaluate the station configuration error for Eikonal tomography. First, for each
virtual source, we compute the synthetic travel times for the same station configuration assuming
a homogeneous phase speed υ
syn
. Then we apply Eikonal tomography to the synthetic travel
times and obtain an estimated 2-D phase speed map with υ
i
inv
representing the local phase speed
at the i-th grid cell. The station configuration error δ
i
at the i-th grid point is defined as
δ
i
= υ
syn
−υ
i
inv
υ
syn
. (6)
Here we use a threshold of 0.025 to further truncate regions with poor station coverage (Fig.
A1.6c).
27
Figure 1.9 shows the resulting phase velocity maps for 7s Rayleigh wave using 4 different
stations as the virtual source. Colors and arrows represent the estimated local phase speed and
propagation direction, respectively. Patterns that match well with the surface geological feature,
such as low velocity zones in LA basin, Salton Trough, and near faults, are consistently observed
in all the maps. However, there are also differences found between these maps. Part of these
differences can be explained by azimuthal anisotropy. As shown in Lin et al. (2009), for a
weakly anisotropic medium (Smith and Dahlen, 1973), the 2-psi anisotropy yields a simple
relation between wave speed and the azimuth
υ ψ
( )
=υ
iso
2ψ
1+
A
2ψ
2
cos 2 ψ −θ
2ψ
( )
⎡
⎣
⎤
⎦
⎧
⎨
⎪
⎩
⎪
⎫
⎬
⎪
⎭
⎪
, (7)
where ψ is the angle of the wave propagation and υ
iso
2ψ
is the isotropic wave speed. A
2ψ
and
θ
2ψ
represent the peak-to-peak amplitude and fast axis of the anisotropy (Fig. A1.6),
respectively. In general, the phase velocity map based on individual effective source (Fig. 1.9) is
noisy due to erroneous phase travel times that we are unable to remove completely using the
quality selection criteria. However, previous Eikonal tomography studies (e.g. Lin et al., 2009,
2013) showed that these errors could be significantly suppressed through stacking.
Since a non-negligible azimuthal anisotropy effect is observed in this region (Fig. A1.6), to
avoid biases in the stacking process, for each location, we first weight each slowness
measurement s
i
inversely proportional to the number (n
i
) of measurements that has an azimuth
differs from the target measurement by less than 20°:
′ s
i
=
1
n
i
s
i
. (8)
Let normalization coefficients η=
1
n
i
i=1
N
∑
and ξ =
1
n
i
2
i=1
N
∑
, then the weighted mean slowness, s
0
,
and the corresponding standard deviation, σ
s
0
, are given by:
s
0
=
1
η
′ s
i
i=1
N
∑
, (9a)
σ
s
0
2
=
ξ
η
3
−η⋅ξ
⋅
1
n
i
s
i
− s
0
( )
2
i=1
N
∑
, (9b)
28
where N is the total number of slowness measurements available at the location from
different virtual sources.

Figure 1.10. Isotropic phase velocities (a-c) and the corresponding of uncertainty distributions
(d-f) of Rayleigh waves at 3s, 7s, 11s. The white lines in (a-c) indicate the boundaries between
distinctive geologic provinces (Hauksson and Meier, 2018).

1.4.2 Isotropic phase and group speed maps
Figures 1.10 and 1.11 show the resulting stacked isotropic phase and group speed maps for
3s, 7s, and 11s Rayleigh waves with corresponding uncertainty distributions. The isotropic speed
maps for Love waves are shown in Figs. A1.7&A1.8. Azimuthal anisotropy can also be derived
for each location at different periods based on equation 7 (e.g. Fig. A.6), but in this paper, we
only focus on the isotropic phase and group velocities. The anisotropy results will be discussed
in a separated study.
The isotropic phase and group speed maps at 3s (top left of Figs. 1.10, 1.11, A1.7, and A1.8)
agree well with the surface geology (e.g. white lines). Low velocity zones (LVZ) are observed at
Southern Central Valley, basins (e.g. LA basin, Ventura basin, San Bernardino basin), Salton
Trough, and complex fault junctions (e.g. SGP, the Trifurcation area in SJFZ). Higher velocities
(~3 km/s) are seen in regions such as the Peninsular Ranges. The LVZ below the basins and
Salton Trough show a flower-type structure (e.g. Zigone et al., 2015), width decreases with depth
or period. Clear velocity contrasts are found across surface traces of the major faults (e.g. SJF,
29
SAF). The Peninsular Ranges have the fastest velocity values of the entire map for all periods.
Consistent with the increasing histogram peak velocity with period shown in Fig. 1.8, faster
velocities are observed for the isotropic phase and group maps at longer periods. The obtained
Rayleigh wave phase and group speed maps are generally similar to results from previous studies
(Barak et al., 2015; Zigone et al., 2015; Roux and Ben-Zion, 2017) in the overlapping area, with
our phase and group velocities being slightly slower beneath basins (e.g. LA basin and Ventura
basin) and showing sharper velocity contrast across the SAF southeast to the SGP at shorter
periods (< 7s).

Figure 1.11. Same as Fig. 1.10 for isotropic Rayleigh wave group velocity.

The uncertainties, calculated using equation 9b, provide an estimate of the variability in
velocities derived from different virtual sources. Larger uncertainties are observed at shorter
periods. This may indicate the quality of isotropic speed maps is lower or the azimuthal
anisotropy (e.g. Fig. A1.6) is larger at the shorter period compared to those at the longer period.
In addition, while the spatial distribution of uncertainties is similar, larger values are also
observed for the group speed than the phase speed, as phase dispersion is intrinsically smoother
and more stable than group dispersion (eq. 3). In addition, we find larger uncertainty values at
the edge of the model (e.g. south of Salton Trough, Southern Central Valley), and this is due to a
poor azimuthal and station coverage. The uncertainty distributions for Love wave (Fig.
A1.7&A1.8) show a similar pattern but with generally larger values, as both the data quality is
30
poorer and the depth sensitivity is larger at shallower depths (Fig. A1.3) for Love wave than
Rayleigh wave.

1.5 Shear wave velocity inversion
1.5.1 Methodology
In section 1.4.2, 2-D isotropic phase and group speed maps for Rayleigh and Love waves at
different periods are derived. These period dependent isotropic phase and group speed maps can
be used to infer the Vs structure. In this study, we adopt the iterative 1-D Vs inversion scheme of
Herrmann (2013) to construct our final 3D shear wave velocity model for the SC plate boundary
region. We use the Southern California Community Velocity Models CVM-H15.1 as our
reference starting model.
In each of these 1-D Vs inversions, to avoid overshooting, we use a damping factor of 50 in
the first 3 iterations, and then lower it to 5 for another 20 iterations. In the process, we fix the
Vp/Vs ratio and Moho depth, and use the differential smoothing constraint to prevent unrealistic
(e.g. large jumps, oscillating-like) shape in the final inverted 1-D Vs profile. Once the final
inverted Vs profile is obtained, we correct the topography effect by simply subtracting the
elevation value from the depth of the Vs profile and assemble all the corrected 1-D Vs profiles to
construct the final pseudo-3D Vs model for the entire region. To evaluate the performance of
every 1D Vs inversion, we define misfit χ as
χ =
1
n
υ
i
obs
−υ
i
pred
( )
σ
i
obs
⎡
⎣
⎤
⎦
2
i=1
n
∑
, (10)
where n is the number of input dispersion data points, υ
i
obs
and υ
i
pred
denote the input and model
predicted dispersion wave speed for the i-th data point, σ
i
obs
is the corresponding data
uncertainty. A χ value less than 1 indicates, on average, the model predicted dispersion curves
fit the input dispersion curves within the corresponding uncertainties. Therefore, we set a
threshold of 2 to reject inverted Vs profiles with poor data fitting. In addition, we can also
compute the χ value for the initial model and compare it with that of the final model to estimate
the general variance reduction.
31

Figure 1.12. (a) Illustration of the iterative 1-D Vs inversion of Herrmann (2013) at an example
grid cell in San Gorgonio Pass located at -117°, 34°. The left panel shows the comparison
between the Rayleigh wave group (in blue) and phase (in red) velocity dispersion measurements
(solid circles) and the best fitting results (solid curves). The error bar indicates the uncertainty
estimated from eikonal tomography (eq. 9b). Rayleigh wave dispersion curves of the initial
model are also displayed as dashed curves. The χ misfit values for both the initial and best fitting
1-D Vs profiles are indicated at the top left corner. The black and red curves in the right panel
denote the initial (CVM-H15.1) and best fitting 1-D Vs models. An estimation of non-uniquness
of the 1-D inversion is illustruated by the gray shaded area given by Fig. 1.13d. The depth
dependent width of the gray shaded area is indicative of the inversion uncertainty and shown as
the gray curve. (b) Same as (a) for using CVM-S4.26 as the initial model. The blue curve in the
right panel represents the best fitting 1-D Vs profile obtained in (a).

Due to the limited period range (i.e. 2.5-16s) of the input surface wave dispersion curves, the
result of the inverted models can be somewhat sensitive to the reference starting model used, and
this sensitivity can vary with depth. Figure 1.12 shows the comparison between the inverted 1-D
Vs profiles obtained using CVM-H15.1 and CVM-S4.26 as the reference starting model at an
example grid point near the SGP. Despite the differences in the reference starting model, the
32
resulting misfits from both inversions are almost the same, and the inverted Vs profiles are also
consistent between 3 km and 25 km, suggesting the inverted result in this depth range is well
constrained by our data. However, large discrepancies are observed in the top 3 km and below
Moho depth, suggesting the Vs values are not well constrained beyond the 3-25 km depth range,
and thus are heavily biased by the initial model. In order to further quantify how the inverted Vs
values are constrained at different depths, we use an improved Neighborhood Algorithm (NA)
developed by Wathelet (2008) to assess the non-uniqueness of the 1-D Vs inversion. By
exploring the physical parameter space (i.e. layer thickness, Vp, Vs, and density), the NA
method can find a collective of models that fits the dispersion data within a given misfit range.
Here we use the variability of all these Vs models as a function of depth to infer the uncertainty.
Figure 1.13 shows an example of such exploration using NA at an example grid cell in the
SGP. In the example, we parameterize the 1-D model with 7 layers (6 layers with linearly
increasing Vs and flexible thickness + a bottom half space). After 200 iterations with 40,200
models being tested, we then select the series of 1-D Vs profiles with misfits less than 1.5 times
the lowest misfit of all models. The survived 1-D Vs profiles are shown as the gray shaded area
in Fig. 1.12 and the uncertainties at each depth can be estimated as the corresponding width of
the gray shaded area. Consistent with what we observed in the comparison between inverted
results using different initial models (Fig. 1.12a), the uncertainty shows much larger values at
shallow depth (i.e. top 3 km) and below 25 km. Consistent observations are seen at several other
grid points that are located at some of the representative geologic provinces (Fig. A1.9), thus in
later sections, we only focus our discussions on the depth range of 3-25 km. For the Vs structures
in the top 3 km, one can defer to either the Vs model presented in Berg et al. (2018), in which
shallow structure is better constrained by H/V ratio, or the geotechnical layer added in the CVM-
H15.1 or CVM-S4.26.
As the Vs depth sensitivity kernels for Rayleigh wave dispersions are different from those of
Love wave dispersions (Fig. A1.3), joint inversion of both Rayleigh and Love wave dispersion
curves could imply stronger constraint in the Vs inversion. However, Rayleigh and Love wave
velocity dispersions are sensitive to different Vs structures, V
SV
and V
SH
, respectively. Thus,
following Zigone et al (2015), we only perform the Vs inversion separately for Rayleigh and
Love waves to account for differences both in data quality and (apparent) radial anisotropy (Fig.
A1.10). Considering the quality of the isotropic Rayleigh wave phase and group speed maps is
33
much higher than those for Love wave, we only focus on the discussion of the Vs model from
jointly inverting the Rayleigh wave phase and group dispersion curves, while the Vs model
derived from Love wave data can be found in the supplementary material.

Figure 1.13. Illustration of Neighborhood Alogrithm (Wathelet, 2008) inversion results. The 1-D
Vp and Vs profiles explored in the inversion are colored according to their misfit, and those with
misfit values less than 1.46 are shown in (a) and (b). The corresponding group and phase velocity
dispersion curves are displayed in (c) and (d). Models with misfit larger than 1.5 times the
minimum misfit value (i.e. 0.41) are discarded, and the minimum and maximum of all the
acceptable 1-D Vs profiles at different depth depict the gray shaded area shown in Fig. 1.12.

1.5.2 3-D shear wave velocity model
34
Figure 1.14 summarizes results associated with the linearized 1-D Vs joint inversions
using Rayleigh wave phase and group dispersion curves. The misfit corresponding to the CVM-
H15.1 (Fig. 1.14a) and CVM-S4.26 (Fig. 1.14c) suggest that the surface wave dispersions
predicted by the initial models are, in general, inconsistent (i.e. large χ values) with the final
results obtained from the Eikonal tomography (Figs. 1.14b&d). The misfit distribution for the
final inverted Vs models is similar regardless of which initial model is used. We prefer the final
Vs model inverted using CVM-H15.1 as the initial model since: 1) Topography is considered in
CVM-H15.1 but ignored in CVM-S4.26; 2) The misfit histogram suggests a slightly smaller χ
value for CVM-H15.1 than CVM-S4.26; 3). Also, the 1-D Vs profiles of CVM-H15.1 are
generally simpler than those of CVM-S4.26 (e.g. Fig. 1.12).

Figure 1.14. Histograms of probability (in gray; PDF) and accumulative (blue curve; CDF)
density distributions for χ misfit. (a) χ misfit values computed for CVM-H15.1 following
equation 10 for all available grid cells. (b) Same as (a) for the best fitting Vs model using CVM-
H15.1 as the initial model. (c) Same as (a) for CVM-S4.26. (d) Same as (a) for the besting fitting
Vs model using CVM-S4.26 as the initial model. The corresponding spatial distributions of the
χ

misfit values are shown in Fig. A1.11.
35
Map views of the final Vs model and differences from the initial model CVM-H15.1 at
depths of 3 km, 5 km, 7 km, 10 km, 15 km, and 20 km are shown in Figs. 1.15&1.16. The largest
differences are observed underneath basins, such as the Salton Trough, and part of ECSZ for
depth between 3 km and 10 km. This is consistent with the fact that the χ misfit values are
significantly reduced in these regions (Fig. A1.11). In general, we observe the following
prominent features in relation to the surface geology (e.g. white lines):

Figure 1.15. Top panels show map view of the final inverted Vs model at 3 km (left), 5 km
(middle), and 7 km (right) depths. To show the Vs model better for regions near SJF and SAF,
the bottom left inset illustrates the Vs model in a different color palette. Model CVM-H15.1 is
used as the initial model here, and the bottom panels illustrate the differences in Vs between the
final and initial Vs models at 3 km, 5 km, and 7 km depths. White lines in the top panels
illustrate the boundaries of different geologic provinces (Hauksson and Meier, 2018).

1). Consistent with the surface geology and previous tomographic imaging results (e.g. Tape
et al., 2010; Lee et al., 2014; Berg et al., 2018), the Southern Central Valley is characterized by
LVZ and it changes rapidly to high velocity Sierra Nevada foothills. 2). The Ventura and LA
basins are well confined and highlighted by the slowest velocities of the entire map in the top 3-7
km. 3). A LVZ is observed within the junction southeast to SGP between SJF and SAF. 4). Clear
LVZ with a width of ~ 20 km is observed centered on the section of SAF southeast of the SGP.
This is consistent with results obtained by Share et al. (2019). Different from the initial model
(Fig. A1.12), this LVZ likely reflecting fault damage is still clearly observed up to a depth of 5
km. A LVZ is also observed around the SJF in the top 5 km. 5). Consistent with the observation
36
of velocity contrast reversal across the SAF northwest and southeast of the SGP (Share and Ben-
Zion, 2016), we observe a clear north-south orienting fast velocity block that cuts through the
SGP at 5 km and 7 km depth, leading to a flipping of the velocity contrast polarity across the fast
velocity anomaly. 6) The Vs around the ECSZ is much slower after the inversion, particularly for
the region north to SGP for depths from 5 km to 10 km, which corresponds well with the area
that has large damage volume in Ben-Zion and Zaliapin (2019). 7) The Salton Trough is imaged
with a well-defined shape of LVZ extended to depth 7 km in our final Vs model. Compared to
the initial model, the velocities are much slower (~ 0.3 km/s) in the top 3-7 km.

Figure 1.16. Same as Fig. 1.15 at depths of 10 km (left), 15 km (middle), and 20 km (right).

8) Different from the initial model, comparing Vs at 10 km and 15 km for our final Vs model
at the south SAF (Fig. 1.16), we see a clear shift in the velocity contrast interface location
indicating a northeast dipping fault plane. 9) The highest velocities are observed in the
Peninsular Ranges, and a sharp velocity contrast from west to east at greater depth (7-15 km;
white vertical line in the insets of Fig. 1.16) that corresponds to the Hemet step-over (Marliyani
et al., 2013) is observed much clearly in the final model. 10) Velocity contrasts across major
faults (e.g. SAF, SJF) previously imaged in other tomography (e.g. Fang et al., 2016; Share et al.,
2019) and fault zone head wave studies (Share and Ben-Zion, 2016, 2018) are observed clearly
in the final Vs model.
37

Figure 1.17. Cross sections of the final inverted Vs model at locations indicated as blue lines in
Fig. 1.1. Localities of major faults, basins, and geomorphic provinces are labeled on the top
topography curve. The red dashed lines beneath LA basin at profile CC’ denote a linear low
velocity zone that is likely associated with the Puente Hills blind-thrust system (Shaw et al.,
2002). In addition, a deep low velocity anomaly outlined by the black dashed circle at profile
DD’ may be related to the large damage volume estimated in Ben-Zion and Zaliapin (2019). The
black dashed lines at cross sections DD’, EE’, and FF’ denote the potential fault planes of SAF
or SJF.

Figure 1.17 presents Vs profiles for six cross-sections crossing the SAF at different locations
from the north (AA’) to south (FF’). In each profile, Vs structures between 3 km and 20 km are
displayed. These cross-sections show the following features: 1) The LA basin is deeper and
larger in our final Vs model that the initial model. 2) Beneath the LA basin, there is a low
velocity fault plane like block dipping towards the northeast. 3) The SJF is identified as a near
vertical dipping fault in DD’ centered on a localized LVZ. 4) Pronounced deep (15-20 km) low
velocity body is found beneath SGP, which is likely linked to large damage volume indicated in
the study of Ben-Zion and Zaliapin (2019). 5) The south SAF is found to be a localized fault
associated with a velocity contrast interface and dipping to the northwest. 6) Localized LVZ with
much slower (~0.2-0.3 km/s) velocities extend to 7 km beneath the ECSZ. These features are
seen much more clearly in our final Vs model compared to the CVM-H15.1 (Fig. A1.13).
The final Vs model using Love wave data is shown in Figs. A1.14-A1.19. Consistent features
in the Vs structures are observed, including as misfit histogram and distribution (Figs.
A1.14&A1.15), localized LVZ in the top 5 km related to fault damage (Fig. A1.16), fault like
structure with a dipping angle of ~30º beneath LA basin, and prominent low velocity body
38
located between 15-20 km depth beneath SGP (Fig. A1.17). Some discrepancies are found
between the derived Vs models from Rayleigh and Love waves, such as the observation of
dipping SAF is clear in the Rayleigh wave results (EE’ & FF’ in Fig. 1.17) but is missing in
those of Love waves (Fig. A1.18). These differences appear to be quite large particularly below a
depth of 7 km, which may be caused by the difference in data quality and the existence of
apparent radial anisotropy. A systematic analysis of the differences between the Rayleigh and
Love wave results will be a topic of a future study.

1.6 Discussion & Conclusions
We obtain 3-D tomographic images of S wave velocities with a grid size of 0.05º in the SC
plate boundary region using Eikonal tomography (Lin et al., 2009; section 1.4) and 1-D
linearized Vs inversion scheme (Herrmann, 2013; section 1.5). The study employed one year of
continuous seismic data recorded on more than 300 stations in SC. The preprocessing steps
discussed in Zigone et al. (2015) are first utilized to compute reliable daily ANC for every
station pair throughout the year (section 1.2), and then phase and group travel time dispersion
relations are extracted automatically from the surface waves reconstructed from all the one-year
stacked cross-correlation function (section 1.3). The Eikonal tomography allows for rapid
derivation of statistically robust and reliable isotropic Rayleigh and Love wave phase and group
velocities between 2.5s and 16s. The final 3-D Vs model is inferred by jointly inverting the
resulting isotropic phase and group dispersion curves through a series of 1-D linearized Vs
inversions at all the grid points.
The study incorporates the following methodological improvements. 1). An automatic
surface wave dispersion-picking algorithm based on frequency-time analysis is developed. To
maximize the number of measurements at shorter periods (e.g. ≤ 3s) while simultaneously
minimizing false detections, we perform the automatic picking procedure not only on the
symmetric signals but also on both the positive and negative time lags of the correlations.
Comparisons between results obtained from different components of the correlation functions
help filtering out erroneous measurements identified through inconsistency. 2). The
determination of phase travel time dispersion picking employs model-predicted travel time
dispersion curve from the CVM-S4.26 to avoid cycle skipping (N in equation 2). 3). The group
travel time dispersion picked using the method of several earlier studies (e.g. Barak et al., 2015;
39
Zigone et al., 2015) is found to be sensitive to noise and has larger uncertainty. We therefore
derive group traveltime dispersion using the obtained phase travel time dispersion following
equation 3, which improves the accuracy of the measurements (Figs. 1.6d&1.7a).
4). In addition to the quality control criteria developed in Lin et al. (2009), we introduce
station configuration error to identify regions that have unreliable gradient estimates due to poor
data coverage (section 1.4.1; Fig. A1.5). This further improves the quality of the final stacked
velocities. 5). We use both the phase and group dispersion curves to invert for Vs structures,
which yields better inversion results (section 1.5.1). 6). The resolvability of the iterative 1-D Vs
inversion is determined typically qualitatively using depth sensitivity kernels of surface wave
velocity (Fig. A1.3; e.g. Zigone et al., 2015). Here we use a neighborhood algorithm (Wathelet
2008; Zigone et al., 2019) to evaluate the non-uniqueness of the inversion quantitatively. The
resulting uncertainties show small values at the depth range of ~3-25km, at which the Vs profiles
inverted using different initial models are consistent. This suggests that although Rayleigh wave
group velocities at 3s have sufficient sensitivity and can constrain the Vs structures in the top 3
km (e.g. Barak et al., 2015; Zigone et al., 2015), the Vs values in the top 3 km from such 1-D
inversion are non-unique and likely be biased by the initial model.
The resulting tomographic model of Vs using Rayleigh wave data is consistent overall with
previous inferences on the large-scale velocity structure in SC (e.g. Tape et al., 2010; Allam et
al., 2012; Lee et al., 2014; Zigone et al., 2015; Barak et al., 2015; Fang et al., 2016; Berg et al.,
2018; Share et al., 2019). However, we find a large discrepancy between the surface wave
dispersions obtained in this study and those predicted by the SCEC community velocity models,
particularly for regions inside the basins and around fault zones (Figs. A1.11&A1.15). The
surface wave imaging results derived in the current study provide likely better results on these
structural features owing to the denser data and methodology improvements included in the
current study. In addition, we observed several important features not included in the initial
models including a reversed polarity of the velocity contrast across the segments of SAF that are
southeast and northwest to the SGP and a northeast dipping SAF southeast to the SGP.
Comparisons between our model results and the distribution of rock damage estimated in
Ben-Zion and Zaliapin (2019) from the background seismicity yields a good correlation between
the LVZ and large estimated damage volumes in the ECSZ (Fig. 1.15) and at depths of ~15-20
km beneath the SGP (DD’ of Fig. 1.17). The low velocity anomalies in the ECSZ seem to
40
coincide with the rupture zones of the M6.1 Joshua Tree, M7.3 Landers, and M6.3 Big Bear
earthquakes happened in 1992. The large damage volume beneath the SGP may be related
(Lyakhovsky and Ben-Zion, 2009) to a significant change in Moho geometry below the South-
Central Transverse Ranges (Zhu and Kanamori, 2000; Ozakin and Ben-Zion, 2015).
Compared to the initial models, our final Vs model better characterizes the fault zones in the
upper crust, which are illuminated by LVZ centered on the fault surface traces in the top 3-7 km.
Interestingly, the LVZ underneath the SJF is less significant compared to initial models,
particularly beneath the trifurcation area, suggesting the area may not be as localized as the south
SAF. In addition, we observe a low velocity strip beneath the LA basin dipping northward to ~15
km depth with an angle of ~30º (CC’ of Fig. 1.17). The estimated surface location, dipping
direction and angle, and depth range of the low velocity strip coincide with features of the Puente
Hills blind-thrust system imaged by Shaw et al. (2002). The estimated ~3-5% velocity reduction
of the low velocity strip compared to the surrounding structures (CC’ of Fig. 1.17) is consistent
with the fact that the Puente Hills blind-thrust system is capable of generating Mw6.0+
earthquake (Shaw et al., 2002).
In conclusion, the noise-based Eikonal tomography results significantly improve the Vs
structures in the top 3-20 km for the southern California plate boundary region, particularly near
fault surface traces. To obtain reliable Vs structures in the top 3 km, either surface wave velocity
dispersion at higher frequencies (e.g. > 1 Hz; Lin et al., 2013) or joint inversion with other types
of measurements such as Rayleigh wave H/V ratio (e.g. Berg et al., 2018) are required. The
observed differences between imaging results using data from Rayleigh and Love waves suggest
the existence of significant apparent radial anisotropy, which may be caused by transverse
isotropy (e.g. Moschetti et al., 2010a) or 3D structural effects (e.g. Levshin and Ratnikova, 1984).
The radial and 2-psi azimuthal anisotropy can provide additional information on crustal
properties (e.g. Moschetti et al., 2010b; Lin et al., 2011) in the study region and will be the
subject of a future study.

41
2. Internal structure of the San Jacinto fault zone at Jackass Flat
from data recorded by a dense linear array

2.0 Preamble
This chapter was published as: Qiu, H., Ben-Zion, Y., Ross, Z. E., Share, P. E., and Vernon,
F. L., 2017, Internal structure of the San Jacinto fault zone at Jackass Flat from data recorded by
a dense linear array: Geophysical Journal International, v. 209, no. 3, p. 1369-1388.
I was the main author, collected the earthquake waveforms recorded from the fault zone
linear array, performed the fault zone imaging analysis, and wrote the manuscript. Yehuda Ben-
Zion advised on the outline of the manuscript and took a lead on the interpretation of the imaging
results. Zachary E. Ross and Pieter-Ewald Share helped developing the analysis tools. Frank L.
Vernon provided the access to the data and general guidance. All authors contribute to paper
revision.
This study was supported by the National Science Foundation (grant EAR-162061) and the
U.S. Department of Energy (awards DE- SC0016520 and DE-SC0016527). We thank the Anza-
Borrego Desert State Park for permission to collect data on Jackass Flat. The seismic instruments
were provided by the Incorporated Research Institutions for Seismology (IRIS) through the
PASSCAL Instrument Center at New Mexico Tech. Data collected are available through the
IRIS Data Management Center. The facilities of the IRIS Consortium are supported by the
National Science Foundation under Cooperative Agreement EAR-1261681 and the DOE
National Nuclear Security Administration.

2.1 Introduction
Fault zone structures contain important information on various aspects of earthquake and
fault mechanics ranging from long-term evolutionary processes to brittle rock rheology and
dynamic stress fields operating during the occurrence of earthquakes (e.g. Ben-Zion and
Sammis, 2003; Xu et al., 2012; Rowe and Griffith, 2015). Some elements in the core structure of
fault zones can control (and reflect past) properties of earthquake ruptures. Specifically,
bimaterial interfaces separating different crustal blocks can influence significantly the location,
mode, speed and propagation direction of earthquake ruptures (e.g., Weertman, 1980; Ben-Zion
42
and Andrews, 1998; Brietzke and Ben-Zion, 2006; Ampuero and Ben-Zion, 2008; Brietzke et al.,
2009; Shlomai and Fineberg, 2016), and hence generated ground motion (e.g. Olsen et al., 2006;
Kurzon et al., 2014). A low velocity fault zone layer can produce oscillations of rupture velocity
and slip rate (e.g., Ben-Zion and Huang, 2002; Huang et al., 2014) and promote rupture
propagation (Weng et al., 2016). Asymmetric rock damage across the seismogenic fault may
indicate repeating occurrence of large earthquake ruptures with statistically preferred
propagation direction (e.g., Ben-Zion and Shi, 2005; Dor et al., 2006, 2008; Lewis et al., 2005).

Figure 2.1. (a) Location map for the San Jacinto fault zone (SJFZ) with 28995 earthquakes used
in this study (circles). The green triangle marks the location of JF and TR fault zone arrays that
are the focus of this paper. Other fault zone arrays are marked by grey triangles. The town of
Anza and surface traces of faults are shown by a green square and black lines, respectively.
Events analyzed for fault zone trapped waves (FZTW), fault zone head waves (FZHW) and P
delay times are contained within the black, purple and red boxes, respectively. Data from the
events marked by red circles are illustrated in later figures. (b) Depth section of hypocenters
projected along the cross-session AA’

in (a). The locations are poorly constrained for along-fault
distances larger than ∼50 km because of reduced network converge.
43
The surface expressions of fault zones often have high geometrical complexity with
hierarchical damage zones and multiple surface traces. This surface complexity and the diffuse
character of low magnitude seismicity (Fig. 2.1) make it difficult to identify the main
seismogenic fault. In the present paper we attempt to clarify the location at seismogenic depth
and internal structure of the Clark fault in the trifurcation area of the San Jacinto fault zone
(SJFZ), about 30 km southeast of Anza, California. The SJFZ is the most seismically active fault
zone in southern California (Hauksson et al., 2012) and accommodates a large portion of the
plate motion in the region (e.g., Fay and Humphreys, 2005; Lindsey and Fialko, 2013).
Paleoseismic and historic records indicate that the SJFZ is capable of large (M
W
> 7.0)
earthquakes (e.g., Petersen and Wesnousky, 1994; Rockwell and Seitz, 2008; Onderdonk et al.,
2013) that pose significant seismic hazard to large urban areas in southern California.
Allam and Ben-Zion (2012), Allam et al. (2014a) and Zigone et al. (2015) derived
earthquake- and noise-based tomographic models for the region around the SJFZ. These studies
imaged with nominal resolution of 1-2 km several key structural features including overall
velocity contrasts across the SJFZ and damage zones at different locations. Here we image
internal components of the Clark fault in the trifurcation area with spatial resolution ranging
from a few 100 m wide damage zone down to a local bimaterial interface. The imaging employs
earthquake data recorded by a dense linear array crossing the Clark fault at the Jackass flat (JF)
site, and a more sparse adjacent array (TR). These arrays are part of a PASSCAL deployment
starting in 2010 of 5 linear arrays crossing the SJFZ at different locations (triangles in Fig. 2.1)
and additional near-fault sensors (Vernon and Ben-Zion, 2010). The JF array crosses the surface
trace of the Clark fault in relatively flat topography and it has 9 broadband seismometers with
~20-30 m spacing between instruments in the center and a total aperture of ~400 m. The TR
array southwest of the JF array is situated about 100 m higher with 4 broadband seismometers
having ~ 0.5-1 km station spacing, and it does not cross the surface trace of any large fault (Fig.
2.2).
The tomographic results for the area indicate a large-scale lithology contrast across the
seismogenic fault at JF, with slower and faster P wave velocities on the SW and NE sides,
respectively (Fig. 2.3). The TR array is on the slower SW side of the fault, while the JF array
appears to cross the boundary between the two blocks although this is not clear given the 1-2 km
resolution of the velocity model. In the following sections we provide detailed results on the
44
velocity structure below the JF and TR arrays using local and teleseismic earthquake data and
several types of analyses. In section 2.2 we describe the data used in this study. Section 2.3
presents results based on delay times of P waves at stations across the arrays, and analyses of
fault zone head and trapped waves recorded by the JF array. The results indicate that the main
seismogenic fault is close to the SW end of the JF array, and the existence of a significant
shallow local low velocity zone NE of the fault. The latter produces a local reversal to the overall
large-scale velocity contrast consistent with persistent preferred propagation direction of
earthquake ruptures to the NW.

2.2 Data
The data analyzed in this work consist of waveforms of ~29,000 local and 9 teleseismic
earthquakes recorded during 2012-2014 with 200 Hz sampling rate by the JF and TR arrays
located in the Jackass Flat area (Fig. 2.2). The local earthquakes are taken from a catalog for the
SJFZ region utilizing the Anza network and nearby stations of the regional southern California
network and several local deployments (White et al., 2019). The corresponding waveforms are
extracted from the YN data set available at the IRIS Data Management Center (Vernon and Ben-
Zion, 2010). The black and purple rectangles in the top panel of Figure 2.1 show events used for
analysis of fault zone trapped and head waves, respectively, while the red rectangle indicates
events used for slowness analysis. Arrival times of P and S body waves generated by the local
events are found with the automated phase picking algorithms of Ross and Ben-Zion (2014) and
Ross et al. (2016). Initial detection of candidate fault zone head and trapped waves are made with
the automated algorithms of Ross and Ben-Zion (2014, 2015). The 9 analyzed teleseismic events
and corresponding waveforms are taken from the SCSN catalog and Southern California
Earthquake Data Center (SCEDC, 2013). These events have M > 6 and depth > 60 km. They are
selected for analysis because their recorded waveforms have sufficient Signal-to-Noise Ratio
(SNR) in a frequency band 0.01-0.5 Hz used in our analysis. The travel time of the teleseismic
phases at the JF and TR arrays are predicted using the global 1D velocity model iasp91 (Kennett
and Engdahl, 1991). The waveforms, arrival times and detected fault zone waves are analyzed in
detail below.
45

Figure 2.2. Location map for the JF and TR arrays at the SE end of the Clark fault. The JF array
has nine stations separated by ∼20–30 m in the center with a total aperture of ∼400 m, and is
located on relatively flat topography (top right inset). The TR array to the SW has four stations
separated by 0.5–1 km, with elevation ∼100 m higher than the JF array.
2.3 Analysis
Four types of studies involving different signals and scales are conducted to image the
internal structure of the SJFZ at the Jackass Flat area: teleseismic delay time analysis, local P
wave delay time analysis, and analyses associated with fault zone head waves (FZHW) and fault
zone trapped waves (FZTW). The analyses are described below starting with large-scale
structural features and progressing to inner fault zone components.
46

Figure 2.3. Average P-wave velocity over the depth range 1–7 km based on tomographic results
of Allam and Ben-Zion (2012). Red triangles indicate stations of the JF and TR arrays, as well as
a nearby CL array deployed in 1999.

2.3.1 Teleseismic arrival time delay analysis
2.3.1.1 Methodology
The selected teleseismic waves sample the crustal structure with near-vertical incidence
angles at lower frequencies than the local seismic waves. The relative arrival time delays of
teleseismic phases across an array can be used to infer variations in crustal velocity structure
below the stations (e.g. Cochran et al., 2009). If the array across a fault has large enough
aperture, the delays can further be used to quantify the average velocities of the different crustal
blocks separated by the fault (Ozakin et al., 2012). The JF and TR arrays have small apertures
and are located (Figure 2.3) within a broad low velocity damage zone that extends over the top 5
km or so of the crust (Allam and Ben-Zion, 2012; Allam et al., 2014a; Zigone et al., 2015). This
configuration is not suitable for imaging with teleseismic arrivals the deep structure across the
fault, but can be used to image variations within the damage zone and clarify the location of the
seismogenic fault.
5.7 5.8 5.9 6.0 6.1 6.2
Vp (km/s)
CL array
JF array
TR array
47

Figure 2.4. Teleseismic data recorded at the JF and TR arrays. (a) Red star and circles denote
locations of events analyzed in the main text and Supporting Information, respectively, while
green circles indicate additionally examined events. (b) Velocity seismograms at the JF and TR
arrays with the P and pP phases generated by the event denoted by the red star in (a), and
template produced by averaging the seismograms recorded by the JF array (bottom trace). The
red and blue dashed lines mark predicted arrival times based on the IASP91 global velocity
model. (c) Spectrogram of the template trace with amplitudes ≥35 per cent of the maximum
value (colors). The P phase energy is mostly below 0.2 Hz, while the pP phase has one-energy
amplitude below 0.2 Hz and another between 0.2 and 0.5 Hz.

We examine delay times of teleseismic P and pP phases across the stations, focusing on
frequency ranges with significant recorded spectral energy. For a low-frequency band, where the
structure is being sensed more broadly by the teleseismic waves, the delay times between
different stations are suitable for analyzing larger-scale velocity variations and can be applied to
an array with large station spacing. In the high-frequency case, the delay times are more sensitive
to the local velocity variations, and can only be derived robustly with small station separations.
As the frequency content of a teleseismic phase varies between events, we study the teleseismic
delay time separately for each individual teleseismic earthquake.
48

Figure 2.5. (a) Tapering window construction for teleseismic P and pP phases at 0.01–0.2 Hz.
Raw template trace on top is first bandpass filtered at the target frequency range (middle black
seismogram). Then, an envelope function is computed (red curve in the middle) and a tapering
window is built using the two nearest local minima (dashed lines) around the peak amplitudes
corresponding to the P and pP phases. The template traces before (purple) and after (black)
tapering are shown in the bottom. (b) Comparisons between velocity seismograms of teleseismic
P and pP phases before and after tapering. (c) Same as (b) for 0.01–0.5 Hz. (d) Same as (b) for
0.2–0.5 Hz.

Figures 2.4 and 2.5 illustrate the analysis of teleseismic P and pP phases generated by a deep
event with M = 6.5 (additional details are given in the next subsection). Since the station spacing
within the JF array is much smaller than the TR array, we generate a template trace by stacking
all waveforms observed at the JF array (Fig. 2.4b, bottom trace). Examining the dominant
frequency bands in the spectrogram of the template trace (Fig. 2.4c) indicates two frequency
bands (0.01-0.2Hz and 0.2-0.5Hz) with sufficient energy. Next we filter the waveforms at each
frequency band and produce shorter waveforms around each teleseismic phase. We select
tapering windows around the phases using two steps. First we filter the template trace to the
target frequency range (e.g., 0.01 – 0.2 Hz in Figure 2.5a) and compute a corresponding
envelope function (middle red curve in Fig. 2.5a). Then we locate the peak amplitude, t
peak
, of the
P or pP phase in the envelope function, find the nearest local minimums, t
min1
, t
min2
, in the
envelope function before and after the peak, and construct a tapering window such that the signal
between t
min1
and t
min2
remains unchanged (bottom red curve in Fig. 2.5a).
49
For the i-th station, the time shift
i
T
~
of a given phase is derived by cross-correlating the
truncated waveform around the target phase with that of the template trace. To obtain delay time
T
i
associated with the local velocity structure, we have to correct for expected delays due to
topography and the geometry of the stations and incoming plane wave. This is done by applying
the following two-step corrections on
i
T
~
.
1) For station i, the travel time of a target teleseismic phase is estimated with the TauP toolkit
(Crotwell et al., 1999) assuming the iasp91 global velocity model and that the station is at sea
level. The expected travel time delay at the i-th station due to the geometry of the stations and
incoming plane wave is approximately given by
Δ
!
t
i
=
!
t
i
−[
!
t
i
i=1
n
∑
]/n , (1)
where n is the number of stations considered in the analysis.
2) Assuming the teleseismic phase is traveling with vertical incident angle in the top portion of
the crust, the relative altitude at station i can be written as n h h h
n
i
i i i
/ ] [
1
∑
=
− = Δ with h
i
being the
station elevation. The delay in arrival time caused by topography is given by Δ
!
τ
i
=Δh
i
v
ref
,
where v
ref
is the velocity of the target phase in the top crust. The delay time at station i can then
be estimated as
T
i
≈
!
T
i
−Δ
!
t
i
−Δ
!
τ
i
. (2)
The spatial pattern of the delay time T
i
across the fault zone arrays can be derived for
different frequency ranges using the described steps. Figure 2.6 shows schematically expected
delay time patterns associated with two ideal conceptual fault models. By comparing the
obtained teleseismic delay times with the schematic patterns, we can estimate the basic type of
crustal structure beneath the arrays.

2.3.1.2 Results
Continuing with the example teleseismic event of Figure 2.4, the spectrogram for the P wave
exhibits a clear dominant frequency of ~0.15 Hz and the pP phase has two clear energy peaks,
one below 0.15 Hz and another between 0.2 Hz – 0.5 Hz (Figure 2.4c). We therefore compute
delay time patterns for this event using the overall frequency band 0.01 – 0.5 Hz and two sub-
50
ranges 0.01 – 0.2 Hz and 0.2 – 0.5 Hz. Figures 2.5b-d presents waveforms bandpass filtered in
these frequency ranges. Although the spectrogram in Fig. 2.4c does not show strong signal above
0.2 Hz for the P phase, the teleseismic P phase appears clearly in the waveforms after bandpass
filtering at 0.2 – 0.5 Hz (Fig. 2.5d). As mentioned, we taper the waveform around a given target
teleseismic phase before calculating cross-correlations for time delays. Comparisons of results
with different tapering choices indicate that the P waveforms are less sensitive than the pP
waveforms to details of the tapering for all three frequency bands. We therefore analyze only
delay time patterns associated with the P phases.

Figure 2.6. Schematic diagrams of slowness patterns for two different fault models. (a) A model
with a fault damage zone beneath the array. (b) Slowness pattern expected for the model in (a).
(c) A structure with a velocity contrast across the fault. (d) Slowness pattern expected for the
model in (c).

Figure 2.7a-c shows delay times for teleseismic P waves associated with the three used
frequency bands. Since the station spacing in TR array (~ 1 km) is much larger than the JF array
(~20 m), we only calculate the delay times for TR array at 0.01 – 0.5 Hz. The delay times are
only 3-4 sampling intervals, but the patterns show persistent variations of progressive velocity
reduction from the SW to the NE (TR04 to JF00). The overall delay time pattern across the
arrays remains when the reference velocity used for correcting the topography of the TR array
51
varies in the range 2-6 km/s (different lines in Fig. 2.7a). The same holds for corresponding
results associated with other examined teleseismic events. The detailed variations across the JF
array associated with the overall frequency range 0.01-0.5 Hz (Fig. 2.7a right) are consistent with
the schematic pattern in Fig. 2.6d, suggesting a possible vertical velocity interface separating two
crustal blocks between stations JFS3 and JFN1. The delay time patterns calculated for the JF
array at frequency ranges 0.01 – 0.2 Hz and 0.2 – 0.5 Hz are displayed (Figs. 2.7b, 2.7c) in the
same range outlined by the red dashed box in Fig. 2.7a. The results computed at the lower
frequency 0.01 – 0.2 Hz are similar to those at the overall range 0.01 – 0.5 Hz. However, for the
higher frequency range 0.2 – 0.5 Hz, the delay time reaches a maximum at stations JFN1 and
JFN2, still with generally more delay in the NW (JF00 – JFN4) than the SW side (JFS4 – JFS1).
This pattern is similar to the schematic curve in Fig. 2.6b, suggesting a low velocity zone (LVZ)
beneath the center of the array. These inferences are supported by results based on additional
teleseismic events, two of which shown in Supp. Figs. A2.1, A2.2, and slowness analysis using
~3500 local events enclosed in the red box of Fig. 2.1 (section 2.3.2).
To summarize, P waves travel faster beneath the TR than the JF array at the overall
frequency range of 0.01 – 0.5 Hz, in contrast with the expectation from the larger scale
tomography results of lower seismic velocities to the SW of the SJFZ (Fig. 2.3). The delay times
calculated within the JF array imply the existence of both a LVZ underneath the center of the
array and a potential vertical velocity interface possibly located between stations JFS3 – JFN1.
The reversed local shallow velocity contrast across the fault compared with the large-scale
seismogenic contrast is consistent with observations presented in later sections, observed
trapping structure to the NE of the fault at the nearby CL array (Lewis et al., 2005) and rock
damage asymmetry across the fault based on analysis of geomorphology features (Wechsler et al.,
2009). The local shallow reversal of the velocity contrast across the fault can be explained by
bimaterial ruptures on the SJFZ in the area with preferred propagation direction to the NW (Ben-
Zion and Shi, 2005; Xu et al., 2012).

52

Figure 2.7. (a) Delay times for teleseismic P phase at 0.01–0.5 Hz. The observed variations in P
arrivals are corrected with the traveltime delays due to the geometry of the station and incoming
plane wave and topography. The black dots denote delay times derived using 4 km s
-1

reference
velocity in topographic correction. The grey dashed lines indicate the lower and higher bound of
delay times computed at the TR array using reference velocities of 2 and 6 km s
-1
, respectively.
(b) Delay times for teleseismic P phase at 0.01–0.2 Hz. The delay times (blue dots) are shown in
the same range outlined by the red dashed box in (a) at the JF array. (c) Same as (b) for 0.2–0.5
Hz.

2.3.2 Local P wave Delay Time Analysis (DTA)
2.3.2.1 Methodology
Direct P wave arrival times from local earthquakes recorded at a dense array can provide
detailed information on the local fault zone structure. Similar to teleseismic arrivals, the
variations in direct P arrivals across the array reflect (with higher frequencies sensitive to smaller
scale features) the geometrical effect of wave propagation and heterogeneous velocity structure
below the array. Given the relatively large spacing between stations of the TR array and high
frequencies of the local data, we focus here on data recorded by the JF array only. After
removing the geometrical effect of the array geometry, the remaining delay time pattern across
the array indicates how P wave arrivals are affected by the local velocity structure. The use of
travel times from thousands of local earthquakes provides statistically robust estimates of local
velocity changes across the array.
Delay Time (s)
Northing from JF00 (km)
TR01
TR02
TR03
TR04
Topographic Correction with Vref
6.0 km/s
2.0 km/s
JFS4
JFS3
JFS2
JFS1
JF00
JFN1
JFN4
JFN3
JFN2
(a)
(b)
(c)
Delay Time (samples)
2
-3
2
-3
53

Figure 2.8. Early P-wave velocity seismograms recorded at the JF array for events marked as
red circles in Fig. 2.1. The red triangles denote the automatic P-wave picks. The green and blue
dashed lines indicate the time of the earliest arriving phase and the interpolated P-wave
traveltime across the array, respectively.

We first obtain automatic P picks (red triangles in Figure 2.8) using the algorithm of Ross
and Ben-Zion (2014). Analysis of significant number of local earthquakes in the Parkfield area
indicates that the median and standard deviation of the differences between manual and
automatic P picks are 0.004s and 0.023s, respectively (Ross and Ben-Zion, 2014). These values
imply that the errors in automatic P picks are negligible for the performed statistical analysis of P
wave delay time pattern. Figure 2.8 illustrates different arrival patterns at the JF array with data
of 4 example earthquakes. For events EQ#1 and EQ#2 southwest of the fault (Fig. 2.1), the P
3.0 3.2 3.4
-0.2
-0.1
0.0
0.1
0.2
Distance to station JF00 (km)
JFN4
JFN3
JFN2
JFN1
JF00
JFS1
JFS2
JFS3
JFS4
EQ#1
7.5 7.7 7.9
JFN4
JFN3
JFN2
JFN1
JF00
JFS1
JFS2
JFS3
JFS4
EQ#2
4.0 4.2 4.4
JFN4
JFN3
JFN2
JFN1
JF00
JFS1
JFS2
JFS3
JFS4
EQ#3
7.3 7.5 7.7
JFN4
JFN3
JFN2
JFN1
JF00
JFS1
JFS2
JFS3
JFS4
EQ#4
Time relative to origin time (s)
54
waves arrive earlier at the southwestern most station JFS4 compared to the northeastern most
station JFN4. The trend reverses for events EQ#3 and EQ#4 that are northeast of the fault (Fig.
2.1). These results indicate a directional dependent arrival pattern due to the source-receiver
geometries. However, in all 4 examples the direct P waves arrive later at the central stations JFS1,
JF00, JFS1 (Fig. 2.8). This consistent late arrival pattern, independent of the earthquake back
azimuth, is produced by a local velocity structure beneath the array.
For a linear array with aperture much smaller than the hypocenter distance, the variations in
travel time due to the different receiver locations can be estimated using a 1D velocity model. To
have a representative 1D velocity model, we average horizontally the 3D P wave velocity model
of Allam and Ben-Zion (2012) for the SJFZ region. The results remain essentially the same when
using a combined velocity model based on the results of Allam and Ben-Zion (2012), Allam et
al. (2014a) and Zigone et al. (2015). The 1D velocity model is used to compute predicted P
arrival time
!
t
ij
based on the algorithm of Kissling (1988) for each i-j earthquake-station pair.
Subtracting the predicted arrival at the central station JF00 gives predicted changes
0
~ ~ ~
ij ij ij
t t t − = Δ of the P arrivals from station JF00 due to the different propagation distances
within the array. This does not remove completely the geometrical effect, but averaging results
for all used earthquakes reduces further this source of variations.
The travel times of P waves depend on the hypocentral distance. To compare results
associated with different source-receiver distances we compute the average slowness α
ij
over
the entire ray path between sources i and receivers j in the following way. We first approximate
the propagation distance with the theoretical ray length
!
Δ
ij
0
from each earthquake to the
reference central station JF00 based on the 1D velocity model. Using the observed P wave travel
time t
ij
, the average slowness is estimated as
α
ij
= t
ij
−Δ
!
t
ij
( )
!
Δ
ij
0
, (3)
where Δ
!
t
ij
is the delay time correction for the geometrical effect for earthquake i at station j. The
average slowness α
ij
can still have systematic variations with the source location. For example,
the average slowness tends to be larger for shallow earthquakes than deeper ones. To reduce the
variability due to source location in the statistical analysis, we define relative slowness as
55

∑
=
=
n
j
ij ij
R
ij
n
1
/ α α α , (4)
where n is the number of stations used in the analysis. The dimensionless relative slowness α
ij
R

quantifies the relative relationship between the average slowness estimated at each station and
the array mean slowness value.

Figure 2.9. Illustration of the slowness calculation in local P-wave data. (a) Early P-wave
velocity seismograms for event EQ#3 (see Fig. 2.1 for location) with the red triangles indicating
the automatic P picks. Each seismogram is shifted in time by the delay due to geometrical effect
predicted by a 1-D velocity model. (b) Same as (a) with horizontal axis normalized by the
theoretical ray length between the hypocentre and reference station JF00. The red triangles
denote the slowness estimates with an average value across the array of 0.1746 km s
-1
. (c) Same
as (b) with horizontal axis normalized by the mean average slowness value across the array
calculated in (b). The relative slowness values are given by the red triangles.

2.3.2.2 Results
Figure 2.9 shows for the example event EQ#3 the observed P arrivals, along with the
estimated average and relative slowness for the JF array. To infer the velocity variations
underneath the array, we evaluate statistically the patterns of both average slowness and relative
4.0 4.1 4.2 4.3 4.4
Time (s)
-0.2
-0.1
0.0
0.1
0.2
Distance to station JF00 (km)
JFN4
JFN3
JFN2
JFN1
JF00
JFS1
JFS2
JFS3
JFS4
0.17 0.18
Average slowness (s/km)
JFN4
JFN3
JFN2
JFN1
JF00
JFS1
JFS2
JFS3
JFS4
Array avg:
0.1746 s/km
0.96 1.0 1.04
Relative slowness
JFN4
JFN3
JFN2
JFN1
JF00
JFS1
JFS2
JFS3
JFS4
(a) (b) (c)
56
slowness associated with many events. To filter erroneous P picks and obtain robust P delay
times, we require the P wave SNR to be larger than 30.0 and P wave apparent velocity along the
orientation of the array to be larger than 14.0 km/s. To ensure the 1D velocity model used to
correct for basic propagation effects is a good approximation for the medium the waves travel
through, we also reject P wave travel times not within ±1 s from the prediction based on the
velocity model.

Figure 2.10. (a) Histogram of average slowness for station JF00. (b) Same as (a) for relative
slowness. (c) Mean values of the average slowness (green) and relative slowness (black)
computed for 3493 events at the JF array. Error bars indicate a range of two standard deviations
about each respective mean value. The horizontal dashed line denotes the mean value of average
and relative slowness across the array.

The local P wave DTA analysis for the JF array starts with ~7800 earthquakes within 10 km
from the fault and within along-strike distance from the JF array of 25 km (red box in Fig. 2.1).
Events not satisfying the above quality criteria are dropped, leaving ~3500 events for refined
statistical analysis of average slowness and relative slowness of P wave at stations across the JF
array. Figure 2.10a and 2.10b show histograms of average slowness and relative slowness for
station JF00, respectively. As indicated by the range of horizontal axes, the width of average
57
slowness distribution is much larger than that of the relative slowness. Using similar histograms
for other stations of the JF array, the mean value and standard deviation of the mean for the
average and relative slowness distributions are calculated at each station. Figure 2.10c
summarizes the mean values of the average slowness (green) and relative slowness (black)
together with 2 times their standard deviation as error bars across the JF array. The two curves,
scaled to have similar overall amplitudes, collapse on the same pattern with difference primarily
in the error bars.
The patterns of the two slowness curves in Fig. 2.10c are similar to the one depicted in Fig.
2.6b, with a zone of larger slowness at the central part, possibly indicating a damaged zone at the
center of the array (Fig. 2.6a). No significant step of the type illustrated in Figs. 2.6c and 2.6d is
observed in the local P wave delay time analysis within the JF array. The shape of the slowness
patterns in Fig. 2.10c, obtained using local P arrivals averaged over thousands of events coming
from different azimuths, is consistent with the delayed time curve computed for teleseismic P
phase at 0.2 – 0.5 Hz (Fig. 2.7c). The dominant frequency of P wave generated by the local
earthquake is approximately 10 Hz. The delays of P waves generated by local events and those
produced by teleseismic events at frequencies above 0.2 Hz are likely dominated by the local
shallow low-velocity FZ structure underneath the JF array.

2.3.3 Fault Zone Head Wave (FZHW)
2.3.3.1 Methodology
Fault zone head waves are emergent seismic phases that refract along a bimaterial interface
in the fault zone structure, propagating along the interface with the faster wave speed and
radiated from there to the slower velocity medium (Ben-Zion, 1989, 1990). In a simple model
with two different quarter-spaces in contact, FZHW are the first arrivals at stations on the side
with slower seismic velocity having normal distance from the interface x less than a critical
distance x
c
given by
x
c
=r⋅tan cos
−1
α
s
/α
f
( ) ( )
, (5)
where r is the propagation distance along the bimaterial interface and
f
α ,
s
α are the P wave
velocities of the faster and slower media, respectively. The travel time of FZHW generated by
58
events on the interface and recorded by stations on the slow side can be expressed (Ben-Zion,
1989) as
t
H
= r /α
f
+ x α
s
−2
−α
f
−2
, (6)
while the corresponding travel time of the direct P wave is
s P
x r t α
2 2
+ = . (7)
The moveout between the FZHW and direct P arrival Δt =t
P
−t
H
for a station on the fault is
Δt = r /α
f
⋅η / 1−η ( )
, (8a)

where η is the fractional velocity contrast defined as α
f
−α
s
( )
/α
f
. In cases when η << 1, the
moveout is given approximately by
Δt ≈ r⋅Δα /α
2
, (8b)
where

α

and Δα

denote the average and differential P wave velocities, respectively (Ben-Zion
and Malin, 1991). In addition, the moveout Δt decreases monotonically with increasing station
distance x from the interface.
Observed FZHW and these relations were used in various studies to image properties of
bimaterial fault interfaces at sections of the San Andreas fault in California (e.g., McGuire and
Ben-Zion, 2005; Zhao et al., 2010; Share and Ben-Zion, 2016), the North Anatolian fault in
Turkey (Bulut et al., 2012; Najdahmadi et al., 2016) and other faults (e.g., Hough et al., 1994;
Zhao and Peng, 2008; Yang et al., 2015). FZHW also have Horizontal Particle Motion (HPM)
with a significant fault-normal component, since they are radiated from the fault, in contrast to
the particle motion of the direct P wave that points in the epicenter direction (e.g., Bulut et al.,
2012; Allam et al., 2014b).
Ross and Ben-Zion (2014) developed an automatic detector that searches the early P
waveform for a possible emergent phase before a sharper arrival, with a time separation between
a minimum value (0.065 s representing the width of a narrow P wave wiggle) and a maximum
value (corresponding to a 10% velocity contrast). The method was tested systematically on both
synthetic seismograms generated with the solution of Ben-Zion and Aki (1990) and observed
data in the Parkfield area where FZHW were detected manually before (Zhao et al., 2010). We
begin the analysis of FZHW by running the detection algorithm of Ross and Ben-Zion (2014) on
waveforms generated by >18,000 earthquakes that are within 10 km from the fault (purple box in
59
Fig. 2.1). To validate the FZHW detections, we visually inspect how the moveout between
tentative FZHW picks and direct P arrivals changes across stations of the JF array. Candidate
FZHW without a moveout trend across the JF array are rejected as false detections. This is
followed by HPM analysis to confirm the existence of significant rotation between the
polarization directions of the first and second phases in the early P waveforms.

Figure 2.11. (a) Location map of 18581 earthquakes (circles) used in the FZHW study (outlined
by the purple box in Fig. 2.1). Green and red circles denote candidate events producing FZHW
flagged through a combination of automatic detection and visual inspection. Data from an
FZHW candidate event (red circle) and a reference event (purple circle) are compared in later
figures. The region outlined by the red box is zoomed-in at the bottom left inset. (b) Depth
section of hypocenters projected along the cross-session AA’

in (a).

2.3.3.2 Results
Applying the automatic detection algorithm of Ross and Ben-Zion (2014) on earthquakes that
are within 10 km from the fault (black line AA’ in Fig. 2.11) leads to numerous candidate FZHW.
About 3000 events produce automatic detections in at least one JF station and about 1600 events
60
produce detection in more than one station. The tentative detected phases are inspected visually
and examined with HPM analysis to yield a high quality data set of FZHW.
Waveforms with automatic FZHW picks at multiple stations are aligned on the P picks made
at the onset of the earliest impulsive arrival. Figs. 2.12a and 2.12b show example waveforms at
stations of the JF array without and with FZHW. In Fig. 2.12a the first arrivals are standard
impulsive P body waves (red triangles), while in Fig. 2.12b waveforms at stations JFN4-JFS2
start with emergent tentative FZHW (green squares) followed by larger amplitude sharper P
body waves (red triangles). The moveout between the candidate FZHW and direct P arrivals in
Fig. 2.12b decreases from NE to SW within the JF array, suggesting that the velocity contrast
interface is on the NE side of the array.

Figure 2.12. (a) Velocity seismograms of early P waveform at the JF array generated by the
reference event denoted by the purple circle in Fig. 2.11. The waveforms are aligned on the
automatic P picks (red triangles). (b) Same as (a) for the candidate event shown as the red circle
in Fig. 2.11. The green squares indicate the automatic FZHW picks.

− 0.4 − 0.2 0.0 0.2 0.4 0.6 0.8 1.0
Time relative toP pick (s)
− 0.2
− 0.1
0.0
0.1
0.2
Distance to JF00 (km)
JFN4
JFN3
JFN2
JFN1
JF00
JFS1
JFS2
JFS3
JFS4
Ref example
(a)
− 0.4 − 0.2 0.0 0.2 0.4 0.6 0.8 1.0
Time relative to P pick (s)
− 0.2
− 0.1
0.0
0.1
0.2 JFN4
JFN3
JFN2
JFN1
JF00
JFS1
JFS2
JFS3
JFS4
HW example
(b)
61
To verify that the emergent early arrivals are not part of the direct P waves, we analyze the
HPM of the early P waveforms. This is illustrated in Figs. 2.13a and 2.13b using the example
waveforms of Fig. 2.12 at station JFS1. The top panels of Fig. 2.13a display HPM trajectories for
successive 0.1s windows of 0.8s in the early P waveform centered at the direct P arrival (red
triangle in Fig. 2.13a bottom). The azimuth and amplitude ratio between the successive sub-
windows of the HPM polarization are marked in the top panels. In this example with no FZHW,
the amplitude ratio becomes very large only at the first window that includes the direct P arrival,
and the HPM becomes polarized in the general horizontal direction of the body wave near the
array. In a waveform that includes FZHW, however, two windows with large amplitude ratios
are observed in the HPM trajectories (top panels in Fig. 2.13b). The first window with a large
amplitude change reflects the transition from noise to FZHW (green square in Fig. 2.13b bottom),
and the second window accompanied also by a significant rotation of the polarization direction
reflects the transition from FZHW to direct P phase.
Similar analyses lead to identification of 61 earthquakes, all located NW of the array (green
and red circles in Fig. 2.11), that produce clear FZHW at multiple JF stations. To estimate
properties of the velocity contrast interface that produce these FZHWs, we align the early P
waveforms at station JF00 generated by the 61 earthquakes on the direct P arrival and plot them
with increasing hypocenter distance (Fig. 2.14a). The moveout between the FZHW (green
squares) and direct P arrivals (red triangle) is constant for a wide range of hypocenter distances
(10 – 60 km). The constant moveout implies that the FZHWs observed at the JF array are
generated by a local bimaterial interface, which may be associated with edge of the fault damage
zone or basin, rather than a deep fault interface (Najdahmadi et al., 2016). Assuming that the
length of the interface at the edge of a damage zone is 3.5-5 km, as suggested by the trapped
waves analysis discussed next, leads to an estimated P wave velocity contrast of 10.7–14.6%
(equation 8a) using a 6 km/s average P wave velocity outside the damage zone. The short
propagation distance along the bimaterial interface explains (equation 5) the rapid change in the
moveout between FZHW and direct P arrival across stations of the JF array in Figure 2.12b, and
the lack of FZHW at the SW most stations as shown in Fig. 2.14b for station JFS4.
62

Figure 2.13. (a) Particle motion analysis on displacement seismograms at station JFS1 generated
by the reference event. The early P waveforms centred at the automatic P pick (red triangle) for
horizontal and vertical components are shown in green, black and purple in the bottom panel,
respectively. The signal is cut into six blocks with equal length by the black dashed lines. The
trajectory of horizontal particle motion in each block is depicted as the grey curve and compared
with the backazimuth of the event (red dashed line) in the corresponding top panel. The
trajectory is magnified in the first three windows with the amplification factor specified at the
top right of the panel. For each signal block, eigenvalue ratio and backazimuth are computed
through polarization analysis used in Bulut et al. (2012). (b) Same as (a) for the FZHW candidate
event. The green square denotes the automatic FZHW pick.
63

Figure 2.14. (a) Velocity seismograms recorded at station JF00 for FZHW candidate events
marked by green and red circles in Fig. 2.11. The early P waveforms are aligned with the
automatic direct P picks (red triangles) and plotted with increasing hypocentre distance. Green
squares indicate candidate FZHW picks. Black and green dashed lines indicate zero time lag and
best-fitting straight line crossing the green squares, respectively. A 0.1 s constant moveout
between the green and black dashed lines is observed. (b) Same as (a) for the southwest end
station JFS4. No FZHW are flagged.

64

Figure 2.15. (a) Location map for earthquakes that are analyzed for the presence of FZTW.
Events flagged as generating high-quality FZTW are displayed as red, green and yellow circles.
Seismograms of events denoted by green and red circles are used for moveout analysis (Fig.
2.16). Waveforms from events marked by red circles are inverted for structural parameters (Fig.
2.17 and Figs A2.7 and A2.8, Supporting Information). (b) Depth section of hypocenters
projected along the cross-session AA’

in (a).

2.3.4 Fault Zone Trapped Wave (FZTW)
2.3.4.1 Methodology
Fault zone trapped waves result from constructive interference of critically reflected phases
within a sufficiently uniform low-velocity fault zone layer (e.g. Ben-Zion and Aki, 1990; Igel et
al., 1997; Jahnke et al., 2002). The common type of FZTW is associated with the antiplane S
case and is analogous to surface Love waves of a horizontally layered structure. These waves
follow the direct S wave, have relatively low frequencies, are at least somewhat dispersive and
exist predominantly in the vertical and fault parallel components of ground motion (Ben-Zion
65
and Aki, 1990). Additional types of FZTW are associated with Rayleigh-type resonance or leaky
modes and appear between the P and S body waves (Ellsworth and Malin, 2011). In this paper
we discuss observations and analysis of the more common Love type trapped waves following
the direct S wave. We also observe possible candidates of Rayleigh or leaky type trapped waves,
but these are less consistent (the waveform character and recording stations vary for different
earthquakes) and are left for future work.
The amplitude, frequency content and duration of FZTW are very sensitive to combinations
of internal fault zone properties including the width, velocity reduction, attenuation coefficient
and propagation distance within the fault zone layer (Ben-Zion, 1998). Observations and
modeling of waveforms containing FZTW can provide high-resolution information on these
properties (e.g., Li et al., 1990; Peng et al., 2003; Lewis et al., 2005; Mizuno and Nishigami,
2006; Calderoni et al., 2012), although the modeling approach should account for the significant
trade-offs between model parameters (Ben-Zion, 1998; Jahnke et al., 2002; Lewis and Ben-Zion,
2010).
Following the pre-processing steps suggested by Ross and Ben-Zion (2015), waveforms are
first corrected for the instrument response, rotated to the fault-parallel component, and then
bandpass filtered at 2 – 20 Hz. The detection of potential FZTW starts by applying the automated
algorithm of Ross and Ben-Zion (2015) on the pre-processed data recorded at the JF array. The
data flagged by the automated algorithm are inspected visually to confirm the existence of
resonance mode wave packages after the direct S arrivals that are observed consistently at a
confined spatial range of the array.
Following Ben-Zion et al. (2003) and later studies, waveforms observed across the array that
include FZTW at some stations are inverted for properties of the trapping structure using a
genetic inversion algorithm with a forward kernel based on the 2D analytical solution of Ben-
Zion and Aki (1990) and Ben-Zion (1998). The analytical solution accounts for multiple plane-
parallel fault zone layers between two quarter spaces. However, in the present paper we find that
a model configuration consisting of a single fault zone layer in a half space (Supp. Fig. A2.5) is
sufficient to fit well the observed FZTW. This simple model has 6 key parameters that are
perturbed in the inversion. The fault zone parameters consist of attenuation coefficient (Q
FZ
), S-
wave velocity (β
FZ
), propagation distance (z
s
) within the fault zone, width (W) and location of
the fault zone layer (x) within the array. The host rock parameters are S-wave velocity (β
H
) and
66
attenuation coefficient (Q
H
). The latter is fixed to be 400, which is the average value for
Southern California (Hauksson and Shearer, 2006), to reduce somewhat the number of
parameters. The genetic inversion algorithm explores systematically the trade-offs between these
model parameters. The propagation time of S waves outside the fault zone layer is an additional
parameter used to check consistency of the results with the overall distances between the
generating earthquakes and JF array.

Figure 2.16. (a) Fault-parallel displacement seismograms generated by event TW4 (see location
in Fig. 2.15). The waveforms are pre-processed following the steps illustrated in Fig. A2.6 in the
Supporting Information. The red triangles denote the available automatic S picks. (b) Same as (a)
for event TW2. (c) Displacement seismograms after pre-processing recorded at station JFS1
generated by events marked by green and red circles in Fig. 2.15. The S waveforms are aligned
on the automatic S picks and plotted with increasing hypocentre distance indicated above each
trace. The blue and red dashed lines denote the S and FZTW arrivals, respectively.

To obtain a good estimation of the FZ parameters, 10,000 models with 50 generations and
200 models per generation are tested in each inversion. A successful inversion includes both
-1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5
Hypocenter distance
16.5 17.0 17.5 18.0 18.5 19.0 19.5
15.5 16.0 16.5 17.0 17.5 18.0 18.5 19.0
Time relative to S pick (s)
29.4 km
62.0 km
62.0 km
56.6 km
18.6 km
58.5 km
43.5 km
62.1 km
62.1 km
62.1 km
18.1 km
62.6 km
61.7 km
60.0 km
59.8 km
19.2 km
54.0 km
60.2 km
station: JFS1 (c)
South (y-) to North (y+)
JFN4
JFN3
JFN2
JFN1
JF00
JFS1
JFS2
JFS3
JFS4
TW4 (a)
Time (s)
South (y-) to North (y+)
JFN4
JFN3
JFN2
JFN1
JF00
JFS1
JFS2
JFS3
JFS4
TW2 (b)
67
good waveform fits, and agreement between parameters of the best fitting model and peaks of
the probability density distributions of parameters generated in the inversion process. As in Ben-
Zion et al. (2003) and later related studies, prior to the inversion the fault-parallel seismograms
are integrated to displacement and convolves with 1/t
1/2
to convert a point source response to that
of an equivalent line dislocation source (e.g. Vidale et al., 1985; Igel et al., 2002). Additional
details on the method can be found in Ben-Zion et al. (2003).

2.3.4.2 Results
As discussed in Ross and Ben-Zion (2015), applying the automatic detection algorithm to the
JF array data recorded during 2013 produced a total of 582 candidate FZTW with 90% of the
detections concentrated between JF00-JFS3. This suggests persistent generation of FZTW by a
narrow FZ layer beneath stations JF00-JFS3 of the JF array. To focus on high quality candidate
FZTW, we apply the automatic detection algorithm on 28995 earthquakes that are within 50 km
from the fault using a higher detection threshold than Ross and Ben-Zion (2015). This leads to
detection of trapped waves at stations JF00-JFS3 with an evaluation score > 6.0 generated by 205
events (colored circles in Fig. 2.15). The evaluation score is computed by summing up the Y-
statistics described in Ross and Ben-Zion (2015). The detected phases in this high quality data
set are further visually inspected and evaluated. We note that many detections including those
producing the clearest candidate FZTW (green and red circles) are generated by events located
considerably away from the fault. This implies a relatively shallow trapping structure that can be
excited by wave energy of regional events that enter the FZ layer from below (Ben-Zion et al.,
2003; Fohrmann et al., 2004).
Supp Fig. A2.6 illustrates the preprocessing steps performed on example event TW2. The
direct S wave picks are followed by resonance modes at stations JFS1-JFS3 that are especially
clear on the displacement seismograms. Similar preprocessing and visual inspection indicate
consistent large amplitude resonance phases confined to stations JFS1-JFS3 for all 205
earthquakes producing high quality FZTW. Figs. 2.16 (a) and (b) display fault parallel
displacement waveforms across the JF array generated by example events TW2 and TW4 (see
Fig. 2.15 for locations). The seismograms have clear FZTW at stations JFS1-JFS3 following the
direct S wave picks (red triangle). Waveforms generated by 18 earthquakes (green circles in Fig.
2.15) with clear and representative FZTW are used for further analysis. To estimate the overall
68
scale of the trapping structure, we align the early S waveforms generated by the 18 earthquakes
at station JFS1 on the direct S arrival, and arrange them with increasing hypocenter distance (Fig.
2.16c). The moveout between the direct S arrivals (blue dashed line) and FZTW (red dashed line)
is constant for a wide range of hypocenter distance (20 – 60 km). This constant moveout
supports the inference based on the regional distribution of events generating FZTW that the
trapping structure below the JF array is relatively shallow.

Figure 17. Inversion results modeling stacked displacement seismograms at the JF array
generated by candidate events TW1 and TW4. (a) Comparison between the observed (black) and
synthetic (green) seismograms. (b) Fitness values (green dots) calculated for different FZ
parameters tested in the inversion. The parameters associated with the best-fitting model (black
circles) are displayed in each panel and used to produce the synthetic waveforms shown in (a).
The black curves indicate probability density functions for the inverted model parameters.

69
Fig. 2.17 presents inversion results for the displacement seismograms at the JF array
generated by earthquakes TW1 and TW4. As the distance between the hypocenters of TW1 and
TW4 is ~2 km, the waveforms generated at each station by the two earthquakes are stacked to
increase the SNR. Performing inversions for waveforms generated by each event separately lead
to very similar results to those discussed below. Fig. 2.17a compares synthetic waveforms
associated with the best fitting model (green lines) and observed seismograms (black lines). Fig.
2.17b presents fitness values calculated by the genetic inversion algorithm for different FZ
parameters. The fitness for a given set of parameters is defined as (1+C)/2, where C is the cross-
correlation coefficient between the set of synthetic and observed waveforms at all stations. The
thin curves in Fig. 2.17b give probability density functions for the various model parameters,
calculated by summing the fitness values of the final 2000 inversion iterations and normalizing
the results to have unit sums. The parameter values associated with the best fitting model (black
dots) are very close to the peaks of the probability density functions. The close waveform fits in
Fig. 2.17a and the results in Fig. 2.17b imply convergent and robust inversion results. The best
fitting parameters indicate that the trapping structure below the JF array has width of ~200 m,
length scale of ~5 km, Q
S
of ~20 and S wave velocity reduction of ~35% compared to the host
rock. The ~5 km propagation distance along the waveguide is associated with a combination of
along-strike and vertical components. Assuming for simplicity that the average along-strike and
vertical components are the same (Ben-Zion et al., 2003), the depth of the trapping FZ layer is
~3.5 km. Supp Figs. A2.7, A2.8 show inversion results of FZTW generated by earthquakes TW2
and TW3. The best fitting parameters and peaks of the probability density functions are similar
to those shown in Fig. 2.17.

2.4 Discussion
We use delay times of P waves from local and teleseismic earthquakes, fault zone head
waves and fault zone trapped waves to image the internal structure of the Clark fault at the
Jackass Flat site in the trifurcation area of the SJFZ. The identification of the various phases is
done first with automated detection algorithms (Ross and Ben-Zion, 2014, 2015) enabling
systematic and objective analyses of large data sets (~29,000 events recorded at 9 stations). This
is followed by visual inspection of detected phases augmented by some additional analysis (e.g.,
examining the horizontal particle motion of candidate head waves) to substantiate the validity of
70
the automatic detections. These procedures lead to identification of >60 near-fault events
producing early P arrivals at multiple stations of the JF array with characteristics of FZHW (Fig.
2.11), >200 broadly distributed events producing high quality FZTW at stations JF00-JFS3 (Fig.
2.15), and P body wave arrivals from ~3500 earthquakes in the region around the SJFZ (Fig.
2.1).

Figure 2.18. Google Earth map view of a 2.5 × 2 km region centered at the JF array. The red
triangles denote stations of the JF and CL arrays with additional two stations of the TR array.
The black dashed thick line parallel to the fault strike indicates the estimated location of
seismogenic fault. The blocks to the NE and SW of the fault are marked as fast and slow based
on previous tomography results (Fig. 2.3). A fault zone model along cross-session AA’
integrating the analysis results is shown in Fig. 2.19.

Teleseismic P waves in the frequency range of 0.01–0.5 Hz observed at the JF and TR arrays
indicate faster seismic velocities beneath the TR array compared to the JF array, in contrast with
the expectation from the larger-scale tomography results (Fig. 2.3). The reversal of seismic
velocities near the fault points to a prominent local LVZ underneath the JF array. This is
supported by both the teleseismic P wave delay time derived from the frequency range 0.2–0.5
Hz and statistical analysis of the P wave delay times associated with numerous local events.
Some of the velocity reduction below the JF array is related to the fact that it is located in a local
basin with sediments. However, the basin is part of the fault zone structure (e.g., Sibson, 1986;
Ben-Zion et al., 2003) and it includes an inner low velocity zone that forms a seismic trapping
structure. Assuming that the average P wave velocity of the host rock is ~6 km/s (Fig. 2.3) and
71
that the LVZ beneath the JF array has ~20% velocity reduction in the top 3.5 km (based on the
velocity contrast indicated by the FZHW and some additional reduction in the more intense inner
damage zone producing FZTW), a vertical incident P wave is delayed beneath the JF array by
~0.15 s compared to the TR array. This is sufficient to reverse the travel time variation that may
be produced by the larger-scale deeper structures shown in Fig. 2.3. The delay time results also
suggest that the main seismogenic fault is between stations JF3 and JF4. This is consistent with a
prominent large-scale geomorphologic feature in this section of the fault (Fig. 2.18).
The observation of a constant moveout between FZHW and direct P phases indicates that the
FZHW are generated by a local velocity contrast, inferred to be at the edge of the LVZ on the NE
side of the JF array. This is confirmed by the fact that the observed 0.1 s constant FZHW–P
moveout implies a 10–15% contrast of P wave velocities across the bimaterial interface, which is
larger than any deep velocity contrast interface found in previous studies at various faults (e.g.
Zhao et al., 2010; Yang et al., 2015; Najdahmadi et al., 2016; Share and Ben-Zion, 2016). The
LVZ beneath the JF array contains a narrower damage zone concentrated below stations JFS1-
JFS3 generating clear FZTW. Inversions of sets of waveforms observed across the JF array with
high quality FZTW at stations JFS1-JFS3 indicate (Figs. 2.17, A2.7, A2.8) that the trapping
structure is ~200 m wide and has V
S
reduction of ~30-40%, Q
s
of ~20 and depth of ~3.5 km. The
relatively shallow depth of the trapping structure is supported by (Fig. 2.15) the wide distribution
of events generating FZTW (e.g. Fohrmann et al., 2014) and the constant moveout (Fig. 2.16)
between the direct S and trapped waves (e.g. Ben-Zion et al., 2003; Peng et al., 2003). These
results are similar to properties estimated from analysis of FZ-reflected phases at the JF array
(Yang et al., 2014), as well as analysis of FZTW and FZ-reflected phases at the CL array
deployed across the Clark fault about 1 km to the SE from the JF array in 1999 (Li and Vernon,
2001; Lewis et al., 2005; Yang and Zhu, 2010). We note that Li and Vernon (2001) suggested
that the trapping structure below the CL array extends to a depth of ~18 km, but this has been
disputed by later analyses (Lewis et al., 2005; Yang and Zhu, 2010).
Fig. 2.19 provides a schematic summary of the different imaging results done in this work,
with a focus on the structure below the dense JF array. The seismogenic fault is inferred to be
located between stations JFS3 and JFS4 and is associated with a strongly asymmetric LVZ that
exists primary on the NE side of the fault. A similar asymmetric LVZ was observed at the nearby
CL array (Lewis et al., 2005; Yang and Zhu, 2010), and in larger-scale analysis of
72
geomorphologic features in the area (Wechsler et al., 2009). The offset of rock damage across
the fault to the side with faster velocities at seismogenic depth (Fig. 2.3) may be produced by
preferred propagation direction of earthquake ruptures in the area to the NW (e.g., Ben-Zion and
Shi, 2005; Xu et al., 2012). This is consistent with rupture directivities of small events (Kurzon
et al., 2014) and along-strike asymmetry of stacked aftershock sequences in the area (Zaliapin
and Ben-Zion, 2011).
The results of this paper are based on data recorded across the Clark fault at one location.
The SJFZ has heterogeneous structures with significant along-strike variations of the surface
geology (e.g., Sharp, 1967), seismicity (e.g. Hauksson et al., 2012), interseismic strain rates (e.g.
Wdowinski, 2009; Lindsey and Fialko, 2013) and large-scale seismic velocities (e.g. Allam et al.,
2014a; Zigone et al., 2015). Yang et al. (2014) and Li et al. (2015) documented along-strike
variations of local damage zones and seismic anisotropy using data recorded by the various
linear arrays across the SJFZ (triangles in Fig. 2.1). Several studies are currently underway on
detailed analyses of the type done in this work using data recorded by the other dense linear
arrays.

Figure 2.19. A conceptual model for fault zone structure at the JF site along the cross-session
AA’

in Fig. 2.18. The lower velocities beneath stations JFS3–JFN2 (warmer colors) are indicated
by delay times of teleseismic and local P waves. The red thick line denotes a local bimaterial
interface that generates FZHW (Section 3.3). A fault zone waveguide is depicted by a vertical
layer with NE boundary located at JFN1. The marked S wave Q, velocity reduction, width and
depth of the waveguide are constrained by waveform inversions (Section 3.4). The inferred
seismogenic fault crosses the JF array between stations JFS3 and JFS4.
JFS4 JFS3
JFS1
SW NE
50m
JFS2 JF00
JFN1
JFN2
JFN3 JFN4
Fault surface trace
Velocity contrast interface
Averaged waveguide
Seismogenic fault
Vs reduction 30-40%
Q = ~20
~200 meters
~3-5 km deep
Slow block Fast block
A A’
Slowness
Figure 19
73
3. Internal structure of the San Jacinto fault zone at Dry Wash from
seismic data of a linear array

3.0 Preamble
This chapter is currently in preparation as: Qiu, H., Ben-Zion, Y., and Vernon, F. L., 2019c,
Internal Structure of the San Jacinto Fault Zone at Dry Wash from Data Recorded by a Dense
Linear Array.
I was the main author, collected the earthquake waveforms recorded from the fault zone
linear array, performed the fault zone imaging analysis, and wrote the manuscript. Yehuda Ben-
Zion advised on the outline of the manuscript and took a lead on the interpretation of the imaging
results. Frank L. Vernon provided the access to the data and general guidance. Both co-authors
helped improving the manuscript.
This study was supported by the National Science Foundation (grant EAR-162061) and the
U.S. Department of Energy (awards DE- SC0016520 and DE-SC0016527). We thank the Anza-
Borrego Desert State Park for permission to collect data on Dry Wash. The seismic instruments
were provided by the Incorporated Research Institutions for Seismology (IRIS) through the
PASSCAL Instrument Center at New Mexico Tech. Data collected are available through the
IRIS Data Management Center. The facilities of the IRIS Consortium are supported by the
National Science Foundation under Cooperative Agreement EAR-1261681 and the DOE
National Nuclear Security Administration.

3.1 Introduction
Here we image internal structure of the Clark fault in the trifurcation area of the San Jacinto
fault zone (SJFZ) using the dense linear array crossing the fault surface trace at the Dry Wash
(DW), ~ 6 km northwest to the Jackass Flat (JF) discussed in Chapter 2 (Fig. 3.1c). As shown in
Figs. 3.1a&b, DW array sits on a relatively flat topography within a narrow valley and contains
12 short period seismometers with ~40 m spacing and a total aperture of ~400 m. The closest
stations in the regional seismic networks are DWRPT in the southwest and ALCY in the north,
which are more than 200 m higher in elevation than the DW array.
74

Figure 3.1. Location map for the DW array. The DW array consists of 12 stations separated by
~20 m in the center, has an aperture of ~400 m, crosses the fault from SW to NW and is located
on relatively flat topography. Complementary to the DW array, the reference station DWRPT
extends farther to the southwest. The elevation is ~200 m higher than the DW array.

Previous large-scale tomographic images for the area (e.g. Allam and Ben-Zion, 2012; Allam
et al., 2014a) indicate a large-scale lithology contrast across the seismogenic fault at DW, with
slower and faster P wave velocities on the SW and NE sides, respectively (Fig. 3.2). Different
from the JF site (Qiu et al., 2017), there are several surface traces of the Clark fault indicating a
much complex local fault structure underneath the DW array. In addition, orthogonal structures
are seen between the Clark and Coyote Creek faults both from fault surface traces (Fig. 3.2) and
seismicity patterns (e.g. Hauksson et al., 2012; Ross et al., 2017a; Cheng et al., 2018). Together
with the fact that the distances between the Clark fault and the other two branches, Buck Ridge
and Coyote Creek, are less than 5 km, it suggests the fault zone at DW site is a very complicated
environment, particularly for the analysis of fault zone phases.
In the following sections, we provide detailed results on the velocity structure below the DW
array using local and teleseismic earthquake data. In section 3.2 we describe the data used in this
study. Section 3.3 presents results based on delay times of P waves at stations across the arrays,
and analyses of fault zone head and trapped waves recorded by the DW array. The results
indicate that the main seismogenic fault is at the center of the DW array (between DW11 and
75
DW12), and the lack of clear S-type trapped waves suggest no effective low velocity guided
zone is found, thus indicating a more localized fault zone structure than that at JF.

Figure 3.2. Average P-wave velocity over the depth range 1–7 km based on tomographic results
of Allam & Ben-Zion (2012). Green and yellow triangles indicate stations of the DW array and
the closest nearby regional stations (DWRPT and ALCY), respectively. The blue star, beach ball,
and black arrow represent the epicenter, focal mechanism, and observed rupture directivity of the
most recent Mw > 5 earthquake, the 2016 Borrego Springs event, occurred in the trifurcation
area of the SJFZ.

76

Figure 3.3. Location map for the San Jacinto fault zone (SJFZ) with ~19,000 earthquakes used
in this study (circles). The green triangle marks the location of DW fault zone array that is the
focus of this paper. Other fault zone arrays are marked by grey triangles. The town of Anza and
surface traces of faults are shown by a green square and black lines, respectively. Events
analyzed for fault zone trapped waves (FZTW), fault zone head waves (FZHW) and P delay
times are contained within the black, purple and red boxes, respectively. Data from the events
marked by red and yellow circles are illustrated in later figures. (b) Depth section of hypocenters
projected along the cross-session AA’ in (a)

3.2 Data
The data analyzed in this work consist of waveforms of ~19,000 local and 58 teleseismic
earthquakes recorded during 2015-2016 with 200 Hz sampling rate by the DW array located in
the Dry Wash area (Fig. 3.3). Although the DW array was deployed in 2012, we only analyze the
data recorded starting from 2015, as stations DW11 and DW12 were installed at the end of the
year 2014. Local earthquakes are taken from the catalog of Hauksson et al. (2012). The
77
corresponding waveforms are extracted from the YN data set available at the IRIS Data
Management Center (Vernon and Ben-Zion, 2010). The black and purple rectangles in the top
panel of Figure 3.3 show events used for analyses of fault zone trapped and head waves,
respectively, while the red rectangle indicates events used for slowness analysis. Arrival times of
P and S body waves generated by the local events are found using the automated phase picking
algorithms of Ross and Ben-Zion (2014) and Ross et al. (2016). Initial detection of candidate
fault zone head and trapped waves are made with the automated algorithms of Ross and Ben-
Zion (2014, 2015). The 58 analyzed teleseismic events and corresponding waveforms are taken
from the SCSN catalog and Southern California Earthquake Data Center (SCEDC, 2013). These
events have M > 5.5 and depth > 60 km. They are selected for analysis because their recorded
waveforms have sufficient Signal-to-Noise Ratio (SNR) in a frequency band 0.2-2 Hz used in
our analysis. The travel time of the teleseismic phases at the DW array are predicted using the
global 1D velocity model iasp91 (Kennett and Engdahl, 1991). The waveforms, arrival times and
detected fault zone waves are analyzed in detail below.

3.3 Methods and Results
3.3.1 Teleseismic earthquake delay time
Figure 3.4 shows the location, waveform, spectrogram, and delay times of one example
teleseismic event. Here, we apply the same teleseismic analysis described in section 2.3.1.1,
cross correlating the tapered P arrival waveform (top panel of Fig. 3.5) of each individual station
and the array mean reference trace. The resulting delay time pattern for the example event of the
frequency band 0.2 – 2 Hz is illustrated in Fig. 3.4e, and it is consistent with the schematic
pattern in Fig. 2.6d suggesting a possible vertical velocity interface separating two crustal blocks
between stations DW11 and DW12 (black dashed lines in Fig. 3.4e). Unlike the results of JF
array discussed in Chapter 2, the inferred velocity contrast is consistent with the larger scale
tomography results of lower seismic velocities to the SW of the SJFZ (Fig. 3.2) and does not
dependent on frequency bands that are used (e.g. 0.2 – 0.5 Hz, 0.5 – 1 Hz, and 1 – 2 Hz; Fig.
A3.1).
78

Figure 3.4. Teleseismic data recorded at the DW array. (a) Velocity seismograms at the DW
array with the P phases generated by the event denoted by the red star in (d). The black dashed
lines mark predicted arrival times based on the IASP91 global velocity model. (b) Template
produced by averaging the seismograms recorded by the DW array. The blue curve depicts the
corresponding envelope function. (c) Spectrogram of the template trace with amplitudes ≥35 per
cent of the maximum value (colors). The P phase energy is mostly between 0.2 Hz and 2 Hz. (d)
Red star denotes the location of events analyzed in the main text, while green circles indicate
additionally examined events. (e) Delay times and maximum correlation coefficients obtained by
applying the analysis described in section 2.3.1.1. The topography is corrected by assuming a
velocity of 4 km/s.

In addition to the delay time pattern, we also observe clear waveform changes in the PcP
phase crossing the boundary between DW11 and DW12 (Fig. 3.4a). To further characterize the
observed waveform changes between stations, we calculate the similarity and delay time
matrices by cross correlating the tapered P waveforms recorded by every two stations (Ozakin et
al., 2012) in DW array. In the later discussions, we use the frequency of 0.2 – 0.5 Hz as it yields
the clearest results. Figure 3.5 shows the resulting patterns for the example teleseismic event
indicated in Fig. 3.4d. In general, we see two sub-squares of high values along the diagonal of
the similarity matrix that are well separated by the black dashed lines (the boundary between
DW11 and DW12), suggesting different structural properties (e.g. seismic velocity) across the
seismogenic fault likely cause the observed waveform changes. As to the delay time matrix, we
see small delay times between stations on the same crustal block (i.e. DW01-DW11 for SW;
79
DW12-DW10 for NE), whereas the time difference between stations from different blocks is ~
0.04s with SW being slower than the NW.

Figure 3.5. Delay time analysis following Ozakin et al. (2012). P waveform recorded by each
station is first tapered following section 2.3.1.1 (top panel), and then cross correlates with each
other for every station pair. The obtained maximum cross correlation coefficient and delay time
matrices are shown as left and right bottom panels in color.

By averaging results obtained from all 58 teleseismic events that yield sufficient SNR (green
circles in Fig. 3.4d), Figure 3.6 shows the mean values of delay times achieved using the
methods described in section 2.3.1.1 (top panels) and illustrated in Fig. 3.5 (bottom panels). To
summarize, P waves travel faster beneath the NE than the SW of DW array at frequency range of
0.2–2 Hz, in agreement with the larger scale tomography results of lower seismic velocities to
the SW of the SJFZ (Fig. 3.2). The delay times calculated within the DW array imply the
existence of a potential vertical velocity interface possibly located between stations DW11 and
DW12 that is also responsible for the observed waveform changes across the contrast plane. The
lack of the existence of a LVZ underneath the center of the array suggests a more localized fault
structure compared to that of JF (Qiu et al., 2017).
80

Figure 3.6. The top panel shows the average time delay pattern of Fig. 3.4e over 58 teleseismic
events. Error bars indicate a range of two standard deviations about each respective mean value.
The red dashed curve illustrates the arctan function that fits the observed delay times the best. δT
and X
c
are the best fitting maximum time difference and location of the contrast plane relative to
the array center. The bottom two panels are the average results of Fig. 3.5.

3.3.2 Local earthquake delay time
Following the processing procedures described in section 2.3.2.1, we obtain the relative
slowness patterns of 1845 local seismic events that pass the selection criterion and the results are
shown in Figure 3.7. The delay times of local seismic events show a faster velocity beneath the
SW (DW01 – DW11) than NE (DW12 – DW10) of the DW array, consistent with the
observations using teleseismic earthquakes discussed in section 3.3.1. The delays of P waves
generated by local events and those produced by teleseismic events at frequencies above 0.2 Hz
are likely dominated by a velocity contrast across the fault plane underneath the DW array.

3.3.3 Fault zone head waves (FZHW)
Same as described in section 2.3.3.1, we flag FZHW candidates by using the automatic
detection algorithm of Ross and Ben-Zion (2014) first. The tentative detected phases are then
inspected visually for verification. Figure 3.8 shows the resulting FZHW candidates after visual
confirmation. P waveforms of two example FZHW candidates are shown in Figure 3.9. For all
81
the candidates, we find FZHW signals at all stations in DW array suggesting the velocity contrast
interface is SW to the entire DW array. Combined with the fact that no obvious evidence of a
localized LVZ beneath the DW array from delay time analysis, the graduate change in P-FZHW
moveout across the array (Fig. 3.9) likely indicates a deep fault contrast interface (Najdahmadi et
al., 2016).

Figure 3.7. (a) Mean values of the relative slowness computed for the 865 events SW (red) and
981 events NE (green) of the fault at the DW array. The black curve denotes the relative
slowness for all 1845 events. Error bars indicate a range of two standard deviations about each
respective mean value. (b) Same as (a) for 617 events SE (red) and 1229 events NW (green) of
the fault. (c) Same as (a) for 685, 544, 180 and 437 events in quadrants I (green), II (purple), III
(blue) and IV (red). Despite the variations between relative slowness calculated for events in
different subsets, the shapes of the mean values remain similar. (d) Example local seismic P
waveform and the corresponding automatic P picks (red triangles). The location of the example
event is shown as the yellow circle in Fig. 3.3. (e) Location map and depth section of 1845 local
earthquakes within the red box shown in Fig. 3.3(a), used in the local P-wave DTA to estimate
relative slowness.
82

Figure 3.8. Same as Fig. 3.3 for FZHW results. The red circles denote the FZHW candidate
events used for demonstration in Fig. 3.9, whereas the blue circles depict all FZHW candidates
after visual verification.

To validate the observed candidates and pinpoint the surface location of the velocity contrast
interface, further analyses, such as horizontal particle motion and comparison with reference
regional seismic stations (e.g. DWRPT, ALCY) are needed.

3.3.4 Fault zone trapped waves

83

Figure 3.9. P waveforms of two example FZHW candidate events shown as the red circles in
Fig. 3.8. The red dashed lines and green solid line indicate the predicted P arrival and the
impulsive P arrival times. Clear emergent FZHW are seen between the red and green lines, with
the P-FZHW moveout decreasing from DW01 to DW10, suggesting a velocity contrast interface
SW to the DW01.

Figure 3.10. Example of the clearest S-type FZTW found in the DW array data. The left panel
shows the location of the example event (red star). The P (vertical) and S (fault parallel)
waveforms are shown on the upper and bottom right panels, respectively. A large amplitude
wave package is seen between stations DW02 – DW11 ~0.5 s after the array mean S pick (blue
dashed lines).

84
Similar to the procedures described in section 2.3.4.1, we use the algorithm of Ross and Ben-
Zion (2015) to flag S-type FZTW candidates. Histograms of the resulting detections are shown
in Fig. A3.2. Most of the detections are made on just one DW station per event, and further
visual verifications indicate no clear S-type FZTW (the best example is shown in Fig. 3.10).
However, plenty of P-type FZTW are observed between stations DW04 – DW06 (Fig. 3.11)
suggesting a ~200 m wide trapping structure for P wave underneath the center of DW array.
Analyses of S wave delay times and modeling of P-type FZTW are required to further
understand the above observations, which will be the subject of future studies.

Figure 3.11. Example of P-type FZTW and waveform changes. The event location of this
example is denoted as the red circle in Fig. 3.3.

In addition to FZTW, we also observe clear waveform changes from the local seismic events
crossing the boundary between DW11 and DW12 (right panel in Fig. 3.11) and fault zone
reflected waves (Fig. A3.3).

3.4 Discussion
We use delay times of P waves from local and teleseismic earthquakes, FZHW and FZTW
to image the internal structure of the Clark fault at the DW site in the trifurcation area of the
0.5 – 10 Hz 10 – 20 Hz
~ 200m
Time relative to predicted P arrival
85
SJFZ. The identification of the various phases is done first with automated detection algorithms
(Ross and Ben-Zion, 2014, 2015) enabling systematic and objective analyses of large data sets
(∼19,000 events recorded at 12 stations). This is followed by visual inspection of detected phases
augmented by some additional analysis to substantiate the validity of the automatic detections.
These procedures lead to identification of 15 near-fault events producing early P arrivals at
multiple stations of the DW array with characteristics of FZHW (Fig. 3.8), no obvious S-type
FZTW are observed but plenty of P-type FZTW are seen at stations DW04 – DW06 (Fig. 3.11),
and P body wave arrivals from ∼1800 earthquakes in the region around the SJFZ (Fig. 3.3).
Teleseismic P waves in the frequency range of 0.2–2 Hz observed at the DW array indicate
faster seismic velocities beneath the NE compared to SW of the DW array, consistent the larger
scale tomography results (Fig. 3.2). The lack of prominent local LVZ underneath the DW array
suggests a localized fault structure and is consistent with the absence of clear S-type FZTW.
Assuming that the average P-wave velocity of the host rock is ∼6 km/s (Fig. 3.2) and that the
velocity contrast across the seismogenic fault is ∼5% in the top 10 km, a vertical incident P wave
is delayed with a maximum time of ~0.08 s beneath the DW array. This is consistent with the
best fitting maximum time differences based on the delay times averaged over 58 teleseismic
events (Fig. 3.6). The delay time results also suggest that the main seismogenic fault is between
stations DW11 and DW12. This is consistent with a prominent large-scale geomorphologic
feature in this section of the fault.
The observation of a graduate changing moveout between FZHW and direct P phases
indicates that the FZHW are generated by a deep fault interface, inferred to be SW of the DW
array (Fig. 3.9). Interestingly, we do not see clear S-type FZTW in data recorded by DW array,
but find several good quality P-type FZTW.
The results of this paper are based on data recorded across the Clark fault at one location.
The SJFZ has heterogeneous structures with significant along-strike variations of the surface
geology (e.g. Sharp, 1967), seismicity (e.g. Hauksson et al., 2012), interseismic strain rates (e.g.
Wdowinski, 2009; Lindsey and Fialko, 2013) and large-scale seismic velocities (e.g. Allam et
al., 2014a; Zigone et al., 2015). Similar to Yang et al. (2014) and Li et al. (2015), we also
documented significant differences between the fault structures at DW and JF.
86
4. Temporal changes of seismic velocities in the San Jacinto Fault
zone associated with the 2016 Mw5.2 Borrego Springs earthquake

4.0 Preamble
This chapter is currently in review as: Qiu, H., Hillers, G., and Ben-Zion, Y., 2019b,
Temporal changes of seismic velocities in the San Jacinto Fault zone associated with the 2016
Mw5.2 Borrego Springs earthquake: Geophysical Journal International.
I was the main author, computed two years of daily ambient noise cross-correlations of two
fault zone linear arrays, performed the monitoring analysis, and wrote the manuscript. Gregor
Hillers advised on the outline of the manuscript and interpretation of the monitoring results.
Yehuda Ben-Zion provided the general guidance. Both co-authors helped improving the
manuscript.
We thank the Anza-Borrego Desert State Park for permission to record data on Jackass Flat
and Dry Wash. The University of California Natural Reserve System Steele/Burnand Anza-
Borrego Desert Research Center Reserve DOI: 10.21973/N3Q94F provided field support. Frank
Vernon supervised the data collection. The seismic instruments were provided by the
Incorporated Research Institutions for Seismology (IRIS) through the PASSCAL Instrument
Center at New Mexico Tech. The data are available through the IRIS Data Management Center.
The facilities of the IRIS Consortium are supported by the National Science Foundation under
Cooperative Agreement EAR-1261681 and the DOE National Nuclear Security Administration.
The study was supported by the National Science Foundation (grant EAR-162061) and the
Department of Energy (award DE-SC0016520).

4.1 Introduction
Co-seismic and post-seismic temporal changes of seismic velocities associated with
earthquakes of different sizes and in different tectonic settings have been widely observed using
earthquake waveforms (e.g. Poupinet et al., 1984; Peng and Ben-Zion, 2006; Sawazaki et al.,
2009; Nakata and Snieder, 2012; Roux and Ben-Zion, 2014) and ambient noise cross correlations
(e.g. Brenguier et al., 2008b, 2014; Froment et al., 2013; Obermann et al., 2014; Hobiger et al.,
2016). Recent advances in deployments of dense seismic arrays around fault zones and analysis
87
techniques provide improved opportunities to track the temporal evolution of seismic velocities
within active fault zones. The San Jacinto fault zone (SJFZ) is one of the most seismically active
fault zones in southern California (Hauksson et al., 2012; Ross et al., 2017a), and is capable of
producing large Mw > 7.0 earthquakes (Petersen and Wesnousky, 1994; Onderdonk et al., 2013;
Rockwell et al., 2015). Several Mw > 5 earthquakes occurred in the SJFZ during the past decade.
The 2016 Mw5.2 Borrego Springs (BS) earthquake had a hypocenter at ~ 12 km depth in the
complex trifurcation area of SJFZ (Fig. 4.1a) where faulting is accommodated by three major
fault branches. Orthogonal structures are seen between the Clark and Coyote Creek faults both
from fault surface traces (Fig. 4.1) and seismicity patterns (e.g. Hauksson et al., 2012; Ross et al.,
2017a; Cheng et al., 2018). Analysis of source properties of the 2016 BS event indicate that it
had significant directivity to the northwest (Ross et al., 2017b). This inference is supported by
the strongly asymmetrical distribution of the aftershocks to the northwest of the mainshock
epicenter (Fig. 4.1).
Five linear arrays that consist of 6-12 stations with an average station spacing of ~ 20-40 m
were deployed across different sections of the SJFZ (Vernon and Ben-Zion, 2010). These include
two arrays, at Jackass Flat (JF) and Dry Wash (DW), which cross the main Clark branch of the
SJFZ about 2-3 km in both directions up and down the fault from the epicenter of the 2016 BS
earthquake (Fig. 4.1). In the present paper we use ambient seismic noise recorded by these two
linear arrays to analyze temporal changes of seismic velocities associated with the 2016 BS
event. Fault zone imaging results at these two sites (Qiu et al., 2017, 2019c) suggest that both
arrays cross the main seismogenic fault. However the fault zone structures in the top few
kilometers have different properties below the two arrays. Based on delay time analysis of body
waves, along with observation and modeling of fault zone trapped waves, the JF array to the
southeast of the BS event epicenter is situated in a low velocity zone with ~35% reduction in
shear wave velocity relative to the host rock (Qiu et al., 2017). In contrast, delay time analysis at
the DW array situated to the northwest and lack of clear trapped waves suggest a more localized
fault zone structure (Qiu et al., 2019c). Regional tomographic images of the SJFZ trifurcation
area (e.g. Allam and Ben-Zion, 2012; Allam et al., 2014a; Zigone et al., 2015; Fang et al., 2016)
suggest a velocity contrast across the Clark fault with the southwest side being slower in the top
few kilometers (background gray scale in Fig. 4.1a). This large-scale velocity contrast implies a
persistent preferred propagation direction of earthquake ruptures to the NW (e.g. Ben-Zion and
88
Andrews, 1998; Ampuero and Ben-Zion, 2008) that is consistent with the clear directivity of the
2016 BS event observed in Ross and Ben-Zion (2016) and Ross et al. (2017b). This
configuration of two across-fault linear arrays at the opposite along-strike directions of the
epicenter of the 2016 BS earthquake provides a unique opportunity to investigate how near fault
environment respond to a nearby moderate earthquake with rupture directivity.
Ambient noise correlation (ANC) techniques have been used extensively to approximate the
Green’s function of the subsurface between two stations (e.g. Campillo and Paul, 2003; Shapiro
and Campillo, 2004; Shapiro et al., 2005; Zigone et al., 2019). However, variations in noise
source distribution (e.g. Kedar and Webb, 2005; Stehly et al., 2006) can bias the reconstructed
direct waves. Code Wave Interferometry (CWI) uses coherent scattered arrivals in the coda part
of the cross correlation (e.g. Snieder et al., 2002; Grêt et al., 2006; Sens-Schönfelder and
Wegler, 2006; Brenguier et al., 2008a,b; Hillers et al., 2015a,b) and is less sensitive to biases
associated with changing source properties (Colombi et al., 2014). In this study, we use CWI to
estimate seismic velocity changes at different coda lapse times and in different frequency bands,
which are a first order proxy for the depth extent of the temporal changes (e.g. Obermann et al.,
2016; Yang et al., 2019).
In the first part of the paper we describe the data and preprocessing steps used to obtain
daily cross correlations. We then improve the data quality by using a denoising technique
(Moreau et al., 2017). To measure changes of seismic velocities δv/v, we use a time-domain
analysis (stretching) method, and further validate the resulting measurements with a frequency-
domain (doublet) method. For both methods, we estimate δv/v for different lapse times in the
coda (τ), component of the correlation tensor, and frequency band at each station pair. The
pairwise δv/v time series involving different stations are then averaged over all station pairs
within each array. We fit the resulting array-median δv/v assuming a combination of terms
including a linear trend, co- and post-seismic earthquake effects, and some periodic loadings that
may be related to natural sources such as tides and atmospheric temperature variations (section
4.4.1). The co-seismic and post-seismic velocity changes are highlighted by removing the
modeled linear and periodic terms from the fitted curve, we analyze their dependence on lapse
time, component, frequency band, and array location (section 4.4.2).

89

Figure 4.1. (a) Map of the fault zone linear arrays (DW in green, JF in red) used in this study
and a nearby TR array (yellow). The average P-wave velocity over the depth range 1–7 km based
on tomographic results of Allam and Ben-Zion (2012) is shown as the gray background (gray –
slow & white – fast). The blue star, beach ball, and green arrow indicate the epicenter location,
focal mechanism and rupture directivity of the 2016 Mw5.2 Borrego Springs earthquake.
Aftershocks detected and relocated by Ross et al. (2017a) are colored by depth, whereas the
green stars indicate Mw3.0+ aftershocks. The regional background seismicity during 2015 –
2016 is denoted as black dots. Surface traces of the San Jacinto fault (SJF) in the trifurcation
area, together with the main faults (e.g. San Andreas fault – SAF, Elsinore fault – EF) in
southern California (bottom left inset), are shown as black lines. The top and bottom right insets
contain the zoom in of the station configuration for the DW and JF arrays, respectively. (b)
Depth section of hypocenters projected along the cross-section AA’ with aftershocks colored by
the time relative to the main shock (largest yellow star). The shaded gray area represents a region
90
in which more than 0.5 m of slip occurred during the main shock according to Ross et al.
(2017b). The Mw3+ aftershocks are also shown as stars scaled with respect to the corresponding
magnitude, and the black dots indicate the locations of background seismicity occurred during
2015 – 2016. The green and red triangles denote the locations of JF and DW arrays, respectively.

For all frequencies, we find a strong inverse relationship between the co-seismic velocity
reduction associated with the 2016 BS event and the coda wave window used to measure δv/v.
The results suggest a much larger co-seismic perturbation of seismic velocities within the fault
zone compared to the surrounding host rock. We also observe a substantial difference in the
magnitude of co-seismic velocity reduction and post-seismic recovery time between the low-
frequency band 0.1-0.4 Hz and the other higher frequencies. This frequency dependence may
indicate seismic velocity changes that are associated with different mechanisms and depth
sections, such as earthquake rupture (deep) and seismic waves passing by (shallow). To further
compare seismic velocity changes beneath the DW and JF arrays, we scale the observed δv/v by
the corresponding amplitude of seasonal variations to eliminate the site-specific responses
associated with local fault zone structure. The resulting scaled co-seismic velocity reductions
show systematically larger values at DW array for both low (0.1-0.4 Hz) and high (1-4 Hz)
frequency bands, which is possibly related to the northwestward rupture directivity of the 2016
BS earthquake.

4.2 Data and processing
We use 2 years (2015-2016) of 3-component continuous waveforms with 200 Hz sampling
rate recorded by two linear fault zone arrays located in the trifurcation area of SJFZ. The JF array
consists of 9 broadband stations, and the DW array has 12 short-period sensors. Both arrays with
an aperture of ~400 m are located almost symmetrically to each side of the epicenter of the 2016
BS earthquake (Fig. 4.1) and are capable of recording ground motion from 0.1 Hz to 10 Hz. The
data from the JF array are analyzed up to November of 2016 when the array is taken out.
To retrieve converged Green’s Function estimates from ANC, transient signals generated by
earthquakes and other large episodic sources have to be removed in preprocessing steps (e.g.
Bensen et al., 2007). Here, we adopt the preprocessing procedure of Poli et al. (2012), Zigone et
al. (2015) and Hillers et al. (2015b) to compute daily ANC. This includes large amplitude
transients removal, spectral whitening and amplitude clipping. The 2016 BS main shock was
91
followed by high rate of aftershocks (Fig. A4.1). Trial and error analyses suggested an optimized
preprocessing with a window length of one-hour, 0.1 – 12 Hz frequency band for spectral
whitening, and a threshold of 3.5 times the standard deviation in amplitude truncation. To reduce
computation time, we downsample the data to 40 Hz prior to the preprocessing and compute
daily correlation functions by cross correlating the preprocessed data for all station pairs and
components within each fault zone linear array.

Figure 4.2. (a) Daily ambient noise cross-correlations com put ed f or year 2015 and 2016
between station pair JFS2 – JF00 after lowpass filtering with wiener filter at 4.5 Hz. The black
dashed line indicates June 10, 2016, the origin time of the main shock, whereas the green line
outlines the Julian day that all daily correlation functions computed before it are stacked to
construct the reference trace. The reference trace is shown in (c). The insets illustrate the zoom in
of the reconstructed coda waves within the black boxes around the main shock origin time. (b)
Same as (a) for station pair DW01 – DW03. Different from (a), an obvious difference is
observed in the travel time of scattered arrivals before and after the main shock, particularly in
the insets. (c) Reference traces computed at station pairs JFS2 – JF00 (black) and DW01 –
DW03 (blue). The insets show the zoom in of the waveform between lapse times of 6s and 40s
(dashed black boxes) at both time lags. The red dashed boxes in the insets represent the
correlation time range used to plot the insets in (a) and (b).

Remaining low quality daily correlations (Fig. A4.2) are removed in a two-step process.
First, we remove outliers based on the median absolute deviation (MAD) of the maximum
92
absolute amplitudes of all daily ANC. Correlation functions with maximum amplitudes three
times smaller or four times larger than the MAD from the median are removed. Second, we
construct a reference signal by stacking all daily ANC, and identify outliers through the
correlation coefficient of the reference and each daily traces. Since the central peak between lag
times of ±1s (|τ| < 1s) dominate the correlation, we first truncate the near-zero peak within the
window of |τ| < 1s and then cross correlate waveforms between lag times of ±60s. Daily
correlations are discarded if the maximum cross-correlation coefficient is smaller than 0.5. After
removing these low quality ANC, we rotate the 9-component cross correlation tensor into a
system of fault perpendicular (R), along fault strike (T), and vertical directions (Z).
Previous studies show that filtering (e.g. Hadziioannou et al., 2011, Stehly et al., 2015) and
sub-stacking at each date over a time interval of ±d days (e.g. Brenguier et al., 2008b) can
increase the signal to noise ratio (SNR) of the correlation functions. In this paper, we adopt the
denoising filtering technique developed by Moreau et al. (2017), which significantly improves
the SNR of each individual daily trace without substantially compromising the temporal
resolution. The technique uses a singular-value decomposition based Wiener filter to denoise the
daily ANC matrix. A 2-D low-pass Wiener filter with a size of 10 days and 5 samples (0.125s)
along the date and correlation lag time dimensions was found to yield the best results (Fig. 4.2).
We then stack the filtered daily correlations from the first day of 2015 to the Julian day that is
one month before the 2016 BS earthquake (green line in Fig. 4.2a&b) to generate the reference
trace (Fig. 4.2c). Ending one month before the 2016 BS event is to guarantee that velocity
perturbations caused by the main shock and aftershocks are not included in the reference, and to
not erroneously associate potential changes at the end of the stack period to medium changes.
To determine the part of the reference trace that is dominated by coda wave signals above
the noise level (i.e. coda wave window), we investigate its amplitude decay as a function of lag
time. We assume that, different from the amplitude of direct waves that is governed by both
geometrical spreading and attenuation, the coda waves are only affected by intrinsic and
scattering attenuations, whereas the base noise level should be independent of the lag time. We
further average the envelope functions over all available station pairs for the ZZ component and
analyze the median envelope. Figure 4.3 shows the best fitting amplitude decay of such median
envelopes at different frequencies as a function of lapse time. Based on the amplitude decay
fitting, we find the coda wave window to be 8s ≤ |τ| ≤ 40s for 0.1-0.4 Hz, 8s ≤ |τ| ≤ 35s for
93
0.5-2 Hz, 6s ≤ |τ| ≤ 23s for 0.75-3 Hz, and 5s ≤ |τ| ≤ 22s for 1-4 Hz. Here, we ignore the
frequency band 0.25-1 Hz due to very low SNR levels (Fig. A4.3). Figure 4.4 further validates
that coherent coda arrivals can be extracted from ANC for frequency bands from 0.1-0.4 Hz up
to 1-4 Hz within the determined coda wave window. Signals in either the entire or part of the
coda wave window can be used to measure δv/v time series.

Figure 4.3. (a) Waveform (blue) and the corresponding envelope function (black) of the one-
year stacked reference trace filtered at 0.1-0.4 Hz averaged over all station pairs in DW and JF
arrays for ZZ component. Fold and average is used to produce symmetric signals. The envelope
function is depicted in log scale. The green dashed curves represent an amplitude decay model
that consists of both geometrical spreading and attenuation terms, whereas the red dashed lines
depicted a model that only includes attenuation effect. Based on the fitting of both models, the
black dashed lines outline the identified coda wave windows. Here, we truncate the coda wave
window at lapse time of 40s to ensure sufficient signal amplitude. Zoom in of the coda
waveforms (6s-40s) are shown in red. The envelope function of individual station pair is shown
as the background gray curves. (b) Same as (a) for frequency band 0.5-2 Hz, the coda wave
window is defined as 8s-35s. (c) Same as (a) for frequency band 0.75-3 Hz, the coda wave
window is defined as 6s-23s. (d) Same as (a) for frequency band 1-4 Hz, the coda wave window
is defined as 5s-22s.

4.3 Methodology
4.3.1 Relative Velocity Change Estimates
94
For a scattered arrival reconstructed from ANC, let t be the absolute travel time observed in
the reference trace, and δt denotes the difference in travel time measured between an instance of
a daily correlation and the reference. The relative velocity change (δv/v) is given by the
relationship (e.g. Ratdomopurbo and Poupinet, 1995; Snieder et al., 2002)
δυ
υ
=−
δt
t
. (1)
Instead of using individual scattered arrivals, a coda wave window is commonly used to
stabilize the estimate. Since scattered arrivals with greater travel time t have a correspondingly
longer propagation path, the resulting estimations represent velocity perturbations averaged over
larger areas surrounding the interstation path. Assuming a homogenous medium change and a
linear relation between δv and v, the obtained δv/v can be interpreted as the temporal perturbation
in medium properties for the sampled region.

Figure 4.4. Coherent scattered arrivals observed in coda waves between lapse time 8s and 15s at
both negative (a) and positive (b) for a 120-day period starting from 60 days before the 2016 Mw
5.2 BS earthquake in the 0.1-0.4 Hz frequency band. The scattered arrivals reconstructed from
the ANC are fairly symmetric between both time lags. Red dashed lines at zero-day offset
indicate the occurrence time of the earthquake. (c), (d), same as (a), (b) but for frequency band
0.5-2 Hz and 30 days before and after the event. The reconstructed coda waves are highly
asymmetric. (e), (f), same as (c), (d) for frequency band 1-4 Hz, and strong asymmetry in coda
waves are also observed. Coherent scattered arrivals are observed throughout the time period for
all three frequency bands.
95
To robustly and reliably estimate δt/t from coda waves, the time-domain stretching method
(e.g. Sens-Schönfelder and Wegler, 2006) and the frequency-domain doublet method (e.g.
Poupinet et al., 1984) are often used. The stretching method dilates or compresses coda
waveforms in time domain, whereas the doublet method uses the moving window cross-spectral
analysis (Clarke et al., 2011) to estimate δv/v. Both methods compare a test trace to a reference
trace, and assume that temporal perturbations in coda wave travel times are small. In this paper,
we use the denoised daily correlation as the test trace, and generate the reference trace as
described in section 4.2.
Both methods have their advantages and drawbacks. Spurious velocity change estimates may
be observed using the stretching method when the coda wave amplitude spectrum varies with
time (Zhan et al., 2013). In contrast, the doublet method provides phase-based estimates and is
less sensitive to the temporal variations in noise amplitude spectrum (Hadziioannou et al., 2009;
Colombi et al., 2014). However cycle skipping, in which phases are matched incorrectly between
the coda waves in the test and reference traces, is a potential problem for the doublet method
(Mikesell et al., 2015). Since these two methods behave differently in the presence of
fluctuations, consistency between estimations obtained from both techniques is indicative of the
robustness of the resulting δv/v.
There are several different noise excitation patterns in the SJFZ environment: microseisms
coming excited offshore and along the coastline for the low-frequency band 0.1-0.4 Hz, and
cultural (e.g. Inbal et al., 2018) and other natural sources (Hillers et al., 2013) for higher
frequencies (e.g. > 1 Hz). Although the noise sources are not evenly distributed for frequencies
above 0.1 Hz in the study region (e.g. Hillers et al. 2013), the coda waves in ANC are fairly
symmetric in the microseism band 0.1-0.4 Hz (e.g. Figs. 4.4a,b). This symmetry in reconstructed
coda waves at this frequency band results from interferometry of multiple-scattered ambient
noise in the heterogeneous Southern California plate boundary region. However the coda waves
are less symmetrical at higher frequencies (e.g. Figs. 4.4c-f). To minimize the effect of non-
isotropic noise sources, we use the coda waves from both positive and negative time lags
simultaneously to determine δv/v.
In contrast to δv/v measured from networks with an aperture of several kilometers, the linear
arrays used here are situated within fault zones that may yield larger temporal changes compared
to the surrounding host rock in response to a loading signal. For such localized temporal changes
96
in rock properties, δt/t tends to be larger at early scattered arrivals and decreases with lapse time
(e.g. Obermann et al., 2016), because coda waves arriving at later lapse times averaged the
seismic velocity variations over a larger area including host rock regions that are less affected
compared to the near fault material. As a result, the underlying assumption of the doublet method
that δt/t is independent of the lapse time t and can be estimated through linear regression in the
frequency domain is not fully compatible with the observation. Although the localized medium
changes may also affect the estimation from stretching method, the time domain method yields a
reliable δt/t estimate of the largest amplitude coda arrival.

Figure 4.5. Coherence and δv/v obtained for (a) DW array using stretching method, (b) DW
array using doublet method, (c) JF array using stretching method, and (d) JF array using doublet
method. Each grey curve represents the resulting CC
max
in the top panel and δv/v in the bottom
panel for each station pair within the array. The blue solid and dashed curves denote the array-
median and corresponding uncertainty of all the grey curves. The median curve is further
smoothed via a median filter and shown as the red curve in the bottom panel for better
visualization. Coda waves from ZZ component within the window from 5s to 22s (Fig. 4.4d)
filtered at 1-4 Hz are used. The black and red dashed lines mark the date of June 10, 2016, the
occurrence time of the 2016 Mw 5.2 BS earthquake and the end of the stacking range to
construct a reference trace, respectively.

97
4.3.2 Curve fitting
Figure 4.5 shows example time series of δv/v and the corresponding coherence estimated
using both the stretching and doublet methods at the DW and JF arrays at 1-4 Hz. In this section,
we model the obtained δv/v time series considering three main components: a long-term linear
trend, quasi-periodic variations, and co-seismic velocity reduction followed by a recovery
process associated with the 2016 BS earthquake. To evaluate and quantify these components, we
approximate the δv/v time series with the following form:
f T
( )
= f
lin
T
( )
+ f
per
T
( )
+ f
eq
T −T
0
( )
. (2)
T is the number of days from the first day of analysis period (2015-01-01) to each target date. T
0
is the origin time of the 2016 BS event. Here f
lin
denotes the linear component and is given by
f
lin
T
( )
= A+B ⋅T. (3)
Since the co-seismic and post-seismic changes in seismic velocities are the main focus of this
paper, we only include the semi-annual (e.g. Johnson et al., 2017) and annual variations (e.g.
Hillers et al., 2015b) in the curve fitting. The quasi-perodic component f
per
can be written as:
f
per
T
( )
=Ccos ω T −T
1ω
( )
⎡
⎣
⎤
⎦
+Dcos 2ω T −T
2ω
( )
⎡
⎣
⎤
⎦
, (4)
with C, D as the amplitudes, and T
1ω
, T
2ω
as the phase terms of the annual and semi-annual
velocity variations. Here ω is a fixed value of 2π/year. The third term in equation 2, f
eq
,
represents the functional form of the co-seismic and post-seismic velocity changes combined,
which can be approximated by the summation of a Heaviside function H(T) and a logarithmic
recovery:
( ) ( ) ( ) ( ) T H T E T H E T f
eq
⋅ + ⋅ + ⋅ − = 1 ln
2 1
, (5)
where E
1
represents the co-seismic velocity reduction and E
2
/E
1
constrains the post-seismic
recovery time (Fig. A4.4a). Eq. (5) is similar to the expression used by Hobiger et al. (2012) with
added 1 to regularize the log function near zero time. Laboratory data indicate that log(t) healing
already operates at 1 s after failure (e.g., Nakatani and Scholz, 2004). Since the minimum
timescale in the preformed analysis is 1 day (~8.6×10
4
s), and the examined post-seismic period
includes only about 6 months, the results do not cover the initial 4-5 decades of the healing
process and do not include a long time intervals (in log scale) of the later post-seismic recovery.
In this case it is convenient to utilize an exponential function to model the post-seismic δv/v, and
derive an effective recovery time scale based on that function (Hobiger et al. 2014). Given the
98
fluctuating character of the results (Fig. 4.5) and the short duration of the available data, we use a
simple exponential function to approximate the post-seismic recovery (Taira et al., 2018):
f
eq
T
( )
= −E⋅H T
( )
⋅e
−T T
c
, (6a)
where E and T
c
represent the resolved co-seismic velocity reduction and post-seismic recovery
time (i.e. 2.3T
c
for 90% recovery). As the Heaviside step function has sharp corners, which is
inconsistent with the smooth δv/v curve estimated from the noise-filtered correlations, we further
refine equation 6a by using the arctan function:
f
eq
'
T
( )
= E ⋅ arctan T λ−T
s
2λ
( )
π +1 2
⎡
⎣
⎤
⎦
⋅e
−H (T−T
s
)⋅ T−T
s
( )/T
c
, (6b)
where λ and T
s
(> 0) characterize the shape of the arctan function (Figs. A4.4b-d). Here
e
−H (T−T
s
)⋅(T−T
s
)/T
c
represents the post-seismic recovery that equals to unity at T ≤ T
s
and e
− T−T
s
( ) T
c
at
T > T
s
. Non-recovering velocity reductions after large earthquakes caused by material damage
have been observed and modeled in previous studies (e.g. Peng and Ben-Zion, 2006; Hobiger et
al., 2016). However, this is not included in the curve fitting analysis as the recovery of seismic
velocities hasn’t completed in the accessible 6 month post-seismic period for frequencies > 0.5
Hz (section 4.5).
We estimate A, B, C, D, T
1ω
, T
2ω
, E, and T
c
following the above equations 2, 3, 4, and 6b, and
the corresponding uncertainties for each array-median δv/v using a nonlinear least square curve
fitting algorithm from the SciPy package (e.g. Figs. 4.6 & 4.7). After obtaining the best fitting
parameters, the linear trend (f
lin
) and annual & semi-annual variations (f
per
) are subtracted from
the obtained δv/v to highlight the co- and post-seismic responses (e.g. Figs. 4.8-4.14). To further
evaluate the significance of the earthquake response modeling part in the curve fitting, we
compute a percentage of variance reduction as
δσ = σ
b
−σ
e
( )
σ
b
, (7)
where σ
b
and σ
e
are the variances computed for the curve fitting with and without considering
co- and post-seismic velocity changes
σ
b
= var f
obs
s
T
( )
− f
lin
T
( )
− f
per
T
( )
⎡
⎣
⎤
⎦
σ
e
= var f
obs
s
T
( )
− f
lin
T
( )
− f
per
T
( )
− f
eq
'
T
( )
⎡
⎣
⎤
⎦
. (8)
In section 4.1, we first model all the array-median δv/v time series measured using different
combinations of method, time window, frequency band, component, and array, using the
99
equations developed in this section, but mainly focus on the ZZ component results. The curve
fitting suggests that the δv/v time series can be well explained by a combinations of seasonal
variations and a co-seismic velocity reduction associated with the 2016 BS event followed by a
post-seismic recovery, for all frequency bands, arrays, and components.

Figure 4.6. Array-median δv/v (blue) and the curve fitting results (red) of DW array at ZZ
component using stretching method for frequency bands (a) 0.1-0.4 Hz, (b) 0.5-2 Hz, (c) 0.75-3
Hz, and (d) 1-4 Hz. The gray shaded area illustrates uncertainties of the array-median δv/v. The
frequency dependent coda wave windows indicated by Fig. 3 are used in the curve fitting at
different frequency bands.

We then remove the modeled non-earthquake (non-EQ) responses from the observed array-
median δv/v time series to highlight the co- and post-seismic (EQ) responses. In section 4.4.2.1,
we compare these “highlighted” array-median δv/v time series estimated using successive
moving sub-windows within the coda wave windows indicated in Figure 4.3. In contrast to the
seasonal variations that are independent of lapse time window, the obtained lapse time dependent
co-seismic velocity reduction is consistent with the hypothesis that the fault damage zone is more
susceptible to nearby moderate earthquake than the surrounding intact rock.

4.4 Results
Figure 4.5 shows the resulting δv/v obtained from ZZ correlations filtered at 1-4 Hz using
100
coda waves in the time window of 5-22s (see Figs. A4.5-A4.7 for other frequency bands). To
ensure the reliability of the obtained temporal changes, for each station pair, both the stretching
(Figs. 4.5a,c) and doublet methods (Figs. 4.5b,d) are used. Results of δv/v from both methods
show the same systematic temporal patterns. The results differ in the coherence values (lower for
the doublet method) and in the magnitudes of the δv/v variations and the corresponding
uncertainties. Systematic higher coherence values from the stretching method suggest that the
estimates are taken using the optimally dilated waveforms. Many of the results from the doublet
method used for comparison with stretching method are shown in the supplementary materials.

Figure 4.7. (a) Best fitting parameter C, amplitude of the seasonal variation, as a function of
frequency band. Results obtained using stretching and doublet methods are shown as dashed
(triangles) and solid (circles) lines, whereas estimates for DW and JF array are displayed in red
and blue, respectively. C values at the low-frequency band 0.1-0.4 Hz are shown in a larger scale
(left axis) compared to those at frequencies above 0.5 Hz (right axis). Red ticks on the left y-axis
illustrate the scale on the right. (b) Best fitting parameter T
1ω
, phase term of the seasonal
variations, as a function of frequency. The black dashed lines indicate the Julian day 200, at
which the annual surface temperature in SJFZ (Hillers et al., 2015b) reaches the maximum. (c)
Same as (a) for parameter E (co-seismic velocity reduction). (d) Same as (b) for parameter T
c

(effective post-seismic recovery time).

Considering the small aperture (~400 m) of the two fault zone arrays, the purpose of this
section is to use the array-median δv/v time series obtained from all station pairs within each
array to infer a robust estimation of the seismic velocity perturbation in the underlying fault
101
zone, and analyze how the observed temporal patterns associate with the 2016 BS event and
other natural loads such as temperature variations. However, similar to the Frechet kernel for
finite frequency travel times (e.g. Dahlen et al., 2000), the sampled volume and thus the
sensitivity of coda waves to velocity perturbations in a heterogeneous medium increases both
horizontally and vertically with lapse time (e.g. Obermann et al., 2013, 2016). The sensitivity
halo for coda waves thus likely exceeds the size of fault damage zone, particularly at low
frequencies (e.g. 0.1-0.4 Hz) and longer lapse time (e.g. 30-40s). A study of δv/v measured for
coda waves at different lapse times can provide constrains on rock susceptibility to perturbations
in fault damage zones and the surrounding host rock.

4.4.1 Results of curve fitting
Figure 4.6 illustrates the curve fitting results of the array median δv/v time series estimated
for the DW array using the stretching method (Fig. A4.8 for the doublet method; Figs. A4.9,
A4.10 for the JF array). ZZ component correlations and the frequency dependent coda wave
windows as indicated by Figure 4.3 are used for the demonstration here. The best fitting
parameters, such as amplitudes (C, D) and phase terms (T
1ω
, T
2ω
) of seasonal and semi-annual
components, co-seismic velocity reduction (E), and post-seismic recovery time (T
c
) are indicated
in the figure. The modeled δv/v time series agree well with the measurements for all the used
frequency bands, arrays, components, and lapse time windows. The resulting δσ (Eq. 7) implies
the co- and post-seismic components are essential to fit the observed δv/v time series, and the
fitness generally increases with frequency. In this section, we mainly focus on the best fitting
parameters shown in Figure 4.7.
At 0.1-0.4 Hz, we find large amplitude periodic (30-40 days periodicity) fluctuations in δv/v
between June and October of 2015 (e.g. Fig. 4.6a), particularly at the DW array. Peak-to-peak
amplitude of the oscillation is comparable to the observed co-seismic reduction in δv/v for both
arrays. However, based on the earthquake catalog of Hauksson et al. (2012), there are no Mw >
3.0 earthquakes in the study region (Fig. 4.1a) between June and October of 2015. This large
short-period background fluctuations correlate well with the visible waveform perturbations in
daily ANC computed for the station pair DW01 – DW03 (Fig. 4.2b, between Julian day 100 and
300). However, we do not observe fluctuations on that time scale in the higher frequency bands
between 0.5 Hz and 4 Hz. Combined with the fact that the amplitude of this short period signal
102
decreases with lapse time (section 4.4.2.1; Figs. A4.11a, b), we conclude that this signal in the
microseisms frequency band is caused by spurious variations in the microseisms excitation
pattern (Hillers et al., 2015b).

4.4.1.1 Seasonal variations – parameters C and T
1ω

The quasi-periodic term f
per
consists of two parts: semi-annual and seasonal variations. The
amplitude of semi-annual component (D) is comparable to the standard deviation of the residual
fluctuations. The semi-annual component is thus negligible and we focus on seasonal velocity
changes (C and T
1ω
). Best fitting values of C and T
1ω
estimated from the DW and JF arrays using
both methods are shown in Figs. 4.7a and 4.7b, respectively.
There are two broad possible causes for the derived seasonal variations in seismic velocities:
changes in medium properties, such as crack density or fluid content, or variations caused by
wavefield properties that are considered spurious (Zhan et al., 2013). To exclude the latter
possibility and identify the dominant driving mechanism for seasonal variations in the SJFZ,
Hillers et al. (2015b) performed lapse time analysis and investigated the delay between the
obtained seasonal variations and various environmental records that exhibit seasonal or semi-
annual periodicity, such as atmospheric temperature, wind speed, and precipitation. They
concluded that the δv/v variations in the SJFZ at 0.1-0.4 Hz are likely caused by wavefield
changes. For the frequency ranges of 0.5–2 Hz and 1–4 Hz, however, the results more likely
reflect the medium response. According to the delay between the observed seasonal variations of
δv/v and air temperature, the most plausible primary mechanism was inferred to be thermoelastic
strain (e.g., Berger, 1975; Ben-Zion and Leary, 1986).
Similarly, we perform lapse time analysis with the corresponding results shown in Figs.
A4.11 and S12, and discuss the δv/v signals and the phase-delay derived at the low-frequency
band 0.1-0.4 Hz and other higher frequency bands (> 0.5 Hz) separately. In the following
discussion, we only investigate the delay between the observed seasonal variations and the
annual air temperature. Based on the surface temperature record at Piñon Flats near the SJFZ
(Hillers et al., 2015b), we find that the air temperature reaches maximum at ~200 day of the year.
Therefore, we use day of the year T
max
= 200d as the peak time of the annual surface temperature
at both array locations to roughly estimate the delay between the δv/v signals and the air
temperature.
103
At the low-frequency band 0.1-0.4 Hz, we find a large deviation between C values (Eq. 4)
achieved either from different methods (circles versus triangles in Fig. 4.7a) or arrays (red versus
blue in Fig. 4.7a). The differences in C values derived from the stretching and doublet methods
indicate that the seasonal signals extracted from time- and frequency-domain are inconsistent.
The systematic larger C value at DW array may reflect differences in sensor types and in-situ
conditions between the two arrays. In addition, the C values at 0.1-0.4 Hz yield significantly
larger seasonal variations compared to the other analyzed higher frequency bands (Fig. 4.7a).
These observations are consistent with the conclusion of Hillers et al. (2015b) that seasonal
variations at 0.1-0.4 Hz in the SJFZ region likely reflect changes in microseism excitation pattern
rather than rock responses.
In contrast to C values, very similar T
1ω
values (Eq. 4) are obtained at the low-frequency
band for both methods (Fig. 4.7b): the seasonal variations peak at day of the year T
1ω
≈ 200d and
150d for the DW and JF arrays, respectively. The fact that the obtained T
1ω
values are equal or
smaller than T
max
(i.e. 200d) rules out thermoelastic strain as the dominant mechanism for the
observed seasonal variations, in agreement with our results on C values. However, consistent
seasonal variations are observed for all lapse times at 0.1-0.4 Hz (Figs. A4.11a&b), in contrast to
the fading annual signals at longer lapse times in the low frequency band shown in Hillers et al.
(2015b). The cause for why the seasonal variations measured at longer lapse times do not decay
for 0.1-0.4 Hz in this study is yet unclear.
At higher frequency bands between 0.5 Hz and 4 Hz, the values of C obtained from the
stretching and doublet methods at both arrays are ~0.05% (Fig. 4.7a). This magnitude of C value
obtained from different frequency bands and processing parameters agrees with the seasonal
variations in the SJFZ environment at 0.5-2 Hz reported in Hillers et al. (2015b) (C ≈ 0.1%).
Combined with the fact that consistent annual variations are retrieved for all lapse time windows
(Figs. A4.11 and A4.12), the seasonal variations at these higher frequency bands likely reflect
the medium response at the fault zone sites.
As shown in Fig. 4.7b, consistent T
1ω
values are achieved using both the stretching and
doublet methods. At the DW array, T
1ω
≈ 240d for all three frequency bands, whereas the T
1ω

values are larger and in general decreases with frequency for JF array. Using the approximation
of T
max
= 200d, we find a ~40-day delay between the obtained seasonal δv/v at the DW array and
the annual surface temperature. This average 40-day delay is expected for a thermoelastic effect
104
in homogenous elastic half-space (e.g. Berger, 1975; Tsai, 2011; Richter et al., 2014). At the JF
array, the delay is slightly longer compared to that at the DW array, which may reflect a thicker
soil layer, and decreases with increasing frequencies (~ 120 days at 0.5-2 Hz to ~ 60 days at 1-4
Hz; Fig. 4.7b). We conclude that thermoelastic strain is the likely primary mechanism for the
observed seasonal variations in δv/v at frequency bands between 0.5 Hz and 4 Hz. The
decreasing delay time with increasing frequency at JF array is consistent with the fact that coda
waves at lower frequencies are generally more sensitive to changes in deeper structures,
associated with somewhat thicker effective unconsolidated layer for thermoelastic strain signal
(e.g. Ben-Zion and Leary, 1986; Prawirodirdjo et al, 2006). The lack of such observable trend in
delay time at DW array may reflect different local fault zone structures beneath the arrays (Qiu et
al., 2017, 2019c).
Since no strong lapse time dependence is found in seasonal variations (Figs. A4.11 and
A4.12), in the later section 4.4.2 we subtract the model predicted “non-EQ responses” (e.g. linear
trend and seasonal variations) from the δv/v curve to highlight the co-seismic and post-seismic
responses, and refer to the corrected δv/v curve as the “EQ responses”.

4.4.1.2 Co-seismic velocity reduction – parameter E
The co-seismic velocity reduction, characterized by a near instantaneous drop, stands out
from the background long period variations in both the raw δv/v time series and curve fitting
results (Fig. 4.6). The parameter E (Eq. 6b) indicates the best fitting estimate of co-seismic
velocity reduction before and after T
0
, the time of 2016 BS event. In general, the co-seismic
component is well fitted in all cases (Figs. 4.6 and A4.8 – A4.10), but the results at the DW array
yield less reliable co-seismic estimates due to the poor fitting of a large amplitude high
frequency oscillation right after the event origin time T
0
. This narrow spike in the δv/v curve with
a width of ~ 2 week may reflect influences of the aftershock sequences (Fig. A4.1), which can
include both medium changes and, in high frequencies, a potentially bias associated with the
wavefield variation. Since most of the aftershocks are located beneath the DW array, it is
consistent with the observation that such narrow spike only appears in the DW array results.
Considering the uncertainty and biases in estimates of co-seismic velocity reduction E,
particularly at the DW array, we focus on the comparison of E values derived from the different
methods. At the low frequency band 0.1-0.4 Hz, for both arrays we find the E estimated from
105
stretching method is ~2 times as large as that of doublet method. The values of E are > 1% at
DW array and between 0.2% and 1% at the JF array (Fig. 4.7c). At frequency bands between 0.5
Hz and 4 Hz, ~0.2%-0.3% value of E is observed in results for both arrays and methods (Fig.
4.7c), which is comparable to previously observed co-seismic velocity reduction in a similar
frequency range (e.g. Brenguier et al., 2008b; Hobiger et al., 2016; Taira et al., 2018). Both
observations suggest that the co-seismic reduction in δv/v curve is a robust signal associated with
subsurface velocity decreases at the time of the 2016 BS earthquake.

4.4.1.3 Post-seismic velocity recovery – parameter T
c

In the curve fitting process, the characteristic time T
c
(Eq. 6b) quantifies the post-seismic
response with a larger value indicating a slower recovery. As discussed in section 4.2, we have
less than 6 month of data to constrain the post-seismic δv/v. When the characteristic time T
c
is
comparable to or larger than 6 months (i.e. ~180 days), T
c
could not be well constrained through
curve fitting and often yields large uncertainty (Fig. 4.7d). In addition, the accuracy of the
estimated post-seismic recovery time T
c
also depends significantly on the quality of the best
fitting co-seismic velocity reduction E, which is poorly fitted for the DW array results due to the
high frequency oscillation immediately after the BS event. Thus, in this section, we only discuss
the reliability of T
c
based on the comparison of results obtained from different methods. In
general, both stretching and doublet methods show compatible trends in Fig. 4.7d. At the low
frequency 0.1-0.4 Hz, T
c
is ~10 days for both arrays and methods. However, T
c
estimates at the
higher frequency bands are much larger. At the JF array, both methods suggest T
c
≈ 110 days for
all frequency bands from 0.5-2 Hz to 1-4 Hz. At the DW array, T
c
values are larger than those of
JF array and increase with frequency (~200 days at 0.5-2 Hz to ~400 days at 1-4 Hz) for both
methods.
Figure 4.7d shows that the uncertainty estimates of T
c
are significantly smaller at the JF array
(~ 30 days) than those at the DW array (> 100 days) at frequency bands above 0.5 Hz. Due to
strong ground motion and non-linear effects of the shallow structures, Hobiger et al. (2016)
found non-negligible permanent weakening in seismic velocities after strong earthquakes (Mw >
6.0) in Japan for frequencies above 0.25 Hz. Therefore, missing the non-recovering co-seismic
velocity reduction term in the curve fitting (Eq. 6) for results at DW array may explain the
observed larger values in both T
c
estimates (> 200 days) and the corresponding uncertainties (>
106
100 days), particularly for the 1-4 Hz band. Since the co-seismic velocity reduction for most of
the δv/v curves estimated at frequencies above 0.5 Hz are not fully recovered at the end of the
analyzed period, the estimates on the post-seismic recovery time T
c
from equation 6b are not
well constrained.

4.4.2 Comparative analyses
In section 4.4.1, curve fitting of δv/v time series achieved from different methods and
frequency bands indicate robust co- and post-seismic responses associated with medium changes
caused by the 2016 BS earthquake. In this section, we remove the modeled linear trend and
seasonal variations to highlight the earthquake related variations (i.e. EQ responses) in δv/v. The
EQ responses resolved using different processing parameters including lapse time (section
4.4.2.1), frequency band (section 4.4.2.2), correlation component (section 4.4.2.3), and array
location (section 4.4.2.4) are directly compared below.

Figure 4.8. (a) Array-median coherence (top panels) and EQ responses of δv/v curves (bottom
panels) measured using different coda wave windows (colored curves) at 0.1-0.4 Hz. (b) Same as
(a) for JF array at 0.1-0.4 Hz. (c) Same as (a) for DW array at 1-4 Hz. (d) Same as (a) for JF
array at 1-4 Hz. Stretching method is used here.

0.1 - 0.4 Hz 1 - 4 Hz
107
4.4.2.1 Lapse time analysis
Figure 4.8 shows comparisons of the EQ responses of δv/v and corresponding coherences
obtained using the stretching method and several 10s-long moving time windows within the coda
wave window as indicated by Figure 4.3. We use frequency bands 0.1-0.4 Hz and 1-4 Hz for the
lapse time analysis as these results are clearest. In general, we observe very high consistency in
co- and post-seismic responses from results for all lapse time windows. However, the co-seismic
velocity reduction decreases monotonically with lapse time of the moving window. To further
investigate the rate of decreases in the co-seismic reduction with lapse time, we analyze the E
values (Eq. 6b) estimated through curve fitting, which reflect the co-seismic medium changes, as
a function of lapse time in Figs. 4.9a, b.

Figure 4.9. Lapse-time dependent co-seismic velocity reduction and median coherence value
derived from curves shown in Figure 4.8. (a) Co-seismic velocity reduction as a function of lapse
time window. Results at 0.1-0.4 Hz for DW and JF array are shown in red and blue, respectively.
(b) Same as (a) for 1-4 Hz. (c) Same as (a) for coherence values. (d) Same as (c) for 1-4 Hz.

For the low-frequency band 0.1–0.4 Hz, the rate at which the value of E decreases with lapse
time is much higher at DW than JF array, particularly at early lapse times. The rate is
approximately the same for both arrays at 1-4 Hz. Consistent with the observation of smaller
SNR for coda waves at later lapse time (Fig. 4.3), coherence values generally decrease with lapse
time (Figs. 4.9c, d).
We also compare the post-seismic recovery obtained for different lapse times. Consistent
with the observed δv/v shown in Fig. 4.6a, the co-seismic velocity reduction is fully recovered
with ~20 days after the 2016 BS event for all lapse times at 0.1-0.4 Hz (Figs. 4.8a, b). This rapid
108
recovery likely reflects the fast post-seismic healing process in fault zones due to larger normal
stress at greater depth. At 1-4 Hz, the co-seismic reduction is not fully recovered due to the
limited period of recording (Figs. 4.8c, d). The longer recovery time at higher frequency is also
consistent with the smaller normal stress at shallower depth that that delays the post-seismic
healing process.

4.4.2.2 Frequency band
Figure 4.10 shows the comparison of the EQ responses estimated at different frequency
bands for both DW and JF arrays. As indicated by the best fitting parameter E and T
c
that
quantify the co-seismic velocity reduction with an arctan function and post-seismic recovery
with an exponential function, respectively, in section 4.4.1.2, the co-seismic velocity reduction
and post-seismic recovery at 0.1-0.4 Hz are significantly larger and faster than those at
frequencies above 0.5 Hz, respectively.
To interpret the results, we use the assumption that coda waves at lower frequency are more
sensitive to changes in deeper structures. This assumes a dominance of surface wave. However,
Yang et al. (2019) noticed a significant effect of shallow seismic velocity changes on phase
velocities of long period (e.g. >15s) surface waves. Therefore, one cannot simply project the
observed δv/v at a low frequency band (e.g. 0.1-0.4 Hz) to deeper structure without considering
the measurements at higher frequency bands (e.g. 1-4 Hz). In this section, we provide a rough
estimation on the depth of co-seismic medium changes, at which coda waves in each frequency
band are the most sensitive to, and validate the assumption that the co-seismic δv/v observed at
lower frequencies are representative to seismic velocity perturbations in deeper structures.
To estimate the depth range of the observed frequency dependent co-seismic δv/v reduction,
we need to consider the sensitivity of coda waves to local medium changes at depth. The coda is
a mixture of body and surface waves in most crustal-like materials, and the partitioning between
their contributions to the coda sensitivity depends on the heterogeneity of the medium and the
coda lapse time (e.g. Obermann et al., 2013, 2016). For a first-order estimation, we ignore the
contribution of body-wave sensitivity and assume the observed δv/v involves medium changes at
depth of approximately one third of the dominating wavelength (i.e. h
max
≈ λ/3; using the depth
sensitivity of surface wave phase velocity to shear wave velocity).

109

Figure 4.10. Comparison of the EQ responses of δv/v measured at four different frequency bands
(colors) for the DW (left panel) and JF (right panel) arrays. ZZ component and stretching method
are used. Median filters are applied to smooth and suppress the background fluctuations in the
raw δv/v curve (transparent curves). The resulting smoothed time series of δv/v are shown as
solid curves. The right vertical axis in red denotes the scale of the red curves measured in
frequency band 0.1-0.4 Hz.

At 0.1-0.4 Hz, the dominating frequency is ~0.14 Hz (primary microseisms), and the
corresponding phase speeds in fault zones at the DW and JF sites (e.g. Roux and Ben-Zion,
2017; Berg et al., 2018) are ~3.2 km/s for Rayleigh wave (e.g. ZZ component) and ~3.4 km/s for
Love wave (e.g. TT component). Therefore, the apparent δv/v is most sensitive to medium
changes at the depth of h
max
= ~7-8 km. Similarly, for frequency bands above 0.5 Hz, h
max
is
smaller than 1 km as the corresponding surface wave speeds are generally slower than 2 km/s.
As shown in section 4.1.2, we observe the largest co-seismic velocity reduction at the low-
frequency band 0.1-0.4 Hz (> 0.5%), and the post-seismic recovery is significantly faster, which
rules out the possibility that the observed co-seismic velocity reduction at the low-frequency
band is biased by medium changes in shallow structures.

4.4.2.3 Components
In the case of surface wave reconstructed from ANC, the recording component is indicative
of the type of surface wave (e.g. ZZ – Rayleigh wave and TT – Love wave). The coda, however,
is a mixture of scattered body and surface waves and thus does not have a simple relation
between waveforms recorded at different components.

110

Figure 4.11. Comparison of the smoothed EQ responses of δv/v measured using different
components (colors) for DW array and stretching method at frequency bands (a) 0.1-0.4 Hz, (b)
0.5-2 Hz, (c) 0.75-3 Hz, and (d) 1-4 Hz.

Figure 4.11 illustrates the EQ responses estimated at different components using the
stretching method for DW array (Fig. A4.13 for doublet method; Figs. A4.14, A4.15 for JF
array). Except for the low frequency RT and ZT results, we observe generally comparable rapid
co-seismic velocity reductions followed by a post-seismic recovery associated with the 2016 BS
earthquake for all components in each frequency band. Although the results include some
variations in co- and post-seismic responses estimated at different components, there is no clear
component-related pattern that can be easily interpreted. In general, the comparison of results
obtained from different components indicates higher similarity at frequencies above 0.5 Hz than
that at the 0.1-0.4 Hz. This is understood to the result from coda waves at high frequencies
contain high degree of mixing modes, i.e. multiple scattering, whereas the mixing is not as
complete at the low-frequency band (Hillers et al., 2013), due to the proximity of the study area
to excitation region and the longer smoothing paths.

4.4.2.4 Array location
To compare the co-seismic medium changes at the DW and JF sites using δv/v estimated
from the two arrays, we remove effects of local structures in the observed co-seismic δv/v. Since
111
both arrays are located in the SJFZ and the inter-array distance is only ~6 km, we assume that the
amplitude of seasonal variations (C), either governed by medium or wavefield changes, is similar
if the local conditions are identical at both sites. Therefore, to remove the local response, we
normalize the EQ responses by the corresponding best fitting C value.

Figure 4.12. Comparison of the scaled EQ responses of δv/v measured at 1-4 Hz for the DW (red)
and JF (blue) arrays. (a) The transparent and solid curves represent the raw and smoothed array-
median δv/v obtained using stretching method at ZZ component, respectively. Each δv/v curve is
normalized by its corresponding best fitting C value (seasonal variation amplitude). (b) Same as
(a) for RR (fault perpendicular) component. (c) Same as (a) for TT (fault parallel) component.
(d), (e), (f) same as (a), (b), (c) for doublet method.

Figure 4.12 shows the comparison of such normalized EQ responses estimated for the DW
(blue in Fig. 4.12) and JF (red in Fig. 4.12) arrays at 1-4 Hz. The top and bottom panels represent
results using the stretching and the doublet method, respectively. Robust and systematic larger
co-seismic reductions are observed at DW compared to those at JF for both methods and all three
analyzed components. Larger co-seismic reduction at DW is also clearly observed at the low
frequency band 0.1-0.4 Hz (Fig. 4.13), but less clear at the other two frequency bands 0.5-2 Hz
(Fig. A4.16) and 0.75-3 Hz (Fig. A4.17). This observation has a simple explanation in terms of
the NW rupture directivity of the 2016 BS earthquake, leading to stronger seismic waves
propagating towards the DW array.

112

Figure 4.13. Same as Fig. 4.12 for frequency band 0.1-0.4 Hz.

4.5 Discussion and Conclusions
We monitor seismic velocities in the SJFZ in a 2 years period (2015-2016) including the
2016 Mw5.2 Borrego Springs (BS) earthquake using ambient seismic noise recorded from two
linear arrays with an aperture of ~400 m. The two arrays at sites Dry Wash (DW) and Jackass
Flat (JF) are located ~3 km away to the NW and SE of the BS earthquake epicenter. Clear
seasonal velocity changes (δv/v) together with robust co- and post-seismic changes associated
with the BS event are estimated using Coda Wave Interferometry (CWI). Results from two
independent methods estimating δv/v in both time- and frequency domains are compared to
identify spurious changes. Through curve fitting in section 4.4.1, we evaluate the obtained
seasonal variations as a function of lapse time, which is later removed to highlight the co- and
post-seismic δv/v changes. To infer additional information, we perform in section 4.4.2
comparative analyses using different coda wave windows, frequency bands, components, and
array location.
For seasonal variations (section 4.4.1.1), the observed δv/v are likely spurious and governed
by wavefield changes at the low-frequency band 0.1-0.4 Hz that cannot be completely removed
using Singular Value Decomposition based Wiener filter. At higher frequencies between 0.5-4
Hz, the observed annual variations in δv/v beneath the DW and JF fault zone arrays are
comparable to those reported in Hillers et al. (2015b) for the region around the SJFZ. The delay
between the obtained seasonal variations at these high frequencies and the surface annual
113
temperature record is ~40-100 days, suggesting thermoelastic strain is a likely primary
mechanism. However, the delay pattern observed at the JF array for different frequency bands is
different from that at DW, indicating effects of different structural properties or additional
processes.
The effective co-seismic velocity reduction E is estimated to be ~ 6% and ~ 1% at the DW
and JF arrays for the low-frequency band 0.1-0.4 Hz, respectively, and much smaller ~0.2% for
frequencies between 0.5 Hz and 4 Hz (Fig. 4.7). The observed large differences in the obtained
co-seismic velocity reduction values E between the two (stretching and doublet) methods at 0.1 –
0.4 Hz (section 4.4.1.2) can be explained by a co-seismic velocity perturbation that is
significantly larger below the arrays than in the surrounding areas. The co-seismic velocity
changes are expected to be higher within the fault zone structure (Peng and Ben-Zion, 2006).
Therefore, the δv/v estimated from coda waves should decrease with lapse time (section 4.3.1),
which is incompatible with the underlying assumption for both stretching and doublet methods.
Since the stretching method is performed in the time domain, the resulting co-seismic δv/v is
mainly governed by the velocity perturbation associated with the largest scattered arrival in the
selected time window, which is typically the arrival at the earliest lapse time (Fig. 4.3). The
doublet method, on the other hand, is equivalent to a moving window cross-spectral analysis, in
which δv/v is estimated through a linear regression of δt and t obtained in moving windows. In
the linear regression, δt measured at different lapse time t is equally weighted. Therefore, when
δt/t decreases with lapse time, smaller δv/v are expected for the doublet method compared to
those of the stretching method. The observation that larger discrepancy between the two methods
is found at the low frequency band is thus a strong indicator that the degree of localization in
velocity changes is higher at deeper structures. From the comparative analysis in section 4.4.2,
the co-seismic velocity reduction is smaller at later lapse time (e.g. Fig. 4.8 & 4.9a,b).
Considering both arrays are situated within fault zones, the inversely frequency dependent E
implies that the co-seismic medium changes are localized, i.e. much larger within fault zones
below the arrays compared to the surrounding host rock.
Assuming coda waves at lower frequency are more sensitive to changes in deeper structures,
we can use the EQ responses of δv/v at 0.1 – 0.4 Hz to represent co- and post- seismic medium
changes in deep structure, and infer the perturbations in seismic velocities at shallow depth using
the results at 1 – 4 Hz. Since coda wave sensitivity to velocity perturbations extends over larger
114
region at longer lapse time, the rate at which co-seismic velocity reduction decreases with lapse
time observed for a particular frequency band is proportional to, and thus indicative of, the fault
zone width, if body wave conversion and sensitivity is neglected here. Therefore, the higher rate
at 0.1 – 0.4 Hz for DW array (Fig. 4.9a) suggests a narrower damage zone beneath the site DW
than JF, which is consistent with imaging results of the internal fault zone structures below these
two arrays: ~ 200 m wide low velocity zone below JF array (Qiu et al., 2017) and no clear
evidence of fault damage zone at DW (Qiu et al., 2019c). In an effort to correlate the observed
frequency dependent co-seismic δv/v with depth (section 4.4.2.2), we conclude that the co-
seismic velocity reduction at 0.1-0.4 Hz is likely related to co-seismic slips of the BS event
occurred at a depth of ~12 km, whereas the medium changes at shallow depth (captured by δv/v
estimated at frequencies above 0.5 Hz) may be a result of local ground shaking generated by the
BS earthquake.
For the post-seismic process, logarithmic recovery with time is commonly observed for time
scales from seconds to months after material failure in laboratory experiments (e.g. Dietrich and
Kilgore, 1996; Nakatani and Scholz, 2004; Johnson and Jia, 2005) and earthquakes (e.g.
Brenguier et al., 2008b; Wu et al., 2009, 2010). Such post-seismic recovery processes have been
related to healing of co-seismic fractures, motion of fluids and post-seismic stress relaxation (e.g.
Ben-Zion, 2008; Roux and Ben-Zion, 2014; Aben et al., 2017; Pei et al., 2019). Better
understanding of the recovery time for fault zone rocks after an earthquake provides important
constraints on fault strength and may influence various key questions ranging from long-term
coupled evolution of earthquakes and faults (e.g., Ben-Zion et al., 1999; Lyakhovsky et al.,
2001) to estimation of recurrence times of earthquakes (Gratier and Gueydan, 2007).
From section 4.4.1.3, the magnitude of the effective post-seismic recovery time T
c
(~100-400
days) estimated for frequencies above 0.5 Hz is comparable or larger than the analyzed post-
seismic period (~ 6 month). This indicates T
c
is not well constrained in the curve fitting for these
high frequencies and it is supported by the large estimated uncertainty values shown in Figure
4.7d. The significantly smaller T
c
(~10 days) at 0.1 – 0.4 Hz indicates a much faster recovery
compared to those at frequencies above 0.5 Hz. In addition, T
c
values for DW array at
frequencies above 0.5 Hz are significantly larger than that for JF array, suggesting a slower post-
seismic healing process in the medium under DW. The observed frequency dependent effective
post-seismic recovery time (T
c
) is most likely related to normal stress increase with depth, as
115
lower frequency coda waves are more sensitivity to changes in deeper structures. Additional
analyses are needed to interpret how normal stress affect the recovery process in seismic
velocities at different frequencies and will be a topic of a future study.
Similar to the curve fitting of the effective co-seismic velocity reduction E and post-seismic
recovery time T
c
demonstrated in section 4.4.1 using an exponential function (Eq. 6b), Figs.
A4.18 & A4.19 provide such estimations using a logarithmic function (Eq. 5) but the overall
fitting is less good. Note, following equation 5, the effective co-seismic reduction and post-
seismic recovery time are given by E = E
1
and T
c
= e
1−1 e ( )E
1
E
2
−1, respectively. We find an
inconsistency between the predicted and observed δv/v curve during the first 1-2 weeks after the
2016 BS event (Fig. A4.18), which can be explained by the numerous aftershocks (Fig. A4.1).
Waves generated by small magnitude (e.g. Mw < 1.0) aftershocks are difficult to remove in the
preprocessing step (section 4.2) and thus contaminate the computed daily ANC that can bias the
estimated δv/v. In addition, aftershocks also generate medium changes (e.g. fault slip) that further
contribute to the δv/v estimated in the first two weeks after the BS event.
Considering the large uncertainties in T
c
estimations and the effect of aftershocks, we take
the part of δv/v two weeks after the main shock (left dashed lines in Fig. 4.14), and illustrate the
recovery process with a time axis in log scale. Here we further ignore data after the day of 100
(right dashed lines in Fig. 4.14) in the post-seismic period, as the background fluctuation is too
large. Under the assumptions that the co-seismic velocity reduction is fully recovered if given
sufficient time and equation 5, we can fit the δv/v curve in the selected post-seismic period
(between dashed lines in Fig. 4.14) with a linear line (red dashed lines in Fig. 4.14). The
resulting effective co-seismic velocity reduction at T = 0, E
log
, and time to fully recover the co-
seismic velocity reduction, T
log
, are given by E
1
and e
E
1
E
2
−1 , respectively. The best fitting
values of E
log
and T
log
are shown in Figure 4.14. In general, T
log
yields the same magnitude of
values to the T
c
shown in Fig. 4.6. The comparison of T
log
obtained from the DW (Fig. 4.14) and
JF (Fig. A4.20) arrays suggests that the observation of a slower post-seismic healing process at
DW (compared to JF) is robust. The co-seismic velocity reduction E
log
estimated through the
above process yield much larger values compared to E estimated in section 4.4.1.2, particularly
at the low frequency band 0.1-0.4 Hz (E ≈ 6% in Fig. 4.6 and E
log
≈ 37% in Fig. 4.14). This
suggests that the fitting with the log function includes information on earlier parts of the
116
recovery process compared to the exponential function.

Figure 4.14. (a) EQ responses of δv/v curve after the origin time of the 2016 BS earthquake (T =
0) estimated at 0.1-0.4 Hz for DW array shown in log time scale. The dashed lines outline the
portion of data that show a clear logarithmic recovery and is used to estimate the effective co-
seismic velocity reduction E
log
(= E
1
) and post-seismic velocity recovery time T
log
(=e
E
1
E
2
−1 )
following equation 5 (red dashed lines). (b) Same as (a) for 0.5-2 Hz. (c) Same as (a) for 0.75-3
Hz. (d) Same as (a) for 1-4 Hz.

Our lapse time analysis shows results that are compatible to the 3-D wavefield simulations
conducted in Obermann et al. (2016). However, full quantification requires further wavefield
simulations using models that include seismic velocity perturbations in thin vertical layers (i.e.
fault zone) with various widths and receivers placed on top and around the fault zone layer. In
addition to the velocity changes associated with the discussed thermoelastic strain and the 2016
BS event, other external loading mechanisms such solid Earth tides (Hillers et al., 2015a; Mao et
al., 2019), precipitation and fluid content in the crust (Sens-Schönfelder and Wegler, 2006;
Clements and Denolle, 2018), and variable ocean height (Wang et al., 2017), and internal
deformations caused by rapid (Brenguier et al., 2008b; Hobiger et al., 2012) and slow (Inbal et
al., 2017; Hillers et al., 2019) slip on faults can also induce medium perturbations and hence
seismic velocity variations over a wide range of scales. In the study area and time interval, these
mechanisms likely have a minor impact on the observed first order effects associated with the BS
117
event. The signals induced by this earthquake are emphasized by removing seasonal variations
through curve fitting (section 4.4.1). However, there are still non-negligible background
fluctuations with ~ 30 days or shorter periodicity. The curve fitting results (section 4.4.1) and
comparative analyses (section 4.4.2) can potentially be refined by applying corrections of
meteorology data (Wang et al., 2017).
The observations at DW and JF arrays are in general consistent, but we observe
discrepancies in the obtained δv/v curves (section 4.4.2.4). Most of these can be explained by
differences in the local fault zone structures. For instance, fluctuations in the noise source
distribution and surface temperature, which are responsible for the observed seasonal variations
at 0.1-0.4 Hz and above 0.5 Hz, respectively, are approximately the same at both arrays as the
inter-array distance is only ~6 km. Therefore, the differences in amplitude and phase of seasonal
variations obtained for both arrays are due to distinctive local structures underneath, and thus can
be used to calibrate the variations associated with the site. For simplification, we ignore the
differences in the phase, and to minimize these effects, we rescale the obtained δv/v by the
corresponding amplitude of the seasonal variations – C (C – Eq. 4). At the low-frequency 0.1-0.4
Hz and high frequency 1-4 Hz bands, we observe systematically larger E/C at the DW site than
at JF, which correlates with the NW rupture directivity of the 2016 BS earthquake. This suggests
that the observed co-seismic velocity reduction obtained from δv/v using CWI resolve directivity
effects of the earthquake rupture.
Here we resolved for the first time the response to earthquake slip of the near-fault material
using data from two fault zone line arrays. We emphasize that so far all previous studies that
report and analyze seismic velocity changes associated with earthquakes resolved the response of
the medium adjacent to or surrounding the earthquake fault, but not processes on the fault or in
the fault zone or fault damage zone itself. Monitoring processes on active faults or in the
associated damage zones is challenging because of their relatively tabular structure that is
difficult to resolve using the diffusive halo. Monitoring the properties of fault zone trapped noise
and reconstructed trapped modes (Hillers et al., 2014; Hillers et al., 2016) may provide further
possibilities of in situ fault zone probing, as well as attempts to reconstruct fault-penetrating
direct waves.
118
Discussion
I use a set of seismological techniques and several seismic signals to image the geometry
and material properties of the San Jacinto fault (SJF) and surrounding upper crust in southern
California (SC) with spatial resolution ranging from ~10-15 km to a few 100 m wide damage
zones to a local bimaterial interface. In addition to structural properties, the temporal evolution
of fault zones in response to a nearby moderate earthquake with strong directivity is analyzed,
providing important information on rock susceptibility to large earthquakes in the SJF.
These results not only provide improved information on fault zone properties relevant to
seismic hazard in the study region, but also develop the advanced data analysis tools to
automatically process large array data sets. Specifically, velocity structures of the surrounding
upper crust obtained from surface wave tomography will reveal smooth large-scale velocity
contrasts across fault planes, fault geometry, and structural complexity (Chapter 1). The results
of slowness analysis and fault zone head wave (FZHW) studies, which illuminate sharp fault
segments predicting the length of future earthquake ruptures (magnitude) and the subsequent
likely amount of ground shaking. The properties of damage zones obtained from fault zone
trapped wave (FZTW) analysis in combination with the distribution of bi-material interfaces will
provide information on statistically preferred rupture directivity of future earthquakes (Chapters
2 and 3). An earthquake with strong directivity could significantly amplify ground motion in the
direction of propagation than in the opposite direction. Such directivity effect is also observed in
the temporal variations of seismic velocities within the SJFZ related to nearby moderate
earthquakes. Monitoring seismic velocities will improve our understanding of the evolving
susceptibility of the shallow fault zone environment (Chapter 4). Integrating results of fault zone
structures and the neighboring upper crust from each technique for the SJF, this work can
provide 4D information along one of the most active fault systems in the SC plate boundary
region with an unprecedented resolution that will improve the seismic hazard assessment of the
nearby densely populated cities, such as Los Angeles.
Based on the thesis, some potential directions for future studies are:
(1) In Chapter 1, we modified and generalized the Eikonal tomography method to seismic
network with irregular station spacing. In addition to the final isotropic 3-D shear wave velocity
model, we observe differences between imaging results using data from Rayleigh and Love
waves, suggesting the existence of significant apparent radial anisotropy, which may be caused
119
by transverse isotropy (e.g. Moschetti et al., 2010a) or 3D structural effects (e.g. Levshin &
Ratnikova, 1984). The radial and 2-psi azimuthal anisotropy can provide additional information
on crustal properties (e.g. Moschetti et al., 2010b; Lin et al., 2011) in the study region and will be
the subject of a future study. In addition, we also find the resolved Vs structure often suffers
from biases caused by topography and large surface wave amplitude spatial variations, which are
non-negligible in regions such as basins and mountains. To address the topography and
amplitude variation issues, I plan to adopt a waveform modeling technique based on the
GGRTM method (Chen 1996) and Helmholtz tomography technique, respectively. Furthermore,
the surface wave data is poorly fitted within the LA basin, suggesting the necessity of future
studies to improve the velocity model beneath the basin for better hazard evaluation in one of the
most populated cities.
(2) In Chapter 2, we developed a standard procedure, analyses of delay time, FZHW, and FZTW,
to image internal structures of fault zones using linear dense arrays. In the study, we find that the
FZHW identification still requires a considerable amount of human interaction (i.e. visual
examination) to rule out false detections from the current automatic algorithm. To further
improve the identification process and also be able to detect S-type FZHW with a lower signal-
to-noise ratio, a new automatic detection algorithm that utilizes the supervised machine learning
techniques could be very helpful.
(3) For FZTW, we observe clear P-type FZTW (Chapters 2 and 3) and the waveform modeling
will be the subject of a future study. Besides, with the help of the recent development of wireless
low-cost three-component geophones, we are able to characterize the full trapped wavefield
(Allam et al., 2017) and the observed amplitude distributions of the trapped waves are consistent
from event to event, and are thus controlled by fault zone structure rather than source
properties. By using a grid search based method to fit the analytical solution to the observed data,
we can provide independent estimations of the damage zone width, velocity reduction,
asymmetry, and across-fault velocity contrast, and then compare to those obtained from the
FZTW waveform modeling method used in Chapter 2.
(4) Clear fault zone reflected waves (FZRW) are identified in Chapter 3. The obtained
differential time and amplitude ratio between the direct and reflected waves can also provide
information on sharp material fault interface and basin shape. Modeling of the observed FZRW,
120
complementary to the study of FZHW, can also provide sharp images of contrast interfaces and
will be the subject of a future study.
(5) In Chapter 4, we monitored the temporal changes of seismic velocities in the near fault
environment, and found the co-seismic velocity reduction is significantly larger near fault
compared to the intact host rock. Besides, the observed co-seismic velocity reduction correlates
well with the rupture directivity of the target event. However, in this study, we didn’t quantify
the depth distribution of the frequency dependent co-seismic velocity reduction and discuss the
physical process of the observed post-seismic recovery. Simulations using proper structural
models under reasonable stress condition will be necessary to unravel these questions. In
addition, we also see the effect of aftershocks in the ambient noise based monitoring results.
Comparison to results from monitoring methods that use other types of data, such as earthquake
doublets, is useful to better quantify the aftershock effects. A future study of these topics can
address many key questions, such as a better understanding of the processes of co-seismic rock
failure and the effective recovery in fault damaged zones at different depths.
In conclusion, I plan to improve the existing analysis methods and incorporate them with
machine learning into a standard analyzing tool that can be easily applied to other dense array
data sets, at which fault zone structures are not well understood (e.g. China).

121
A1. Chapter 1 supplementary figures

Figure A1.1. (a) One-year stacked cross correlation at ZZ component as a function of
interstation distance. The correlations are stacked for each 1 km distance bin, with black and
white indicating positive and negative amplitudes. Clear Rayleigh wave between the blue
(moveout of 1.5 km/s) and red (moveout of 4 km/s) dashed lines is observed. The maximum
amplitudes are located at a moveout of ~3 km/s. (b) Same as (a) for TT component (i.e. Love
wave). (c) Average group (blue) and phase (orange) velocity dispersion for Rayleigh wave
computed for CVM-S4.26. Error bars indicate a range of one standard deviation of the
corresponding mean values. (d) Same as (c) for Love wave. The surface waves signals are well
bounded by the time window between moveouts of 1.5 km/s and 4 km/s. This is consistent with
the average phase and group velocity dispersions shown in (c), (d) for CVM-S4.26.

122

Figure A1.2. (a) One year stacked ANC at station pair GSC-SDD. The acasual and causal sides
are depicted in blue and black. The dashed lines outline the tapering and signal windows. (b)
Same as Fig. 1.5c for the acausal side. (c) Same as Fig. 1.5c for the causal side. Large differences
in amplitude as a function of center period is observed, which also contribute to the asymmetry
observed in the one year stacked correlogram.

123

Figure A1.3. Depth sensitivity kernels of the surface wave velocity. (a) Rayleigh wave phase
velocity depth sensitivity kernels with respect to Vs calculated at 3s, 7s, and 11s shown in blue,
red, and black, respectively. The dashed curves indicate the results after normalization. (b) Same
as (a) for Love wave phase velocity. (c) Rayleigh wave group (dashed) and phase (solid) velocity
depth sensitivity kernels with respect to Vs (red) and Vp (gray) calculated at period 7s. (d) Love
wave group (dashed) and phase (solid) velocity depth sensitivity kernels with respect to Vs
calculated at period 7s. (e) The Vs profile used to calculate the depth sensitivity kernels. The Vp
profile used in the forward calculation is computed as 1.732 time the corresponding Vs values,
and a homogenous density profile (2.7 g/cm
3
) is assumed.

(e)
124

Figure A1.4. (a) All the phase traveltime measurements (colored triangles) associated with the
central station GSC (red circle) for 7s Rayleigh wave. The black triangles show the locations of
the receivers that do not have a valid measurement when GSC is the virtual source. (b) Travel
time map constructed using the measurements shown in (a) through minimum curvature
interpolation on a 0.05°×0.05° grid. The solid curves represent the phase fronts of the
interpolated travel time map, while the dashed curves depicted the phase front predicted by the
mean phase velocity of the entire map. (c) Curvature map estimated using finite difference. The
measurements shown in (a) are indicated as dots, with red color dots indicating outlier
measurements with large curvature value and green dots as measurements with outlier gradient
values. The white curves depict the phase front of the interpolated traveltime map shown in (b).
(d) The phase speed map derived using eikonal equation. (e) Same as (c) after removing the
outliers. (f) Same as (d) after removing the outliers shown in (c).

125

Figure A1.5. Illustration of removing regions with poor travel time interpolation. (a) Same as
Fig. A1.4d after removing regions within two wavelengths of the GSC and those identified by
the “three- out of four quadrants” selection criterion (Lin et al., 2009) with a searching radius of
50 km. The triangles denote the discrete travel time measurements used to derive the phase
velocity. (b) Differences in interpolated travel time between two different schemes for phase
front interpolation. The color denotes the difference in seconds. The white curves indicate the
contour of 0.5s difference. (c) Estimated station configuration error (section 4.1) map. The white
curves depicted the contour of an error value of 0.025. (d) Final phase velocity map derived
using virtual source GSC for 7s Rayleigh wave (same as Fig. 1.7b). Same as (a) by further
truncating regions with travel time difference shown in (b) > 0.5 and areas with station
configuration error > 0.025.

126

Figure A1.6. 2-psi anisotropy analysis performed on two example grid cells. The black dots
represent the local phase velocity derived from individual eikonal phase velocity maps as a
function of azimuth. In order to fit the systematic pattern, we average the phase velocities within
each 20 degree azimuth bin, and the mean and standard deviation are displayed as the red dot and
error bar. The best fitting 2-psi azimuthal anisotropy is depicted as the red curve with dashed
lines representing the estimated isotropic phase speed. The best fitting amplitude and fast axis of
the 2-psi azimuthal anisotropy are indicated on the top left, while the location of the grid point is
indicated on the bottom left.

127

Figure A1.7. Same as Fig. 1.10 for Love wave. The gray squares in the top left panel show the
location of the grid points explored in Fig. A1.9.

Figure A1.8. Same as Fig. 1.11 for Love wave.

128

Figure A1.9. Illustration of 1D shear wave velocity inversion using Rayleigh wave phase and
group velocity dispersion curves at grid points (gray squares in Fig. A1.7) located at (a) Ventura
basin, (b) Mojave desert, (c) LA basin, (d) Cajon Pass, (e) Trifurcation area in SJFZ, and (f)
Eastern California Shear Zone (ECSZ).

129

Figure A1.10. (a) Same as Fig. 1.12a for Love wave phase and group dispersion curves joint
inversion. The misfit value (0.51) is larger compared to that of the Rayleigh wave (0.37) shown
in Fig. 1.8a. Same initial model extracted from CVM-H15.1 is used, and the Vs profile inverted
from Rayleigh wave data (blue curve; obtained in Fig. 1.12a) is shown in the right panel for
comparison. (b) Same as (a) for joint inversion results using both Rayleigh and Love wave
dispersion curves. The misfit value is much larger (1.2) here indicating the existence of apparent
radial anisotropy as Rayleigh and Love wave dispersion curves cannot be explained with
sufficient accuracy by using just one 1-D Vs profile.

(a)
(b)
130

Figure A1.11. Spatial distribution of the misfit for the joint inversion of Rayleigh wave phase
and group dispersion curves. (a) Misfit calculated for each grid cell using CVM-H15.1. (b) Same
as (a) for the inverted velocity model using CVM-H15.1 as the initial model. (c) Same as (a) for
CVM-S4.26. (d) Same as (b) for using the initial model, CVM-S4.26.

131

Figure A1.12. Map view of Vs maps extracted from CVM-H15.1 at depths of 3 km, 5 km, 7 km,
10 km, 15 km, and 20 km. The color palette used in each panel is the same as the map view of
the final Vs at the corresponding depth.

132

Figure A1.13. Same as Fig. 1.17 for cross sections of Vs in terms of relative changes (final -
initial) and the initial model, CVM-H15.1 (right column).

133

Figure A1.14. Same as Fig. 1.14 for joint inversion of Love wave phase and group dispersion
curves.

(a) (b)
(c) (d)
134

Figure A1.15. Same as Fig. A1.11 for joint inversion of Love wave phase and group dispersion
curves.

135

Figure A1.16. Same as Fig. 1.10 for Love wave joint inversion results.

Figure A1.17. Same as Fig. 1.11 for Love wave joint inversion results.

136

Figure A1.18. Same as Fig. 1.17 for Love wave joint inversion results.

137

Figure A1.19. Same as Fig. A1.13 for Love wave joint inversion results.

138
A2. Chapter 2 supplementary figures

Figure A2.1. (a) Same as Fig 2.4(a). The red circle denotes event #1 analyzed in the Supporting
Information; the event information is shown in the title. (b) Spectrogram for the analyzed event
(same as Fig. 2.4c). (c) Same as Fig. 2.5(b) for the analyzed event. (d) Same as Fig. 2.7(a) for the
analysed event. (e) Same as Fig. 2.7(b) for the analyzed event.

Event #1 P phase 0.01 -0.2 Hz
Time to origin (s)
Time to origin (s)
Northing from JF00 (km)
TR04
TR03
TR02
TR01
JF
JFS4
JFS3
JFS2
JFS1
JF00
JFN1
JFN2
JFN3
JFN4
(a) (b)
(c)
(d)
(e)
Delay time (s)
139

Figure A2.2. Same as Fig. A2.1 in the Supporting Information for event # 2 analyzed in the
Supporting Information.

JFS4
JFS3
JFS2
JFS1
JF00
JFN1
JFN2
JFN3
JFN4
(a) (b)
(c) (d)
Event #2 P phase 0.2-0.5 Hz
Time to origin (s) Northing from JF00 (km)
Time to origin (s)
Delay time (s)
140

Figure A2.3. (a) Location map and depth section of 3493 local earth- quakes within the red box
shown in Fig. 2.1(a), used in the local P-wave DTA to estimate relative slowness. (b) Same as
(a) with events SW and NE of the fault denoted by red and green circles, respectively. (c) Same
as (a) with events SE and NW of the JF array denoted by red and green circles, respectively. (d)
Same as (a) with green (I), purple (II), red (III) and blue (IV) circles representing events west,
north, east and south of the JF array, respectively.

141

Figure A2.4. (a) Mean values of the relative slowness computed for the 1434 events SW (red)
and 2059 events NE (green) of the fault at the JF array. Same as Fig. 2.10(c); the black curve
denotes the relative slowness for all 3493 events. Error bars indicate a range of two standard
deviations about each respective mean value. (b) Same as (a) for 1132 events SE (red) and 2361
events NW (green) of the fault. (c) Same as (a) for 637, 1724, 797 and 335 events in quadrants I
(green), II (purple), III (blue) and IV (red). Despite the variations between relative slowness
calculated for events in different subsets, the shapes of the mean values remain similar.

142

Figure A2.5. A model for a vertical low-velocity fault zone layer in a half-space used in the
inversion of fault zone trapped waves (Ben-Zion and Aki 1990; Ben-Zion 1998). The source is
an SH line dislocation with coordinates (x
S
, z
S
). Width, quality factor and shear wave velocity are
marked by W, Q and β, respectively. Subscripts H and FZ denote the half-space and fault zone
layer.

143

Figure A2.6. Pre-processing steps for inversion of trapped waves. (a) Velocity seismograms at
the JF array generated by event TW2 (see Fig. 2.15 for location) are corrected for instrument
response and rotated to the fault-parallel direction. (b) Bandpass filtering at 2–20 Hz. (c)
Integrating velocity to displacement seismograms. (d) Convolving waveforms shown in (c) with
1/t
1/2
. Red triangles denote the automatic S picks and stations in bold indicate waveforms with
clear candidate trapped waves following the S arrival.

144

Figure A2.7. Inversion results of displacement seismograms at the JF array generated by
candidate event TW3. Instead of the simpler model shown in Fig. A2.5 in the Supporting
Information, different shear wave velocities are allowed for the quarter spaces left (QS SW) and
right (QS NE) to the vertical fault zone layer. (a) Comparison between the observed (black) and
synthetic (green) seismograms. (b) Fitness values (green dots) calculated for different FZ
parameters tested in the inversion. The parameters
145

Figure A2.8. Same as Fig. 2.17 for event TW2.

146

A3. Chapter 3 supplementary figures

Figure A3.1. Delay time patterns of the example teleseismic event shown in Fig. 3.4 for (a) 0.2 –
0.5 Hz, (b) 0.5 – 1 Hz, and (c) 1 – 2 Hz. As teleseismic P wave is much sensitive to shallower
structure at higher frequency and thus has slower velocity, here we use 2 km/s in (a), 1 km/s in
(b), and 0.5 km/s in (c) for topography correction.

(a)
(b)
(c)
147

Figure A3.2. S-type FZTW candidate identified from the automatic detection algorithm of Ross
et al. 2015.

148

Figure A3.3. P (vertical; top panel) and S (fault parallel – bottom panel; fault normal – bottom
panel) waveforms of the event shown in Fig. 3.10. Clear fault zone reflected waves are observed
after both P and S waveforms at the vertical component.

149
A4. Chapter 4 supplementary figures

Figure A4.1. Magnitude-time distribution (black dots) of the seismicity that occurred within the
trifurcation area of the SJFZ shown in Figure 4.1 for the period 2 days before and 15 days after
the occurrence time of the 2016 BS earthquake. The beach ball indicates the focal mechanism of
the main shock. Red stars denote aftershocks with Mw > 3.0. The daily frequency of seismicity
is shown in the blue.

Figure A4.2. Same as Fig. 4.2 for daily cross correlations before removing outlier correlations
and Wiener filtering.

150

Figure A4.3. Same as Fig. 4.3a for frequency band 0.25-1 Hz. We do not include this frequency
band due to insufficient SNR.

151

Figure A4.4. (a). Synthetic co- and post-seismic δv/v curves calculated using equation 5 for a
selection of E
1
/E
2
values. (b) Synthetic co- and post-seismic velocity change predicted using
equation 6a (in black) and 6b (in red). The difference between these two curves is shown in blue.
(c) Same as (a) using equation 6b for the parameter λ. (d) Same as (c) for the parameter T
s
.

152

Figure A4.5. Same as Fig. 4.5 but for 0.1-0.4 Hz using a coda wave window of 8-40s. Different
from Fig. 4.5, large periodic fluctuations (~ 6% peak to peak variation & ~40-day periodicity) in
δv/v is observed between May and October in 2015.

153

Figure A4.6. Same as Fig. 4.5 but for 0.5-2 Hz and using a coda wave window of 8-25s.

154

Figure A4.7. Same as Fig. 4.5 but for 0.75-3 Hz and using a coda wave window of 6-23s.

155

Figure A4.8. Same as Fig. 4.6 for results using doublet method. Although the amplitudes of
seasonal and co-seismic components are smaller compared to those using stretching method, the
temporal patterns are similar regardless of the method used.

156

Figure A4.9. Same as Fig. 4.6 for JF array.

157

Figure A4.10. Same as Fig. A4.9 for results using doublet method.

158

Figure A4.11. (a) Same as Fig. 4.8a before removing the long term variations including a linear
trend, annual and semi-annual variations. (b) Same as (a) but for JF array. (c), (d) Same as (a),
(b) for frequency band 0.5-2 Hz. A very high consistency in seasonal variations is observed for
all lapse time windows. Same observations are seen in results from doublet method (not shown
here).

159

Figure A4.12. Same as Fig. A4.11 for frequency bands (a), (b) 0.75-3 Hz and (c), (d) 1-4 Hz.
Almost the same seasonal variation is observed for all lapse time windows. Same observations
are seen in results from doublet method (not shown here).

160

Figure A4.13. Same as Fig. 4.11 using the doublet method. Here only results for components ZZ
(vertical-vertical in black), RR (fault parallel – fault parallel in blue), and TT (fault normal –
fault normal in red) are shown.

161

Figure A4.14. Same as Fig. 4.11 for the JF array.

162

Figure A4.15. Same as Fig. A4.14 for doublet method.

163

Figure A4.16. Same as Fig. 4.12 for the frequency band 0.5-2 Hz.

Figure A4.17. Same as Fig. 4.12 for frequency band 0.75-3 Hz.

164

Figure A4.18. Same as Fig. 4.6 using equation 5. Here E = E
1
and T
c
= e
1−1 e ( )E
1
E
2
−1.

165

Figure A4.19. Same Fig. A4.18 for JF array.

166

Figure A4.20. Same Fig. 4.14 for JF array.

167
References
Aben, F., Doan, M., Gratier, J., and Renard, F., 2017, Experimental postseismic recovery of
fractured rocks assisted by calcite sealing: Geophysical Research Letters, v. 44, no. 14, p.
7228-7238, doi: 10.1002/2017gl073965.
Aki, K., and Richards, P., 2009, Quantitative Seismology: University Science Books, Mill
Valley.
Allam, 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, no. 2, p. 1181-1196, doi: 10.1111/j.1365-246x.2012.05544.x.
Allam, A., Ben-Zion, Y., Kurzon, I., and Vernon, F., 2014a, Seismic velocity structure in the Hot
Springs and Trifurcation areas of the San Jacinto fault zone, California, from double-
difference tomography: Geophysical Journal International, v. 198, no. 2, p. 978-999, doi:
10.1093/gji/ggu176.
Allam, A., Ben-Zion, Y., and Peng, Z., 2014b, Seismic Imaging of a Bimaterial Interface Along
the Hayward Fault, CA, with Fault Zone Head Waves and Direct P Arrivals: Pure and
Applied Geophysics, v. 171, no. 11, p. 2993-3011, doi: 10.1007/s00024-014-0784-0.
Allam, A. A., Qiu, H., Lin, F. C., and Ben-Zion, Y., 2017, Damage Zone Resonance:
Observation and Analytical Fitting of Fault Zone Trapped Wave Normal Modes Using
Large-N Near-fault Seismic Arrays: Abstract presented at 2017 Fall Meeting, AGU, New
Orleans, Louisiana, 11-15 December.
Ampuero, J., 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, no. 2, p. 674-692, doi: 10.1111/j.1365-246x.2008.03736.x.
Barak, S., Klemperer, S., and Lawrence, J., 2015, San Andreas Fault dip, Peninsular Ranges
mafic lower crust and partial melt in the Salton Trough, Southern California, from
ambient-noise tomography: Geochemistry, Geophysics, Geosystems, v. 16, no. 11, p.
3946-3972, doi: 10.1002/2015gc005970.
Barmin, M., Ritzwoller, M., and Levshin, A., 2001, A Fast and Reliable Method for Surface
Wave Tomography: Pure and Applied Geophysics, v. 158, no. 8, p. 1351-1375, doi:
10.1007/pl00001225.
168
Bensen, G., Ritzwoller, M., Barmin, M., Levshin, A., Lin, F., Moschetti, M., Shapiro, N., and
Yang, Y., 2007, Processing seismic ambient noise data to obtain reliable broad-band
surface wave dispersion measurements: Geophysical Journal International, v. 169, no. 3,
p. 1239-1260, doi: 10.1111/j.1365-246x.2007.03374.x.
Ben-Zion, Y., and Leary, P., 1986, Thermoelastic strain in a half space covered by
unconsolidated material: Bulletin of the Seismological Society of America, v. 80, no. 4,
p. 971-994.
Ben-Zion, Y., 1989, The response of two joined quarter spaces to SH line sources located at the
material discontinuity interface: Geophysical Journal International, v. 98, no. 2, p. 213-
222, doi: 10.1111/j.1365-246x.1989.tb03346.x.
Ben-Zion, Y., 1990, The response of two half spaces to point dislocations at the material
interface: Geophysical Journal International, v. 101, no. 3, p. 507-528, doi:
10.1111/j.1365-246x.1990.tb05567.x.
Ben-Zion, Y., and Aki, K., 1990, Seismic radiation from an SH line source in a laterally
heterogeneous planar fault zone: Bulletin of the Seismological Society of America, v. 80,
no. 4, p. 971-994.
Ben-Zion, Y., and Malin, P., 1991, San Andreas Fault Zone Head Waves Near Parkfield,
California: Science, v. 251, no. 5001, p. 1592-1594, doi: 10.1126/science.251.5001.1592.
Ben-Zion, Y., 1998, Properties of seismic fault zone waves and their utility for imaging low-
velocity structures: Journal of Geophysical Research: Solid Earth, v. 103, no. B6, p.
12567-12585, doi: 10.1029/98jb00768.
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, no. 4, p.
1085-1094.
Ben-Zion, Y., Dahmen, K., Lyakhovsky, V., Ertas, D., and Agnon, A., 1999, Self-driven mode
switching of earthquake activity on a fault system: Earth and Planetary Science Letters, v.
172, no. 1-2, p. 11-21, doi: 10.1016/s0012-821x(99)00187-9.
Ben-Zion, Y., and Huang, Y., 2002, Dynamic rupture on an interface between a compliant fault
zone layer and a stiffer surrounding solid: Journal of Geophysical Research, v. 107, no.
B2, doi: 10.1029/2001jb000254.
169
Ben-Zion, Y., and Sammis, C., 2003, Characterization of Fault Zones: Pure and Applied
Geophysics, v. 160, no. 3, p. 677-715, doi: 10.1007/pl00012554.
Ben-Zion, Y., Peng, Z., Okaya, D., Seeber, L., Armbruster, J.G., Ozer, N., Michael, A.J., Baris,
S., and Aktar, M., 2003, A shallow fault-zone structure illuminated by trapped waves in
the Karadere–Duzce branch of the North Anatolian Fault, western Turkey: Geophysical
Journal International, v. 152, no. 3, p. 699-717.
Ben-Zion, Y., and Shi, Z., 2005, Dynamic rupture on a material interface with spontaneous
generation of plastic strain in the bulk: Earth and Planetary Science Letters, v. 236, no. 1-
2, p. 486-496, doi: 10.1016/j.epsl.2005.03.025.
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, no. 4, doi: 10.1029/2008rg000260.
Ben-Zion, Y., and Zaliapin, I., 2019, Spatial variations of rock damage production by
earthquakes in southern California: Earth and Planetary Science Letters, v. 512, p. 184-
193, doi: 10.1016/j.epsl.2019.02.006.
Berg, E., Lin, F., Allam, A., Qiu, H., Shen, W., and Ben-Zion, Y., 2018, Tomography of
Southern California Via Bayesian Joint Inversion of Rayleigh Wave Ellipticity and Phase
Velocity From Ambient Noise Cross-Correlations: Journal of Geophysical Research:
Solid Earth, v. 123, no. 11, p. 9933-9949, doi: 10.1029/2018jb016269.
Berger, J., 1975, A note on thermoelastic strains and tilts: Journal of Geophysical Research, v.
80, no. 2, p. 274-277, doi: 10.1029/jb080i002p00274.
Boué, P., Roux, P., Campillo, M., and de Cacqueray, B., 2013, Double beamforming processing
in a seismic prospecting context: GEOPHYSICS, v. 78, no. 3, p. V101-V108, doi:
10.1190/geo2012-0364.1.
Brenguier, F., Shapiro, N., Campillo, M., Ferrazzini, V., Duputel, Z., Coutant, O., and
Nercessian, A., 2008a, Towards forecasting volcanic eruptions using seismic noise:
Nature Geoscience, v. 1, no. 2, p. 126-130, doi: 10.1038/ngeo104.
Brenguier, F., Campillo, M., Hadziioannou, C., Shapiro, N., Nadeau, R., and Larose, E., 2008b,
Postseismic Relaxation Along the San Andreas Fault at Parkfield from Continuous
Seismological Observations: Science, v. 321, no. 5895, p. 1478-1481, doi:
10.1126/science.1160943.
170
Brenguier, F., Campillo, M., Takeda, T., Aoki, Y., Shapiro, N., Briand, X., Emoto, K., and
Miyake, H., 2014, Mapping pressurized volcanic fluids from induced crustal seismic
velocity drops: Science, v. 345, no. 6192, p. 80-82, doi: 10.1126/science.1254073.
Brietzke, G.B., and Ben-Zion, Y., 2006, Examining tendencies of in-plane rupture to migrate to
material interfaces: Geophysical Journal International, v. 167, no. 2, p. 807-819.
Brietzke, G., Cochard, A., and Igel, H., 2009, Importance of bimaterial interfaces for earthquake
dynamics and strong ground motion: Geophysical Journal International, v. 178, no. 2, p.
921-938, doi: 10.1111/j.1365-246x.2009.04209.x.
Bulut, F., Ben-Zion, Y., and Bohnhoff, M., 2012, Evidence for a bimaterial interface along the
Mudurnu segment of the North Anatolian Fault Zone from polarization analysis of P
waves: Earth and Planetary Science Letters, v. 327-328, p. 17-22, doi:
10.1016/j.epsl.2012.02.001.
Calderoni, G., Di Giovambattista, R., Vannoli, P., Pucillo, S., and Rovelli, A., 2012, Fault-
trapped waves depict continuity of the fault system responsible for the 6 April 2009 MW
6.3 L'Aquila earthquake, central Italy: Earth and Planetary Science Letters, v. 323-324, p.
1-8, doi: 10.1016/j.epsl.2012.01.003.
Campillo, M., and Paul, A., 2003, Long-Range Correlations in the Diffues Seismic Coda:
Science, v. 299 no. 5606, p. 547-549, doi:10.1126/science.1078551.
Chen, X., 1996, Seismogram synthesis for multi-layered media with irregular interfaces by
global generalized reflection/transmission matrices method. III. Theory of 2D P-SV
case: Bulletin of the Seismological Society of America, v. 86, no. 2, p. 389-405.
Cheng, Y., Ross, Z., and Ben-Zion, Y., 2018, Diverse Volumetric Faulting Patterns in the San
Jacinto Fault Zone: Journal of Geophysical Research: Solid Earth, v. 123, no. 6, p. 5068-
5081, doi: 10.1029/2017jb015408.
Clarke, D., Zaccarelli, L., Shapiro, N., and Brenguier, F., 2011, Assessment of resolution and
accuracy of the Moving Window Cross Spectral technique for monitoring crustal
temporal variations using ambient seismic noise: Geophysical Journal International, v.
186, no. 2, p. 867-882, doi: 10.1111/j.1365-246x.2011.05074.x.
Clements, T., and Denolle, M., 2018, Tracking Groundwater Levels Using the Ambient Seismic
Field: Geophysical Research Letters, v. 45, no. 13, p. 6459-6465, doi:
10.1029/2018gl077706.
171
Cochran, E., Li, Y., Shearer, P., Barbot, S., Fialko, Y., and Vidale, J., 2009, Seismic and
geodetic evidence for extensive, long-lived fault damage zones: Geology, v. 37, no. 4, p.
315-318, doi: 10.1130/g25306a.1.
Colombi, A., Chaput, J., Brenguier, F., Hillers, G., Roux, P., and Campillo, M., 2014, On the
temporal stability of the coda of ambient noise correlations: Comptes Rendus
Geoscience, v. 346, no. 11-12, p. 307-316.
Crotwell, H.P., Owens, T.J., and Ritsema, J., 1999, The TauP Toolkit: Flexible seismic travel-
time and ray-path utilities: Seismological Research Letters, v. 70, no. 2, p. 154-160.
Dahlen, F. A., Huang, S.-H., and Nolet, G., 2000, Fr'echet kernels for finite-frequency travel
times - I. Theory: Geophysical Journal International, v. 141, p. 157-174.
Dieterich, J. H., and Kilgore, B. D., 1996, Imaging surface contacts: Power law contact
distributions and contact stresses in quarts, calcite, glass, and acrylic plastic:
Tectonophysics, v. 256, p. 219–239, doi:10.1016/0040-1951(95)00165-4.
Dor, O., Rockwell, T.K., and Ben-Zion, Y., 2006, Geological observations of damage
asymmetry in the structure of the San Jacinto, San Andreas and Punchbowl faults in
Southern California: a possible indicator for preferred rupture propagation direction: Pure
and Applied Geophysics, v. 163, no. 2-3, p. 301-349.
Dor, O., Yildirim, C., Rockwell, T.K., Ben-Zion, Y., Emre, O., Sisk, M., and Duman, T.Y.,
2008, Geological and geomorphologic asymmetry across the rupture zones of the 1943
and 1944 earthquakes on the North Anatolian Fault: possible signals for preferred
earthquake propagation direction: Geophysical Journal International, v. 173, no. 2, p.
483-504.
Ellsworth, W.L., and Malin, P.E., 2011, Deep rock damage in the San Andreas Fault revealed by
P-and S-type fault-zone-guided waves: Geological Society, London, Special Publications,
v. 359, no. 1, p. 39-53.
Fang, H., Zhang, H., Yao, H., Allam, A., Zigone, D., Ben-Zion, Y., Thurber, C., and
van der Hilst, R., 2016, A new algorithm for three-dimensional joint inversion of body
wave and surface wave data and its application to the Southern California plate boundary
region: Journal of Geophysical Research: Solid Earth, v. 121, no. 5, p. 3557-3569, doi:
10.1002/2015jb012702.
172
Fay, N., and Humphreys, E., 2005, Fault slip rates, effects of elastic heterogeneity on geodetic
data, and the strength of the lower crust in the Salton Trough region, southern California:
Journal of Geophysical Research: Solid Earth, v. 110, no. B9, doi:
10.1029/2004jb003548.
Fohrmann, M., Igel, H., Jahnke, G., and Ben-Zion, Y., 2004, Guided Waves from Sources
Outside Faults: An Indication for Shallow Fault Zone Structure?: Pure and Applied
Geophysics, v. 161, no. 11-12, doi: 10.1007/s00024-004-2553-y.
Froment, B., Campillo, M., Chen, J. H., and Liu, Q. Y., 2013, Deformation at depth associated
with the 12 May 2008 MW 7.9 Wenchuan earthquake from seismic ambient noise
monitoring: Geophysical Research Letters, v. 40, no. 1, p. 78-82,
doi:10.1029/2012gl053995.
Gouédard, P., Yao, H., Ernst, F., and van der Hilst, R. D., 2012, Surface wave eikonal
tomography in heterogeneous media using exploration data: Geophysical Journal
International, v. 191, no. 2, p. 781-788, doi:10.1111/j.1365-246X.2012.05652.x.
Gratier, J.P., and Gueydan, F., 2007, Effect of fracturing and fluid–rock interaction on seismic
cycles: Tectonic Faults: Agents of Change on a Dynamic Earth, v. 95, p. 319e356.
Grêt, A., Snieder, R., and Scales J., 2006, Time-lapse monitoring of rock properties with coda
wave interferometry, Journal of Geophysical Research, v. 111, no. B3,
doi:10.1029/2004jb003354.
Hadziioannou, C., Larose, E., Coutant, O., Roux, P., and Campillo, M., 2009, Stability of
monitoring weak changes in multiply scattering media with ambient noise correlation:
Laboratory experiments: The Journal of the Acoustical Society of America, v. 125, no. 6,
p. 3688-3695, doi: 10.1121/1.3125345.
Hadziioannou, C., Larose, E., Baig, A., Roux, P., and Campillo, M., 2011, Improving temporal
resolution in ambient noise monitoring of seismic wave speed: Journal of Geophysical
Research, v. 116, no. B7, doi: 10.1029/2011jb008200.
Hauksson, E., and Shearer, P.M., 2006, Attenuation models (QP and QS) in three dimensions of
the southern California crust: Inferred fluid saturation at seismogenic depths: Journal of
Geophysical Research: Solid Earth, v. 111, no. B5.
173
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, no. 5, p. 2239-2244, doi:10.1785/0120120010.
Hauksson, E., and Meier, M., 2018, Applying Depth Distribution of Seismicity to Determine
Thermo-Mechanical Properties of the Seismogenic Crust in Southern California:
Comparing Lithotectonic Blocks: Pure and Applied Geophysics, v. 176, no. 3, p. 1061-
1081, doi: 10.1007/s00024-018-1981-z.
Herrmann, R. B., 2013, Computer Programs in Seismology: An Evolving Tool for Instruction
and Research: Seismological Research Letters, v. 84, no. 6, p. 1081-1088, doi:10.1785/
0220110096.
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, no. 7, p. 2182-2197, doi:10.1002/ggge.20140.
Hillers, G., Campillo, M., Ben-Zion, Y., and Roux, P., 2014, Seismic fault zone trapped noise:
Journal of Geophysical Research: Solid Earth, v. 119, no. 7, p. 5786-5799, doi:
10.1002/2014jb011217.
Hillers, G., Retailleau, L., Campillo, M., Inbal, A., Ampuero, J., and Nishimura, T., 2015a, In
situ observations of velocity changes in response to tidal deformation from analysis of the
high-frequency ambient wavefield: Journal of Geophysical Research: Solid Earth, v. 120,
no. 1, p. 210-225, doi: 10.1002/2014jb011318.
Hillers, G., Ben-Zion, Y., Campillo, M., and Zigone, D., 2015b, Seasonal variations of seismic
velocities in the San Jacinto fault area observed with ambient seismic noise: Geophysical
Journal International, v. 202, no. 2, p. 920-932, doi: 10.1093/gji/ggv151.
Hillers, G., and Campillo, M., 2016, Fault zone reverberations from cross-correlations of
earthquake waveforms and seismic noise: Geophysical Journal International, v. 204, no.
3, p. 1503-1517, doi: 10.1093/gji/ggv515.
Hillers, G., Campillo, M., Brenguier, F., Moreau, L., Agnew, D. C., and Ben-Zion, Y.,
2019, Seismic velocity change patterns along the San Jacinto fault zone following the
2010 M7.2 El Mayor-Cucapah and M5.4 Collins Valley earthquakes: Journal of
Geophysical Research (in review).
174
Hobiger, M., Wegler, U., Shiomi, K., and Nakahara, H., 2012, Coseismic and postseismic elastic
wave velocity variations caused by the 2008 Iwate-Miyagi Nairiku earthquake, Japan:
Journal of Geophysical Research: Solid Earth, v. 117, no. B9, doi:
10.1029/2012jb009402.
Hobiger, M., Wegler, U., Shiomi, K., and Nakahara, H., 2014, Single-station cross-correlation
analysis of ambient seismic noise: application to stations in the surroundings of the 2008
Iwate-Miyagi Nairiku earthquake: Geophysical Journal International, v. 198, no. 1, p. 90-
109, doi: 10.1093/gji/ggu115.
Hobiger, M., Wegler, U., Shiomi, K., and Nakahara, H., 2016, Coseismic and post-seismic
velocity changes detected by Passive Image Interferometry: comparison of one great and
five strong earthquakes in Japan: Geophysical Journal International, v. 205, no. 2, p.
1053-1073, doi: 10.1093/gji/ggw066.
Hough, S.E., Ben-Zion, Y., and Leary, P., 1994, Fault-zone waves observed at the southern
Joshua Tree earthquake rupture zone: Bulletin of the Seismological Society of America,
v. 84, no. 3, p. 761-767.
Huang, Y., Ampuero, J., and Helmberger, D., 2014, Earthquake ruptures modulated by waves in
damaged fault zones: Journal of Geophysical Research: Solid Earth, v. 119, no. 4, p.
3133-3154, doi: 10.1002/2013jb010724.
Igel, H., Ben-Zion, Y., and Leary, P.C., 1997, Simulation of SH-and P-SV-wave propagation in
fault zones: Geophysical Journal International, v. 128, no. 3, p. 533-546.
Igel, H., Jahnke, G., and Ben-Zion, Y., 2002, Numerical simulation of fault zone guided waves:
accuracy and 3-D effects: Pure and Applied Geophysics, v. 159, p. 2067-2083.
Inbal, A., Ampuero, J. P., and Avouac, J. P., 2017, Locally and remotely triggered aseismic slip
on the central San Jacinto Fault near Anza, CA, from joint inversion of seismicity and
strainmeter data: Journal of Geophysical Research: Solid Earth, v. 122, no. 4, p. 3033-
3061, doi:10.1002/2016jb013499.
Inbal, A., Cristea‐Platon, T., Ampuero, J., Hillers, G., Agnew, D., and Hough, S., 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.
175
Jahnke, G., Igel, H., and Ben-Zion, Y., 2002, Three-dimensional calculations of fault-zone-
guided waves in various irregular structures: Geophysical Journal International, v. 151,
no. 2, p. 416-426.
Johnson, C., Fu, Y., and Bürgmann, R., 2017, Stress Models of the Annual Hydrospheric,
Atmospheric, Thermal, and Tidal Loading Cycles on California Faults: Perturbation of
Background Stress and Changes in Seismicity: Journal of Geophysical Research: Solid
Earth, v. 122, no. 12, p. 10,605-10,625, doi: 10.1002/2017jb014778.
Johnson, P.A., and Jia, X., 2005, Nonlinear dynamics, granular media and dynamic earthquake
triggering: Nature, v. 473, p. 871–874.
Kedar, S., and Webb, F. H., 2005, The Ocean's Seismic Hum: Science, v. 307,
doi:10.1126/science.1108380.
Kennett, B.L.N., and Engdahl, E.R., 1991, Travel times for global earthquake location and phase
identification: Geophysical Journal International, v. 105, no. 2, p. 429-465.
Kissling, E., 1988, Geotomography with local earthquake data: Reviews of Geophysics, v. 26,
no. 4, p. 659-698.
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, no. 11, p. 3045-3081.
Lee, E.-J., Chen, P., Jordan, T. H., Maechling, P. B., Denolle, M. A. M., 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:
Solid Earth, v. 119, no. 8, p. 6421-6451, doi:10.1002/2014jb011346.
Lévěque, J. J., Rivera, L., and Wittlinger, G., 1993, On the use of the checker-board test to assess
the resolution of tomographic inversions: Geophysical Journal International, v. 115, no.
1, p. 313-318.
Levshin, A., and Ratnikova, L., 1984, Apparent ansiotropy in inhomogeneous media:
Geophysical Journal International, v. 76, no. 1, p. 65-69, doi: 10.1111/j.1365-
246x.1984.tb05022.x.
Lewis, M.A., Peng, Z., Ben-Zion, Y., and Vernon, F.L., 2005, Shallow seismic trapping structure
in the San Jacinto fault zone near Anza, California: Geophysical Journal International, v.
162 no. 3, p. 867-881.
176
Lewis, M.A., and Ben-Zion, Y., 2010, Diversity of fault zone damage and trapping structures in
the Parkfield section of the San Andreas Fault from comprehensive analysis of near fault
seismograms: Geophysical Journal International, v. 183, no. 3, p. 1579-1595.
Li, Y., Leary, P., Aki, K., and Malin, P., 1990, Seismic Trapped Modes in the Oroville and San
Andreas Fault Zones: Science, v. 249, no. 4970, p. 763-766, doi:
10.1126/science.249.4970.763.
Li, Y., Wu, Q., Pan, J., and Sun, L., 2012, S-wave velocity structure of northeastern China from
joint inversion of Rayleigh wave phase and group velocities: Geophysical Journal
International, v. 190, no. 1, p. 105-115.
Li, Y.G., and Vernon, F.L., 2001, Characterization of the San Jacinto fault zone near Anza,
California, by fault zone trapped waves: Journal of Geophysical Research: Solid Earth, v.
106, no. B12, p. 30671-30688.
Li, Z., Peng, Z., Ben-Zion, Y., and Vernon, F. L., 2015, Spatial variations of Shear-wave
Anisotropy Near the San Jacinto Fault Zone in Southern California: Journal of
Geophysical Research: Solid Earth, v. 120, p. 8334-8347, doi: 10.1002/2015JB012483.
Lin, F. C., Ritzwoller, M. H., Townend, J., Bannister, S., and Savage, M. K., 2007, Ambient
noise Rayleigh wave tomography of New Zealand: Geophysical Journal International, v.
170, no. 2, p. 649-666.
Lin, F.-C., Moschetti, M. P., and Ritzwoller, M. H., 2008, Surface wave tomography of the
western United States from ambient seismic noise: Rayleigh and Love wave phase
velocity maps: Geophysical Journal International, v. 173, no. 1, p. 281-298,
doi:10.1111/j.1365-246X.2008.03720.x.
Lin, F.-C., Ritzwoller, M. H., and Snieder, R., 2009, Eikonal tomography: surface wave
tomography by phase front tracking across a regional broad-band seismic array:
Geophysical Journal International, v. 177, no. 3, p. 1091-1110, doi:10.1111/j.1365-
246X.2009.04105.x.
Lin, F.-C., and Ritzwoller, M. H., 2011, Helmholtz surface wave tomography for isotropic and
azimuthally anisotropic structure: Geophysical Journal International, v. 186, no. 3, p.
1104-1120, doi:10.1111/j.1365-246X.2011.05070.x.
177
Lin, F.-C., Li, D., Clayton, R. W., and Hollis, D., 2013, High-resolution 3D shallow crustal
structure in Long Beach, California: Application of ambient noise tomography on a dense
seismic array: Geophysics, v. 78, no. 4, p. Q45-Q56, doi:10.1190/GEO2012-0453.1.
Lindsey, E.O., and Fialko, Y., 2013, Geodetic slip rates in the southern San Andreas Fault
system: Effects of elastic heterogeneity and fault geometry: Journal of Geophysical
Research: Solid Earth, v. 118, no. 2, p. 689-697.
Lyakhovsky, V., Ben-Zion, Y., and Agnon, A., 2001, Earthquake cycle, fault zones, and
seismicity patterns in a rheologically layered lithosphere: Journal of Geophysical
Research: Solid Earth, v. 106, no. B3, p. 4103-4120, doi: 10.1029/2000jb900218.
Lyakhovsky, V., and Ben-Zion, Y., 2009, Evolving geometrical and material properties of fault
zones in a damage rheology model: Geochemistry, Geophysics, Geosystems, v. 10, no.
11, doi:10.1029/2009gc002543.
Mao, S., Campillo, M., Hilst, R., Brenguier, F., Stehly, L., and Hillers, G., 2019, High Temporal
Resolution Monitoring of Small Variations in Crustal Strain by Dense Seismic Arrays:
Geophysical Research Letters, v. 46, no. 1, p. 128-137, doi: 10.1029/2018gl079944.
Marliyani, G., Rockwell, T., Onderdonk, N., and McGill, S., 2013, Straightening of the
Northern San Jacinto Fault, California, as Seen in the Fault-Structure Evolution of the
San Jacinto Valley Stepover: Bulletin of the Seismological Society of America, v. 103,
no. 3, p. 2047-2061, doi: 10.1785/0120120232.
McGuire, J., and Ben-Zion, Y., 2005, High-resolution imaging of the Bear Valley section of the
San Andreas Fault at seismogenic depths with fault-zone head waves and relocated
seismicity: Geophysical Journal International, v. 163, no. 1, p. 152-164.
Mikesell, T., Malcolm, A., Yang, D., and Haney, M., 2015, A comparison of methods to estimate
seismic phase delays: numerical examples for coda wave interferometry: Geophysical
Journal International, v. 202, no. 1, p. 347-360, doi: 10.1093/gji/ggv138.
Mizuno, T., and Nishigami, K., 2006, Deep structure of the Nojima Fault, southwest Japan,
estimated from borehole observations of fault-zone trapped waves: Tectonophysics, v.
417, no. 3, p. 231-247.
Moreau, L., Stehly, L., Boué, P., Lu, Y., Larose, E., and Campillo, M., 2017, Improving ambient
noise correlation functions with an SVD-based Wiener filter: Geophysical Journal
International, v. 211, no. 1, p. 418-426, doi: 10.1093/gji/ggx306.
178
Moschetti, M., Ritzwoller, M., Lin, F., and Yang, Y., 2010a, Crustal shear wave velocity
structure of the western United States inferred from ambient seismic noise and
earthquake data: Journal of Geophysical Research, v. 115, no. B10, doi:
10.1029/2010jb007448.
Moschetti, M., Ritzwoller, M., Lin, F., and Yang, Y., 2010b, Seismic evidence for widespread
western-US deep-crustal deformation caused by extension: Nature, v. 464, no. 7290, p.
885-889, doi: 10.1038/nature08951.
Najdahmadi, B., Bohnhoff, M., and Ben‐Zion, Y., 2016, Bimaterial interfaces at the Karadere
segment of the North Anatolian Fault, northwestern Turkey: Journal of Geophysical
Research: Solid Earth, v. 121, no. 2, p. 931-950, doi: 10.1002/2015jb012601.
Nakata, N., and Snieder, R., 2012, Estimating near-surface shear wave velocities in Japan by
applying seismic interferometry to KiK-net data: Journal of Geophysical Research: Solid
Earth, v. 117, no. B1, p. n/a-n/a, doi: 10.1029/2011jb008595.
Nakatani, M., and Scholz, C., 2004, Frictional healing of quartz gouge under hydrothermal
conditions: 1. Experimental evidence for solution transfer healing mechanism: Journal of
Geophysical Research: Solid Earth, v. 109, no. B7, doi: 10.1029/2001jb001522.
Obermann, A., Planès, T., Larose, E., Sens-Schönfelder, C., and Campillo, M., 2013, Depth
sensitivity of seismic coda waves to velocity perturbations in an elastic heterogeneous
medium: Geophysical Journal International, v. 194, no. 1, p. 372-382, doi:
10.1093/gji/ggt043.
Obermann, A., Froment, B., Campillo, M., Larose, E., Planès, T., Valette, B., Chen, J., and Liu,
Q., 2014, Seismic noise correlations to image structural and mechanical changes
associated with theMw7.9 2008 Wenchuan earthquake: Journal of Geophysical Research:
Solid Earth, v. 119, no. 4, p. 3155-3168, doi: 10.1002/2013jb010932.
Obermann, A., Planès, T., Hadziioannou, C., and Campillo, M., 2016, Lapse-time-dependent
coda-wave depth sensitivity to local velocity perturbations in 3-D heterogeneous elastic
media: Geophysical Journal International, v. 207, no. 1, p. 59-66, doi:
10.1093/gji/ggw264.
Olsen, K., Day, S., Minster, J., Cui, Y., Chourasia, A., Faerman, M., Moore, R., Maechling, P.,
and Jordan, T., 2006, Strong shaking in Los Angeles expected from southern San
179
Andreas earthquake: Geophysical Research Letters, v. 33, no. 7, doi:
10.1029/2005gl025472.
Onderdonk, N., Rockwell, T., McGill, S., and Marliyani, G., 2013, Evidence for Seven Surface
Ruptures in the Past 1600 Years on the Claremont Fault at Mystic Lake, Northern San
Jacinto Fault Zone, California: Bulletin of the Seismological Society of America, v. 103,
no. 1, p. 519-541, doi: 10.1785/0120120060.
Ozakin, Y., Ben-Zion, Y., Aktar, M., Karabulut, H., and Peng, Z., 2012, Velocity contrast across
the 1944 rupture zone of the North Anatolian fault east of Ismetpasa from analysis of
teleseismic arrivals: Geophysical Research Letters, v. 39, no. 8, doi:
10.1029/2012gl051426.
Ozakin, Y., and Ben-Zion, Y., 2015, Systematic Receiver Function Analysis of the Moho
Geometry in the Southern California Plate-Boundary Region: Pure and Applied
Geophysics, v. 172, no. 5, p. 1167-1184, doi:10.1007/s00024-014-0924-6.
Pei, S., Niu, F., Ben-Zion, Y., Sun, Q., Liu, Y., Xue, X., Su, J., and Shao, Z., 2019, Seismic
velocity reduction and accelerated recovery due to earthquakes on the Longmenshan
fault: Nature Geoscience, v. 12, no. 5, p. 387-392, doi: 10.1038/s41561-019-0347-1.
Peng, Z., Ben-Zion, Y., Michael, A.J., and Zhu, L., 2003, Quantitative analysis of seismic fault
zone waves in the rupture zone of the 1992 Landers, California, earthquake: evidence for
a shallow trapping structure: Geophysical Journal International, v. 155, no. 3, p. 1021-
1041.
Peng, Z., and Ben-Zion, Y., 2006, Temporal Changes of Shallow Seismic Velocity Around the
Karadere-Düzce Branch of the North Anatolian Fault and Strong Ground Motion: Pure
and Applied Geophysics, v. 163, no. 2-3, p. 567-600, doi: 10.1007/s00024-005-0034-6.
Petersen, M.D., and Wesnousky, S.G., 1994, Fault slip rates and earthquake histories for active
faults in southern California: Bulletin of the Seismological Society of America, v. 84, no.
5, p. 1608-1649.
Poli, P., Pedersen, H. A., and Campillo, M., 2012, Noise directivity and group velocity
tomography in a region with small velocity contrasts: the northern Baltic shield:
Geophysical Journal International, v. 192, no. 1, p. 413-424, doi: 10.1093/gji/ggs034.
Poupinet, G., Ellsworth, W., and Frechet, J., 1984, Monitoring velocity variations in the crust
using earthquake doublets: An application to the Calaveras Fault, California: Journal of
180
Geophysical Research: Solid Earth, v. 89, no. B7, p. 5719-5731, doi:
10.1029/jb089ib07p05719.
Prawirodirdjo, L., Ben-Zion, Y., and Bock, Y., 2006, Observation and modeling of thermoelastic
strain in Southern California Integrated GPS Network daily position time series: Journal
of Geophysical Research: Solid Earth, v. 111, no. B2, p. n/a-n/a, doi:
10.1029/2005jb003716.
Qin, L., Ben-Zion, Y., Qiu, H., Share, P. E., Ross, Z. E., and Vernon, F. L., 2018, 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,
no. 1, p. 98-114.
Qiu, H., Ben-Zion, Y., Ross, Z. E., Share, P. E., and Vernon, F. L., 2017, Internal structure of the
San Jacinto fault zone at Jackass Flat from data recorded by a dense linear
array: Geophysical Journal International, v. 209, no. 3, p. 1369-1388.
Qiu, H., Ben-Zion, Y., and Lin, F.-C., 2019a, Eikonal tomography of the Southern California
plate boundary region: Journal of Geophysical Research: Solid Earth (in review).
Qiu, H., Hillers, G., and Ben-Zion, Y., 2019b, Temporal changes of seismic velocities in the San
Jacinto Fault zone associated with the 2016 Mw5.2 Borrego Springs earthquake:
Geophysical Journal International (in review).
Qiu, H., Ben-Zion, Y., and Vernon, F. L., 2019c, Internal Structure of the San Jacinto Fault Zone
at Dry Wash from Data Recorded by a Dense Linear Array (in prep).
Ratdomopurbo, A., and Poupinet, G., 1995, Monitoring a temporal change of seismic velocity in
a volcano: application to the 1992 eruption of Mt. Merapi (Indonesia): Geophysical
Research Letters, v. 22, no. 7, p. 775-778, doi:10.1029/95GL00302.
Richter, T., Sens-Schönfelder, C., Kind, R., and Asch, G., 2014, Comprehensive observation and
modeling of earthquake and temperature-related seismic velocity changes in northern
Chile with passive image interferometry: Journal of Geophysical Research: Solid Earth,
v. 119, no. 6, p. 4747-4765, doi: 10.1002/2013jb010695.
Ritzwoller, M. H., Lin, F. C., and Shen, W., 2011, Ambient noise tomography with a large
seismic array: Comptes Rendus Geoscience, v. 343, no. 8-9, p. 558-570.
Rockwell, T.K., and Seitz, G.G., 2008, Observations of mode-switching from long paleoseismic
records of earthquakes on the San Jacinto and San Andreas faults: implications for
181
making hazard estimates from short paleoseismic records: Abstract presented at 33rd
International Geological Congress, Oslo, Norway, 6-14 August.
Rockwell, T., Dawson, T., Young Ben-Horin, J., and Seitz, G., 2015, A 21-Event, 4,000-Year
History of Surface Ruptures in the Anza Seismic Gap, San Jacinto Fault, and
Implications for Long-term Earthquake Production on a Major Plate Boundary Fault:
Pure and Applied Geophysics, v. 172, no. 5, p. 1143-1165, doi: 10.1007/s00024-014-
0955-z.
Ross, Z.E., and Ben-Zion, Y., 2014, Automatic picking of direct P, S seismic phases and fault
zone head waves: Geophysical Journal International, v. 199, no. 1, p. 368-381.
Ross, Z.E., and Ben-Zion, Y., 2015, An algorithm for automated identification of fault zone
trapped waves: Geophysical Journal International, v. 202, no. 2, p. 933-942.
Ross, Z. E., White, M. C., Vernon, F. L., and Ben‐Zion, Y., 2016, An improved algorithm for
real‐time S‐wave picking with application to the (augmented) ANZA network in
Southern California: Bulletin of the Seismological Society of America, v. 106, no. 5, p.
2013-2022.
Ross, Z., Hauksson, E., and Ben-Zion, Y., 2017a, Abundant off-fault seismicity and orthogonal
structures in the San Jacinto fault zone: Science Advances, v. 3, no. 3, p. e1601946, doi:
10.1126/sciadv.1601946.
Ross, Z., Kanamori, H., and Hauksson, E., 2017b, Anomalously large complete stress drop
during the 2016 M w 5.2 Borrego Springs earthquake inferred by waveform modeling
and near-source aftershock deficit: Geophysical Research Letters, v. 44, no. 12, p. 5994-
6001, doi: 10.1002/2017gl073338.
Roux, P., and Ben-Zion, Y., 2014, Monitoring fault zone environments with correlations of
earthquake waveforms: Geophysical Journal International, v. 196, no. 2, p. 1073-1081,
doi: 10.1093/gji/ggt441.
Roux, P., and Ben-Zion, Y., 2017, Rayleigh phase velocities in Southern California from
beamforming short-duration ambient noise: Geophysical Journal International, v. 211, no.
1, p. 450-454, doi:10.1093/gji/ggx316.
Rowe, C., and Griffith, W., 2015, Do faults preserve a record of seismic slip: A second opinion:
Journal of Structural Geology, v. 78, p. 1-26, doi: 10.1016/j.jsg.2015.06.006.
182
Sawazaki, K., Sato, H., Nakahara, H., and Nishimura, T., 2009, Time-Lapse Changes of Seismic
Velocity in the Shallow Ground Caused by Strong Ground Motion Shock of the 2000
Western-Tottori Earthquake, Japan, as Revealed from Coda Deconvolution Analysis:
Bulletin of the Seismological Society of America, v. 99, no. 1, p. 352-366, doi:
10.1785/0120080058.
SCEDC, 2013, Southern California Earthquake Data Center. Caltech Dataset
doi:10.7909/C3WD3xH.
Sens-Schönfelder, C., and Wegler, U., 2006, Passive image interferometry and seasonal
variations of seismic velocities at Merapi Volcano, Indonesia: Geophysical Research
Letters, v. 33, no. 21, doi: 10.1029/2006gl027797.
Sethian, J., 1996, A fast marching level set method for monotonically advancing fronts.:
Proceedings of the National Academy of Sciences, v. 93, no. 4, p. 1591-1595, doi:
10.1073/pnas.93.4.1591.
Shapiro, N. M., and Campillo, M., 2004, Emergence of broadband Rayleigh waves from
correlations of the ambient seismic noise: Geophysical Research Letters, v. 31, no. 7,
doi:10.1029/2004gl019491.
Shapiro, N. M., 2005, High-Resolution Surface-Wave Tomography from Ambient Seismic
Noise: Science, v. 307, no. 5715, p. 1615-1618, doi: 10.1126/science.1108339.
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, doi:10.1002/2016GL070774.
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, no. 2, p. 819-832.
Share, P.-E., and Ben-Zion, Y., 2018, A bimaterial interface along the northern San Jacinto fault
through Cajon Pass: Geophysical Research Letters, v. 45, no. 11, p. 622–11,631, doi:
10.1029/2018GL079834.
Share, P. E., Guo, H., Thurber, C. H., Zhang, H., and Ben-Zion, Y., 2019, Seismic Imaging of
the Southern California Plate Boundary around the South-Central Transverse Ranges
Using Double-Difference Tomography: Pure and Applied Geophysics, v. 176, no. 3, p.
1117-1143.
183
Sharp, R.V., 1967, The San Jacinto fault zone in the Peninsular Ranges of southern California:
Geological Society of America Bulletin, v. 78, p. 705-730.
Shaw, J. H., Plesch, A., Dolan, J. F., Pratt, T. L., and Fiore, P, 2002, Puente Hills Blind-Thrust
System, Los Angeles, California: Bulletin of the Seismological Society of America, v.
92, no. 8, p. 2946-2960.
Shaw, J., Plesch, A., Tape, C., Suess, M., Jordan, T., Ely, G., Hauksson, E., Tromp, J., Tanimoto,
T., Graves, R., Olsen, K., Nicholson, C., Maechling, P., and Rivero, C. et al., 2015,
Unified Structural Representation of the southern California crust and upper mantle:
Earth and Planetary Science Letters, v. 415, p. 1-15, doi: 10.1016/j.epsl.2015.01.016.
Shlomai, H., and Fineberg, J., 2016, The structure of slip-pulses and supershear ruptures driving
slip in bimaterial friction: Nature Communications, v. 7, no. 1, doi:
10.1038/ncomms11787.
Sibson, R. H., 1986, Rupture interaction with fault jogs: Earthquake Source Mechanics, v. 37, p.
157-167.
Silver, P., Daley, T., Niu, F., and Majer, E., 2007, Active Source Monitoring of Cross-Well
Seismic Travel Time for Stress-Induced Changes: Bulletin of the Seismological Society
of America, v. 97, no. 1B, p. 281-293, doi: 10.1785/0120060120.
Smith, M. L., and Dahlen, F. A., 1973, The azimuthal dependence of Love and Rayleigh wave
propagation in a slightly anisotropic medium: Journal of Geophysical Research, v. 78, no.
17, p. 3321-3333, doi:10.1029/JB078i017p03321.
Smith, W. H. F., and Wessel, P., 1990, Gridding with continuous curvature splines in tension,
Geophysics, v. 55, no. 3, p. 293-305.
Snieder, R., Grêt, A., Douma, H., and Scales, J., 2002, Coda Wave Interferometry for Estimating
Nonlinear Behavior in Seismic Velocity: Science, v. 295, no. 5563, p. 2253-2255, doi:
10.1126/science.1070015.
Stehly, L., Campillo, M., and Shapiro, N., 2006, A study of the seismic noise from its long-range
correlation properties: Journal of Geophysical Research, v. 111, no. B10, doi:
10.1029/2005jb004237.
Stehly, L., Froment, B., Campillo, M., Liu, Q., and Chen, J., 2015, Monitoring seismic wave
velocity changes associated with the Mw 7.9 Wenchuan earthquake: increasing the
184
temporal resolution using curvelet filters: Geophysical Journal International, v. 201, no.
3, p. 1939-1949, doi: 10.1093/gji/ggv110.
Taira, T., Nayak, A., Brenguier, F., and Manga, M., 2018, Monitoring reservoir response to
earthquakes and fluid extraction, Salton Sea geothermal field, California: Science
Advances, v. 4, no. 1, p. e1701536, doi: 10.1126/sciadv.1701536.
Tape, C., Liu, Q., Maggi, A., and Tromp, J., 2010, Seismic tomography of the southern
California crust based on spectral-element and adjoint methods: Geophysical Journal
International, v. 180, no. 1, p. 433-462, doi:10.1111/j.1365-246X.2009.04429.x.
Tsai, V.C., 2011, A model for seasonal changes in GPS positions and seismic wave speeds due
to thermoelastic and hydrologic variations: Journal of Geophysical Research, v. 116,
doi:10.1029/2010JB008156.
Vernon, F. and Ben-Zion, Y., 2010, San Jacinto Fault Zone Experiment. International Federation
of Digital Seismograph Networks: Other/Seismic Network, doi:10.7914/SN/YN_2010.
Vidale, J., Helmberger, D.V., and Clayton, R.W., 1985, Finite-difference seismograms for SH
waves: Bulletin of the Seismological Society of America, v. 75, no. 6, p. 1765-1782.
Wang, Q., Brenguier, F., Campillo, M., Lecointre, A., Takeda, T., and Aoki, Y., 2017, Seasonal
Crustal Seismic Velocity Changes Throughout Japan: Journal of Geophysical Research:
Solid Earth, v. 122, no. 10, p. 7987-8002, doi: 10.1002/2017jb014307.
Wathelet, M., 2008, An improved neighborhood algorithm: Parameter conditions and dynamic
scaling: Geophysical Research Letters, v. 35, no. 9, doi:10.1029/2008gl033256.
Wechsler, N., Rockwell, T.K., and Ben-Zion, Y., 2009, Analysis of rock damage asymmetry
from geomorphic signals along the trifurcation area of the San-Jacinto Fault:
Geomorphology, v. 113, p. 82-96.
Weng, H., Yang, H., Zhang, Z., and Chen, X., 2016, Earthquake rupture extents and coseismic
slips promoted by damaged fault zones: Journal of Geophysical Research: Solid Earth, v.
121, no. 6, p. 4446-4457, doi: 10.1002/2015jb012713.
Wdowinski, S., 2009, Deep creep as a cause for the excess seismicity along the San Jacinto fault:
Nature Geoscience, v. 2, no. 12, p. 882–885.
Weertman, J., 1980, Unstable slippage across a fault that separates elastic media of different
elastic constants: Journal of Geophysical Research: Solid Earth, v. 85, no. B3, p. 1455-
1461.
185
White, M. C. A., Ross, Z. E., Vernon, F. L., and Ben-Zion, Y., 2019, A detailed earthquake
catalog for the San Jacinto fault zone region in southern California: Journal of
Geophysical Research: Solid Earth, doi: 10.1029/2019JB017641.
Wielandt, E., 1993, Propagation and structural interpretation of non-plane waves: Geophysical
Journal International, v. 113, p. 45-53.
Wu, C., Peng, Z., and Ben-Zion, Y., 2009, Non-linearity and temporal changes of fault zone site
response associated with strong ground motion: Geophysical Journal International, v.
176, no. 1, p. 265-278, doi: 10.1111/j.1365-246x.2008.04005.x.
Wu, C., Peng, Z., and Ben-Zion, Y., 2010, Refined thresholds for non-linear ground motion and
temporal changes of site response associated with medium-size earthquakes: Geophysical
Journal International, v. 182, no. 3, p. 1567-1576, doi: 10.1111/j.1365-
246x.2010.04704.x.
Xu, H., Luo, Y., Chen, C., and Xu, Y., 2016, 3D shallow structures in the Baogutu area,
Karamay, determined by eikonal tomography of short-period ambient noise surface
waves: Journal of Applied Geophysics, v. 129, p. 101-110.
Xu, S., Ben-Zion, Y., and Ampuero, J.P., 2012, Properties of inelastic yielding zones generated
by in-plane dynamic ruptures—II. Detailed parameter-space study: Geophysical Journal
International, v. 191, no. 3, p. 1343-1360.
Yao, H., and van der Hilst, R. D., 2009. Analysis of ambient noise energy distribution and phase
velocity bias in ambient noise tomography, with application to SE Tibet: Geophysical
Journal International, v. 179, no. 2, p. 1113-1132.
Yang, C., G. Li, F. Niu and Y. Ben-Zion (2019). Significant effects of shallow seismic and stress
properties on phase velocities of Rayleigh waves up to 20 s Pure Appl. Geophys., 176,
doi: 10.1007/s00024-018-2075-7.
Yang, H., and Zhu, L., 2010, Shallow low-velocity zone of the San Jacinto fault from local
earthquake waveform modelling: Geophysical Journal International, v. 183, p. 421-432.
Yang, H., Li, Z., Peng, Z., Ben-Zion, Y., and Vernon, F., 2014, Low-velocity zones along the
San Jacinto Fault, Southern California, from body waves recorded in dense linear arrays:
Journal of Geophysical Research: Solid Earth, v. 119, no. 12, p. 8976-8990, doi:
10.1002/2014jb011548.
186
Yang, W., Peng, Z., Wang, B., Li, Z., and Yuan, S., 2015, Velocity contrast along the rupture
zone of the 2010 Mw6. 9 Yushu, China, earthquake from fault zone head waves: Earth
and Planetary Science Letters, v. 416, p. 91-97.
Zhu, L., and Kanamori, H., 2000, Moho depth variation in southern California from teleseismic
receiver functions: Journal of Geophysical Research: Solid Earth, v. 105, no. B2, p. 2969-
2980, doi:10.1029/1999jb900322.
Zaliapin, I., and Ben-Zion, Y., 2011, Asymmetric distribution of aftershocks on large faults in
California: Geophysical Journal International, v. 185, p. 1288-1304, doi: 10.1111/j.1365-
246X.2011.04995.x.
Zhan, Z., Tsai, V., and Clayton, R., 2013, Spurious velocity changes caused by temporal
variations in ambient noise frequency content: Geophysical Journal International, v. 194,
no. 3, p. 1574-1581, doi: 10.1093/gji/ggt170.
Zhao, P., and Peng, Z., 2008, Velocity contrast along the Calaveras fault from analysis of fault
zone head waves generated by repeating earthquakes: Geophysical Research Letters, v.
35, no. 1.
Zhao, P., Peng, Z., Shi, Z., Lewis, M.A., and Ben-Zion, Y., 2010, Variations of the velocity
contrast and rupture properties of M6 earthquakes along the Parkfield section of the San
Andreas fault: Geophysical Journal International, v. 180, no. 2, p. 765-780.
Zigone, D., Ben-Zion, Y., Campillo, M., and Roux, P., 2015, Seismic Tomography of the
Southern California Plate Boundary Region from Noise-Based Rayleigh and Love
Waves: Pure and Applied Geophysics, v. 172, no. 5, p. 1007-1032, doi: 10.1007/s00024-
014-0872-1.
Zigone, D., Ben-Zion, Y., Lehujeur, M., Campillo, M., Hillers, G., and Vernon, F. L., 2019,
Imaging subsurface structures in the San Jacinto fault zone with high frequency noise
recorded by dense linear arrays: Geophysical Journal International, v. 217, p. 879–893,
doi: 10.1093/gji/ggz069.

Abstract (if available)

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

PDF

Multi-scale imaging of major fault zones in Southern California

PDF

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

PDF

Detailed properties of seismic waveforms and earthquake sources in Southern California

PDF

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

PDF

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

PDF

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

PDF

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

PDF

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

PDF

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

PDF

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

PDF

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

PDF

Stress-strain characterization of complex seismic sources

PDF

Interaction between dynamic ruptures and off-fault yielding characterized by different rheologies

PDF

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

PDF

Seismicity distribution near strike-slip faults in California

PDF

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

PDF

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

PDF

Direct observation of fault zone structure and mechanics in three dimensions: a study of the SEMP fault system, Austria

PDF

Observations and modeling of dynamically triggered high frequency burst events

Asset Metadata

Creator Qiu, Hongrui (author)

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

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Geological Sciences

Publication Date 07/18/2019

Defense Date 04/17/2019

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

Tag body waves,coda waves,earthquake dynamics,guided waves,interface waves,OAI-PMH Harvest,rheology and friction of fault zones,seismic interferometry,seismic noise,surface waves

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Ben-Zion, Yehuda (committee chair), Pilch, Krzysztof (committee member), Sammis, Charles (committee member), West, Joshua (committee member), Wu, Francis (committee member)

Creator Email hongruiq@usc.edu,qiuhonrui@gmail.com

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

Unique identifier UC11660102

Identifier etd-QiuHongrui-7555.pdf (filename),usctheses-c89-182946 (legacy record id)

Legacy Identifier etd-QiuHongrui-7555.pdf

Dmrecord 182946

Document Type Dissertation

Format application/pdf (imt)

Rights Qiu, Hongrui

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