Resolving co- and early post-seismic slip variations of the 2021 M_W 7.4 Madoi earthquake in east Bayan Har block with a block-wide distributed deformation mode from satellite synthetic aperture radar data

1.   Introduction

On 21 May 2021 at 18:04:13 (UTC time), an earthquake with a moment magnitude (M_W) of 7.4 (hereafter Madoi earthquake) jolted the northeast Bayan Har block at 98.46°E, 34.65°N in the Tibetan Plateau (Figure 1). As one of the most active tectonic regions in the Tibetan Plateau, the Bayan Har block is well known for its intensive tectonic deformation (e.g., Wang M and Shen Z-K, 2020) and frequent large earthquakes (e.g., Lasserre et al., 2005; Wang Q et al., 2011; Elliott et al., 2010; Sun JB et al., 2018; Li X et al., 2020). The 2021 Madoi earthquake was the most powerful one in the Bayan Har block and its surroundings since the 2008 M_W 7.8 Wenchuan earthquake, leading to strong regional seismic intensity of up to X (scaled in modified Mercalli intensity); the ground surface was heavily fractured and buildings suffered extensive damage (Pan JW et al., 2021; Li ZM et al., 2021). In addition, the Madoi earthquake was followed by more than 1000 aftershocks within 6 days, including 252 M > 5 events (Wang WL et al., 2021; Figure 1c).

Figure  1.  Tectonic setting of the study area. (a) Focal mechanisms of M_W ≥ 6.8 earthquakes recorded by GCMT in/around the Tibetan Plateau. The Bayan Har block is highlighted by a light red filled polygon. Thick cyan arrows characterize the clockwise rotation of the crust. (b) Thin black lines are regional faults; thick blue lines indicate the east Kunlun fault (EKLF) and the Jiangcuo fault (JCF). Green and yellow rectangles outline the spatial coverage of the Sentinel-1 data used in this study. Red dots are the aftershock relocations during 2021/05/21 and 2021/05/29 from Wang WL et al., (2021). Black arrows are the GNSS-measured coseismic horizontal deformation (Li ZC et al., 2021); purple arrows are deformations predicted by the coseismic slip model (Figure 3a) inverted in this study.

DownLoad: Full-Size Img PowerPoint

Clockwise rotation of the crust is one of the most prominent kinematic features of the east Tibetan Plateau. The Bayan Har block, however, is located right at a transition zone where the crustal movement switched from the north-south direction to the north-east/east-west direction (Cyan arrows in Figure 1a), accommodating a significant portion of the eastward extrusion of the plateau. As a response to this shear strain buildup, all of the M_W ≥ 6.8 earthquakes recorded instrumentally since 1976 in the interior of the Tibetan Plateau have been associated with the Bayan Har block (Figure 1a). Among them are the 2001 M_W 7.8 Hoh Xil earthquake which occurred in the block’s north boundary, and the 2010 M_W 6.9 Yushu earthquake in its south boundary. The Madoi earthquake is located about 70 km south of the east Kunlun fault (Figure 1b) which has been of concern for dozens of years because of its high seismic hazard (e.g., Van der Woerd et al., 2002; He WG et al., 2006; Kirby et al., 2007; Diao FQ et al., 2019; Wang M and Shen Z-K, 2020; Zhu LY et al., 2021), and it fills the seismic gap in the northeast Bayan Har block where no M_W ≥ 6.8 earthquake had been recorded instrumentally by the Global Centroid Mo- ment Tensor (GCMT) catalog (Figure 1a). An earthquake of magnitude similar to that of the Madoi earthquake is unusual and thus provides us a great opportunity to gain a better understanding of the kinematics of large earthquakes, and how the extrusion deformation across the east Bayan Har block is accommodated. Beyond that, the Madoi earthquake also enables us to assess the likelihood of future earthquakes on surrounding faults, from the spatiotemporal evolution of its stress transfer (e.g, Li SS et al., 2020).

Coseimic surface deformation of the Madoi earthquake was preliminarily measured by GNSS (Global Navigation Satellite System) (Li ZC et al., 2021), but data from GNSS stations in the epicentral area were too limited to reveal detailed coseismic ground motion (Figure 1b) or to invert for elaborate fault rupture description. Unlike previous M_W ≥ 6.8 earthquakes occurring in the Bayan Har block, the Madoi earthquake is the first large event from which the entire coseismic surface deformation can be precisely and timely mapped, from both ascending and descending satellite synthetic aperture radar (SAR) data with a short temporal baseline (e.g, 6 days). Here we use the spaceborne Sentinel-1 data to characterize the co- and early post-seismic surface deformation and subsurface fault slip distributions of the Madoi earthquake. Then the inferred co- and post-seismic slip distributions are used to shed light on fault slip heterogeneities and to explore spatiotemporal evolution of the early afterslip. We further discuss how shear strain is accommodated across the east Bayan Har block based on a comprehensive analysis of the GNSS velocity profiles, strain rates, and magnetotelluric data. Finally, we evaluate the co- and post-seismic stress loadings on regional faults and identify fault segments of high seismic hazard.

