Bayesian inversion and quantitative comparison for bilaterally quasi-symmetric rupture processes on a multisegment fault in the 2021 M_w7.4 Maduo earthquake

SUMMARY

On 2021 May 21, the |$M_{w}7.4$| Maduo earthquake occurred in the southern Qinghai Province, China. This earthquake ruptured approximately 160 km along the Kunlunshankou–Jiangcuo fault, an east–west trending fault located in the middle of the Bayan Har Block. The seismogenic fault exhibits an apparent simple geometry, characterized by fault branches to the east and a splay fault to the west. Despite the apparent simplicity of the fault’s structure, a noteworthy level of variability and inconsistencies persist in the representations of fault geometry in published rupture models of the earthquake. Our study employs a Bayesian approach to elucidate both the fault geometry and kinematic rupture parameters of the earthquake. We use 3-D displacements obtained from synthetic aperture radar images and teleseismic data to quantify the rupture kinematics. We then conduct three separate finite-fault inversions using individual data sets, and perform a subsequent joint inversion for a comprehensive analysis. Additionally, we employ teleseismic back projection to complement the kinematic analysis of the earthquake rupture. Our results unveil a nearly symmetric bilateral rupture in the east–west direction, characterized by an average rupture speed of approximately 2.7 km s⁻¹. The rupture to the east displays a heightened level of complexity, manifested in at least five discernible stages, whereas the rupture to the west is comparatively simpler. The eastward rupture directly triggered the southern branch of the bifurcating fault, with a notable delay of approximately 3 s on the northern branch. Several studies have presented coseismic slip models for the earthquake. An analysis of variability among 10 slip models, including our preferred model, highlights that fault geometry and inversion strategy (e.g. fault discretization, smoothing factor) contribute to considerable variability in both slip magnitude and slip extent on the fault, despite similar data types being used in the inversions. Furthermore, the finite-fault model acquired through slip inversion plays a crucial role in calculating Coulomb failure stress change (⁠|$\Delta {\rm CFS}$|⁠) transmitted from the source fault to neighbouring receiver faults. Understanding how the variability in slip models influences |$\Delta {\rm CFS}$| calculations is essential for conducting comprehensive analyses in seismic hazard studies. Our findings highlight that discrepancies in fault geometry contribute to the variance of |$\Delta {\rm CFS}$| in the regions delineating positive and negative stress change. Meanwhile, variability in slip magnitude substantially impacts the variability of |$\Delta {\rm CFS}$| in the vicinity of the source fault. Furthermore, our analysis of |$\Delta {\rm CFS}$| calculations using our preferred slip models indicates that a major event on the Maqin-Maqu segment, a well-recognized seismic gap on the East Kunlun Fault (EKF), could potentially be advanced in time.

1 INTRODUCTION

An |$M_{w}7.4$| earthquake struck the Maduo County, Qinghai Province, Southern China on 2021 May 21. This event occurred in the northeastern part of the Tibetan Plateau within the Bayan Har Block (BHB) that separates the Qaidam-Qilian block (QQB) in the north from the |$\sim$|1600 km East Kunlun Fault (EKF) and from the Qiangtang block (QTB) in the south by Xianshuihe Fault (Fig. 1a). The two seismically active block boundary faults have hosted several large historic earthquakes with magnitudes |$M_{w}\ge 7.0$| in the recent past, including the 1937 |$M_{w} 7.8$| Huashixia earthquake, the 2001 |$M_{w} 7.8$| Kunlun earthquake, the 2008 |$M_{w} 7.1$| Yatian earthquake and the 2008 |$M_{w} 7.9$| Wenchuan earthquake (Fig. 1a). Unlike these block boundary earthquakes, the left-lateral strike-slip |$M_{w} 7.4$| Maduo earthquake initiated on and ruptured the secondary Kunlunshankou-Jiangcuo fault (KLSK-JCF) (Li et al. 2021; Xie et al. 2022; Zhang et al. 2022c) in the eastern-central area of the BHB (Fig. 1a). It partially overlaps with the western segment of a previously recognized but un-named fault (see HimaTibetMap data base; Taylor & Yin 2009; Styron et al. 2010). Field investigations show that the Maduo earthquake caused notable ground shaking and damaged structures located at the eastern end of the coseismic surface ruptures (Zhang et al. 2022c).

Unlike the surrounding large boundary faults that are primarily responsible for large earthquakes (Fig. 1a), no major earthquakes have occurred within the BHB near the JCF since the 1947 |$M_{w} 7.3$| Dari earthquake. The estimated low fault slip-rate (⁠|$1.2 \pm 0.8$| mm yr⁻¹) (Zhu et al. 2021b) and limited fault length make the faults within the BHB less likely to generate large earthquakes. The block-bounding fault zone of the Tibetan plateau exhibits exceptional kinematic features and demonstrates higher strain rates (⁠|$\gt 20$| nanostrain yr⁻¹) (Wang et al. 2019; Wang & Shen 2020) than the estimated on JCF of |$\lt 20$| nanostrain yr⁻¹ (Zhao et al. 2021). However, the Maduo earthquake highlights that regions with low and difficult-to-quantify interseismic strain rates, below the resolution of current instruments, still have the potential to generate large earthquakes.

Investigating the spatiotemporal rupture evolution of large earthquakes is crucial to understand the physical processes that result in the surface deformation and the potential correlation between the rupture processes and the seismotectonic conditions, which is particularly important for seismic hazard assessment. However, it is well known that finite-fault rupture models for a single event inferred from different data sets and/or by different research teams show distinct discrepancies (Mai et al. 2016). Therefore, it is important to quantify the variations of inferred rupture models and their corresponding effects on, for example, the estimated Coulomb failure stress change (⁠|$\Delta {\rm CFS}$|⁠) in the region.

The robustness of finite-fault models is critical to earthquake source modelling as it determines the reliability and accuracy of the inferred rupture characteristics. A primary concern in finite-fault inversion is the non-uniqueness of solutions, particularly if relying solely on teleseismic data sets, where different finite-fault models may translate to fairly similar seismograms (López-Comino et al. 2015). Therefore, combining various data sets is crucial to reduce the non-uniqueness of the inversion solutions and help to better constrain the earthquake source models (Weston et al. 2014). Additionally, the variability in slip models is influenced by inversion approaches, modelling assumptions and parameter choices (see e.g. Weston et al. 2012; Mai et al. 2016). Consequently, the derived rupture models can be viewed as samples from a broad spectrum of plausible models or as distinct representations from different parameter spaces. These variations result in different solutions for the same seismic event. Furthermore, a synthetic study by Frietsch et al. (2019) reveals that different inversion methodologies can lead to disparate solutions. Specifically, simultaneous inversion of fault parameters effectively recovers input parameters, while sequential inversion introduces substantial errors (Frietsch et al. 2019).

To thoroughly investigate the spectrum of potential rupture scenarios of an earthquake, it is imperative to establish the distribution of all credible source models, considering the trade-offs and uncertainties associated with each rupture parameter. To achieve this, various nonlinear optimization approaches have been developed, such as simulated annealing (e.g. Hartzell et al. 1996; Bouchon et al. 2000; Delouis et al. 2002; Liu & Archuleta 2004), neighbourhood searches (e.g. Vallée & Bouchon 2004) and genetic algorithms (e.g. Emolo & Zollo 2005). However, previous publications then often relied on a single best-fitting model or the mean model of a limited number of models. As a result, additional uncertainty quantification is required to accurately assess model variability and model robustness. In light of this, Bayesian inference techniques were proposed to infer kinematic finite-fault rupture parameters (Monelli & Mai 2008; Minson et al. 2013; Dettmer et al. 2014; Vasyura-Bathke et al. 2020) and dynamic source parameters (Gallovič et al. 2019), along with their associated uncertainty estimates.

