science.sciencemag.org/content/366/6463/346/suppl/DC1
Supplementary Materials for
Hierarchical interlocked orthogonal faulting in the 2019 Ridgecrest
earthquake sequence
Zachary E. Ross*, Benjamín Idini, Zhe Jia, Oliver L. Stephenson, Minyan Zhong, Xin
Wang, Zhongwen Zhan, Mark Simons, Eric J. Fielding, Sang
-Ho Yun, Egill Hauksson,
Angelyn W. Moore, Zhen Liu, Jungkyo Jung
*Corresponding author. Email: zross@gps.caltech.edu
Published 18 October 2019,
Science
366
, 346 (2019)
DOI:
10.1126/science.aaz0109
This PDF file includes:
Materials and Methods
Figs. S1 to S19
Tables S1 to S4
References
Submitted Manuscript: Confidential
2
Materials and Methods
Template matching and seismicity relocation
We used all of the continuous waveform data available for EH and HH channels within 80
km of the Ridgecrest mainshock from the Southern California Earthquake Data Center
(doi:10.7909/C3WD3xH1). The data were filtered between 2-15 Hz with a zero phase
Butterworth filter. All of the earthquakes listed in the Southern California Seismic Network
catalog from 2019/07/04 - 2019/07/25 within 60 km of the mainshock were sel
ected as template
events. The template events were relocated individually with a 3D velocity model (
33
) and then
pairwise relocated with HypoDD (
34
) using just the phase data from the SCSN. We applied a
template matching detection algorithm to the continuo
us waveform data starting on the first day
of the sequence. The procedure closely followed that of Ross et al. (
35
). This resulted in 111,918
detections. Then, we cross
-correlated each of these detections with the 100 nearest template
events to measure dif
ferential times. We used seismograms 1.5 s long, starting 0.25 s before the
arrival time. When catalog picks were unavailable we used 1D arrival predictions for selecting
the window start time. The maximum source-receiver distance was 80 km, and we only selected
event pairs that were within 5 km hypocentral distance. The detections were then relocated with
GrowClust (
36
), where we set a minimum correlation threshold of 0.75 and required a minimum
of 8 differential times per event pair for relocation. In tot
al 46,512 events were able to be
precisely relocated.
To estimate the relative errors on the relocated solutions, we performed a bootstrap analysis
by resampling the data 100 times and performing a new inversion for each resample (
36
). For the
horizontal a
nd vertical errors, the 90
th
percentile is 98 m and 340 m, respectively.
Multiple-subevent inversion for the M
w
6.4 and M
w
7.1 earthquakes
Here, we apply a subevent inversion method to estimate rupture process of the Mw 6.4 and
Mw 7.1 earthquake sequence. In the method, we treat large earthquake as a series of subevents of
varying location, timing and point source focal mechanism. This flexible yet simple source
parameterization allows us to constrain first
-order rupture complexity of large earthquakes
robustly. Mostly important for the complex Ridgecrest earthquakes, the subevent method does
not need to assume any particular fault geometry and rupture sequence, while still capturing
majority of the moments by inverting regional and teleseismic waveforms.
Subevent methods have been developed since the 1980s (
37
–
41
). The subevent event method
applied here based on a multiple point source method that has been successfully applied on deep
earthquakes (
40
,
42
). To resolve the rupture processes of the Ridgecrest
events, we incorporate a
variety of dataset, including teleseismic P and SH waves and near
-field full waveforms for
improved spatial and temporal resolutions. The ground phase velocity of the near field Rayleigh
and Love waves (~ 3 km/s) are much less than that of the teleseismic body waves (~ 10
-20
km/s), making it possible to distinguish rupture on closely located fault segments. We optimize
some of the source parameters (locations, centroid times and durations) nonlinearly, while
inverting for the rest of the source parameters (moment tensors) linearly. For the nonlinear
parameters, we generate Markov chains with a Metropolis
-Hasting algorithm, in which the
proposal models are generated by sampling through one of the nonlinear parameters while
keeping th
e other nonlinear parameters at their current values (
43
).
In each step of forward calculation, we linearly invert for the subevent moment tensors by
extending the approach
used by Minson and Dreger (
44
) to multiple subevents. We generate 72
Markov Chains with random first samples, and finally keep the 24 best fitting chains, to
Submitted Manuscript: Confidential
3
eliminate the dependency of the inversion on the initial values. Subevent moment tensors are
constrained to be deviatoric, with
no isotropic components. We apply a bounded uniform prior
probability density function for all non-linear parameters in the inversion. We also adopt a
penalty term to accommodate the moment
-duration scaling relationship observed for large
earthquakes (
45
) by rejecting models of extremely sharp or flat source time functions. We set the
data error to be 10% of the final misfit with the optimal models; the true data errors of
seismograms are trivial, but additional prediction errors may be introduced due to nonlinear
wave propagation effects (e.g. inaccurate velocity model).
We use teleseismic (epicentral distances from 30 to 90 degrees) P wave records of 33
stations, teleseismic SH wave records of 39 stations, regional (epicentral distances from 50 to
150 km) Rayleigh waves of 31 stations and regional Love waves of 30 stations. The data are
selected from all available Global Seismic Network broadband stations and the Southern
California Earthquake Data Center strong motion stations for good quality and azimuthal
coverage. We remove the instrument response and linear trends of the waveforms, and rotate the
two horizontal components to the radial and transverse components. We apply filter bands of
0.01-0.2 Hz for teleseismic body waves and 0.02-0.2 Hz for regiona
l surface waves, and allow
time shifts up to 2.0 s for P waves and 5.0 s for SH and surface waves to account for path
complexities and picking errors. The calculation of Green’s functions is based on the propagator
matrix method with plane wave approximati
on (
37
) for the teleseismic body waves, and the
frequency-wavenumber integration method
(
46
) for regional surface waves. The source side
velocity model is based on a combination of a 4
-layer 1D elastic model (Table S1; (
47
)) and
iasp91 model (
48
). We start
with one subevent and iteratively increase the number of subevents
until the main features of the waveform are well fitted and the sum of subevent moments agrees
with the long period moment.
The aftershock locations and the InSAR interferograms illuminate
the geometry of the faults,
which provide important constraints on our subevent models. For the M 6.4 event, we anchor the
location of the third subevent by the (arbitrarily chosen) maximum surface offset point on the
NE
-SW trending conjugate fault observed from the unwrapped interferograms, while the first
and second subevents are allowed to move freely in the whole space. A qualitative assessment of
the reliability of subevent locations is shown by the ensemble misfits and spanning areas in Fig.
S1. For the M 7.1 event, we fix the last subevent at the location where a slip patch is observed
from the static slip model, and search for the locations of other subevents along a NW
-SE
trending plane derived from the surface rupture and aftershocks. The depths of the M 7.1
subevents are searched for between 1-12 km, consistent with the distributions of static slip and
aftershocks with depth. The distribution of the Markov Chain ensemble is shown in Fig. S3.
Geodetic Observations
—GPS
We used geodetic observations
from ground stations of the Global Positioning System
(GPS) networks covering the area within 500 km of the earthquakes. Most of the stations are part
of the Network of the Americas operated by UNAVCO. We processed GPS data using the JPL
processing softwa
re GIPSY
-OASIS (https://gipsy-oasis.jpl.nasa.gov/) in Precise Point
Positioning mode (
49
) using JPL Final orbits and clocks. We estimated daily position time series
for all stations. Troposphere delays and gradients were estimated every 300 seconds, and si
ngle
station ambiguity resolution (
50
) was performed. Mean positions derived from daily positions of
15 days before Mw6.4 and 3 days after Mw7.1 earthquakes are then used to estimate total
Submitted Manuscript: Confidential
4
coseismic offsets and uncertainties of the earthquake sequence foll
owing the approach in Liu et
al (
51
).
Geodetic Observations
—Satellite Imaging
We used data from two synthetic aperture radar (SAR) satellite systems to extract
information about the earthquakes. We used SAR data acquired by the Japan Aerospace
Exploration Agency (JAXA) Advanced Land Observation Satellite
-2 (ALOS
-2) and by the
Copernicus Sentinel
-1A and -1B satellites operated by the European Space Agency (ESA). We
performed three types of analysis of the SAR data to extract information about the static gro
und
displacement during the earthquakes and to map the locations where the fault ruptures reached
the surface, including interferometric SAR (InSAR) analysis, interferometric coherence and
coherence change used to calculate a Damage Proxy Map (DPM), and SA
R pixel matching or
pixel offset tracking (described below). We calculated interferograms from the ALOS
-2 data and
interferograms, DPM, and SAR pixel offsets from the Sentinel
-1 data.
The JAXA ALOS
-2 SAR uses an L
-band (24 cm wavelength) radar and ALOS
-2 s
cenes we
used scenes from ascending path 65 that were acquired before (2018-04-16) and after (2019-07-
08) the earthquake sequence (Table S2). The pre
-event scenes were acquired in stripmap mode
(SM3), and the post
-event scene was acquired in ScanSAR mode (
Figure 3). The stripmap-
ScanSAR processing was done using additional modules (
52
) of the InSAR Scientific
Computing Environment (ISCE) (
53
). The interferogram was processed with the Shuttle Radar
Topography Mission (SRTM) version 3 digital elevation model
at 1
-arcsecond spacing and had 1
sample in range and 2 lines in azimuth averaged (1 by 2 looks) for the phase unwrapping with
SNAPHU (
54
).
The Copernicus Sentinel
-1 SAR uses an S
-band (5.6 cm wavelength) radar and two tracks
cover the Ridgecrest area, ascending track 64 and descending track 71. Scenes we used are listed
in Table S2. Sentinel
-1 SAR data was acquired in Terrain Observation by Progressive Scans
(TOPS) mode (250 km swath) and processed with the ISCE software, using the SRTM elevation
data. Inte
rferograms were processed with averaging of 3 range samples and 1 azimuth line to do
the phase unwrapping with SNAPHU.
The Sentinel
-1 coseismic interferograms extend 5 days and 11 days after the M7.1 earthquake,
with the first post
-quake acquisition on 10 July for track 64 and 16 July for track 71. This means
that the interferograms include a small amount of postseismic deformation. Preliminary analysis
of the Sentinel
-1 data acquired in the month and a half after the first acquisitions shows that the
posts
eismic deformation rate is moderate, with a maximum rate of about 1 mm/day in the radar
line-of-sight. Unless the postseismic deformation rate decreased rapidly in the first few days,
then the contribution of postseismic deformation to the coseismic measur
ements is likely 1 cm or
less.
Offset Field Estimation
We estimate offset fields from coregistered Sentinel
-1 (S1) SAR amplitude images using a
feature tracking algorithm (
55
), which performs dense cross
-correlation matching operation to
estimate both ran
ge and azimuth offsets. We perform feature tracking on one 6-day pair of
ascending track 64 (A64 pair) and one 12-
day pair of descending track 71 (D71 pair), and apply
a median filter that corresponds to a size about 1km x 1km on the obtained offset fields
. The
spatial resolution of both offset fields is about 200 m x 200 m. Compared with interferograms
which also measure displacement in the radar line-of-sight (LOS) direction, range offset fields
Submitted Manuscript: Confidential
5
are of lower spatial resolution but are able to give measure
ments at near
-fault area where
interferograms decorrelate. Besides LOS offsets, offset fields include azimuth offsets measured
in the satellite along
-track direction, which provide additional constraint to the fault model.
Azimuth offset fields is noisier than range offset fields, because azimuth pixel size (14 m) is
much larger than range pixel size (2.3 m).
Static Inversion
Data sources:
We used GPS data, ALOS
-2 ascending, and Sentinel
-1 descending coseismic
interferograms
together with Sentinel
-1 azim
uth offsets for the static inversion (described
above).
Satellite-data downsampling and error models:
For computational efficacy, each interferogram
is downsampled following a data
-resolution approach which attempts to provide locally averaged
data at a scale commensurate with the ability of the model to produce variable predictions at that
location (
56
). A data covariance matrix is constructed from the standard deviation of the data
contained in the blocks emerging from the decimation and a correlation function. The correlation
function is derived from the region in the images where the deformation caused by the
earthquake is negligible. We compute an empirical covariogram (
56
) over this region from which
we extract the correlation function as a function
of distance and use this to estimate the full data
covariance matrix. Errors in the interferogram phase are largely due to atmospheric water vapor
variations that have a strong spatial correlation, so this procedure captures that. Errors in the
pixel offs
ets are due to the local cross
-correlation quality and have little effects from atmospheric
water vapor.
GPS:
We combine the GPS measurements from the M
w
6.4 and M
w
7.1 earthquakes to obtain a
measurement consistent with the deformation recovered from SA
R. The data covariance is set to
be diagonal with an independent variance for each station.
Fault geometry:
We manually infer a fault trace from the damage proxy map and the wrapped
interferograms. The resulting trace is consistent the complex deformation observed at surface
(Fig. 5), but contrasts with the localization of deep seismicity (> 5 km) in an almost s
traight
plane (Fig. 3). As a result, we set a straight fault plane below 5 km depth that follows the deep
seismicity and connects to the inferred fault trace using appropriate dipping angles (Fig. 5).
We selected a depth
-varying patch size consistent with
the loss of resolution that geodetic data
provides as with increasing depth. Our preferred fault geometry has two segments with
subparallel strands. Since proximal parallel strands will have strong trade-offs, we only project
one of the subparallel strands throughout the seismogenic zone, limiting the other to the
shallowest few kilometers. It is not possible to discriminate which of these two strands
accommodate slip at depth but we assume here that the strand with the maximum surface offset
corresponds t
o the dominant fault at depth.
The elastic model and prediction errors:
We use the 4
-layer 1D elastic model shown in Table S1
for Southern California (
47
). Our imperfect knowledge of the elastic structure may lead to
overfitting given resulting errors in our model predictions. We account for prediction errors
Submitted Manuscript: Confidential
6
through a prediction covariance matrix, C
p
, using an approach based on perturbation analysis
(
57
). C
p
is added to the data covariance matrix, Cd, to form a total misfit covariance, C
x
= C
p
+
C
d
, which is used in the inverse problem.
C
p
depends on an assumed covariance for the elastic
structure, C
u
, here assumed to be diagonal and to correspond to a 10% uncertainty on log(mu),
where mu is the shear modulus. (
58
,
59
).
Bayesian sampling:
Our results in Fig. 4 represent the mean value of an ensemble of models
obtained from sampling the posterior PDF using a modified version of CATMIP, an MCMC-like
algorithm
(
60
). This procedure is free of any a priori smoothing, a common practice in
geophysical inversions that may lead to inconsistent results depending on the arbitrary selection
of the smoothing parameters (
60
). The algorithm is embedded into AlTar, a fully
parallelized
software suite designed for sampling of large inverse problems (
58
,
59
). We select
the priors based on the known tectonics of the region. Strike
-slip priors are set as bounded
uniform distributions between -2 and 12 m oriented in a direction
consistent with plate motion
(e.g., dextral for the M
w
7.1 event and sinistral for the M
w
6.4 event). Dip slip priors use a
Gaussian prior with 0.5 m as standard deviation and centered at zero.
Garlock creep analysis
Individual interferograms are often cont
aminated by atmospheric noise, making it hard to
extract small signals from them accurately. In order to resolve creep on the Garlock fault, we
take all available data for Sentinel 1 (a C
-band SAR satellite operated by the European Space
Agency) ascending track 64 between June 1
st
and August 27
th
and use this to form a time series
of ground displacement. We then fit a step function at the time of the earthquake (
61
). The fitted
step function will contain coseismic deformation, as well as any motion that has
occurred after
the earthquake but on timescales much smaller than the time series. The step function fit
mitigates atmospheric noise and allows much smaller offset signals to be observed.
This step function is plotted in Figure 6, along with the gradient shading of the step function,
allowing us to resolve areas of high gradient in the deformation field. The gradient shading
shows a zone of high deformation gradient to the south of the Mw 7.1 mainshock along a section
of the Garlock fault approximately 30 km long, with several zones of very sharp deformation
gradients. As there was minimal seismicity near the Garlock fault, we infer this deformation to
be due to creep on the Garlock fault.
To explore the magnitude of the creep we take profiles across the Garlock where the
deformation gradient is sharpest. In order to separate out the motion on the Garlock fault from
the long wavelength deformation field due to the earthquake we take profiles across the Garlock
fault then remove a polynomial from the profiles
. The residuals are plotted in Figure 6.
The overall width of the deformation profile is controlled by the deepest extent of the creep;
examining the profiles indicates that creep is not likely to extend below the upper few hundred
meters. The gradient of
the deformation profile on the fault is determined by the upper extent of
the creep, with increasing deformation gradient as creep comes nearer to the surface. Profiles
suggest observable surface creep in a few small sections, with much larger zones where creep
has not made it to the surface. Profile AA’ shows a deformation profile consistent with near
surface slip that hasn’t broken the surface. A maximum line of sight offset of 20 mm at or very
near the surface is observed in profile BB’. Profile CC’ shows what appears to be two strands of
fault offset, with the largest offset off the main strand of the Garlock.
Submitted Manuscript: Confidential
7
Burst boundaries in Sentinel
-1 interferograms can create apparent discontinuities in the
calculated deformation field. One such burst boundary, s
ubparallel to the Garlock fault. is
labelled in Figure 6. Care should be taken to distinguish this from an actual tectonic signal.
Submitted Manuscript: Confidential
10
Thickness (km)
Vs (km/s)
Vp (km/s)
Density (g/cm3)
2.5
2.6
4.34
2.75
3.0
3.5
5.88
2.77
24.5
3.6
6.30
2.90
0.0
4.5
7.74
3.10
Table S1
.
The elastic model for Southern California (
47
) used in the computation of static Green ́s
Functions and prediction errors.
Submitted Manuscript: Confidential
11
Sensor
Track number
(Ascending/Descending)
Preseismic
Acquisition
Date
Postseismic
acquisition date
ALOS
-
2*
A65
20180416***
20190708
Sentinel
-
1**
A64
20190704****
20190710
Sentinel
-
1
D71
20190704
20190716
Table S2
*Advanced Land Orbiting Satellite 2, L-band SAR satellite operated by the Japanese Aerospace
Exploration Agency
** C-band SAR satellite operated by the European Space Agency
***ALOS A65 has not been acquired regularly, meaning it was necessary to use a scene taken
more than a year earlier
**** Currently Sentinel 1 A64 is acquired every six days (using Sentinel A and B of the
constellation) and D71 every 12 days
Submitted Manuscript: Confidential
12
Centroid
time (s)
Duration
(s)
Longitude
(°)
Latitude
(°)
Depth
(km)
Mrr (10
26
dyne
-cm)
Mtt (10
26
dyne
-cm)
Mpp (10
26
dyne
-cm)
Mrt (10
26
dyne
-cm)
Mrp (10
26
dyne
-cm)
Mtp
(10
26
dyne
-cm)
E1
2.68
6.42
-
117.512
35.724
9.76
0.0025
-
0.1831
0.1806
0.0256
0.0344
0.0198
E2
6.24
6.40
-
117.510
35.682
5.88
0.0117
-
0.2265
0.2147
0.0110
-
0.0052
-
0.0061
E3
9.05
6.08
-
117.552
35.633
4.48
0.0047
-
0.2331
0.2284
0.0163
-
0.0267
0.0419
Table S3.
Subevent model parameters for the M 6.4 earthquake.
Submitted Manuscript: Confidential
13
Centroid
time (s)
Duration
(s)
Longitude
(°)
Latitude
(°)
Depth
(km)
Mrr (10
26
dyne
-cm)
Mtt (10
26
dyne
-cm)
Mpp (10
26
dyne
-cm)
Mrt (10
26
dyne
-cm)
Mrp (10
26
dyne
-cm)
Mtp (10
26
dyne
-cm)
E1
7.00
9.26
-
117.627
35.794
8.20
-
0.0316
-
2.6543
2.6860
0.1987
0.0525
0.7137
E2
10.28
6.73
-
117.572
35.741
4.37
0.0280
-
0.8067
0.7787
0.2993
0.2412
0.1385
E3
14.53
4.65
-
117.540
35.709
5.86
0.1354
-
0.4079
0.2724
0.0307
0.2180
0.0689
E
4
19.88
5.70
-
117.444
35.616
4.00
0.0388
-
0.5098
0.4710
0.0173
0.1461
0.0861
Table S4.
Subevent model parameters for the M 7.1 earthquake.
Submitted Manuscript: Confidential
14
Fig. S1.
Subevent locations of the M 6.4 eve
nt revealed by the ensemble of Markov Chain samples.
Circles in cyan show the optimal subevent locations corresponding to Fig. 4. The dot clouds
indicate Markov Chain samples of different horizontal locations for three subevents. Their color
show the corre
sponding data misfit. Faults are indicated by the black lines. Stars indicate the
hypocenters of the M 6.4 and M 7.1 events.
Submitted Manuscript: Confidential
15
Fig. S2.
Waveform fits for the preferred subevent model of the M 6.4 earthquake. Data and synthetic
waveforms are shown in black and red, respectively. The numbers below ea
ch trace are the
azimuth and distance in degrees. (A) P waves in velocity. (B) P waves in displacement. (C) SH
waves in displacement. (D) Rayleigh waves in velocity. Traces are aligned by the event origin
time. (E) Love waves in velocity. Traces are aligned by the event origin time.
Submitted Manuscript: Confidential
16
Fig. S3.
Subevent locations of the M 7.1 event. Circles in cyan show the optimal subevent locations
corresponding to Fig. 4. (A) A NW
-SE trending vertical plane (N40°E) where the subevent
locations are searched. The plane is indicated by the dashed blue line. Backg
round gray lines
show the faults. Stars indicate the hypocenters of the M 7.1 and M 6.4 events. (B) The ensembles
of Markov Chain samples for 4 subevents of the M 7.1 mainshock. Samples (dot clouds) are
plotted on the cross
-section of the NW-SE trending pl
ane in (A). Their color show the
corresponding data misfit.
Submitted Manuscript: Confidential
17
Fig. S4.
Same as Fig. S2 but for the M 7.1 event.
Submitted Manuscript: Confidential
18
Fig. S5.
Seismicity cross
-section along NW-SE fault plane. Earthquakes that occurred prior to the Mw
7.1 mainshock are colored black, while events afterward are colored red. M>5 events are
indicated by blue stars.
Submitted Manuscript: Confidential
19
Fig. S6.
Seismicity cross
-section along SW-NE fault pl
ane. Earthquakes that occurred prior to the Mw
7.1 mainshock are colored black, while events afterward are colored red. Mw 6.4 event is
indicated by blue star. A large gap in seismicity is present with dimensions roughly 6x6 km (blue
dashed region). Note
the vertical lineations indicating orthogonal faulting throughout the profile.
Submitted Manuscript: Confidential
20
Fig. S7.
Evolution of seismicity projected a
long strike of M
w
7.1 mainshock. The M
w
6.4 event ruptured a
~6 km segment of the NW fault (Fig. 5) and left behind a gap of aftershocks. The NW edge of
the seismic activity expanded over 34 hours through the occurrence of a series of M > 4
earthquakes, le
ading to the nucleation of the M
w
7.1 earthquake.
Submitted Manuscript: Confidential
21
Fig. S8.
Offset fields for Sentinel
-1 data.
Submitted Manuscript: Confidential
22
Fig. S9.
Summary of information for the Garlock swarm. Upper panel shows a map view of seismicity.
Brown lines indicate faults. Lower panel shows time history of magnitude and event counts.
Submitted Manuscript: Confidential
23
Fig. S10.
The top two panels show log of the
seismicity density between the Mw 6.4 and Mw 7.1 (left,
blue color scale) and after the Mw 7.1 (right, red color scale) in map view. The lines show a
simplified surface fault geometry of the Mw 6.4 and 7.1 earthquakes combined. The main
rupture during the
Mw 7.1 is represented by the thick line between A and A’. The bottom two
panels show log of the seismicity density from the top two panels projected onto the line between
A and A’ extended vertically down to 15 km, only considering seismicity within 5 km of the
surface. Black dashed regions enclose areas of low event density within each seismicity
distribution. Again left is between the Mw 6.4 and the Mw 7.1 and right is post Mw 7.1.