The N-S direction strike-slip activities in the Pamir hinterland under oblique convergence: the 2015 and 2023 earthquakes

SUMMARY

The prominent Pamir plateau holds considerable significance in comprehending the processes of Asian continental collisional orogeny. However, due to harsh natural conditions and low seismic activity within the Pamir hinterland, our understanding of this region remains deficient. Recent major events and the accumulation of geodetic observations present a rare opportunity for us to get insights into the tectonic activities and orogenic processes occurring in this region. First, employing Sentinel-1 and Advanced Land Observation Satellite (ALOS)-2 Synthetic Aperture Radar (SAR) images, we acquire coseismic displacements associated with the most recent earthquakes in 2015 and 2023. Subsequently, we conduct the source models inversion with the constraints of surface displacements based on a finite-fault model. Our results reveal displacements ranging from −0.8 to 0.8 m for the 2015 M_w 7.2 Tajik earthquake and −0.25 to 0.25 m for the 2023 M_w 6.9 Murghob event, respectively. The optimal three-segment model for the 2015 event ruptured a fault length of 89 km with a surface rupture extending 59 km along the Sarez–Karakul fault (SKF), characterized predominantly by left-lateral strike-slip motion, with a maximum slip of 3.5 m. Meanwhile, our preferred uniform slip model suggests that the 2023 event ruptured an unmapped fault in the southern Pamir region with a strike angle of 31° and a dip angle of 76.8°. The distributed slip model indicates that the 2023 event ruptured a fault length of 32 km, resulting in an 8 km surface rupture. This event is characterized by left-lateral strike slip, with a peak slip of 2.2 m. Secondly, the Coulomb stress calculations demonstrate that the 2023 event was impeded by the 2015 event. Finally, interseismic Global Positioning System data revel a relative motion of 3.4–5.7 mm yr⁻¹ in the N-S component and 3.2–3.8 mm yr⁻¹ in the E-W component along the SKF in the Pamir hinterland, respectively. These N-S direction strike-slip activities and slip behaviours support an ongoing strong shear and extension in the Pamir regime, which is a response to the oblique convergence between the Indian and Eurasian plates.

Satellite geodesy, Seismic cycle, Earthquake source observations, Continental neotectonics, Dynamics and mechanics of faulting

1 INTRODUCTION

The salient Pamir plateau, also known as western Himalayan syntaxis, situated at the northwestern end of the Himalayan–Tibetan Orogen, is encompassed by four tectonic units: the Indian plate to the south, the south Tienshan Orogen to the north, the Tajik Basin to the west and the Tarim Basin to the east (e.g. Burtman & Molnar 1993; Fig. 1). As a result of the Indo-Eurasian plate collision during the late Palaeocene (e.g. Najman et al. 2017; Chen et al. 2018), the Pamir experienced 300–900 km N-S shortening in its interior since the Cenozoic (Sobel et al. 2013; Thompson Jobe et al. 2017). The ongoing convergence of the Indo-Eurasian plates gives rise to the formation of the Pamir, characterized by prominent topography including a thick crust measuring ∼60–70 km, an average high elevation of 4000–5000 m, and a rare intracontinental deep subduction (Schneider et al. 2013; Sippl et al. 2013; Bloch et al. 2021). According to the Global Centroid-Moment-Tensor (GCMT) earthquake catalogue, the Pamir region has experienced intense seismic activity, including a shallow seismicity zone along the Pamir arcuate margin and a vigorous ‘S’ shape intermediate-depth seismicity zone (> 200 km) beneath the Hindu Kush (Mechie et al. 2012; Sippl et al. 2013; Schneider et al. 2019). Kufner et al. (2016) proposed an ‘indenter model’ to account for the Pamir's protrusion hundreds of kilometres north of the Hindu-Kush and Tibetan plateau, which effectively explains the observed surface features. As one of the most active regions in Central Asia, the Pamir serves as a natural laboratory for comprehending the growth and continental collisional orogenesis of the Himalayan–Tibetan Orogen (Thompson Jobe et al. 2017; Chen et al. 2018; Perry et al. 2019).

Figure 1.

Tectonic setting of the Himalayan–Tibetan Orogen and its surrounding active blocks. (a) Distribution of Asian orogenic belts associated with the collision of Indian and Eurasian plate. The large and small dashed rectangular indicates the region plotted in panels (b) and (c), respectively. (b) Topography and major active faults, and seismic activities in the Pamir and its adjacent region. The coloured beach balls depict the earthquakes (M ≥ 4.5) recorded by the GCMT catalogue (since 1978) in different depths (< 40 km with blue colour, 40–80 km with yellow colour and >80 km with green colour). Black lines depict the major faults ans structures in the Pamir, involving, from south to north, MBT, MMT, MKT, RPS, TTS, MPT and SKF. Uppercase letters of N, C and S represent the northern, central and southern tectonic units of Pamir separated by PRS and TTS (Schurr et al. 2014). The star indicates the location of the 2023 M_w 6.9 Murghob event. (c) A detailed active structure in the Pamir hinterland. Beach balls represent the major recent earthquakes in the Pamir hinterland and its northern margin. Various coloured lines indicate the major and secondary faults. Circle points indicate the aftershocks within one month of the 2023 event from USGS. Light yellow zones indicate some domes in the central Pamir, and azure zones indicate the magmatic rocks formed during the Jurassic to Cenozoic (Zanchetta et al. 2018). White rectangles depict the footprints of the SAR images used for the 2023 event in this study, while these for the 2015 event have not been plotted that can be found in the Sangha et al. (2017). Green square point represents a major city of Murghob near the epicentre.

