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.