SEISMOGENIC ANALYSIS OF THE 2023 MW6.9 MURGHOB EARTHQUAKE, TAJIKISTAN, BASED ON MULTI-TRACK INSAR AND MODELLING

2.1. InSAR Processing

We utilised Sentinel Application Platform (SNAP), an open-source tool to process Sentinel-1 images, downloadable from the European Space Agency's (ESA) official website. One of SNAP's key advantages is its user-friendly graphical user interface (GUI), simplifying image processing. The processing stages include co-registration, interferogram formation, phase unwrapping, and terrain correction. Co-registration is a crucial step that aligns images by utilizing statistical values to achieve sub-pixel accuracy in spatial matching between them. This process consists of two main steps: the Apply Orbit Information step, which employs precise orbit files, and the Back Geocoding step, which uses the Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) with a 1 arc-second resolution.

An interferogram is generated by multiplying the reference image with the complex conjugate of the secondary image, where the amplitudes of both images are combined, and the resulting phase reflects the phase difference between the two images. The interferometric phase at each pixel of the SAR image is determined by the difference in the travel paths of the radar signals between the two SAR acquisitions to the resolution cell under consideration. The calculated interferogram exhibits phase variations ([TeX:] \mathrm{\Delta}\varnothing insar$\mathrm{\Delta }⌀insar$) due to several contributing factors, including flat Earth phase ([TeX:] \mathrm{\Delta}\varnothing orb)$\mathrm{\Delta }⌀orb)$, topographic phase ([TeX:] \mathrm{\Delta}\varnothing$\mathrm{\Delta }⌀$topo), atmospheric phase ([TeX:] \mathrm{\Delta}\varnothing atm$\mathrm{\Delta }⌀atm$), phase noise ([TeX:] \mathrm{\Delta}\varnothing noise$\mathrm{\Delta }⌀noise$), and surface deformation ([TeX:] \mathrm{\Delta}\varnothing defo$\mathrm{\Delta }⌀defo$). These components are represented mathematically, as shown in Equation 1.

[TeX:] \mathrm{\Delta}\varnothing insar = \ \mathrm{\Delta}\varnothing orb + \ \mathrm{\Delta}\varnothing topo + \ \mathrm{\Delta}\varnothing atm + \ \mathrm{\Delta}\varnothing noise\ + \ \mathrm{\Delta}\varnothing\ defo$\mathrm{\Delta }⌀insar=\mathrm{\Delta }⌀orb+\mathrm{\Delta }⌀topo+\mathrm{\Delta }⌀atm+\mathrm{\Delta }⌀noise+\mathrm{\Delta }⌀defo$ (1)

Factors other than surface deformation (∆∅defo) are eliminated to extract the displacement or vertical change value. The process begins with subtracting the Flat Earth Phase, which is calculated using orbit metadata and removed from the interferometric phase. Subsequently, the Topographic Phase is eliminated by differencing the interferogram data with the Digital Elevation Model (DEM). The atmospheric phase, comprising tropospheric and ionospheric components (Furuya, et al., 2017; Setiawan & Furuya, 2021; Setiawan & Furuya, 2024), is considered. However, for certain cases, such as large earthquakes, the phase changes due to atmospheric effects are significantly smaller than the deformation signal (He et al., 2023). To further reduce noise, multilooking and filtering are applied. Goldstein Phase Filtering (Goldstein & Werner, 1998) filters the data. Since the phase in interferometry is inherently ambiguous and confined to the range of 2π, a phase unwrapping process is necessary to connect the interferometric phase with the actual topographic height. Additionally, a multi-looking process is performed after interferogram formation to optimise phase unwrapping. The multi-looking parameters used are 6 for range and 2 for azimuth. The Statistical-Cost, Network-Flow Algorithm for Phase Unwrapping (SNAPHU) (Chen & Zebker, 2001) is utilised for phase unwrapping. After the unwrapping process, Phase-to-Displacement Conversion and terrain correction are performed. The Phase-to-Displacement Conversion transforms the image data, measured initially in radians, into metric units (Equation 2). Meanwhile, the Terrain Correction process involves geocoding the image by correcting geometric distortions using the Digital Elevation Model (DEM), resulting in a georeferenced product.