Present geodetic observations reveal a crustal shortening rate of 24 mm yr⁻¹ spanning a 300–400 km N-S distance across the Pamir regime, including a shortening rate of ∼10 mm yr⁻¹ at the Hindu Kush and ∼12 mm yr⁻¹ at the Main Pamir Thrust (MPT) system (Li et al. 2012; Ischuk et al. 2013; Zhou et al. 2016; Metzger et al. 2020). Along the Pamir's leading edge, strong crustal shortening and lateral extrusion are localized on some major thrust and strike-slip faults, such as the Karakorum fault to the west, exhibiting a left-lateral strike-slip rate of 6.9–10.8 mm yr⁻¹ (Robinson 2009; Cowgill 2010), the Darvaz fault to the east, displaying a right-lateral strike-slip rate of 10–15 mm yr⁻¹ (Furuya & Satyabala 2008), and the MPT system to the north, with a thrust-slip rate of 12 mm yr⁻¹ (Zhou et al. 2016). In the Pamir hinterland, there is only one major NNE-SSW striking Sarez–Karakul fault (SKF, Fig. 1). By utilizing local seismic records, Schurr et al. (2014) described the strike-slip and extensional tectonics within the Pamir interior, suggesting minimal seismicity and internal deformation in this region. To elucidate the prevailing deformation pattern of the Pamir, various kinematic models have been proposed, including the synorogenic extension model (Brunel et al. 1994), the radial thrusting model (Thomas et al. 1994), the oroclinal bending of the Pamir–Nanga Parbat syntaxis system (Robinson et al. 2004), the gravitational collapse model (Thiede et al. 2013) and the indenter model (Kufner et al. 2016). By using receiver functions and seismic tomography, some recent investigations reveal a complex dynamic process underlying the Pamir indenter, involving the interaction between northern and western Pamir (Schneider et al. 2019), intermediate-depth seismicity and its presence (e.g. Bloch et al. 2021; Xu et al. 2021).

Intense seismic activities characterize the Pamir arcuate margin, yet only a limited number of moderate-size major shallow events have been recorded instrumently in Pamir hinterland. The 2015 M_w 7.2 Tajik earthquake, being the first major event extensively observed by geodetic measurements in the Pamir hinterland, garnered attention from scholars worldwide. Based on the analysis of geodetic and seismic waveform data, previous studies suggest that the 2015 event ruptured the NNE-SSW striking SKF spanning a length of 80 km, predominantly characterized by left-lateral strike-slip motion (e.g. Metzger et al. 2017; Sangha et al. 2017; Elliott et al. 2020). In addition, Bloch et al. (2023) investigated the foreshock, main shock and aftershock sequences in the Pamir during the 2015–2017 by using the local seismic network and suggested the facilitated fluid migration caused by the main shock should trigger the following M_w  > 6 earthquakes. These source models provide valuable insights into the comprehension of tectonic activities in the Pamir hinterland. Due to the harsh natural conditions (i.e. extreme topographic relief, and broad permanent ice cover) and political unrest, geodetic observations within the Pamir hinterland remain very sparse, impeding a thorough understanding of this region.

On February 23 at 0:37 UTC (5:37 local time), a strong M_w 6.8 earthquake struck the west of Murghob, Tajikistan, in the Pamir hinterland. The event was perceptible in distant locations such as the city of Dushanbe, which serves as the capital of Tajikistan, as well as various other regions within the central Asian nation (Fig. 1c). The U.S. Geological Survey (USGS) estimated that over 2000 people were exposed to this main shock. Fortunately, no casualties or structural damage have been reported. In the month subsequent to the main shock, numerous aftershocks ranging from M 4.0–5.0 were observed predominantly in the N-S striking direction. The preliminary focal mechanism results from both USGS and GCMT suggested that the 2023 event was characterized by a nearly pure strike-slip fault rupture with a high dip angle, while their respective source parameters exhibited significant disparities. The occurrence of the 2023 event renews our focus on the Pamir hinterland, thereby raising numerous crucial issues regarding the tectonic activities within this prominent orogenic belt. For instance, what are the detailed rupture features of this 2023 event? If this 2023 event also ruptured a strike-slip fault oriented approximately in N-S striking direction, are there any distinctions compared to the 2015 event? Given the proximity of their epicentres, could there be any influence of the 2015 M_w Tajik earthquake on the occurrence of the 2023 event?

To address these inquiries, comprehensive geodetic observations and robust source models are indispensable. In this study, we analyse the coseismic deformation and source models of both the 2015 and 2023 events using L-band ALOS-2 and C-band Sentinel-1 SAR images. Furthermore, we discuss potential fault interactions between the 2015 Tajik event and the 2023 Murghob event. Finally, we integrate available interseismic Global Positioning System (GPS) observations to estimate the slip rate along the SKF striking direction in the Pamir hinterland, aiming to explore how fault activities modulate the interior shortening of orogenic belt and their possible implications in central Asia.

2 TECTONIC SETTING

The Pamir plateau stands out prominently within the Himalaya–Tibetan orogenic belt due to its northward convex shape and deflection, resulting from the amalgamation of various tectonic terranes along-strike by a westward prolongation of the Himalaya–Tibetan Plateau during the India-Asia convergence (Robinson 2009; Sobel et al. 2013; Li et al. 2018, 2020). It is separated by a series of major arcuate faults and sutures, which include, from south to north, the Main Boundary Thrust fault (MBT), the Indus Suture (Main Mantle Thrust fault, MMT), the Shyok Suture (i.e. the Main Karakoram Thrust fault, MKT), the Rushan–Pshart Suture (RPS) zone, the palaeotethys Tanymas Thrust Suture (TTS) zone and the MPT system (including the Main Frontal Thrust and the Pamir Frontal Thrust, i.e. MFT and PFT), respectively. Additionally, it is bounded by two major marginal strike-slip faults, the Darvaz fault on the west flank margin and the Karakorum fault on the east flank margin (Fig. 1b). The Pamir regime generally originates from several Mesozoic–Palaeozoic terranes, and has been subdivided into three major geological units by the TTS and RPS, that is, the northern Pamir, central Pamir and southern Pamir (Burtman & Molnar 1993; Robinson 2009). The northern Pamir constitutes the outer margin of the salient that consists of the Asian affinity terranes, and the central and southern Pamir consist of the Cimmerian Gondwanan terranes (Zanchetta et al. 2018). The Pamir plateau, as a whole, displays extreme elevations, an over-thickened crust, and inherited structures that potentially influence its contemporary deformation, fault occurrence, local tectonic features and seismic activity (Elliott et al. 2020).