2.   Data and Inversion Configuration

2.1   Data

In this study, to characterize the kinematic responses of the Earth’s surface to the 2021 Madoi earthquake, we apply a combination of interferometry and SAR pixel offset tracking methods to Sentinel-1 SAR data from the satellite’s ascending (path 99) and descending (path 106) tracks (Table S1). Detailed InSAR and pixel offset tracking processing steps are described in the supporting information (Text S1). In addition, we process three pairs of Sentinel-1 data (Table S1) to capture the early 6-day and 24-day postseismic surface deformation associated with the earthquake, following the same InSAR procedure as used in the coseismic case. The mapped co- and post-seismic deformation fields allow us to reveal rupture details of this earthquake as well as the spatiotemporal evolution of the early afterslip, given the limited ground-based geodetic measurements (e.g, GNSS) and the sparsely distributed seismic network near the epicenter due to the harsh environment. The coseismic horizontal displacements measured at six GNSS stations near the epicentral area (Li ZC et al., 2021) are not included in our inversions (Figure 1b); instead, we use them as an independent dataset to examine the reliability of the coseismic slip model that was inverted from the SAR data.

2.2   Inversion Configuration

2.2.1   Fault model

The pixel offsets along the range directions provide valuable constraints on the near field deformation and help reveal details of how the fault ruptured the Earth’s surface, especially in regions where InSAR observations decorrelate due to possible surface ruptures of an earthquake. Therefore, we use the range offset maps to determine the location (i.e., fault traces) of the earthquake fault. By far, the dip direction and dip angles are the remaining two parameters of the seismogenic fault that we needed to solve. In the following inversions, based on the coseismic InSAR and range offset observations, we test two fault model scenarios — dipping to northeast and southwest respectively — as suggested by focal mechanisms from GCMT and the United States Geological Survey (USGS). For each scenario, we simply assume that the fault models share the same dip orientation, and vary the dip angle from 50° to 90°. The optimal dip angle is determined when the root-mean-square of fitting residuals between the data and model predictions reaches the minimum.

2.2.2   Inversion method

We use the CosInv software (Wang S et al., 2017, 2019, 2020), which can be freely accessed from http://doi.org/10.5281/zenodo.5098387, to conduct the geodetic inversions. In the sensitivity analysis of dip directions, we discretize the surface plane of each candidate fault model with a grid patch of 5 km × 5 km, and include both the coseismic InSAR and range offset observations. We weight the influence of the InSAR and range offset tracking observations on inversions based on the variance of each dataset, estimated from the nondeforming area (Table S1). The Green’s functions between slip on each sub-fault element and the down-sampled observations (Figure S2; Text S2) are calculated using the elastic half-space dislocation model (Okada, 1985). A left-lateral strike-slip component is allowed only in the strike slip direction with guidance from the USGS and GCMT focal mechanism solutions. To avoid non-physical oscillation in the slip distribution, we adopt a modified Laplacian smoothing operator (Wang S et al., 2017) to constrain the fault motion; this operator is especially suitable for earthquakes with potential surface ruptures (Jiang GY et al., 2013). The optimal smoothing factor is determined from the trade-off curve defining the balance relation between the data misfit and slip roughness (Figure S4a). Once the optimal dip direction and dip angle are resolved, we re-discretize the fault interface with the small patch size of 2 km × 2 km to recover both the co- and post-seismic slip distributions in detail. The above inversion method is applied also to the postseismic slip inversion.

3.   Results

3.1   Co- and Post-seismic Surface Observations

As documented in Figure 2, the Madoi earthquake produced a highly complex coseismic surface deformation field along the mapped fault traces. The maxima observed in the coseismic range offset and InSAR displacement maps are found to be ~1.5 m and ~1.0 m respectively along the line-of-sight (LOS) direction. The coseismic observations from descending and ascending tracks present markedly different fringe patterns on two sides of the fault traces: the north side moves toward the ascending track satellite, while far away from the descending track satellite. This suggests that the earthquake is dominated by left-lateral strike-slip motion. Ground displacement discontinuities are clearly seen in the range offset maps (Figures 2g and 2j), suggesting that slip during the earthquake has reached the Earth’s surface. In addition to the main rupture along fault segment 1, displacement discontinuities are found also along a secondary fault branch (segment 2) near the easternmost edge of the main rupture. This implies that at least two fault segments were involved in the earthquake, which is evidenced from the coherence maps with a cabined lower coherence stripe along the identified fault traces (Figure S1). We note that the aftershock relocations are slightly offset (by ~4 km) from our mapped fault traces (Figures 2h, 2k and S1), which may be due to the limited number and poor spatial coverage of seismic stations, as pointed out by Wang WL et al. (2021). We thus do not use the aftershock relocations to constrain our fault model.

