University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content SEISMIC VELOCITY STRUCTURE AND SPATIOTEMPORAL PATTERNS OF
MICROSEISMICITY IN ACTIVE FAULT ZONES OF SOUTHERN CALIFORNIA
by
Malcolm Charles Adan White
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulllment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(GEOLOGICAL SCIENCES)
August 2021
Copyright 2021 Malcolm Charles Adan White
Dedication
I dedicate this work to my greatest supporter, my wife Sara.
ii
Acknowledgements
Thank you to my wife, Sara, for your unwavering support, for your sacricial service of our family, and for
your patience when this work consumed our evenings and weekends. Thank you to my parents, Donald
and Ruth, for providing me with the innumerable advantages that enabled me to pursue an intellectually
stimulating profession and for bearing with me when I questioned my desire to do so in times of vexation.
Thank you to my church for being there to support me independent of my line of work. Thank you to my
advisor, Yehuda, for working with me in an unorthodox arrangement. Thank you to my co-authors, Frank,
Hongjian, and Nori, for many fruitful discussions.
Finally, I give my thanks to Jesus Christ, the A and
in whom we live and move and have our being.
You drew me up from the pit of destruction, out of the miry bog, and set my feet upon a rock, making my
steps secure. To you be glory, majesty, dominion, and authority, before all time and now and forever.
Great are the works of the Lord, studied by all who delight in them. Full of splendor and majesty is his
work, and his righteousness endures forever. —Psalm 111:2-3
iii
TableofContents
Dedication ii
Acknowledgements iii
ListofTables viii
ListofFigures ix
Abstract xii
Chapter1: Introduction 1
Chapter2: AutomatedcatalogingofmicroseismicityintheSanJacintoFaultZone 5
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.1 Seismic stations around the San Jacinto fault zone . . . . . . . . . . . . . . . . . . . 8
2.2.2 Choosing and downloading the data . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2.3 Velocity model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.4 Analyst reviewed control dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3.2 Phase I—Detecting and locating earthquakes . . . . . . . . . . . . . . . . . . . . . . 12
2.3.2.1 Detecting P-waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3.2.2 Detecting S-waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3.2.3 Associating phase detections with events . . . . . . . . . . . . . . . . . . 15
2.3.2.4 Locating events individually . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3.3 Phase II—Relocating earthquakes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3.3.1 Cross-correlating waveforms . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3.3.2 Relocating with a double-dierence algorithm . . . . . . . . . . . . . . . 17
2.3.4 Phase III—Calculating event sizes and post-processing . . . . . . . . . . . . . . . . 17
2.3.4.1 Quality control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3.4.2 Calculating local magnitude, M
L
. . . . . . . . . . . . . . . . . . . . . . . 18
2.3.4.3 Calculating moment magnitudes, M
w
. . . . . . . . . . . . . . . . . . . . 18
2.3.4.4 Estimating minimum magnitude of catalog completeness, M
C
. . . . . . 19
2.3.4.5 Analyzing similar-event chains . . . . . . . . . . . . . . . . . . . . . . . 19
2.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.4.1 Basic features of the catalog . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
iv
2.4.2 Structural features in the catalog . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.4.2.1 General structural features . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.4.2.2 Hot Springs area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.4.2.3 Trifurcation area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.4.2.4 Anza Seismic Gap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.5.1 Characteristics of the catalog . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.5.2 Characteristics of the fault zone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
Chapter3: Thetheoreticalbasisfortraveltimetomographyinanisotropic,linearlyelastic
continuum 41
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.2 Elastodynamic equations of motion in a linearly elastic, isotropic continuum . . . . . . . . 43
3.3 The eikonal equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.4 Fermat’s principle and traveltime tomography . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.5 Fundamentals of traveltime tomography . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
Chapter4: Solvingtheeikonalequation 57
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.2.1 Deriving the eikonal equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.2.2 Statement of the problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.2.3 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.2.4 The Fast Marching Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.2.5 Coordinate systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.2.6 Order of the nite-dierence operators . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.2.7 Ray tracing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.3.1 Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.3.2 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.3.3 Speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.3.4 Generality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.3.5 Extensibility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.3.5.1 Tracking secondary phases . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.3.5.2 Locating earthquakes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.4.1 Comparing with FM3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.4.2 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.4.3 Jupyter Notebook Examples and Documentation . . . . . . . . . . . . . . . . . . . 85
4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.6 Data and Resources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.7 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
Chapter5: TraveltimetomographyfromanautomatedcatalogofRidgecrestaftershocks 87
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
5.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
v
5.3.1 Detecting earthquakes and measuring phase arrival times . . . . . . . . . . . . . . 94
5.3.2 Associating arrivals with earthquakes . . . . . . . . . . . . . . . . . . . . . . . . . 95
5.3.3 Locating earthquakes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
5.3.4 Estimating local event magnitudes,M
L
, and magnitude of catalog completeness,
M
C
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
5.3.5 Traveltime tomography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
5.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
5.4.1 Earthquake catalog . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
5.4.1.1 Eastern Little Lake Fault . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5.4.1.2 Southern Little Lake Fault . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.4.1.3 The Garlock Fault and Fremont Valley . . . . . . . . . . . . . . . . . . . 104
5.4.1.4 Coso Range . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.4.1.5 Argus Mountains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
5.4.2 Velocity models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
5.4.2.1 Model residuals, resolution, and robustness . . . . . . . . . . . . . . . . . 105
5.4.2.2 Surface geology and shallow velocities . . . . . . . . . . . . . . . . . . . 107
5.4.2.3 Horizontal slices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
5.4.2.4 Structural variability parallel to the ELLF . . . . . . . . . . . . . . . . . . 114
5.4.2.5 Fault-normal sections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
5.4.3 Key structural features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
5.4.3.1 Central ELLF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
5.4.3.2 Northwest ELLF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
5.4.3.3 Southeast ELLF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
5.4.3.4 Garlock Fault at Fremont Valley . . . . . . . . . . . . . . . . . . . . . . . 126
5.5 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
Chapter6: UpdatingthecatalogfortheSanJacintoFaultZone 132
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
6.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
6.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
6.3.1 Detecting earthquakes and picking phase arrivals . . . . . . . . . . . . . . . . . . . 136
6.3.2 Associating phase arrivals with earthquake sources . . . . . . . . . . . . . . . . . . 136
6.3.3 Estimating absolute event locations . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
6.3.4 Relocating event ensembles with double dierences . . . . . . . . . . . . . . . . . . 141
6.3.5 Calculating local magnitudes,M
L
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
6.3.6 Quality control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
6.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
6.4.1 Catalog overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
6.4.2 Selected structural features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
6.4.2.1 Hot Springs Area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
6.4.2.2 Trifrucation Region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
6.4.2.3 Cahuilla Swarm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
6.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
6.5.1 Prior probabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
6.5.2 Derivative studies and further work . . . . . . . . . . . . . . . . . . . . . . . . . . 152
6.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
Chapter7: Discussion 154
vi
Bibliography 156
AppendixA:Detectingandpickingphasearrivals 170
A.1 Measuring P-wave arrival times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
A.2 Measuring S-wave arrival times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
vii
ListofTables
2.1 Parameters for detecting P-waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2 Phase detection rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3 Catalog quality metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.4 Similar-event chains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.1 PyKonal execution time for problems of various grid sizes. . . . . . . . . . . . . . . . . . . 76
5.1 Traveltime tomography parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
viii
ListofFigures
1.1 Kernel density estimate of seismic activity in Southern California since 2008. . . . . . . . 2
2.1 Network map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Network history . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Automated processing procedure schematic . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.4 Frequency-magnitude distribution model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.5 Seismicity overvieww map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.6 Seismicity detection rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.7 Location error maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.8 Magnitude of completeness maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.9 Map of similar-event locations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.10 Vertical transect (fault parallel) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.11 Vertical transect (fault perpendicular; Hot Springs area) . . . . . . . . . . . . . . . . . . . . 29
2.12 Vertical transect (fault perpendicular; Trifurcation area) . . . . . . . . . . . . . . . . . . . . 30
2.13 Similar-event chain example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.14 Vertical transect (fault parallel; S-wave velocity background) . . . . . . . . . . . . . . . . . 32
2.15 Shatter networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.1 Essential elements of the FMM algorithm for solving the eikonal equation implemented
by PyKonal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.2 ISO coordinate-system convention adopted for spherical coordinates. . . . . . . . . . . . . 66
ix
4.3 Comparison of fraveltime elds and absolute errors computed for a point source using
Cartesian and spherical coordinates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.4 Comparison of fraveltime elds and absolute errors computed for a line source using
Cartesian and spherical coordinates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.5 Convergent behavior of traveltimes computed in a complex 2D velocity model using
PyKonal in Cartesian coordinates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.6 Convergence to analytical solution of ray paths computed using PyKonal. . . . . . . . . . 75
4.7 Convergent behavior of traveltimes computed in a complex 2D velocity model using
PyKonal in spherical coordinates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.8 Tracking reected waves using PyKonal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.9 Comparison of absolute errors for traveltime elds computed using PyKonal and FM3D. . 83
5.1 Network map. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.2 Seismicity overview map. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
5.3 Frequency-magnitude distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
5.4 Event detection rate timeline. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5.5 Model residual distributions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
5.6 Example checkerboard tests. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
5.7 Surface geology and shallow velocity structure. . . . . . . . . . . . . . . . . . . . . . . . . 110
5.8 Horizontal slices ofV
p
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
5.9 Horizontal slices ofV
s
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
5.10 Horizontal slices of
Vp
=Vs model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
5.11 Fault-parallel transects ofV
p
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
5.12 Fault-parallel transects ofV
s
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
5.13 Fault-parallel transects of
Vp
=Vs model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
5.14 Fault-normal transects ofV
p
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
5.15 Fault-normal transects ofV
s
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
x
5.16 Fault-normal transects of
Vp
=Vs model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
5.17 Structural characteristics of the central segment of the Eastern Little Lake Fault. . . . . . 123
5.18 Structural characteristics of the northwest terminus of the Eastern Little Lake Fault. . . . 124
5.19 Structural characteristics of the southeast segment of the Eastern Little Lake Fault. . . . . 125
5.20 Structural characteristics of the Garlock Fault at the Fremont Valley. . . . . . . . . . . . . 126
6.1 Station map. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
6.2 Example error model used to estimate hypocenter coordinates using MCMC sampling. . . 140
6.3 Example 2D and 1D marginalized posterior sample distributions for hypocenter
coordinates sampled using MCMC sampling. . . . . . . . . . . . . . . . . . . . . . . . . . 142
6.4 Histograms of hypocenter-parameter uncertainties for all cataloged events. . . . . . . . . 144
6.5 Overview map of cataloged seismicity in our focus region. . . . . . . . . . . . . . . . . . . 145
6.6 Comparative histograms of events cataloged per month in our focus region for this study,
White et al. (2019), and Hauksson et al. (2012). . . . . . . . . . . . . . . . . . . . . . . . . . 146
6.7 Map of cataloged seismicity near the Hot Springs Area. . . . . . . . . . . . . . . . . . . . . 147
6.8 Map of cataloged seismicity in the Trifurcation Region. . . . . . . . . . . . . . . . . . . . . 149
6.9 Maps of cataloged seismicity in a sub-region of the Trifurcation Region. . . . . . . . . . . 150
6.10 Map of cataloged seismicity near Cahuilla, California. . . . . . . . . . . . . . . . . . . . . 151
A.1 (a) 3C seismic waveform data with three registered P-wave arrivals (vertical dashed blue
lines) and (b) derived statistics used to compute the (c) characteristic function (CF) and
dynamic threshold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
xi
Abstract
I develop and apply automated processing procedures to derive seismicity catalogs from raw seismograms
for two highly active fault zones in Southern California: (a) the San Jacinto Fault Zone (SJFZ) and (b)
the region surrounding the 2019 Ridgecrest, California, earthquake sequence. For the Ridgecrest region,
I also derive 3D models of seismic P- and S-wave velocity (V
p
andV
s
, respectively) structure, and their
ratio (
Vp
=Vs). The SJFZ catalog comprises more than 2.510
5
earthquakes in a 410
3
km
2
region between
1 January 2008 and 1 January 2021. The Ridgecrest catalog comprises more than 9.510
4
earthquakes in
a four month period including the rst three months of the aftershock sequence. Jointly interpreting spa-
tiotemporal seismicity patterns, geologic data, and velocity models provides new insights into subsurface
structure of and physical processes operating in the fault zones. To build the presented catalogs and ve-
locity models, I (a) develop and implement a new algorithm to detect earthquakes and pick phase arrivals,
(b) implement the Fast Marching Method to solve the eikonal equation in 3D media, (c) develop and im-
plement two new algorithms for locating earthquakes in 3D media, and (d) modify a recently developed
traveltime tomography method to adaptively resolve structural features based on data distribution. The
entire framework developed is modular, enabling components to be replaced as methods advance, and
can be applied to a wide range of network geometries. The presented methodological developments, data
products, and interpretations will benet various future studies of seismotectonics and crustal structure
in Southern California and elsewhere.
xii
Chapter1
Introduction
The boundary between the Pacic and North American tectonic plates transects Southern California and
threatens its over 23 million residents and infrastructure with the potential to unleash devastating earth-
quakes. To eectively mitigate the associated seismic hazard and plan for the future, we need a compre-
hensive understanding of the subsurface structure of the various fault zones it comprises and the physical
processes that operate within them. But, the vast majority of seismic activity in Southern California oc-
curs between 4 and 16 km below sea level, and we are conned to observations at or very near Earth’s
surface. The San Andreas Fault Observatory at Depth (SAFOD), for example, reaches a maximum depth
of 3.2 km below the surface. So how are we to infer the physical properties of and processes operating in
the subsurface suciently to prepare for the inevitable if we cannot observe either of them directly?
Since the early 1900s, seismologists have sought to answer this question by studying the potentially
destructive seismic waves—elastically transmitted physical disturbances that propagate through Earth’s
interior and over its surface—themselves to infer properties of their earthquake sources and the media
through which they propagate. The San Andreas Fault (SAF), spanning roughly 1200 km from just east
of the Salton Sea (about 80 km north of Mexico) to Cape Mendocino (on the Californian coast, 180 km
south of the Oregon border), is nominally considered the plate boundary and has thus justiably received
considerable attention from the public and scientists alike. Although the SAF is thelargest fault, and thus
1
Figure 1.1: Kernel density estimate of seismic activity in Southern California since 2008. The thick blue
line marks the San Andreas Fault, the nominal boundary between the oceanic Pacic plate to the west and
the continental North American plate to the east. Thin, white lines indicate other known faults.
capable of producing the largest earthquakes, the San Jacinto Fault Zone (SJZF) is the most seismically
active fault zone in Southern California, making it a critical component of the plate-boundary system.
Furthermore, diverse fault zones in Southern California have historically hosted large (magnitudeM > 6)
earthquakes, including, most recently, a pair of earthquakes (moment magnitudesM
w
6:4 andM
w
7:1) near
the town of Ridgecrest, California in July 2019. Evidently, the plate-boundary system is complex and its
many components must be understood holistically and integrated into a unied framework.
What is the subsurface structure of the highly active SJFZ? What can we learn about future large
earthquakes from the Ridgecrest earthquake pair? These are key questions that this thesis aims to address.
2
This thesis answers these questions by analyzing microseismicity (small earthquakes) in the SJFZ and
during the aftershock sequence following the Ridgecrest earthquake pair. Microearthquakes are a par-
ticularly powerful source of information sheerly because of their quantity; they are exponentially more
abundant than their larger counterparts. And, because faults are, by denition, where earthquakes oc-
cur, careful analysis of microseismicity can yield detailed images of fault-zone structure, and studying
their spatial and temporal (spatiotemporal) distributions can reveal insights into the nature of the physical
processes that generate them.
The research presented in this thesis comprises two main case studies: investigating fault-zone struc-
ture revealed by (a) a decade of microseismicity in the SJFZ, and (b) three months of intense aftershock
activity following the Ridgecrest earthquake pair. In the rst case, we develop and apply an automated
processing procedure to nine years of archived seismic data from the SJFZ to build a database of over
100 000 earthquakes. We subsequently analyze spatiotemporal patterns of seismicity in this database to
infer physical properties of the SJFZ. In the second case, we adapt our automated procedure from the rst
case to build a database for the aftershock sequence of the Ridgecrest earthquake pair and then use this
database to create 3D models of seismic wavespeeds, which serve as proxies for other material properties
(e.g., rigidity, temperature, and uid content). Jointly analyzing the resultant database and velocity models
leads us to a set of hypotheses for explaining the cessation of the main earthquakes, which will help guide
future research seeking cessation patterns of large earthquakes.
Prior to the work presented in this thesis, all available decade-scale catalogs of seismic activity in
the SJFZ were based primarily on labor-intensive manual analysis by seismic analysts, and were thus
relatively incomplete—meaning they lack many of the microearthquakes we recovered in our research. We
engineered an automated procedure to derive a catalog of microseismicity for the SJFZ from raw seismic
waveform data that can be readily modied and adapted for other study regions, as we demonstrate in the
Ridgecrest study. In the Ridgecrest study, we contribute the most comprehensive analysis of the aftershock
3
sequence to date and the most detailed models of seismic velocity structure for the region. The catalogs
and models obtained will benet future studies in the SJFZ and Ridgecrest regions.
This document is divided into ve main content chapters (excluding this Introduction and the Discus-
sion chapter). Chapter 2 introduces our preliminary study of microseismicity in the SJFZ (White et al. 2019).
Chapter 3 presents the theoretical basis for ray-based tomography. Chapter 4 presents a critical piece of
software that we developed based on theory in the preceding chapter and in support of subsequent work.
Chapter 5 presents the Ridgecrest study. Chapter 6 presents ongoing work to update and extend the SJFZ
catalog with additional data and rened methods. Finally, Chapter 7 oers a brief, summarizing discussion
of the thesis. A preface precedes each of Chapters 2-6 to unify what were originally independent projects
into a single narrative. I recommend that the reader start with this Introduction, the preface to each chap-
ter, and the Discussion to orient themselves to the overall trajectory of the thesis before engaging with
any one chapter. Otherwise, the unifying narrative may remain unclear.
4
Chapter2
AutomatedcatalogingofmicroseismicityintheSanJacintoFaultZone
Preface
The content of this chapter comes from an article published in the Journal of Geophysical Research: Solid
Earth in 2019 under the titleAdetailedearthquakecatalogfor theSanJacintofault-zoneregioninsouthern
California (White et al. 2019). It presents an automated processing procedure for deriving a seismicity
catalog from raw seismic waveform data and application of it to nine years of data (2008 through 2016) to
obtain a catalog for the SJFZ with over 100 000 earthquakes. We interpret spatiotemporal microseismicity
patterns to infer (a) the geometry of the brittle-to-ductile transition zone—a region over which the Earth’s
crust undergoes a categorical change in mechanical properties (i.e., from brittle to plastic); (b) that diverse
deformation processes occur within the fault zone (i.e., modeling the fault zone with a small number of
planar surfaces is inadequate); and (c) that, in certain instances, a sequence of similar events that links two
dissimilar events exists, demonstrating that the diverse processes occurring within the fault zone are not
all categorically dierent (e.g., the spectrum of observed deformation patterns is continuous, not discrete).
The procedure developed is modular, making it suitable to subsequent modication and update, as will be
demonstrated in Chapters 5 and 6.
The work presented in this chapter is the predecessor of ongoing work, which is presented in 6, to
extend and update the catalog developed in this chapter using additional data and rened methods.
5
Abstract
We develop an automated processing procedure to derive a new catalog of earthquake locations, magni-
tudes, and potencies, and analyze nine years of data, between 2008 and 2016 in the San Jacinto fault zone
region. Our procedure accounts for detailed 3D velocity structure using a probabilistic, global-search loca-
tion inversion and obtains high-precision relative event locations using dierential travel times measured
by cross-correlating waveforms. The obtained catalog illuminates spatiotemporal seismicity patterns in
the fault zone with observations for 108,800 earthquakes in the magnitude range -1.8 to 5.4. Inside a focus
region consisting of an 80-km by 50-km rectangle oriented parallel to the main fault trace, we estimate a
99% detection rate of earthquakes with magnitude 0.6 and greater and detect and locate about 60% more
events than those present in the Southern California Seismic Network catalog (Hutton et al. 2010). The
results provide the most complete catalog available for the region during the analyzed period, and include
both deeper events and very shallow patches of seismicity not present in the Hauksson et al. (2012) catalog.
The seismicity exhibits a variety of complex patterns, from which we interpret mechanical properties of
the host rock, and the fraction of event pairs with waveforms having cross-correlation coecients 0.95
is only about 3%, indicating diverse processes operating in the fault zone.
2.1 Introduction
Catalogs of earthquake locations, occurrence times and sizes are foundational for numerous studies of
earthquake physics, crustal dynamics and structural seismology. They dene the existence and geometry
of seismogenic faults at depth, provide information on spatiotemporal seismicity patterns, and are the
basis of many derivative studies including body-wave tomography and analyses of seismic hazard. More
detailed earthquake catalogs directly increase the available information on brittle-failure processes, crustal
structures and related topics.
6
The San Jacinto fault zone (SJFZ) is the most seismically-active fault zone in Southern California
(Hauksson et al. 2012; Ross et al. 2017; Sanders and Kanamori 1984) during the modern instrumental pe-
riod, and has ruptured in fteen M>7 earthquakes over the past 4,000 years (Rockwell et al. 2015). The
fault zone accommodates 10-21 mm/yr of the plate motion in Southern California, and has accumulated 25
km of slip (Fialko 2006; Sharp 1967). The SJFZ has signicant scientic and societal relevance because its
high seismic activity can provide detailed information on seismogenic processes in a complex active plate
boundary region, while posing signicant seismic hazard to large urban areas.
Since 1932, the Southern California Seismic Network (SCSN) has been the authoritative agency for
monitoring seismicity in Southern California; the agency currently operates over 400 seismic stations
and regularly updates a catalog of analyst-reviewed earthquake parameters (Hutton et al. 2010), referred
to here as SCSN_catalog_2010. A series of studies over the past two decades (Hauksson 2000; Hauksson
and Shearer 2005; Hauksson et al. 2012; Lin et al. 2007; Shearer et al. 2005) relocated seismicity from the
SCSN catalog from 1981, using 3D velocity information, a source-specic station-term relocation method
and a double-dierence relocation based on waveform cross-correlations. These studies culminated in a
waveform-relocated catalog (Hauksson et al. 2012), referred to here asHYS_catalog_2011, which has been
extended to later years and provides the standard earthquake catalog for Southern California.
The SCSN catalog covers the entire Southern California region, and omits some data from more focused
temporary local networks that make negligible dierence for routine processing on a regional scale. The
SJFZ Experiment Network (Vernon and Ben-Zion 2010) includes 35 stations near the central portion of the
SJFZ (Figure 2.1) that allow derivation of a more detailed earthquake catalog in that region. In the present
study, we use data from the SJFZ Experiment and additional four networks operating near the SJFZ to
obtain a high-resolution earthquake catalog for the region marked by a black, dashed rectangle in Figure
2.1. In the following sections 2.2 and 2.3 we describe in detail the data and methods used to derive the
catalog. The resulting catalog is used in section 2.4 to illuminate seismogenic structures in the SJFZ, and
7
117.5°W 117°W 116.5°W 116°W 115.5°W
33°N
33.5°N
34°N
34.5°N
San Diego
Anza
Palm Springs
Network
AZ
CI
PB
SB
YN
Figure 2.1: Map of seismic stations (inverted triangles color-coded by operating network; see legend at
bottom left) used in this study, surface traces of known Quarternary faults (thin black lines), and our focus
region (black, dashed rectangular outline). The cities of San Diego, Anza, and Palm Springs are shown for
reference. The region mapped in the large-scale map corresponds to the solid, red rectangular outline in
the small-scale map (top right).
the observed features are discussed in section 2.5 in relation to tomographic results and some mechanical
models for the region.
2.2 Data
2.2.1 SeismicstationsaroundtheSanJacintofaultzone
Five seismic networks operating in the SJFZ region with complementary capabilities are used in this study
(Figures 2.1 and ??): the Anza Network (AZ; Vernon 1982), the SJFZ Experiment Network (YN; Vernon
8
and Ben-Zion 2010), the Plate Boundary Observatory Borehole Network (PB), the Southern California
Seismic Network (CI; California Institute of Technology and United States Geological Survey Pasadena
1926; Southern California Earthquake Center 2013), and the UC Santa Barbara Engineering Seismology
Network (SB). The AZ network is the longest running network of those that focus on the SJFZ, operating
since 1982 and archiving continuous waveforms since 1998. The YN network is the densest network of
those focusing on the SJFZ, with 65 stations deployed as either stand-alone stations or in fault-normal linear
arrays. The PB network contributes high-quality data from eight high-frequency geophones installed in
boreholes between 5 and 230 m deep. The CI network is the most extensive network, covering the whole
of Southern California and providing a regional backbone for this study. The SB network operates two
down-hole arrays in the SJFZ, from which we use data from one instrument per array.
2.2.2 Choosinganddownloadingthedata
Because our aim is to catalog SJFZ microseismicity, and the AZ and YN networks target this region, we
build an aggregate network called here the Augmented Anza Network (AZ+) using stations from these as
the core (Figure 2.1). We include data from all ve networks previously mentioned, and choose stations to
evenly cover the fault zone area. The SCSN archived continuous waveforms at maximum sample rates of
40 samples-per-second (s.p.s.) prior to 2008, so we begin our analysis in 2008 when 100 s.p.s. continuous
SCSN data archives became available.
Many stations have co-located instruments (e.g., broadband seismometers and strong-motion accelerom-
eters), and digital data are sometimes archived at multiple sample rates. Because this study aims to detect
small seismic events, we preferentially select a single data stream for each station based on two criteria:
i) instrument type and ii) sample rate. Data from broadband seismometers are most preferred, followed
in order of decreasing preference by intermediate-period seismometers, high-frequency geophones, and
strong-motion accelerometers. The highest sample-rate data available (up to 250 s.p.s.) is preferred, and no
9
Figure 2.2: The number of stations used in this study as a function of time. The colored, solid lines represent
the number of stations stratied by operating network (see legend; right middle) and are scaled relative to
the left vertical axis. The black, dashed line represents the total number of stations and is scaled relative
to right vertical axis. The increasing number of stations between 2008 and 2014 is chiey the result of new
station installation. The rapid expansion of the YN network starting in late 2011 is particularly important
for this study.
10
data with sample rate lower than 40 s.p.s. is used in this study. Wherever possible, we use low-preference
data streams to ll gaps in preferred data streams that are longer than 24 hours.
Waveform data from 116 stations comprise the nal 12.5 TB dataset. Most of the data were downloaded
from archives at UC San Diego; however, some gaps were lled with data downloaded with the Seismogram
Transfer Program (STP; used for data from the CI network) and the BREQ_FAST web resource provided
by the Incorporated Research Institutions for Seismology’s Data Management Center (IRIS-DMC; used for
all other networks). All data from the CI network used in this study can be downloaded using STP, and all
other data can be downloaded from the IRIS-DMC using BREQ_FAST.
2.2.3 Velocitymodel
The three-dimensional velocity structure in the region is modeled with a hybrid framework derived from
(Fang et al. 2016), referred to as FANG16, and the Southern California Earthquake Center’s (SCEC’s) Com-
munity Velocity Model Harvard version 15.1, referred to as CVM-H15.1. The FANG16 model was obtained
directly from its authors and CVM-H15.1 was obtained using the SCEC UCVM software. We combine
these two models because FANG16 resolves detailed structure within and around the SJFZ, but lacks reso-
lution along the edges of our study area; CVM-H15.1 complements FANG16 by resolving regional velocity
structure well, but lacks the detail of FANG16 inside the focus region of our study.
The hybrid velocity model (Supplementary Figures S1, S2) is dened on the same set of grid nodes as
FANG16 and is derived by computing a weighted average of FANG16 and CVM-H15.1. Complementary
weight functions are dened using a 2D Gaussian distribution oriented parallel to the fault zone (Supple-
mentary Figure S1) so that FANG16 is heavily weighted in the SJFZ and CVM-H15.1 is heavily weighted
near the edges of the study area. The uppermost layer of FANG16 is located at -1.5 km (i.e., 1.5 km above
sea level) and CVM-H15.1 is not dened for negative depths, so we set the -1.5 km layer for CVM-H15.1
equal to 0.9 times the 0.0 km layer. An additional constant velocity layer is added to the nal model at 30.0
11
km depth using velocities from the IASPEI91 model (Kennett and Engdahl 1991); 6.5 km/s and 3.75 km/s
for P- and S-waves respectively.
2.2.4 Analystreviewedcontroldataset
To validate the performance of dierent processing parameters, we compare intermediate results with an
unpublished, 40-day, analyst-reviewed catalog of phase arrival times and earthquake locations observed on
data of the AZ+ network, referred to as ANZA_review. This catalog contains 1,200 events, 36,608 P-wave
arrival times, and 32,590 S-wave arrival times, each reviewed by an analyst with the goal of registering
every readily-identiable phase arrival and locating all associated earthquakes.
2.3 Methods
2.3.1 Overview
The processing procedure used in this study is presented in three main phases (Figure 2.3), each broken
down into individual steps in the following subsections: Phase I—Detecting and locating earthquakes;
Phase II—Relocating earthquakes; and Phase III—Calculating event sizes (potencies, magnitudes) and post-
processing. For reproducibility, we strive to use open-source or free software; major software components
and notes for accessing them are specied in Supplementary Table S2.
2.3.2 PhaseI—Detectingandlocatingearthquakes
2.3.2.1 DetectingP-waves
A short-term average to long-term average ratio (STA/LTA) detection algorithm (Allen 1982) targets P-
wave arrivals by operating on vertical-component data with sample-rate dependent parameters (Table
2.1). If the sample rate is greater than or equal to 80 s.p.s., waveform data are preprocessed with a 10 Hz
12
Phase I
Detect P-waves
dbdetect
Phase II
Phase III
Associate detections with
events
dbgrassoc
Detect S-waves
dbshear
Locate earthquakes
individually
NonLinLoc
Cross-correlate waveforms
ddcc
Double-difference relocate
events
GrowClust
Calculate magnitudes
dbevproc, dbmw
Quality control
Figure 2.3: Schematic diagram outlining the automated processing procedure. Steps are divided into three
phases: Phase I—Detecting and locating earthquakes; Phase II—Relocating earthquakes; and Phase III—
Calculating event sizes and post-processing. Each step gives a pithy description of its purpose, and the
name of software components used are provided in bold. Components that we develop or contribute to
are indicated in green. Program names beginning with “db” come from the Antelope software package
(BRTT Inc.) and all others are open-source. All software used was obtained for free; either under an
academic license, or as open-source software.
13
Table 2.1: Sample-rate dependent processing parameters used for detecting P-waves. The extremely nar-
row window used for 80 s.p.s. data is more sensitive to the targeted short-duration, impulsive signals
from microseismicity than the wider window used for lower sample-rate data. The parameters were cho-
sen to ensure STAs were computed using at least 10 samples, however, this criteria was relaxed for the
case of 80 s.p.s. data streams, of which there were only 4 used.
80 s.p.s. < 80 s.p.s.
Preprocessing lter 10 Hz high-pass 4 Hz high-pass
Short-term average window length 0.1 s 0.25 s
Long-term average window length 10.0 s 10.0 s
high-pass lter and window lengths of 0.1 s and 10.0 s are used for the STA and LTA respectively. Data with
a lower sample rate are preprocessed with a 4 Hz high-pass lter and window lengths of 0.25 s and 10.0 s
are used for the STA and LTA respectively. The STA/LTA algorithm registers a detection when the ratio
of the two averages (each calculated as the root-mean-square (RMS) signal amplitude) exceeds an onset
threshold of 5 and remains above a secondary threshold of 2 for at least 20 s. The long-term average is held
constant from the time the ratio exceeds the onset threshold until it falls below the secondary threshold.
We chose these parameters for detecting P-waves after comparing results from two tests using dierent
parameters selected based on experience. In the rst test, we preprocessed all data uniformly, regardless
of sample rate, with a 1-10 Hz band-pass lter and used window lengths of 1.0 s and 10.0 s for the STA and
LTA respectively. In the second test, we processed data as described in the preceding paragraph. Validating
the results of both tests againstANZA_review showed that the chosen sample-rate dependent parameters
yield the greatest number of and most precise detections.
2.3.2.2 DetectingS-waves
Each detected P-wave triggers a search for the corresponding S-wave using a multi-step algorithm de-
scribed by Ross et al. (2016b). This algorithm is based on ve statistical quantities, each derived from
three-component data in a sliding window: i) the signal rectilinearity, ii) the signal incidence angle, iii)
the STA/LTA ratio, iv) the signal kurtosis, and v) the rst dierence of the kurtosis. We preprocess wave-
forms with a 3-20 Hz band-pass lter. The algorithm derives signal rectilinearity and incidence angle from
14
the covariance matrix (Jurkevics 1988) in a 3 s sliding window and uses them to construct a polarizing
lter designed to isolate energy with high rectilinearity and near-vertical incidence angle. The algorithm
then applies this lter to preprocessed, horizontal-component data and registers an initial detection using
an STA/LTA algorithm with parameters like those for detecting P-waves. A procedure based on tandem
analysis of the kurtosis and kurtosis rate (-of-change) signals then renes the initial detection. We use the
same algorithm and parameters as those presented in Ross et al. (2016a) in an earlier application to data
from the SJFZ region. The full set of used parameters is provided in the Supplementary Table S5.
2.3.2.3 Associatingphasedetectionswithevents
Phase detections are associated with events by examining the spatial coherency of the detections. Speci-
cally, the procedure checks that a minimum of 5 unique stations have detections at times that agree, within
a tolerance of 1.5 s, with the arrival times modeled by a 1D travel-time calculator. For each subset of 20
consecutive detections, the procedure searches a 3D spatial grid for a tentative source location that passes
this check. If such a source location is found, the procedure creates a unique event ID and associates it with
the coherent phase detections. If more than one location passes this initial check, the procedure associates
the arrivals with the location that minimizes the RMS of the dierence between observed and modeled
arrival times.
We chose the minimum ve-station threshold after comparing results from three tests, where we varied
the minimum station threshold between four and six. In these tests, events were tagged as genuine if their
location matched (within a prescribed space-time tolerance) an event in ANZA_review, and were tagged
as false if they did not. A ve-station threshold recovers 81% of the events in ANZA_review, with a 34%
false-detection rate, and increasing the station threshold to 6 reduces the event-recovery rate to 75%. Most
of the false detections with the adopted ve-station threshold are culled in later processing steps.
15
2.3.2.4 Locatingeventsindividually
After associating phase detections with events, the events are individually located using the algorithm
NonLinLoc (Lomax et al. 2009), which implements a probabilistic, global-search for locating earthquakes
in heterogeneous 3D velocity structure. We use our hybrid 3D velocity model and the nite dierence
method of Podvin and Lecomte (1991) for solving the Eikonal equation to calculate travel-time-lookup
tables for each station. All events are located using these travel-time-lookup tables with default algorithm
parameters for locating local events (Supplementary Table S7). We report 68% (one-sigma) condence
intervals produced byNonLinLoc as absolute location uncertainties; these are estimated assuming arrival-
time uncertainties can be modeled by a Gaussian distribution with mean 0.1 s and 0.2 s for P- and S-waves,
respectively.
2.3.3 PhaseII—Relocatingearthquakes
2.3.3.1 Cross-correlatingwaveforms
To relocate clusters of events, we measure dierential travel times for pairs of closely spaced events by
cross-correlating waveforms. Each event is cross-correlated with each of its 200 nearest neighbors that oc-
curred later in time (correlating only with later events avoids redundant calculations), and each component
of data is correlated independently for each unique station-phase pair with at least one arrival associated
with the event pair.
To cross-correlate a pair of traces, the travel time of the arrival from the earlier event (or the later event
if no arrival registered for the earlier one) is calculated by subtracting the origin time from the arrival time.
This travel time is then added to the origin time of the later event (or earlier event when appropriate) to
obtain the arrival time corresponding to zero travel-time dierence between the events. The traces are
then band-pass ltered between 2 and 15 Hz, windowed around these arrival times (0.25 s before and
0.75 s after for P-waves, 0.25 s before and 1.25 s after for S-waves) and cross-correlated with a maximum
16
oset of0.5 s. The oset that produces the maximum absolute cross-correlation coecient is recorded
as the dierential travel time (i.e. double dierence), and measurements with a coecient less than 0.6 are
discarded.
We use the Adaptable Seismic Data Format (Krischer et al. 2016) to store waveform data, instead of a
conventional format like miniSEED or SAC, for its superior scalability with respect to I/O resource con-
sumption. This is an important technical consideration for cross-correlating a data set of this size; us-
ing a methodology designed to work with miniSEED data overwhelmed the I/O resources of the high-
performance computing cluster we use, making the calculation untenable.
2.3.3.2 Relocatingwithadouble-dierencealgorithm
The dierential travel times obtained by cross-correlating waveforms are input, along with an average 1D
velocity model, to theGrowClust (Trugman and Shearer 2017) double-dierence algorithm for relocating
earthquakes. A minimum of six dierential travel time are required to relocate each event pair. Aside
from the minimum number of dierential travel times per event pair, the used parameters (Supplementary
Table S8) are as those provided by the example distributed with the source code. In an earlier stage of
this project, theHypoDD algorithm was used. It was replaced byGrowClust because limits on computer
memory and the number of tunable parameters makeHypoDD cumbersome for a data set of this volume.
2.3.4 PhaseIII—Calculatingeventsizesandpost-processing
2.3.4.1 Qualitycontrol
To obtain a catalog that maximizes the ratio of genuine to false events detected, all events not relocated
by GrowClust are validated by a ve-point quality-control test, and retained only if they satisfy all ve
criteria. Events relocated byGrowClust bypass the quality-control test because it is unlikely that a false
17
event will correlate suciently well with another event for GrowClust to relocate it. All events in the
nal catalog are either relocated byGrowClust or satisfy the following:
1. event has a horizontal error (reported byNonLinLoc) less than 3 km,
2. event has a vertical error (reported byNonLinLoc) less than 3 km,
3. event is recorded by at least one station within 10 km epicentral distance,
4. event has at least three S-wave arrival detections,
5. event has an S-to-P phase-picks ratio of at least 1/4.
2.3.4.2 Calculatinglocalmagnitude,M
L
Local magnitudes, M
L
, are calculated using Richter’s original method (Richter 1935), as implemented by
thedbevproc program and accompanying moduleMlrichter.pm of the Antelope software package. The
maximum signal amplitude of horizontal component data is measured, on integrated velocity data with
simulated Wood-Anderson instrument response, in a window that starts 10.0 s before the predicted P-wave
arrival and is twice the predicted S-minus-P time in length. The noise level is calculated as 1.414 times the
standard deviation of 10.0 s of data preceding the predicted P-wave arrival time, and maximum amplitude
measurements are discarded if they are less than 3 times this value.
2.3.4.3 Calculatingmomentmagnitudes,M
w
Moment magnitudes, M
w
, are estimated using a spectral-based method implemented by the dbmw pro-
gram (Ross et al. 2016b) and parameters listed in the Supplementary Table S10. This method estimates
scalar seismic potencies and moments from the low-frequency asymptote of observed displacement spec-
tra, and then uses the scaling relation of Hanks and Kanamori (1979) to get the moment magnitude.
18
2.3.4.4 Estimatingminimummagnitudeofcatalogcompleteness,M
C
The minimum magnitude of completeness for the catalog, M
C
, is estimated using a method based on the
work of Ogata and Katsura (1993). This involves modeling the observed frequency-magnitude distribution
using an exponentially-modied Gaussian PDF (Figure 2.4) dened by Eqn. (2.1):
f
M
(m;;;) exp