In the Pamir hinterland, it is characterized by mountainous terrain, Cenozoic domes, and pull-apart basins. The late Cenozoic tectonic activities, involving crustal extension collapse, doming, and exhumation, are generally related to three major types of tectonic units (Sobel et al. 2013; Rutte et al. 2017; Fig. 1c). (1) Large-scale gneiss domes are broadly distributed in the Pamir hinterland, such as the Sarez dome, Msukol dome and Yazgulem dome in central Pamir. These E-W extended domes are regarded as the upper crustal syntectonic extension, indicative of ceased N-S extension in the last 2 Ma and ongoing E-W extension (Robinson et al. 2004; Ischuk et al. 2013; Stübner et al. 2013; Zanchetta et al. 2018). (2) The reactivated RPS zone, oriented E-W, hosts several isolated slices of serpentinite that corroborate the closure of the Rushan–Pshart Ocean, and it was reworked during the Cenozoic shortening (Zanchetta et al. 2018). (3) The NNE-SSW striking SKF system, developed along the transtensional Karakul graben, intersects the Muji fault on the northern Pamir thrust system and extended southward to the Yashykul Lake. It cuts across these convergent shear zones and divides the northward convergent Pamir into the eastern and western parts (Kulikova et al. 2016). Considering the slip behavious of the SKF, Cowgill (2010) suggested that the deflection structures and deformation features on the west and east flanks of the Pamir are asymmetric. In addition to these major tectonic units, numerous secondary shallow faults have been inferred based on geological landform data (Schurr et al. 2014). Note that lots of potentially active faults can be found along the NNE-SSW striking direction in the southern Pamir. These shear and extensional structures support a complicated convergent and extensional evolution in the Pamir hinterland (e.g. Searle & Hacker 2019), resulting from either deep slab subduction or gravitational collapse, or both of them (e.g. Sippl et al. 2013; Zanchetta et al. 2018).

Recent earthquake catalogue and geodetic observations reveal shallow tectonic activities in the Pamir hinterland at a low level. The 1911 M_s 7.3 Sarez earthquake was the largest event recorded instrumentally in the region, and it triggered the largest single non-volcanic landslide in the past century (Kulikova et al. 2016). Schurr et al. (2014) speculated the 1911 event was caused by either the reactivation of the E-W trending RPS zone or the sinistral transtensional SKF system. Due to limited global seismological observatories at that time, the focal mechanism of the 1911 event was enigmatic and its location was with largely uncertain (Fig. 1c). After 100 yr, the 2015 M_w 7.2 event has been the second major event occurred in this region, and also being the first major event well documented through a global seismic network and geodetic observations (Kulikova et al. 2016; Sangha et al. 2017; Elliott et al. 2020). Previous studies suggested that the 2015 event ruptured the SKF with a length of ∼ 80 km, and its slip reached the surface (e.g. Metzger et al. 2017; Sangha et al. 2017; Elliott et al. 2020). In comparison to the 2015 events, the 2023 Murghob earthquake was presumably inferred as a similar strike-slip event occurring near the RPS zone in southern Pamir, but with an epicenter located further south.

3 DATA SET AND METHODS

3.1 SAR images and processing

Benefiting from the burgeoning remotely sensed satellites and their improved image quality (e.g. high spatial resolution in large-scale, short revisit-time and optimal coherence), imaging geodesy observations have become a routine avenue to capture surface displacement, proving successful in global earthquake monitoring and response. In this study, both the available L-band ALOS-2 Phased Array type L-band Synthetic Aperture Radar (PALSAR) and C-band Sentinel-1 SAR images are collected to derive the coseismic displacements associated with the 2015 and 2023 events, respectively. For the 2015 event, we selected the SAR images referring to Sangha et al. (2017). Given that one of the ALOS-2 pairs used in the Sangha et al. (2017) covers only a small part of the coseismic deformation region, it has been excluded from our study. The other four InSAR pairs, including both ascending and descending Sentinel-1 and ALOS-2 images, have been used in this study (detailed information listed in the Supporting Information, Table S1). These pairs for the 2015 event exhibit a short perpendicular baseline and temporal baseline. For the 2023 event, the selected Sentinel-1 images include two ascending pairs and one descending pair, while the ALOS-2 images comprise two ascending pairs. As shown in Fig. 1(c), each pair of Sentinel-1 images provides complete coverage of the epicentral region. The two ALOS-2 images pairs with the same path orbit number cover a small area, but they capture complementary deformation patterns in the epicentral region. Note that most of image pairs used for the 2023 event have a long temporal baseline of 300–400 d, as this region predominantly exhibits low coherence due to its steep topography and glacier coverage (detailed information listed in the Supporting Information, Table S2).

All of these images are processed with the two-pass differential Interferometric Synthetic Aperture Radar (InSAR) method based on the commercial GAMMA platform (Wegnüller et al. 2016). During the interferometric processing, the 1-arcsec Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) is used to remove the topographic phase contribution (Farr et al. 2007), and a powerful spectral filter is used to improve the signal-to-noise ratio (SNR). We unwrap each interferogram with the Statistical-cost, Network-flow Algorithm for Phase Unwrapping (SNAPHU) algorithm, followed by a transformation from the radar coordinate frame to the geographic coordinate of the SRTM DEM. Finally, we correct all potential long-wavelength errors (i.e. orbital and atmospheric errors) in the interferograms by applying a linear ramp correction. To evaluate the precision of these final interferograms, 1-D covariance function was used to calculate their standard deviation (e.g. Parsons et al. 2005). We also use the pixel offset tracking method to investigate the range and azimuth offsets of the 2023 event, while no available displacement signals have been retrieved due to its low SNR. Given that each interferogram generally comprises millions of data points, conducting inversion computations becomes nearly infeasible. Thus, subsampling with a suitable point data set is imperative for source modelling. The computing efficiency and robustness of inversion were balanced by subsampling each interferogram based on its displacement gradient using the Quadtree algorithm (Jónsson et al. 2002). This approach retained sufficient information within about 500–1500 points to describe the deformation features.

3.2 Source modelling

