Stress triggering and future seismic hazards implied by four large earthquakes in the Pamir from 2015 to 2023 revealed by Sentinel-1 radar interferometry

SUMMARY

The M_w 6.8 Murghob earthquake is the third earthquake in an M_w 6.4+ sequence occurring in the Pamir initiated by the 2015 Sarez M_w 7.2 earthquake. It is of great significance to investigate their interactions and to assess future seismic hazards in the region. In this paper, we use Sentinel-1 radar interferometric data to retrieve coseismic deformation, invert for the slip distributions of the four events, and then investigate their interactions. The cumulative Coulomb failure stress changes (ΔCFS) suggest that the 2023 Murghob earthquake was promoted by the three prior earthquakes in the sequence. Pre-stress from historical earthquakes is a key factor in explaining the triggering mechanism of the two 2016 M_w 6.4+ earthquakes. Stress loading and unloading effects on major faults in the region indicate that future attention should be paid in (1) the segment of the Sarez-Karakul fault north of the Kokuibel Valley, (2) the segment of the Sarez-Murghab thrust fault west of the Sarez-Karakul fault and (3) the east segments of the Pamir thrust fault system, all with a large positive ΔCFS.

Radar interferometry, Asia, Earthquake hazards, Earthquake interaction, forecasting and prediction, Earthquake source observations

1 INTRODUCTION

Large earthquakes induce regional stress changes that can alter regional seismicity activity, for instance, promoting or delaying the occurrence of an earthquake. Therefore, it is essential to investigate the coseismic stress changes induced by historical and recent earthquakes for the assessment of future seismic hazards. Estimates of Coulomb failure stress changes (ΔCFS) are the most common and direct way to relate past and future seismic events, and can be calculated to investigate earthquakes interaction and assess future seismic hazards. The slip distribution is the necessary input data used to calculate the ΔCFS resulting from large earthquakes and is directly related to the accuracy of the ΔCFS calculations, thus it is essential to acquire precise and reliable fault parameters and slip distributions (Weston et al. 2011; Zhu et al. 2021).

Since the Sarez M_w 7.2 earthquake occurred on 7 December 2015, seismicity has increased in the Pamir region of Central Asia. Evidence is the successive Sary-Tash M_w 6.4 earthquake on 26 June 2016, the Muji M_w 6.6 earthquake on 25 November 2016 and the Murghob M_w 6.8 earthquake on 23 February 2023 (Fig. 1). All the four earthquakes occurred in sparsely populated high-altitude areas (with an average altitude of 5000 m) with harsh topographic and climatic conditions. The sparse distribution of seismic station poses a challenge to investigate the seismogenic tectonics and local Coulomb stress evolution. With its notable features for monitoring surface deformation at all-times, in all weather conditions with a wide coverage and high spatial resolution, space-borne Interferometric Synthetic Aperture Radar (InSAR) has been proven by numerous studies to be a valuable means to study and understand seismic faults in remote areas (e.g. Li et al. 2011; Salvi et al. 2012; Elliott et al. 2016; Ha et al. 2022; Xiao et al. 2022; Liu et al. 2023; Yu et al. 2024).

Figure 1.