2

2 +
2

exp (m)

m

+
2

!
(2.1)
where

m(+
2
)

is a Gaussian CDF with mean

+
2

and standard deviation, and is the
decay-rate of the exponential component of the distribution. The exponential component in Eqn. (2.1)
models the classical frequency-magnitude distribution following Gutenberg-Richter statistics, and the
Gaussian CDF component acts as a “thinning operator”, representing network sensitivity as a function
of magnitude. We estimate parameters,, and by numerically maximizing the log-likelihood function
(i.e., using the maximum-likelihood estimator implemented byscipy.stats.rv_continuous.t (Jones et al.
2001)), and dene the minimum magnitude of catalog completeness as the 99
th
percentile of the Gaussian
PDF underlying .
Spatial variations in M
C
are analyzed by assigning each node of a grid with 0.1 node-spacing the value
of M
C
estimated for all events within 0.1 epicentral distance of that node; nodes with fewer than 200 events
within 0.1 are not assigned an M
C
value.
2.3.4.5 Analyzingsimilar-eventchains
To investigate the degree of similarity between events we build similar-event chains based on waveform
similarity using a simple chaining rule like the one used by Peng and Ben-Zion (2005): if event A is similar
to event B, and event B is similar to event C, then events A, B, and C form a similarity chain. Events A and
C will not necessarily satisfy the criteria used to assess similarity, but are considered indirectly similar.
19
Figure 2.4: Example t of an exponentially-modied Gaussian PDF (blue curve; Eqn. 2.1) to an observed
frequency-magnitude distribution (black X’s; data taken fromHYS_catalog_2011 for a 0.1 square centered at
33.5 , 116.5 ) using maximum-likelihood estimation. The distribution is decomposed into two components:
i) an exponential component (orange, dashed line) representing ideal Gutenberg-Richter statistics, and ii) a
Gaussian CDF component (green, dashed line) representing network sensitivity. We dene the magnitude
of completeness as the 99
th
percentile of the PDF underlying the Gaussian CDF; i.e. the event magnitude
with a 99% probability of being detected.
20
Table 2.2: SJFZ_catalog_2008-2016 detection counts at dierent steps in the processing procedure. as-
sociating, QC, and DD relocating refer to the procedures described in Sections 2.3.2.3, 2.3.4.1 and 2.3.3.2
respectively.
Before associating After associating After QC After DD
relocating
“P” detections 11,391410 2,159,717 1,683,038 1,575,515
“S” detections 4,286,019 1,666,272 1,268,476 1,195,924
Events 0 160,867 108,800 102,719
Two events are considered similar in this context if the maximum cross-correlation coecient (for either
P- or S-waves; measured as described in Section 2.3.3.1) exceeds a chosen threshold (we analyze results for
various thresholds between 0.70 and 0.95) at a minimum of ve unique stations.
2.4 Results
2.4.1 Basicfeaturesofthecatalog
The performed analysis provides a new catalog, namedSJFZ_catalog_2008-2016, of earthquake parameters
for 160,867 events. During the quality-control stage, we discard 52,067 of these events as either spurious or
poorly-constrained, leaving 108,800 high-quality locations (Figure 2.5) with 105,762 (local/Richter) mag-
nitudes, 97,709 seismic potencies (scaled and reported as moment magnitudes), 1,683,038 P-wave arrival
times, and 1,268,476 S-wave arrival times (Table 2.2). About 63% of the retained events are inside the focus
region, with 16 events for every 10 events inSCSN_catalog_2010 for the same region (Figure 2.6), and 40%
of the retained events are within the Trifurcation area (TR).
Aftershocks of the three largest events in the focus region during the examined period (red stars in
Figure 2.6) increase considerably the event detection rates, but additional spikes of seismicity rates that
are temporally uncorrelated with signicant events are also observed. Visual inspections of the results
indicate that most of these features are produced by earthquake swarms. A short-lived (7-day) swarm
3 km NNW of Cahuilla produces the largest such uncorrelated spike (early 2014), a long-lived (100-day)
21
Figure 2.5: Map of 108,800 cataloged events in and around the focus region, with color showing depth
(scale at bottom). Stars represent moderate events (M 3.5; scale at bottom left) taken from the USGS
catalog—USGS catalog data are used for moderate events because our processing procedure is tuned for
detecting microseismic events. Surface traces of vertical transects A-A
0
, B-B
0
, and C-C
0
—plotted in Figures
2.10, 2.11, and 2.12 respectively—are shown as solid green lines; the center of each transect is marked by
a blue +. Locations of Clark (CL), Hot Springs (HS), Coyote Creek (CC), and Buck Ridge (BR) faults are
annotated with white arrows. Towns of Palm Springs, Anza, and Cahuilla are marked by white circles with
black edge. The Hot Springs area (surrounding center of transect B-B
0
) is characterized by deep seismicity
(> 13 km) primarily distributed inside two elongate volumes, which trend sub-parallel to the trend of and
are displaced to the NE of the main fault trace—few moderate, and no M 4, events occurred here. The
Trifurcation area (surrounding the center of transect C-C
0
) is characterized by networks of discrete surfaces
with conjugate orientations at intermediate depths (1-13 km) and deep seismicity (> 13 km) resembling
structure in the Hot Springs area—deep seismicity is chiey located between the Clark and Buck Ridge
fault traces. The largest events inside the focus region occurred in the Trifurcation area.
22
Figure 2.6: Events detected inside the focus region as a function of time in one-year bins (top panel) and
one-week bins (bottom panel) forSCSN_catalog_2010 (grey bars) andSJFZ_catalog_2008-2016 (black bars).
Over 16 events are detected inSJFZ_catalog_2008-2016 for every 10 events inSCSN_catalog_2010. Moderate
events (M 3.5; taken from the USGS catalog) are shown as white stars in the bottom panel (scaled relative
to the right vertical axis). The three largest events to occur inside the focus region (2010 M
w
5.4 Borrego
Springs, 2013 M
L
4.7 Anza Borrego, and 2016 M
w
5.2 Borrego Springs; shown as red stars) are observed in
the bottom panel as a brief spike in event detection rates followed by a decay pattern typical of aftershock
sequences. Other spikes that are temporally uncorrelated with signicant mainshocks (e.g. early 2014 and
early 2015) are visually inspected and conrmed to be swarm activity.
swarm3 km SW of the Coyote Creek (CC) fault produces the second to largest such spike (early 2015),
and additional short-lived swarms are observed near the Coyote Creek and Clark (CL) faults. It is also
interesting to note the variable event rates that follow the occurrence of earthquakes with M
L
4, likely
reecting variable stress conditions.
Absolute event locations are well-constrained, both horizontally and vertically, inside the focus region
(Figure 2.7, Table 2.3), but vertical uncertainties increase rapidly outside the footprint of the core network
coverage. Contrary to conventional assumption, the median vertical uncertainties inside the focus region
23
Table 2.3: Basic characteristics of SJFZ_catalog_2008-2016 in dierent geographical regions: # of events is
the total number of events in the catalog, Hor. error and Vert. error are median absolute horizontal and
vertical location errors respectively (68% condence intervals reported by NonLinLoc), and M
C
is the
magnitude of completeness estimated by the 99
th
percentile of the Gaussian component of the best-tting
exponentially-modied Gaussian PDF.
Study area Focus region HS TR
# of events 108,800 68,642 14,814 43,006
Hor. error 1.4 km 1.2 km 1.3 km 1.1 km
Vert. error 1.5 km 1.1 km 1.1 km 1.1 km
M
C
0.8 0.6 0.7 0.5
are smaller than horizontal, though the opposite is true outside the focus region. The magnitudes of
detected events range from -1.8 to 5.4, and the catalog is complete above M
L
0.6 inside the focus region
of this study. The SCSN_catalog_2010 is complete above M
L
0.9 inside the focus region, so SJFZ_catalog_-
2008-2016 outperforms SCSN_catalog_2010 in the focused region. On the other hand, SCSN_catalog_2010
outperforms the derived SJFZ_catalog_2008-2016 in the outer region (Figure 2.8).
Event pairs inSJFZ_catalog_2008-2016 generally exhibit weak waveform similarity (Table 2.4); at least
6 of every 10 events in the catalog are suciently unique (in terms of waveform signature) that they
are omitted from every similarity chain. Fewer than 3% of events are considered highly similar (maximum
cross-correlation coecient0.95 at 5 or more stations), and, regardless of the cross-correlation coecient
threshold chosen, a relatively small number of similarity chains contain most of the similar events. These
similarity chains form highly-localized structures in the core fault zone surrounded by diuse clouds of
dissimilar events, are strongly correlated with spatiotemporal event clusters, and give useful information
on the principle seismogenic structures in the fault zone (Figures 2.9, 2.10, 2.11, and 2.12).
2.4.2 Structuralfeaturesinthecatalog
2.4.2.1 Generalstructuralfeatures
We divide the SJFZ into four distinct, quasi-tabular regions based on observed characteristic patterns of
seismicity (Figures 2.10, 2.11, and 2.12). Three characteristic patterns are dened for this purpose: i) seismic
24
Figure 2.7: Spatial distribution of median absolute location errors (68% condence intervals reported by
NonLinLoc). The side length of each cell is 0.025 and the cell color represents the median value for
all events within 0.05 epicentral distance of the cell’s center; cells with fewer than 50 associated events
are assigned a null value. Inside the focus region, the median horizontal and vertical errors are 1.2 km
and 1.1 km, respectively, and outside of the focus region they increase to 1.8 km and 2.2 km, respectively.
Inside the focus region—where station density is high, many events are deep, and the velocity model is
well constrained—the median vertical error is smaller than the median horizontal error. Outside the focus
region, the median vertical error is larger than the median horizontal error. The dashed rectangular outline
indicates our focus region, as in Figure 2.1.
25
Figure 2.8: Spatial distribution of minimum magnitude of catalog completeness for SJFZ_catalog_2008-
2016 (left) andHYS_catalog_2011 (right) estimated by the 99
th
percentile of the Gaussian component of the
best-tting (in the maximum-likelihood sense) exponentially-modied Gaussian PDF (Eqn. 2.1). The side
length of each cell is 0.1 and the cell color represents the magnitude of completeness for all events within
0.1 epicentral distance from the cell’s center; cells with fewer than 200 associated events are assigned a
null value. SJFZ_catalog_2008-2016 andHYS_catalog_2011 are complete above M
L
0.6 and 0.9, respectively,
inside the focus region and above M
L
1.9 and 1.5, respectively, outside the focus region (and within the
mapped area).
26
Table 2.4: Summary of similar events in SJFZ_catalog_2008-2016: CC thresh is the cross-correlation coe-
cient threshold applied for chain linking,N
sim
is the total number of events belonging to similarity chains,
N
chain
is the total number of similarity chains, N
chain5
is the number of similarity chains with at least 5
events,N
sim5
is the numbe of events belonging to similarity chains with at least 5 events, andPercent is the
percent of events in the nal catalog that belong to similarity chains. The number and length of similarity
chains increases as the cross-correlation threshold decreases, until a threshold 0f 0.85 is reached, at which
point the length of chains continues to increase, but the number of chains begins to decrease: relaxing the
linking criteria beyond the critical threshold of 0.85 causes formerly independent chains to merge.
CC thresh N
sim
N
chain
N
chain5
N
sim5
Percent
0.95 2,880 511 121 1918 2.6
0.90 8,425 870 265 6,951 7.7
0.85 14,359 920 280 12,772 13.2
0.80 20,550 849 259 19,073 18.9
0.75 27,049 682 204 25,831 24.9
0.70 33,953 488 143 33,089 31.2
Figure 2.9: Map view of events belonging to similarity chains (red points) obtained for cross-correlation
coecient thresholds of 0.9 (left) and 0.8 (right); isolated events are shown as green points. Similar events
tend to occur near the center of diuse clouds of surrounding dissimilar events and form compact lin-
eations. The relatively-compact lineations formed by similarity chains (compared to diuse surround-
ing events) may be an artifact of the double-dierence relocation method, or may be evidence of highly-
localized zones of repeated failure accommodating signicant deformation.
27
Figure 2.10: Vertical transects (see also Figures 2.11 and 2.12) along A-A
0
(Figure 2.5) of all seismicity within
15 km of the nominal fault plane collapsed onto the plane (by orthogonal projection). In the top panel,
events are color-coded by origin time (scale at top), signicant events (M 4; taken from USGS catalog)
are shown as stars (magnitude scale at the right of panel) color-coded by origin time, and boundaries of the
four mechanically-distinct regions we interpret are outlined by thin black dash-dotted lines and labeled
R-I through R-IV. Events belonging to similarity chains with cross-correlation coecient thresholds of 0.9
(middle panel) and 0.8 (bottom panel) are shown as red points and independent events are shown as green
points (bottom two panels). Vertical dashed lines indicate the extent of the Hot Springs area (HS), the
Anza Seismic Gap (ASG), and Trifurcation area (TR), as we dene them. Regions I and IV are characterized
by an absence of seismicity, Region II is characterized by dense space-time clusters of events with similar
waveforms distributed on discrete quasi-planar surfaces, and Region III is characterized in this projection
by a ribbon of spatiotemporally diuse seismicity.
28
Figure 2.11: Vertical transects as in Figure 2.10 along B-B
0
through the Hot Springs area (Figure 2.5). The
vertical projections of the Clark (CL) and Hot Springs (HS) faults’ surface traces are marked by dashed
vertical lines. There are no events M> 4 in this transect, and most seismicity along this transect is in
Region III. Region III seismicity exhibits greater waveform similarity than corresponding seismicity in TR
(Figure 2.12), and, from this perspective, appears to form limbs of volumetrically distributed seismicity
which converge at low angles near the top of the region.
29
Figure 2.12: Vertical transects as in Figure 2.10 along C-C
0
through the Trifurcation area (Figure 2.5). Verti-
cal projections of the Coyote Creek (CC), Clark (CL), and Buck Ridge (BR) faults’ surface traces are marked
by dashed vertical lines (the location of the Coyote Creek fault is extrapolated by continuing the mapped
segment of the fault to the point of intersection with the transect). Region II seismicity along this tran-
sect reveals numerous discrete surfaces that were active during brief periods of intense seismic activity—
quintessential event clusters—and these event clusters exhibit high waveform similarity (middle and bot-
tom panels). The lateral extent of Region II seismicity narrows with increasing depth and converges to a
9 km wide zone of diuse Region III seismicity at13 km depth.
30
Figure 2.13: P-wave waveforms (top two panels) and event location maps (bottom two panels) for the
largest similarity chains obtained with cross-correlation coecient thresholds of 0.9 (left two panels) and
0.8 (right two panels). Similar P-wave waveforms observed at station PB.B082 (inverted blue triangle
on location maps) are cross-correlated and aligned with a template event trace (marked by a red star
on location maps; trace index 0 on waveform plots), and then sorted in order of decreasing similarity
(similarity with template decreases as trace index increases in waveform plots). Event markers in location
maps are color-coded by maximum cross-correlation coecient with respect to the template trace. A
modest reduction of the similarity-chain cross-correlation threshold from 0.9 to 0.8 results in a continuous
chain of similar events linking events that are separated by>10 km epicentral distance and span all three
major fault traces. Despite the continuity of the similarity chain in this instance, waveforms from event
pairs with large inter-event separation are only moderately similar (maximum cross-correlation coecient
0:5).
31
Figure 2.14: Vertical transect along A-A
0
(Figure 2.5) with earthquakes plotted as magenta points and
the FANG16 S-wave velocity model as background. Seismogenic structures are more fully illuminated in
SJFZ_catalog_2008-2016 than inHYS_catalog_2011 (Supplementary Figure S7), and the depth extent of the
seismogenic zone between 20 and 40 km horizontal oset is greater in SJFZ_catalog_2008-2016 than in
HYS_catalog_2011.
quiescence, ii) spatiotemporally localized, discrete manifolds of seismicity, and iii) spatiotemporally diuse
clouds of seismicity.
Region I consists of the uppermost1 km of crust and is nearly devoid of seismicity, although a note-
worthy cluster of events is observed near the SE end of the Trifurcation area. This shallow cluster is
correlated spatially with a strong bilateral fault-parallel S-wave-velocity gradient in FANG16-S (Figure
2.14) and the SW slope of the Santa Rosa Mountains.
Region II extends from the bottom of Region I to13-15 km depth—the lower bound of this region
varies along strike—and contains 62% of the observed seismicity in and around the SJFZ (an 80-km long,
30-km wide rectangular swath oriented parallel to the main fault trace plotted in Figure 2.10). Seismicity in
this region clusters tightly in space and time, forming angular networks of discrete surfaces with conjugate
32
Figure 2.15: Map of seismicity in the Trifurcation area showing intermediate-depth “shatter networks”—
networks of quasi-planar fracture surfaces with conjugate orientations. Shatter networks are distributed
sub-orthogonally to the trend of the main fault trace and produce most of the similar waveforms observed.
orientations along with signicant gaps of seismicity (Figures 2.10, 2.12, and 2.15). All seven M4 events
that occurred during the observation period nucleated in this region.
Region III populates a3-6 km ribbon beneath Region II with abundant, spatiotemporally diuse, and
low-magnitude earthquakes (37% of observed seismicity in the fault zone, no M4 events, and only four
3M4 events). The downwards bulge near -15 km horizontal oset in Figure 2.10 is correlated somewhat
with Thomas Mountain and the associated topographic ridge (Figure S7). Basement lithology may control
this feature (Magistrale and Zhou 1996), and it may also be related to the dip of the Moho in this region
(Lewis et al. 2000; Ozakin and Ben-Zion 2015; Zhu and Kanamori 2000).
Region IV includes everything beneath Region III and is almost completely devoid of seismicity.
33
The lower bound of the seismogenic zone generally shallows toward the SE, with the deepest seis-
micity observed in the Hot Springs area (HS). Separating the deep seismicity of HS to the NW from the
generally shallower seismicity in the highly-productive TR area to the SE is the Anza Seismic Gap (ASG) of
microseismicity (Sanders and Kanamori 1984). The ASG area is associated with a highly-localized section
of the SJFZ with paleoseismic evidence of moderate and large events (Rockwell et al. 2015), a high-velocity
contrast across the fault (Allam and Ben-Zion 2012; Share et al. 2018), and a strong overlying velocity
gradient in FANG16-S (Figure 2.14).
2.4.2.2 HotSpringsarea
The Hot Springs area is predominantly populated by seismicity in Region III (Figure 2.11); 78% of HS
seismicity is in Region III. The deepest earthquakes observed in the SJFZ, and the only M3.5 earthquakes
(four events) observed in Region III, are in HS. The earthquakes in HS Region III include anastomosing
NE-dipping limbs of diusely distributed seismicity that converge at low angles near the top of the region
(Figure 2.11, Supplementary Movie S1) where seismicity rates quickly drop. Region III seismicity in HS
exhibits stronger waveform similarity than corresponding Region III seismicity in TR (Figures 2.10, 2.11,
and 2.12). Region II seismicity is sparse in HS and is mostly located10 km SW of the main fault segments.
2.4.2.3 Trifurcationarea
The Trifurcation area (Figure 2.12) is the most seismically active portion of the SJFZ, with abundant
seismicity in Regions II and III; 66% of core fault zone seismicity is in TR. Seismicity in TR is asym-
metrically distributed between regions II (75%) and III (24%) and depicts “shatter networks”—principal
sets of commonly-oriented discrete fracture surfaces co-existing with sub-orthogonal conjugate surfaces—
distributed along a trend sub-orthogonal to the main fault trace (Figure 2.15). Region II seismicity in TR
includes spatiotemporally-dense event clusters, comprised of events with highly-similar waveforms, that
34
are dipping to the NE with angles that become shallower below about 10 km. Region III seismicity in TR,
however, is laterally distributed across a9 km swath below the three major faults in the area and exhibits
greater waveform diversity than the overlying Region II seismicity.
Certain events occurring on conjugate surfaces in TR produce highly-similar waveforms (cross-correlation
threshold0.9), and similarity chains with a threshold of 0.8 linking events spanning>10 km and all three
major TR faults are observed (Figure 2.13). Waveforms produced by events from such chains with large
separation, however, show only moderate similarity (maximum cross-correlation coecient of0.5), and
the coherent energy of these coincides with the highest-amplitude signal. Coherent, high-amplitude sig-
nals like this may provide information on regional processes that unify events belonging to similarity
chains.
2.4.2.4 AnzaSeismicGap
The Anza Seismic Gap (Figure 2.10), as considered here, is a region nearly devoid of seismic activity be-
tween TR and HS. It extends9.5 km in the along-fault direction in Region II, spans the entire seismogenic
zone, and narrows to4 km at the base of Region III. Its SE boundary is vertical and is met with intense
intermediate-depth TR activity while its NW boundary dips steeply to the SE and is met with modest
shallow-to-intermediate-depth HS activity. Region III seismicity rates are signicantly reduced outside
ASG on either side.
2.5 Discussion
2.5.1 Characteristicsofthecatalog
The developed new catalog complements an ongoing suite of investigations targeting the SJFZ: It provides
an enhanced regional seismological framework to contextualize highly-detailed studies of internal SJFZ
structures (Qiu et al. 2017; Share et al. 2017; Qin et al. 2018); it provides an improved data set for repeating
35
tomographic inversions at higher resolution (Allam and Ben-Zion 2012; Allam et al. 2014a; Fang et al.
2016); and provides a more complete catalog for testing dierent mechanical models of the deforming
lithosphere (e.g., Ben-Zion and Lyakhovsky 2006; Wdowinski 2009; Cheng et al. 2018). The resolving
power and coverage of the derived results bridge the gap between the internal fault-zone studies at specic
sites, and the lower-resolution high-coverage tomographic inversions. At present,SJFZ_catalog_2008-2016
is the most complete set of observed seismic activity in the SJFZ for 2008-2016, and continuing operation
of key stations in the SJFZ Experiment Network will allow the catalog to be extended beyond 2016.
We compare SJFZ_catalog_2008-2016 with the currently best regional catalog, HYS_catalog_2011, for
Southern California. HYS_catalog_2011 is derived fromSCSN_catalog_2010, so the catalog derived here has
the same relative increase in the number of events (about 60% more events thanHYS_catalog_2011). Signif-
icant along-strike structural features of the two catalogs are generally similar (Figure 2.14 and Supplemen-
tary Figure S3); however, their expression inSJFZ_catalog_2008-2016 is fuller. Despite general similarities,
seismicity shallows toward the SE more rapidly in HYS_catalog_2011. This bears important implications
on the depth of the seismogenic zone, and thus maximum likely rupture dimensions in the SJFZ. Our re-
sults are likely to be more reliable in our focus region (rectangle in Figure 2.1) for two reasons: We locate
events using i) increased station density, and ii) a higher resolution velocity model. SJFZ_catalog_2008-2016
provides a better framework than HYS_catalog_2011 for most studies inside the SJFZ; however, HYS_cat-
alog_2011 remains the standard regional catalog around the SJFZ. Additionally, we recommend using the
HYS_catalog_2011 if results are needed prior to 2008, or if high-accuracy (i.e. analyst-reviewed) arrival
times are needed.
SJFZ_catalog_2008-2016 is a fully-automated catalog targeting microseismic events, and it is impractical
to manually review each phase detection or even each event. This will become increasingly relevant as time
passes; seismic networks will mature, technology will advance, and methodology will improve, resulting in
larger earthquake catalogs. AlthoughSJFZ_catalog_2008-2016 is comparable to (and in some regards better
36
than) SCSN_catalog_2010 and HYS_catalog_2011, it is a categorically-dierent catalog because of its fully-
automated nature. Despite extensive eorts to tune our processing procedure, the employed automated
detection algorithms are not as good as trained human analysts and can make imprecise or erroneous
observations. An occasional secondary phase (e.g. reection or P-to-S converted phase) may trip the
S-wave detector. Fault zone head wave arrivals at near-fault stations (e.g., Ben-Zion 1989; Allam et al.
2014b; Ross and Ben-Zion 2014) or short bursts of noise may be mislabeled as a P-wave arrival. STA/LTA
algorithms tend to pick arrivals later than the true arrival time, and this type of error is greater for emergent
arrivals than impulsive. By comparing 16,597 P- and 12,015 S-wave arrival times fromSJFZ_catalog_2008-
2016 with corresponding arrival times inANZA_review, we estimate that the average errors associated with
determining phase arrival times are 0.07 s and 0.14 s for P- and S-waves, respectively—average errors are
estimated by computing the mean arrival-time residual relative to analyst-reviewed observations. New
phase-picking techniques using machine-learning algorithms (e.g., Meier et al. 2019) may approach the
accuracy of trained human analysts and be incorporated into future automated workows of the type
used in this work.
The absolute event-location uncertainty for individual events in SJFZ_catalog_2008-2016 are likely
greater than in SCSN_catalog_2010 and its derivative HYS_catalog_2011, because of arrival-time errors.
However, the relative event locations in SJFZ_catalog_2008-2016 may be more precise than in HYS_cat-
alog_2011, because similar automated procedures were used for cross-correlating waveforms from both
catalogs, and SJFZ_catalog_2008-2016 leverages increased station density to constrain locations. Further-
more, if we treat absolute arrival-time errors as symmetrically distributed (this is an approximation that
ignores the eect of skewness), they add random errors to event locations, i.e. we would expect them to
increase the diusivity of event clusters without imparting a systematic shift to the cluster centroid. If the
cluster centroid is accurate, as expected since we use a detailed 3D velocity model, the double-dierence
procedure should (at least partially) correct random location errors by pulling events back to their relative
37
position in the cluster. If the location of an event is accurate relative to its cluster centroid, and the cluster
centroid is accurate in an absolute sense, then the event location is also accurate in an absolute sense.
The percent of events producing highly-similar waveforms in the SJFZ (Figures 2.9, 2.10, 2.11, and 2.12)
is far lower than observed for seismicity associated with highly-localized sections of large displacement
faults. Most notably, these include the creeping section of the San Andreas fault (e.g., Nadeau et al. 2003;
Rubin and Gillard 2004; McGuire and Ben-Zion 2005) and the Calaveras fault in northern California (e.g.,
Scha et al. 2002). Clusters of events with highly-similar waveforms in other locations such as the Mar-
mara region of the North Anatolian fault were assumed to reect aseismic creep (Bohnho et al. 2017). It
is possible that some of the highly-similar chains of similar events in the SJFZ are associated with some
aseismic creep, although the geometrical character of the similar-event clusters observed here are signi-
cantly more complex than in the studies mentioned above. The chains of events producing highly-similar
waveforms in the SJFZ may be associated with some forms of earthquake swarms as reported for events in
Japan (e.g., Ito 1990). Some although not all the similar-event clusters overlap with reports of earthquake
swarms in the SJFZ (e.g., Vidale and Shearer 2006).
The weak waveform similarity of events in SJFZ implies that many of the new events that we detect will
not be recovered by common template-matching techniques because it is unlikely that the waveforms from
these events will be similar to any template waveforms in existing catalogs. Common template-matching
techniques are particularly well-suited for detecting very small events with waveforms that closely resem-
ble waveforms of previously detected (template) events, but a diverse library of initial templates is needed
to maximize the ecacy of template-based processing. Diversifying the available library of templates in
SJFZ is a key contribution of SJFZ_catalog_2008-2016. The existence of continuous similarity chains that
link distant dissimilar events (Figure 2.13) suggests that a recursive adaptation to the template-matching
technique, as in Frank and Abercrombie (2018), may be well-suited for the SJFZ.
38
The automated procedure that we developed to deriveSJFZ_catalog_2008-2016 begins with raw wave-
form data and network meta-data. Because the procedure does not depend on any a priori information
about seismicity (e.g. template waveforms) in the target region, it is applicable to diverse data sets and is
particularly well-suited for detailed reconnaissance of new eld areas.
2.5.2 Characteristicsofthefaultzone
The brittle portion of a fault zone, as reected by microseismicity, is generally correlated with heat ow
and rheological properties of the rocks (e.g., Sibson 1982; Shimamoto 2003; Ben-Zion 2008). The absence
of seismicity in most of Region I is consistent with similar observations at many other locations. This
is typically interpreted as reecting velocity-strengthening frictional properties in the top few km (e.g.,
Blanpied et al. 1991; Rice 1993), but can also reect relatively thick fault zone at shallow depth (Hillers
et al. 2006) and the general low conning pressure inhibiting the shallow crust from accumulating high
elastic strain energy. There is, however, a noteworthy cluster of shallow events at the SE end of TR, which
correlates with the Clark and Buck Ridge faults, a lateral S-wave velocity gradient in FANG16-S (Figure
2.14), and the SW slope of the Santa Rosa Mountains. The increased conning pressure at shallow depths
caused by the Santa Rosa Mountains may suciently increase the brittle strength of the underlying crust
to allow detectable events to nucleate.
Region II corresponds with the geodetically-inferred locking depth along the SJFZ (Fialko 2006; Wdowin-
ski 2009). Strain energy in this region is released mostly through brittle fracture, although the region in the
Anza Seismic Gap is a clear exception, and it is interesting to note that Inbal et al. (2017) detected triggered
aseismic slip below and around the ASG. The seismicity in Region II is associated with seismic quiescence
punctuated by brief events of intense deformation on networks of quasi-planar surfaces. All three M>4
events occur within 2 km of the lower boundary of Region II.
39
Region III is where the brittle-strength and the ductile-strength of the comprising rocks are similar
and produce a brittle-to-ductile transition zone (Sibson 2002; Kohlstedt et al. 2004; Abolfathian et al. 2018).
Strain energy accumulated in Region III is released during both brittle and ductile failure events. Ac-
celerated rock healing and partial ductile failures limit the rate of brittle failure events and development
of highly-localized fracture surfaces, thus producing diuse clouds of distributed seismicity. Since some
strain energy in this region is released aseismically, only four M3 events occur in Region III (all inside
HS).
The ductile strength of rocks in Region IV is suciently lower than the brittle strength that hypocenters
do not occur here, although rupture zones and early aftershocks of large earthquakes may occur in this
region (e.g., Rolandone et al. 2004; Ben-Zion and Lyakhovsky 2006; Jiang and Lapusta 2017).
2.6 Conclusions
We integrate nine years of data (2008-2016) from ve seismic networks operating near the San Jacinto
Fault Zone, then develop and apply an automated processing procedure to derive a new catalog—called
SJFZ_catalog_2008-2016—of earthquake locations, magnitudes, and potencies from raw waveform data.
Our catalog contains 108,800 earthquakes in the magnitude range -1.8 to 5.4, each located using a detailed
3D velocity model, and is complete above M
L
0.6 inside our focus region. The obtained catalog has the
most detailed information on seismicity in the focused study region between 2008 and 2016. The events
contained in the catalog may be used as templates to detect additional events (e.g., Gibbons and Ringdal
2006; Shelly et al. 2016). This can increase the size of the catalog by a factor of 10 or more (e.g., Ross et al.
2017; Ross et al. 2019) and may be the subject of a follow-up work. The derived catalog illuminates seis-
mogenic structures in the SJFZ in unprecedented detail and includes both deeper and shallower seismicity
than in HYS_catalog_2011. The brittle-to-ductile transition zone is evident in our results as a3-6 km
ribbon of spatiotemporally-diuse seismicity at the bottom of the seismogenic zone.
40
Chapter3
Thetheoreticalbasisfortraveltimetomographyinanisotropic,linearly
elasticcontinuum
Preface
The content of this chapter is adapted from a term paper written for a course on Advanced Seismology.
Although its contents are the least novel—the key results (i.e., the elastodynamic equations of motion and
the eikonal equation) obtained are independently derived here, though well known from much earlier
work by others—it is my personal favorite (Section 3.2 in particular) because of its mathematical beauty;
it is unmuddied by the complexities of data “from the wild.” I begin with the Euler-Lagrange equations of
the Lagrangian formalism for classical eld theory and, using the elastic potential density for an isotropic,
linearly elastic continuum under innitesimal strain, derive the fundamental elastodynamic equations of
motion (EOMs) that form the theoretical basis for the traveltime tomography performed in Chapter 5. As-
suming that the physical properties of the continuum vary suciently gradually, these EOMs are shown to
reduce to the familiar homogeneous wave equation, and, further assuming that the frequency of vibrations
is suciently large, the wave equation is shown to reduce to the eikonal equation. We solve the eikonal
equation numerically using the Fast Marching Method in Chapter 4. I also lay out in this chapter some of
the fundamentals that are assumed for the traveltime tomography in Chapter 5.
41
3.1 Introduction
A primary goal of the geophysicist is to create models of the spatiotemporal distribution of physical pa-
rameters in the Earth that permit predictive statements to be made with reasonable accuracy. For example,
one must know the speed at which seismic waves travel through the Earth with some accuracy if he wishes
to predict the amount of time those waves take to travel from point A (e.g., an earthquake hypocenter) to
point B (e.g., a densely populated urban center). Of course possessing such models is only a necessary con-
dition for making predictive statements, not a sucient one. One also needs the mathematical equipment
to solve the appropriate predictive equations for the given model(s) and problem(s) at hand. In this report,
I focus on a narrow class of such Earth models—viz., models of seismic velocity structure—and I elucidate
key mathematical concepts involved in one particular method for deriving them. Inevitably, the problems
of making predictions for a given model (the forward problem) and deriving more accurate models (the
inverse problem) are inextricably linked, and I must cover the former in some detail to clarify the latter.
This report focuses on deriving models of seismic velocity structure using the method of traveltime
tomography and is constructed as follows. In Section 3.2, I derive the fundamental elastodynamic equation
of motion in an isotropic linearly elastic continuum and then show that, under the assumption of slowly
varying material properties, the equation of motion can be decomposed into a set of scalar wave equations.
In Section 3.3, I derive the eikonal equation using a Rytov ansatz for the wave equation and taking the
high-frequency limit. In Section 3.4 I review Fermat’s principle as a justication for linearized traveltime
tomography. In Section 3.5 I describe a traveltime tomography method based on ray theory and Poisson-
Voronoi subspace projections of the inverse problem (Fang et al. 2020).
42
3.2 Elastodynamic equations of motion in a linearly elastic, isotropic
continuum
In this section, I use the Euler-Lagrange equations from the Lagrangian formalism for classical eld theory
derive the fundamental equation of motion governing vibrations in an linearly elastic, isotropic continuum—
i.e., the elastodynamic equation of motion.
For a continuum with mass density = (r) the non-relativistic kinetic-energy density is
T =
1
2
_ u
i
_ u
i
(3.1)
in whichu
i
= u
i
(r) is the scalar displacement eld in thei
th
direction of a Cartesian coordinate sys-
tem. The Einstein summation convention is implied for repeated indices throughout this chapter unless
otherwise noted.
The elastic potential density is
V =
1
2
C
ijkl
"
ij
"
kl
(3.2)
in whichC
ijkl
=C
ijkl
(r) is the stiness tensor, and"
ab
="
ab
(r) is the innitesimal strain tensor:
"
ab
=
1
2
(@
b
u
a
+@
a
u
b
): (3.3)
43
The shorthand notation@