A finite-fault model with nine parameters in elastic half-space (i.e. the longitude, latitude and depth of the fault location, the length and width of the fault plane size, the strike angle and dip angle of the fault plane orientation, and the slip and rake angle of the fault plane motion) is commonly used to describe the coseismic surface displacement resulting from the subsurface fault rupture of an earthquake (Okada 1985). Resolving these parameters, with the constraint of surface displacements, poses a high-dimensional nonlinear inversion problem. Two primary strategies for source modelling have been extensively employed. The first strategy involves a traditional two-step inversion procedure, commencing with a nonlinear inversion predicated on a uniform slip model (USM), succeeded by a linear inversion utilizing a distributed slip model (DSM). The objective of the USM is to determine the fault geometry (the fault location, fault plane size and its orientation) by assuming uniform slip across the rectangular fault plane, rendering it suitable for modelling simple rupture scenarios such as single planar fault ruptures. Various nonlinear optimizations, such as simplex search algorithm, simulated annealing, and Bayesian posterior probability density search, have been carried out to search a global optimal solution of nine parameters of the USM with the constraints from geodetic observations (e.g. Clarke et al. 1997; Cervelli et al. 2001; Bagnardi & Hooper 2018). By fixing the fault geometry derived from the USM, the fault plane can be discretized into numerous smaller patches, enabling subsequent employment of DSM to invert for detailed variable rupture slip vectors on each patch. The second strategy is specifically designed for distributed slip inversion. When the slip rupture exhibits a large scale and/or highly complex pattern, it becomes inadequate to simply approximate its fault geometry with a single planar surface. Under such circumstances, multisegment or curve plane fault geometry would be more suitable. To reduce potential complexities in fault geometry, certain fault geometry parameters can be fixed or constrained by using a priori information, such as relocated aftershocks, surface rupture traces, topography and landforms, and local seismic reflection profiles (e.g. Wang et al. 2009; Elliott et al. 2020). Subsequently, the DSM on each of these structures can be inverted to determine a final optimal slip solution model by minimizing fitting residual.

The 2015 event was caused by a complex large-scale curve subsurface slip, and some previous studies proposed a multisegmented model to invert its source modelling (e.g. Metzger et al. 2017; Sangha et al. 2017; Elliott et al. 2020). These studies constrained some of the fault geometry tightly (e.g. fault striking and length) by using the surface rupture trace from the imaging geodesy observation, and then searched the other source parameters. For example, a three-segment fault model was constrained by ALOS-2 and Sentinel-1 interferograms in Sangha et al. (2017), while a nine-segment fault model was constrained by Sentienl-1 interferograms and optical images offsets in Elliott et al. (2020). Based on these previous studies and our coseismic interferograms, a three-segment fault model is proposed to approximate the curve strike-slip rupture in this study. The fault location, strike angle and length of these segments are fixed, and adjacent segments are connected with each other at their updip ends (detailed information listed in Supporting Information, Table S3). Our three-segment fault model includes a southwest segment of 47.5 km, a central segment of 17 km, and a northeast segment of 37.6 km, and their strike angles are 208.8°, 231.5° and 215.2°, respectively. The maximum slip is confined to < 3.5 m for all these three segments. To determine the dip angle of each segment fault, the initial dip angle of each segment ranges from 60° to 90° with an interval of 1°, followed by a grid search method to determine the final optimal model from the candidate DSM models through SDM (steepest descent method) inversion (Supporting Information, Fig. S1).

For the 2023 event, we employ the classic two-step inversion procedure. First, we utilize a Bayesian inverse approach developed by Bagnardi & Hooper (2018) for the USM inversion. This Bayesian method considers the model fitting residual as well as the posterior probability density of the source parameters, and the initial parameters can be defined in a wide range (detailed information listed in Supporting Information, Table S4), thereby enabling a much more robust parameter estimation as well as completely search in underestimated parameter space. A total of 1 × 10⁶ iterations has been implemented in the Bayesian inversion. Results of the first 1 × 10⁴ iterations are considered as a burn-in period and have been excluded from analysis, and then the final optimal model is estimated based on these remaining candidate models using a posterior probability distribution statistic (Supporting Information, Fig. S2 and Table S5). Based on the fixed fault geometry structure obtained from the USM, an SDM developed by Wang et al. (2009) is used to invert the variable slip distribution. To avoid the influence of edge effects, the fault plane has been extended to a size of 40-km length × 32-km width,with its upper boundary reaching the surface. The fault plane was discretized into a total of 320 subpatches, each measuring 2 km × 2 km. In addition, the slip magnitude and rake angle in the DSM inversion range from 0–3.5 m, and −60° to 60°, respectively.

3.3 Coulomb stress calculation

To evaluate the earthquake interactions and their effects on the local seismicity in the Pamir hinterland, we employ a Coulomb failure criterion with an index of |${{\Delta }}{\rm CFS}$| to analyse the coseismic stress field changes caused by two shallow moderate earthquakes (i.e. the 2015 Tajik event and the 2023 Murghob event). By using the Coulomb v3.4 package (Toda et al. 2005), we calculate the static Coulomb stress change induced by the 2015 event on the 2023 event, as well as the static Coulomb stress change induced by both the 2015 and 2023 events on a receiver fault similar to the 2023 event, respectively. The source models of the 2015 and 2013 events are from the inverted DSM in this study. When the 2023 event is considered as receiver fault, its fault nodal plane (31°/76.8°/−9.8°) is constructed based on its USM. For all the stress change calculations, a friction coefficient of 0.4 is used.

3.4 Interseismic GPS data integration

GPS has been recognized as a revolutionary tool for investigating the crustal motion, local strain field and fault tectonic activities. It has been rapidly developed during the past three-decade years (Kreemer et al. 2014). Benefit from the establishment of the crustal movement observation network of China (CMONOC), as well as the implementation of some local GPS networks, the GPS spatial density of the Himalayan–Tibetan Orogen and its neighbourhood are greatly increased (e.g. Kreemer et al. 2014; Wang & Shen 2020; Li et al. 2022). However, the current GPS measurements in Pamir remain limited in terms of density due to restricted accessibility and a lack of a stable campaign and/or continuous GPS network. Over the past few decades, several episodic campaign GPS measurements in different parts of the Pamir have been carried out (e.g. Kreemer et al. 2014; Zhou et al. 2016; Perry et al. 2019). By integrating these different GPS data sources, we can enhance the spatial density of GPS measurements and facilitate an investigation into orogenic motion in the Pamir hinterland. Referring to the most recent large-scale GPS horizontal velocity results proposed by Wang & Shen (2020), we collect these published GPS measurements in the Pamir and align them with a same stable Eurasia reference frame by using a rotation Euler pole. The final integrated GPS horizontal measurements are used to analyse the slip rate in the Pamir hinterland. By delineating two GPS profiles across the SKF strike direction in Pamir hinterland, we estimate the slip rates parallel and perpendicular to the target fault stike direction.