Figure  2.  Co- and post-seismic displacements due to the Madoi earthquake. (a) Observed coseismic InSAR displacements from ascending track 99. (b) Modeled displacements for (a). (c) Fitting residual between (a) and (b). (d−f) are the same as (a−c), but for descending track 106. (g−i) and (j−l) are the same as (a−c) and (d−f) respectively, but for range offset displacements. Inset panel is the zoom-in of the displacement, as outlined by the dashed rectangle. (m−o) and (p−r) are the same as (a−c) and (d−f) respectively, but for the early 24-day postseismic deformation. Note that the coseismic InSAR displacements are rewrapped with a fringe interval of 0 to 0.1 m to highlight the complex deformation of the earthquake. Pink lines are the surface traces identified from range offset maps. Black dots are aftershock relocations (Wang WL et al., 2021).

DownLoad: Full-Size Img PowerPoint

One predominant feature of the coseismic observations is the apparent along-strike variations, as highlighted in the re-wrapped interferograms (Figures 2a and 2d), among which interferometric fringes along east parts (98.5°E−99.5°E) of the fault traces are characterized with more dense and broader distributed fringes than along the west parts. This implies that large coseismic slip heterogeneities could exist on the seismogenic fault and that coseismic strain release may have been highly segmented during the earthquake.

Figures S7a, 2m, and 2p show the early 6-day and 24-day postseismic InSAR deformation respectively, mainly localized along the central part of the coseismic fault surface rupture and clearly showing the fault location. The spatial correlation between the post- and co-seismic deformation suggests that the postseismic motion is also dominated by horizontal movement. The postseismic displacement reaches up to 10 mm in LOS during the time period of 2021/05/26−2021/06/01, and increases to 25 mm during the time period of 2021/05/26−2021/06/19.

3.2   Co- and Post-seismic Slip Distributions

Our analysis on the optimal dip orientation of the seismogenic fault shows that the fitting residual of the data reaches a minimum (~14.4 cm) when the dip angle is 83° for the northeast dipping candidate fault model (Figure S3). However, for the scenario in which dipping is to the southwest, the optimal dip angle cannot be obtained even under the extreme case of dip angle equal to 90°. Thus, we conclude that the seismogenic fault of the Madoi earthquake is dipping to the northwest with a dip angle of 83°. The inverted coseismic slip model shows three major features: (1) The earthquake has broken the Earth’s surface with geodetic-inferred surface ruptures consistent with the coseismic offsets reported from field work (Pan JW et al., 2021; Li ZM et al., 2021; Figure 3a). (2) Major fault slip appears to be concentrated at shallow depths, no more than 12 km, with a maximum slip of ~7.2 m observed at ~7 km depth. (3) Fault ruptures show apparent segmentations with five geodetic asperities (Asperity 1−Asperity 5) identified along the fault. Our optimal coseismic slip model can explain not only the prominent ground motions mapped from the SAR data (Figure 2) but also the GNSS data (Figure 1b). Relatively large fitting residuals in the SAR-mapped coseismic deformation are found only close to the fault traces, which is possibly due to: (1) the complexity of the near-field deformation; (2) the fact that the coseismic deformation includes some early 5-day postseismic deformation; (3) the complex of the real Earth model (Dong J et al., 2021), and (4) the existence of inelastic deformation in the near-field (Pan JW et al., 2021; Li ZM et al., 2021) that cannot be modeled by elastic models. The optimal slip distribution pro- duces a geodetic seismic moment of 1.71 × 10²⁰ N·m, corresponding to an M_W 7.4 event with a shear modulus of 30 GPa.

Figure  3.  (a) Optimal coseismic slip distribution for the Madoi earthquake. Black dots are the aftershock relocations from 2021/05/21 to 2021/05/29 (Wang WL et al., 2021). Red outlined circles are the coseismic offset from field work (Pan JW et al., 2021; Li ZM et al., 2021). To give a better version of slip on all segments, we make a shift to show segment 2. (b) Early 24-day (2021/05/26–2021/06/19) afterslip, inferred from postseismic observations shown in Figures 2m and 2p. Green polygons outline the significant 6-day afterslip distribution, as shown in Figure S7b.

DownLoad: Full-Size Img PowerPoint