is used throughout this chapter to represent partial dierentiation with respect
to the
th
independent spacetime coordinate; greek-letter indices represent arbitrary spacetime coordi-
nates whereas roman-letter indices are restricted to representing spatial coordinates only.
For isotropic media, the stiness tensor is highly symmetric and can be written as
C
ijkl
=
ij

kl
+ (
ik

jl
+
il

jk
) (3.4)
in which
ij
is the Kronecker delta and = (r) and = (r) are Lamé parameters (of the rst and
second kind, respectively).
So, the Lagrangian density for an isotropic continuum is
L =TV (3.5)
=
1
2
_ u
i
_ u
i

1
2
C
ijkl
"
ij
"
kl
(3.6)
=
1
2
_ u
i
_ u
i

1
8
C
ijkl
(@
j
u
i
+@
i
u
j
) (@
l
u
k
+@
k
u
l
); (3.7)
and using the symmetry of the stiness tensor (3.4), this becomes
L =
1
2
_ u
i
_ u
i

1
2
@
i
u
i
@
k
u
k

1
8
(@
l
u
k
+@
k
u
l
) (@
l
u
k
+@
k
u
l
)

1
8
(@
k
u
l
+@
l
u
k
) (@
l
u
k
+@
k
u
l
)
(3.8)
=
1
2
_ u
i
_ u
i

1
2
@
i
u
i
@
j
u
j

1
4
(@
i
u
j
@
i
u
j
+ 2@
i
u
j
@
j
u
i
+@
j
u
i
@
j
u
i
)
(3.9)
44
The Euler-Lagrange equations for a scalar eld' are
@L
@'
=@

@L
@ (@

')

: (3.10)
The displacement eld in thei
th
direction thus obeys
@L
@u
i
=@

@L
@ (@

u
i
)

: (3.11)
Now, I establish a few lemmas that will facilitate expanding the Euler-Lagrange equations.
Lemma3.2.1.
@
k

@
@ (@
k
u
l
)

1
2
@
i
u
i
@
j
u
j

= (@
l
+@
l
)@
i
u
i
(3.12)
Proof.
@
k

@
@ (@
k
u
l
)

1
2
@
i
u
i
@
j
u
j

=
1
2
@
k

@
@ (@
k
u
l
)
(@
i
u
i
@
j
u
j
)

(3.13)
=@
k
(
kl
@
i
u
i
) (3.14)
=@
l
(@
i
u
i
) (3.15)
= (@
l
+@
l
)@
i
u
i
(3.16)
45
Lemma3.2.2.
@
k

@
@ (@
k
u
l
)

1
4
@
i
u
j
@
i
u
j

=
1
2
(@
k
+@
k
)@
k
u
l
(3.17)
Proof.
@
k

@
@ (@
k
u
l
)

1
4
@
i
u
j
@
i
u
j

=
1
4
@
k

@
@ (@
k
u
l
)
(@
i
u
j
@
i
u
j
)

(3.18)
=
1
4
@
k
(2@
k
u
l
) (3.19)
=
1
2
(@
k
+@
k
)@
k
u
l
(3.20)
Lemma3.2.3.
@
k

@
@ (@
k
u
l
)

1
2
@
i
u
j
@
j
u
i

= (@
k
+@
k
)@
l
u
k
(3.21)
Proof.
@
k

@
@ (@
k
u
l
)

1
2
@
i
u
j
@
j
u
i

=
1
2
@
k

@
@ (@
k
u
l
)
(@
i
u
j
@
j
u
i
)

(3.22)
=
1
2
@
k
(2@
l
u
k
) (3.23)
= (@
k
+@
k
)@
l
u
k
(3.24)
Substituting the Equation (3.9) into Equation (3.11) and using Lemmas (3.2.1), (3.2.2), and (3.2.3) to
simplify, it follows that
46
u
l
= (@
l
+@
l
)@
i
u
i
+ (@
k
+@
k
)@
k
u
l
+ (@
k
+@
k
)@
l
u
k
; (3.25)
or,
u
l
=@
l
@
i
u
i
+@
k
(@
k
u
l
+@
l
u
k
) + (@
l
)@
i
u
i
+@
k
(@
k
u
l
+@
l
u
k
): (3.26)
Equation (3.26) is the general form of the fundamental elastodynamic equation of motion for a linearly
elastic isotropic continuum in index notation.
Next, I will write Equation (3.26) in vector notation. First, I rewrite it as
u
l
= ( + 2)@
l
@
i
u
i
@
k
(@
l
u
k
@
k
u
l
) + (@
l
)@
i
u
i
+@
k
(@
k
u
l
+@
l
u
k
): (3.27)
Then, using
@
i
u
i
=r u; (3.28)
@
l
' = (r')
l
; (3.29)
47
@
k
(@
l
u
k
@
k
u
l
) = (r (r u))
l
; (3.30)
and
@
k
u
l
= (ru)
kl
; (3.31)
Equation (3.27) can be written
u = ( + 2)r (r u)r (r u) + (r u)r +

ru + (ru)
T

r: (3.32)
Equation (3.32) is the vector-notation equivalent of Equation (3.26). Taking the one-dimensional Fourier
transform of Equation (3.32) along the time axis, each time derivative on the left-hand side transforms to
multiplication byi!, and letting
2
=
+2

and
2
=

gives
u =

2
!
2
r (r u) +

2
!
2
r (r u) (r u)
r
!
2

ru + (ru)
T

r
!
2
: (3.33)
If we assume that spatial gradients of the Lamé parameters and are small at the frequency of interest—
i.e.,
r
!
;
r
!
;—then we can neglect the last two terms and write
48
u

2
!
2
r (r u) +

2
!
2
r (r u); (3.34)
or, taking the inverse Fourier transform along the time axis,
u
2
r (r u)
2
r (r u): (3.35)
This assumption simplies the elastodynamic equation of motion at the expense of being able to accurately
describe motion in a medium with material discontinuities (for whichr;r!1).
Now I will show that Equation (3.35) describes a decoupled pair of components of motion—viz., com-
pressional and shear motion—each of which can be described by the scalar wave equation. First, I dene
new variables for the compressional strain =r u and shear strain
=r u. With this change of
variables, Equation (3.35) becomes
u
2
r
2
r
: (3.36)
Taking the divergence of Equation (3.36) yields

2
r
2
; (3.37)
whereas taking the curl yields
49

2
rr
(3.38)
=
2
r
2

: (3.39)
Equation (??) is the scalar wave equation and Equation (3.39) is three scalar wave equations written in
vector notation.
In this section I showed how the elastodynamic equation of motion for a isotropic linearly elastic
continuum results from using a non-relativistic kinetic-energy density and the appropriate elastic potential
density in the Lagrangian formalism for classical eld theory. Then I showed that, under the assumption
of slowly varying (in space) material properties, the elastodynamic equation of motion decouples into a set
of four scalar wave equations describing two components of motion (one describing compressional strain
and three describing shear strain). In the next section, I will start with the scalar wave equation to derive
the eikonal equation under a high-frequency approximation.
3.3 Theeikonalequation
In this section, I derive the eikonal equation—a relationship between traveltimes and wavespeed based on
a high-frequency approximation of the wave equation—following the derivation in Rawlinson et al. (2008)
as presented in White et al. (2020).
The homogeneous scalar wave equation can be written as
r
2
(r;t) =
1
v
2
(r)
@
2
(r;t)
@t
2
; (3.40)
a second-order linear partial dierential equation in space and time with general solutions of the form
50
(r;t) =A (r)e
i!((r)t)
; (3.41)
where , r,t,v,A,!, and represent the wavefunction, position vector, time, wave velocity, wave ampli-
tude, angular frequency, and the spatially dependent traveltime, respectively.
Substituting (3.41) into (3.40) gives

A!
2
v
2
=

r
2
A!
2
Ajrj
2

i!

2rAr +Ar
2

;
(3.42)
which, after factoring out and separating the real and imaginary parts, yields a pair of equations:
r
2
A
A
!
2
jrj
2
=
!
2
v
2
; (3.43)
and
2rAr +Ar
2
= 0: (3.44)
Dividing (3.43) by!
2
and taking the high-frequency limit!!1 gives the eikonal equation:
jrj
2
=
1
v
2
(3.45)
The high-frequency approximation made to obtain (3.45) is conventionally held to be valid when the wave-
length of the propagating wave is much shorter than the scale of velocity structure heterogeneity, although
discontinuities are permitted (Rawlinson et al. 2008).
51
Expression (3.44) is called the transport equation and can be used to compute the amplitude of the
propagating wave as a function of space if the traveltime eld is known (e.g., Buske and Kastner (2004)
and Qian and Symes (2002a)).
Sethian (1996) introduced the Fast Marching Method, an unconditionally stable nite-dierence solver
for the eikonal equation, and White et al. (2020) introduced PyKonal, a convenient Python package imple-
menting the Fast Marching Method in two and three dimensions for spherical and Cartesian coordinates.
3.4 Fermat’sprincipleandtraveltimetomography
In this section I briey review Fermat’s principle as a justication for traveltime tomography.
Fermat’s principle is a specialization of Hamilton’s principle of least action which states that the evo-
lution of a mechanical system satises
S = 0 (3.46)
in whichS represents the variation of the action,S, with respect to variation of the system trajectory
through conguration space. The action is dened as
S
Z
t
2
t
1
Ldt; (3.47)
in whichL is the system’s Lagrangian, andt
1
andt
2
are xed endpoints in time. Fermat’s principle can
be written as
52
= 0 (3.48)
in which