[TeX:] d = \ - \frac{\lambda\mathrm{\Delta}\varnothing_{d}}{4}$d=-\frac{\lambda \mathrm{\Delta }{⌀}_d}{4}$ (2)

Where d is the surface displacement, λ is the wavelength, [TeX:] \nu$\nu$ is Poisson’s ratio, [TeX:] \varnothing_{d}${⌀}_d$is the phase different.

2.2. Earthquake Modelling

Earthquake modelling in this study is conducted using Pyrocko (https://git.pyrocko.org/pyrocko), an open-source software designed for applications in seismology and geophysics. Built on the Uniplexed Information Computing System (UNIX), Pyrocko offers diverse toolboxes for analyzing and processing seismic and geophysical data. It is particularly suited for studying Earth's movements, earthquakes, and geodynamics (Heimann et al., 2017).

Kite, one of the toolboxes available in Pyrocko, can be used to process data from the Sentinel Application Platform (SNAP). Specifically, Kite is designed to process and analyze surface displacement data derived from InSAR (Isken et al., 2017). Within Kite, a data down-sampling process is performed to facilitate modelling. The quadtree sampling methodology efficiently manages high-resolution Interferometric Synthetic Aperture Radar (InSAR) data. This approach subdivides the dataset into a hierarchical arrangement of tiles based on spatial variability, as illustrated in Figure 3. The epsilon parameter modulates the sensitivity of tile refinement, thereby ensuring that regions exhibiting significant variability are sampled at elevated resolutions, while uniform areas are represented by a smaller number of larger tiles. Parameters such as minimum tile size and maximum tile size establish the constraints for tile dimensions, thereby achieving a balance between resolution and computational efficiency. This adaptive technique facilitates effective data processing while maintaining critical spatial details, particularly in regions characterized by complex surface displacement patterns.

Figure 3. Downsampling of InSAR data using the Kite module. The figure illustrates how interferogram data points are reduced to representative subsets, allowing efficient modelling while preserving the main deformation signal.

Earthquake modelling in this study is performed using the Grond (Heimann et al., 2018) toolbox, with data imported from the Kite toolbox. To conduct modelling utilizing Pyrocko, it is essential to acquire two distinct scenes, specifically those characterized by ascending and descending trajectories within the earthquake-affected region. The modelling employs a rectangular finite fault model, requiring the input of various parameters. Eight iterations are conducted, progressively increasing: 250, 500, 750, 1000, 3125, 6250, 12500, and 25000 iterations. The primary input parameters for the model include depth values, east_shift, north_shift, length, width, strike, dip, and rake. Calculations are performed using the Bayesian Bootstrap Optimization (BABO) statistical approach. BABO is a global optimization algorithm designed for multi-objective functions through directed searches. It enables the formation of a set of objective function, identifying multiple minima, irregularly shaped minima, and evaluating parameter trade-offs and uncertainties in terms of probability.

2.3. Multi-Track InSAR Processing

The multi-track InSAR approach combined multiple Sentinel-1 tracks to differentiate surface displacement's horizontal (dH) and vertical (dV) components. The dH and dV maps were derived from ascending and descending track data and thoroughly processed to ensure precise deformation characterizations. This approach facilitates enhanced discernment of horizontal and vertical movements that might remain obscured in the Line-of-Sight (LOS) data alone.

It is important to note that the multi-track technique implemented in this analysis explicitly maps the vertical (dV) and horizontal (dH) displacement dimensions. However, most SAR observation (including Sentinel-1) has a near-polar orbit, so the North-South displacement cannot be observed and produces a negligible phase change (Equation 3) (Fuhrmann & Garthwaite, 2019; Fujiwara et al., 2000). Nonetheless, complete 3D displacement from InSAR observation may be obtained if right-looking and left-looking images are available (Wright et al., 2004), or by combining with Global Navigation Satellite System (GNSS) observation.

[TeX:] \binom{v_{\text{asc~}}}{v_{\text{des~}}} = \begin{pmatrix} - sin\theta_{\text{asc~}}cos\alpha_{\text{asc~}} & cos\theta_{\text{asc~}} \\ - sin\theta_{\text{des~}}cos\alpha_{\text{des~}} & cos\theta_{\text{des~}} \end{pmatrix}\binom{v_{E}}{v_{U}}\ \ $\left(\genfrac{}{}{0ex}{}{v_{\mathrm{\text{asc }}}}{v_{\mathrm{\text{des }}}}\right)=\left(\begin{array}{cc}\hfill -sin{\theta }_{\mathrm{\text{asc }}}cos{\alpha }_{\mathrm{\text{asc }}}\hfill & \hfill cos{\theta }_{\mathrm{\text{asc }}}\hfill \\ \hfill -sin{\theta }_{\mathrm{\text{des }}}cos{\alpha }_{\mathrm{\text{des }}}\hfill & \hfill cos{\theta }_{\mathrm{\text{des }}}\hfill \end{array}\right)\left(\genfrac{}{}{0ex}{}{v_E}{v_U}\right)$ (3)

3.1. InSAR Result

The interferometric process using Synthetic Aperture Radar (SAR) images produces a spectral pattern known as fringes, where each concentric circle represents half the wavelength of the radar signal. These fringes were converted into unwrapped phase imagery using the SNAPHU algorithm and subsequently transformed into displacement measurements in meters. This process provides a detailed depiction of surface deformation within the research area.

Analysis of the ascending track data (see Figure 4a) reveals displacement values ranging from -0.21 to -0.04 meters. Similarly, the descending track data (see Figure 4b) indicates displacement values ranging from -0.13 to 0.02 meters. Particularly, the area surrounding the earthquake's epicenter exhibits a distinct pattern, with some regions showing significant subsidence while others experience uplift. This phenomenon is attributed to active fault movement in the region, particularly during the Murghob earthquake, highlighting the dynamic tectonic activity shaping the Earth's surface in this area.

Figure 4. Unwrapped interferograms reveal coseismic deformation from the Murghob earthquake: (a) ascending and (b) descending Sentinel-1 tracks. Changes in colour gradients point to displacements in the line of sight, predominantly with maximum alterations occurring near the epicentral region. Negative measurements reveal ground displacement moving away from the satellite (subsidence), while positive measurements reflect an upward shift (uplift).

3.2. Earthquake Modelling Result

We used the standard epsilon value of 0.04 m in the quadtree meshing method. This value was chosen to preserve minor deformation signals while preventing excessive patch subdivision from noise. A lower epsilon value would increase the number of quadtree patches, potentially introducing noise, whereas a higher value could smooth deformation features excessively, reducing spatial resolution.

We adjusted the prior distributions for source parameters (strike, dip, rake, depth, and slip) from their default settings for the Bayesian inversion process. These adjustments were informed by data from the USGS and GCMT assessments of the 2023 Mw 6.9 Murghob earthquake and relevant regional tectonic factors. This use of informative priors ensured that the search space accurately represented credible fault geometries for the research area. Although we did not perform a sensitivity analysis, the posterior results matched catalog solutions, suggesting that the selected priors and the standard epsilon value did not significantly bias the inversion results.

The modelling stage using Grond produces a series of images (see Figure 5 and Figure 6) critical to the research analysis. These include the observed, modelled, and residual images, each serving a specific purpose in understanding the fault dynamics. The observed image represents the raw data, visualizing all available information without modifications or processing. This image serves as a baseline, capturing the unaltered characteristics of the earthquake's surface deformation. The modelled image depicts the data after being processed through the modelling software. It provides a refined representation of the fault system by incorporating the modelled parameters, offering insight into the earthquake's underlying mechanisms. Lastly, the residual image highlights the discrepancies between the observed and modelled data. It represents the portions of the observed data that the model could not fully explain, pointing to areas where additional complexities or variations remain unaccounted for in the current modelling framework.

The iterative process is critical in producing an effective model. We performed the iteration process eight times, adjusting the iteration count parameter at each step. This iterative approach is vital because it allows the model to be refined through repetition, improving the accuracy of the model's representation of the observed phenomenon. We decided to halt the iterations after reaching 25,000 iterations, as the modelling results at that point were deemed sufficient and aligned with our research objectives. In addition to generating the images, the modelling process yields the best-fit parameters. Our result displayed a calculated fault strike of 128°, with a dip of 67° and a rake of 177°.

Figure 5. Modelling results for the ascending track: (A) observed displacement from InSAR, (B) best-fit synthetic model, and (C) residuals (observed–modelled). The residuals are significantly reduced, indicating the model adequately reproduces the observed deformation.

Figure 6. Modelling results for the descending track: (A) observed displacement from InSAR, (B) best fit synthetic model, and (C) residuals (observed–modelled). The residuals indicate a strong correlation between observations and model prediction, akin to the ascending track

3.3. Multi-Track InSAR Result

Figure 7a and Figure 7b represent the horizontal (dH) and vertical (dV) displacements, respectively, derived from multi-track InSAR analysis of the 2023 Murghob earthquake in Tajikistan, offering invaluable insight into the surface deformation caused by this significant seismic event. Figure 7 shows horizontal displacement (dH), with values ranging from -0.53 to 0.12 meters, while Figure 7 depicts vertical displacement (dV) with values between -0.13 and 0.32 meters. Both displacement fields were generated by combining multiple Sentinel-1 radar tracks, allowing for precise separation of horizontal and vertical movement components that would otherwise be hidden in line-of-sight (LOS) data alone. The ability to isolate these movements is crucial for understanding the full extent of tectonic shifts during the earthquake.

Figure 7. Studying the multi-track InSAR results uncovers (A) east to west horizontal movement and (B) vertical changes in position. The horizontal dimension points to substantial east–west ground activity in the fault sector, whereas the vertical dimension stresses uplift and subsidence dynamics near the epicenter. The data corresponds with strike-slip faulting tied to the seismic phenomenon.

4.1. The essential of Digital Elevation Model (DEM) selection in InSAR processing

Utilizing an external DEM to correct topographic phase is crucial in InSAR processing (Hanssen, 2002; Furuya, 2011). A higher spatial resolution DEM can be beneficial in rugged terrain (such as in our study area), since DEM height errors cause residual fringes in the interferogram (Lazecký et al., 2015). The high variable topography may be smoothed out in a coarse DEM, which later can cause unwrapping error. However, in relatively flat areas, the difference between using 30 m or 90 m DEM in InSAR processing could be minimal, since the elevation-induced phase is small (Du et al., 2017).

4.2. Fault Parameters and Tectonic Complexity of the 2023 Murghob Earthquake

The fault parameters obtained in our study show some differences compared to the results of Chen et al. (2024); Shi et al. (2023), although all studies agree on the dominant right-lateral strike-slip event of the earthquake. In our study, the calculated fault strike was 128°, with a dip of 67° and a rake of 177°, reflecting a high dip, almost vertical fault geometry. On the other hand, Chen et al. (2024) observed a strike of 124°, a dip of 86°, and a rake of -147.35°, while Shi et al. (2023) reported a strike of 28.08°, a dip of 72.89°, and a rake of -1.85° (see Table 2). These variations in fault parameters may arise due to differences in the data sources, modelling techniques, and assumptions used by each study.

While our approach integrates 30 m resolution Digital Elevation Models (DEM) and both ascending and descending track InSAR data to achieve detailed spatial resolution of surface deformation, we do not aim to claim methodological superiority over prior studies by institutions such as GFZ (Germany), IPGP (France), USGS, or the Global CMT (gCMT) project. Instead, our study offers a complementary perspective—focusing on the near-surface deformation field and possible complex rupture behaviour. InSAR-based modelling is particularly sensitive to surface displacements and provides a fine-scale view of fault slip distribution, while seismic-based methods such as waveform or moment tensor inversion (e.g. USGS, gCMT) are more effective in characterizing source parameters at depth. Thus, both approaches offer valuable and complementary insight into the rupture process.

In terms of fault geometry, our results align more closely with the models proposed by Chen et al. (2024) and GFZ, especially regarding strike and dip values. This alignment may be attributed to the shared reliance on geodetic data and similar inversion methodologies. In contrast, the USGS and gCMT models suggest different orientations and rake angles, possibly generalizing the event as a pure strike-slip rupture. These differences highlight how methodological approaches—InSAR versus seismic inversion—can influence fault parameter interpretations across studies.

Table 2. Parameters comparison with our study

The differences in fault geometry and slip characteristics further underline the complexity of the earthquake's seismogenic structure. While prior studies (Chen et al., 2024; Shi et al., 2023) focus on the dominant strike-slip component with minor variations in dip and rake angles, our model reveals additional complexity by identifying normal faulting contributions in the northern rupture zone, particularly along the Sarez-Karakul Fault system and the Sarez-Murghab Thrust system. This mixed-slip behaviour, which is less emphasized in other studies, may reflect a more intricate rupture process. Analogous patterns of mixed-mode rupture have been documented in other tectonically complex settings, such as the 2008 Wenchuan (Lei & Zhao, 2009) and 2011 Tohoku earthquakes (Yagi & Fukahata, 2011), where rupture propagation involved multiple fault systems and mechanisms.

The use of InSAR in our study allows direct measurement of surface deformation at high spatial resolution, making it well suited to identify distributed deformation and possible secondary faulting. However, we recognize the importance of combining this geodetic information with seismic data to better capture the full rupture process from surface to depth.

All three studies—our own, Chen et al. (2024), and Shi et al. (2023)—recognize the involvement of multiple fault zones in the rupture. Prior studies. emphasize the Sarez-Karakul and Aksu Murghab fault systems as major contributors (Chen et al., 2024), while Shi et al. (2023) focus more heavily on the Sarez-Karakul zone. Our findings are consistent with these, but introduce added complexity by suggesting the involvement of a blind right-lateral strike-slip fault along the eastern segment of the Aksu Murghab Fault system. Additionally, we observe indications of normal and minor left-lateral strike-slip components in the northern rupture area, which were not as clearly detected in other models. These differences further emphasize the tectonic intricacy of the Central Asia collision zone, where oblique convergence between the Indian and Eurasian plates leads to varied faulting styles (Avouac & Tapponnier, 1993; Singh & Ghosh, 2019; Styron et al., 2011; Tapponnier et al., 2001; Tapponnier, et al., 1981).

Previous seismic events may also help explain the rupture behaviour of the 2023 earthquake. The 2015 Mw7.2 Murghab earthquake occurred near the epicentral area and involved left-lateral strike-slip and normal faulting in the northern segment (Sangha et al., 2017). It likely altered the local stress field, with Coulomb stress transfer potentially promoting failure along the Sarez-Karakul Fault system and other regional structures in 2023 (Li et al., 2021; Li, et al., 2020; Wang et al., 2017). Similar triggering effects have been observed in other large strike-slip earthquakes, such as the 1999 Izmit (Reilinger et al., 2000) and 1992 Landers events (Massonnet et al., 1993).

In conclusion, all major studies agree that the 2023 Tajikistan earthquake was predominantly a right-lateral strike-slip event. However, there are essential differences in fault parameters, geometries, and inferred fault system contributions. Our model, which incorporates additional normal faulting and thrust components, adds new dimensions to understanding this rupture. These differences reflect the complexity of the regional tectonics and the complementary nature of diverse modelling techniques. These differences reflect not only the complexity of the regional tectonics but also the complementary nature of diverse modelling techniques. We emphasize that our results are not meant to supersede existing models, but rather to enhance the seismogenic interpretation of this event. We advocate for the continued integration of high-resolution InSAR, seismic data, and tectonic history to develop more complete and accurate assessments of earthquake behaviour in tectonically complex regions like the Pamir and Hindu Kush (Ischuk et al., 2013; Schurr et al., 2014; Sippl et al., 2013).

4.3. Insights into Fault Dynamics and Tectonic Complexity from Multi-Track InSAR Displacements

Vertical displacement can be derived from a single pair of InSAR based on Line of Sight (LoS) and satellite incidence angles (Cigna & Tapete, 2022). Nonetheless, the conversion must be carefully taken and assume the movement is only in the vertical direction. Thus, any movement in the horizontal component can have huge inaccuracies (Fuhrmann & Garthwaite, 2019). In the region with two SAR tracks (ascending and descending) available, pseudo 3D (2.5 D) can be used (Fujiwara et al., 2000).

The vertical displacement (dV) image from multi-track InSAR highlights a predominant upward movement in certain areas, with the positive displacement (green) being more dominant than the negative displacement (red) (see Figure 7). This vertical movement suggests the occurrence of uplift, which is consistent with normal faulting or thrust faulting processes that often accompany lateral fault movement in compressional settings. The Sarez-Karakul normal fault and other thrust and strike-slip fault components in the region likely contributed to these vertical shifts.

On the other hand, the horizontal displacement (dH) image reveals significant westward and eastward motions, with red zones indicating negative displacements (westward movement) and green zones showing positive displacements (eastward movement) (see Figure 7). This pattern suggests that the earthquake involved lateral fault movements, typical of strike-slip faulting, where two blocks of the Earth's crust move horizontally past each other. The distribution of displacements in these directions aligns with the region's fault system, specifically the Aksu Murghab Fault system, known for its strike-slip characteristics.

The tectonic complexity of the Tajikistan region is further underscored by the various fault systems that interact in the area. Studies such as those by Chen et al. (2024; Shi et al. (2023) have highlighted the importance of the Sarez-Karakul Fault system and the Aksu Murghab Fault system as primary contributors to the 2023 earthquake. However, our study adds further detail by suggesting that the rupture also likely involved a blind right-lateral strike-slip fault, primarily driven by the Aksu Murghab Fault system in the east. Additionally, we observe influences of normal faulting and left-lateral strike-slip faulting in the northern part of the region, which were less emphasized in previous studies.

The region’s tectonic setting is further complicated by its position within the broader Central Asia collision zone, where the Indian and Eurasian plates interact. Various faulting styles characterize this collision zone, including oblique collision, strike-slip faulting, and thrust faulting. The role of these different fault types in shaping the region’s tectonic landscape has been well-documented in numerous studies (Avouac & Tapponnier, 1993; Singh & Ghosh, 2019). The interactions between these fault systems contribute to the complex seismic activity observed in Tajikistan and surrounding regions. Understanding the contributions of various fault zones to earthquake rupture is essential for improving seismic hazard assessments and developing more effective disaster mitigation strategies in the region.

Additionally, the analysis of horizontal and vertical displacements through multi-track InSAR provides a detailed understanding of the seismic deformations caused by the 2023 Murghob earthquake. The horizontal displacements reveal significant strike-slip faulting, while the vertical displacements suggest a complex combination of uplift and normal faulting. These findings underscore the importance of further research into the fault systems in Tajikistan to better understand the tectonic forces shaping the region and to improve earthquake preparedness and risk assessment.