For the 2021 Maduo earthquake, different types of data sets have been used to estimate slip models: for example, coseismic slip models constrained by 3-D displacements derived from synthetic aperture radar (SAR) images (e.g. Chen et al. 2021; He et al. 2021; Jin & Fialko 2021; Zhao et al. 2021; Wang et al. 2022a; Liu et al. 2022b), static GNSS offsets and SAR images (e.g. Guo et al. 2022; Hong et al. 2022; Li et al. 2022b) or using a combination of multiple data sets including teleseismic P waves, SH waves, SAR images, static GNSS offsets and high-rate GNSS (e.g. Chen et al. 2022; Wang et al. 2022b). However, most published models rely on conventional optimization approaches that only yield a single solution to best fit the used observation data within the selected parameter ranges, therefore, they lack rigorous uncertainty quantification. The uncertainty in slip models is of critical importance to understand the implications of rupture-model variations, for example on the reliability of |$\Delta {\rm CFS}$| calculation, which is used to characterize the stress conditions at which failure can occur and to investigate stress interactions between earthquakes (Woessner et al. 2006, 2012; Dutta et al. 2018).

To investigate the coseismic slip of the 2021 Maduo earthquake and its uncertainty, we use geodetic data and teleseismic waveform recordings, either jointly or separately, to perform Bayesian inversions. Before doing so, we first infer the fault geometry parameters simultaneously using geodetic and teleseismic data to better constrain the solution. We further use two global seismic arrays to perform back projection (BP) to track the spatiotemporal evolution of high-frequency radiation generally associated with the rupture front propagation (e.g. Ishii et al. 2005; Krüger & Ohrnberger 2005).

Furthermore, we compare our slip models with previous studies and conduct a comprehensive variability analysis of different slip models of the Maduo earthquake. We also investigate how the slip variability persisting among the slip models and the corresponding uncertainties affect |$\Delta {\rm CFS}$| on surrounding faults due to the main shock of the Maduo earthquake.

2 DATA AND METHOD

Utilizing a sequential Monte Carlo sampling algorithm, we follow a Bayesian approach to infer both the seismogenic fault geometry and the associated rupture parameters of the Maduo earthquake. This approach allows us to estimate the geometry uncertainty of the seismogenic faults that hosted the earthquake. Below, we present the data sets used to constrain the model parameters, the detailed parametrization steps and then the results of the fault and rupture parameter inference.

Multiple data sets, including geodetic and seismic data from local, regional or teleseismic sources, offer distinct wavelength and frequency resolution for capturing the rupture process, whereby each data set exhibits different sensitivity to the Earth’s structure. Frietsch et al. (2019) suggest that constraining the downdip segmentation when inferring the fault geometry is challenging without employing local/regional seismic data. Therefore, such data is necessary to detect the rupture of both shallow and deep parts of a potential listric fault (Frietsch et al. 2021). Teleseismic data, due to their distance, are less sensitive to lateral heterogeneities in the Earth’s crust. Considering these factors, along with the accuracy of the 3-D Earth’s subsurface model and the computational constraints of incorporating 3-D heterogeneities in finite-fault inversion, most source studies using teleseismic data adopt a 1-D global Earth model.

Our analysis of the Maduo earthquake’s fault geometry utilize 3-D surface displacements and teleseismic P-waveform data sets, with the addition of SH waves for kinematic finite-fault inference. Local/regional seismic data was excluded from both inversions due to its potential sensitivity to 3-D local/regional lateral heterogeneity compared to teleseismic data, necessitating precise local/regional 3-D velocity models which we did not have access to. A potential downdip segmentation of fault geometry was not considered due to the unavailability of local/regional seismic data sets. We also employ the back projection method (Ishii et al. 2005; Krüger & Ohrnberger 2005) to analyse the spatiotemporal evolution of the high-frequency seismic radiation. This approach complements rupture models obtained from kinematic finite-fault inversion, providing further insights into the earthquake rupture process.

2.1 Geodetic data

We utilize ascending and descending Sentinel-1 and Advanced Land Observation Satellite-2 (ALOS-2) SAR data, specifically four pairs of SAR images, to quantify the coseismic surface deformation caused by the Maduo earthquake. Our approach involves differential interferometric SAR (DInSAR) (Gabriel et al. 1989), pixel-offset tracking (Michel et al. 1999), multiple aperture interferometry (MAI) (Bechor & Zebker 2006) and burst overlap interferometry (BOI) techniques (Grandin et al. 2016). As a result, we obtain line-of-sight (LOS) displacements from DInSAR and azimuthal (AZI) displacements from MAI and BOI. Pixel-offset tracking allowed us to derive both an independent set of LOS and AZI displacements. We then use the strain model and variance component estimation method (Liu et al. 2018, 2022a) to combine all surface displacement observations to produce an accurate representation of the 3-D ground surface displacements (i.e. east–west, north–south and vertical displacements; see Fig. 2).

The 3-D surface displacement images clearly reveal a distinct and contrasting offset along and across the fault line, indicating a well-defined surface rupture extending toward the northwest near the Jiangcuo fault. Fig. 2 shows dominant left-lateral displacement of the ruptured fault, with maximum slip of 1.7 and 2.4 m in the east and west directions, respectively. The 3-D displacements reveal noticeable variations in the horizontal and vertical displacement amplitudes along the fault. The northern side of the fault experiences larger displacements compared to the southern side. Notably, along the fault trace, the epicenter region displays a relatively smaller fault offset than areas further to the west and east.

We downsample the full-resolution 3-D ground displacement fields (Fig. 2) for computational efficiency using the 2-D quadtree algorithm (Jónsson et al. 2002) implemented in the KITE software (Isken et al. 2017) to ensure that more refined sampling near the rupture zone with higher coseismic strain release and coarser sampling farther away. This strongly reduces the number of data points from millions to 1312, 319 and 374 for the E–W, N–S and vertical displacements, respectively (Figs 5a, d, g). We then calculate the covariance of each displacement component by employing an approach proposed by Sudhaus & Jónsson (2009) and use the estimated observation error variances to weight each component for the slip model inversion.

2.2 Teleseismic data

We use broad-band teleseismic waveform data from the Australian National Seismograph Network (AU) and the European Array (EU) (see Fig. 3a) to perform back projection. We download teleseismic recordings from the Incorporated Research Institutions for Seismology (IRIS) web service and remove the instrument response. To homogenize the time series, we downsample the corrected waveforms to a uniform 25 Hz sampling rate. We then apply cross-correlation for 10 s of filtered (0.5–1 Hz) waveforms around the P-wave arrival to align the P-wave train. This process enhances the coherency of the seismic traces and corrects the effects of 3-D structure heterogeneities that are not considered in the theoretical traveltimes calculation using the global 1-D velocity model. Only stations with higher average cross-correlation coefficients (⁠|$\ge 0.9$|⁠) are used. We then perform the back projection in the frequency band 0.5–1 Hz.

For the Bayesian inversions of the fault geometry and the spatiotemporal slip evolution, we select 13 teleseismic stations with a wide azimuth coverage, within epicentral distance range |$30^{\circ }-60^{\circ }$| (Supplementary Fig. S1a). The downloaded raw waveforms are transformed to ground displacements and rotated into the radial (R), transverse (T) and vertical (Z) components. For fault geometry inference, we apply a two-pass Butterworth filter (order 4), frequency range of 0.01 to 0.1 Hz. Specifically, we only utilize the vertical (Z) component to perform fitting of a 80 s time window of P waveforms, extending from 20 s before to 60 s after the first P-wave arrival. For kinematic finite-fault inversion, we use P- and SH-waves recordings to constrain the rupture parameters (Figs 3b and c); here, we consider a higher frequency range of 0.05–0.5 Hz. We only use the vertical (Z) and transversal (T) components to fit the corresponding P and SH waveform, respectively, over a 100 s time window covering 30 s before and 70 s after the first P-wave arrival.

2.3 Back projection method

Since it was introduced by Ishii et al. (2005) and Krüger & Ohrnberger (2005), the BP method has been used in many seismic-source studies to track the spatiotemporal evolution of seismic radiation sources during earthquake ruptures (e.g. Meng et al. 2012; Li & Ghosh 2017; Mai et al. 2023; Zhang et al. 2023; Vasyura-Bathke et al. 2024). The source of high-frequency seismic radiation being tracked could be attributed from many different physical processes, such as a rapid change in rupture velocity or slip (e.g. Meng et al. 2011), heterogeneity in slip rate or variation in fault geometry (Madariaga et al. 2006; Li et al. 2022a). The different frequency ranges used in BP may track different sources. In particular, radiation of low-frequency waves (0.1–0.5 Hz) mostly align with areas of large coseismic slip (e.g. Melgar et al. 2016; Yin et al. 2016), while radiated high-frequency waves (0.5–1 Hz) align with the edges of large slip areas and positive stress change (e.g. Yin & Denolle 2019).