Z
B
A
ds
v
; (3.49)
A andB are xed spatial endpoints, andds is an element of path length betweenA andB; andv are
phase traveltime and velocity, respectively, as above. The variation is the change in the traveltime with
respect to a change in the path taken between xed endpointsA andB.
In other words, Fermat’s principle dictates that the path a seismic wave travels between pointsA (e.g.,
the source location) and B (e.g., the receiver location) is stationary with respect to time. Thus, to rst
order, small perturbations of the path from the true path have no eect on the traveltime and we are free
to make small adjustments to the velocity structure along the raypath (which in turn perturbs the raypath
by a small but negligible amount) to obtain better agreement between our velocity model and observed
data.
3.5 Fundamentalsoftraveltimetomography
In this section I review the inverse theory behind traveltime tomography and develop the method of
Poisson-Voronoi parameterization following Fang et al. (2020).
The fundamental equation in linear geophysical inverse theory is
53
Gm = d (3.50)
in which G is anNM matrix dening a linear relationship between the model parameters and the data
predictions, m is anM 1 column vector of model parameters, and d is anN 1 column vector of ob-
served data. Because errors in the observed data prevent an exact solution to Equation (3.50), approximate
solutions of the form
m
0
= arg min
m
kGm dk (3.51)
are often sought in whichkk is typically taken to be the`
2
(or Euclidean) norm. In practice, Equation (3.51)
is virtually always poorly conditioned—i.e., small perturbations of d can produce large perturbations in
~ m—because the data sampling is heterogeneous. One way to stabilize the inverse problem is to reformulate
it as
m
0
= arg min
m

kGm dk
2
+kLmk
2
+kImk
2

(3.52)
in which and are scale factors, L is a smoothing matrix, and I is the identity. The additional terms in
Equation (3.52) are smoothing and regularization constraints imposed to obtain smooth “minimum-norm”
solutions.
Alternatively, smoothing and regularization constraints can be applied implicitly by formulating the
inverse problem using a Poisson-Voronoi parameterization as follows:
54
1. Randomly sampleK points from a probably density function (PDF) dened over the spatial domain
of the model.
2. Dene a matrix P2 R
MK
such that elementP
ij
is one when the point at which thei
th
model
parameter is discretized is closest to thej
th
randomly sampled point and zero otherwise.
3. Compute = GP, the sensitivity matrix for a reparameterization of the model usingK Voronoi
cells dened by the randomly sampled points in Step 2.
4. Solve the inverse problem
~ m
0
= arg min
~ m
k ~ m dk (3.53)
in which ~ m is the reparameterized model.
5. Compute m
0
= P ~ m
0
, the model parameters on the original regular grid.
6. Repeat Steps 1–5n times to obtain a set ofn such models.
7. Compute a nal aggregate model by averaging each model parameter over the set ofn models.
The sensitivity matrix, , in this reformulation is better conditioned than G, and the resultant ~ m
0
are
thus stabilized without introducing subjective smoothing and regularization constraints as in Equation
(3.52). Smoothing and regularization are instead introduced implicitly by the averaging and model repa-
rameterization, respectively. The only subjective parameters involved in this approach are the number of
Voronoi cellsK and the number of models to average overn.
In traveltime tomography based on ray theory, the data vector d is taken to be the traveltimes of
observed phase arrivals (computed as the time dierence between observed phase arrival time and earth-
quake origin time assuming the latter is known), the model parameters m are the slowness values at points
55
on a regular grid, and the sensitivity of thej
th
data observation to thei
th
model parameter,G
ij
, is the arc
length of thej
th
ray path through the voxel whose center is thei
th
node of the grid on which the model
was discretized. Ray paths and thus arc lengths can be computed numerically using an approach based on
solving the eikonal equation as in White et al. (2020).
Conclusion
In this chapter I lay out the foundation for seismic traveltime tomograpy based on ray theory. I derive
the equation of motion for an isotropic and linearly elastic continuum and then contextualize the approx-
imation of this equation using a set of decoupled scalar wave equations. I then make a high-frequency
approximation to obtain the eikonal equation as a further simplication of the relationship between seis-
mic wavespeed and traveltimes. I argue that Fermat’s principle provides a basis for linearized inverse
theory in the case of traveltime tomography and then briey illustrate an approach for solving such an
inverse problem. This chapter provides the physical rationale supporting research currently being done to
infer material properties of the subsurface in the region surrounding theM
w
6:4 andM
w
7:1 earthquake
pair that occurred near Ridgecrest, CA in July 2019.
56
Chapter4
Solvingtheeikonalequation
Preface
The content of this chapter comes from an article published in Seismological Research Letters in 2020
under the title PyKonal: A Python Package for Solving the Eikonal Equation in Spherical and Cartesian
CoordinatesUsingtheFastMarchingMethod (White et al. 2020). Having derived the eikonal equation from
(Lagrangian) rst principles in Chapter 3, we now solve it numerically using the Fast Marching Method
(FMM). Our FMM implementation is designed to fulll a need within the seismological community for a
freely available, accurate, stable, fast, general, extensible, and easy-to-use 3D traveltime calculator and ray
tracer. It is used as a core computational engine in Chapters 5 and 6.
57
Abstract
This article introduces PyKonal: a new open-source Python package for computing traveltimes and trac-
ing raypaths in 2D or 3D heterogeneous media using the Fast Marching Method for solving the Eikonal
equation in spherical and Cartesian coordinates. Compiled with the Cython compiler framework,PyKonal
oers a Python application program interface (API) with execution speeds comparable to C or Fortran
codes. Designed to be accurate, stable, fast, general, extensible, and easy to use, PyKonal oers low- and
high-level API functions for full control and convenience, respectively. A scale-independent implementa-
tion allows problems to be solved at micro, local, regional, and global scales, and precision can be improved
over existing open-source codes by combining dierent coordinate systems. The resulting code makes
state-of-the-art computational capabilities accessible to novice programmers and is ecient enough for
modern research problems in seismology.
4.1 Introduction
How long do seismic waves generated at one point take to reach another? What path does the energy
take to get there? These basic questions are key for seismologists’ ability to locate earthquakes, image
subsurface structure, and pursue many other fundamental studies. The vast literature documenting dif-
ferent approaches to answering these questions dwarfs the number of corresponding practical tools freely
available to seismologists. Here, we present PyKonal: a new open-source Python package for computing
traveltimes and tracing raypaths using the Fast Marching Method (FMM) for solving the eikonal equation.
Most numerical approaches for computing traveltimes and raypaths can be classied into one of two
main categories: ray-based methods, which solve the kinematic ray equations; and grid-based methods,
which usually solve the eikonal equation using nite dierences or use Dijkstra-like network algorithms
58
(Dijkstra 1959) to determine the shortest path between two points. Ray-based methods compute travel-
times along individual raypaths and provide no explicit information about traveltimes along alternative
raypaths. They are favorable for their speed when the number of raypaths that the user is interested in
is small, but they suer from inaccuracy and instability in regions with signicant velocity heterogeneity.
Grid-based methods, on the other hand, compute the entire traveltime eld—i.e. from the source to every
point on a grid—and then use the gradient of the traveltime eld to determine individual raypaths. While
less ecient for determining individual raypaths, they are preferable for their stability and accuracy in
complex media and are more ecient when traveltimes are needed at a large number of points.
The kinematic ray equations can be formulated as an initial-value problem, which can be solved using
shooting methods (e.g., White 1989), or as a boundary-value problem, which can be solved using bending
methods (e.g., Julian and Gubbins 1977; Thurber and Ellsworth 1980; Um and Thurber 1987). See Rawlinson
et al. (2008) for a thorough review of ray- and grid-based methods.
A plethora of grid-based methods for solving the eikonal equation have been proposed since the late
1980s, each with dierent computational complexity and stability properties, originating primarily in the
elds of computational mathematics and physics. Most grid-based eikonal solvers implement the concept
of a narrow-band computational front (a limited region of the computational domain where information
is strategically propagated from regions where the solution is known to regions where it is unknown).
Sweepingmethods are one class of methods, stemming from the work of Zhao (2004), that omit the concept
of a narrow band, which, despite their computational eciency, nd little use in seismological applications
because of their instability in heterogeneous media.
Grid-based eikonal solvers implementing a narrow band generally use either an expanding box or the
expanding wavefront as the computational front. Reshef and Koslo (1986) were probably the rst to
propose a grid-based method for solving the eikonal equation in seismological applications, and Vidale
(1988) contributed signicantly to the expanding-box formulation. The most severe limitation of Vidale’s
59
method is that it breaches the causality of the eikonal equation under certain conditions and is thus un-
stable. A number of eorts have been directed at resolving this instability and improving the eciency
and generality of that expanding-box formulation (e.g., Vidale 1990; van Trier and Symes 1991; Podvin
and Lecomte 1991; Schneider et al. 1992; Hole and Zelt 1995; Afnimar and Koketsu 2000). High-order
essentially non-oscillatory (ENO) nite-dierence schemes (Harten and Osher 1987) and their weighted
extensions (WENO; Liu et al. 1994) have been implemented in the expanding-box framework for seismo-
logical applications (e.g., Kim and Cook 1999; Qian and Symes 2002a; Qian and Symes 2002b; Buske and
Kastner 2004); however the post-processing proposed to overcome the instability inherent in expanding-
box methods (e.g., Kim and Cook 1999) increases algorithmic complexity.
Qin et al. (1992) were apparently the rst to replace the expanding box with a computational front that
conforms approximately to the shape of the advancing wavefront. This expanding-wavefront formulation
always honors the causality of the eikonal equation, and is thus unconditionally stable, but does so at
the cost of increased computational expense because a sorted list of narrow-band values must always be
maintained. In an independent eort, Sethian (1996) combined the stability of the expanding-wavefront
formulation with the computational eciency of heap-sort technology to develop the FMM, which simul-
taneously oers unconditional stability and computational eciency. Sethian and Popovici (1999) intro-
duced the FMM to the seismological community and Rawlinson and Sambridge (2004a) generalized the
method for multiple reection and transmission phases.
While previous work established a solid practical foundation for estimating traveltimes and raypaths
used to advance the understanding of the Earth’s interior, improving precision over and providing greater
exibility than existing software can lead to better results. We implement the FMM here because it is
relatively easy to code (reducing the likelihood of bugs), unconditionally stable, and computationally e-
cient. The presented new Python package,PyKonal, can improve on various seismic imaging applications.
60
PyKonal can solve the eikonal equation in two or three dimensions using Cartesian or spherical coordi-
nates and combinations thereof. With possible applications to locating earthquakes, performing seismic
tomography, computing raypath parameters, and Kirchho depth migration, PyKonal is designed to be
accurate, stable, fast, general, extensible, and easy to use.
In the following sections we outline the theory and implemented methodology (Section 4.2), present
the performance of the implementation with respect to the design criteria above (Section 4.3), and com-
pare with existing tools and discuss potential extensions and use cases (Section 4.4). Concluding remarks,
including discussion of ongoing work, are nally oered (Section 4.5).
4.2 Method
In this section, we review the derivation of the eikonal equation following Rawlinson et al. (2008), state
the problem solved by the FMM, and describe how we implement it.
4.2.1 Derivingtheeikonalequation
The homogeneous scalar wave equation,
r
2
(r;t) =
1
v
2
(r)
@
2
(r;t)
@t
2
; (4.1)
is a second-order linear partial dierential equation in space and time with general solutions of the form
(r;t) =A (r)e
i!((r)t)
; (4.2)
in which , r, t, v, A, !, and represent the wavefunction, position vector, time, wave velocity, wave
amplitude, angular frequency, and the spatially dependent traveltime, respectively.
Substituting (4.2) into (4.1) gives
61

r
2
A!
2
jrj
2
i!

2rAr +Ar
2

=
1
v
2
!
2
; (4.3)
which, after factoring out and separating the real and imaginary parts, yields a pair of equations:
r
2
A!
2
jrj
2
=
!
2
v
2
; (4.4)
and
2rAr +Ar
2
= 0: (4.5)
Dividing (4.4) by!
2
and taking the high-frequency limit!!1 gives the eikonal equation:
jrj
2
=
1
v
2
(4.6)
The high-frequency approximation made to obtain (4.6) is conventionally held to be valid when the wave-
length of the propagating wave is much shorter than the scale of velocity structure heterogeneity, although
discontinuities are permitted (Rawlinson et al. 2008).
Expression (4.5) is called thetransportequation and can be used to compute the amplitude of the prop-
agating wave as a function of space if the traveltime eld is known (e.g., Buske and Kastner 2004; Qian
and Symes 2002b).
4.2.2 Statementoftheproblem
Given the velocity eld,v (r), we aim to determine the traveltime eld, (r), by solving (4.6) under various
boundary conditions.
62
Consider a parametric surface inR
3
, (t), representing a zero-phase wavefront that propagates per-
pendicular to itself with velocityv assumed to be a strictly positive function of space, and let (r) be the
time it takes to reach a point with position vector r from an arbitrary boundary value, (t = 0). The
eikonal equation (4.6) approximates the relationship between (r) andv (r) at high frequencies, and (t)
is the surface dening the wavefront at timet—i.e. the set of r such that (r) =t. The evolution of (t)
can be tracked by solving (4.6) for (r), givenv (r) and appropriate boundary conditions. The problem is
to solve (4.6) for (r), givenv (r) and (t = 0).
4.2.3 Solution
Sethian (1996) developed an ecient numerical method—the FMM—for approximating a solution to (4.6),
subject to an entropy condition: the wavefront can only cross each point once. Imposing this condition
implies that the solution obtained comprises only rst arrivals; the wavefront propagates along the minimal
path in terms of time between any two points.
The next section presents the FMM as implemented herein. See the original paper by Sethian (1996)
for a detailed derivation and proof that the method produces a valid solution of the eikonal equation.
4.2.4 TheFastMarchingMethod
Evaluating (4.6) at a discrete number of points gives

(r)
ijk

2
=
1
v
2
ijk
; (4.7)
where
f
ijk
f (i
1
;j
2
;k
3
); (4.8)
63
withijk being indices belonging to the set of integers, and
1
,
2
, and
3
being discretization intervals
along the
1
,
2
, and
3
axes of a general orthogonal coordinate system spanningR
3
, respectively.
To satisfy the entropy condition imposed on the approximate solution to (4.6) obtained by the FMM,
the gradient operator must be approximated in such a way that each component of the traveltime gradient
is always non-negative. Doing so ensures that information always ows in the same direction as the
propagating rst-arrival wavefront. Dierent entropy-satisfying approximations have been proposed (e.g.,
Osher and Sethian 1988; Rouy and Tourin 1992) and, as in (Sethian and Popovici 1999), we choose the
scheme from Rouy and Tourin (1992):
jr
ijk
j
3
X
a=1
max

1
h
a
D
a
ijk
;
1
h
a
D
+a
ijk
; 0

: (4.9)
Substituting (4.9) into (4.7) gives
3
X
a=1
max

1
h
a
D
a
ijk
;
1
h
a
D
+a
ijk
; 0

2

1
v
2
ijk
0; (4.10)
where h
a
represents the scale factor along the
a
axis for the used coordinate system, and D
a
ijk
and
D
+a
ijk
represent, respectively, backward and forward nite-dierence approximations to the derivative of
atijk along the
a
axis.
Equation (4.10) has anupwind dierence structure—a causality relationship implying that the value of

ijk
is fully determined by neighbouring values,
lmn
, such that
lmn
<
ijk
. In other words, information
propagates in one direction—from lesser values of to greater—and unknown values of can be derived
from neighboring known values by solving (4.10). Thus, starting with known values
ijk
= 0, the domain
of known values can be expanded until it includes all values by iteratively solving (4.10). This amounts to
solving a quadratic equation in
ijk
after expansion. The eciency of the FMM rests on carefully choosing
64
foreach ijk in griddo

ijk
1 Appendijk to Unknown
end
foreach ijk2 (t = 0)do

ijk
0 Transferijk from Unknown to Trial
end
while length(Trial)> 0do
f
ijk index of node with smallest value of in Trial Transfer
f
ijk from Trial to Known foreach
Neighbour,lmn, of
f
ijk do

lmn
solution of Eqn (??) if

lmn
<
lmn
then

lmn

lmn
if lmn2 Unknownthen
Transferlmn from Unkown to Trial
end
end
end
end
Figure 4.1: Essential elements of the FMM algorithm for solving the eikonal equation (4.6) implemented
by PyKonal.
the order in which to update the unknown values, as presented in the following paragraph and in Figure
4.1.
Initialize three empty sets: Known, Trial, and Unknown. Assign all nodes to Unknown. Set
ijk
= 0
for all nodes on the initial wavefront and transfer them from Unknown to Trial. Observe that the upwind
structure of (4.10) implies that the node in Trial with the smallest value of cannot change (if multiple
values are equally the smallest, choose one at random—this will not aect the solution); transfer it from
Trial to Known and temporarily label it Active. Update the value of each neighbor of Active that is not
in Known by using values in Known to solve (??) and transfer any of them that are in Unknown to Trial.
Iteratively transfer the node in Trial with the smallest value of to Known and update its neighbors until
all nodes are in Known.
Any sorting algorithm can be used to determine the node inTrial with the smallest value of. Follow-
ing Sethian (1996), we use a heap-sort algorithm, which hasO(N logN) worst-case performance where
N is the total number of grid nodes. In practice, the computational eciency of the FMM is nearly linear
inN.
65
x
y
z
r

Figure 4.2: ISO coordinate-system convention adopted for spherical coordinates, where x, y, and z are
standard Cartesian coordinates, and,, and are radial, polar, and azimuthal coordinates, respectively,
with2 [0;1],2 [0;], and2 [0; 2).
4.2.5 Coordinatesystems
The method has been presented in a general coordinate system thus far, but a specic choice must be made
in practice and each coordinate system has strengths and weaknesses, which are illustrated in Section 4.3.
Choosing a particular set of coordinates is tantamount to dening the scale factors, h
n
, in (4.10). In
Cartesian coordinates, the scale factors are h
x
= h
y
= h
z
= 1 and in spherical coordinates they are
h

= 1, h

= , and h

= sin, where the International Organization for Standardization’s (ISO)
convention is adopted for spherical coordinates (Figure 4.2).
Spherical coordinates are ideal for tracking spherical wavefronts, whereas Cartesian coordinates are
ideal for tracking planar wavefronts (Section 4.3). Seismologists often treat sources of seismic energy as
point sources, making spherical coordinates suitable at small and intermediate distances. Wavefronts from
suciently distant sources are often assumed to be planar, for which Cartesian coordinates are ideal. When
wavefronts must be tracked over great distances and the sphericity of the Earth must be accounted for,
66
spherical coordinates are again naturally suited to the problem. It is often desirable to combine dierent
coordinate systems when solving a particular problem.
Implementing the FMM solution of the eikonal equation in a general coordinate system makes inte-
grating dierent systems easy. To combine dierent coordinate systems, we rst solve the eikonal equation
in one, then map the solution into a second and solve the eikonal equation in the remaining unknown re-
gion. This procedure can be iteratedadinnitum, but combining two or three coordinate systems is likely
sucient for most problems in seismology.
The unconditional stability of the FMM owes to a proper entropy-condition satisfying approximation
for the gradient and has only been proven in the case of rst-order-accurate nite-dierence operators.
Strictly speaking, transitioning between coordinate systems may destabilize the FMM by violating the
entropy condition. Transitioning from the rst coordinate system to the second assumes that the wavefront
does not propagate from the second back into the rst anywhere. This is likely of little concern in the vast
majority of practical use cases; however, users should be aware of this feature and use appropriate caution.
One further caveat users should be aware of is the singularity in the gradient operator at the origin of
spherical coordinate systems. The eikonal equation is unsolvable there because the gradient operator is
undened. Thus, PyKonal will raise an error if the user attempts to place a node at the origin of a spherical
grid. Grid nodes can be placed suciently close to the origin that this has little eect in practice, but
any wavefront propagating through the origin will be distorted. Future updates to PyKonal may include
implementing an unstructured-mesh update scheme (Sethian 1999) at the origin of spherical coordinate
systems to mitigate the singularity there.
67
4.2.6 Orderofthenite-dierenceoperators
The accuracy of the FMM depends partially on the order of the nite-dierence operators used where
generic operators are specied in (4.10). The simplest but most inaccurate scheme uses only rst-order
operators:
D

1
ijk
=D

1
1
ijk

ijk

i1jk

1
; (4.11)
D
+
1
ijk
=D
+
1
1
ijk

i+1jk

ijk

1
; (4.12)
and similar denitions in
2
and
3
, where
n
m
indicates ann-th order nite-dierence operator along the

m
axis (e.g.D
+
1
1
ijk
andD

2
3
ijk
represent the nite-dierence approximation of the derivative of computed
using a rst-order forward operator along the
1
axis and a second-order backward operator along the
3
axis, respectively). More accurate approximations can be obtained by using higher-order operators, e.g.,
second-order operators:
D

1
ijk
=D

2
1
ijk

i+2jk
4
i+1jk
+ 3
ijk
2
1
; (4.13)
D
+
1
ijk
=D
+
2
1
ijk

i2jk
4
i1jk
+ 3
ijk
2
1
; (4.14)
and similar denitions in
2
and
3
.
The causality of (4.10) implies that using higher-order operators is only permissible when the trav-
eltime eld at all nodes involved in an update is monotonically decreasing away from the node being
updated. For example, a second-order update along thex-axis using a backward nite-dierence operator
68
is only permissible when
i2jk

i1jk

ijk
. As in Sethian (1999), we implement second-order op-
erators whenever known values satisfy the causality condition and rst-order operators otherwise, which
means that the update scheme is of mixed order.
4.2.7 Raytracing
Energy propagates perpendicular to the wavefronts in isotropic media, so the gradient of the traveltime
eld is always tangent to the raypath between two points. Thus the raypath can be expressed as a para-
metric curve, R (s), satisfying
dR (s)
ds
=r; (4.15)
wheres represents distance along the curve. Equation (4.15) is a rst-order ordinary dierential equation,
which can be solved by integrating overds using the Runge-Kutta method. PyKonal implements a simple
rst-order Runge-Kutta method (Euler’s method).
Because the source location necessarily coincides with the global minimum of the traveltime eld, the
gradient of the traveltime vanishes there and we cannot integrate (4.15) from source to receiver. Instead,
we integrate (4.15) in the reverse direction, from receiver to source, which means that we nd the path of
steepest descent (in terms of traveltime) between the receiver and source.
The fact that the global minimum of the traveltime eld coincides with the source location presents
a convenient condition for stopping integration of (4.15). Because we are following the path of steepest
traveltime descent, the traveltime at each point along the raypath should be lower than at every point that
precedes it. Thus, the ray will satisfy (R (s
1
)) < (R (s
2
)) for alls
1
> s
2
. If the ray overshoots the
source location, it will start traveling “uphill” and will violate this condition. Thus, we stop integrating as
soon as the ray encounters a point associated with a traveltime value exceeding that of the immediately
preceding point.
69
4.3 Results
We now demonstrate that ve of our six aims in building this tool are achieved: The code is accurate,
stable, fast, general, and extensible. Users may judge whether our tool satises our sixth aim: easy to use.
4.3.1 Accuracy
Consider the trivial case of waves emanating from a point source withv = 1 km=s everywhere (Figure
4.3). Spherical wavefronts expand outward from the source, reaching a given point after an amount of
time that is in direct proportion to the distance from the source. Comparing this analytical solution to a
numerical solution computed using Cartesian coordinates reveals numerical anisotropy. As discussed by
Alkhalifah and Fomel (2001), errors are greatest in the near-source region—where wavefront curvature is
large relative to the node spacing—and in regions where the wavefront is oblique to the coordinate axes.
In spherical coordinates, however, the numerical solution is exact up to computational precision because
the radial axis is always orthogonal to the wavefront and the partial derivative of the traveltime eld
with respect to azimuth (and polar angle in 3D) is always zero. This simple test suggests that, even for
more complicated velocity elds, spherical coordinates centered on a point source are more accurate than
Cartesian coordinates in the near-eld region where wavefront curvature is high.
If the point source in the above example is replaced by a line (in the 2D case) or a plane (in the 3D case)
source, the opposite holds and Cartesian coordinates are more accurate (Figure 4.4). In this case, they-
axis is always orthogonal to the wavefronts leading to an accurate solution, whereas spherical coordinates
suer from numerical anisotropy.
Much work has been devoted to reducing numerical anisotropy introduced by nite-dierence schemes,
particularly in the eld of hyperbolic partial dierential equations (see Sescu (2015) for a thorough review).
Adapting such approaches to our context seems plausible, however the preceding two examples elucidate
an important feature of the FMM for solving the eikonal equation: judicious design of the computational
70
Figure 4.3: (a,b) Traveltime elds and (c,d) absolute errors for a point source at (x;y) = (0; 0) computed
in (a,c) Cartesian coordinates and (b,d) spherical coordinates. (b,d) Portions of the plots are blank because
the radial axis of the spherical grid used extends to a maximum of 25 km.
71
Figure 4.4: (a,b) Traveltime elds and (c,d) absolute errors for a line source aty = 0 computed in (a,c)
Cartesian coordinates and (b,d) spherical coordinates. (b,d) Portions of the plots are blank because the
radial axis of the spherical grid used extends to a maximum of 25 km.
72
grid can yield highly accurate solutions while eectively minimizing numerical anisotropy. Because nu-
merical anisotropy can be eectively reduced in this manner, we avoid endeavoring here to implement
complicated nite-dierence schemes to further reduce it.
4.3.2 Stability
To be useful for solving modern problems in seismology, any ray-tracing tool must be numerically stable
in highly complex velocity structures. The FMM is unconditionally stable (Sethian and Popovici 1999)—
meaning that solutions converge to a stable solution that is consistent with the exact solution in the limit
that the node intervals go to zero—making it well-suited to any problem in seismology where the eikonal
equation is valid. To demonstrate this convergent behavior, we consider a point source at the surface of
the Marmousi2 velocity model (Versteeg 1994)—a standard model in exploration seismology comprising
complex structures with lateral and vertical gradients, discontinuities, inversions, folded layers, pinch outs,
and lenses— (Figure 4.5). Solutions converge towards a stable solution as the node intervals are successively
halved. Unconditional stability results from proper choice of an entropy-satisfying gradient approximation
(Expression 4.9), and although it has only been proven in the case of rst-order accurate nite-dierence
operators, these results suggest that the algorithm remains stable with higher-order operators for complex
practical use cases.
The unconditional stability of the FMM implies that raypaths obtained by integrating (4.15) can be made
arbitrarily accurate by making the computational grid suciently dense (Figure 4.6). Errors accumulate
along the integration path, so they are least at the beginning of the integration path (near the receiver) and
greatest at the end (near the source). This error accumulation is exacerbated by relatively large traveltime
errors in the source region. Additionally, the criteria used for terminating integration permits inaccurate
raypaths to overshoot the source location; however, the amount of overshoot decreases as grid density
increases.
73
Figure 4.5: (a) The Marmousi2 velocity model. (b) Traveltime eld computed on a Cartesian grid with 6801
and 1401 nodes along thex andy axis, respectively, for a source located at the top left corner (x;y) = (0; 0).
Black dashed lines indicate wavefronts at 0.25 s intervals. (c-i) Dierence, t, at coincident nodes between
traveltime eld in (b) and traveltime eld computed on a grid decimated by a factor,d, specied in the top
left corner of each plot.
74
Figure 4.6: Rays traced through a medium with linear velocity gradient v (z) = (4:5 + 0:25z) [km=s]
using dierent grid densities. The densest grid contains 4096 and 1024 nodes in thex andy directions,
respectively, and the grid density is iteratively decimated by a factor of 2 whered represents the decimation
factor relative to the densest grid (e.g. the grid ford = 8 has
4096
8
= 512 and
1024
8
= 128 nodes in the
x andy directions, respectively). Decimation factors corresponding to each ray are given in the legend
at right. The analytical solution for the raypath is shown as a dashed black line. The black stars and
inverted black triangles represent the source and receiver locations, respectively. (a-c) Zoomed in regions
(22x magnication; axes ticks are at 0:2 km intervals) near the (a) source, (b) midpoint, and (c) receiver. (d)
The entire raypath.
75
Table 4.1: Execution time for problems of various grid sizes. Simulations were run on a Dell Optiplex 9020
workstation with a quad-core 3:4 GHz Intel Core i7-4770 CPU.
Grid size Number of nodes Execution time [s]
2 2 2 2
3
1.410
−4
4 4 4 2
6
1.110
−4
8 8 8 2
9
6.510
−4
16 16 16 2
12
2.810
−3
32 32 32 2
15
2.810
−2
64 64 64 2
18
2.510
−1
128 128 128 2
21
2.6
256 256 256 2
24
3.710
1
512 512 512 2
27
4.410
2
4.3.3 Speed
Having established the accuracy of our FMM implementation for the simplest cases (Section 4.3.1) and its
stability in the complex case of the Marmousi2 model (Section 4.3.2), we turn to its execution speed. We
consider a suite of 3D problems of various sizes with a random velocity model and a source at a corner of
the grid. The execution time scales nearly linearly with grid size (Table 4.1) making large problems (> 10
6
nodes) tractable and small problems (< 10
6
nodes) trivial.
4.3.4 Generality
In Section 4.3.2, we demonstrated that our implementation of the FMM in Cartesian coordinates is applica-
ble at local scales where the curvature of the Earth is negligible. Here we demonstrate that our implemen-
tation in spherical coordinates works well at global scales and the property of stable convergence holds.
Note that it is the exact same code base that solves the eikonal equation in both Cartesian and spherical
coordinates, eliminating the need for multiple implementations—only the scaling factors in (4.10) change
between operating modes, and this is done automatically at run time.
76
Figure 4.7: (a) A 2D slice of the MITP2008 velocity model shown as percent P-wave deviation from the
ak135 reference model (Kennett et al. 1995). (b) Traveltime eld computed on a spherical grid with 1024
and 2048 nodes along the and axis, respectively, for a source located at 1444 km depth. The white star
with black outline and black dashed lines indicate the source location and wavefronts at 20 s intervals,
respectively. (c-i) Dierence, t, at coincident nodes between traveltime eld in (b) and traveltime eld
computed on a grid decimated by a factor,d, specied in the top left corner of each plot.
77
We repeat the exercise from Figure 4.5, but this time we use a 2D slice from the Li et al. (2008) global
mantle P-wave model (MITP2008) in spherical coordinates. As in the case of Cartesian coordinates, the
traveltime eld converges to a stable solution as the node intervals decrease (Figure 4.7).
4.3.5 Extensibility
4.3.5.1 Trackingsecondaryphases
The FMM is limited to tracking rst arrivals, but secondary reected arrivals contain much useful informa-
tion for seismologists. These secondary arrivals can be tracked using the multi-stage approach developed
by Rawlinson and Sambridge (2004b) and extended to 3D spherical coordinates by de Kool et al. (2006). We
now show that our code can readily be extended to implement such a multi-stage approach for tracking
secondary arrivals.
To show this, we track waves reected from the lower boundary of a rectangular section (Figure 4.8).
First, we compute the downgoing traveltime eld from the source to all points. Then, we re-initialize
the traveltime eld and set the traveltime at each node along the bottom edge to the traveltime of the
incident downgoing wave. The upgoing reected waves are then tracked to all points. This relatively
simple example serves as a proof of concept: PyKonal can be extended to solve more complex problems than
tracking rst arrivals. Tracking multiple reections and refractions from complex undulating surfaces,
however, requires additional functionality that PyKonal does not currently oer. FM3D (Rawlinson and
Sambridge 2004b; de Kool et al. 2006) is a more exible and mature software package with respect to this
class of problems.
4.3.5.2 Locatingearthquakes
Seismic analysts locate earthquakes on a daily basis, and typically use 1D velocity models to represent
the 3D Earth structure because simple algorithms that can be incorporated into relevant workows are
78
Figure 4.8: Reected waves propagating through a velocity model with random velocity perturbations.
The white star with black outline indicates the source location. (a,b) The black dashed lines represent
wavefronts at 0.25 s intervals for (a) downgoing wavefronts and (b) wavefronts reected by the lower
boundary near 12 km depth. (c) The black lines represent raypaths reected by the lower boundary and
arriving at 5 km intervals at the surface.
79
not available. Research seismologists are left to determine more accurate locations that account for 3D
structure, which could be obtained in the rst place with the necessary tools. As a nal demonstration, we
show that it is straightforward to extend the PyKonal algorithm to locate earthquakes.
To show this, we relocate a set of 148 well-constrained events in Southern California, taken from the
catalog of White et al. (2019), using a simple grid search followed by dierential evolution optimization
(Storn and Price 1997) in a small neighborhood around the optimal grid point. Dierential evolution is
a genetic algorithm that performs a stochastic global search to optimize complex multivariate functions.
Each of the chosen events is in the focus region of the study by White et al. (2019), is associated with
at least 128 arrivals with root-mean-square residual (RMS)< 1 s, and was located using the NonLinLoc
software package (Lomax et al. 2009) for probabilistic non-linear earthquake location in 3D heterogeneous
Earth models and a 3D velocity model derived from the results of Fang et al. (2016).
To relocate events we use the source-receiver reciprocity and compute traveltime elds for each re-
ceiver by treating them as sources. We use the same velocity model as in White et al. (2019) and a spherical
grid with 64, 128, and 128 nodes in the radial, polar, and azimuthal directions, respectively. The origin time
is estimated for each node-arrival pair by subtracting the node-to-station traveltime from the observed ar-
rival time, and the node that minimizes the standard deviation of origin time estimates for all arrivals
provides an initial estimate of the hypocenter location. The average of the origin times estimated for the
best tting node provides an initial origin time value. The dierential evolution algorithm implemented
by the scipy.optimize.dierential_evolution function in the SciPy Python package (Jones et al. 2001)
then renes this initial location by searching a small region around the initial location for the value that
minimizes the RMS between observed and synthetic arrival times.
Relocated events are consistent with NonLinLoc locations (Figure ??), but have lower RMS in all 148
cases suggesting that this rudimentary algorithm consistently nds small adjustments that better t the
data. More sophisticated location algorithms can be designed using PyKonal to solve the forward problem.
80
4.4 Discussion
4.4.1 ComparingwithFM3D
PyKonal builds on the foundation laid by FM3D (de Kool et al. 2006) in four important ways: (a) it can
improve accuracy by using hybrid coordinate systems; (b) it accounts for the azimuthal periodicity inherent
in spherical coordinates to allow wavefronts to propagate across the = 0 plane, making many global-
scale problems easier to solve; (c) it solves both 2D and 3D problems; and (d) it oers a simple Python API
to accommodate diverse use cases.
Repeating the simple point-source error analysis from Section 4.3.1 for identical problems solved using
FM3D and PyKonal demonstrates the improved accuracy achieved by PyKonal (Figure 4.9). We consider
a point source in homogeneous velocity structure and compare errors associated with FM3D and two
dierent congurations of PyKonal. In all three cases, the computational grid comprises 256, 5, and 256
nodes along the radial, polar, and azimuthal axes, with 10 km, 0.125°, and 0.125° node intervals, respectively.
Both FM3D and PyKonal can signicantly reduce the error throughout the computational domain by using
a rened source grid. The rened source grid for FM3D in this test spans 20 coarse-grid nodes along
each coordinate axis and the node intervals are decreased by a factor of 10; we chose these parameters
because they provide a reasonable tradeo between execution time and accuracy. The rened source grid
for PyKonal consists of a spherical grid centered on the source with a grid renement factor of 5 and
extending 40 coarse grid nodes in the radial direction. The accuracy of FM3D in this test approaches that
of PyKonal in the region far from the source, but required over an order of magnitude more execution time
(41.80 s compared to 1.12 s) to do so on our workstation (a Dell Optiplex 9020 workstation with a quad-core
3:4 GHz Intel Core i7-4770 CPU). Although PyKonal can achieve better accuracy with less execution time
than FM3D, and is thus preferable in certain use cases, we note that FM3D oers more exibility as it can
81
track multiple reection and transmission phases as well as teleseismic phases propagating through local
3D structure.
Wavefronts that cross over the = 0 plane need to be carefully treated to solve problems at global
scales. If the periodicity across this plane is neglected, the obtained solutions will, in the best-case scenario,
only be valid in a hemisphere centered on the source. By properly accounting for this periodicity, our tool
accommodates global-scale problems without any extra eort required of the user.
Python is a popular high-level interpreted programming language that suers from slow execution as a
result of the fact that it checks data types at run time. Languages that check data types during compilation,
like C, C++, and Fortran, execute quickly but are slow for developments. We leverage the high-level
scripting capabilities of Python while maintaining C-like speeds by producing compiled C extensions for
Python using the Cython compiler framework. The code functions as any pure Python package would, but
with the speed of C. This makes the presented tool user friendly while standing up to computationally-
demanding tasks.
4.4.2 Applications
A robust ray tracer has many applications; we highlight the four uses that we had in mind while developing
this tool: locating earthquakes, body-wave tomography, computing take-o angles, and Kirchho depth
migration.
Locating earthquakes typically involves nding the set of spacetime coordinates that minimizes a mist
function between observed and synthetic arrival times. In seismically active fault zones where highly
damaged rocks are juxtaposed against one or two dierent host rocks (e.g., Ben-Zion 2003; Fang et al.
2016) and other areas where wavespeed varies signicantly, 3D heterogeneities need to be considered to
accurately locate events. In Section 4.3.5.2, we showed that developing an inversion framework to locate
earthquakes usingPyKonal to solve the forward problem is straightforward. More sophisticated algorithms
82
Figure 4.9: Absolute errors for traveltime elds computed through a homogeneous velocity model (v =
1 km=s) using three dierent methods: (a) PyKonal with a single coarse grid; (b) FM3D with a rened
source grid and a coarse far-eld grid; and (c)PyKonal with a rened source grid and coarse far-eld grid.
83
than the one presented here can be developed. One potential application of PyKonal is to incorporate it
into a probabilistic location algorithm that uses a Markov Chain Monte Carlo approach to includeapriori
information about the data and model. Such an algorithm may provide more reliable location estimates
and associated uncertainties.
Similar to locating earthquakes, body-wave tomography typically involves minimizing a mist func-
tion between observed and synthetic arrival times. In addition to the optimal spacetime coordinates of
the sources, the optimal velocity structure is sought in body-wave tomography. In active-source surveys,
where the source locations are known, only the velocity model is sought. Again, being able to accurately
synthesize traveltimes in the presence of 3D heterogeneity is crucial for deriving accurate velocity models.
Inverting rst-motion polarity data for source mechanisms depends on the take-o angles at which
observed rays leave the focal sphere, and small dierences in take-o angles can produce signicantly
dierent responses at distant observation points. Because PyKonal can reduce traveltime errors in the
near-source region (Figure 4.9), more accurate raypaths, and thus take-o angles, can be computed. More
accurate take-o angles should yield more accurate focal mechanisms and lead to better resolution of
rupture kinematics.
Kirchho depth migration (Schneider 1978) is a method for transforming seismic data from the time
domain to the depth domain that nds wide use in seismic image processing as a means of resolving
structural complexities by backpropagating scattered energy to the scatterer’s position in the image. Mi-
grating data in this way requires computing traveltime elds that account for 3D structural heterogeneity.
PyKonal can be integrated easily into imaging workows to facilitate clearer and more accurate images of
subsurface structures.
84
4.4.3 JupyterNotebookExamplesandDocumentation
To promote reproducibility and provide examples of how to use the code, we include as supplementary
material Jupyter notebooks that can be used to recreate Figures 4.3-4.8. We encourage interested readers
to begin familiarizing themselves with the tool by reproducing the gures in this article. After running
the example notebooks, users will benet from API documentation, which is being actively developed and
is available at https://malcolmw.github.io/pykonal- docs.
4.5 Conclusions
PyKonal is a new open-source Python package for computing traveltimes and tracing raypaths in 2D and
3D based on the Fast Marching Method for solving the eikonal equation (Sethian 1996). It is designed to
be easy enough for novice programmers to use and ecient enough to solve complex research problems
in seismology without sacricing generality or extensibility. The most recent release of the code can
be downloaded from https://github.com/malcolmw/pykonal/releases. Users are encouraged to submit bug
reports via GitHub or via email to the rst author at malcolm.white@usc.edu.
We are now using PyKonal in a exible tomographic method requiring minimal a priori information
based on Fang et al. (2020), and are extending that method to include a non-linear event relocation step
similar to the one in Section 4.3.5.2. We are applying the new methods to complementary data sets (Hutton
et al. 2010; White et al. 2019) to derive integrated hierarchical images of the crust in Southern California
with focus on the San Jacinto Fault Zone and region surrounding the 5 July 2019M
w
7.1 Ridgecrest earth-
quake sequence.
85
4.6 DataandResources
TheMITP2008 velocity model used in this paper is available via the cited reference (Li et al. 2008). Figures
3-10 were made using Matplotlib (Hunter 2007).
4.7 Acknowledgements
The study was supported by the U.S. Department of Energy (award DE-SC0016520). The manuscript ben-
etted from useful comments by N. Rawlinson and an anonymous referee.
86
Chapter5
TraveltimetomographyfromanautomatedcatalogofRidgecrest
aftershocks
Preface
The content of this chapter comes from an article currently under review (as of July 28, 2021) at Geophysical
Journal International under the title Detailed traveltime tomography and seismic catalog around the 2019
M
w
7.1 Ridgecrest, California, earthquake using dense rapid-response seismic data (White et al. 2021) and
forms the second main case study of this thesis. We update the automated processing procedure from
Chapter 2 to use more advanced algorithms and extend it to perform traveltime tomography based on the
theory and software presented in the preceding two chapters. We apply the revised procedure to data
from a network of instruments that were densely deployed in response to the Ridgecrest earthquake pair.
The resultant catalog comprises over 95 000 earthquakes during the rst three months of the aftershock
sequence (including one month of foreshocks), and traveltime tomography yields detailed 3D models ofV
p
,
V
s
and
Vp
=Vs structure. Jointly interpreting aftershock seismicity, velocity structure, and surface geology
leads us to three retrospective hypotheses:
(a) Rigid crust beneath the Coso Mountains formed a barrier that the northwest-propagating rupture
front of theM
w
7.1 earthquake could not penetrate, causing it to arrest.
87
(b) Compliant crust associated with the Garlock Fault dissipated stress concentrated at the southeast-
propagating rupture front of theM
w
7.1 earthquake, causing it to arrest.
(c) Isolated seismicity beneath the Fremont Valley accommodated transtensional strain induced on the
Garlock Fault by the main events in a pull-apart basin structure.
88
Abstract
We derive a detailed earthquake catalog andV
p
,V
s
, and
Vp
=Vs models for the region around the 2019M
w
6.4
and M
w
7.1 Ridgecrest, California earthquake sequence using data recorded by rapid-response, densely
deployed sensors following the Ridgecrest mainshock and the regional network. The new catalog spans a
four-month period, starting on 1 June 2019, and it includes nearly 95,000 events detected and located with
iterative updates to our velocity models. The nal models correlate well with surface geology in the top 4
km of the crust and spatial seismicity patterns at depth. Joint interpretation of the derived catalog, velocity
models, and surface geology suggests that (a) a compliant low-velocity zone near the Garlock Fault arrested
theM
w
7.1 rupture at the southeast end; (b) a sti high-velocity zone beneath the Coso Mountains acted as
a strong barrier that arrested the rupture at the northwest end; and (c) isolated seismicity on the Garlock
Fault accommodated transtensional-stepover strain triggered by the main events. The derived catalog and
velocity models can be useful for multiple future studies, including further analysis of seismicity patterns,
derivations of accurate source properties (e.g., focal mechanisms), and simulations of earthquake processes
and radiated seismic waveelds.
5.1 Introduction
On 4 July 2019, aM
w
6.4 earthquake ruptured a previously unmapped fault near the town of Ridgecrest,
California. Thirty-four hours later, aM
w
7.1 earthquake ruptured another nearby unmapped fault during
the largest earthquake to rattle southern California since the 1999 M7.1 Hector Mine earthquake. The
earthquake science community responded quickly, and eld teams led by the United States Geological
Survey (USGS) and the Southern California Earthquake Center (SCEC) deployed 480 seismic sensors to
record seismic data generated by the subsequent aftershock sequence (Steidl et al. 2019; Catchings et al.
2020; Cochran et al. 2020). The data collected by these sensors complement a wide array of additional
89
in-situ eld observations (e.g., Brandenberg et al. 2020; DuRoss et al. 2020; Floyd et al. 2020; Mattioli et al.
2020; Ponti et al. 2020) and remote-sensing observations (e.g., Donnellan et al. 2020; Fielding et al. 2020;
Hudnut et al. 2020; Jin and Fialko 2020; Magen et al. 2020; Pierce et al. 2020) of various kinds. The results
from these data highlighted the signicant structural and earthquake complexities in the area.
Multiple studies have investigated details of the Ridgecrest earthquake sequence by examining spa-
tiotemporal seismicity patterns observed in earthquake catalogs (e.g., Lee et al. 2020a; Lee et al. 2020b; Lin
2020; Liu et al. 2020; Lomax 2020; Ross et al. 2019; Shelly 2020). Previous such analyses, however, only
cover about one month of the aftershock sequence or less and do not include results based on the rapid-
response dense deployment of seismographs. Although the standard earthquake catalog for this region
(Hauksson et al. 2012) oers extensive records of the Ridgecrest sequence, it also does not use the dense
data that were collected in response to theM
w
6.4 andM
w
7.1 events. Furthermore, no study until now
has combined fully 3D procedures for locating earthquakes and updating velocity models. To improve
the information on seismicity patterns and crustal structures in the region around the Ridgecrest earth-
quake sequence, in this study, we derive a local earthquake catalog and velocity models (V
p
,V
s
, and
Vp
=Vs)
using rapid-response dense-deployment seismic data and data from the SCSN. The obtained new catalog
has roughly 95,000 earthquakes covering one month of foreshocks and three months of the aftershock
sequence and is used to perform detailed traveltime tomography analysis for seismic velocity structures
in the area.
To derive an earthquake catalog that is independent ofapriori observations, we processed raw wave-
form data recorded by 152 seismic sensors over a four-month period, starting one month before the main-
shocks, using an automated processing procedure. We then iteratively updated event locations and velocity
models using fully 3D methods based on the Fast Marching Method for solving the eikonal equation (Fang
et al. 2020; White et al. 2020). Joint interpretation of the derived seismicity andV
p
,V
s
, and
Vp
=Vs models,
along with the surface geology, suggests that (a) compliant crust associated with the Garlock Fault arrested
90
rupture of theM
w
7.1 event propagating to the southeast; (b) rigid crust beneath the Coso Mountains acted
as a strong barrier that arrested the rupture propagation to the northwest; and (c) isolated seismicity on
the Garlock Fault accommodated transtensional-stepover strain triggered by the main events.
The data used in this paper are described in Section 5.2. A general description of the processing proce-
dure is provided in Section 5.3, and more detailed technical material is presented in Appendix A. The main
results obtained are presented in Section 5.4, and we discuss our interpretation of the results in Section 5.5.
The derived seismic catalog and velocity models can be useful for multiple future studies in the region.
5.2 Data
Rapid response teams, led by the USGS and SCEC, started progressively deploying 480 seismic sensors in
the region surrounding theM
w
6.4 andM
w
7.1 earthquakes on 7 July (Steidl et al. 2019; Catchings et al.
2020; Cochran et al. 2020). The network codes assigned by the International Federation for Digital Seismo-
graph Networks (FDSN) to these deployments are 3J (Steidl et al. 2019), GS (Albuquerque Seismological
Laboratory and United States Geological Survery 1980), and ZY (California Institute of Technology and
United States Geological Survey Pasadena 1926). Data from the 3J network are accessible through the In-
corporated Research Institutions for Seismology Data Management Center (IRIS DMC), data from the ZY
network are available through the Southern California Earthquake Data Center (SCEDC; Southern Cali-
fornia Earthquake Center 2013), and data from the GS network are available through both the IRIS DMC
and the SCEDC.
The 3J network comprises 461 three-component (3C) nodal seismometers, deployed in fourteen sub-
arrays. Two sub-arrays form rectangular grids, eleven form fault-perpendicular linear arrays, and the
remaining one is a sparse array of three nodal seismometers that were co-located with broadband seis-
mometers. In this study, we use data from the two rectangular grids (Figure 5.1). Forty-seven nodal seis-
mometers, centered on the mainshock ruptures, operating between 14 July and 10 September, and deployed
91
Figure 5.1: Map of stations used in this study, color coded by FDSN network code. Data from rapid-
response deployments used in this study include fteen broadband seismometers (GS and ZY networks;
gold triangles) and seventy-eight nodal seismometers (3J network; red triangles). An additional fty-ve
permanent broadband stations (CI, NN, and SN networks; blue triangles) provide regional coverage of the
study area. Beach balls show focal mechanisms ofM
w
6.4 andM
w
7.1 mainshocks, as determined by the
USGS. Solid black lines indicate Quaternary fault traces, and the black dashed line delineates the focus
region targeted in this study.
92
with nominal inter-station spacing of 10 km, make up the rst such array, which we herein refer to as 3J.R.
Thirty-one nodal seismometers, centered on an isolated swarm of seismicity on the Garlock Fault to the
southwest of the mainshocks, operating between 8 August and 10 September, and deployed with nominal
inter-station spacing of 5 km, make up the second rectangular array, which we herein refer to as 3J.G. All
data from the 3J.R and 3J.G arrays were sampled at 500 s
1
(see Catchings et al. (2020) for details).
The GS network added ten new stations (Cochran et al. 2020), which began recording as early as 7
July 2019. All ten sites were occupied by strong-motion accelerometers, and six sites were co-located with
broadband seismometers. In this study, we use data from the six broadband seismometers, all of which
were sampled at 100 s
1
.
The ZY network added nine stations (Cochran et al. 2020), each with a co-located broadband seismome-
ter and strong-motion accelerometer, which also began recording as early as 7 July 2019. Here, we use data
from all nine broadband seismometers, which were sampled at 100 s
1
for all but one seismometer, which
was sampled at 200 s
1
.
The permanent regional network, operated by Caltech and assigned the CI network code (California
Institute of Technology and United States Geological Survey Pasadena 1926; Southern California Earth-
quake Center 2013), provides continuity of data coverage before and after the rapid-response deployments,
along with a coarse regional coverage of our study area. In this study, we use data from 55 stations main-
tained by the CI network, plus three stations contributed by the Nevada Seismic Network (network code
NN; University of Nevada Reno 1971), and one station contributed by the Southern Great Basin Network
(network code SN; University of Nevada Reno 1980). All data from the permanent regional network are
recorded by broadband seismometers and are archived with a uniform sample rate of 100 s
1
at the SCEDC.
During the procedure for associating detected phase arrivals with earthquake sources, we use a 1D
velocity model based on the model of Zhang and Lin (2014). We derive updated velocity models from
traveltime tomography for various initial models, including the SCEC Community Velocity Model (CVM)
93
version H15.1 (CVMH; Shaw et al. 2015), SCEC CVM version S4.26 (CVMS; Lee et al. 2014), Modied
Hadley-Kanamori 1D model (HK1D; Hauksson 2000), and a regional model in development (FANG3D),
derived using phase arrivals downloaded from the SCEDC and the same traveltime tomography method
used here. We use the HK1D model as a starting model instead of the 1D model used for associating
phase arrivals to test model robustness by starting with a reasonable model with minimal a priori infor-
mation about local structure. Data for CVMH, CVMS, and HK1D were extracted using the SCEC Unied
Community Velocity Model software (Small et al. 2017).
5.3 Methods
In this section, we describe key elements of the workow used to derive an earthquake catalog and tomo-
graphic models from raw waveform data, proceeding through each processing step in order.
5.3.1 Detectingearthquakesandmeasuringphasearrivaltimes
Detecting earthquakes in continuous waveform data is the rst critical procedure in our automated work-
ow, and this problem is intricately connected to that of measuring the arrival times of seismic phases (P-
and S-waves); we treat the two problems in tandem. The description that follows is conceptual, and the
reader is referred to Appendix A for a detailed technical explanation.
We detect earthquakes on a per-station basis by applying a dynamic threshold to a characteristic func-
tion derived from 3C waveform data that targets potential P-wave arrivals. We target P-waves because
they arrive rst and are often clearer than the later arriving S-waves, which can be obfuscated by coda
from the preceding P-wave. The characteristic function combines measures of signal energy (via the ratio
of short- to long-term amplitude averages, or the STA/LTA), the signal kurtosis, and the ratio of hori-
zontally to vertically polarized energy. Our algorithm registers a candidate P-wave arrival whenever the
94
characteristic function exceeds its own average over the preceding 5 s by a factor of 6 or more. After reg-
istering a candidate P-wave arrival, the algorithm measures its arrival time using the Akaike Information
Criterion (AIC; Akaike 1974; Maeda 1985) to determine the sample that optimally divides a small window
of data (1 s centered on the intersection of the characteristic function and threshold) into noise and signal
segments. We process the entire data set in this way to obtain a comprehensive set of candidate P-wave
arrivals.
After detecting potential P-waves and measuring their arrival times throughout the entire data set,
we probe seismograms during the time interval between successive pairs of P-wave arrivals for S-wave
arrivals, assuming that at most one S-wave arrival exists between each successive pair of P-wave arrivals.
Instead of detecting potential S-waves and then measuring their arrival times, we reverse the procedure.
Assuming that an S-wave exists, we estimate its arrival time using the AIC method and then perform
retrospective tests to validate its authenticity, discarding the measurement if it fails these tests. These
retrospective tests comprise a signal-to-noise ratio threshold and a horizontal-to-vertical amplitude ratio
threshold.
5.3.2 Associatingarrivalswithearthquakes
Having detected and measured arrival times of potential P- and S-waves, we associate them with earth-
quake sources using the Rapid Earthquake Association and Location algorithm (REAL; Zhang et al. 2019)
and the 1D average prole of the 3D velocity model for the Coso geothermal area from Zhang and Lin
(2014). We test various algorithm parameters and nd that our results are primarily sensitive to the thresh-
old for the number of associated phases; in the nal analysis, we retain only events associated with at least
eight P-wave arrivals, four S-wave arrivals, and sixteen total arrivals (P and S).
95
5.3.3 Locatingearthquakes
Earthquake hypocenters comprise four coordinates—three spatial and one temporal—which can be inferred
from observed phase arrival times. We locate earthquakes by nding the hypocenter coordinates, h
0
, that
minimize the`
2
-norm of the residual vector between observed and synthetic phase arrival times. That is
h
0
arg min
h
kd p (h)k
2
(5.1)
where d and p are vectors of observed and synthetic phase arrival times for a given event, respectively, and
h is an arbitrary hypocenter four-vector. Determining h
0
for a given set of observations using Equation
(5.1) requires two things: (a) a method to compute synthetic phase arrival times for an arbitrary hypocenter
(the forward problem); and (b) a method to search for h
0
in the space of possible values of h (the inverse
problem).
A convenient software for solving the forward problem is provided by PyKonal (White et al. 2020),
which we use as a core computational engine for both locating earthquakes and updating the velocity
model. Using PyKonal to compute p(h), we obtain an approximate solution to Equation (5.1) using the
Dierential Evolution Optimization algorithm (Storn and Price 1997) implemented by the SciPy Python
package (Virtanen et al. 2020).
5.3.4 Estimatinglocaleventmagnitudes,M
L
,andmagnitudeofcatalogcompleteness,
M
C
To estimate the local magnitude,M
L
, for each event, we measure maximum peak-to-peak S-wave ampli-
tude on horizontal-component data with simulated Wood-Anderson response (A in Equation (5.2) below)
and use the equations of Hutton and Boore (1987, their Equations (1) and (3)) with no station corrections
96
(i.e.,S = 0 in their Equation (1)). That is, we compute the local magnitude from an individual amplitude
measurement made on an instrumentd kilometers from the event hypocenter as
M
L
= log
10
(A) + 1:11 log
10