In this study, we simply rate afterslip as the only relevant mechanism that contributes most to the observed early postseismic motion, considering that (1) the early postseismic deformation is localized along the fault traces; (2) poroelastic rebound contributes most to the vertical displacement, as it predicts a similar spatial deformation pattern for both the descending and ascending viewing geometries (Figure S8), but the early postseismic motion of the Madoi earthquake is dominated by horizontal displacement as discussed in Section 3.1; and (3) poroelastic rebound produces a deformation that spatially differs from postseismic observations (Figure S8). Assuming that afterslip is the only mechanism that drove the observed postseismic deformation, we find that major afterslips are distributed along the central part of the coseismic fault surface ruptures (Figures 3b and S7b). The 24-day afterslip model exhibits a maximum slip of 7.2 cm, and yields a geodetic seismic moment of 1.67 × 10¹⁸ N·m (corresponding to an M_W 6.1 event), which is only ~0.98% of that released by the coseismic slip. Afterslip seems to have evolved, spatially and temporally: motion on discrete fault patches during the time period of 2021/05/26−2021/06/01 (Figure S7a) evolves towards diffuse motion on shallow fault patches (< 12 km in depth) during the time period of 2021/05/26−2021/06/19 (Figure 3b).

4.   Discussion

4.1   Variable Co- and Post-seismic Slip Models

Different from the recent super-shear strike-slip M 7.5 Palu earthquake in Indonesia (e.g, Socquet et al., 2019), the coseismic rupture of the Madoi earthquake exhibits apparent variation on the fault along strike. Previous studies have attributed such a variation to variable frictional properties and/or interseismic strain partitioning on the fault (e.g., Barnhart et al., 2015; Bassett et al., 2016). Additionally, fault geometry, e.g., with bends or step-overs, could also exert a predominant role in controlling rupture propagation (e.g., Lozos et al., 2011; Sangha et al., 2017; Wang S et al., 2017; Zhang Y et al., 2021). Such a lateral coseismic slip variation is found also in the 2001 M_W 7.8 Hoh Xil earthquake that occurred on the east Kunlun fault (Lasserre et al., 2005). The crucial difference between the two cases is that the Madoi earthquake shows an apparent shallow slip deficit of ~64% in the top upper crust (Figure S5), while this was not resolved for the latter. Constrained by the entire coseismic deformation from InSAR and the range offset observations, the Madoi earthquake’s coseismic shallow slip deficit is a reliable feature, independent from the smoothing factor used in the inversion (Figure S6). The shallow slip deficit in the Madoi earthquake is consistent with that of the M_W 7.2 Baja earthquake (60%) which, similarly, was confirmed from geodetic data collected from observation points providing good coverage of the affected region (e.g., Xu X et al., 2016). Regardless of whether the slip variation along the strike and downdip is related to the velocity-strengthening property of the upper crust, the low level of shear stress, or the complexity of fault geometry, the finding in this study suggests that the seismogenic fault has relatively low structural maturity.

The maximum afterslip of the Madoi earthquake in the first 6-day period (from 2021/05/26 to 2021/06/01) is 5 cm, but it increases only to 7.2 cm in the first 24-day period (from 2021/05/26 to 2021/06/19), suggesting an afterslip decay, as has previously been documented extensively (Marone et al., 1991). The early afterslip is seen also to have evolved in space, from motion on discrete fault patches 11 days following the earthquake (Figure S7b) to wide motion on fault patches shallower than 12 km in depth one month after the shock (Figure 3b); the final result is an obvious spatial overlap of the post- and co-seismic slips (Figure 3). Although this spatial overlap is in contrast to the common understanding that afterslip should be distributed around the seismogenic layer (Scholz, 1998), it can be well explained by non-steady-state frictional behavior (Helmstetter and Shaw, 2009) or residual stress heterogeneities on the fault that can drive continued slow slip on low-slip rupture zone patches (Wang S et al., 2019) as depicted in Figure 3.

Postseismic slip, normally continuing for a period of a few months or more, is widely observed subsequent to large earthquakes (e.g., Ryder et al., 2007; Wen YM et al., 2012). The whole spatiotemporal evolution of the afterslip can be recovered by using geodetic time series data (e.g., InSAR and GNSS), but detailed analysis is beyond the scope of the study; it is expected to be conducted in the future.

4.2   Block-wide Distributed Deformation in the East Bayan Har Block

Based on the most complete GNSS data available so far in the Tibetan Plateau (Wang M and Shen Z-K, 2020; Figure 4a), we calculate the second invariant strain rate (Figure 4b; Auxiliary Text S3), a full strain field of the horizontal strain tensors, using the velocity interpolation for strain rate (VISR) code (Shen Z-K et al., 2015). Unsurprisingly, a prominent rate of 120 nanostrain/yr appears across the Himalayas; a rate of 30 nanostrain/yr is found across the north boundary of the Tibetan Plateau. Notably, a significant strain rate, 30−40 nanostrain/yr, is diffusely distributed across the east Bayan Har block, suggesting a block-wide crustal extrusion in the region (Figure 4b). Viscous flow in the mid-to-lower crust is proposed as a mechanism responsible for the eastward extrusion of the Tibetan Plateau (e.g., England and Molnar, 1997; Royden et al., 1997). Seismic and magnetotelluric data have revealed a mid-to-lower crustal zone of low-velocity and high electric conductivity beneath the east Bayan Har block (Huang SY et al., 2020; Zhan Y et al., 2021). Low shear velocity and/or high electric conductivity can suggest a viscous or partially melted material at depth. Unlike the kinematics of rigid block models (Tapponnier et al., 1982; Avouac and Tapponnier, 1993), viscous models predict that shear deformation tends to distribute over a block-wide zone rather than focus on some large-scale faults for which high slip rates are often expected.