As the BP method has poor depth resolution (Ishii et al. 2005; Li et al. 2024), we use a horizontal source plane at a fixed depth of 10 km. The source region is bounded by |$97.13^{\circ }$| to |$100.13^{\circ }$| in longitude and |$33.4^{\circ }$| to |$26.4^{\circ }$| in latitude. We discretize the plane into a grid of potential seismic-radiation sources with uniform |$0.015^{\circ }$| spacing in both latitude and longitude. Different from Ishii et al. (2005), which used stacked beams, here we use a semblance term, a normalization of beam power by the energy summed over all considered seismograms (Roessler et al. 2010). The semblance converges to |$1/N$| (where N is the number of seismic stations) for uncorrelated noise, but noticeably higher for coherence signals.

2.4 Bayesian inversion to infer the fault geometry and kinematic finite-fault model

The Bayesian Earthquake Analysis Tool (BEAT; Vasyura-Bathke et al. 2019, 2020) has been applied for both fault geometry and finite-fault model inference. One advantage of BEAT is that it can estimate the posterior probability density function (PDF) of each model parameter |${\bf m}$| (⁠|$P\left({\bf m}|{\bf d}_{obs}\right)$|⁠), based on the Bayesian Theorem (Bayes & Price 1763) under the assumption of Gaussian distributed residuals |${\bf r}_{k}\left({\bf m}\right)$|⁠, expressed as:

where |$P\left({\bf m}\right)$| is a prior PDF, |$L\left({\bf m},\sigma _{k}\right)$| is a likelihood function and |${\bf r}_{k}\left({\bf m}\right)={\bf d}_{k,obs}-{\bf d}_{k}\left({\bf m}\right)$| denotes the residual, with |${\bf d}_{k}\left({\bf m}\right)$| as the predicted synthetic of data set k (i.e. teleseismic data set and coseismic 3-D surface deformation). |$C_{k}$| is the covariance matrix describing noise statistics of residual errors of each data set, and |$\sigma _{k}$| is a hierarchical scaling parameter to estimate the standard deviation of the residuals. We use an advanced sequential Monte Carlo (SMC) algorithm (Del Moral et al. 2006; Vasyura-Bathke et al. 2020) to sample the posterior PDF. It samples from the simulation of a series of ‘transitional’ intermediate PDFs from the prior to posterior as shown in eq. (1). The algorithm also enables likelihood tempering (or annealing), combining importance sampling and resampling, and parallel computing. This sampling method is particularly useful for exploring a high-dimensional and multimodal posterior PDF (e.g. Neal 1996; Del Moral et al. 2006, 2007; Herbst & Schorfheide 2014; Nguyen et al. 2016; Amaya et al. 2021).

Within Bayesian inference, parameter and uncertainty estimates are notably influenced by data errors originating from selections made in Earth velocity models and source parametrization. To implement the method effectively, assumptions on noise statistics (⁠|$C_{k}$| in eq. 2) are necessary. This involves specifying distribution types and corresponding parameters that require appropriate assumptions or estimations. We employ data covariance estimation to mitigate errors arising from both measurement inaccuracies and errors in theoretical assumptions. This type of covariance is independent of the uncertainties in the assumed velocity model. In Bayesian joint inversion, estimating |$C_{k}$| is particularly important considering the contributions of multiple data sets to the solutions. Therefore, we estimate a specific covariance matrix for each data set. These matrices provide reliable statistical measures of how significant each data set contributes to the solution. In addition, these statistical measures also provide more objective weights to quantify the uncertainty.

2.4.1 Fault geometry inference

For large crustal earthquakes, a well-constrained fault geometry is required for accurate coseismic slip inversion. However, the geometry for any given event could be inconsistently estimated when analysed with different data sets and different approaches, which in turn, may have large uncertainties, especially for complex rupture events. In our study, we first use the aforementioned high-resolution coseismic deformation data to define the surface fault traces. The surface rupture traces show a linear delineation in the central portion of the coseismic rupture that separates two different deformation characteristics, with one splay fault in the west, and a bifurcation in the east. Thus we divide the fault traces into four segments, shown as F-1, F-2, F-3 and F-4 in Fig. 1(b). F-1 is the westernmost segment, stretching 30 km in the strike direction of 274.5|$^{\circ }$|⁠. F-2 is the longest central segment, extending 100 km in a strike direction of 285.2|$^{\circ }$|⁠. F-3 and F-4 are the bifurcation faults in the east, with the 35 km F-3 to the north in the strike direction of 265.4|$^{\circ }$|⁠, and the 20 km F-4 to the south in a strike direction of 105.5|$^{\circ }$|⁠. More segments could be used to represent the rupture traces better, but this would increase the computational cost required for the Bayesian source inversion without adding more information due to the band-width limitation of the seismic data.

We treat each fault segment as a rectangular element with uniform slip. For each fault segment, we estimate its geometry and fault slip kinematics parameters: dip, rake, width and slip. In addition, assuming a half-cosine source time function (STF; Lay et al. 2010), we estimate nucleation points and their corresponding nucleation times and slip duration of each fault segment. We do not estimate orientation, length and depth of each fault segment in the inversion, but use a priori information from surface data to constrain these three parameters. The strike orientation and length of fault segments are determined from the rupture traces mapped from the 3-D surface displacements, while the top of all fault segments is set to 0 km considering surface-rupture traces of the earthquake. Each segment’s strike and length are shown in Fig. 1(b); for all four segments, there are eight geometrical parameters (dip and width), eight slip parameters (slip orientation and total slip) and eight parameters of source time function (rupture start time and duration) and eight parameters of nucleation points. To partially address the inaccuracies in the Green’s function resulting from the lateral heterogeneity of Earth structures, we also estimate a time-shift for the P waveform at each seismic station (Mustać et al. 2020; Vasyura-Bathke et al. 2021). We perform a joint Bayesian inversion utilizing 3-D surface displacements and teleseismic P-wave data with the SMC sampling algorithm considering the full data error variance–covariance matrix of both the geodetic and seismic data sets (with the seismic noise measured within a 5 s long time window prior to the P-wave arrival) to estimate the hierarchical scaling for both surface displacements and seismic data.

For the inversion, we adopt a prior PDF to infer the fault geometric parameters of each segment, assuming that fault parameters follow a uniform distribution for a pre-defined value range. For example, we set the prior dip angle ranging from 70|$^{\circ }$| to 110|$^{\circ }$| (⁠|$\lt 90^{\circ }$| indicates a north-dipping segment, and |$\gt 90^{\circ }$| a south-dipping one and reverse for F-4). This dip angle range should sufficiently cover the likely dip angles, given that the E–W surface displacements exhibits almost symmetric behaviour across the fault surface traces, indicating a nearly vertical fault plane. For fault width, we define the prior range between 0 and 30 km considering the aftershock distribution (Wang et al. 2021) and the upper crustal thickness of the BHB block (Wang et al. 2008). The solution of uniform slip of all fault segments is assumed to be between 0.1  and 8 m. Furthermore, for the rake angle parameter of each segment, we define the range value between |$-90^\circ$| and |$90^\circ$|⁠.

2.4.2 Kinematic finite-fault inference