d
100

+ 0:00189 (d 100) + 3: (5.2)
We computeM
L
for both horizontal components of every station that registered an S-wave arrival and
report the median value as the event magnitude. We also report the median absolute deviation as a measure
of uncertainty.
Following White et al. (2019), we estimate the magnitude of catalog completeness, M
C
, as the 99
th
percentile of the Gaussian component of an exponentially modied Gaussian probability density function
(pdf). That is, we model the observed frequency-magnitude distribution of the catalog using Equation
(5.3):
f
M
(m;;;) exp

2

2 +
2

exp (m)

m

+
2

!
; (5.3)
in which is a Gaussian cumulative distribution function (cdf), and, , and are model parameters
estimated using maximum-likelihood estimation. The 99
th
percentile of the Gaussian pdf corresponding
to is reported as the magnitude of catalog completeness.
5.3.5 Traveltimetomography
Starting with an initial model (i.e., one of CVMH, CVMS, HK1D, or FANG3D), we iteratively derive model
updates via traveltime tomography formulated using the method of Poisson-Voronoi subspace projections
(Fang et al. 2020). This method is used here with three key changes: (a) cell centers are sampled ran-
domly from a non-uniform distribution; (b) cells are contracted along the vertical axis; and (c) arrival time
observations are randomly sampled in proportion to non-uniform weights.
97
Drawing Voronoi cell centers from a non-uniform distribution permits model resolution to vary as
the data allow and structures demand. To obtain variable model resolution, we distribute Voronoi cell
centers vertically in proportion to the vertical velocity gradient and horizontally in proportion to the
station density. To concentrate Voronoi cells where vertical velocity gradients are strongest (e.g., in the
top few kilometers of the crust), their depth coordinates are drawn from a pdf dened by the smoothed
gradient magnitude of the average 1D velocity prole. To concentrate Voronoi cells where station coverage
is densest, their horizontal coordinates are drawn from a pdf dened by the kernel density estimate for the
station distribution with bandwidth determined using Scott’s rule (Scott 2015) as implemented in the SciPy
Python package byscipy.stats.gaussian_kde (Virtanen et al. 2020). This automated choice of bandwidth
eectively increased Voronoi-cell density beneath the 3J.R and 3J.G arrays without over-concentrating the
cells beneath individual stations.
Neglecting discontinuities between lithologic units, we assume that vertical velocity gradients are
larger than their horizontal counterparts, particularly in the top few kilometers of the crust, where the
low conning pressure and surface processes enhance the generation of rock damage. To improve resolu-
tion of vertical velocity gradients, we deform Voronoi cells by compressing them along the vertical axis.
Doing so increases the sensitivity of the inverse problem to vertical velocity gradients at the expense of
decreased sensitivity to horizontal gradients. We dene theaspectratio of a Voronoi cell as the ratio of its
vertical to horizontal dimensions (a high aspect ratio implies a attened pancake-like Voronoi cell).
In the ideal case, raypaths sample the entire model space uniformly. In reality, however, ray cover-
age is strongly heterogeneous. Because our formulation of the tomographic inverse problem comprises
multiple sub-problems, each using a random subset of the observed data, we can sample the observations
in a way that homogenizes ray coverage for each sub-problem. To homogenize ray coverage, we sample
observations in proportion to a weight that we assign each: i.e.,w
i
= e
P (
i
;
i
;z
i
;
i
;
i
)
, in whichw
i
is
the weight assigned to thei
th
ray,
i
,
i
, andz
i
are the latitude, longitude, and depth coordinates of the
98
associated event, respectively,
i
and
i
are the event-to-station azimuth and distance, respectively, and
P () is the joint pdf representing the probability of sampling an observation with the given parameters
from the observed data set. We compute these weights using a 5D kernel density estimate to estimateP (),
with all ve coordinates normalized to lie in [0; 1] and a bandwidth of 0.1. This choice of bandwidth for
normalized data was empirically determined to eectively homogenize ray sampling.
With the algorithm of Fang et al. (2020) modied as described above, we iteratively rene the initial
models by gradually increasing the number of Voronoi cells and cycling through a set of three aspect ratios
(Table 5.1). The rst three iterations use 64 Voronoi cells and an aspect ratio of 1 for the rst iteration,
4 for the second, and 8 for the third. The next three iterations double the number of Voronoi cells (128)
and proceed through the same three aspect ratios (1, 4, and 8). We progressively double the number
of Voronoi cells in this way until a maximum number of 1024 is reached. Each iteration comprises 128
random realizations for each phase (P and S), and each random realization uses 32,648 randomly sampled
observations, which helps mitigate biases due to dierences between the number of observations for each
phase (e.g., having more P- than S-wave observations). The median result is used to update the models
at each iteration, and we relocate all events in the catalog using the updated model after each model
update. We derive the
Vp
=Vs model by directly dividing theV
p
model by theV
s
at each grid node. Model
uncertainties (
P
and
S
for theV
p
andV
s
models, respectively) are estimated by the standard deviation
of the random realizations at each grid node and propagated to the
Vp
=Vs model via Equation (5.4):

P=S
=
s

@
@V
p
V
p
V
s

2

2
P
+

@
@V
s
V
p
V
s

2

2
S
(5.4)
99
Table 5.1: The number of Voronoi cells and aspect ratio used for each iteration of the velocity model update.
Iterations # of Voronoi cells Aspect ratio
1,2,3 64 1,4,8
4,5,6 128 1,4,8
7,8,9 256 1,4,8
10,11,12 512 1,4,8
13,14,15 1024 1,4,8
5.4 Results
5.4.1 Earthquakecatalog
The derived earthquake catalog has a total of 95 533 events from 1 June 2019 through 30 September 2019,
with 94 839 of the events inside our focus region (Figure 6.5), associated with 2:510
6
P- and 1:510
6
S-wave arrival times. The catalog is complete aboveM
L
1:81 (Figure 5.3) and exhibits sharp increases in
the daily event detection rate at the times of theM
w
6:4 andM
w
7:1 mainshocks, reecting the increased
seismicity rate of the aftershock sequences (Figure 5.4). A second sharp increase in the daily detection rate
coincides with the rst phase of the dense deployment (3J.R) on 14 July 2019, followed by a gradual decrease
corresponding to the decaying aftershock intensity. The second phase of the dense deployment (starting
8 August 2019) does not produce a corresponding increase in the daily event detection rate because many
of the additional sensors were deployed in the 3J.G array on the Garlock Fault, which sustained relatively
modest seismicity rates.
The detected seismicity can generally be categorized into ve distinct groups associated with (a) the
eastern Little Lake Fault (ELLF; Plesch et al. 2020), (b) the southern Little Lake Fault (SLLF), (c) the Fremont
Valley, (d) the Coso Range, and (e) the Argus Mountains and Panamint Valley—although it is dicult to
dierentiate seismicity on the ELLF from that on the SLLF where the two faults intersect. We briey
describe prominent characteristics of the spatial distribution of seismicity in each of these groups in the
remainder of this section. Seismicity of the ELLF and beneath the Fremont Valley are discussed in greater
detail in Section 5.4.3.
100
Figure 5.2: Seismicity map (events color-coded by hypocentral depth) showing locations of 94,855 earth-
quakes located in this study. Solid black lines represent cross sections shown in Figures 5.11-5.15. Various
geographical features are labeled, and the locations of theM
w
6.4 andM
w
7.1 mainshocks (small yellow
and large red star, respectively) are shown for reference.
101
Figure 5.3: Observed (black crosses) and modeled (blue curve) frequency-magnitude distribution of the
obtained catalog. The vertical, dashed, red line marks the estimated magnitude of catalog completeness,
M
C
, and the shaded, red area indicates the incomplete portion of the catalog.
102
Figure 5.4: Number of events per day in the catalogs of the SCSN (Hauksson et al. (2012); blue curve), Ross
et al. (2019; orange curve), Shelly (2020; green curve), and this study (red curve). Event counts are scaled
according to the vertical axis on the left-hand side of the plot. Timeline of the number of stations used
(black dashed line), scaled according to the vertical axis on the right-hand side of the plot.
5.4.1.1 EasternLittleLakeFault
The majority of observed seismicity is associated with the ELLF. Seismicity in this group collapses to a
relatively narrow volume (5 km wide between 4 and 16 km deep) over a lateral interval of approximately
20 km, where the main rupture cuts through the Quaternary alluvium of the Indian Wells Valley (Figures
6.5 and 5.7). Events along this segment of the ELLF occur predominantly between 2 and 10 km depth
and dene a single simple lineation in map view. To the southeast of this localized segment, where the
main fault intersects granodiorites of the Argus Mountains, seismicity diuses through a complex net-
work of southwest-/northeast-trending structures that intersect the ELLF at nearly right angles. Further
to the southeast, seismicity gradually shallows towards the Garlock Fault and trifurcates into three discrete
branches before abruptly deepening near the Garlock Fault. Whereas three distinct branches of seismicity
characterize the southeast end of this main group of aftershocks, seismicity at the northwest end, near the
foot of the Coso Range, diuses through a portion of the crust subtended by the ELLF and an10-km-long,
103
orthogonally oriented feature that delineates the northwest terminus of the aftershocks associated with
the ELLF.
5.4.1.2 SouthernLittleLakeFault
Much of the seismicity associated with the SLLF is dicult to dierentiate from that associated with the
ELLF because the two intersect, and we observe the most diuse spatial distribution of seismicity around
the intersection. One distinct cluster of seismicity, however, occurs near the southwest end of the SLLF
between 4 and 8 km depth. Events associated with the SLLF span20 km along the fault and are anked on
the northwest and southeast by multiple parallel lineations that crosscut the ELLF (the complex network
of intersecting features referred to in the paragraph above) and give insight into stepover kinematics.
5.4.1.3 TheGarlockFaultandFremontValley
In addition to the deep seismicity at the southeast end of theM
w
7.1 rupture, the Garlock Fault hosts an
isolated swarm of deep seismicity (8 to 12 km deep) beneath the Fremont Valley to the southwest of the
main rupture. This area was targeted by the second phase of the dense deployment (i.e., the 3J.G sub-array).
The Fremont Valley occupies the region between two strands of the Garlock Fault that are separated by a
stepover, and this swarm of deep seismicity occurred entirely within this stepover region, primarily along
two northeast-/southwest-trending lineations that apparently dip to the southeast. Although the number
of events here is relatively small (2,500), these events are valuable because they illuminate structure on
the enigmatic Garlock Fault that would otherwise remain obscure.
5.4.1.4 CosoRange
A conspicuous seismic gap characterizes the Coso Range south of the Coso Volcanic Field, and to the
northwest of this seismic gap, there lies a roughly 20 by 10 km patch of shallow (< 6 km deep) seismicity.
104
Activity in this region was contemporaneous with that in the main rupture area and forms three substruc-
tures: (a) a very shallow fan of seismicity at the southern end; (b) four pockets of seismicity at the north
end of the shallow fan; and (c) a linearly distributed cluster north of the three pockets of seismicity.
5.4.1.5 ArgusMountains
A nal group of earthquakes forms isolated pockets distributed over a broad region approximately parallel
to the Garlock Fault in the Argus Mountains and further east, starting roughly 15 km to the east of the
diuse seismicity at the northwest end of the ELLF. The distinct clusters of shallow events that compose
this area of seismicity exhibit a weak east/west trend and seem to occur preferentially in regions where
the surface geology is igneous (Figure 5.7), such as in the Argus Mountains, in contradistinction to regions
where the surface geology is sedimentary, such as in the Searles and Panamint valleys.
5.4.2 Velocitymodels
We derive 3D models ofV
p
,V
s
, and
Vp
=Vs structure using the entire earthquake catalog described above
and four dierent starting models (CVMH, CVMS„ HK1D, and FANG3D). In this section, we focus on our
preferred model, which was derived from the CVMH starting model; additional results are available in
the Supplementary Material. The derived models correlate well with surface geology and observed spatial
seismicity patterns. In the following, we describe model residuals (Figure 5.5), resolution (Figure 5.6),
robustness, and salient model features in relation to surface geology and spatial seismicity patterns. In
Section 5.5, we discuss our interpretations of the models.
5.4.2.1 Modelresiduals,resolution,androbustness
Each of the starting models produces a nal model with comparable distributions of residuals between
observations and model predictions (Figure 5.5), with the 1D starting model yielding the largest mist.
Although our preferred model has larger residuals for S-wave data than the other two 3D starting models,
105
Figure 5.5: Final distribution of (a) P-wave and (b) S-wave arrival-time residuals for each of the initial
models. The top, middle, and bottom horizontal line for each distribution marks the 95
th
percentile, mean,
and 5
th
percentile, respectively.
we prefer it for its smoothness, consistency with surface geology, and its relatively modest uncertainty
(particularly in comparison with the CVMS starting model). The mean residual was reduced from 0.081
s (with 0.334 s standard deviation) in the initial model to 0.007 s (with 0.295 s standard deviation) in the
preferred model for P-waves and from -0.127 s (with 0.479 s standard deviation) to -0.014 s (with 0.387 s
standard deviation) for S-waves.
To assess model resolution, we perform a series of checkerboard tests, in which5% anomalies are su-
perimposed on a 1D background model, and we present an illustrative example here (Figure 5.6); additional
results are available to the interested reader in the Supplementary Material.
106
The recovered checkerboard models indicate that, as expected, model resolution is best within the
footprint of the dense deployment and degrades progressively outside of it. The recovered models resolve
internal detail of ne-scale anomalies (10 km horizontal and4 km vertical extent) as deep as 16 km in
key areas of the focus region; however, the top 3 km of the crust (i.e., 2 km below sea level) are poorly
resolved throughout the model, as are regions deeper than 16 km. The relative northeast/southwest dis-
placement between seismicity on the main faults (ELLF, SLLF, and Garlock) and the dense deployment
arrays (3J.R and 3J.G) imparts an overall northeast/southwest trend to the region of the model that is well
resolved.
Major model features interpreted below are generally consistent for dierent starting models where
checkerboard tests indicate good resolution, with the exception of the low-velocity zone (LVZ) below 6 km
depth associated with the Coso Volcanic Field (discussed below), which is absent in the models derived from
the HK1D and, to a lesser extent, the FANG3D model. This model feature is, however, believed to be realistic
because an associated LVZ to the SE below 8 km depth exists in all derived models and the oversimplied
velocity structure in the model derived from the HK1D starting model is dicult to interpret in light of
the heterogeneity we expect the Coso Volcanic Field to introduce. Thus, although the geometrical details
and amplitudes of anomalies depend on the starting model, their rst-order characteristics are robust and
justify modest interpretation.
5.4.2.2 Surfacegeologyandshallowvelocities
Correlating our velocity models with surface geology provides a rst-order method of validating our mod-
els of the shallow crust, and we indeed nd strong correlations with ourV
p
,V
s
, and
Vp
=Vs models (Figure
5.7). Large-scale features can be generally described as (a) high-velocity regions with
Vp
=Vs > 1:70 coin-
ciding with the predominantly igneous composition of the southeastern Sierra Nevada, Coso, Argus, El
107
Figure 5.6: (a) Map view of a horizontal slice at 4 km depth through the true and recovered models of a
checkerboard test for P-wave resolution. Thick, black lines show the surface traces of vertical transects
(b)CACA
0
and (c)C1C1
0
.
108
Paso, and Rand mountains, and (b) low-velocity regions with
Vp
=Vs < 1:70 coinciding with alluvium in the
Indian Wells, Panamint, Fremont, and Searles valleys.
The highestV
p
andV
s
values in the shallow crust are found in the Sierra Nevada Mountains; however,
station coverage is poor in this region, and our inversion cannot constrain much of this structure beyond
the initial model, resulting in high-velocity “bullseye” features. Comparable V
p
values characterize the
Coso, Argus, El Paso, and Rand mountains, but values in the Coso and Argus mountains are notably lower
than those in the Sierra Nevada Mountains—possibly because the lithology associated with the Sierra
Nevada Mountains extends to greater depths, making the structure easier to resolve. Alluvium overlies
granodiorite at the intersection between the ELLF and the Argus Mountains, which manifests as a region
of locally reducedV
p
relative to the higherV
p
associated with exposed granodiorite north and south of the
intersection. These lowV
p
values may also be related to rock damage associated with the SLLF and ELLF
(Qiu et al. 2021).
The lowest V
p
and V
s
values in the shallow crust are found in the Indian Wells Valley, and the as-
sociated basin structure in our models coincides with the shape of the basin, as inferred from the surface
geology. Our models show a steep velocity gradient at the edge of the basin at the foot of the Sierra Nevada
Mountains but a more gradual transition zone at the Coso and Argus mountains.
The Garlock Fault cuts through two low-velocity regions. The rst, at the SE end of the ELLF, correlates
with sedimentary geology of the Searles Valley and may have played a key role in terminating theM
w
7:1
rupture. The second coincides with the fault stepover in the Fremont Valley and the associated swarm of
deep seismicity there. Both of the features also coincide with reduced
Vp
=Vs, and they are separated by a
region of high
Vp
=Vs that correlates well with rhyolitic surface geology. We discuss these features in greater
detail in Section 5.5.
109
Figure 5.7: Horizontal slices of (a) P-wave velocity and (b) S-wave velocity at 2 km below sea level, (c) map
of surface geology compiled by the USGS with legend below, and (d) elevation map.
110
Figure 5.8: Horizontal slices through our V
p
model at various depths (referenced relative to sea level).
Panels in the top row are scaled according to the top color bar, and panels in the second row are scaled
according to the bottom color bar. Seismicity within500 m of each depth slice are shown as white dots.
111
Figure 5.9: Horizontal slices through ourV
s
model, as discussed in Figure 5.8.
112
Figure 5.10: Horizontal slices through our
Vp
=Vs model, as discussed in Figure 5.8.
113
5.4.2.3 Horizontalslices
Our velocity models generally decorrelate with surface geology below 4 km depth (Figures 5.8-5.10). Promi-
nentV
p
andV
s
low-velocity zones (LVZs) are observed near the northwest end of the ELLF below 6 km
depth; however, they do not coincide exactly. TheV
p
LVZ is only faintly visible at 6 km depth, is most
prominent at 10 km depth, and migrates toward the southeast with increasing depth, localizing directly
beneath the northwest terminus of the ELLF. TheV
s
LVZ, on the other hand, is most prominent near 6
km depth and shows little lateral migration with depth, suggesting laterally varying uid content. A re-
gion of low
Vp
=Vs occurs at 10 km depth, where theV
p
LVZ localizes beneath the northwest terminus of
the ELLF and grows to envelope the entire ELLF by 14 km depth. A similar
Vp
=Vs feature is observed as
shallow as 6 km deep beneath the Fremont Valley. This low
Vp
=Vs feature beneath the Fremont Valley is
the result of a colocatedV
p
LVZ without a correspondingV
s
LVZ. In this instance, the low
Vp
=Vs feature is
most prominent at 10 km depth. This leads us to an important observation of this study: deep earthquake
source regions associated with the Ridgecrest earthquake correlate with regions of lowV
p
and
Vp
=Vs. To
the east of the Fremont Valley, aV
s
LVZ below 6 km depth with high
Vp
=Vs between 0 and 14 km depth
coincides with rhyolitic surface geology and a segment of the Garlock Fault, along which Barnhart et al.
2019 observed aseismic surface creep.
5.4.2.4 StructuralvariabilityparalleltotheELLF
Vertical cross sections taken parallel to the ELLF (00
0
through 44
0
; Figures 6.5 and 5.11-5.13) elucidate
variable seismicity patterns and velocity structure along the strike of the ELLF and with depth.
The majority of observed seismicity (77%) occurred within 5 km of the primary trace of the ELLF
(section 2 2
0
), and theM
w
7.1 mainshock occurred in a region with a fault-parallel gradient inV
s
below
6 km depth (Figure 5.12). Southeast of the mainshock, where most aftershocks occurred, ourV
s
model
has stable values near3.7 km s
1
, consistent with V
s
expected of granitic to dioritic composition in
114
that depth range (Christensen 1996). OurV
p
model shows a maximum value of6.3 km s
1
in the same
region, consistent with velocities expected for the granitic end of the granite-diorite range of composition,
although the spatial extent of this feature in ourV
p
model diers from that in theV
s
model, which produces
a correspondingly low
Vp
=Vs region and may be evidence of faulting, fracturing, or temperature dierences.
To the northwest of the mainshock,V
s
decreases to 3.55 km s
1
andV
p
is limited to6.0 km s
1
in
a wedge-shaped structure that pinches out to the northwest;V
p
beneath the wedge is as low as5.5 km
s
1
. V
p
within the wedge is consistent with that of a granitic gneiss, and decreasedV
p
below the wedge
can be explained by increased temperature or shearing. AverageV
p
andV
s
values between 4 and 12 km
depth beneath the Coso Volcanic Field (near -40 km horizontal oset in Figures 5.11 and 5.12) of5.8 and
3.5 km s
1
, respectively, are consistent with various lithologies including basaltic, which is a favourable
interpretation because it is consistent with nearby surface geology.
5.4.2.5 Fault-normalsections
Vertical cross sections taken perpendicular to the ELLF (AA
0
throughKK
0
; Figures 6.5 and 5.14-5.16)
elucidate variable seismicity patterns and velocity structures across the ELLF and with depth.
A lateral gradient from the higher velocity crust beneath the Sierra Nevada Mountains to the lower
velocity crust beneath the Coso Range characterizes the northwesternmost fault-normal proles (sections
AA
0
andBB
0
). The velocity gradient in the top 4 km of the crust reverses (low on the southwest
and high on the northeast) near the Coso Volcanic Field (section CC
0
), with velocities below 4 km
depth gradually transitioning to a more homogeneous distribution in sectionsDD
0
throughFF
0
.
Strong lateral velocity gradients in the top 4 km decrease in sectionsGG
0
throughKK
0
, except for
low-velocity basins in sectionsII
0
andJJ
0
, and the cross-fault velocity contrast reverses polarity
southeast of sectionFF
0
(along which the projected hypocenter of theM
w
7.1 event occurred), with
the northeast side of the fault being higher velocity than the southwest side.
115
Figure 5.11: Northwest-/southeast-trending, fault-parallel vertical slices (traces 0 0
0
through 4 4
0
on
Figure 6.5) of ourV
p
model at 10 km intervals. Seismicity within5 km of each vertical plane is projected
onto the plane and shown as white dots. The large and small stars on section 2 2
0
mark the projected
locations of theM
w
7.1 andM
w
6.4 mainshocks, respectively. Solid white lines indicate contours of relative
uncertainty,
P
. The single- and cross-hatch patterns indicate regions where 1:75%
P
< 2:5% and
2:5%
P
< 5%, respectively, and
P
< 1:75% in unhatched regions.
116
Figure 5.12: Vertical slices through ourV
s
model, as discussed in Figure 5.11. As in Figure 5.11, the single-
and cross-hatch patterns indicate regions where 1:75%
S
< 2:5% and 2:5%
S
< 5%, respectively,
and
S
< 1:75% in unhatched regions.
117
Figure 5.13: Vertical slices through our
Vp
=Vs model, as discussed in Figure 5.11. the single- and cross-
hatch patterns indicate regions where 2:5%
P=S
< 5% and 5%
P=S
< 10%, respectively, and