Figure  4.  Block-wide distributed deformation in the east Bayan Har block. (a) GNSS velocity fields with respect to the Eurasian plate (Wang M and Shen Z-K, 2020). (b) Second invariant horizontal strain rate. The color is saturated between 0 and 40 nanostrain/yr to highlight the strain rate in the interior of the Tibetan Plateau. Black contours highlight the strain rate between 25 and 40 nanostrain/yr. (c) GNSS velocity components along profiles shown in panel (a). Green highlights the region of block-wide deformation in the east Bayan Har block.

DownLoad: Full-Size Img PowerPoint

In addition, previous studies show that the inferred shallow depth (<12 km) of the slip distribution of major earthquakes and their associated shallow aftershocks can argue against assignment of shear deformation alone to the elastic strain associated with slips on large-scale faults typically acting as the boundaries of rigid blocks (Zhang PZ et al., 2004). Instead, shear deformation driven by crustal extrusion in terms of viscous flow is most likely absorbed through a block-wide crust. Crustal shortening across the east Bayan Har block is very limited, as evidenced by the fact that crustal velocities parallel to the selected profiles are consistent (Figure 4c).

We note that the block-wide shear deformation shows a remarkable difference from west to east, such as the difference in rate from ~13 mm/yr along profile A-A' to ~6 mm/yr along profile B-B'. This shear deformation variation could reflect a lateral heterogeneity of viscous flow beneath the region. A recent magnetotelluric study in the east Banyan Har block images a major low resistivity zone of less than 10 Ωm at depths from 20 to 50 km near the epicentral area, whereas the thickness of such a low resistivity layer is decreased to ~10 km at ~102°E, and restored to ~30 km at ~104°E where the Longriba fault developed (Zhan Y et al., 2021). The finding of this thickness difference in the low resistivity layer precisely supports the lateral heterogeneity of viscous flow beneath the east Bayan Har block, and reveals a variable diffuse deformation mode in the region. The high electric conductivity at mid-to-lower crust appears to be the means by which viscous driving stress is transferred to the upper brittle seismogenic layer, making the east Bayan Har block a hot spot to host large earthquakes.

In addition, geodetic observations over decadal timescales and geological studies from millennial-scale geomorphic markers both document a decreased slip rate on the east Kunlun fault from west to east (Kirby et al., 2007; Lin A and Guo J, 2008; Diao FQ et al., 2019; Duvall and Clark, 2010; Zhu LY et al., 2020; Zhu Y et al., 2021; Zhu LY et al., 2021), terminating within the thickened crust of the plateau at its easternmost tip. Therefore, any extrusion of the Tibetan lithosphere accomplished by slip along the Kunlun fault must be absorbed by internal deformation of the surrounding plateau (Kirby et al., 2007). Large earthquakes occurring in the interior of the Bayan Har block, e.g., the 2021 Madoi earthquake, bear this out.

4.3   Seismic Hazard on Regional Faults Imposed from the Madoi Earthquake

Given that elastic deformation is widespread in the east Bayan Har block and that the region is characterized by broadly diffused active faults (Figure 1b), the stress transfer from the Madoi earthquake could increase the potential for failure in the rock of nearby faults. Using the PSCMP/PSGRN software (Wang RJ et al., 2006) with crustal structure from CRUST1.0 (Laske et al., 2013), we calculated Coulomb failure stress changes from both the co- and post-seismic (due to viscoelastic relaxation of lower crust and upper mantle) deformation at a depth of 8 km, which is approximately the depth at which major coseismic slip of the Madoi and Hoh Xil earthquakes occurred, to evaluate the stress transfer effect. Detailed configurations for the calculation are presented in Supplementary Materials Text S4.