For finite-fault inference (FFI), we use the geometry of fault segments obtained from the maximum a posteriori (MAP) solution of the fault geometry inference (see MAP solutions corresponding to SC-1 in Supplementary Table S1 for more details). For fault width, we extend the fault to 30 km in downdip direction for all segments. The top edges of all the faults are assumed to be at the Earth surface according to the coseismic surface displacements and field investigations (Pan et al. 2021; Ren et al. 2022; Xie et al. 2022). The fault model is divided into 222 sub-faults of a uniform area of 5×5 km|$^2$|⁠. Different from the conventional slip inversion, which assumes constant rupture speed and rise time, the Bayesian inversion in the BEAT allows for solving the rupture parameters (strike-slip, dip-slip, rupture speed and rise time) on each sub-fault, including the hypocentre location. All rupture parameters are sampled within their pre-defined prior boundary values with equal probability under a uniform distribution assumption. For slip parameters, earthquake scaling relations of McGarr & Fletcher (2003) and Thingbaijam et al. (2017) suggest a maximum slip of |$\sim$|10 m for a magnitude |$M_{w}7.4$| event. Therefore, we set the slip ranges to [−0.05,10] m for strike-slip and [−10.0,10.0] m for dip-slip. This enables us to appropriately infer the sinistral fault movement with varying slip along both strike and dip directions. The lower and upper boundaries of prior rupture parameters are presented in Supplementary Table S2. We apply Laplacian smoothing for the slip distribution during the inversion (Fukuda & Johnson 2008; Vasyura-Bathke et al. 2020).

To compute synthetic seismograms and static displacements, we use a layered elastic half-space model by combining the CRUST2 Earth model (Bassin et al. 2000) for shallower crustal structure and the AK135 (Kennett et al. 1995) for the deeper structures. The synthetic teleseismic waveforms are generated by convolving the sub-fault slip-rate function of a sampling rate of |${\rm d}t = 0.5$| s with the elementary Green’s Functions (GFs) (Heimann et al. 2019). The GFs are pre-calculated using a spherical Earth model code (QSSP; Wang et al. 2017) for each sub-fault and station pair. Additionally, we calculate the coseismic GFs and resulting static displacements using the PSGRN/PSCMP program (Wang et al.2006).

3 RESULTS

After describing the data and methods we use to estimate fault geometry and rupture kinematics of the Maduo earthquake, we now present our results in detail before discussion further implications of the findings.

3.1 Geometry of fault segments

Estimating the precise fault geometry is the initial step before constraining the kinematic rupture process. Starting with the dip of each fault segment, we observe that the 1-D marginal posterior distributions of fault dip (Fig. 4) shows a high probability at 80|$^{\circ }$|–85|$^{\circ }$| (north dipping) for fault-segment F-1, and 83|$^{\circ }$|–87|$^{\circ }$| (north dipping) for segment F-2. For segments F-3 and F-4, we find 108|$^{\circ }$|–110|$^{\circ }$| (or 70|$^{\circ }$|–72|$^{\circ }$| dipping to the south) and 90|$^{\circ }$|–107|$^{\circ }$| (or 73|$^{\circ }$|–90|$^{\circ }$| dipping to the north), respectively. Thus, the MAP solutions of fault-segment dip show nearly vertically dipping fault segments, F-1, F-2 and F-4 (dipping to the north by 84.2|$^{\circ }$|⁠, 85.2|$^{\circ }$| and 87.6|$^{\circ }$|⁠). Fault segment F-3, on the other hand, is found dipping |$\sim 71^{\circ }$| to the south. In terms of fault width, we find that F-1 and F-4 have larger fault widths, extending to 21 and 29 km, respectively. In contrast, F-2 and F-3 show smaller fault width of |$\sim$|16 and |$\sim$|11 km, respectively. Figs 5(b), (e), (h) summarize the inferred fault model from the MAP solution projected to the ground surface. Regarding the inferred uniform slip on each segment, we find 2.5 m of slip on segment F-2, which is a bit higher than the 2.1 m estimated for segment F-1. Segments F-3 (1.3 m slip) and F-4 (0.9 m slip) show less average fault-segment displacement. Combining the fault dimensions and slip values, the estimated segment-specific seismic moment magnitudes |$M_{w}$| are 6.9, 7.3, 6.7 and 6.7 for F-1, F-2, F-3 and F-4, respectively, with a total |$M_{w}$| of 7.4. The corresponding seismic waveform fitting from the MAP solution of these fault parameters is shown in Fig. 6.

The marginal posterior PDF of correlations between segment-specific geometric fault parameters quantifies how reliable the data constrains the model parameters. The 1-D marginal posteriors show a multimodal PDF of width for fault segments F-1, F-2 and F-4 (Fig. 4). Such multimodal posteriors may arise from the nature of the SMC sampling algorithm, which utilizes a sequence of intermediate distributions and a tempering procedure to transition smoothly from a simple initial distribution to a more complex target distribution. This sampling algorithm is advantageous for handling high-dimensional and multimodal posterior distributions in case where a standard Monte Carlo technique may lead to slow convergence or poor estimates (e.g. Neal 1996; Del Moral et al. 2006, 2007; Nguyen et al. 2016). For segment F-2, this marginal PDF of correlations also shows a negative trade-off between fault dip and rake angle, and between fault width and slip, but it is still restricted to small parameter space. We observe little trade-off for other parameters, indicating that the data used for inferring fault geometry provides adequate constraints.

3.2 Rupture process of the 2021 |$M_{w}7.4$| Maduo earthquake

We adopt the inferred fault model geometry of the MAP solution to estimate the kinematic finite-fault rupture parameters, using the teleseismic data (P and SH waves), and the 3-D ground surface deformation data derived from the InSAR data. For a complete analysis, we carry out the FFI using first each data set separately before conducting a joint inversion. This approach aims to provide insights into the kinematic source model for this earthquake, assessing also the strength and limitations of each data set in constraining the rupture evolution.

The two slip models inferred from single data sets, teleseismic waveforms or geodetic data, show good agreement in their slip distributions (Fig. 7). Both models exhibit small amounts of slip around the hypocentre but large slip (⁠|$\gt 2$| m) at depths above 15 km, with a concentrated high slip patch on the eastern F-2 near the bifurcating fault. The peak slip from the geodetic data slip inversion is 4.9 m, slightly larger than the 4.7 m from the teleseismic data inversion. The standard deviation in sub-fault slip, obtained from a single data set (Supplementrary Figs S2a and b) indicates that the InSAR data effectively constrains the slip in shallower parts of the fault whereas teleseismic data are more sensitive to slip occurring at greater depths (Yong et al. 2014; Steinberg et al. 2020). The discrepancies between these two slip models may also arise because only fault slip is estimated using geodetic data, while additional temporal rupture parameters are inferred in the teleseismic inversion. López-Comino et al. (2015) report that teleseismic data only weakly constrain the solution, particularly in case of complex earthquake ruptures involving more heterogeneous slip patterns.

The joint-inversion rupture model exhibits analogous features to the slip models of the single data set inversions (Fig. 8). It reveals prevailing left-lateral slip, primarily distributed at shallow depths of less than 15 km, and multiple high- and low-slip patches, particularly in the eastward section of segment F-2. Naturally, areas of relatively low slip occur near asperities (i.e. areas of large slip), particularly at shallow depths (0–10 km) along the fault rupture. Based on two years of post-seismic surface displacements, the shallow low-slip areas appear to experience a relatively high amount of afterslip (Liu et al. 2024). Considering the faulting mechanism, the study by Qi et al. (2011) suggests that such peak-and-trough like slip patterns require spatial variations in rock strength, frictional properties or pre-earthquake loading on the fault. These variations then cause the stress levels in low-slip regions to remain below the fault’s failure threshold during an earthquake rupture. Essentially, these factors controlling slip behaviour may be attributed to spatial variations in fault friction, which explains the occurrence of multiple high-slip areas along the fault strike under the assumption of rate-and-state frictional behaviour (Kaneko et al. 2010; Qi et al. 2011).

The peak slip of 5.2 m is located close to the bifurcating fault in the southeast and extends to the surface, leading to more substantial surface displacements in that region. Field measurements further support this finding, indicating general coseismic surface offsets of approximately 1–2 m, with the maximum offset reaching |$\sim$|2.7 m (Ren et al. 2022). The standard deviation of fault slip in the joint inversion appears to be smaller at shallower depths than that from inversion using only seismic data (Supplementrary Fig. S2c), primarily due to the stronger constraints imposed by the geodetic data. In contrast, at greater depths the standard deviation of slip becomes smaller than that from the inversion using only geodetic data (Supplementrary Fig. S2), reflecting the influence of teleseismic data, which exhibits relatively higher sensitivity to fault slip occurring at increased depths. Additionally, Supplementrary Fig. S3 reveals the notable pattern that aftershocks tend to cluster in low-slip areas, while regions with larger slip appear to have fewer aftershocks. For further analysis in this paper, we refer to this joint-inversion rupture model as the preferred model.