P=S
< 2:5% in unhatched regions.
118
Seismicity patterns are simplest where lateral velocity structure is relatively homogeneous (e.g., sec-
tions EE
0
and FF
0
), whereas seismicity is more complex and broadly distributed where lateral
velocity gradients exist (e.g., sectionsDD
0
andGG
0
throughII
0
).
5.4.3 Keystructuralfeatures
5.4.3.1 CentralELLF
TheM
w
6.4 earthquake nucleated in a region with a strongV
s
gradient across the SLLF at 12 km depth
and propagated primarily to the southwest towards a LVZ (Figure 5.17). Similarly, theM
w
7.1 earthquake
nucleated in a region with strongV
s
gradients along and across the ELLF at 4 km depth. The entire central
segment of the ELLF is characterized by cross-faultV
s
contrasts (slower on the southwest) in the top 4 km,
with the steepest gradients to the northwest of theM
w
7.1 and southeast of theM
w
6.4 at 4 km depth.
5.4.3.2 NorthwestELLF
At the northwest terminus of the ELLF, a cross-faultV
s
contrast is observed below 6 km depth, with the
low-velocity side to the northeast (Figure 5.18). Above 6 km depth, thisV
s
contrast reverses. Similarly,
relatively low V
s
to the north of the fault terminus (below 6 km depth) is replaced by relatively high
velocities above 6 km.
5.4.3.3 SoutheastELLF
The southeast terminus of the ELLF exhibits a strong fault-parallelV
s
gradient that decreases toward the
southeast in the shallow crust (top 4 km), where seismicity trifurcates and shallows (Figure 5.19). Below 6
km depth, highV
s
localizes to a “bullseye” that encloses the trifurcation region.
119
Figure 5.14: Southwest-/northeast-trending fault-normal vertical slices (tracesAA
0
throughKK
0
on
Figure 6.5) through ourV
p
model, as discussed in Figures 5.11-5.12. Relative-uncertainty contours are as
in Figure 5.11.
120
Figure 5.15: Southwest-/northeast-trending fault-normal vertical slices through ourV
s
model, as discussed
in Figure 5.14. Relative-uncertainty contours are as in Figure 5.12.
121
Figure 5.16: Southwest-/northeast-trending fault-normal vertical slices through our
Vp
=Vs model, as dis-
cussed in Figure 5.14. Relative-uncertainty contours are as in Figure 5.13.
122
Figure 5.17: (a) Seismicity in the region surrounding the central segment of the ELLF and (b)-(f) horizontal
slices of theV
s
model at 2, 4, 6, 9, and 12 km depth.
123
Figure 5.18: (a) Seismicity in the region surrounding the northwest terminus of the ELLF and (b)-(f) hori-
zontal slices of theV
s
model at 2, 4, 6, 8, and 10 km depth.
124
Figure 5.19: (a) Seismicity in the region surrounding the southeast terminus of the ELLF and (b)-(f) hori-
zontal slices of theV
s
model at 2, 4, 6, 8, and 10 km depth.
125
Figure 5.20: (a) Seismicity in the region of the Fremont Valley and (b)-(f) horizontal slices of theV
s
model
at 2, 4, 6, 8, and 10 km depth.
5.4.3.4 GarlockFaultatFremontValley
The lowV
s
of sediment in the Fremont Valley is clearly apparent in the shallow crust (top 4 km; Figure
5.20). Moderately reducedV
s
is also observed at 8-10 km depth in a region that is bound on the north and
south by mapped strands of the Garlock Fault and to the east and west by detected seismicity. Bounding
seismicity to the east and west form roughly northeast-/southwest-trending lineations that apparently dip
to the southeast.
126
5.5 DiscussionandConclusions
We processed four months of raw waveform data recorded starting one month before the M
w
6.4 and
M
w
7.1 mainshocks using automated procedures to derive an earthquake catalog with over 95,000 events.
We then perform 3D, ray-based tomography modeling using data from this catalog to relocate earthquakes
and obtain detailed models of seismic P- and S-wave speeds (V
p
andV
s
) and their ratio
Vp
=Vs. Our derived
seismic catalog contributes to the growing list of earthquake catalogs for the 2019 Ridgecrest earthquake
sequence (Lee et al. 2020a; Lee et al. 2020b; Lin 2020; Liu et al. 2020; Lomax 2020; Ross et al. 2019; Shelly
2020), but our catalog is unique in at least three ways: (a) it is the only catalog derived using recordings
from the rapid-response dense-deployment of seismometers (Catchings et al. 2020); (b) it is the temporally
longest catalog, focusing exclusively on the Ridgecrest event sequence; and (c) it is the only catalog built
using fully 3D location methods and velocity model updates. Our catalog contributes the largest number
of new phase arrival times to the standard catalog from the SCSN. Furthermore, the events in our catalog
can be used as templates to detect many more events that produce similar waveforms (e.g., Shelly et al.
2016; Ross et al. 2017).
The velocity structure in the upper 8 km of the crust along the segment of the ELLF southeast of the
M
w
7.1 epicenter is characterized by relatively highV
p
andV
s
beneath the Argus Mountains and Span-
gler Hills and lowV
p
andV
s
near the southeastern end of the rupture, where it abuts the Garlock Fault.
The network of orthogonal branches of seismicity along this segment of theM
w
7.1 fault have been noted
in other studies (e.g., Ross et al. 2019) and are reminiscent of the “shatter networks” observed near the
Trifurcation Area of the San Jacinto Fault Zone (White et al. 2019). These complex structures may result
from failure of highly brittle crust, a hypothesis consistent with high shear moduli that explains the ob-
served elevatedV
p
andV
s
. Strain diuses across three sub-parallel strands of seismicity, where the highest
seismic velocities along the ELLF are observed (immediately northwest of the Garlock Fault). The abrupt
deepening of seismicity near the Garlock Fault may result from a highly compliant (i.e., low shear moduli
127
and thus low-velocity) crust associated with the Garlock Fault that absorbs strain and reduces the stress
at shallow depths below the level needed to produce brittle failure of detectable magnitude. This highly
compliant portion of the crust, which we interpret as due to a combination of weak sedimentary lithology
at the surface and fault damage at depth, likely played a key role in terminating the southeast rupture of
the ELLF during theM
w
7.1 mainshock by readily deforming and dissipating the stress.
Northwest of theM
w
7.1 epicenter, the ELLF penetrates the Indian Wells Valley, and aftershock activity
localizes along a single linear (in map view) zone. Cross-fault velocity contrasts as large as 10.5% forV
p
and
13.8% forV
s
in the upper 4 km of crust along this segment of the ELLF suggest that a bimaterial interface
at the edge of the Indian Wells Valley may have played a key role in localizing the rupture along this
fault segment (e.g., Brietzke and Ben-Zion 2006). Seismicity at the northwest end of the ELLF terminates
along a zone that trends orthogonal to the ELLF, near the transition between Indian Wells Valley and
the Coso Mountains. Apparently, the rupture on the ELLF did not have suciently concentrated stress
to penetrate the more rigid (higher velocity) Coso Mountains and was diverted to failures in the weaker
(slower) sediments of the Indian Wells Valley. The high angle between ruptures in the Indian Wells Valley
and the Ridgecrest area is expected of ruptures that decelerate rapdily upon encountering a barrier (Xu and
Ben-Zion 2013). We hypothesize that opposing structures caused the arrest of theM
w
7.1 rupture of the
ELLF at either end: (a) a strong barrier to the northwest (i.e., the Coso Mountains) that the rupture could
not penetrate, and (b) a compliant buer that dissipated stress to the southeast (i.e., the Garlock Fault).
The proposed hypothesis, which is necessarily retrospective, is primarily intended to be explanatory,
not predictive. Our aim is not to suggest a predictive framework for earthquake cessation, but rather
to propose potential mechanisms that explain the cessation of this single earthquake. Only after many
similar observations would we hope to extract patterns with predictive power. We suggest then, that the
hypothesis be used not to make predictions about where a future earthquake rupture might arrest, but to
128
guide future research looking for earthquake cessation patterns that may eventually yield such predictive
capacity.
It is interesting to note that where seismicity trifurcates at the southeastern end of the ELLF, it pen-
etrates a region with greaterV
s
(and thus greater rigidity) than the region that we hypothesize acted as
a rigid barrier at the northwestern end. This indicates that our hypothesis is, at best, an incomplete ex-
planation. The variables that aect the ability of a rupture front to penetrate a rigid body are manifold:
e.g., the amount of stress concentrated along the rupture front, the rupture orientation with respect to the
boundary of the rigid body and any internal planes of weakness, the geometry of the rigid body, and the
geometry of the rupture front. Speculating on the eect of these variables in the present case is beyond
the scope of this paper, as they are either unknown or uncertain, but their eects might be systematically
tested in a laboratory or computational environment.
The isolated swarm of seismicity in the stepover region along the Garlock Fault form northeast-
/southwest-trending lineations that extend between the two main branches of the Garlock Fault. Given the
right-lateral sense of motion on the ELLF during the main rupture, these lineations likely represent accom-
modation of transtensional stresses induced on the Garlock Fault by theM
w
7.1 mainshock in a pull-apart
basin structure.
Our earthquake tomography model has low resolution in the shallow crust, where important velocities
are needed for ground motion estimation (Seyhan and Stewart 2014; Juarez and Ben-Zion 2020), and our
results are also unable to resolve internal fault detail in the immediate vicinity of the Ridgecrest rupture.
Improved results for the shallow crust can be obtained using active sources (e.g., Louie 2001; Schuster 2009)
or the ambient seismic noise (e.g., Lin et al. 2013; Zigone et al. 2019). Detailed imaging of the rupture zone,
on the other hand, is best accomplished using seismic phases that propagate along and across the fault zone
and has been targeted in a separate study (Qiu et al. 2021). The velocity structure in our study area has
various low-velocity zones at depth that are likely produced by a combination of lithological variations,
129
uids, temperature, and rock damage. Low
Vp
=Vs in deep earthquake source regions are consistent with the
observations of Lin and Shearer (2009), which they interpreted as due to the presence of uids (i.e., water).
Using the derived catalog to estimate earthquake-induced rock damage (Ben-Zion and Zaliapin 2019) can
improve estimates of rock types and uid content in various locations.
DataandAcknowledgments
The derived seismic catalog andV
p
,V
s
, and
Vp
=Vs velocity models can be obtained at https://data.mendeley.
com/datasets/x8v5wkbj6r and https://data.mendeley.com/datasets/gv33tgvt5f, respectively. All gures were
created using Matplotlib (Hunter 2007). Numerous open-source Python packages were critical in complet-
ing this work, including but not limited to NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Pandas
(McKinney 2010), ObsPy (Beyreuther et al. 2010), and IPython/Jupyter (Perez and Granger 2007). Focal
mechanism data for theM
w
6.4 andM
w
7.1 mainshocks shown in Figure 5.1 were obtained from the USGS
website at https://earthquake.usgs.gov/earthquakes/eventpage/ci38443183/moment- tensor (last accessed 14
December 2020) and https://earthquake.usgs.gov/earthquakes/eventpage/ci38457511/moment- tensor (last
accessed 14 December 2020), respectively. The geologic data shown in Figure 5.7 were downloaded from
the USGS website (https://mrdata.usgs.gov/geology/state/state.php?state=CA ; last accessed 1 July 2020).
Fault traces of Quaternary faults in Figures 5.1, 6.5, 5.7-5.9, and 5.17-5.20 were obtained from the USGS
and California Geological Survey, Quaternary fault and fold database for the United States (https://www.
usgs.gov/natural- hazards/earthquake- hazards/faults ; last accessed 1 December 2019). Traces of surface
ruptures from the 2019 Ridgecrest sequence shown in Figures 6.5, 5.7-5.9, and 5.17-5.20 were obtained from
the Supplementary Material accompanying Ponti et al. (2020). Locations of theM
w
6.4 andM
w
7.1 main-
shocks from Lomax (2020) were used in Figures 6.5 and 5.11-5.17. Data from the 3J network are described
in Catchings et al. (2020). Data from the GS and ZY networks are described in Cochran et al. (2020). Data
130
from the 3J and GS networks and can be accessed through the facilities of IRIS Data Services, and speci-
cally the IRIS Data Management Center. IRIS Data Services are funded through the Seismological Facilities
for the Advancement of Geoscience (SAGE) Award of the National Science Foundation under Cooperative
Service Agreement EAR-1851048. Data from the ZY network can be accessed though the Southern Cali-
fornia Earthquake Data Center Southern California Earthquake Center (2013). Research for this paper was
supported by the National Science Foundation (grants EAR-1945781 and EAR-1722561) and the Southern
California Earthquake Center (based on NSF Cooperative Agreement EAR-1600087 and USGS Cooperative
Agreement G17AC00047). We thank Frederik Tilmann, an anonymous referee, and Editor Ana Ferreira for
useful comments that helped us to improve the manuscript.
131
Chapter6
UpdatingthecatalogfortheSanJacintoFaultZone
Preface
The content of this chapter comes from an article currently in preparation (as of July 28, 2021) for submis-
sion to Seismological Research Letters under the title Catalog Update (2008-2020): A Detailed Earthquake
CatalogfortheSanJacintoFault-ZoneRegioninSouthernCalifornia. It presents ongoing work to extend and
update the catalog for the San Jacinto Fault Zone presented in Chapter 2 using additional data and rened
methods. The catalog in Chapter 2 spanned 2008 through 2016; we are now extending the catalog through
2020. We also make two key changes to the processing procedure from Chapter 5: a) we use a deep-learning
model to detect earthquakes and pick phase arrivals, and b) we implement a new Markov-Chain Monte
Carlo (MCMC) sampling approach to locate earthquakes. The benets of the new MCMC-based location
algorithm are that it a) accounts for observational and prediction errors in a general and exible way, and
b) it provides robust uncertainty estimates. Using the updated procedure and additional data, we detect
and locate more than 2.5 times as many events as the standard regional catalog. These newly detected
earthquakes reveal previously unseen internal details of fault-zone structure.
132
6.1 Introduction
6.2 Data
Five dierent networks contribute data from 136 seismic stations that we use in this study (Figure 6.1).
Data from 68 stations from theCI network (California Institute of Technology and United States Geological
Survey Pasadena 1926; Southern California Earthquake Center 2013) provide coarse regional coverage of
the study area. Data from 29 stations from theAZ network (Vernon 1982), 28 stations from theYN network
(Vernon and Ben-Zion 2010), 8 stations from thePB network, and 1 station from theSB network provide
dense coverage targeting our focus region.
Data from the CI network are downloaded from the Southern California Earthquake Data Center’s
(SCEDC; Southern California Earthquake Center 2013) archive hosted by Amazon Web Services. All data
from the CI network used in this study were archived with a sample rate of 100 s
−1
, except two stations (SVD
and DGR, which which were sampled at 80 s
−1
until 10 October 2008 and 4 December 2009, respectively,
after which they were sampled at 100 s
−1
). Fifty-eight CI stations used in this study contribute data from
three-component (3C), high-gain, broadband seismometers. Nine contribute data from 3C, strong-motion
accelerometers and single-component (1C), high-gain, short-period seismometers; we preferentially select
3C data for these stations, resorting to 1C data when 3C data is unavailable, if possible. A single CI station
used in this study (SKY) contributes only 1C data (from a high-gain, short-period seismometer).
Data for all stations from the AZ, YN, PB, and SB networks are downloaded from the Incorporated
Research Institutions for Seismology Data Management Center (IRIS DMC).
Nineteen AZ stations contribute data from 3C, broadband, high-gain seismometers, variably sampled
at 40, 100, or 200 s
−1
. Ten contribute data from 3C, strong-motion accelerometers, variably sampled at 200
or 250 s
−1
.
133
Figure 6.1: Locations of 136 seismic stations used in this study. Gold squares indicate the locations of
stations contributing data from strong-motion accelerometers to this study, inverted, red triangles indicate
stations contributing data from high-gain seismometers, and blue circles indicate stations contributing data
from both types of instruments.
134
Twenty-three YN stations contribute data from 3C, high-gain seismometers (11 short-period and 12
broadband), uniformly sampled at 200 s
−1
. Five contribute data from strong-motion accelerometers, uni-
formly sampled at 200 s
−1
.
All eight PB stations contribute data from 3C, high-gain seismometers installed in boreholes between
125 and 230 m below the surface and variably sampled at 100 or 200 s
−1
.
Finally, the SB network contributes 3C data from a single, strong-motion accelerometer installed in a
borehole 100 m below the surface.
The 3D P- and S-wave velocity models used in this study are taken directly from White et al. (2019),
which were derived by combining results from Fang et al. (2016) and Shaw et al. (2015).
6.3 Methods
In this chapter, we further modify the basic procedure from White et al. (2019) and White et al. (2021).
Instead of the statistical methods previously used to detect earthquakes and estimate (pick) phase arrival
times, we use a deep-learning model that performs these tasks simultaneously (EQTransformer; Mousavi
et al. 2020). And, instead of using the NonLinLoc algorithm (Lomax et al. 2009) or Dierential Evolution
Optimization (DEO; Storn and Price 1997) to locate events as in White et al. (2019) and White et al. (2021),
respectively, we implement an MCMC sampling approach that takes into account observational and pre-
diction errors and yields uncertainty estimates for hypocenter coordinates. As in White et al. (2019), we
cross-correlate waveforms to measure dierential arrival times for relocating event ensembles using the
GrowClust double-dierence algorithm (Trugman and Shearer 2017). In the remainder of this section, we
describe key elements of the revised procedure.
135
6.3.1 Detectingearthquakesandpickingphasearrivals
To detect earthquakes and pick P- and S-wave arrivals, we use the EQTransformer algorithm (Mousavi
et al. 2020), a deep-learning model for simultaneously solving these coupled problems using a single en-
coder to encode high-level representations of seismograms and three separate decoders to decode these
(representations) into probability sequences representing the (a) presence of an earthquake signal, (b) on-
set of P-wave energy, and (c) onset of S-wave energy. Although the data set used to train the EQTrans-
former model (Mousavi et al. 2019) comprises recordings of earthquakes from diverse tectonic settings
and geographical locations at a broad range of hypocentral distances, all were from seismometers (i.e.,
no recordings from accelerometers). Despite this shortcoming, we nd that the model generalizes to data
recorded by accelerometers well. Furthermore, because the training data set comprises many recordings
from earthquakes in Southern California, it performs well with minimal or no hyperparameter tuning. In
this chapter, we use the model with default hyperparameter values distributed by its authors.
6.3.2 Associatingphasearrivalswithearthquakesources
As in White et al. (2021), we use the Rapid Earthquake Association and Location algorithm (REAL; Zhang
et al. 2019) to associate phase arrivals picked by EQTransformer with their earthquake sources. We derive
a 1D velocity model from the 3D model used in White et al. (2019) and, using this 1D model, register events
that associate with at least four P-wave arrivals and six total (P- and S-wave) arrivals. This low association
threshold maximizes the number of detected events at the expense of many false positives, which will need
to be discarded during subsequent processing steps.
6.3.3 Estimatingabsoluteeventlocations
To estimate absolute locations for individual earthquakes, we implement an MCMC sampling procedure
to approximate the posterior distribution for hypocenter coordinates and report themaximumaposteriori
136
(MAP) estimate as the event location and mean absolute deviation (from the MAP estimate) as associated
uncertainties. To do so, we use Bayes’ Theorem as follows.
Bayes’s Theorem can be stated as
P

m

d

=
P

d

m

P (m)
P (d)
(6.1)
in which m and d represent vectors of model parameters (three spatial and one temporal hypocenter
coordinates) and observed data (phase arrival times), respectively,P (mjd) is the probability that model
parameters m describe the model that generated observations d (the posterior distribution),P (djm) is
the probability of observing d if m is the generative model (the likelihood), andP (m) andP (d) are the
probabilities of the model parameters and data observations being sampled from the set of possible values
(the prior and marginal probabilities), respectively.
To let the data dominate the inversion procedure, we impose minimalapriori information on hypocen-
ter coordinates by choosing a non-informative prior; i.e., a uniform distribution over a nite region of
Euclidean space and time,
:
P (m) =P (x;y;z;t) =
8
>
>
>
>
<
>
>
>
>
:
for (x;y;z;t)2

0 otherwise;
(6.2)
in which is an appropriate constant such that P integrates to unity, and x, y, z, and t are Cartesian
hypocenter coordinates and time, respectively, or
137
P (m) =P (;;;t) =
8
>
>
>
>
<
>
>
>
>
:

2
sin for (;;;t)2

0 otherwise
(6.3)
in which,, and are spherical hypocenter coordinates (radius, polar angle, and azimuth, respectively).
As is common practice in Bayesian analyses, we further simplify the problem by neglecting the marginal
probability,P (d), so that, parameterizing the model using spherical coordinates, we have
P

m

d

/P

d

m

2
sin: (6.4)
It is worth noting here that there is some controversy regarding how to dene prior probabilities,
with no single choice satisfying all philosophical perspectives (e.g., frequentists and strict Bayesians). We
acknowledge the need for further comment on our choice and defer such discussion to Section 6.5.
To dene the likelihood function,P (djm), we model the vector of residuals, r, between observations
and predictions as the sum of a random observational error (i.e., pick error) and a random prediction
error (i.e., traveltime error generated by inaccurate models of velocity structure): i.e., r = "
obs
+ "
pred
,
in which "
obs
and "
pred
are independent, random vectors representing the observational and prediction
errors, respectively. Then, letting d =
^
d (m) + r, in which
^
d (m) represents the vector of predicted
arrival times for hypocenter m, we have
P

m

d

/P

^
d (m) + r

m

(6.5)
=P

^
d (m)

m

~P

r

m

; (6.6)
138
in which~ represents the convolution operator. Because
^
d (m) is deterministic, its PDF is a Dirac-delta
function, and we thus have
P (mjd)/P (rjm) (6.7)
=P

"
obs
+"
pred

m

(6.8)
=P

"
obs

m

~P

"
pred

m