Tectonic setting of the study region. (a) Topographic map of the collision of the Indian Plate with the Eurasian Plate. The black dashed boxed area corresponds to the coverage of (b). (b) Tectonic setting of the Pamir Plateau. The red beachballs represent the four earthquakes in this study, and the grey beachballs represent regional historical strong earthquakes. The purple dashed box corresponds to (c). (c) Locally enlarged tectonic setting of four earthquakes in the Pamir Plateau from 2015 to 2023. Note that the faults in the map are from the GEM Global Active Faults Database (Styron & Pagani 2020; https://blogs.openquake.org/hazard/global-active-fault-viewer/), and locally corrected based on faults mapped by Schurr et al. (2014) and Metzger et al. (2017). AMFS: Aksu-Murghab fault system; DF: Darvaz fault; KF: Karakorum fault; KSF: Kongur Shen fault; MF: Muji Fault; PTS: Pamir thrust system; RE: Reshun-Hunza fault system; SKF: Sarez-Karakul fault system; SMTS: Sarez-Murghab thrust system and TTS: Tanymas thrust fault system.

A number of previous researches have been carried out on the three strong seismic events that occurred in the complex tectonic setting of the Pamir region from 2015 to 2016. Sangha et al. (2017), Metzger et al. (2017), Elliott et al. (2020) and Jin et al. (2022) investigated the fault geometry and slip distribution of the 2015 Sarez earthquake by combining InSAR, SAR and optical offsets, seismology, optical remote sensing and field surveys. In the fault modelling, Sangha et al. (2017) and Metzger et al. (2017) divided the fault into three segments for modelling, Jin et al. (2022) added a parallel fault near the Lake Karakul, and Elliott et al. (2020) subdivided it into nine segments for modelling. Although they used different fault segments, at least three fault segments were required to model the event. The epicentre of the 2016 Sary-Tash earthquake was located in the Trans Alai Range (at a maximum altitude of ∼7100 m); Although permanent snow coverage impeded coherent Sentinel-1 interferograms, He et al. (2018b) studied the seismogenic fault of the event using Sentinel-1 satellite radar observations. The 2016 Muji earthquake was extensively studied by Bie et al. (2018), Feng et al. (2017), He et al. (2018a) and Ma et al. (2018), who determined the source parameters and slip distribution of the event using InSAR and Global Navigation Satellite System (GNSS) observations. The source mechanism solutions for the 2023 Murghob earthquake vary considerably from different agencies (Table 1). Shi et al. (2023) used the Markov chain Monte Carlo method to globally search for fault geometry of the event based on Sentinel-1 satellite radar observations, and they concluded that the seismogenic fault is a steeply SE-dipping fault (∼73° dip and 28° strike; Table 1).

Table 1.

Source parameters of the 2023 Murghob earthquake. Note that the location of the model determined by Yong Zhang (pers. Comm.), Shi et al. (2023) and our study represents the location of the maximum slip on the corresponding fault plane for the distributed slip model. USGS: United States Geological Survey. GCMT: Global Centroid Moment Tensor. IPGP: Institute de Physique du Globe de Paris. GFZ: German Research Center for Geosciences (last accessed on 2 October 2023).

Source	Location			Fault plane			Magnitude
	Longitude (°)	Latitude (°)	Depth(km)	Strike (°)	Dip (°)	Rake (°)	(Mw)
USGS (2023)	73.23	38.055	9	203	57	−21	6.9
				305	72	−146
GCMT (2023)	73.22	38.16	14.5	120	87	−178	6.8
				30	88	−3
IPGP (2023)	73.21	38.073	8	123	71	−162	6.8
				27	73	−20
GFZ (2023)	73.29	38.06	10	122	73	−156	6.8
				25	67	−18
Yong Zhang	73.275	38.178	10.5	30	86	−4	6.9
Shi et al. (2023)	73.236	38.142	13.3	28.1	72.9	−1.9	6.83
This paper	73.182	38.105	9.7	33.1	78.0	−5.8	6.84

Some previous studies have been carried out on the interaction relationships between the first three events in the sequence. Jin et al. (2022) determined the fault geometry and slip distribution of the 2015 Sarez earthquake by combining InSAR from Sentinel-1 and Advanced Land Observing Satellite-2 (ALOS-2) and GNSS observations, then calculated the ΔCFS generated by the event on the fault planes of the two 2016 earthquakes. They found that the 2015 Sarez earthquake slightly prevented the nucleation of the 2016 Sary-Tash earthquake with a ΔCFS of −0.003 bar and promoted the nucleation of the 2016 Muji earthquake with a ΔCFS of +0.034 bar, and concluded that static or quasi-static Coulomb stress transfer hardly explains the triggering relationships between the Sarez event and the pair of M_w 6.4+ earthquakes. Bloch et al. (2022a) calculated the ΔCFS distribution generated by the 2015 Sarez earthquake on the fault of the two 2016 M_w 6.4+ earthquakes to be +0.04 and -0.19 bar, respectively, using the fault slip distribution model from Metzger et al. (2017) who determined the distributed slip model by combining Sentinel-1 InSAR, GNSS and post-seismic field observations. They concluded that the static ΔCFS effect was too small to explain the occurrence of the two 2016 earthquakes, and thought that fluids migrated through the damaged fault zones and might trigger the subsequent earthquakes. Probably due to the different source models constrained by different geodetic observations or stress calculation methods, Jin et al. (2022) and Bloch et al. (2022a) came to nearly opposite ΔCFS values, thus it is essential to determine unified and reliable fault models using unified geodetic observation constraints and inversion methods. In addition, whether the recent 2023 Murghob earthquake, occurring ∼40 km east of the 2015 Sarez earthquake, was triggered by the three prior strong events in the sequence and how they affect future seismicity in the region deserves in-depth investigation.

In this paper, we used Sentinel-1 radar observations to obtain the coseismic surface displacements of the four earthquakes in the Pamir region, determined the fault geometries and refined slip distributions based on rectangular dislocation models in an elastic half-space (Okada 1985). Our efforts refine the interaction of the four large earthquakes occurring between 2015 and 2023, and the future seismicity based on the Coulomb failure criterion.

2 TECTONIC SETTING

The Pamir Plateau is a result of the extrusion and collision of the Indian Plate into the Eurasian Plate and is located on far to the north of the collision (Robinson et al. 2004; Molnar & Tapponnier 1975; He et al. 2018a; Fig. 1a). After the collision at ∼50 Ma, it formed by at least 300 km of crustal shortening between the northern and southern Pamir Plateau (e.g. Fu et al. 2010; Burtman & Molnar 1993; Molnar & Tapponnier 1975). As the Pamir Plateau moved northwards, it disconnected the Tarim and Afghan-Tajik basins, and the last vestige of this disappearing crust is now the Alai Valley, located between the Northern Pamirs and the Southern Tien Shan (Fig. 1b; Schurr et al. 2014). The Trans Alai Range at a maximum altitude of ∼7100 m is the southern boundary of the Alai Valley, and is the front edge of the Pamir Thrust System (PTS), which accommodates most of the convergence of the Pamirs with the Tien Shan at about 13–15 mm yr⁻¹ (Zubovich et al. 2016, 2010), as quantified by the GNSS velocity field. It accounts for more than 1/3 of the total convergence of the Indian and Eurasian plates at this longitude (∼38 mm yr⁻¹), causing some of the highest strain rates worldwide at the PTS in the entire India–Asia collision zone (Reigber et al. 2001; Molnar & Stock 2009; DeMets et al. 2010; Zubovich et al. 2010). To the east of the Pamir Plateau, the Karakorum Fault extends northwards to the Muji Graben (Fig. 1c), which forms the eastern boundary of the Pamir Plateau. Within the Pamir Plateau, the Sarez-Karakul Fault (SKF), an active NNE-trending left-lateral slip fault system with a tensile component, divides the interior of the Pamir Plateau, linking the Karakul and Sarez lakes (Fig. 1c), accommodating the left-lateral shear and east–west extension of the Pamir Plateau (Schurr et al. 2014).

The Pamir region is dominated by reverse and strike-slip movements (Burtman & Molnar 1993; Molnar & Ghose 2000), specifically extrusion and dextral strike-slip shearing on the northern margin, and sinistral strike-slip and tensile extension in the interior of the Pamir Plateau (Schurr et al. 2014). According to the USGS earthquake catalogue, there have been five shallow-source strong events (M_w 6.4+) in the region since 2000 (Fig. 1b), including an M_w 6.7 earthquake in 2008 in Nura, Xinjiang, which occurred less than 50 km away from the two 2016 earthquakes (Sippl et al. 2014; Teshebaeva et al. 2014). The other four strong earthquakes since 2015 are the main focus of this study, being plotted with red beach balls in Fig. 1.

3 METHODOLOGY

3.1 InSAR data

We collected Sentinel-1A/B satellite radar images covering the four earthquakes (see Table 2 for details) and retrieved the coseismic deformation fields using differential InSAR by the open-source GMTSAR software (Sandwell et al. 2011; Xu et al. 2020). Then we used the Precise Orbit Ephemerides and Shuttle Radar Topography Mission 1-arcsec (∼30 m) digital elevation model to correct for orbital and topographic effects, respectively (Farr et al. 2007; Weiss et al. 2020). Interferograms are multilooked by a factor of 10 in range and 2 in azimuth and filtered by the adaptive phase filter (Goldstein & Werner 1998) with a filter strength of 0.5 to reduce phase noise. Given the fact that snow coverage destroys so much of the signal, we masked the poorly coherent areas with a threshold of 0.1 before phase unwrapping to mitigate the effects of unwrapping errors (Song et al. 2022; Zhang et al. 2023). Finally, we retrieved the line of sight (LOS) surface deformation from the interferometric fringe map using the SNAPHU unwrapping method (Chen & Zebker 2001).

Table 2.

Parameters of the selected interferograms and the numbers of downsampling points. Note that the abbreviations ‘A’ and ‘D’ represent ‘Ascending’ and ‘Descending’; the abbreviations ‘Perp.B’ represents perpendicular baseline; the column ‘Points’ represents the number of downsampled data sets used for the modelling.

Event	Track	Primary scene	Secondary scene	Perp.B	Points
		(YYYYMMDD)	(YYYYMMDD)	(m)
Sarez	A100	20151206	20151230	−18.434	5944
	D005	20151118	20151212	135.339	6924
Sary-Tash	A100	20160615	20160709	−11.599	847
	D107	20160616	20160710	−60.751	1552
Muji	A027	20161113	20161207	−95.518	2664
	D107	20161125	20161219	80.373	2297
Murghob	A100	20230203	20230323	22.294	1649
	D005	20230221	20230305	135.949	1608

To mitigate the effects of long-wave atmosphere error, all interferograms were corrected by the Generic Atmospheric Correction Online Service for InSAR (GACOS, Yu et al. 2018a, b, 2017). To reduce the spatial correlation and to enhance the computational efficiency, we downsampled the coseismic interferograms of the four earthquakes using a quadtree downsampling method (Jónsson et al. 2002), and the number of downsampling data sets that went into the modelling was shown in Table 2.

3.2 Modelling methodology

We used the PSOKINV software package that applies the Multipeak Particle Swarm Optimization (MPSO) algorithm (Feng & Li 2010) to invert the fault geometry and slip distribution in two steps (Fukahata & Wright 2008). In the first step, the fault geometries are determined in a least-squared sense assuming that the rectangular fault slips uniformly in an elastic half-space. In the inversion, the initial, minimum and maximum values of the fault parameters are set by considering the source mechanisms given by different earthquake agencies and the InSAR deformation field characteristics, and the optimal fault geometry is acquired by a non-linear search using the MPSO algorithm (see Table 3). In the second step, based on the fault parameters determined by non-linear inversion, the length (along strike) and width (along dip) of the rectangular fault plane are suitably extended and discretized into numbers of fault patches (e.g. 1 km × 1 km) to allow for spatial variation of slip on each patch, and then an iterative grid search method is used to estimate the spatially variable slip distribution.

Table 3.

Source geometric parameters of the four earthquakes in this study. Note that the Sarez S1, Sarez S2 and Sarez S3 represent northeastern segment, central segment and southwestern segment of the 2015 Sarez earthquake fault model, respectively. The location (Longitude, Latitude, Depth) is the fault top centre from the non-linear inversion. Dip represents the optimal dip angle determined using the L-curve.

Event	Longitude (°)	Latitude (°)	Depth (km)	Strike (°)	Dip (°)	Rake (°)	Length (km)	Width (km)	Magnitude (Mw)
Sarez S1	73.141	38.571	1.830	210.6 ± 0.3	87	16.9 ± 0.7	27.1 ± 0.4	17.2 ± 1.4	7.22
Sarez S2	73.016	38.437	0.752	232.8 ± 0.4	64	−8.4 ± 0.3	14.3 ± 0.3	10.9 ± 0.3
Sarez S3	72.885	38.303	0.547	210.6 ± 0.2	82	−8.1 ± 0.2	26.8 ± 0.4	10.7 ± 0.2
Sary-Tash	73.592	39.421	7.937	274.7 ± 2.7	76	125.3 ± 5.3	15.5 ± 0.8	4.0 ± 1.0	6.36
Muji	74.248	39.211	2.591	105.5 ± 0.3	73	−172.0 ± 0.9	35.7 ± 1.2	8.7 ± 2.7	6.59
Murghob	73.208	38.167	2.136	33.1 ± 0.9	78	−5.0 ± 0.6	25.0 ± 0.1	10.1 ± 0.5	6.84

Following Feng et al. (2013), we used logarithmic trade-off curves plotting the model roughness and residuals to update the optimal dip angle and determine the smoothness factor of the distributed slip model. Additionally, a Monte Carlo estimation of correlated noise is used to estimate errors in fault geometry parameters determined by the MPSO non-linear search (Funning et al. 2005; Wright et al. 2003). The observations are added to Gaussian noise to generate 100 perturbed data sets, and each set is modelled separately to calculate the uncertainties in the fault parameters determined by the inversion (see Table 3). Meanwhile, we also used the Monte Carlo method (Biggs et al. 2006; Parsons et al. 2006) to estimate the uncertainty of our slip distribution model. The same parameter is used to invert each of these 100 perturbed datasets to derive 100 slip distribution models. Then, we calculated the standard deviation of all these slip distribution models to quantify the slip uncertainties at each fault patch.

4 SOURCE MODELLING AND SLIP DISTRIBUTION

4.1 The 2015 Sarez M_w 7.2 earthquake

The initial event in the sequence of four earthquakes occurred near Sarez Lake in the Pamir Plateau, where there was also a strong earthquake (M 7.3) in 1911 (Fig. 1; Kulikova et al. 2016; Elliott et al. 2020). It caused a ∼37 km discontinuous surface rupture that extends northward from the Sarez Lake to the Kokuibel River valley (Fig. 1c; Metzger et al. 2017; Elliott et al. 2020), roughly along the left-lateral SKF within the plateau. Combined with the USGS focal mechanism (USGS 2015) and InSAR coseismic deformation (Fig. 2), the event can be initially identified as a left-lateral slip event occurring on the SKF. To obtain unified and reliable fault slip models, we redefined the slip distribution model of this earthquake. In the modelling, we firstly used single and double faults to conduct modelling trials and found that all of them had difficulty in retrieving the coseismic deformation field. Then, we divided the fault into three subfaults based on the Sentinel-1 SAR offsets, and used the MPSO algorithm to non-linearly search for the best-fitting fault geometry.

Figure 2.

Coseismic modelling of the 2015 Sarez earthquake. (a) Observed coseismic interferograms, modelled interferograms and the residuals from Track 100A and 005D. Note that the red boxes indicate the surface projections of the fault planes, and the solid blue lines represent the intersection between the fault planes and the surface. The labels ‘los’ and ‘Az’ are abbreviations for range and azimuth directions. (b) Coseismic slip distribution in 3-D with a topographic map. (c) Fault plane slip model of the 2015 Sarez earthquake and its uncertainty estimated using Monte Carlo method (Parsons et al. 2006).

The modelling results suggest that at least three subfaults with different fault geometries (defined as the northeastern, central and southwestern segments, from north to south in the map view) are required (Fig. 2b). The non-linear search results indicate a near-vertical northeastern segment and westward dipping faults in the central and southwestern segments, with dips of 69.2 ± 0.6° and 80.0 ± 0.4°, respectively (Table 3). The northeastern and southwestern segments strike almost identically at 210.6 ± 0.3° and the central segment at 232.8 ± 0.4°. Furthermore, the rake angles (−8.4 ± 0.3° and −8.1 ± 0.2°, respectively) of the central and southwestern segments of the main slip indicate that the fault nature is left-lateral strike-slip with a tensile component, further confirming that the seismogenic fault is the left-lateral strike-slip SKF. In addition, we also searched for optimal dip angles and smoothness factors for the three segments, indicating that the optimal models dip 87° southeast, 64° northwest and 82° northwest and the smoothing factors are the same for all segments (Fig. S1). The linear inversion results show that the majority of slip localized mainly on the central and southwestern segments, with a maximum sliding of 5.5 m occurring at a depth of 8.8 km in the southwestern segment (Figs 2b and c); The maximum slip uncertainty is ∼0.18 m, which is less than 4 per cent of the maximum slip (5.5 m), suggesting that our slip model is reliable. The estimated moment release is 7.46 × 10¹⁹ N·m, corresponding to a moment magnitude of M_w 7.22. Although we have used a relatively simplified fault model, we have also modelled almost similar coseismic displacement (Fig. 2a), with an overall explained ratio [defined as (1−abs(residuals)/observations] of 96.2 per cent.

The dips and orientations of the three fault segments of Sangha et al. (2017) and Metzger et al. (2017) are slightly different (Table 4). Compared to Sangha et al. (2017), our results are more consistent with Metzger et al. (2017) in the fault orientations, and also agree with Elliott et al. (2020) with respect to close fault segmentation. In addition, compared to the steeper dips of Sangha et al. (2017) and Metzger et al. (2017), our model is consistent with Elliott et al. (2020), especially for the central segment. The seismic moment released by the fault slip is slightly smaller than the results from Sangha et al. (2017), but in closer agreement with Metzger et al. (2017) and the GCMT catalog (GCMT 2015; Table 4).

Table 4.

Parameters of different fault segments of the 2015 Sarez earthquakes.Note that the Sarez S1, Sarez S2 and Sarez S3 represent northeastern segment, central segment and southwestern segment of the 2015 Sarez earthquake fault model, respectively.

Model	Fault segment	Length (km)	Strike (°)	Dip (°)	Orientation	Moment (× 1019 N·m)
Metzger et al. (2017)	Sarez S1	23.8	205	89.3	SE	6.2
	Sarez S2	18.5	227	81.8	NW
	Sarez S3	23.5	217	87.7	NW
Sangha et al. (2017)	Sarez S1	26.55	210.1	83	NW	13.7
	Sarez S2	12.56	234.6	80	NW
	Sarez S3	41.22	215.4	89	SE
This study	Sarez S1	27.1	210.6	87	SE	7.5
	Sarez S2	14.3	232.8	64	NW
	Sarez S3	26.8	210.6	82	NW
GCMT (2015)	-	-	-	-	-	7.8

4.2 The 2016 Sary-Tash M_w 6.4 earthquake

The second event in the sequence of four earthquakes occurred near the Trans-Alai Range on the border of the Pamir Plateau and the Southern Tien-Shan extrusion collision (Figs 1b and c). According to the USGS focus mechanism (USGS 2016a), the seismogenic fault of the 2016 Sary-Tash event may be a shallowly south-dipping reverse fault (node I: 55°/35°/79°, corresponding to strike, dip and rake) or a steeply north-dipping reverse fault (node II: 248°/56°/97°). Since it is difficult to determine the location and orientation of faults from InSAR displacement maps (Fig. 3a), we consider two fault scenarios, a steeply north-dipping fault model (FM1) or a shallowly south-dipping fault model (FM2), to model this earthquake.

Figure 3.

Coseismic modelling of the 2016 Sary-Tash earthquake. (a) Observed coseismic interferograms, modelled interferograms and the residuals from Track 100A and 107D. Other features as in Fig. 2. (b) Coseismic slip model of the 2016 Sary-Tash earthquake and its uncertainty. (c1) Spatial distribution of aftershocks of the 2016 Sary-Tash earthquake. The purple dots represent aftershocks, from Bloch et al. (2022b), the blue and red boxes are the ground projections of the FM1 and FM2 fault planes respectively, and the cyan and yellow lines represent the intersection between the fault planes and the surface of FM1 and FM2. (c2) Aftershock profiles with the projected lines of the fault planes of FM1 (blue line) and FM2 (red line).

The results show that the surface displacements simulated by both models are in high agreement with the observed coseismic displacement. The FM1 has a strike of 85.8 ± 5.4°, a dip of 25 ± 0.4°and a rake angle of 72.5 ± 2.4°, with an overall explained ratio of 90.7 per cent and an RMS of 1.8 cm, and the FM2 has a strike of 274.7 ± 2.7°, a dip of 72.7 ± 2.8° and a rake angle of 125.3 ± 5.3°, with an overall explained ratio of 92.6 per cent and an RMS of 1.6 cm. Although the fit residuals of the FM2 are smaller, the seismogenic fault geometry cannot be determined using InSAR observations alone.

To further analyse which model is more reasonable, we used the aftershock relocation results provided by Bloch et al. (2022b) to draw aftershock profiles (Fig. 3c2) and projected the inversion-derived fault planes of FM1 and FM2 onto the aftershock profiles, as shown by the blue and red lines in Fig. 3(c2). As can be seen from the aftershock profiles, the FM1 cuts across the area of concentrated aftershocks and concentrates more aftershock activities in its footwall. In general, aftershocks are concentrated in the vicinity of strong coseismic slip areas, and the FM1 is difficult to explain reasonably for aftershocks at depth. Therefore, we prefer the north-dipping FM2 as the main fault plane for the event, because under this model, the aftershocks are mainly concentrated in the hanging wall of the fault, and the aftershocks are concentrated at depths of 6–16 km, which is consistent with the depth range of the main slip obtained from the inversion of the FM2 model. Similarly, He et al. (2018b) also found that both a steeply north-dipping fault and shallowly south-dipping fault can explain the InSAR surface displacement characteristics. Considering a better fit to the InSAR data, focal mechanism solutions and regional seismogenic thickness, they preferred the north-dipping fault geometry as the source for the event, which agrees with our solution based on the aftershock profiles.

Based on the above analysis, the 2016 Sary-Tash M_w 6.4 earthquake had a north-dipping fault with a dip of 72.7 ± 2.8°, a near east–west strike of 274.7 ± 2.7° and a rake angle of 125.3 ± 5.3°, indicating that this earthquake was a reverse fault event with a right-lateral strike-slip component. The fault geometry is generally consistent with (He et al. 2018b, ∼67° dip, 266° strike, and 126° rake). The optimal dip angle determined using the L-curve test is 76° (Fig. S2). The slip distribution shows that the main slip occurred at depths ranging from 6 to 16 km, with the maximum slip of 1.6 m occurring at a depth of 10.8 km (Fig. 3b).The maximum slip uncertainty is ∼0.17 m (Fig. 3b), which is about 10 per cent of the maximum slip. The moment released by the distributed slip was 3.85 × 10¹⁸ N·m, corresponding to a moment magnitude of Mw 6.36, which is slightly lower than (He et al. 2018b, 6.61 × 10¹⁸ N·m), but closer to the USGS seismic moments (USGS 2016a; 4.55 × 10¹⁸ N·m).

4.3 The 2016 Muji M_w 6.6 earthquake

The third event in the sequence of four earthquakes occurred in the Muji graben on the northeastern boundary of the Pamir Plateau, ∼60 km away from the 2016 Sary-Tash earthquake (Figs 1b and c). Field survey materials demonstrate limited localized surface ruptures near the Muji Fault (two ∼10-m-long ruptures; Chen et al. 2016), suggesting that this earthquake was a blind fault event that occurred on the Muji Fault. According to the focal mechanism from USGS (USGS 2016b), the earthquake might occur on a near north–south trending left-lateral strike-slip fault plane (node I: 199°/84°/14° for strike/dip/rake) or a near east–west trending right-lateral strike-slip fault plane (node II: 107°/76°/174°). Combining the characteristics of the InSAR coseismic deformation fields of the ascending and descending tracks and the tectonic background of the regional faults, the parameters of node II were chosen as a reference and constrained the fault strike to be between 80° and 130°, the dip to be between 60° and 90° and the rake to be between 150° and 200° for the non-linear search.

The model parameters converged to a strike of 105.5°, a dip of 73.5°, and a rake of −172°, indicating that the earthquake was a dextral slip event with a normal fault component and its seismogenic fault is south-dipping. The optimal dip angle determined using the L-curve test is 73° (Fig. S3). The distributed slip pattern indicates that the seismogenic fault plane had two asperities (Figs 4 b and c), with slip occurring mainly in the 2–15 km depth range; The maximum slip occurred at a depth of 4.8 km, on the east asperity, with a slip of 1.3 m; The maximum slip uncertainty is ∼0.15 m, which is located in the deeper western asperity (Fig. 4c). The seismic moment released from the fault plane was 8.7 × 10¹⁸ N·m, corresponding to a moment magnitude of M_w 6.59. The distributed slip model reached an overall explained ratio of 92 per cent with an RMS of 1.5 cm.

Figure 4.

Coseismic modelling of the 2016 Muji earthquake. (a) Observed coseismic interferograms, modelled interferograms and the residuals from Track 027A and 107D. Other features as in Fig. 2. (b) Coseismic slip distribution in 3-D with a topographic map. (c) Fault plane slip model of the 2016 Muji earthquake and its uncertainty.

The fault geometry of our model (Table 3) is almost consistent with Bie et al. (2018, 106.4°/70°/−176° for strike/dip/rake) and Feng et al. (2017, 105.5°/80°/−161°) who also inverted the fault model by using the Sentinel-1 InSAR observations constraints, but differs slightly from the fault model of He et al. (2018a) derived from the combined InSAR and GNSS inversion. The seismic moment released by our fault model is also consistent with Bie et al. (2018, 8.27 × 10¹⁸ N·m) and Feng et al. (2017, 9.87 × 10¹⁸ N·m), slightly smaller than that of He et al. (2018a, 1.35 × 10¹⁹ N·m).

4.4 The 2023 Murghob M_w 6.8 earthquake

The latest event in the sequence of four earthquakes occurred ∼67 km west of Murghob, Tajikistan, which is less than 45 km away from the 2015 M_w 7.2 Sarez left-lateral slip seismic event (Figs 1b and c). Although different agencies have given different focal mechanisms for the event, most of them (e.g. GCMT, GFZ and IPGP) agreed that the earthquake either ruptured on a SEE-trending right-lateral strike-slip fault or ruptured on an NNE-trending left-lateral strike-slip fault. Similar to the 2016 Sary-Tash earthquake, the fault geometry of this event could not be uniquely determined from the InSAR surface displacement characteristics (Fig. 5). Therefore, we modelled the two fault scenarios, a SEE-trending right-lateral strike-slip fault plane (FP1) and an NNE-trending left-lateral strike-slip fault plane (FP2). Considering that the two nodes of the updated USGS focal mechanism are controversial with other agencies only in terms of fault orientation, we constrain the range of fault dips to [30°, 150°] in our non-linear search for FP1 and FP2 so as to include the case of the two nodes of the USGS.

Figure 5.

(a) Observed coseismic interferograms, modelled interferograms and the residuals from Track 100A and 005D for the model FP1. (b) As (a), but for the model FP2. Other features as in Fig. 2.

The modelling results show that the FP1 has a strike of 127 ± 0.2°, a dip of 85.0 ± 0.2° and a rake angle of −178.0 ± 0.1°, with an overall explained ratio of 88.9 per cent and an RMS of 2.7 cm (Fig. 5a), and the FP2 has a strike of 33.1 ± 0.9°, a dip of 78.0 ± 2.2° and a rake angle of −5.8 ± 0.6°, with an overall explained ratio of 93.4 per cent and an RMS of 2.1 cm (Fig. 5b). It is obvious that the FP2 model exhibits smaller fitting residuals for both the ascending and descending InSAR interferograms (Figs 5a and b). Considering that the interior of the Pamir Plateau is dominated by left-lateral tensile movements (Bloch et al. 2022a; Schurr et al. 2014), we prefer FP2 as the rupture model for the 2023 Murghob earthquake.

Based on the above analysis, the 2023 Murghob earthquake had a SE-dipping fault with a dip of 78.0 ± 2.2°, a near NNE strike of 33.1 ± 0.9° and a rake angle of −5.8 ± 0.6°, indicating that this earthquake was a left-lateral strike-slip fault event with a tensile component, which is consistent with Shi et al. (2023) and the slip model (Fig. S4) of Yong Zhang (personal communication). The optimal dip angle determined using the L-curve test is 78° (Fig. S5). The distributed slip results show that the slip occurred mainly in the depth range of 2–15 km, with the maximum slip of 4.0 m occurring at a depth of 9.7 km (Figs 6a and b). The maximum slip is slightly larger than that of Shi et al. (2023), and its depth is slightly shallower than that of Shi et al. (2023), which may be attributed to different downsampling strategies and model smoothing parameters. The maximum slip uncertainty is ∼0.28 m (Fig. 6b), which is less than 10 per cent of the maximum slip (4 m), suggesting that our slip model is reliable. The seismic moment released from the fault plane was 2.0 × 10¹⁹ N·m, corresponding to a moment magnitude of M_w 6.84, which is consistent with Shi et al. (2023) and slightly lower than that of Yong Zhang determined by seismic wave inversion (Table 1).

Figure 6.

(a) Coseismic slip distribution in 3-D with a topographic map. (b) Fault plane slip model of the 2016 Muji earthquake and its uncertainty.

5 DISCUSSIONS ON LOCAL STRESS EVOLUTION

As one of the regions with the highest strain rates globally, the Pamir region is of great interest for the seismic community (e.g. Bloch et al. 2022a; Jin et al. 2022; Xiong et al. 2019). According to the Coulomb failure criterion, seismic activities are promoted when ΔCFS is positive and inhibited when ΔCFS is negative, which has been observed by numerous past events (e.g. Harris 1998; King et al. 1994; Stein 1999; Yan et al. 2018; Zhu & Wen 2010). In this study, we calculated static ΔCFS based on the unified fault slip models of these four seismic events, using the Coulomb 3.4 software (Lin & Stein 2004). The effective friction coefficient was set to 0.4 as suggested by previous studies (Jin et al. 2022; Xiong et al. 2019).

5.1 Stress evolution induced by the recent four major earthquakes

Considering the centroid depths of the four events, we calculated uniformly the ΔCFS at 10 km to investigate their interactions, and the resultant ΔCFS evolution is shown in Fig. 7. First, using the Sarez event as the source and the modelled fault geometries of the three subsequent events as the receivers, we calculated the static stress perturbations from the Sarez event on several subsequent events (Figs 7a–c). The resultant ΔCFS indicates the fault plane of the Sary-Tash earthquake is located entirely in the positive ΔCFS region, with a maximum ΔCFS of +0.15 bar, suggesting the Sarez event might promote the rupture of the Sary-Tash event. However, the fault plane of the Muji earthquake is located entirely in the ΔCFS shadow zone, indicating the Sarez event might delay the Muji earthquake. The majority of the fault plane of the Murghob earthquake is located in the stress-loading region of the Sarez earthquake, with a maximum ΔCFS of +2.73 bar, revealing a strong triggering effect of the Sarez earthquake on the Murghob earthquake. Secondly, we also calculated the static ΔCFS for the Sarez and Sary-Tash events on the two earthquakes subsequently (Figs 7d and e), and found that the Muji earthquake was still located in the ΔCFS shadow zone, reaching a similar conclusion with using only the Sarez event as the source. Finally, we calculated the ΔCFS generated cumulatively by the first three events (i.e. Sarez, Sary-Tash and Muji) on the Murghob event (Fig. 7f). Not surprisingly, the cumulative stress effects of the first three events still promoted the Murghob earthquake.

Figure 7.

Static ΔCFS evolution of the four earthquakes at the depth of 10 km.

The calculated ΔCFS of the two 2016 M_w 6.4+ earthquakes is consistent in sign with Bloch et al. (2022a), but differ from Jin et al. (2022), which may be attributed to the different source models derived from different geodetic constraints or calculation methods. According to Shan et al. (2013) and Shao et al. (2016), since steady tectonic loading can cause a stress change of ∼0.1 bar in a short period of time (Stein et al. 1997), the ΔCFS is considered insignificant when it is less than +0.1 bar. Thus, the calculated ΔCFS of the two 2016 M_w 6.4+ earthquakes in this study are too small (∼0.1 bar) to seem to be sufficient to explain the triggering mechanism of the two 2016 events, due to the tremendous decay of the ΔCFS over distances of 150 km. Similarly, Bloch et al. (2022a) and Jin et al. (2022) have come to the same conclusion as us.

Were the 2016 M_w 6.4+ events triggered by the 2015 Sarez earthquake, and if so, what was the triggering mechanism? Bloch et al. (2022a) argues that fluid migration in the damage fault zone induced by preceding earthquake may have played an important role in triggering the two subsequent events, but Jin et al. (2022) suggests that delayed dynamic triggering may be another way to explaining the causal relations (if any) between the 2015 Sarez earthquake and the subsequent M_w 6.4+ events. Earthquake triggering is a complex process, according to Mildon et al. (2019) and Liu et al. (2023), Coulomb pre-stress calculated for faults is an ignored yet vital factor for earthquake triggering.

To address the above-mentioned question, we calculated the pre-stress generated by historical strong earthquakes prior to the 2015-2023 earthquake sequence. We collected the fault models (Table S1) of the 2008 Nura earthquake modelled by Teshebaeva et al. (2014) and a simplified rectangular model of the 1974 Markansu earthquake [from Xiong et al. (2019) and Jackson (1979)] calculated from empirical equations (Wells & Coppersmith 1994). They were used as source faults along with the fault model of the 2015 Sarez earthquake obtained in this study and the fault geometry of the two 2016 earthquakes as receiver faults to calculate the ΔCFS on the fault of the two 2016 M_w 6.4+ earthquakes. The resulting ΔCFS is shown in Fig. 8. It is evident that when the coseismic slips of the 2008 Nura and 1974 Markansu earthquakes are taken into account, half of the 2016 Sary-Tash earthquake fault is located in the positive ΔCFS region, with a maximum ΔCFS of +6.32 bar. Also, it can be found that the seismogenic fault of the Muji earthquake is located entirely in the region of the positive ΔCFS generated by historical earthquakes, with a maximum ΔCFS of +2.52 bar (Fig. 8). If we consider the effects of regional historical earthquakes, the static stress transfer can be used to reasonably explain the triggering mechanism of the two 2016 M_w 6.4+ earthquakes. Therefore, when investigating earthquake triggering mechanisms, pre-stresses generated by historical earthquakes may play a key role and are not to be ignored.

Figure 8.

Static ΔCFS at the depth of 10 km generated by preceding earthquakes on the fault of the two 2016 M_w 6.4+ earthquakes.

5.2 Future seismic hazards

Post-seismic deformation plays an important role in calculating ΔCFS as well. In particular, on a long timescale, whether or not to take into account the post-seismic effects may lead to different conclusion (e.g. Chéry et al. 2001; Pollitz et al. 2003; Xiong et al. 2019). However, during a short time period (e.g. 3–5 yr), when afterslip and poroelastic rebound mechanisms dominate surface deformation, the ΔCFS from afterslip or poroelastic rebound is small and localized compared to coseismic ΔCFS, and shall not significantly change the amount of total ΔCFS (e.g. Huang et al. 2014; Wang et al. 2017; Liu et al. 2023). In addition, due to the poor observation conditions of the Pamir Plateau, it is difficult to detect small signals in the post-seismic period (Jin et al. 2022). Similarly, interseismic tectonic signals with smaller deformations are more difficult to obtain in this region. Therefore, given that all these events in this study occurred in recent years and the poor observation conditions for inter- and post-seismic deformation in the region, we only calculated coseismic ΔCFS in this study.

We collected the fault geometries (Table 5) from historical literature and GEM database (Styron & Pagani 2020) for some major faults in the region. Then, we calculated the ΔCFS distribution generated by the recent four earthquakes and the two historical M_w 6.6+ earthquakes mentioned above, on these major fault planes in the Pamir region. Note that the receiver faults were discretized into 5 km × 2 km patches extending downdip to 30 km from the surface.

Table 5.

Fault geometries of the main faults in the Pamir region.

Fault plane	Abbr.	Orientation	Dip (°)	Rake (°)	Source
Aksu-Murghab fault system	AMFS	-	90	180	Schurr et al. (2014)
Darvaz fault	DF	E	80	−30	GEM
Kongur Shan fault	KSF	W	75	−90	Brunel et al. (1994)
Muji fault	MF	S	75	−172	This paper
Pamir thrust system	PTS	S	45	90	Schurr et al. (2014)
Reshun-Hunza fault system	RE	N	50	90	GEM
Sarez-Karakul fault system	SKF	W	85	−5	This paper
Sarez-Murghab thrust system	SMTS	S	45	90	Schurr et al. (2014)
Tanymas thrust fault system	TTS	N	45	90	Schurr et al. (2014)

The resulting ΔCFS distribution is presented in Fig. 9. In the interior of the Pamir Plateau, one apparent pattern is that the unloaded fault segments are mostly distributed near the intersection area of the NNE-trending SKF fault and the E-trending TTS fault and SMTS fault, which means that the occurrence of historical earthquakes unloaded the stresses and somewhat reduced the future earthquake risk in the intersection region. The loaded fault segments are the southern segment of the SKF fault (area marked by grey dashed lines circle in Fig. 9), the segment of the SKF fault north of the Kokuibel Valley (area marked by blue dashed lines circle in Fig. 9), the segments of the SMTS fault east and west of the SKF fault, the TTS fault east of the SKF fault, and the eastern segment of the AMFS fault. In particular, the segment of the SKF fault north of the Kokuibel Valley (Fig. 1c), and the segment of the SMTS fault west of the SKF fault, with maximum stress loadings of more than +5 bar, are of concern for seismic hazards. On the northern edge of the Pamir Plateau, stress-loaded fault segments are concentrated on the east segments of the PTS fault (area in the green dashed circle in Fig. 9) and localized section of the Muji Fault (Fig. 9). In addition, the occurrence of historical earthquakes also unloaded the stresses of the northern part of the Kongur Shan fault in the northeastern Pamir Plateau and the left-lateral Darvaz fault in the western Pamir Plateau.

Figure 9.

The ΔCFS distribution generated by the recent four earthquakes and the two historical M_w 6.6+ earthquakes on major fault planes in the Pamir region. The topographic map corresponding to the ΔCFS distribution is plotted in the upper left-hand corner, and the locally enlarged map of the ΔCFS distribution is plotted in the upper right-hand corner.

The effect of post-seismic deformation on ΔCFS calculations is usually small in the short-term post-seismic period. To verify this, we used the post-seismic kinematic slip model for the 2015 Sarez earthquake from Jin et al. (2022) as a source fault along with the coseismic slip distribution models of the four earthquakes to calculate the static ΔCFS on the major fault planes in the region (Fig. S6). Adding the post-seismic slip barely alters the coseismic ΔCFS distributions (Fig. S6), suggesting that the post-seismic effect is small compared to that of coseismic slip on regional Coulomb stress evolution in the early 5 yr after the 2015 Sarez earthquake. In addition, the post-seismic seismic moment (2.04 × 10¹⁸ N·m) is less than 3 per cent of the seismic moment released by the coseismic rupture (9.09 × 10¹⁹ N·m; Jin et al. 2022).

6 CONCLUSIONS

Coseismic displacements of the four earthquakes from 2015 to 2023 in Pamir region have been derived from Sentinel-1A/B satellite radar observations. We determined the fault geometries and slip distribution models of the four events based on elastic half-space dislocation model using unified InSAR observations as constraints. Based on the fault slip models, we investigated the interactions of the four strong seismic events and assessed future seismic hazards.

The 2015 Sarez event promoted the rupture of the 2016 Sary-Tash earthquake, but delayed the occurring of the 2016 Muji earthquake, if we calculate the ΔCFS evolution only using the four fault models obtained in this study as the source. When we consider the pre-stresses generated by the 1974 Markansu earthquake and the 2008 Nura earthquake, the preceding earthquakes had a significant triggering effect on the two 2016 M_w 6.4+ earthquakes. Cumulatively, the three earthquakes prior to the 2023 Murghob event had a strong triggering effect on the 2023 fault, with a stress loading of more than +2 bar on the 2023 fault plane.

According to the stress loading and unloading effects on major fault planes in Pamir region generated cumulatively by the preceding earthquakes, we hold that future attention should be paid in (1) the segment of the Sarez-Karakul fault north of the Kokuibel Valley, (2) the segment of the Sarez-Murghab thrust fault west of the Sarez-Karakul fault and (3) the east segments of the Pamir thrust fault system.

Acknowledgements

The authors are grateful to Kosuke Heki and the other two anonymous reviewers for their constructive comments and valuable suggestions. This work was supported by the National Natural Science Foundation of China (41941019), the Shaanxi Key Science and Technology Innovation Team Project (Ref. 2021TD-51), the Shaanxi Province Geoscience Big Data and Geohazard Prevention Innovation Team (2022), the Fundamental Research Funds for the Central Universities, CHD (Refs. 300102260301, 300102261108 and 300102262902). In addition, we are very grateful to Yong Zhang for providing the rupture model of the 2023 Murghob earthquake and Wasja Bloch for providing the aftershock relocation data of the 2016 Sary-Tash earthquake.

DATA AVAILABILITY

Sentinel-1 SAR images were downloaded from European Space Agency. GMTSAR software package was downloaded from https://topex.ucsd.edu/gmtsar/. GACOS data were downloaded from the Generic Atmospheric Correction Online Service for InSAR (http://www.gacos.net/). The PSOKINV software we used to conduct the geodetic modelling in this study can be freely accessed from https://github.com/wpfeng/psokinv (accessed on 6 November 2023). Most figures were plotted using the Generic Mapping Tools version 6.2 (https://www.generic-mapping-tools.org).