The total seismic moment of the preferred slip model is 1.9 × |$10^{20}$| Nm (Fig. 8d), equivalent to a moment magnitude |$M_{w}=7.46$|⁠. It is slightly larger than the GCMT and USGS catalogue estimates (⁠|$M_{w}7.4$| and |$M_{w}7.34$|⁠, respectively). However, this is not unexpected since some post-seismic or aftershocks deformation may be included in the InSAR data collected 2-14 d after the main shock (He et al. 2021; Jin & Fialko 2021).

The comprehensive analysis of the kinematic finite-fault model and back projection images reveals an almost symmetrical bilateral rupture propagating at an average subshear velocity of 2.7 km s⁻¹ toward the northwest and southeast (Fig. 8a). The presence of multiple disconnected high-slip patches that correlate with rupture velocity variations interpreted from the FFI rupture-time contours and the clustered semblance high-frequency radiations from back projection indicate multiple stages of the rupture process, suggesting different dynamic failure episodes, as the rupture propagates away from the hypocentre to the east.

Fig. 8(a) illustrates that in the first 4–5 s, the rupture propagates radially away from the hypocentre. During this initial phase, the moment-rate function continuously increases and strong high-frequency energy radiation occurs (Figs 8b and c). Next, bilateral pulse-like rupture propagation occurs at nearly constant rupture speed for a few seconds. This corresponds to the first plateau in the moment rate function and minor seismic radiation seen in the back projection semblance. The rupture velocity starts to accelerate to the east at |$\sim$|7 s and decelerates at |$\sim$|12 s, interpreted from the narrowed and expanded interval of rupture time contour. In addition, we also observe large-slip areas that ruptured within this time range and up to 17 s. This variation in rupture velocity and the presence of large slip areas correlate with clustered high-frequency radiation observed in BP images. Considering the frequency range of the used P-wave recordings, we hypothesize that this combination contributes to stronger high-frequency seismic radiation within this time range in the eastward rupture.

Following a deceleration phase, just after 12 s, the rupture undergoes a rupture-speed increase, reaching an average speed of about 3 km s⁻¹ as it breaks the bifurcating faults. The clustering of radiation seen in the semblance plot suggests that the change in fault geometry likely results in the release of high-frequency seismic radiation followed by an abrupt stopping phase. In contrast, the rupture to the west exhibits less variation, reflected in a relatively constant interval of rupture time contours along strike. The seismic energy imaged by the back projection appears less clustered and more evenly distributed on the western F-2 segment. The fault geometry change at the splay fault in the west seems to be responsible for the clustered semblance around this location at |$\sim$|19–24 s.

To summarize, our kinematic finite-fault inversion and back projection analysis of the Maduo earthquake reveal an overall subshear rupture velocity and nearly symmetrical bilateral propagation, but substantial local rupture-speed variations. Additionally, the FFI-model portrays a continuous rupture propagation from F-2 to F-4, but a 3 s delay from F-2 to F-3. We attributed this delay to a change in fault orientation of F-3, making it less optimally oriented within the regional stress field and thereby requiring higher dynamic stress to be accrued before rupture is initiated on F-3.

4 DISCUSSION

In the following, we discuss specific aspects of the rupture process and fault geometry of the Maduo earthquake as inferred from our detailed source study. We also quantitatively compare our finite-fault model with previously published solutions and examine the implied stress changes due to this earthquake on nearby faults. Computing the change in |$\Delta {\rm CFS}$| due to the spatially heterogenous slip on the fault allows estimating if this earthquake brought any nearby fault closer to failure, which in turn may have implications for time-dependent seismic hazard in the region.

4.1 Fault geometry parameters

The near-linear surface-fault trace of the main rupture segments of the Maduo earthquake suggests a simple fault structure featuring a dominant straight fault, terminated by two splay faults where rupture stopped. However, detailed analysis shows it comprises several segments with different strike and dip directions. In terms of fault segmentation, several studies diverge in their modelling of the eastern fault section. Some studies conceptualized the eastern segment as a continuous splay fault (e.g. Yue et al. 2022; Wang et al. 2022) denoted as F-2 to F-3 in our study, without the southern branch. In contrast, others interpreted it as two bifurcating segments (e.g. Jin & Fialko 2021; Hong et al. 2022) (here, northern segment F-3 and southern segment F-4).

Detailed fault geometry investigations, incorporating analyses of aftershock distribution and near-field geodetic observations, reveal dip-angle variations along the fault. Our fault models show that segments F-1 and F-2 dip to the north, matching estimates based on near-field geodetic observations (e.g. Jin & Fialko 2021; Zhao et al. 2021). However, the aftershock-based geometry shows contrasting south- and north-dipping directions for segments F-1 and F-2, respectively (e.g. Wang et al. 2021; He et al. 2022). This discrepancy may be attributed to at least two factors, namely limited resolution of surface-rupture observations and inherent uncertainties in aftershock locations. The latter may reach several kilometres, depending on seismic-network configuration and station-event distances (Bondár et al. 2004). Previous studies also show inconsistent dip angles for segment F-3 to the east. Some models derived from geodetic near-field observations suggest a north dipping segment (e.g. Chen et al. 2021; Zhao et al. 2021; He et al. 2022) but aftershock data suggest a near-vertical structure (He et al. 2021). Different from those models, our inferred model reveals that the northeastern segment F-3 dips to the south at an angle of 71|$^\circ$|⁠, consistent with a previous study that links field observations and relocated aftershocks (Ren et al. 2022). Furthermore, from geodetic estimates, Wei et al. (2022) also infer a south dipping segment, but at a steeper angle (80|$^\circ$|⁠). For the southeastern segment, F-4, the aftershock data reveal a near-vertical dip to the north (e.g. He et al. 2021; Xu et al. 2023), which is consistent with our results.

Our fault geometry results are based on combined seismic and geodetic data and thus provide more constraints than previous studies. Furthermore, the inversion for fault parameters for all segments is conducted simultaneously to reduce likely trade-offs between fault parameters of different segments, particularly the well-known trade-off between fault dimension, fault depth and slip (Frietsch et al. 2019). Next, we summarize our sensitivity tests using different seismic data sets, which help to assess the effects of variations in seismic-station distributions on the solution.

4.2 Effect of different seismic data set distribution to fault geometric parameter inference

Using a similar inversion strategy, sampling parameters and prior distributions of fault parameters, we select another 20 stations, along with geodetic data, to perform a new joint inversion for fault geometry inference (Supplementary Fig. S1). We refer to this scenario as SC-2, while the previous one is SC-1. Compared to SC-1, the MAP solution of SC-2 (Supplementary Table S1) shows slightly shallower dip angles for segments F-1 (from |$84.2^{\circ }$| to |$83^{\circ }$|⁠) and F-2 (from |$85.2^{\circ }$| to |$83.8^{\circ }$|⁠), while segment F-3 remains unchanged and segment F-4 has a slightly steeper dip (from |$92.4^{\circ }$| to |$91.3^{\circ }$|⁠). These differences are small, however, suggesting that our estimated fault parameters are stable and robust. Our Bayesian inversion, using both the SC-1 and SC-2 station selections, yields generally similar fault geometries (Supplementary Fig. S4 and Table S1), and they both agree on the MAP solution for all segments. As shown in Fig. 4 and Supplementary Fig. S5 for SC-1 and SC-2, respectively, the marginal posterior PDFs consistently show multimodal distributions for the fault width of segments F-1, F-2 and F-4, but their MAP solutions are consistent. Overall, both the Bayesian inversion results and the marginal posterior PDF analysis show similar results for SC-1 and SC-2, implying that the estimated fault geometry is insensitive to the station selection as long as the full azimuthal range is covered.