(6.9)
With an appropriate model forP (rjm), this formulation is suitable for approximating the posterior dis-
tribution using MCMC sampling.
We develop station- and phase-dependent models of observational errors by tting a set of Cauchy
distributions to the histograms of residuals between automated picks produced by EQTransformer and
manual picks made by SCSN analysts for each station and phase (Figure 6.2). Because we are ignorant of
the true distribution of traveltime errors produced by inaccurate velocity models, we assume they follow
a normal distribution that depends only on the phase (i.e., are station independent). To determine the
mean and standard deviation of these two models, we t them to the histograms of residuals between
observations and arrival times predicted using preliminary 3D locations obtained with the DEO method
from White et al. (2021). This is an overestimate of the prediction error because it assumes that the entire
mist between observations and predictions can be explained by the prediction error alone (i.e., assuming
there is no observational error). Modeling the observational and prediction errors with Cauchy and normal
distributions, respectively, yields a Voigt prole for the total residual, which is dominated by the traveltime
errors.
With the likelihood function dened using the error model described above, we approximate the pos-
terior distribution for the hypocenter coordinates (e.g., Figure 6.3) using the MCMC sampler implemented
139
Figure 6.2: Example (a) observational, (b) prediction, and (c) total residual models for P-wave arrivals
observed at station AZ.BZN. (a) The solid black curve is the histogram of residuals between automated
and manual P-wave picks at station AZ.BZN, and the orange curve is the best-tting Cauchy distribution.
(b) The blue curve is the normal distribution used to model station-independent traveltime errors. (c)
The green curve is the Voigt prole used to model the distribution of total residuals for P-wave arrivals
observed at station AZ.BZN.
140
by the emcee Python package (Foreman-Mackey et al. 2013). We report the MAP estimate as hypocen-
ter coordinates and 2.5 times the mean absolute deviation of MCMC samples from the MAP estimate as
associated uncertainties.
6.3.4 Relocatingeventensembleswithdoubledierences
After estimating absolute hypocenter coordinates for each event individually, we cross-correlate wave-
forms from each event with waveforms from its 200 nearest-neighbor events to measure dierential ar-
rival times. Waveforms are rst bandpass ltered between 3 and 20 Hz, and then three-component (3C)
template waveforms are extracted for each registered arrival (starting 0.25 s before the arrival and ending
0.75 s after it for P-waves, and starting 0.25 s before the arrival and ending 1.25 s after it for S-waves). Each
3C template is then cross-correlated with 3C test waveforms from each of its 200 nearest-neighbor events
that have registered arrivals for the same station and phase; test waveforms are all 1 s long and are cen-
tered on the arrival. Dierential arrival times with maximum cross-correlation coecient less than 0.75
are discarded and remaining results are averaged for each event, station, and phase triplet. Observations
are discarded for event pairs with fewer than six dierential arrival times at this point, and the remaining
observations are used (with the same 1D model used to associate arrivals with sources) to relocate events
using GrowClust. We perform 10 bootstrap samples to estimate relative location uncertainties.
6.3.5 Calculatinglocalmagnitudes,M
L
To estimate the local magnitude for each event, we use the same method as White et al. (2021). See Section
5.3.4 of Chapter 5 for details.
141
Figure 6.3: Example 2D (o-diagonal subplots) and 1D (diagonal subplots) marginalized posterior sample
distributions for hypocenter coordinates sampled using MCMC sampling. Solid, orange lines and red
squares indicate MAP parameter estimates. Dashed, blue lines (diagonal subplots) indicate uncertainty
estimates (2.5 times the mean absolute deviation) for each coordinate; MAP and uncertainty estimates for
each coordinate are indicated by diagonal subplot titles.
142
6.3.6 Qualitycontrol
We apply quality control metrics that are similar to those used in White et al. (2019), however, we mod-
erately relax the depth-error threshold because uncertainties reported in White et al. (2019) were 68%
condence intervals, whereas we now report 95% condence intervals. Thus, every event retained in the
nal catalog was relocated using GrowClust or is associated with latitude, longitude, depth, and time un-
certainty less than 0.03°, 0.03°, 5 km, and 1 s, respectively.
6.4 Results
6.4.1 Catalogoverview
Before applying quality control metrics to discard events with large uncertainties, the catalog comprises
437 995 events (Figure 6.4). The nal catalog, after quality control, comprises 252 007 events, 179 421 of
which were relocated by GrowClust. Inside our focus region, we detect 163 052 events (Figre??)with me-
dian absolute latitude, longitude, depth, and time uncertainties of 0.011°, 0.014°, 2.2 km, and 0.29 s, respec-
tively. The nal catalog contains 65 % more events than the White et al. (2019) catalog and 152 % more than
the Hauksson et al. (2012) catalog inside our focus region between 2008 and 2016 (Figure 6.6).
6.4.2 Selectedstructuralfeatures
In the remainder of this section, we present selected structural features observed in the catalog. The results
presented are intended merely to acquaint the reader with broad-scale structures and general character-
istics of the catalog—they should not be construed as exhaustive. We present the catalog here as a data
product for derivative studies and leave more detailed analysis and interpretation of spatiotemporal seis-
micity patterns for future work.
143
Figure 6.4: Histograms of uncertainties associated with hypocenter parameters for events in the catalog
before applying quality control metrics. The vertical, red, dashed lines indicate quality control thresh-
olds. Green-shaded regions indicate the subset of events retained after quality control; red-shaded regions
indicate the subset of events discarded.
144
Figure 6.5: Overview map of cataloged seismicity in our focus region.
145
Figure 6.6: Histograms of events detected per month in our focus region for this study (blue), White et al.
(2019, orange), and Hauksson et al. (2012, green).
146
Figure 6.7: Map of seismicity from our catalog in the Hot Springs Area.
6.4.2.1 HotSpringsArea
The Hot Springs Area (Figures 6.5 and 6.7) occupies the northwestern portion of and comprises the deepest
seismicity in our focus region. The majority of seismicity here is distributed along the main Clark fault
branch below 12 km depth, and an abrupt depth transition (deeper to the NW) along this primary fault
branch coincides with diuse limbs splaying with E/SE trend toward the Hot Springs Fault. Deep seismicity
associated with the Hot Springs Fault delineates two discrete structures, parallel to the surface trace of the
fault but oset from one another. O-fault seismicity in this region is relatively shallow (<10 km).
6.4.2.2 TrifrucationRegion
The Trifurcation Region (Figures 6.5, 6.8, and 6.9) occupies the SE portion of our focus region, is the
most seismically active region in the SJFZ (and Southern California), and juxtaposes the Anza Seismic
Gap—a region immediately NW of the Trifurcation Region that is nearly devoid of seismic activity. Three
147
fault strands (viz., the Coyote Creek, Clark, and Buck Ridge strands) converge toward the northwest here,
activating a complex network of entwined substructures.
Seismicity in the southeastern Trifurcation Region is generally deep (>12 km) and is predominantly
distributed between or near the Coyote Creek and Clark fault branches. Although most seismicity here is
diuse and distributed roughly parallel to the surface traces, a few faint, deep lineations cross-cutting the
primary trend are observed.
At the point of convergence in the Trifurcation Region, seismicity manifests two rst-order signatures:
a) deep, diuse seismicity subtended by the Coyote Creek and Buck Ridge faults, and b) shallow networks
of discrete, conjugate fractures west of the Coyote Creek and north of the Buck Ridge faults.
The majority of shallow seismicity to the north of the Buck Ridge fault occurred following the 2013
M
L
4.7 Buck Ridge earthquake and has since been progressively localizing to a single discrete structure,
which is likely the principal subsurface trace of the Buck Ridge fault. Shallow seismicity to the west of the
Coyote Creek fault, on the other hand, appears to have occurred in punctuated periods of intense activity
on individual structures. At least one such structure shows evidence of having reactivated.
Two periods of heightened activity characterize the deep, diuse seismicity subtended by the Coyote
Creek and Buck Ridge faults; one coinciding with the 2016M
L
5.2 Borrego Springs earthquake, and the
other coinciding with the 2020 M
L
4.9 Anza earthquake. The spatial distribution of events associated
with the 2020 Anza earthquake are more compact than those associated with the 2016 Borrego Springs
earthquake.
6.4.2.3 CahuillaSwarm
The Cahuilla Swarm, a four-year period of increased seismic activity near the town of Cahuilla, California,
that started in 2016, reveals some of the most detailed structure observed in our catalog (Figure 6.10). An
east-northeast dipping structure is clearly evident, and internal fault braiding is observed on sub-kilometer
148
Figure 6.8: Map of seismicity from our catalog in the Trifurcation Region.
scale. The resolution of this structure in our catalog is comparable to that of the focused study by Ross
et al. (2020), but we achieve similar resolution with much greater spatial and temporal coverage.
6.5 Discussion
6.5.1 Priorprobabilities
In this chapter we use Bayes’ Theorem to estimate hypocenter coordinates, and our choice for the prior
probability merits discussion. First, note that we could have chosen a prior that reects ourapriori knowl-
edge about the distribution of seismicity (i.e., an informative prior). For example, we could have chosen
a prior to reect the fact that the depth distribution of seismicity in the SJFZ roughly follows a normal
distribution centered at12 km. It is our view, however, that such a choice would impart an unjustied
149
Figure 6.9: Maps of seismicity from our catalog in a small sub-region of the Trifurcation Region, color-
coded by (a) depth and (b) time.
150
Figure 6.10: Map of seismicity from our catalog near the town of Cahuilla, California, where the Cahuilla
Swarm occurred.
151
bias on all hypocenters irrespective of the constraining observations. Instead, we wish to maximally allow
the data to speak for itself.
We do so by encoding in our prior the fact that, within a nite region of Euclidean space and time,
no contained subregion is a priori more likely to host a hypocenter than another contained subregion
of equal volume. This is achieved using a uniform distribution over Cartesian hypocenter coordinates
and time, which transforms to the distribution over spherical hypocenter coordinates and time given by
Equation (6.3). Note that a simple uniform distribution over spherical coordinates (e.g., Smith et al. 2021)
concentrates probability density near the origin and poles of the coordinate system and will thus bias
hypocenters toward those regions.
Strictly speaking then, the hypocenters herein obtained do not correspond to Maximum Likelihood
Estimates (MLEs), which those who espouse the frequentist viewpoint might prefer, because the chosen
prior is non-uniform over model parameters. In the present case, however, the dierence between the
MAP and MLE for a given hypocenter’s coordinates is negligible as the radial and polar-angle intervals
over which the prior is non-zero are suciently narrow that the prior is locally uniform and the likelihood
dominates the posterior.
6.5.2 Derivativestudiesandfurtherwork
The derived catalog is suitable for multiple derivative studies. We have presented a preliminary analysis of
spatiotemporal seismicity patterns here, but the manifold structures illuminated by the catalog merit more
thorough analysis—such analysis is beyond the scope of the present work and left for a followup study.
Traveltime tomography using our new catalog is likely to be a high-yield endeavor because we detect
many new events and associated phase arrivals, and, because we process data from many stations omitted
from the standard regional catalog (Hauksson et al. 2012), many new ray paths are sampled. Because our
152
analysis is independent ofapriori template events, it is also likely that we have detected diverse new event
classes that can be used as templates to further improve catalog completeness via template matching.
In this study, we modeled prediction errors using a pair of phase-dependent normal distributions tted
to the residuals between observations and predictions made using preliminary hypocenter coordinates.
Instead of tting the distributions in this way, the location and scale parameters (means and variances,
respectively) could be treated as additional model parameters to be inverted for. Estimating parameters
for prediction error models using MCMC sampling would not only be more accurate than those used, thus
making hypocenter estimates more robust, but would also be useful for quantifying uncertainty associated
with inaccurate velocity models. Such analysis is beyond the scope of the present work, but may form a
useful followup study.
6.6 Conclusions
In this chapter, we further rene the automated processing procedure rst presented in Chapter 2 and
improve the SJFZ catalog by applying it to thirteen years of data (2008-2020) to derive a new catalog with
over a quarter million events (2.5 times the standard regional catalog in our focus region). The updated
processing procedures includes an MCMC sampling approach for estimating hypocenter parameters; this
newly developed method has the benets of providing uncertainty estimates and being able to account for
arbitrarily detailed models of observation and prediction errors. Selected structural features observed in
the catalog are described to familiarize the reader with its general characteristics; however, it is presented
primarily as a data product for use in derivative studies, without detailed analysis or interpretation. The
catalog can be obtained at https://data.mendeley.com/datasets/7ywkdx7c62.
153
Chapter7
Discussion
This thesis aims at addressing key questions about the subsurface geometry of and physical processes
operating in active fault zones of Southern California. We developed an automated processing procedure
to derive earthquake catalogs and models of seismic velocity structure from raw seismic waveform data
and applied variations of this procedure to data from two main study areas: a) the San Jacinto Fault Zone
(SJFZ), and b) the region around the 2019 Ridgecrest earthquake sequence.
The automated processing procedure developed is a general, modular framework that can be modi-
ed and applied to a wide range of tectonic environments and network geometries. Novel technical de-
velopments presented in this thesis, apart from the catalogs and velocity models themselves, include (a)
implementation of a robust numerical solver for the eikonal equation in 3D media, (b) an algorithm used
to detect earthquakes and pick phase arrivals in the Ridgecrest study, (c) two algorithms for locating earth-
quakes in 3D media, and (d) modication of a new traveltime tomography method to adaptively resolve
structural features based on data distribution. The modular nature of the framework allowed integration
of addtional algorithmic advances made by others (e.g., the REAL algorithm used to associate phase picks
with their earthquake sources in Chapters 5 and 6, and the deep-learning model used to detect earthquakes
and pick phase arrivals in Chapter 6).
154
Applying our procedure to data from the SJFZ produced a microseismicity catalog with 2.5 times as
many events as the standard regional catalog. Newly detected events provide import constraints on fault-
zone geometry and reveal previously unobserved structural features.
Applying our procedure to data from the Ridgecrest region produced the most comprehensive cata-
log currently available for the Ridgecrest aftershock sequence. From the aftershock catalog, we derive
high-resolution models of seismic velocity structure, and jointly interpreting the spatial distribution of
aftershocks, seismic velocity structure, and surface geology leads us to retrospective hypotheses for ex-
plaining the cessation of theM
w
7.1 mainshock.
Catalogs for both the SJFZ and Ridgecrest regions and seismic velocity models for Ridgecrest will
benet future studies in these regions. The processing procedure also provides means of gaining new
insights into other fault zones using existing data and is currently be applied to data from the Dead Sea
Tranform Fault and Banda Arc Subduction Zone regions.
155
Bibliography
N. Abolfathian, P. Martínez-Garzón, and Y. Ben-Zion (2018). “Spatiotemporal Variations of Stress and
Strain Parameters in the San Jacinto Fault Zone”. In: Pure and Applied Geophysics 176, pp. 1145–1168.
doi: 10.1007/s00024- 018- 2055- y.
Afnimar and K. Koketsu (2000). “Finite dierence traveltime calculation for head waves travelling along
an irregular interface”. In: Geophysical Journal International 143.3, pp. 729–734.doi:
10.1046/j.1365- 246X.2000.00269.x.
H. Akaike (1974). “A new look at the statistical model identication”. In: IEEE Transactions on Automatic
Control 19.6, pp. 716–723.doi: 10.1109/TAC.1974.1100705.
Albuquerque Seismological Laboratory and United States Geological Survery (1980).USGeologicalSurvey
Networks. International Federation of Digital Seismograph Networks.doi: 10.7914/SN/GS.
T. Alkhalifah and S. Fomel (2001). “Implementing the Fast Marching Eikonal Solver: Spherical Versus
Cartesian Coordinates”. In: Geophysical Prospecting 49.2, pp. 165–178.doi:
10.1046/j.1365- 2478.2001.00245.x.
A. Allam and Y. Ben-Zion (2012). “Seismic velocity structures in the southern California plate-boundary
environment from double-dierence tomography”. In: Geophysical Journal International 190.2,
pp. 1181–1196.doi: 10.1111/j.1365- 246X.2012.05544.x.
A. Allam, Y. Ben-Zion, I. Kurzon, and F. Vernon (2014a). “Seismic velocity structure in the Hot Springs
and Trifurcation areas of the San Jacinto fault zone, California, from double-dierence tomography”.
In: Geophysical Journal International 198.2, pp. 978–999.doi: 10.1093/gji/ggu176.
A. Allam, Y. Ben-Zion, and Z. Peng (2014b). “Seismic Imaging of a Bimaterial Interface Along the
Hayward Fault, CA, with Fault Zone Head Waves and Direct P Arrivals”. In: Pure and Applied
Geophysics 171.11, pp. 2993–3011.doi: 10.1007/s00024- 014- 0784- 0.
R. Allen (1982). “Automatic phase pickers: Their present use and future prospects”. In: Bulletin of the
Seismological Society of America 72.6, pp. 173–186.
W. D. Barnhart, G. P. Hayes, and R. D. Gold (2019). “The July 2019 Ridgecrest, California, Earthquake
Sequence: Kinematics of Slip and Stressing in Cross-Fault Ruptures”. In: Geophysical Research Letters
46.21, pp. 11859–11867.doi: 10.1029/2019GL084741.
156
Y. Ben-Zion (1989). “The response of two joined quarter spaces to SH line sources located at the material
discontinuity interface”. In: Geophysical Journal International 98.2, pp. 213–222.
— (2003). “Appendix 2 Key formulas in earthquake seismology”. In: International Geophysics. Vol. 81.
PART B, pp. 1857–1875.doi: 10.1016/S0074- 6142(03)80304- 2.
— (2008). “Collective behavior of earthquakes and faults: Continuum-discrete transitions, progressive
evolutionary changes, and dierent dynamic regimes”. In: Reviews of Geophysics 46.4, pp. 1–70.doi:
10.1029/2008RG000260.
Y. Ben-Zion and V. Lyakhovsky (2006). “Analysis of aftershocks in a lithospheric model with seismogenic
zone governed by damage rheology”. In: Geophysical Journal International 165.1, pp. 197–210.doi:
10.1111/j.1365- 246X.2006.02878.x.
Y. Ben-Zion and I. Zaliapin (2019). “Spatial variations of rock damage production by earthquakes in
southern California”. In: Earth and Planetary Science Letters 512, pp. 184–193.doi:
10.1016/j.epsl.2019.02.006.
M. Beyreuther, R. Barsch, L. Krischer, T. Megies, Y. Behr, and J. Wassermann (2010). “ObsPy: A Python
Toolbox for Seismology”. In: Seismological Research Letters 81.3, pp. 530–533.doi:
10.1785/gssrl.81.3.530.
M. L. Blanpied, D. A. Lockner, and J. D. Byerlee (1991). “Fault stability inferred from granite sliding
experiments at hydrothermal conditions”. In: Geophysical Research Letters 18.4, pp. 609–612.doi:
10.1029/91GL00469.
M. Bohnho, C. Wollin, D. Domigall, L. Küperkoch, P. Martínez-Garzón, G. Kwiatek, and P. E. Malin
(2017). “Repeating Marmara Sea earthquakes: Indication for fault creep”. In: Geophysical Journal
International 210.1, pp. 332–339.doi: 10.1093/gji/ggx169.
S. J. Brandenberg, J. P. Stewart, P. Wang, C. C. Nweke, K. Hudson, C. A. Goulet, X. Meng, C. A. Davis,
S. K. Ahdi, M. B. Hudson, A. Donnellan, G. Lyzenga, M. Pierce, J. Wang, M. A. Winters, M.-P. Delisle,
J. Lucey, Y. Kim, T. W. Gallien, A. Lyda, J. S. Yeung, O. Issa, T. Buckreis, and Z. Yi (2020). “Ground
Deformation Data from GEER Investigations of Ridgecrest Earthquake Sequence”. In: Seismological
Research Letters 91.4, pp. 2024–2034.doi: 10.1785/0220190291.
G. B. Brietzke and Y. Ben-Zion (2006). “Examining tendencies of in-plane rupture to migrate to material
interfaces”. In: Geophysical Journal International 167.2, pp. 807–819.doi:
10.1111/j.1365- 246X.2006.03137.x.
S. Buske and U. Kastner (2004). “Ecient and accurate computation of seismic traveltimes and
amplitudes”. In: Geophysical Prospecting 52.4, pp. 313–322.doi: 10.1111/j.1365- 2478.2004.00417.x.
California Institute of Technology and United States Geological Survey Pasadena (1926). Southern
California Seismic Network.doi: 10.7914/SN/CI.
157
R. D. Catchings, M. R. Goldman, J. H. Steidl, J. H. Chan, A. A. Allam, C. J. Criley, Z. Ma,
D. S. Langermann, G. J. Huddleston, A. T. McEvilly, D. D. Mongovin, E. M. Berg, and Y. Ben-Zion
(2020). “Nodal Seismograph Recordings of the 2019 Ridgecrest Earthquake Sequence”. In:
Seismological Research Letters 91.6, pp. 3622–3633.doi: 10.1785/0220200203.
Y. Cheng, Z. E. Ross, and Y. Ben-Zion (2018). “Diverse Volumetric Faulting Patterns in the San Jacinto
Fault Zone”. In: Journal of Geophysical Research: Solid Earth 123.6, pp. 5068–5081.doi:
10.1029/2017JB015408.
N. I. Christensen (1996). “Poisson’s ratio and crustal seismology”. In: Journal of Geophysical Research:
Solid Earth 101.B2, pp. 3139–3156.doi: 10.1029/95JB03446.
E. S. Cochran, E. Wolin, D. E. McNamara, A. Yong, D. Wilson, M. Alvarez, N. van der Elst, A. McClain,
and J. Steidl (2020). “The U.S. Geological Survey’s Rapid Seismic Array Deployment for the 2019
Ridgecrest Earthquake Sequence”. In: Seismological Research Letters 91.4, pp. 1952–1960.doi:
10.1785/0220190296.
M. de Kool, N. Rawlinson, and M. Sambridge (2006). “A practical grid-based method for tracking multiple
refraction and reection phases in three-dimensional heterogeneous media”. In: Geophysical Journal
International 167.1, pp. 253–270.doi: 10.1111/j.1365- 246X.2006.03078.x.
E. W. Dijkstra (1959). “A note on two problems in connexion with graphs”. In: Numerische Mathematik
1.1, pp. 269–271.doi: 10.1007/BF01386390.
A. Donnellan, G. Lyzenga, A. Ansar, C. Goulet, J. Wang, and M. Pierce (2020). “Targeted High-Resolution
Structure from Motion Observations over the Mw 6.4 and 7.1 Ruptures of the Ridgecrest Earthquake
Sequence”. In: Seismological Research Letters 91.4, pp. 2087–2095.doi: 10.1785/0220190274.
C. B. DuRoss, R. D. Gold, T. E. Dawson, K. M. Scharer, K. J. Kendrick, S. O. Akciz, S. J. Angster,
J. Bachhuber, S. Bacon, S. E. K. Bennett, L. Blair, B. A. Brooks, T. Bullard, W. P. Burgess, C. Chupik,
M. DeFrisco, J. Delano, J. F. Dolan, E. Frost, N. Graehl, E. K. Haddon, A. E. Hatem, J. L. Hernandez,
C. Hitchcock, K. Hudnut, J. Thompson Jobe, R. Koehler, O. Kozaci, T. Ladinsky, C. Madugo,
D. S. McPhillips, C. Milliner, A. Morelan, B. Olson, J. Patton, B. Philibosian, A. J. Pickering, I. Pierce,
D. J. Ponti, G. Seitz, E. Spangler, B. Swanson, K. Thomas, J. Treiman, F. Valencia, A. Williams, and
R. Zinke (2020). “Surface Displacement Distributions for the July 2019 Ridgecrest, California,
Earthquake Ruptures”. In: Bulletin of the Seismological Society of America 110.4, pp. 1400–1418.doi:
10.1785/0120200058.
H Fang, R. D. van der Hilst, M. V. de Hoop, K. Kothari, S. Gupta, and I. Dokmanić (2020). “Parsimonious
Seismic Tomography with Poisson Voronoi Projections: Methodology and Validation”. In:
Seismological Research Letters 91.1, pp. 343–355.doi: 10.1785/0220190141.
H. Fang, H. Zhang, H. Yao, A. Allam, D. Zigone, Y. Ben-Zion, C. Thurber, and R. van der Hilst (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”. In: Journal of Geophysical Research:
Solid Earth 121.5, pp. 3557–3569.doi: 10.1002/2015JB012702.
Y. Fialko (2006). “Interseismic strain accumulation and the earthquake potential on the southern San
Andreas fault system”. In: Nature 441.7096, pp. 968–971.issn: 0028-0836.doi: 10.1038/nature04797.
158
E. J. Fielding, Z. Liu, O. L. Stephenson, M. Zhong, C. Liang, A. Moore, S. Yun, and M. Simons (2020).
“Surface Deformation Related to the 2019 Mw 7.1 and 6.4 Ridgecrest Earthquakes in California from
GPS, SAR Interferometry, and SAR Pixel Osets”. In: Seismological Research Letters 91.4,
pp. 2035–2046.doi: 10.1785/0220190302.
M. Floyd, G. Funning, Y. Fialko, R. Terry, and T. Herring (2020). “Survey and Continuous GNSS in the
Vicinity of the July 2019 Ridgecrest Earthquakes”. In: Seismological Research Letters 91.4,
pp. 2047–2054.doi: 10.1785/0220190324.
D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman (2013). “emcee : The MCMC Hammer”. In:
Publications of the Astronomical Society of the Pacic 125.925, pp. 306–312.doi: 10.1086/670067.
W. B. Frank and R. E. Abercrombie (2018). “Adapting the matched-lter search to a wide-aperture
network: An aftershock sequence and an earthquake swarm in connecticut”. In: Bulletin of the
Seismological Society of America 108.1, pp. 524–532.doi: 10.1785/0120170190.
S. J. Gibbons and F. Ringdal (2006). “The detection of low magnitude seismic events using array-based
waveform correlation”. In: Geophysical Journal International 165.1, pp. 149–166.doi:
10.1111/j.1365- 246X.2006.02865.x.
T. C. Hanks and H. Kanamori (1979). “A moment magnitude scale”. In: Journal of Geophysical Research
84.B5, p. 2348.doi: 10.1029/JB084iB05p02348.
C. R. Harris, K. Millman Jarrod, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser,
J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane,
J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser,
H. Abbasi, C. Gohlke, and T. E. Oliphant (2020). “Array programming with NumPy”. In: Nature
585.7825, pp. 357–362.doi: 10.1038/s41586- 020- 2649- 2.
A. Harten and S. Osher (1987). “Uniformly High-Order Accurate Nonoscillatory Schemes. I”. In: SIAM
Journal on Numerical Analysis 24.2, pp. 279–309.doi: 10.1137/0724022.
E. Hauksson (2000). “Crustal structure and seismicity distribution adjacent to the Pacic and North
America plate boundary in southern California”. In: Journal of Geophysical Research: Solid Earth
105.B6, pp. 13875–13903.doi: 10.1029/2000JB900016.
E. Hauksson and P. Shearer (2005). “Southern California hypocenter relocation with waveform
cross-correlation, part 1: Results using the double-dierence method”. In: Bulletin of the Seismological
Society of America 95.3, pp. 896–903.issn: 00371106.doi: 10.1785/0120040167.
E. Hauksson, W. Yang, and P. M. Shearer (2012). “Waveform Relocated Earthquake Catalog for Southern
California (1981 to June 2011)”. In: Bulletin of the Seismological Society of America 102.5,
pp. 2239–2244.doi: 10.1785/0120120010.
G. Hillers, Y. Ben-Zion, and P. M. Mai (2006). “Seismicity on a fault controlled by rate- and
state-dependent friction with spatial variations of the critical slip distance”. In:JournalofGeophysical
Research: Solid Earth 111.1, pp. 1–23.doi: 10.1029/2005JB003859.
159
J. A. Hole and B. C. Zelt (1995). “3-D nite-dierence reection travel times”. In: Geophysical Journal
International 121.2, pp. 427–434.doi: 10.1111/j.1365- 246X.1995.tb05723.x.
K. W. Hudnut, B. A. Brooks, K. Scharer, J. L. Hernandez, T. E. Dawson, M. E. Oskin, J. Ramon
Arrowsmith, C. A. Goulet, K. Blake, M. L. Boggs, S. Bork, C. L. Glennie, J. C. Fernandez-Diaz,
A. Singhania, D. Hauser, and S. Sorhus (2020). “Airborne Lidar and Electro-Optical Imagery along
Surface Ruptures of the 2019 Ridgecrest Earthquake Sequence, Southern California”. In: Seismological
Research Letters 91.4, pp. 2096–2107.doi: 10.1785/0220190338.
J. D. Hunter (2007). “Matplotlib: A 2D Graphics Environment”. In: Comput. Sci. Eng. 9.3, pp. 90–95.doi:
10.1109/MCSE.2007.55.
K. Hutton, J. Woessner, and E. Hauksson (2010). “Earthquake Monitoring in Southern California for
Seventy-Seven Years (1932-2008)”. In: Bulletin of the Seismological Society of America 100.2,
pp. 423–446.doi: 10.1785/0120090130.
L. K. Hutton and D. M. Boore (1987). “The Ml Scale in Southern California”. In: Bulletin of the
Seismological Society of America 77.6, pp. 2074–2094.
A. Inbal, J. P. Ampuero, and J. P. Avouac (2017). “Locally and remotely triggered aseismic slip on the
central San Jacinto Fault near Anza, CA, from joint inversion of seismicity and strainmeter data”. In:
Journal of Geophysical Research: Solid Earth 122.4, pp. 3033–3061.doi: 10.1002/2016JB013499.
A. Ito (1990). “Earthquake swarm activity revealed from high-resolution relative hypocenters - clustering
of microearthquakes”. In: Tectonophysics 175.1-3, pp. 47–66.doi: 10.1016/0040- 1951(90)90129- V.
J. Jiang and N. Lapusta (2017). “Connecting depth limits of interseismic locking, microseismicity, and
large earthquakes in models of long-term fault slip”. In: Journal of Geophysical Research: Solid Earth
122.8, pp. 6491–6523.doi: 10.1002/2017JB014030.
Z. Jin and Y. Fialko (2020). “Finite Slip Models of the 2019 Ridgecrest Earthquake Sequence Constrained
by Space Geodetic Data and Aftershock Locations”. In: Bulleting of the Seismological Sociecty of
America 110.4, pp. 1660–1679.doi: 10.1785/0120200060.
E. Jones, T. Oliphant, and P. Peterson (2001). SciPy: Open source scientic tools for Python – Reference
Guide, 0.14.url: http://www.scipy.org/.
A. Juarez and Y. Ben-Zion (2020). “Eects of Shallow-Velocity Reductions on 3D Propagation of Seismic
Waves”. In: Seismological Research Letters 91.6, pp. 3313–3322.doi: 10.1785/0220200183.
B. R. Julian and D. Gubbins (1977). “Three-dimensional seismic ray tracing”. In: Journal of Geophysics 43,
pp. 95–114.
A. Jurkevics (1988). “Polarization Analysis of Three-Component Array Data”. In: Bulletin of the
Seismological Society of America 78.5, pp. 1725–1743.
B. L. N. Kennett and E. R. Engdahl (1991). “Traveltimes for global earthquake location and phase
identication”. In: Geophysical Journal International 105, pp. 429–465.
160
B. L. N. Kennett, E. R. Engdahl, and R. Buland (1995). “Constraints on Seismic Velocities in the Earth from
Traveltimes”. In: Geophysical Journal International 122.1, pp. 108–124.doi:
10.1111/j.1365- 246X.1995.tb03540.x.
S. Kim and R. Cook (1999). “3-D traveltime computation using second-order ENO scheme”. In:Geophysics
64.6, pp. 1867–1876.doi: 10.1190/1.1444693.
D. L. Kohlstedt, B. Evans, and S. J. Mackwell (2004). “Strength of the lithosphere: Constraints imposed by
laboratory experiments”. In: Journal of Geophysical Research: Solid Earth 100.B9, pp. 17587–17602.
doi: 10.1029/95jb01460.
L. Krischer, J. Smith, W. Lei, M. Lefebvre, Y. Ruan, E. S. de Andrade, N. Podhorszki, E. Bozdağ, and
J. Tromp (2016). “An Adaptable Seismic Data Format”. In: Geophysical Journal International 207.2,
pp. 1003–1011.doi: 10.1093/gji/ggw319.
E.-J. Lee, P. Chen, T. H. Jordan, P. B. Maechling, M. A. M. Denolle, and G. C. Beroza (2014). “Full-3-D
tomography for crustal structure in Southern California based on the scattering-integral and the
adjoint-waveeld methods”. In:JournalofGeophysicalResearch:SolidEarth 119.8, pp. 6421–6451.doi:
10.1002/2014JB011346.
E.-J. Lee, W.-Y. Liao, D. Mu, W. Wang, and P. Chen (2020a). “GPU-Accelerated Automatic Microseismic
Monitoring Algorithm (GAMMA) and Its Application to the 2019 Ridgecrest Earthquake Sequence”.
In: Seismological Research Letters 91.4, pp. 2062–2074.doi: 10.1785/0220190323.
E.-J. Lee, D. Mu, W. Wang, and P. Chen (2020b). “Weighted Template-Matching Algorithm (WTMA) for
Improved Foreshock Detection of the 2019 Ridgecrest Earthquake Sequence”. In: Bulletin of the
Seismological Society of America 110.4, pp. 1832–1844.doi: 10.1785/0120200020.
J. L. Lewis, S. M. Day, H. Magistrale, J. Eakins, and F. Vernon (Apr. 2000). “Regional crustal thickness
variations of the Peninsular Ranges, southern California”. In: Geology 28.4, pp. 303–306.doi:
10.1130/0091- 7613(2000)028< 0303:RCTVOT> 2.3.CO;2.
C. Li, R. D. van der Hilst, E. R. Engdahl, and S. Burdick (2008). “A New Global Model for P Wave Speed
Variations in Earth’s Mantle”. In:Geochemistry,Geophysics,Geosystems 9.5.doi: 10.1029/2007GC001806.
F.-C. Lin, D. Li, R. W. Clayton, and D. Hollis (2013). “High-resolution 3D shallow crustal structure in
Long Beach, California: Application of ambient noise tomography on a dense seismic array”. In:
Geophysics 78.4.issn: 0016-8033.doi: 10.1190/geo2012- 0453.1.
G. Lin (2020). “Waveform Cross-Correlation Relocation and Focal Mechanisms for the 2019 Ridgecrest
Earthquake Sequence”. In: SeismologicalResearchLetters 91.4, pp. 2055–2061.doi: 10.1785/0220190277.
G. Lin and P. M. Shearer (Sept. 2009). “Evidence for water-lled cracks in earthquake source regions”. In:
Geophysical Research Letters 36.17, p. L17315.issn: 0094-8276.doi: 10.1029/2009GL039098.
G. Lin, P. M. Shearer, and E. Hauksson (2007). “Applying a three-dimensional velocity model, waveform
cross correlation, and cluster analysis to locate southern California seismicity from 1981 to 2005”. In:
Journal of Geophysical Research 112.B12, B12309.issn: 0148-0227.doi: 10.1029/2007JB004986.
161
M. Liu, M. Zhang, W. Zhu, W. L. Ellsworth, and H. Li (2020). “Rapid Characterization of the July 2019
Ridgecrest, California, Earthquake Sequence From Raw Seismic Data Using Machine-Learning Phase
Picker”. In: Geophysical Research Letters 47.4, pp. 1–9.doi: 10.1029/2019GL086189.
X.-D. Liu, S. Osher, and T. Chan (1994). “Weighted Essentially Non-oscillatory Schemes”. In: Journal of
Computational Physics 115.1, pp. 200–212.doi: 10.1002/fld.3889.
A. Lomax (2020). “Absolute Location of 2019 Ridgecrest Seismicity Reveals a Shallow Mw 7.1 Hypocenter,
Migrating and Pulsing Mw 7.1 Foreshocks, and Duplex Mw 6.4 Ruptures”. In: Bulleting of the
Seismological Society of America 110.4, pp. 1845–1858.doi: 10.1785/0120200006.
A. Lomax, A. Michelini, and A. Curtis (2009). Earthquake Location, Direct, Global-Search Methods.
Springer Berlin Heidelberg, pp. 1–33.doi: 10.1007/978- 3- 642- 27737- 5.
J. N. Louie (2001). “Faster, Better: Shear-Wave Velocity to 100 Meters Depth from Refraction Microtremor
Arrays”. In: Bulletin of the Seismological Society of America 91.2, pp. 347–364.doi: 10.1785/0120000098.
N. Maeda (1985). “A Method for Reading and Checking Phase Time in Auto-Processing System of Seismic
Wave Data”. In: Zisin (Journal of the Seismological Society of Japan. 2nd ser.) 38.3, pp. 365–379.doi:
10.4294/zisin1948.38.3_365.
Y. Magen, A. Ziv, A. Inbal, G. Baer, and J. Hollingsworth (2020). “Fault Rerupture during the July 2019
Ridgecrest Earthquake Pair from Joint Slip Inversion of InSAR, Optical Imagery, and GPS”. In:
Bulletin of the Seismological Society of America 110.4, pp. 1627–1643.doi: 10.1785/0120200024.
H. Magistrale and H.-W. Zhou (1996). “Lithologic Control of the Depth of Earthquakes in Southern
California”. In: Science 273.5275, pp. 639–642.doi: 10.1126/science.273.5275.639.
G. S. Mattioli, D. A. Phillips, K. M. Hodgkinson, C. Walls, D. J. Mencin, B. A. Bartel, D. J. Charlevoix,
C. Crosby, M. J. Gottlieb, B. Henderson, W. Johnson, D. Maggert, D. Mann, C. M. Meertens,
J. Normandeau, J. Pettit, C. M. Puskas, L. Rowan, C. Sievers, and A. Zaino (2020). “The GAGE Data
and Field Response to the 2019 Ridgecrest Earthquake Sequence”. In: Seismological Research Letters
91.4, pp. 2075–2086.doi: 10.1785/0220190283.
J. McGuire and Y. Ben-Zion (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”. In:
Geophysical Journal International 163.1, pp. 152–164.doi: 10.1111/j.1365- 246X.2005.02703.x.
W. McKinney (2010). “Data Structures for Statistical Computing in Python”. In: Proceedings of the 9th
Python in Science Conference. Ed. by S. van der Walt and J. Millman, pp. 56–61.doi:
10.25080/Majora- 92bf1922- 00a.
M.-A. Meier, Z. E. Ross, A. Ramachandran, A. Balakrishna, S. Nair, P. Kundzicz, Z. Li, J. Andrews,
E. Hauksson, and Y. Yue (2019). “Reliable Real-Time Seismic Signal/Noise Discrimination With
Machine Learning”. In: Journal of Geophysical Research: Solid Earth 124.1, pp. 788–800.doi:
10.1029/2018JB016661.
162
S. M. Mousavi, W. L. Ellsworth, W. Zhu, L. Y. Chuang, and G. C. Beroza (2020). “Earthquake
transformer—an attentive deep-learning model for simultaneous earthquake detection and phase
picking”. In: Nature Communications 11.1, p. 3952.doi: 10.1038/s41467- 020- 17591- w.
S. M. Mousavi, Y. Sheng, W. Zhu, and G. C. Beroza (2019). “STanford EArthquake Dataset (STEAD): A
Global Data Set of Seismic Signals for AI”. In: IEEE Access 7, pp. 179464–179476.doi:
10.1109/ACCESS.2019.2947848.
R. Nadeau, W. Foxall, P. A. Johnson, T. V. McEvilly, and M. Antolik (2003). “Seismological studies at
Parkeld III: microearthquake clusters in the study of fault-zone dynamics”. In: International Journal
of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 31.6, p. 271.doi:
10.1016/0148- 9062(94)90077- 9.
Y. Ogata and K. Katsura (1993). “Analysis of temporal and spatial heterogeneity of magnitude frequency
distribution inferred from earthquake catalogues”. In: Geophysical Journal International 113.3,
pp. 727–738.doi: 10.1111/j.1365- 246X.1993.tb04663.x.
S. Osher and J. A. Sethian (1988). “Fronts propagating with curvature-dependent speed: Algorithms based
on Hamilton-Jacobi formulations”. In: Journal of Computational Physics 79.1, pp. 12–49.doi:
10.1016/0021- 9991(88)90002- 2.
Y. Ozakin and Y. Ben-Zion (2015). “Systematic Receiver Function Analysis of the Moho Geometry in the
Southern California Plate-Boundary Region”. In: Pure and Applied Geophysics 172.5, pp. 1167–1184.
doi: 10.1007/s00024- 014- 0924- 6.
Z. Peng and Y. Ben-Zion (2005). “Spatiotemporal variations of crustal anisotropy from similar events in
aftershocks of the 1999 M7.4 Izmit and M7.1 Duzce, Turkey, earthquake sequences”. In: Geophysical
Journal International 160.3, pp. 1027–1043.doi: 10.1111/j.1365- 246X.2005.02569.x.
F. Perez and B. E. Granger (2007). “IPython: A System for Interactive Scientic Computing”. In:
Computing in Science and Engineering 9.3, pp. 21–29.doi: 10.1109/MCSE.2007.53.
I. Pierce, A. Williams, R. D. Koehler, and C. Chupik (2020). “High-Resolution Structure-From-Motion
Models and Orthophotos of the Southern Sections of the 2019 Mw 7.1 and 6.4 Ridgecrest Earthquakes
Surface Ruptures”. In: Seismological Research Letters 91.4, pp. 2124–2126.doi: 10.1785/0220190289.
A. Plesch, J. H. Shaw, Z. E. Ross, and E. Hauksson (2020). “Detailed 3D Fault Representations for the 2019
Ridgecrest, California, Earthquake Sequence”. In: Bulletin of the Seismological Society of America
110.4, pp. 1818–1831.doi: 10.1785/0120200053.
P. Podvin and I. Lecomte (1991). “Finite dierence computation of traveltimes in very contrasted velocity
models: a massively parallel approach and its associated tools”. In: Geophysical Journal International
105.1, pp. 271–284.doi: 10.1111/j.1365- 246X.1991.tb03461.x.
163
D. J. Ponti, J. L. Blair, C. M. Rosa, K. Thomas, A. J. Pickering, S. Akciz, S. Angster, J.-P. Avouac,
J. Bachhuber, S. Bacon, N. Barth, S. Bennett, K. Blake, S. Bork, B. Brooks, T. Bullard, P. Burgess,
C. Chupik, T. Dawson, M. DeFrisco, J. Delano, S. DeLong, J. Dolan, A. Donnellan, C. DuRoss,
T. Ericksen, E. Frost, G. Funning, R. Gold, N. Graehl, C. Gutierrez, E. Haddon, A. Hatem, J. Helms,
J. Hernandez, C. Hitchcock, P. Holland, K. Hudnut, K. Kendrick, R. Koehler, O. Kozaci, T. Ladinsky,
R. Leeper, C. Madugo, M. Mareschal, J. McDonald, D. McPhillips, C. Milliner, D. Mongovin,
A. Morelan, S. Nale, J. Nevitt, M. O’Neal, B. Olson, M. Oskin, S. Padilla, J. Patton, B. Philibosian,
I. Pierce, C. Pridmore, N. Roth, D. Sandwell, K. Scharer, G. Seitz, D. Singleton, B. Smith-Konter,
E. Spangler, B. Swanson, J. Thompson Jobe, J. Treiman, F. Valencia, J. Vanderwal, A. Williams, X. Xu,
J. Zachariasen, J. Zimmerman, and R. Zinke (2020). “Documentation of Surface Fault Rupture and
Ground-Deformation Features Produced by the 4 and 5 July 2019 Mw 6.4 and Mw 7.1 Ridgecrest
Earthquake Sequence”. In: SeismologicalResearchLetters 91.5, pp. 2942–2959.doi: 10.1785/0220190322.
J. Qian and W. W. Symes (2002a). “An adaptive nite-dierence method for traveltimes and amplitudes”.
In: Geophysics 67.1, pp. 167–176.doi: 10.1190/1.1451472.
— (2002b). “An adaptive nite-dierence method for traveltimes and amplitudes”. In: Geophysics 67.1,
pp. 167–176.doi: 10.1190/1.1451472.
F. Qin, Y. Luo, K. B. Olsen, W. Cai, and G. T. Schuster (1992). “Finite-dierence solution of the eikonal
equation along expanding wavefronts”. In: Geophysics 57.3, pp. 478–487.doi: 10.1190/1.1443263.
L. Qin, Y. Ben-Zion, H. Qiu, P.-E. Share, Z. E. Ross, and F. L. Vernon (2018). “Internal structure of the San
Jacinto fault zone in the trifurcation area southeast of Anza , California , from data of dense seismic
arrays”. In: 213.1, pp. 98–114.doi: 10.1093/gji/ggx540/4757070.
H. Qiu, Y. Ben-Zion, R. D. Catchings, M. R. Goldman, A. A. Allam, and J. H. Steidl (2021). “Detailed
seismic imaging of the Mw 7.1 Ridgecrest earthquake rupture zone from data recorded by dense
linear arrays”. In: Journal of Geophysical Research: Solid Earth In press.
H. Qiu, Y. Ben-Zion, Z. E. Ross, P.-E. Share, and F. L. Vernon (2017). “Internal structure of the San Jacinto
fault zone at Jackass Flat from data recorded by a dense linear array”. In: Geophysical Journal
International 209.3, pp. 1369–1388.doi: 10.1093/gji/ggx096.
N. Rawlinson, J. Hauser, and M. Sambridge (2008). “Seismic ray tracing and wavefront tracking in
laterally heterogeneous media”. In: Advances in Geophysics 49.07, pp. 203–273.doi:
10.1016/S0065- 2687(07)49003- 3.
N. Rawlinson and M. Sambridge (2004a). “Multiple reection and transmission phases in complex layered
media using a multistage fast marching method”. In: Geophysics 69.5, pp. 1338–1350.doi:
10.1190/1.1801950.
— (2004b). “Wave Front Evolution in Strongly Heterogeneous Layered Media Using the Fast Marching
Method”. In: Geophysical Journal International 156.3, pp. 631–647.doi:
10.1111/j.1365- 246X.2004.02153.x.
M. Reshef and D. Koslo (1986). “Migration of common-shot gathers”. In: Geophysics 51.2, pp. 324–331.
doi: 10.1190/1.1442091.
164
J. R. Rice (1993). “Spatio-temporal Complexity of Slip on a Fault Rate- dependent friction”. In: Journal of
Geophysical Research 98.93, pp. 9885–9907.doi: 10.1029/93JB00191.
C. F. Richter (1935). “An instrumental earthquake magnitude scale”. In: Bulletin of the Seismological
Society of America 25, pp. 1–32.
T. K. Rockwell, T. E. Dawson, J. Young Ben-Horin, and G. Seitz (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”. In: Pure and Applied Geophysics 172.5,
pp. 1143–1165.issn: 0033-4553.doi: 10.1007/s00024- 014- 0955- z.
F. Rolandone, R. Bürgmann, and R. M. Nadeau (2004). “The evolution of the seismic-aseismic transition
during the earthquake cycle: Constraints from the time-dependent depth distribution of aftershocks”.
In: Geophysical Research Letters 31.23, pp. 1–4.doi: 10.1029/2004GL021379.
Z. E. Ross and Y. Ben-Zion (2014). “Automatic picking of direct P, S seismic phases and fault zone head
waves”. In: Geophysical Journal International 199.1, pp. 368–381.doi: 10.1093/gji/ggu267.
Z. E. Ross, Y. Ben-Zion, M. C. A. White, and F. L. Vernon (2016a). “Analysis of earthquake body wave
spectra for potency and magnitude values: implications for magnitude scaling relations”. In:
Geophysical Journal International 207.2, pp. 1158–1164.doi: 10.1093/gji/ggw327.
Z. E. Ross, E. S. Cochran, D. T. Trugman, and J. D. Smith (2020). “3D fault architecture controls the
dynamism of earthquake swarms”. In: Science 368.6497, pp. 1357–1361.issn: 0036-8075.doi:
10.1126/science.abb0779.
Z. E. Ross, E. Hauksson, and Y. Ben-Zion (2017). “Abundant o-fault seismicity and orthogonal structures
in the San Jacinto fault zone”. In: Science Advances 3.3, e1601946.doi: 10.1126/sciadv.1601946.
Z. E. Ross, D. T. Trugman, and P. M. Shearer (2019). “Searching for Hidden Earthquakes in Southern
California”. In: Science.
Z. E. Ross, M. C. A. White, F. L. Vernon, and Y. Ben-Zion (2016b). “An Improved Algorithm for Real-Time
S -Wave Picking with Application to the (Augmented) ANZA Network in Southern California”. In:
Bulletin of the Seismological Society of America 106.5, pp. 2013–2022.doi: 10.1785/0120150230.
E. Rouy and A. Tourin (1992). “A Viscosity Solutions Approach to Shape-From-Shading”. In: SIAM
Journal on Numerical Analysis 29.3, pp. 867–884.doi: 10.1137/0729053.
A. M. Rubin and D. Gillard (2004). “Aftershock asymmetry/rupture directivity among central San
Andreas fault microearthquakes”. In: Journal of Geophysical Research: Solid Earth 105.B8,
pp. 19095–19109.doi: 10.1029/2000jb900129.
C. O. Sanders and H Kanamori (1984). “A seismotectonic analysis of the Anza seismic gap, San Jacinto
fault zone, southern California”. In: Journal of Geophysical Research 89.B7, pp. 5873–5890.issn:
01480227.doi: 10.1029/JB089iB07p05873.
165
D. P. Scha, G. H. R. Bokelmann, G. C. Beroza, F. Waldhauser, and W. L. Ellsworth (2002).
“High-resolution image of Calaveras Fault seismicity”. In: Journal of Geophysical Research: Solid Earth
107.B9, ESE 5–1–ESE 5–16.doi: 10.1029/2001JB000633.
W. A. Schneider (1978). “Integral Formulation for Migration in Two and Three Dimensions”. In:
Geophysics 43.1, pp. 49–76.doi: 10.1190/1.1440828.
W. A. Schneider, K. A. Ranzinger, A. H. Balch, and C. Kruse (1992). “A dynamic programming approach
to rst arrival traveltime computation in media with arbitrarily distributed velocities”. In: Geophysics
57.1, pp. 39–50.doi: 10.1190/1.1443187.
G. T. Schuster (2009). Seismic Interferometry. Cambridge: Cambridge University Press.isbn:
9780511581557.doi: 10.1017/CBO9780511581557.
D. W. Scott (2015). Multivariate Density Estimation: Theory, Practice, and Visualization. Wiley Series in
Probability and Statistics. Wiley.isbn: 9780471697558.
A. Sescu (2015). “Numerical anisotropy in nite dierencing”. In: Advances in Dierence Equations 2015.1,
p. 9.doi: 10.1186/s13662- 014- 0343- 0.
J. A. Sethian (1996). “A Fast Marching Level Set Method for Monotonically Advancing Fronts”. In:
Proceedings of the National Academy of Sciences 93.4, pp. 1591–1595.doi: 10.1073/pnas.93.4.1591.
— (1999). “Fast Marching Methods”. In: SIAM Review 41.2, pp. 199–235.doi: 10.1137/S0036144598347059.
J. A. Sethian and A. M. Popovici (1999). “3-D Traveltime Computation Using the Fast Marching Method”.
In: Geophysics 64.2, pp. 516–523.doi: 10.1190/1.1444558.
E. Seyhan and J. P. Stewart (2014). “Semi-Empirical Nonlinear Site Amplication from NGA-West2 Data
and Simulations”. In: Earthquake Spectra 30.3, pp. 1241–1256.doi: 10.1193/063013EQS181M.
P.-E. Share, Y. Ben-Zion, Z. E. Ross, H. Qiu, and F. L. Vernon (2017). “Internal structure of the San Jacinto
fault zone at Blackburn Saddle from seismic data of a linear array”. In: Geophysical Journal
International 210.2, pp. 819–832.doi: 10.1093/gji/ggx191.
P.-E. Share, H. Guo, C. H. Thurber, H. Zhang, and Y. Ben-Zion (2018). “Seismic Imaging of the Southern
California Plate Boundary around the South-Central Transverse Ranges Using Double-Dierence
Tomography”. In: Pure and Applied Geophysics 176, pp. 1117–1144.doi: 10.1007/s00024- 018- 2042- 3.
R. V. Sharp (1967). “San Jacinto Fault Zone in the Peninsular Ranges of Southern California”. In:
Geological Society of America 78.6, pp. 705–730.doi: 10.1130/0016- 7606(1967)78.
J. H. Shaw, A. Plesch, C. Tape, M. P. Suess, T. H. Jordan, G. Ely, E. Hauksson, J. Tromp, T. Tanimoto,
R. Graves, K. Olsen, C. Nicholson, P. J. Maechling, C. Rivero, P. Lovely, C. M. Brankman, and
J. Munster (2015). “Unied Structural Representation of the southern California crust and upper
mantle”. In: Earth and Planetary Science Letters 415, pp. 1–15.doi: 10.1016/j.epsl.2015.01.016.
166
P. Shearer, E. Hauksson, and G. Lin (2005). “Southern California hypocenter relocation with waveform
cross-correlation, part 2: Results using source-specic station terms and cluster analysis”. In: Bulletin
of the Seismological Society of America 95.3, pp. 904–915.issn: 00371106.doi: 10.1785/0120040168.
D. R. Shelly (2020). “A High-Resolution Seismic Catalog for the Initial 2019 Ridgecrest Earthquake
Sequence: Foreshocks, Aftershocks, and Faulting Complexity”. In: Seismological Research Letters 91.4,
pp. 1971–1978.doi: 10.1785/0220190309.
D. R. Shelly, W. L. Ellsworth, and D. P. Hill (2016). “Fluid-faulting evolution in high denition:
Connecting fault structure and frequency-magnitude variations during the 2014 Long Valley Caldera,
California, earthquake swarm”. In: Journal of Geophysical Research: Solid Earth 121.3, pp. 1776–1795.
doi: 10.1002/2015JB012719.
T. Shimamoto (2003). “The Mechanics of Earthquakes and Faulting”. In: Journal of Structural Geology.
2nd ed. Vol. 15. 6. New York: Cambridge University Press, pp. 810–811.isbn: 9780511818516.doi:
10.1016/0191- 8141(93)90067- k.
R. Sibson (2002). “Earthquakes and Rock Deformation in Crustal Fault Zones”. In: Annual Review of Earth
and Planetary Sciences 14.1, pp. 149–175.doi: 10.1146/annurev.earth.14.1.149.
R. H. Sibson (1982). “Fault zone medels, heat ow, and the depth distribution of earthquakes in the
continental crust of the United States”. In: Bulletin of the Seismological Society of America 72.1,
pp. 151–163.
P. Small, D. Gill, P. J. Maechling, R. Taborda, S. Callaghan, T. H. Jordan, K. B. Olsen, G. P. Ely, and
C. Goulet (2017). “The SCEC Unied Community Velocity Model Software Framework”. In:
Seismological Research Letters 88.6, pp. 1539–1552.doi: 10.1785/0220170082.
J. D. Smith, Z. E. Ross, K. Azizzadenesheli, and J. B. Muir (2021). “HypoSVI: Hypocenter inversion with
Stein variational inference and Physics Informed Neural Networks”. In: arXiv: 2101.03271.url:
http://arxiv.org/abs/2101.03271.
Southern California Earthquake Center (2013). SCEDC.doi: 10.7909/C3WD3xH1.
J. Steidl, R. Catchings, and A. Allam (2019). RAMP deployment of 3C nodal for July Searles Valley 2019
Earthquake. International Federation of Digital Seismograph Networks.doi: 10.7914/SN/3J\_2019.
R. Storn and K. Price (1997). “Dierential Evolution – A Simple and Ecient Heuristic for Global
Optimization over Continuous Spaces”. In: Journal of Global Optimization 11.4, pp. 341–359.doi:
10.1023/A:1008202821328.
C. H. Thurber and W. L. Ellsworth (1980). “Rapid Solution of Ray Tracing Problems in Heterogeneous
Media”. In: Bulletin of the Seismological Society of America 70.4, pp. 1137–1148.
D. T. Trugman and P. M. Shearer (2017). “GrowClust: A Hierarchical Clustering Algorithm for Relative
Earthquake Relocation, with Application to the Spanish Springs and Sheldon, Nevada, Earthquake
Sequences”. In: Seismological Research Letters 88.2A, pp. 379–391.doi: 10.1785/0220160188.
167
J. Um and C. Thurber (1987). “A Fast Algorithm for Two-Point Seismic Ray Tracing”. In: Bulletin of the
Seismological Society of America 77.3, pp. 972–986.
University of Nevada Reno (1971). Nevada Seismic Network. International Federation of Digital
Seismograph Networks.doi: 10.7914/SN/NN.
— (1980). Southern Great Basin Network. International Federation of Digital Seismograph Networks.doi:
10.7914/SN/SN.
J. van Trier and W. W. Symes (1991). “Upwind nite-dierence calculation of traveltimes”. In: Geophysics
56.6, pp. 812–821.doi: 10.1190/1.1443099.
F. L. Vernon (1982). ANZA Regional Network. San Diego.doi: 10.7914/SN/AZ.
F. L. Vernon and Y. Ben-Zion (2010). San Jacinto Fault Zone Experiment Network.doi:
10.7914/SN/YN{\_}2010.
R. Versteeg (1994). “The Marmousi Experience: Velocity Model Determination on a Synthetic Complex
Data Set”. In: The Leading Edge 13.9, pp. 927–936.doi: 10.1190/1.1437051.
J. E. Vidale (1988). “Finite-Dierence Calculation of Travel Times”. In: Bulletin of the Seismological Society
of America 78.6, pp. 2062–2076.
— (1990). “Finite-Dierence Calculation of Traveltimes in Three Dimensions”. In: Geophysics 55.5,
pp. 521–526.doi: 10.1190/1.1442863.
J. E. Vidale and P. M. Shearer (2006). “A survey of 71 earthquake bursts across southern California:
Exploring the role of pore uid pressure uctuations and aseismic slip as drivers”. In: Journal of
Geophysical Research: Solid Earth 111.5.doi: 10.1029/2005JB004034.
P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski,
P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov,
A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, I. Polat, Y. Feng, E. W. Moore, J. VanderPlas,
D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald,
A. H. Ribeiro, F. Pedregosa, and P. van Mulbregt (2020). “SciPy 1.0: fundamental algorithms for
scientic computing in Python”. In: Nat. Methods 17.3, pp. 261–272.doi: 10.1038/s41592- 019- 0686- 2.
S. Wdowinski (2009). “Deep creep as a cause for the excess seismicity along the San Jacinto fault”. In:
Nature Geoscience 2.12, pp. 882–885.doi: 10.1038/ngeo684.
D. J. White (1989). “Two-Dimensional Seismic Refraction Tomography”. In: Geophysical Journal
International 97.2, pp. 223–245.doi: 10.1111/j.1365- 246X.1989.tb00498.x.
M. C. A. White, Y. Ben-Zion, and F. L. Vernon (2019). “A Detailed Earthquake Catalog for the San Jacinto
Fault-Zone Region in Southern California”. In: Journal of Geophysical Research: Solid Earth 124.7,
pp. 6908–6930.doi: 10.1029/2019JB017641.
168
M. C. A. White, H. Fang, R. D. Catchings, M. R. Goldman, J. H. Steidl, and Y. Ben-Zion (2021). “Detailed
traveltime tomography and seismic catalogue around the 2019 M w7.1 Ridgecrest, California,
earthquake using dense rapid-response seismic data”. In: Geophysical Journal International 227.1,
pp. 204–227.doi: 10.1093/gji/ggab224.
M. C. A. White, H. Fang, N. Nakata, and Y. Ben-Zion (2020). “PyKonal: A Python Package for Solving the
Eikonal Equation in Spherical and Cartesian Coordinates Using the Fast Marching Method”. In:
Seismological Research Letters 91.4, pp. 2378–2389.doi: 10.1785/0220190318.
S. Xu and Y. Ben-Zion (2013). “Numerical and theoretical analyses of in-plane dynamic rupture on a
frictional interface and o-fault yielding patterns at dierent scales”. In: Geophysical Journal
International 193.1, pp. 304–320.doi: 10.1093/gji/ggs105.
M. Zhang, W. L. Ellsworth, and G. C. Beroza (2019). “Rapid Earthquake Association and Location”. In:
Seismological Research Letters 90.6, pp. 2276–2284.doi: 10.1785/0220190052.
Q. Zhang and G. Lin (2014). “Three-dimensional Vp and Vp/Vs models in the Coso geothermal area,
California: Seismic characterization of the magmatic system”. In: Journal of Geophysical Research:
Solid Earth 119.6, pp. 4907–4922.doi: 10.1002/2014JB010992.
H. Zhao (2004). “A fast sweeping method for Eikonal equations”. In: Mathematics of Computation 74.250,
pp. 603–628.doi: 10.1090/S0025- 5718- 04- 01678- 3.
L. Zhu and H. Kanamori (2000). “Moho depth variation in southern California from teleseismic receiver
functions”. In: Journal of Geophysical Research: Solid Earth 105.B2, pp. 2969–2980.doi:
10.1029/1999JB900322.
D. Zigone, Y. Ben-Zion, M. Lehujeur, M. Campillo, G. Hillers, and F. L. Vernon (2019). “Imaging
subsurface structures in the San Jacinto fault zone with high-frequency noise recorded by dense
linear arrays”. In: Geophysical Journal International 217.2, pp. 879–893.doi: 10.1093/gji/ggz069.
169
AppendixA
Detectingandpickingphasearrivals
Bold-face lowercase letters represent vectors throughout this article. Consider some time series, x, with
n samples, x =hx
0
;x
1
;:::;x
n1
i. Then x denotes the mean value of x, x [k] x
k
denotes the sample
with indexk, and x [i :j]hx
i
;x
i+1
;:::;x
j1
i denotes the segment of x with indices running fromi to
j 1. A superscript letter attached to the label of a time series indicates that it corresponds to a particular
data component—viz., Z, N, or E—or phase—viz., P or S. Square brackets [ ] enclose integer indices of time
series, parentheses ( ) enclose function arguments, and curly bracesfg indicate algebraic grouping.
A.1 MeasuringP-wavearrivaltimes
Measuring the arrival times of the seismic phases used to locate events is the rst critical step in the
workow we use to build an earthquake catalog. S-wave arrival times are typically harder to measure
than those for P-waves because they arrive later and are interfered with by the P-wave coda, so we begin
our processing by targeting P-wave arrivals. The existence of a P-wave must rst be detected before its
arrival time can be measured; the following describes our approach to solving the two-part problem of
detecting the existence of and measuring the arrival time of P-waves.
A characteristic function (CF) based on three simple and intuitive statistics—the ratio of the signal
variance in a short time window to the same in a longer window (STA/LTA), the signal kurtosis in a
170
Figure A.1: (a) 3C seismic waveform data with three registered P-wave arrivals (vertical dashed blue lines)
and (b) derived statistics used to compute the (c) characteristic function (CF) and dynamic threshold.
sliding window, and the ratio of the signal variance on the vertical channel to the mean of the same on the
horizontal channels (V/H)—coupled with an adaptive threshold robustly detects P-wave arrivals (Figure
A.1).
The STA/LTA, s
P
, is computed for vertical-component data, z, using short- and long-term window
lengths of 0.25 s and 4 s, respectively:
s
P
i
=
var (z [i 0:25r :i])
var (z [i 4r :i])
; (A.1)
in whichr is the sampling rate and var () is the sample variance of the argument:
var (x)
1
n 1
n1
X
i=0
fx
i
xg
2
: (A.2)
171
The sliding-window kurtosis for the vertical component, k
P
, is computed using a window length of 5
s:
k
P
i
=

4
(z[i 5r :i])
fvar (z [i 5r :i])g
2
; (A.3)
in which
4
() is the fourth central moment of the argument:

4
(x)
1
n 1
n1
X
i=0
fx
i
xg
4
: (A.4)
The V/H ratio, v, is computed using a window length of 0.5 s:
v
i
=
2 var (z [i 0:5r :i])
var (n [i 0:5r :i]) + var (e [i 0:5r :i])
: (A.5)
in which n and e represent North-South and East-West component data, respectively.
Finally, the CF for detecting P-waves, c
P
, is dened as the product of these three statistics:
c
P
i
s
P
i
k
P
i
v
i
; (A.6)
and a threshold, t
P
, is dened as six times the RMS of the CF in the preceding 5 s:
t
P
i
6
v
u
u
t
1
5r
5r
X
j=0
n
c
P
ij
o
2
: (A.7)
The algorithm registers a detection at every index wherec
P
i
>t
P
i
and thins clusters by retaining only
the detection with the highest signal-to-noise ratio (SNR) in each 1.5 s window. Then, the arrival time for
each detected P-wave is measured using the Akaike Information Criterion (AIC; Akaike 1974) by extracting
a 3 s window centered on each detection and registering the sample index,i

, corresponding to the global
minimum of the AIC for that window:
172
i

= arg min
i
(i log (var (z [0 :i])) +fnig log (var (z [i :n])))
(A.8)
A.2 MeasuringS-wavearrivaltimes
Having determined the onset time of all candidate P-wave arrivals, the algorithm targets corresponding
S-wave arrivals using a similar procedure, assuming, however, that at most one S-wave arrival exists in
the time between two successive P-wave arrivals.
The STA/LTA in the S-wave case, s
S
, is dened using a 0.5 s short-term window and an expanding
long-term window:
s
S
i
=
var (n [i 0:5r :i]) + var (e [i 0:5r :i])
var (n [0 :i]) + var (e [0 :i])
: (A.9)
The short-term window is longer for S-waves than P-waves because S-waves tend to have longer duration.
The expanding window used to dene the long-term average serves to favor earlier arrivals over later when
multiple high-energy signals occur within the segment of data being processed.
The kurtosis for the S-wave case, k
S
, is dened as the average kurtosis in a 0.5 second window across
the horizontal channels:
k
S
i
=
1
2

4
(n[i 5r :i])
fvar (n [i 5r :i])g
2
+

4
(e[i 5r :i])
fvar (e [i 5r :i])g
2

(A.10)
The V/H ratio is inverted when dening the S-wave characteristic function, c
S
to give
c
S
i

s
S
i
k
S
i
v
i
; (A.11)
173
and the index,i

arg max
i
c
S
i
, coinciding with the maximum value of c
S
provides an initial measure-
ment of the arrival time. The AIC is then computed for a window centered on i

with length equal to
the time interval between the preceding P-wave arrival and i

. The index corresponding to the global
minimum of the AIC is registered as the S-wave arrival time if the SNR at that index is greater than two.
174

Abstract (if available)

Abstract I develop and apply automated processing procedures to derive seismicity catalogs from raw seismograms for two highly active fault zones in Southern California: (a) the San Jacinto Fault Zone (SJFZ) and (b) the region surrounding the 2019 Ridgecrest, California, earthquake sequence. For the Ridgecrest region, I also derive 3D models of seismic P- and S-wave velocity (Vₚ and Vₛ, respectively) structure, and their ratio (Vₚ/Vₛ). The SJFZ catalog comprises more than 2.5×10⁵ earthquakes in a 4×10³ km² region between 1 January 2008 and 1 January 2021. The Ridgecrest catalog comprises more than 9.5×10⁴ earthquakes in a four month period including the first three months of the aftershock sequence. Jointly interpreting spatiotemporal seismicity patterns, geologic data, and velocity models provides new insights into subsurface structure of and physical processes operating in the fault zones. To build the presented catalogs and velocity models, I (a) develop and implement a new algorithm to detect earthquakes and pick phase arrivals, (b) implement the Fast Marching Method to solve the eikonal equation in 3D media, (c) develop and implement two new algorithms for locating earthquakes in 3D media, and (d) modify a recently developed traveltime tomography method to adaptively resolve structural features based on data distribution. The entire framework developed is modular, enabling components to be replaced as methods advance, and can be applied to a wide range of network geometries. The presented methodological developments, data products, and interpretations will benefit various future studies of seismotectonics and crustal structure in Southern California and elsewhere.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

PDF

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

PDF

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

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

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

PDF

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

PDF

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

PDF

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

PDF

Stress-strain characterization of complex seismic sources

PDF

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

PDF

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

PDF

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

PDF

Structural clustering analysis of CVMS-4.26: a 3D seismic velocity model for southern California

PDF

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

PDF

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

PDF

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

PDF

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

PDF

Detailed properties of seismic waveforms and earthquake sources in Southern California

PDF

Crustal structure and seismotectonic patterns of the Northern Hikurangi margin, New Zealand, from a dense, short-period seismic deployment

Asset Metadata

Creator White, Malcolm Charles Adan (author)

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

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Geological Sciences

Degree Conferral Date 2021-08

Publication Date 07/31/2021

Defense Date 04/28/2021

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

Tag crustal structure,Earthquakes,faults,OAI-PMH Harvest,seismology,seismotectonics,tomography

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Ben-Zion, Yehuda (committee chair), Haldar, Justin (committee member), Sammis, Charles (committee member), Vidale, John (committee member)

Creator Email malcolm.white@usc.edu,malcolmw@mit.edu

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

Unique identifier UC15670637

Legacy Identifier etd-WhiteMalco-9953

Document Type Dissertation

Format application/pdf (imt)

Rights White, Malcolm Charles Adan

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 author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

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

Repository Email cisadmin@lib.usc.edu