4 COSEISMIC OFFSETS AND SOURCE MODELS

4.1 The 2015 M_w 7.2 Tajik event

In Figs 2(a)–(d), both the ascending and descending Sentinel-1 and ALOS-2 interferograms clearly exhibit symmetrical opposite displacements in line of sight (LOS) on each side of the SKF, indicating a strike-slip rupture with high dip angle for the 2015 Tajik event (positive LOS displacements indicate the ground motion towards the satellite, and this definition is consistently used hereinafter). The displacements of these interferograms range from −0.8-to 0.8 m, and certain discrepancies among them can be attributed to the satellite observing geometries. With reference to the epicentral location determined by the USGS (Fig. 1), these displacements initiated from the Yashykul lake on the southern boundary of the Sarez dome, and then propagated northward along the SKF over a distance of 60–70 km (Fig. 2). The Sentienel-1 interferograms in Figs 2(a) and (b) exhibit a larger uncorrelated region along the SKF, resulting from the unwrap failure. Different from the C-band Sentinel-1 images, the L-band ALOS-2 interferograms in Figs 2(c) and (d) exhibit much higher coherence that can capture complete coseismic displacements in near filed of this 2015 event. Therefore, some previous studies that solely relied on Sentinel-1 images for the source modelling of this 2015 event would lack some crucial near field information. Note that we used the same InSAR pairs as Sangha et al. (2017), while the final interferograms in this study exhibit significantly enhanced quality, such as clear and complete fringes, owing to an improved registration accuracy. The standard deviations for T005D, T100A, P057D and P163A are 1.5, 2.1, 2.7 and 3.1 cm, respectively (Supporting Information, Table S1). For simplicity, the subsampled points of these interferograms are assigned equal weight in the source modelling.

Figure 2.

Coseismic displacement in LOS for the 2015 Tajik event. The Sentienl-1 interferograms are shown in (a) and (b) for the descending track T005D and ascending track T100A, respectively. The ALOS-2 interferograms are shown in (c) and (d) for the descending track P057D and ascending track P163A, respectively. Along the dashed AA′ line, detailed observed displacements are shown in the corresponding insert map. Positive LOS displacements indicate the ground motion toward the satellite, and this definition has been used in the following Figs 4, 5 and 7.

The DSM of our three-segment fault model is shown in Fig. 3. The final dip angles for each segment of our three-segment fault model are 87°, 70° and 87°, respectively (the root-mean-square error of each DSM with candidate dip angle are shown in the Supporting Information, Fig. S1). The fault geometry structure of the northern and central segments in this study exhibits similarities to that proposed by Elliott et al. (2020), whereas the southern segment aligns with the results of Sangha et al. (2017). The main reasons for these differences can be attributed to two factors. First, the fault surface trace segmentation can be tightly constrained by comprehensive near-field observations from ALOS-2 interferograms in this study or the optical POT images in Elliott et al. (2020) (i.e. SPOT 6–7 and Landsat-8 images). Secondly, both ALOS-2 and Sentinel-1 images have been utilized for source modelling in this study and Sangha et al. (2017), while only Sentinel-1 images were used in the work of Elliott et al. (2020). The significant slips exceeding 0.5 m cover an area of 22-km length × 30-km width within the northern segment, exhibiting a ‘pear’ shape, 17-km length × 12-km width in the central segment and 36-km length × 18-km width in the southern segment, respectively, suggesting a total subsurface slip length of 65 km. The surface rupture spans approximately 51 km (slip > 0.5 m) with a peak slip of 3.5 m can be found in Fig. 3, corresponding to the uncorrelated region along the SKF. If considering slip with > 0.2 m, this 2015 event ruptured a fault length of 89 km and surface rupture of 59 km, which contrasts previous estimations by some studies that suggested a fault length of 80 km and surface rupture of 35–40 km (e.g. Metzger et al. 2017; Sangha et al. 2017). Note that the 2015 event is dominated by left-lateral strike slip, while some normal slip components or thrust slip components exist in different segment along the SKF. The total released seismic moment yields to 8.1 × 10¹⁹ N·m, equivalent to M_w 7.2, similar to that of GCMT and USGS. The synthetic and residual displacements derived from DSM depicted in Fig. 4 exhibit a well-fitting between the observation and model, except for some large misfits near the eastern side of the central segment (slip uncertainty and resolution testing of the 2015 event are shown in the Supporting Information, Figs S3 and S4).

Figure 3.

Coseismic slip distribution of the three-segment structure for the 2015 event. Arrows show the sense of motion of the hangingwall and contours are meters of coseismic slip.

Figure 4.

(a)–(d) Synthetic and (e)–(h) residual displacement for interferograms of T005D, T100A, P057D and P163A, based on the DSM in Fig. 3. Red line indicates fault line projected on the surface along its dip plane upward, and black rectangular frame is the slip plane in subsurface. Along the dashed AA′ line, detailed synthetic displacements are shown in the corresponding insert map. Statistical histograms indicate the distribution of residuals in each panel.

4.2 The 2023 M_w 6.9 Murghob event