Our results show that the coseismic stress exerts an increase of more than 0.43 bar on the east part of the Tuosuo Lake fault segment for a length of 50 km, and 0.16 bar for 200 km on the west part of the Maqên−Maqu segment (Figure 5a), results that are consistent with calculations based on a preliminary slip model from USGS (Li YJ et al., 2021a). The postseismic process applies additional load to the above two parts, due to viscoelastic relaxation, increasing the stress on both to ~0.5 bar in 5 years, to ~0.75 bar in 20 years, and to ~1.0 bar in 50 years. These calculations suggest that the two sections now have heightened seismic hazard. In particular, close attention should be paid to the seismic hazard on the Maqên–Maqu segment because it passes through populated areas and surface ruptures are not recorded in its history (Van der Woerd et al., 2002; He WG et al., 2006; Xu XW et al., 2006; Wen XZ et al., 2007; Guo J et al., 2007). Considerable co- and post-seismic loadings of larger than 0.5 bar are also found on the western segments of the south margin fault of the Gadê and the Dari fault, indicating presumably high seismic hazard on those two faults. In contrast, stress loadings on the main Bayan Har fault and the Wudaoliang–Chomsakongma fault are to some extent limited. We also calculated co- and post-seismic stresses at depths of 5 and 15 km, and found that the stress increases and associated spatial patterns are consistent with the 8 km depth results (Figure S10), highlighting the reliability of the above analysis.

Figure  5.  Co- (a) and post-seismic (b−d) Coulomb failure stress changes on regional faults from the Madoi earthquake. The beach balls are the regional M_W ≥ 6.0 historical earthquakes; the Madoi earthquake is highlighted in red. 1 Kunlun fault; 2 South margin fault of Gadê; 3 Dari fault; 4 Main Bayan Har fault; 5 Wudaoliang−Chomsakongma fault.

DownLoad: Full-Size Img PowerPoint

5.   Conclusions

In this study, we have inverted high-quality co- and early post-seismic satellite radar observations to reveal subsurface fault slip variabilities related to the 2021 Madoi M_W 7.4 earthquake. Our results show the following: (1) The coseismic deformation of the Madoi earthquake is very complex, with a surface rupture of ~170 km in length. (2) The coseismic slip distribution exhibits strong heterogeneities with a maximum slip of ~7 m and a geodetic seismic moment of 1.71 × 10²⁰ N·m. (3) Assuming that after- slip is the only mechanism that drives the observed postseismic deformation, the 24-day (2021/05/26–2021/06/19) afterslip model exhibits a maximum slip of 7.2 cm and a geodetic seismic moment of 1.67 × 10¹⁸ N·m. (4) Analysis of the GNSS profiles, strain rate, and magnetotelluric data suggests a block-wide crustal extrusion in the east Bayan Har block. (5) Co- and post-seismic stress loadings from the Madoi earthquake exacerbate seismic hazard on the Maqên–Maqu and Tuosuo Lake fault segments of the Kunlun fault.

Data Availability Statement