The MAP solutions of SC-1 and SC-2 show comparable fits for both the 3-D surface displacement (Fig. 5 and Supplementary Fig. S6) and seismic waveforms (Fig. 6 and Supplementary Fig. S7), indicating a robust solution. However, we observe that some seismic stations located in near coastal areas, such as AU.MANU and AU.GIRL, exhibit a relatively poor fit in the low-frequency range (0.01–0.1 Hz). This discrepancy may stem from low-frequency noise in the observational data, potentially caused by oceanic and coastal waves (Hasselmann 1963; Longuet-Higgins 1950; Bromirski & Duennebier 2002), and the uniform-slip-per-segment source models of scenarios SC-1 and SC-2. Spatially varying finite-fault rupture parameters do improve the fit for these stations in the 0.05–0.5 Hz range, documented by an increase in VR |$\sim$|58 per cent and |$\sim$|15 per cent for AU.MANU and AU.GIRL, respectively. In addition, incorporating geodetic data in the joint inversion penalizes the fit to the seismic data, whilst different noise levels of each data set also affects the respective misfits. Furthermore, the varying non-zero time-shifts observed at each station suggest lateral variations in the velocity–density structure that are not captured in the 1-D Earth model used in our inversion.

4.3 Near-field surface displacements and rupture speed

Different from previous studies (e.g. Yue et al. 2022; Xu et al. 2023), we apply a Bayesian kinematic FFI in which we allow all rupture parameters to be spatially variable, sampled from appropriately specified a priori distributions. However, due to the high computational cost of the Bayesian sampling process, we use a relatively large sub-fault size of 5 km by 5 km compared to other studies (see Supplementary Table S3).

Furthermore, the surface-displacement misfits (Supplementary Fig. S8) reveal that the horizontal motions are better explained than the vertical ones, as seen by VR values of 97.4 per cent and 70 per cent for E–W and N–S displacements, respectively, compared to only 31.5 per cent for the vertical component. Notable residuals near the surface-rupture traces suggest that the assumption of purely elastic deformation in our modelling may be insufficient. Those large residuals may originate from inelastic off-fault deformation (e.g. Wollherr et al. 2019) and the simplified fault geometry adopted that obviously does not fully capture the detailed fine-scale geometry. High residuals near the fault trace may also stem from averaging the slip values over the sub-fault dimensions of 5 × 5 km|$^{2}$|⁠. Addressing these aspects is important for refining the accuracy and reliability of our surface displacement predictions. Additionally, the seismic waveform fits in Supplementary Fig. S9 show the match between model-generated waveforms and observed seismic data. Overall, approximately 60 per cent of the seismic data for both P-wave and SH-wave phases are well modelled, with VR exceeding 50 per cent. However, the amplitudes of the primary seismic pulse and early P- and SH-wave phases are frequently biased, leading to lower VR values (e.g. for stations IU.DAV, II.ERM and AU.GIRL).

Some previous studies on the Maduo earthquake suggested supershear rupture propagation towards the east based on back projection results (e.g. Yue et al. 2022; Zhang et al. 2022b). Our rupture models, in contrast, reveal an average subshear rupture speed of  2.7 km s⁻¹, propagating in a quasi-symmetric bilateral manner in the northwest and southeast directions, in agreement with Wei et al. (2022). Their study also suggests a direct triggering of both fault branches (F-3 and F-4 in our study) based on back projection results, however, our kinematic inversion indicates a |$\sim$|3 s delay for rupture on F-3. This delay may be attributed to the change in fault geometry, as F-2 and F-4 both have a steep dip angle to the north, while segment F-3 dips to the south with a lower dip angle.

4.4 Variability in slip models of the 2021 |$M_{w}7.4$| Maduo earthquake

To assess the variability of slip models published for the Maduo earthquake, we conduct a systematic quantitative analysis based on our work (Model-1a, b and c) and nine other published studies (Model-2 to Model-10, Fig. 9) (He et al. 2022; Jin & Fialko 2021; Zhao et al. 2021; Hong et al. 2022; Lyu et al. 2022; Wei et al. 2022; Yue et al. 2022; Wang et al. 2022; Xu et al. 2023). Note that the labelling of these models does not imply any ranking.

Robust metrics for quantitatively comparing rupture models are needed to identify models with common features. Comparing multiple slip models of the same earthquake, derived from various inversion approaches, data sets, assumptions on the Earth’s structure and model-space parametrization, reveals shared characteristics and may offer insights to better understand the detailed earthquake rupture process. In addition, the variability in inferred slip models for the same earthquake can be regarded as an epistemic uncertainty in a seismic-hazard context, which thus may be important for future simulation-based seismic hazard studies for the region.

To perform a quantitative comparison and capture the source–parameter variability in terms of location, amount and extent of slip on the fault, we employ the multidimensional scaling (MDS) approach of Razafindrakoto et al. (2015). Prior to the MDS analysis, we examine slip-model variations with respect to maximum and average fault slip. In addition, we also calculate the first-order centroid location of the slip models defined as |$C\left(c_{i},c_{j},c_{k}\right)$|⁠, where |$c_{i},c_{j},$| and |$c_{k}$| are described in eqs (3) (Razafindrakoto et al. 2015).

where |$U_{i}^{x}$|⁠, |$U_{j}^{x}$| and |$U_{k}^{x}$| are the slip projection at point x of the fault along the |$i-$|⁠, |$j-$| and |$k-$|axis, and n represents the total number of points. A summary of seismic moment and maximum slip for Models 1 to 10 is presented in Supplementary Table S3.

Examining the centroid locations for Models 1–10 (Fig. 10), we find that they are distributed within a depth range of 2 to 11 km. Notably, Model-4 (Jin & Fialko 2021) exhibits the shallowest slip centroid, whereas our models (Model-1a, b and c) and Model-2 (He et al. 2022) are located at greater depths (Fig. 10b). Intriguingly, the centroids of nine models form an elliptical cluster around a latitude of |$34.6^{\circ }$| and a longitude of |$98.49^{\circ }$|⁠, implying variability in slip magnitude both laterally and with depth.

In Fig. 11, average depth-dependent slip profiles, averaged along strike, for each individual fault segment are shown. For segments F-1 and F-2, we observe low variability between all models, which effectively show similar average slip patterns down to a depth of |$\sim$|10 km. However, for segments F-3 and F-4, the vertical slip profiles show substantial variations, which we attribute to different fault geometries and fault-model discretizations adopted in the source-inversion studies for these segments.

Moreover, the variation in average slip on segments F-3 and F-4 may be attributed to insufficient constraints during the inversion process, particularly considering that segment F-4 is considerably shorter than other segments and is positioned in close proximity to segment F-3. This spatial arrangement could result in a dominant influence of segment F-3 in fitting the observational data during the inversion. While some models indicate a shallow slip deficit, our models do not capture this feature, which aligns with our use of coarser and larger sub-fault size. Supplementary Figs S10 and S11 provide detailed comparisons of the average slip and shallowest slip across horizontal profiles, respectively.

Subsequently, we broaden the scope of the slip-model comparison with a more comprehensive assessment of the spatial distribution of slip by projecting and then re-griding different slip models onto the same overall fault geometry, specifically the Models 1 fault geometry in our study. For models lacking the southern segment F-4 in the east [Model-5 and Model-8 from Wang et al. (2022) and Yue et al. (2022), respectively], our comparison and analysis focus specifically on the variability of slip for the shared fault segments.

Razafindrakoto et al. (2015) highlights that the MDS point clouds depend on the chosen comparative metric, with the grey-scale exhibiting sensitivity to sub-fault geometry and small-scale slip variation. Aligning with the definition of large-slip and very-large-slip asperities by Mai et al. (2005), we convert the slip magnitudes into three colour intensity levels, as defined in Razafindrakoto et al. (2015): low-to-moderate slip (⁠|$U\le \frac{1}{3}U_{\rm max}$|⁠, U represents the slip amplitude of each patch), large-slip (⁠|$\frac{1}{3}U_{\rm max}\lt U\lt \frac{2}{3}U_{\rm max}$|⁠) and very-large-slip (⁠|$U\ge \frac{2}{3}U_{\rm max}$|⁠). We name such transformed slip information as slip intensity level.