The interferograms presented in Figs 5(a)–(d) provide complete coverage of the epicentral deformation region associated with the 2023 event. The two descending Sentinel-1 interferograms shown in Figs 5(a) and (b) reveal deformation ranging from −0.25 to 0.25 m, exhibiting a similar shape of ‘butterfly’. These displacements in a four-quadrant distribution include two large lobes and two small lobes on the southeast and northwest sides, respectively. The displacements in Fig. 5(c) range from −0.18 to 0.18 m, displaying a symmetric deformation but oriented oppositely with respect to the N-S striking axis. Considering the different observing geometry of SAR satellite, the symmetric displacements in both descending and ascending interferograms suggest that the LOS displacements are mainly contributed by the horizontal components, with minimal contribution from vertical components. This pattern commonly corresponds to a strike-slip rupture event characterized by a high dip angle. The entire ALOS-2 interferogram includes a left and a right pair, and their displacements range from −0.5 to 0.29 m (Fig. 5d). The insert map in each interferogram exhibits the coseismic displacement along profile AA′, revealing a symmetrical logarithmic decrement as distance increases from the epicentral location. In comparison to the L-band ALOS-2 interferogram, some uncorrelated regions can be found in the epicentral location in all three C-band Sentinel-1 interferograms, which may be attributed to either a glacier cover or phase unwrapping failure. The complete near-field displacement of the ALOS-2 interferogram suggests the fault may rupture along an NNE striking. The standard deviations for T005D, T107D, T100A, P163Al and P163Ar are 0.6, 0.6, 1.9, 4.6 and 2.6 cm, respectively (Table S2, Supporting Information), suggesting that the precision of the Sentinel-1 interferograms is nearly three times higher than that of the ALOS-2 images. Given the total subsampled points of Sentinel-1 images is about three times than that of the ALOS-2 images, we also assign equal weights for all these InSAR measurements in the following source modelling.

Figure 5.

Observed, synthetic and residual displacements in LOS for the 2023 Tajik event based on the USM. Three Sentienl-1 interferograms are shown in (a)–(c) for the descending track T005D, the descending track T107D and the ascending track T100A, respectively. The ALOS-2 interferograms consist of two subpairs (left and right) are shown in (d) for the ascending track P163A. Their relative synthetic and residual displacement based on the optical USM in this study are shown in (e)–(h) and (i)–(l), respectively. Along the dashed AA′ line, detailed observed and synthetic displacements are shown in the corresponding insert map. Statistical histograms indicate the distribution of residuals in each panel, respectively.

The optimal USM solution suggests that the 2023 event is characterized by an NNE left-lateral strike slip with a high dip angle, which does not align with any known fault in this region. The 2023 event ruptured a fault with a strike angle of 31°, a dip angle of 76.8° and a rake angle of −9.8°. The uniform slip is confined by a rectangular plane with a size of 17-km length × 16-km width (detailed information can be found in the Supporting Information, Table S4). The uniform slip consists of a 1.9 m strike-slip and a −0.33 m dip-slip, indicating that the rupture is dominated by strike slip with a minor normal slip. Given that the updip boundary of the USM at a shallow depth of 3.1 km and a significant total uniform slip of 1.93 m, the rupture of this 2023 event may reach the surface. In comparison to the focal mechanisms of both USGS and GCMT, our USM does not support either of the two nodal planes. Note that the first nodal plane of these seismic waveform solutions also exhibits a left-lateral strike-slip rupture with a high dip angle, but in the opposite dipping direction. It is understandable that the seismological data is hard to distinguish the dipping direction when the rupture occurred on a strike-slip fault with a high dip angle. If considering an elastic shear module of 30 GPa, the USM released a total released seismic moment of 1.69 × 10¹⁹ N·m, which is smaller than that reported by GCMT and USGS. The forward displacements and corresponding residuals are shown in Figs 5(e)–(h) and (i)–(l), respectively. Regarding the Sentinel-1 residuals, only some large misfits (within −4 cm) can be found near the hypocentral region, such as a region with light blue coloūr on the west side, and these discrepancies can be attributed to the model error. As for the ALOS-2 residuals, there are also some large misfits (∼6–7 cm) in the southwest far field, which can be attributed to atmospheric or/and ionospheric error. However, these residuals generally exhibit a random pattern with most of them are less than < 4 cm, indicating that the InSAR observations are well explained by the USM.

As shown in Fig. 6, the variable slip distribution renders only one simple asperity characterized by an elliptic pattern. The significant slip exceeding 0.5 m covers an area of 32-km length × 29-km width, and it is confined to a depth of 0–28 km with a peak slip of 2.2 m. Similar to the 2015 event, the slip vectors of the 2023 event almost exhibit a pure left-lateral strike slip, except for some normal slip components around the max slip region and a little thrust slip close to the surface. Note that the main slip has reached the surface with a length of 8 km, suggesting a surface rupture resulting from this event. Therefore, the uncorrelated region of the interferograms near the causative fault can be attributed to the an unwrapping failure caused by a phase gradient exceeding 2|$\pi $| between adjacent pixels in the interferogram when there exists a shallow fault rupture. The DSM yields a geodetic moment of 2.9 × 10¹⁹ N·m, equivalent to M_w 6.9, which is larger than that inferred from the USM, GCMT and USGS. The synthetic LOS displacements and their corresponding residuals based on the DSM are shown in Fig. 7 (slip uncertainty and resolution testing of the 2015 event are shown in the Supporting Information, Figs S5 and S6). The Sentinel-1 residuals are within 3 cm (Figs 7e–g), and large residuals observed in the ALOS-2 interferograms (Fig. 7h) are mainly from the systematic offsets between the left and right pairs. The large misfits near the epicentral region in the USM (i.e. Figs 5i–l) are absent in the DSM (i.e. Figs 7e–h), suggesting that the model error has been greatly reduced.

Figure 6.

Coseismic slip distribution for the 2023 event in this study. Arrows show the sense of motion of the hangingwall and contours are meters of coseismic slip.

Figure 7.

(a)–(d) Synthetic and (e)–(h) residual displacement for interferograms of T005D, T107D, T100A and P163A, based on the DSM in Fig. 6. Red line indicates fault line projected on the surface along its dip plane upward, and black rectangular frame is the slip plane in subsurface. Statistical histograms indicate the distribution of residuals in each panel.

5 DISCUSSION

5.1 Rupture slip behaviours and fault interactions

The high topography and thick crust in the Pamir hinterland are believed to be a result of the Indo-Eurasian plate collision, yet the active deformation in the region and its relationship with the orogenic process are still poorly known due to extreme field conditions and rough terrain. Robust earthquake source models would provide crucial insights to enhance our comprehension. Currently, only three moderate shallow earthquakes (M_w|$\ge $| 6.0) have been recorded instrumently in the Pamir hinterland, indicating very low seismic activities compared to the Pamir arcuate margin. The 1911 M_w 7.3 Sarez earthquake went undetected by reliable seismic seismological observatories, let alone precise geodetic observations. Inheriting from the recent 2015 M_w 7.2 Tajik earthquake, this 2023 M_w 6.9 Murghob earthquake has been the second major event in the Pamir hinterland that can be investigated through imaging geodesy observations. InSAR observations derived from multisource SAR satellites as well as multi-observing geometries greatly reduce the effects of incoherence in the interferograms and enable a complete observation of near-field surface displacements, which imply much tighter constraints for robust estimations of source model inversions for these two recent events.

