Supplementary Materials
Geodetic observations and three-dimensional co-seismic
displacements
Both Sentinel-1 and ALOS-2 (advanced land observation satellite-2) SAR
data (see Table S1 for details) are used to obtain the co-seismic
displacement observations of Maduo earthquake, where the DInSAR
(differential interferometric SAR)[Gabriel et al. , 1989], POT
(pixel offset tracking)[Michel et al. , 1999], MAI (multiple
aperture interferometry)[Bechor and Zebker , 2006], and/or BOI
(burst overlap interferometry)[Grandin et al. , 2016] methods
are employed to process these SAR images based on the GAMMA software.
The SRTM (shuttle radar topographic mission) 1-arc-second DEM (digital
elevation model) data is used here for correcting the topographic
component and for assisting in co-registration process. We apply the
multi-look operation of 30×8, 8×20, and 5×32 (range×azimuth) for
Sentinel-1, ascending strip-map-mode ALOS-2, and descending ScanSAR-mode
ALOS-2 SAR images, respectively, and the final resolution of
displacement observations is about 100m×100m. Particularly for the
descending ScanSAR-mode ALOS-2 data, there are five separate beams, and
each beam is separately processed based on the DInSAR, MAI and POT
methods, then mosaiced together in the geographical coordinate system
based on the overlaps between adjacent beams.
The DInSAR interferograms are filtered by the adaptive filter based on
local fringe spectrum[Goldstein and Werner , 1998], and
unwrapped by the minimum cost flow method[C W Chen and Zebker ,
2002]. For the L-band ALOS-2 data, due to its longer wavelength, the
corresponding interferograms are contaminated with the ionospheric
delays. For the ascending strip-map-mode ALOS-2 SAR images with small
coverage, we correct the contained ionospheric delay by first degree
polynomials. For the descending ScanSAR-mode ALOS-2 SAR images, we apply
a range split spectrum method[Gomba et al. , 2016] to mitigate
the ionospheric delay[J Liu et al. , 2021]. Different from the
DInSAR method that can only obtain displacements along the LOS
(line-of-sight) direction, the POT method can obtain displacements along
both the LOS and azimuth directions. Before the POT process of
Sentinel-1 data, the SAR images should be deramped[Wegmüller et
al. , 2016]. To obtain the pixel offset, we employ a matching window
size of 128×128 pixels with the oversampling factor of two for all SAR
images, and use a second-order polynomial to fit the possible ramp based
on the signals in the far-field regions[J Hu et al. , 2012].
Due to the significant ionospheric disturbs, the POT azimuth
displacements from the ALOS-2 images can hardly recognize any valuable
displacement signals. The signal-to-noise ratio (SNR) of POT
displacements is generally proportional to the spatial resolution of SAR
images, and may also be affected by the spatial baseline of SAR images
from the comparison between the ascending/descending Sentinel-1 POT
azimuth displacements. Similar to the POT azimuth observations of ALOS-2
data, the MAI observations are also seriously contaminated by the
strip-shape ionospheric disturbs. Due to the limited azimuth doppler
bandwidth and spatial resolution, the MAI method is not used for the
Sentinel-1 data. For the BOI method, it is only applicable for the
Sentinel-1 data in the burst overlap area, whose size is about
86km×1.5km (range×azimuth) accounting for about 7.5% of a burst along
the azimuth direction. For a target burst overlap area, four SLC (single
look complex) images can be obtained from the two successive bursts of
the primary and secondary SAR images, and the BOI azimuth displacement
is derived by double differential operation of these four SLC images.
Therefore, sixteen displacement observations can be obtained, among
which the azimuth displacement observations from the ALOS-2 data contain
serious ionospheric delays and the descending Sentinel-1 POT azimuth
observation contains serious decorrelations.
To calculate the three-dimensional (3D) co-seismic displacements, a
InSAR-based Strain Model and Variance Component Estimation method
(SM-VCE)[J Liu et al. , 2018; J Liu et al. , 2019a] is
used to combine the remaining 11 displacement observations. The strain
model represents the geophysical relationship between adjacent points’
3D displacements[Guglielmino et al. , 2011], and can be used
to establish the observing function between the 3D displacements at a
target point and its surrounding observations within a window. Compared
with the traditional methods that calculate the 3D displacements of a
target point based only on observations on this point[Jung et
al. , 2011; Wright et al. , 2004], the SM-VCE method can
incorporate context observations in the calculation. Besides, the
relative weight of different observations in the SM-VCE method is
determined by the well-known variance component estimation method in a
posteriori way[S W Hu and Xiao , 2016], which is more accurate
than the weight determined by some a priori information about the
observation noise.
Since the observations within a window are used to establish the
observing function, in the final 3D displacement field, it is possible
for the SM-VCE method to fill in some gaps in the observations. Some
band-shaped signals are observed in the north-south displacement field,
which can be attributed to the fact that the BOI observations are only
available in the burst overlap regions and more sensitive to the
north-south component than the east-west. In addition, one key point for
SM-VCE method is the window size when establishing the observing
function based on the strain model. In the window, the displacement
gradient is assumed to be constant[Guglielmino et al. , 2011].
If the window size is too large, this assumption would be violated,
while if too small, there are insufficient data points to suppress the
high-frequency noise (e.g., the random noise in the Sentinel-1 POT
azimuth observations). In this paper, the window size of 31×31 pixels is
empirically determined after comparing with the result of other window
sizes. Besides, when calculating 3D displacements for a target point
near the fault rupture region, heterogeneous observations at both sides
of the fault would be included in the window, which is undesirable and
would deteriorate the 3D displacements. In this paper, the fault traces
are manually mapped from the ascending Sentinel-1 POT LOS observation,
and then used to discard the observations within the window across the
fault.
Finite fault inversion
We derive the FFM with a simulated anealing inversion method [Ji
et al. , 2002], which allows the joint inversion of seismic waveforms
and static deformation data. Since seismic and static data provide
complementary constraints on the kinematic rupture process, joint
inversion can greatly suppress trade-off among model parameters.
The seismic waveform dataset includes teleseismic P and SH waves,
regional broadband and nearby high-rate GPS waveform records. We
download teleseismic waveform data from IRIS (http://www.iris.edu)
recorded by the Global Seismological Network (GSN) that has the most
uniform spatial coverage to provide best teleseismic (30° - 90°) dataset
for FFM inversion. We selected 42 P-waves and 39 SH-waves based on the
data quality and azimuthal coverage. Instrument responses are removed
from the raw waveform data and converted to displacement after a
low-pass filtering with corner frequency of 1.0 Hz. Although 1Hz is
relatively high frequency for deterministic inversion, the waveforms are
dominated by relatively low frequency energy. We also use
three-component 1Hz high-rate GPS waveform observaations recorded by 10
nearest stations (Fig.S4), see Gao et al. [2021] for the
detailed GPS data processing and data availability. The sampling rate of
the GPS data is 1 Hz, so in the inversion we filter these data to 0.02
to 0.4 Hz. Besides, three-component displacement waveform records of 7
regional broadband stations are also included in the inversion (top
traces in Fig.S3). These broadband stations are part of the national
broadband seismic network. Both the regional and teleseismic Green’s
functions are calcualted with a Tibet crustal velocity model proposed byL Zhu and Helmberger [1996] using a frequency-wavenumber
integration method[L P Zhu and Rivera , 2002].
We use 6 down-sampled InSAR and SAR datasets in the joint inversion. We
did not inlcude all the datasets in Fig.S1 in the inversion, only the
images that have high data quality are used. For InSAR data, we solve
for quadratic ramps in the inversion to correct for orbital errors. The
static green’s functions at free surface are calculated by using the
same 1D velocity model[L Zhu and Helmberger , 1996] as for the
dynamic Green’s function.
To parameterize the finite fault model, we divide the rectangle fault
plane into smaller sub-faults and solve for the slip amplitude and
direction, risetime and rupture velocity on each sub-fault. For each
parameter, we specify the bounds and a discretization interval. The
bounds and increments are selected based on trial-and-error. We define
the misfit function as:
Ewf + WI×EI +
WS×S + WM×M
where Ewf is the waveform misfit. EI is
the geodetic misfit, S is a normalized, second derivative of slip
between adjacent patches (smoothing), M is a normalized seismic moment,
and WI, WS and WM are
the relative weighting for the geodic misfit, smoothing, and seismic
moment, respectively. In the joint inversion, the relative weighting is
realized by normalizing the misfit of each dataset by the misfit derived
by inverting the individual dataset. Based on our previous experiences
(e.g. [Wei et al. , 2011]) we set the weight for the geodetic
misfits double that of the waveform misfits. A simulated annealing
algorithm[Ji et al. , 2002] is adopted to find the best
fitting model parameters in the joint.
We set the sub-fault size to be 3 km × 2 km, constrain the rake angles
to be between -45° to 45°, and allow rupture velocity to vary between
3.0- 2.0 km/s, as guided by the back-projection results. We also assume
the risetime to be an arc of the cosine function and to vary between 1.0
and 9.0 s with 0.8 s steps. The slip amplitude can change from 0 to 6 m.
Note that the selection of rupture speed range is narrowed down by the
BP results, which can reduce the trade-offs between rupture speed and
risetime.
We use 10 fault segments (Table S2) to mimic the primary variations of
strike as revealed in the 3D deformation (Fig.S2). The dip direction and
dip angle of these fault segments are determined by fitting the geodetic
observations and the aftershock distribution. We then conduct joint
inversion of waveform data and the geodetic data. Our preferred finite
fault model is presented in Fig.2a and Fig.3, with the waveform fittings
presented in Fig.S3-5 and the geodetic data fitting shown in Fig.S6-8.
As shown in these figures, both datasets are very well fitted. To
further verify the quality of waveform fitting, we also convert the
horizontal components of the high-rate GPS waveform to velocity and
filter them to 0.02 to 0.4 Hz, which still show good agreements between
data and the synthetics (Fig.S4a). At broadband, both the dynamic
waveforms and the static offsets in these GPS data are also well-fitted
(Fig.S4b).
Multiple Point Source Inversion
To constrain the first order complexity of the rupture process of the
Maduo mainshock, we conduct a multiple point source (MPS)
inversion[Q Shi et al. , 2018] using nearby high-rate GPS
waveforms, regional broadband waveforms and teleseismic body waves. The
same high-rate GPS data used in the FFM inversion are adopted here. We
download the teleseismic broadband data from the Incorporated Research
Institutions for Seismology (IRIS) data centre and regional broadband
data (within the epicentre distance of 800 km) from China Earthquake
Networks Centre (CENC). We remove the instrument response for the
downloaded data and select 61 regional broadband stations and 49
teleseismic stations with relatively even azimuthal coverage based on
their coordinates (Fig. S9) and data quality (e.g., signal to noise
ratio, whether the waveform is clipped).
We conduct MPS inversion using the Markov Chain Monte Carlo (MCMC)
sampling algorithm that is proposed by [Q Shi et al. , 2018].
Considering the different site condition and 3D velocity structure, we
apply station-dependent filtering to the regional and teleseismic
waveforms. At the same time, we apply broader band for the high-rate GPS
waveforms. The determination of frequency ranges is initially based on
the similar degree of complexities of the observations and synthetics.
We then conduct the MPS inversion in an iterative fashion and further
refine the frequency ranges with updated MPS solutions. We start the
inversion using three point-sources, and then gradually increase to
seven sources. Here, we show the optimal waveform fitting is achieved
using six point-sources.
The rupture process of the Mw7.4 earthquake is well represented by six
point-sources (Figure 2 and Figure S10-12). The statistics of waveform
cross-correlation coefficients are shown in Figure S13. The first
subevent (Mw7.07) is located near the epicentre of the earthquake, only
about 7 km to the southeast. The following subevents show a bilateral
rupture processes, with two subevents (Mw6.95 and Mw6.55) located to the
west and the other three (Mw6.93, Mw6.84 and Mw6.83) located to the east
of the first subevent. We add together the moment-weighted moment tensor
of the six subevents to calculate the cumulative moment tensor, which
shows a relative small compensated linear vector dipole (CLVD)
component. It should be noted that the CLVD percentage in the cumulative
moment tensor is mainly contributed by the oblique north-dipping
subevent that is located near the west end of the fault.
Back-projection
We back-project teleseismic high frequency waveform data recorded by the
Australian (AU) and European (EU) arrays (Fig. S14) to image the rupture
process of the Maduo earthquake. Traditionally, beginning portion of the
array P-waves is align to calibrate the 3D velocity structure and this
calibration is applied to the rest of the rupture area. However, due to
3D velocity structure at the source side, the calibration obtained for
the beginning portion may not be sufficient for later part of the
rupture that are located further away from the epicenter. To handle this
issue and improve the BP resolution, we apply a new travel-time
calibration strategy [Zeng et al. , 2022] to the MUltiple
SIgnal Classification BP method[Meng et al. , 2011; Zeng
et al. , 2020] to better resolve the high-frequency (HF) radiation of
the 2021 Mw 7.4 Maduo earthquake. To understand the performance of the
new calibration and derive an uncertainty estimation, we use the
conventional calibration and the new calibration strategy to relocate
M>5 earthquakes, mostly aftershocks, that are located near
the surface rupture. The BP locations of these events are then compared
with the relocation results from [W Wang et al. , 2021] (Fig.
S15). We first note that in both conventional and calibrated BP
relocation results, events located to the west/east of the mainshock
hypocenter have large mis-relocation in AU/EU arrays, respectively.
Considering that in the calibrated BP we assume velocity heterogeneities
bias travel-time error linearly, this mis-relocation pattern change
around the mainshock epicenter may indicate some structural
heterogeneities here. Compared with relocation from conventional BP,
relocation from the calibrated BP using EU data can better match the
catalog location in west of the epicenter. Relocation from the
calibrated BP using AU data has similar resolution as conventional BP.
Considering the different mis-location pattern we observed, therefore,
to reach higher resolution, we use AU and EU arrays to image the east
and west propagation ruptures, respectively. And the average
mis-location results in a lower bound uncertainty estimation of 7.5 and
4.5 km for calibrated BP using AU and EU array data.