Sentinel-1 SAR images and GNSS data can be freely accessed from the Copernicus Open Access Hub (https://scihub.copernicus.eu) and Wang M and Shen Z-K (2020) respectively. Aftershock locations are available through Wang WL et al. (2021). Coseismic offsets determined in the field work are available through Pan JW et al. (2021) and Li ZM et al. (2021). VISR code used to calculate the strain field is available through Shen Z-K et al. (2015).

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Acknowledgments

We thank Editor Wei Leng and two anonymous reviewers for their comments and suggestions. This work is supported by the Natural Science Foundation of Jiangsu Province (Grant No. SBK2020043202) and by Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University (No.19-01-08). Figures were prepared using Generic Mapping Tools (Wessel et al., 2013). Sentinel-1 interferograms are processed with software GAMMA (Wegnüller et al., 2016).

Supplementary Materials for “Resolving co- and early post-seismic slip variations of the 2021 M_W 7.4 Madoi earthquake in east Bayan Har blockwith a block-wide distributed deformation mode fromsatellite synthetic aperture radar data”

Text S1   InSAR and Pixel Offset Tracking Processing

We used Gamma software (Wegnüller et al., 2016) to process and align the single look complex (SLC) image pairs; the topographic phase component was removed based on a 30-m-resolution digital elevation model (Farr et al., 2007). We applied an adaptive filter (Goldstein and Werner, 1998) to the differential interferograms with a small pixel window of 32 × 32 by two times, so that the phase noise could be reduced without distorting the signal. We then unwrapped the interferograms using a minimum cost flow algorithm (Chen and Zebker, 2000). We further detrended the interferograms by fitting a quadratic surface plane in the nondeforming area to account for possible orbital and ionospheric errors (e.g., Wang S et al., 2020). Finally an empirical linear model (e.g., Wang S et al., 2020), regarding topography and phase, was used to mitigate the atmospheric phase delays.

Considering the possibility of surface rupture of the Madoi earthquake, InSAR observations could encounter near-field decorrelation. Thus, we used pixel offset tracking to capture coseismic surface displacements in the near field, based on pixel-by-pixel cross correlation of SAR intensity (Michel et al., 1999; De Zan, 2014). Here, we calculate only the range offset displacements because the theoretical accuracy of Sentinel-1 azimuth offsets (~1.40 m) is insufficient to map the coseismic deformation caused by such a rough east-west-oriented strike-slip event as the Madoi earthquake for which the accuracy of the range offset can be ~0.23 m. During the calculation of range offset, we set the search window size of SLC offsets to be 300 × 60 pixels to specify a square searching range. We also set the window size for calculating the cross-correlation function between pixels and the acceptance threshold of cross-correlation coherence to 32 × 32 pixels and 0.1, respectively. A median filter with a window size of 9 × 9 pixels is further used to reduce noise in the estimated offsets.

Text S2   Co- and Post-seismic Observations Downsampling

High-resolution InSAR and range offset displacement maps are preferable for examining the rupture details of an earthquake, especially for constraining near surface ruptures. However, inputting all of them to the modeling process is to some extent impractical. To reduce redundancy and facilitate the co- and post-seismic modeling, the original high-resolution co- and post-seismic ground deformation can be downsampled without degrading the contained information as they are highly correlated spatially. Accordingly, we apply a uniform downsampling method to the near-, medium- and far-field observations respectively, with different configuration parameters (Table S2). Different from the quad-tree and model resolution-based downsampling schemes, our method here not only overcomes the problem of oversampling data in areas with large phase noise, but does so without requiring an initial slip model. Figure S2 shows our downsampled observations for the co- and post-seismic surface deformation, which are seen to preserve the surface rupture information in the near-field regions in great detail.

Text S3   Strain Rate Calculation

Geodetic estimates of crustal strain rates could give us a more specific idea about how fast the lithosphere is being deformed in response to tectonic plate movement. They can also provide us detailed information on how the continents deform, e.g., a diffusely distributed strain rate usually means the crust is deforming in accord with a continuum model. To examine how the Bayan Har block deforms, we here interpolate the horizontal GNSS velocity from Wang M and Shen Z-K (2020) to derive a smoothly distributed strain rate field, using VISR software developed by Shen Z-K et al. (2015). A priori information, in the form of weighted smoothing, must be used to facilitate the solution. The degree of smoothing reflects the spatial scale over which strain rates would be averaged and is closely related to the spatial density and coverage of the GNSS sites from which the data are obtained. At a given site, the horizontal velocity field in its vicinity can be approximated by a bilinear function model consisting of translation, rotation, and strain rates. These unknown parameters can be resolved through a least squares inversion when the fitting residual of the GNSS velocity reaches a minimum value. To determine the optimal smoothing distance D, we introduce a quantity W_t to denote the threshold of the sum of the reweighting coefficients $W_t={\text{Σ}}_{i}{L}_{i}\times{V}_{i}$ of the data. ${L}_{i}$, the distance dependent weighting, is assumed to be a Gaussian function ${L}_{i}=\rm{exp}\left(-\Delta R^2/{D}^{2}\right),$ in which $\Delta R$ is the distance between an interpolation site and a GNSS site. ${V}_{i}$, the spatial weighting, depends on the area of the Voronoi cell. The spatial smoothing parameter D varies for each interpolation site. More information about VISR can be found in Shen Z-K et al. (2015). In this study, we set W_t = 24; the calculated results of the spatial smoothing distance and the second invariant of horizontal strain rates are shown in Figures S9 and 4b respectively.

Text S4   Lithospheric Structure Used for the Calculation of Co- and Post-seismic Coulomb Failure Stress Changes

While each event produces a net reduction of stress on the involved fault, events can also result in stress increases on nearby faults (King et al., 1994). With further tectonic loading it seems logical that such sites of stress rise should be the foci of future events. To examine the impact of the Madoi earthquake on potential failures in the regional faults, we calculated both the co- and post-seismic Coulomb failure stress changes at 8 km depth using PSCMP/PSGRN software (Wang RJ et al., 2006). To represent the lithospheric structure, we used a three-layer model consisting of an elastic upper crust layer overlying a lower crust and a half-space upper mantle. The effective friction coefficient was assumed to be 0.4, a typical value commonly used for strike-slip faults (e.g., Liu G et al., 2019). The thickness of the elastic layer was set to 20 km based on the coseismic slip distribution of the Madoi earthquake. We assumed a simple linear Maxwell rheology for the lower crust and upper mantle, as have previous studies (e.g., Liu G et al., 2019; Diao et al., 2019). Crustal velocity parameters were taken from CRUST1.0 (Table S3). The viscosity of the lower curst was set to 3.0 × 10¹⁸ Pa·s based on interseismic studies on the Kunlun fault (Diao FQ et al., 2019); that of the corresponding upper mantel was set to 1.0 × 10²⁰ Pa·s, as deduced from the postseismic modeling of the Hoh Xil earthquake (Wen YM et al., 2012).

  S1.  Sentinel-1 data used to measure the co- and post-seismic ground deformation triggered by the Madoi earthquake.

Deformation Type	Track	Reference Date (YYYMMDD)	Repeat Date (YYYYMMDD)	σ† (cm)
Coseismic InSAR	99 (Ascending)	20210520	20210526	10
	106 (Descending)	20210520	20210526	9
Coseismic range offset	99 (Ascending)	20210520	20210526	16
	106 (Descending)	20210520	20210526	9
Postseismic InSAR	99 (Ascending)	20210526	20210619	2.2
	106 (Descending)	20210526	20210601	0.9
	106 (Descending)	20210526	20210619	1.7
Notes: †σ is the standard deviation estimated from pixels by masking out the median- and near-field regions shown in Figure S2.

 | Show Table

DownLoad: CSV

  S2.  Block size used in the downsampling InSAR and range offset observations.

	Near field	Medium field	Far field
InSAR	64	128	512
Range offset	32	64	256

 | Show Table

DownLoad: CSV

  S3.  Lithospheric structure used for the calculation of co- and post-seismic Coulomb failure stress changes.

Layer	Depth (km)	VP (km/s)	VS (km/s)	Density (km/m3)	$\eta$ (Pa·s)
Upper crust	0.0	6.0000	3.5200	2720.0	0.0E+00
	20.0	6.0000	3.5200	2720.0	0.0E+00
Lower crust	20.0	6.6000	3.8200	2850.0	3.0E+18
	60.0	6.6000	3.8200	2850.0	3.0E+18
Upper mantle	60.0	8.1900	4.5400	3370.0	1.0E+20
Note: VP and VS are the seismic velocities of P and S waves respectively. $\eta$ represents the viscosity of the lithosphere.

 | Show Table

DownLoad: CSV

  S1.  (a) and (b) are the coherence maps of the coseismic InSAR interferograms from the ascending and descending tracks respectively. Red dots are the aftershock relocations from 2021/05/21 to 2021/05/29 (Wang WL et al., 2021). Dashed green lines are the identified surface traces of the seismogenic fault of the Madoi earthquake. We note that coherence exhibits lower values along the identified fault traces.

DownLoad: Full-Size Img PowerPoint

  S2.  Downsampled observations for the co- and post-seismic ground deformation, with track number (e.g., T099) and time period (e.g., 2021/05/20–2021/05/26) of the interferograms indicated

.

DownLoad: Full-Size Img PowerPoint

  S3.  Testing two scenarios of fault models with variation of dip orientations.

DownLoad: Full-Size Img PowerPoint

  S4.  Trade-off curve between the misfit of co- and post-seismic observations and slip roughness for the optimal seismogenic fault dipping to the northeast with a dip angle of 83°. (a) and (b) are results for the co- and post-seismic slip models respectively.

DownLoad: Full-Size Img PowerPoint

  S5.  Variable slip features of the Madoi earthquake. (a) Distribution of the cumulative coseismic slip. Horizontal axis represents the normalized cumulative slip computed by integrating slip along strike. This panel is used to highlight the observed shallow slip deficit. (b) Seismic moment released along strike.

DownLoad: Full-Size Img PowerPoint

  S6.  (a) Strike- and dip-slip components of the optimal coseismic slip model, with the total slip distribution shown in Figure 3a. Note that only the strike-slip component is allowed on segment 2. (b) Total slip distribution estimated with a smoothing factor of 0.7. (c) Same as (b) but with a smoothing factor of 19.3.

DownLoad: Full-Size Img PowerPoint

  S7.  Observed early 6-day postseismic InSAR LOS displacements from descending track 109, with primary and secondary images acquired at 2021/05/26 (YYYY/MM/DD) and 2021/06/01 respectively. (b) Modeled displacements from the best fitting afterslip model shown in panel (d) for (a). (c) Fitting residual between (a) and (b). Pink lines represent the surface trace identified from the range offset maps. Black dots are the aftershock relocations (Wang WL et al., 2021).

DownLoad: Full-Size Img PowerPoint

  S8.  Predicted surface LOS displacements due to poroelastic rebound calculated by differencing the coseismic displacements modeled under fully drained and undrained conditions respectively (Jónsson et al., 2003). The Poisson’s ratio for the drained and undrained conditions are assumed to be 0.25 and 0.28 respectively.

DownLoad: Full-Size Img PowerPoint

  S9.  Estimated smoothing distance with W_t = 24. Gray dots denote locations of the GNSS sites used in this study.

DownLoad: Full-Size Img PowerPoint

  S10.  (a−d) Same as Figures 5a-5d, but for the co- and post-seismic Coulomb failure stress changes at 5 km. (e−h) Same as Figures 5a-5d, but for the co- and post-seismic Coulomb failure stress changes at 15 km.

DownLoad: Full-Size Img PowerPoint