Subsequently, we employ MDS to assess similarities and dissimilarities in the 2-D grey-scale metric (Wilson et al. 1997) of all slip models. Note that no reference model is defined and that each model is given an equal weight. Consequently, the MDS configuration facilitates a comparative assessment by capturing and scaling the dissimilarities among all slip models. To quantify the dissimilarity, we define four distinct categories based on the dissimilarity metric: (1) ‘Excellent’, |$\lt $|5 per cent; (2) ‘Good’ (5–20 per cent); (3) ‘Fair’ (20–40 per cent); (4) ‘Poor’ (⁠|$\gt $|40 per cent). These categories are described in a 2-D MDS configuration that visualizes the dissimilarities in a 2-D plot with two axes, dimension 1 and dimension 2, associated to the overall slip intensity level and to the spatial extent of this intensity, respectively (Razafindrakoto et al. 2015). Consequently, the two dimensions encapsulate distinct aspects of spatial slip heterogeneity, with no correlation between them.

In the combined MDS configuration, we explored two specific cases. The first case, with large common slip areas, encompasses the slip across the fault width of 25 km in the down-dip direction (bounded by dashed blue lines in Fig. 9). Conversely, the second case, involving small common areas, is restricted to slip within the upper 15 km (bounded by dashed red lines in Fig. 9). Fig. 12 displays the MDS configuration from all 12 slip models, assessed individually for each segment across the fault width of 15 and 25 km. Additionally, Fig. 13 illustrates the collective MDS configuration for F-1|$\cup$|F-2|$\cup$|F-3.

In Fig. 12, variations in dissimilarities among models across four fault segments are evident. For segment F-1, within a shallower common area (0–15 km), 10 out of 12 models exhibit ‘good-to-fair’ similarities, with dissimilarities particularly noted in our single data set seismic model and the one of Wei et al. (2022). However, for larger common areas (0–25 km), the model from He et al. (2022) stands out with the largest dissimilarities. In segment F-2, both small and common areas consistently show ‘good-to-fair’ similarity values, revealing two groups with different distributions of the slip intensity level. For segment F-3, expanding the area used for MDS reduces the dissimilarities, indicating improved model similarities. In segment F-4, enlarging the common area from 15 to 25 km increases the similarities among some models but accentuates dissimilarity for the model by He et al. (2022), potentially due to its large slip at greater depth, a characteristic also observed in segment F-1.

The 2-D MDS representations of slip variability for segments F-1, F-2 and F-3 account for more than 90 per cent of the overall dissimilarity, a value deemed acceptable according to Razafindrakoto et al. (2015) and Hair et al. (2010). However, the 2-D configuration for segment F-4 explains only 73 per cent, although it still exceeds the minimum acceptable level of 60 per cent (Hair et al. 2010). This suggests that visualizing the dissimilarities in slip intensity levels for segments F-1, F-2 and F-3 in a 2-D configuration is reliable, while for segment F-4, the 2-D configuration appears somewhat less valid but remains within an acceptable range.

Overall, the 2-D MDS of each fault segment of all twelve slip models (including our slip models obtained from a single data set inversion) for both large and small common areas reveal fault segments F-2 and F-3 having ‘good-to-fair’ similarities. This means these slip models capture common features of the rupture along these two fault segments. Most slip models reveal higher similarities in the slip-extension metric, while they are less similar in the slip intensity levels.

The 2-D MDS configuration shown in Fig. 13 illustrates the combined slip intensity level of segments F-1 through F-3, indicating minimal dissimilarities among all models and demonstrating ‘excellent to good’ similarities, except for the model of He et al. (2022). However, with an expansion to a larger common area, the similarities increase, causing the model by He et al. (2022) to move closer to other models.

4.5 Coulomb failure stress change (⁠|$\Delta {\rm CFS}$|⁠) from different slip models of the 2021 Maduo earthquake

Studies show that |$\Delta {\rm CFS}$| due to a ruptured fault may either delay or advance the next earthquake on a nearby fault according to the Coulomb failure criterion (e.g. Harris 1998; Stein 1999; King & Cocco 2001; Freed 2004; Toda & Stein 2020). Therefore, computing |$\Delta {\rm CFS}$| has become a common method to analyse and forecast aftershock occurrence following large earthquakes, although with large uncertainties mainly caused by non-unique slip inversion solutions and unknown receiver fault geometry (Woessner et al. 2006, 2012; Sharma et al. 2020). Its calculation employs coseismic slip models and is based on the expression (Harris 1998; King & Cocco 2001):

where |$\Delta \tau$| and |$\Delta \sigma _{n}$| denote the shear and normal stress changes, respectively, and |$\mu ^{\prime }$| is the effective frictional coefficient. Eq. (4) describes that increasing |$\Delta {\rm CFS}$| (positive value) may push a nearby fault closer to failure, while a negative |$\Delta {\rm CFS}$| may prevent or delay fault failure.

Calculating |$\Delta {\rm CFS}$| requires defining the geometry and faulting mechanism of the target faults onto which the stress perturbations are resolved. In this analysis, we calculate |$\Delta {\rm CFS}$| from our preferred slip model (Model-1c) and nine other slip models (Model-2 to Model-10) to assess the impact of slip model discrepancies, with and without considering fault model variability. Additionally, we propagate the slip-model uncertainties to |$\Delta {\rm CFS}$| by calculating the |$\Delta {\rm CFS}$| for the ensemble of slip models obtained from sampling the posterior PDF (Woessner et al. 2012; Dutta et al. 2018). We assume a uniform half-space with Poisson’s ratio of 0.25, Young modulus of 80 GPa and effective friction coefficient |$\mu ^{\prime }$| of 0.4. Considering the tectonic setting, strain field orientation and the fault map in Fig. 1, we find that large faults surrounding the ruptured segments of the Maduo earthquake have approximately similar orientations as the KLSK-JCF. Given that F-2 is the longest segment with an orientation nearly parallel to these surrounding large faults, we select it as the target fault (strike: 285|$^{\circ }$|⁠, dip: 85|$^{\circ }$|⁠). Additionally, based on the distribution of aftershocks and the location of the large slip, the |$\Delta {\rm CFS}$| calculation is performed at a depth of 7.5 km.

4.5.1 Variability in |$\Delta {\rm CFS}$| from different slip models for the Maduo earthquake

The computed |$\Delta {\rm CFS}$| from different slip models, illustrated in Fig. S12 in the supplementary materials, reveals an overall consistent |$\Delta {\rm CFS}$| map. To further examine the |$\Delta {\rm CFS}$| variability, we calculate the mean, standard deviation and coefficient of variation (CV) of |$\Delta {\rm CFS}$|⁠. The mean of the |$\Delta {\rm CFS}$| derived from the 10 slip models (excluding our slip models from single data set inversion) is depicted in Fig. 14(a). It highlights a positive |$\Delta {\rm CFS}$| along the extension of the KLSK-JCF fault to the west, the eastern part of the Tuosuo Lake segment of the East Kunlun Fault (EKF) and the Maqin-Maqu segment (MMS). This suggests that the Maduo earthquake exerted stress, bringing these segments closer to failure.

The standard deviation of |$\Delta {\rm CFS}$| (Fig. 14b) indicates substantial variability in proximity to the source fault, diminishing with increasing distance from the fault. The CV of |$\Delta {\rm CFS}$| (Fig. 14c) assesses the relative variability around the mean |$\Delta {\rm CFS}$| values. The CV map reveals elevated values near the fault and along the boundaries between positive and negative |$\Delta {\rm CFS}$| lobes, indicating regions where either the standard deviation is high or the mean |$\Delta {\rm CFS}$| is low. It is important to note that these |$\Delta {\rm CFS}$| variabilities stem from differences in both slip models and fault geometry.