Source models of previous earthquakes suggest that active fault slip behaviors in the Pamir hinterland are dominated by left-lateral strike slip primarily oriented along the N-S direction with a high dip angle. As shown in Fig. 1(c), most of active faults (including major and secondary fault) in the Pamir hinterland run from east to west with an imbricated thrust slip feature associated with its contemporary orogenic activities, and the SKF is the only major strike-slip fault system oriented along the NNE-SSW direction. Indeed, the detailed fault slip behaviors within the Pamir interior remained poorly understood until the occurrence of the 2015 event (e.g. Metzger et al. 2017; Bloch et al. 2023). By using the re-analysed high-quality InSAR measurements, we revisit the source model of the 2015 event with a three-segment structure. A detailed slip distribution reveals that the 2015 event ruptured the southern part of the SKF with a total length of 89 km at depths ranging from 0 to 24 km in a SSW direction (Fig. 3). The major slips of this event include at least three asperities, and some slips have reached the surface. These surface ruptures have been attested by discontinuous scarps from optical images (e.g. Elliott et al. 2020). Previous studies generally agree that the 2015 event ruptured a steep multisegment strike-slip seismotectionics, while they exhibit some different detailed structure and slip features (e.g. Metzger et al. 2017; Sangha et al. 2017; Elliott et al. 2020). In all these previous models, the slip distribution in our three-segment structure model is much closer to a nine-segment model from Elliott et al. (2020), excepting some differences in the northern segment (e.g. fault dip angle and shallow slip pattern). The main reason for the discrepancies can be attributed to different near-field surface displacement constraints.

For the 2023 event, it ruptured an unmapped fault with a length of 32 km at depths of 0–28 km in the NNE direction, and its slips are solely composed of only one simple asperity. Aftershocks occurring within one month extend along N-S direction, aligning with the striking direction of the causative fault. Note that our model represents a different nodal plane that does not correspond to either of the two nodal planes determined by USGS and GCMT. A possible reason for this phenomenon could be attributed to the large uncertainties of seismic focal mechanisms in the cases of steep strike-slip event. In comparison to the 2015 event, the 2023 event exhibits some similar slip behaviors, including parallel fault strike, a high dip angle of > 70°, and a pure left-lateral strike slip. Different from a fault segmentation of the 2015 event, the 2023 event ruptured a single fault with an opposite dip trend. Assuming the fault nodal plane solution (strike/dip/rake: 35°/73°/12°) of the 1911 Sarez from Kulikova et al. (2016) is robust, it can be found that the 2023 event exhibit a similar fault nodal plane (strike/dip/rake: 31°/76.8°/−9.8°), even though they should not be repeated events as their epicentral locations are different.

Note that three major events occurred sequentially from west to east within a 100 km range in the Pamir hinterland, near the Rushan–Pshart Suture zone (Fig. 1). With a similar left lateral strike slip behavior, these three parallel causative faults should release a total relative strike-slip rupture of 8–9 m (including slip of 2.8 m for the 1911 event estimated by an empirical formula refer to Wells & Coppersmith 1994; slip of 3.5 m for the 2015 event and 2.2 m for the 2023 event refer to the DSM in this study) along the N-S direction, implying a strong partitioned slip along the N-S direction as a response to differential compression-shortening deformation within the Pamir interior. As shown in Fig. 8, the 2023 event occurred in a negative |${{\Delta }}{\rm CFS}$| region caused by the 2015 event. In addition, significant post-seismic displacements associated with the 2015 event have been only observed at the northeastern end of the rupture, far from this 2023 event (Jin et al. 2022). Therefore, a simple explanation for these fault interactions is that the occurrence of the 2023 event should be not triggered by the 2015 event. Note that Bloch et al. (2023) proposed a fluid triggering model for the 2015–2017 earthquake sequence, which could be a possible triggering explanation of the 2023 event.

Figure 8.

Coulomb failure stress changes at a depth of 10 km caused by the recent local earthquakes in the Pamir hinterland and their effects. (a) Static Coulomb stress changes caused by the 2015 Tajik earthquake imposed on the 2023 event; (b) static Coulomb stress changes caused by both the 2015 M_w 7.3 Tajik and 2023 M_w 6.9 Murghob earthquakes imposed on a receiver fault with strike of 31.3°, dip of 76.8° and rake of 9.8°. Green rectangular frames represent the rupture locations of these events.

5.2 Extensional motions in the Pamir interior

The extreme elevation of the Pamir should generate a large gravitational potential energy (GPE) that drives its extensional motion, collapses, and lateral extrusion (e.g. Jay et al. 2017; Rutte et al. 2017). By using the dynamic simulation, Jay et al. (2017) suggested that E-W extension is prominent in the Pamir interior, corresponding to the formation of lots of late gneiss domes in the Pamir hinterland with a long axis along the E-W direction (Sobel et al. 2013; Rutte et al. 2017). Comparatively, differential N-S extension with low magnitude occurred only around the gneiss domes (Jay et al. 2017), which is consistent with cessation of N-S extension in the last 2 Ma (Robinson et al. 2004; Ischuk et al. 2013; Stübner et al. 2013; Zanchetta et al. 2018). Note that these remarkable domes are exclusively located along the RPS, an alternative explanation for the extension dome could also be related to either the deep subduction of steep slab or backarc extension resulting from slab rollback (Jay et al. 2017). That is, the deep steep sab beneath the Hind Kush can facilitate the arc-along extension (partial melting, Schurr et al. 2014). Assuming that surface deformation and fault occurrence are controlled by a balance of local body forces, it is expected that obvious E-W extensional structure or/and motion would be observed throughout the Pamir interior. However, the only major E-W extensional structure in Pamir is the 250 km long Kongur Extensional system that located on the northeastern Pamir margin. Apart from some small normal displacements have been suggested on the Karakul rift based on the estimation of Quaternary offset stream channels, no other E-W extensional structures or/and motions have been attested within the Pamir interior (Robinson et al. 2004; Schurr et al. 2014).