Given that the location of positive and negative lobes of |$\Delta {\rm CFS}$| is sensitive to the fault strike (Dutta et al. 2018), we extend our analysis to calculate similar statistical measures for the projected slip models. This involves projecting the original slip models from previous studies onto our preferred fault geometry, thereby isolating the effect of different slip distributions. The mean |$\Delta {\rm CFS}$| for projected slip models (Fig. 15a) exhibits subtle differences with |$\Delta {\rm CFS}$| patterns in Fig. 14(a), indicating consistency in the |$\Delta {\rm CFS}$| calculated from both models. However, the standard deviation shows a slightly narrower area with large near-fault variability. This reduction in |$\Delta {\rm CFS}$| variability around the fault source aligns with expectations, as the slip models share a similar fault geometry. Therefore, the variability in |$\Delta {\rm CFS}$| primarily stem from differences in slip distribution, both in strike-slip and dip-slip magnitude. Notably, by removing the variability in fault geometry through the projection of slip to our fault geometry, the CV map in Fig. 15(c) illustrates a variability decrease of |$\Delta {\rm CFS}$|⁠, particularly at the boundaries of positive and negative lobes.

4.5.2 |$\Delta {\rm CFS}$| variability from ensemble solutions of joint seismic and geodetic finite-fault inversion

To assess how the uncertainty in our slip models obtained from Bayesian inference influences the |$\Delta {\rm CFS}$| calculation, we use 135 000 slip-model realizations as an ensemble of representative slip solutions. The mean of |$\Delta {\rm CFS}$| (Fig. 16a) reveals a pattern of |$\Delta {\rm CFS}$| lobes similar to that in Fig. 14. This suggests that the variability in our ensemble solutions does not alter the regions experiencing decreased or increased stress. However, the standard deviation and coefficient of variation (Figs 16b and c) are notably smaller compared to those in Fig. 14. This reduction in variability of |$\Delta {\rm CFS}$| indicates greater consistency among ensemble solutions, because variability of |$\Delta {\rm CFS}$| is primarily driven by slip variability.

4.6 Stress transferred by the |$M_{w}7.4$| Maduo earthquake to surrounding faults

Small variations in |$\Delta {\rm CFS}$| from our ensemble slip solutions allow us to use the MAP solution of coseismic slip (Fig. 8) to investigate the stress transferred by the Maduo main shock to the surrounding faults. For the receiver faults, the dip and rake angles are taken from the historical seismic events (Feng et al. 2022). These receiver faults are assumed to extend down to 30 km uniformly along the dip direction. We then compute the |$\Delta {\rm CFS}$| for a constant effective friction coefficient of 0.4, a commonly assumed |$\mu ^{\prime }$|-value for Coulomb stress transfer analysis. We note that the effective friction coefficient could vary between 0.0 and 0.8, resulting in 20 per cent–40 per cent variation in |$\Delta {\rm CFS}$| according to Parsons (2005).

The |$\Delta {\rm CFS}$| in Fig. 17 shows that the |$M_{w}7.4$| Maduo earthquake leads to an increase in |$\Delta {\rm CFS}$| by 0.1 MPa along the eastern part of the Tuosuo Lake segment on the EKF and along the fault extension of the ruptured segments to the west. A positive |$\Delta {\rm CFS}$| also occurs along the MMS and the eastern extension of the KLSK-JCF. This may advance the occurrence of the next large earthquake on the MMS of the EKF, where the seismic gap along the EKF has been hypothesized (e.g. Shan et al. 2015; Zhu et al. 2021a). This is consistent with previous reports (He et al. 2021, 2022; Guo et al. 2022; Feng et al. 2022). The Maqin-Maqu segment is characterized by the recurrence intervals of 500–700 and approximately 1000 yr, respectively (Li et al. 2011; Shan et al. 2015). Incorporating historical earthquakes from 1879 to 2008, Shan et al. (2015) argued that these earthquakes increased the stress on the segment, potentially leading to an advanced failure for |$\sim$|160 and |$\sim$|250 yr on the Maqin and Maqu segment, respectively. In addition, Li et al. (2011) reported that the most recent large earthquake on the Maqin segment occurred around 514–524  CE, while the Maqu segment last ruptured between |$\sim$|1055 and 1524 CE. Our |$\Delta {\rm CFS}$| calculation (Fig. 17) shows that the Maduo earthquake imposes an increasing stress on the MMS which may lead to accelerated failure on this segment.

5 CONCLUSION

The absence of uncertainty quantification in previous finite-fault inversion studies of the Maduo earthquake motivated us to apply a Bayesian approach to generate and explore an ensemble of solutions from which to then select the most probable ones for detailed analyses. Our study thus represents the first attempt to employ fully Bayesian inversion for inference of fault geometry and rupture parameters across multiple segments, utilizing combined seismic and geodetic data sets. Our primary goal is to advance the understanding of the rupture process of the Maduo earthquake by employing kinematic finite-fault inversion and back projection analyses, at the same time assessing also rupture-parameter uncertainties.

Our findings indicate that the Maduo earthquake ruptured bilaterally at a quasi-symmetric speed, with a relatively smooth forward rupture propagation to the west compared to the more intricate pattern observed towards the east. In contrast to previous studies that suggest an average supershear rupture speed to the east, our results suggest a relatively stable sub-shear rupture at approximately 2.7 km s⁻¹ in both the northwest and southeast directions. Along the bifurcated fault in the eastern region, the rupture appears to propagate continuously from segment F-2 to segment F-4, albeit with a slight delay of approximately 3 s from segment F-2 to segment F-3.

The lack of openly available local/regional seismic records prevented us from adding further constraints on the solutions. Because each data set has a different wavelength and frequency resolution on the physical process of the rupture, the resolution of our kinematic finite-fault models could potentially be improved if regional data were available, like in previous FFI-studies for moderate to large earthquakes (e.g. Delouis et al. 2002; Cirella et al. 2009; Yokota et al. 2011; Kobayashi et al. 2016).

Our thorough exploration of slip model variability highlights consistent findings. Factors such as fault geometry and inversion strategy play a crucial role in determining model discrepancies, even when using similar data sets. The intricacies of constraining the rupture process on adjacent fault branches with only minor spatial separation, coupled with limited near-field data, pose challenges in determining accurate fault parameters, particularly for the eastern bifurcated faults of the Maduo earthquake.

Examining |$\Delta {\rm CFS}$| change maps from different slip models underscores the substantial impact of fault geometry differences on |$\Delta {\rm CFS}$| variability, particularly at the boundaries of positive and negative stress change lobes. Meanwhile, this near-fault variability is more influenced by slip magnitude variability. Compared to the slip variability among models from various studies, our assessment of uncertainty in slip models from an ensemble of solutions in the final stage of our joint Bayesian inversion reveals reduced slip variance, notably minimizing the variability in |$\Delta {\rm CFS}$| near the fault. From |$\Delta {\rm CFS}$| analysis, one could speculate that the Maqin-Maqu segment of the EKF, with a potential seismic gap identified, may have been brought closer to failure as a result of the Maduo earthquake.

ACKNOWLEDGMENTS

We express our gratitude to the anonymous reviewers and the editor, Prof Ana M G Ferreira, for providing constructive comments that contribute to the improvement of the manuscript. The research presented herein received support from the King Abdullah University of Science and Technology (KAUST) under grant BAS/1/1339-01-01.

DATA AVAILABILITY

Real‐time aftershock locations were provided by Wang et al. (2021). Satellite data were downloaded from the Sentinel data hub (https://scihub.copernicus.eu). The moment-tensor solutions shown in red colour in Fig. 1 are taken from Global Centroid–Moment–Tensor (CMT) catalogue (https://www.globalcmt.org). 3-D surface deformation data derived from satellite data are available upon request to the authors. Teleseismic waveforms for back projection were downloaded from the Incorporated Research Institutions for Seismology (IRIS, https://www.iris.edu). Teleseismic waveforms for finite‐fault inversion were obtained from the IRIS, Geoforschungsnetz (GEOFON) at https://geofon.gfz-potsdam.de, and Observatories and Research Facilities for European Seismology (ORFEUS) at https://www.orfeus-eu.org data centres, respectively. Furthermore, this work utilized the following open-source libraries: numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), and pyrocko (https://pyrocko.org/) (Heimann et al. 2019). Plots were generated using matplotlib (Hunter 2007), the Generic Mapping Tools (GMT) (Wessel et al. 2013), and pygmt (https://www.pygmt.org/) (Uieda et al. 2021).

REFERENCES