The recent seismic activities in the Pamir hinterland greatly increase the significance of approximately NNE-SSW striking faults as a response to its interior extension. In this study, our results demonstrate that the 2015 event ruptured the southern half of the SKF in the northern and central Pamir, and the 2023 event ruptured a fault that extends into the southern Pamir (Fig. 9a). The final GPS measurements in Fig. 9(a) show a diffuse and continuous first-order deformation feature with approximately N-S shortening within the Pamir plateau, corresponding to the N-S collisional compression of the Indian–Eurasian plates. In Fig. 9(b), the NEE-trending profile A across the northern and central Pamir shows a relative motion of 5.7 and 3.8 mm yr⁻¹ parallel and perpendicular to the SKF, respectively. Note that the velocity decline across the SKF is steep and stepwise, corresponding to a left-lateral strike-slip fault that regulates the shortening of N-S compression. In Fig. 9(c), a similar profile B across the southern Pamir shows a relative motion of 3.4 and 3.2 mm yr⁻¹ parallel and perpendicular to the SKF, respectively. It can be found that the southern Pamir also experiences a differential N-S shortening motion similar to that observed in the northern and central Pamir, but with a small magnitude. These obvious interseismic strike-slip rates reveal strong shear stress within the Pamir interior, consistent with a maximum shear slip-line field oriented along the NNE direction from Wang & Barbot (2023). When the ongoing N-S compression exceeds the shear capacity of the Pamir regime, the terranes could be torn accompanied by an ‘I’ shape strike-slip faulting and earthquake rupture. These strike-slip shear zones can serve as channels of ongoing lower crustal flow from the partially molten rocks (Sass et al. 2014). In addition, both of these two profiles exhibit ∼3 mm yr⁻¹ relative motion in the E-W direction, suggesting the E-W extension should be ongoing and accompanied with the strike-slip motion within the Pamir hinterland.

Figure 9.

Interseismic GPS velocity and profile analysis in the Pamir plateau. In subfigure a, arrows represent horizontal GPS-derived interseismic velocities with respect to stable Eurasia. The black lines represent the active faults. Red lines indicate the rupture caused by the 2015 and 2023 events. The coloured beach balls depict the earthquakes (M ≥ 6.0) recorded by the GCMT catalogue (since 1978). In subfigures (b) and (c), upper dots are the velocities projected onto the direction parallel to the faults. Lower dots are the velocities projected onto the direction normal to the faults. In the upper dots, positive values represent the northward direction, while negative values represent the southward direction. In the lower dots, positive values represent the eastward direction, while negative values represent the westward direction. Dotted line is the best-fitting line of velocity data. Strips denote acceptable ranges of average velocity components.

As the thrust-slip activities migrate from the interior of the orogeny to its margins, the strike-slip and normal-slip activities become even more active during the late stage of continental collision orogeny, and provide essential information to manifest its growth and dynamic process. Different from the N-S normal-slip faulting in the central Tibet, the Pamir hinterland exhibits ongoing NNE-SSW strike-slip faulting. This difference should be attributed to the direct collision of the Indian and Asian lithospheres occurring beneath the Pamir Plateau rather than under the Tibetan plateau, leading to variations in deep subduction structures between these regions (Zhao et al. 2010) (Fig. 10).

Figure 10.

Schematic map of orogenic pattern in Pamir (referring to Perry et al. 2019).

6 CONCLUSIONS

The coseismic displacements of the most recent 2015 and 2023 earthquakes are derived from both the Sentinel-1 and ALOS-2 images in this study. Our imaging geodesy observations reveal that the 2015 and 2023 events have symmetric deformation patterns with displacement magnitudes ranging from −0.8 to 0.8, and −0.25 to 0.25 m, respectively. By using these complete and detailed near-field geodetic constraints, we invert the fault structure and slip behaviour based on the finite-fault source modelling. Our results indicate that the 2015 event ruptured a three-segment structure with dip angles of 87°, 70° and 87° along the SKF from south to north, and the 2023 event ruptured an unmapped single fault with a dip angle of 76.8° in the southern Pamir. Note that both of these two events exhibit steep shallow left-lateral strike-slip motions along the NNE-SSW striking direction, accompanied with obvious surface ruptures, but with opposite dip trends. The causative faults of these two events indicate that the NNE-SSW striking direction strike-slip activities in the Pamir hinterland can not only occur on the preexisted SKF in northern and central Pamir, but also continuing propagate southward, suggesting an N-S direction shear transect across the Pamir regime. In addition, the 2023 event should not be triggered by the 2015 event based on the static Coulomb stress changes analysis in this study. By analysing the interseismic GPS measurements, we evaluate the strike-slip activities of the SKF in the Pamir hinterland. These GPS profiles indicate that the SKF exhibits a steep stepwise relative motion of 3.4–5.7 mm yr⁻¹ in the N-S component and a relative motion of 3.2–3.8 mm yr⁻¹ in the E-W component, respectively, implying ongoing compression orogeny and interior varied shortening rate in the Pamir hinterland.

ACKNOWLEDGEMENTS

The Sentinel-1 SAR (https://search.asf.alaska.edu/) data were provided by the ESA through their open data policy. Most Figures were plotted using Generic Mapping Tools (GMT). The open AI tools of ChatGPT 3.5 and AIbox of NetEase Youdao dictionary were used to check possible grammatical errors. This work is supported by the National Natural Science Foundation of China (42174004, 41974004 and 41704005). PH is an academic visitor at the University of Leeds, sponsored by the State Scholarship Fund from the China Scholarship Council..

DATA AVAILABILITY

InSAR data can be downloaded from the website of Zenodo (https://doi.org/10.5281/zenodo.8155679). GPS data are available in the flowing publications. The data sets for the Pamir are from Kreemer et al. (2014, https://doi.org/10.1002/2014GC005407), Zhou et al. (2016, https://doi.org/10.1016/j.gr.2016.03.011) and Perry et al. (2019, https://doi.org/10.1029/2018GL080065). The data for the Tien Shan are mainly from Li et al. (2022, https://doi.org/10.1029/2022GL099105). The data for the Tibet is mainly from Wang & Shen (2020, https://doi.org/10.1029/2019JB018774).