Coseismic Source Model of the February 2023 Mw 6.8 Tajikistan Earthquake from Sentinel-1A InSAR Observations and Its Associated Earthquake Hazard

Abstract

On 23 February 2023, an Mw 6.8 earthquake struck the border of Tajikistan and Xinjiang China, the source mechanism remains controversial according to different seismic inversions. To better comprehend the source characteristics and the surface deformation pattern, we used the ascending and descending orbital Sentinel-1A SAR data to obtain the coseismic deformation of this earthquake based on the traditional two-pass differential interferometric synthetic aperture radar (InSAR). The source model is inverted from the InSAR coseismic deformation results. The possible Coulomb Failure Stress (CFS) transfer is analyzed based on the preferred source model. The results illustrate that the earthquake ruptured a blind left-lateral strike-slip fault of strike 28.1° with a maximum slip of 1.53 m and the total geodetic moment is 1.99 × 10¹⁹ N·m (Mw 6.83). The strike direction and the fault characteristics suggest the Seismogenic fault is a secondary fault of the Sarez–Karakul Fault System. The 2015 Mw 7.2 Sarez Earthquake plays a triggering role in the occurrence of the 2023 Tajikistan earthquake. Earthquake hazard on Sarez–Karakul Fault System and Sarez–Murghab Thrust System is enhanced due to the Coulomb stress loaded by the Tajikistan earthquake.

1. Introduction

According to the officially released results of the China Earthquake Network Center (CENC), on 23 February 2023, an M 7.2 earthquake struck eastern Tajikistan, which is only approximately 82 km away from the China–Tajikistan border (Figure 1). The results show that the epicenter of the Tajikistan earthquake was located at (73.29°E, 37.98°N), the event initiated at 08:37:39 Beijing time, and the depth is 10 km. Different institutes or agencies released the moment tensor solutions of this earthquake, however, there are noteworthy discrepancies between them (

Table 1. Source mechanism solutions of the Tajikistan earthquake from different organizations.

Source	Epicenter		Depth (km)	Fault Plane			Mw
	Longitude (°)	Latitude (°)		Strike (°)	Dip (°)	Rake (°)
USGS	73.23	38.06	13.5	203/305	57/72	−21/−146	6.87
gCMT	73.22	38.15	16.9	210/120	88/89	1/178	6.8
GFZ	73.29	38.06	10.0	122/25	73/67	−156/−18	6.8
IPGP	73.208	38.073	8.0	27/123	73/71	−20/−162	6.81
CENC	73.29	37.98	12.0	313/188	52/53	−131/−50	6.8
This study a	73.206	38.153	3.67	28.08	72.89	−1.85	6.83

Note: USGS—U.S. Geological Survey. gCMT—the Global Centroid-Moment-Tensor (CMT). GFZ—GFZ German Research Centre for Geosciences. IPGP—Institut de physique du globe de Paris. CENC—China Earthquake Network Center. ^a The shear modulus used is 33 GPa. The epicenter is the fault top center of the non-linear inversion result.

). The epicenter locations determined by them are similar, conversely, the fault depths, and the other parameters are different especially the fault strike direction. The centroid moment tensor solutions of the CENC show that the first fault planar is NW directional (https://data.earthquake.cn/datashare/report.shtml?PAGEID=earthquake_dzzyjz, accessed on 20 March 2023). The W-phase moment tensor solutions released by the U.S. Geological Survey (USGS, https://www.usgs.gov/programs/earthquake-hazards/earthquakes, accessed on 23 February 2023) and the Global Centroid Moment Tensor Project (gCMT, https://www.globalcmt.org, accessed on 23 February 2023) show that the fault strike direction of the first fault planar is the SSW. Differently, the fault directions of the GFZ German Research Centre for Geosciences (GFZ, https://geofon.gfz-potsdam.de/, accessed on 23 February 2023) and the Institut de physique du globe de Paris (IPGP, http://geoscope.ipgp.fr/index.php/en, accessed on 23 February 2023) are SE and NE. Moreover, the slipping rakes of these organizations are different too.

Geologically, the area where the Tajikistan earthquake occurred belongs to the Western Himalayan Syntaxis, the tectonic setting is like the Eastern Himalayan Syntaxis and is very complex. The collision and the continuous extrusion between the Indian plate and the Eurasian plate started 50 Ma ago and the Indian plate subducts with approximately 38 mm/yr beneath the Eurasian plate [1]. The collision and the extrusion of these two plates caused crustal continuous deformation and many faults in this area [2,3,4]. Interestingly, the properties of these faults are distinct, for instance, the Sarez–Karakul fault is sinistral predominantly, the Aksu–Murghab Fault is a dextral fault, the Tanymas and the Sarez–Murghab faults are thrusting, the small faults at the south of the Sarez–Karakul fault are normal [5]. The Global Strain Rate Map [6] shows that the strain rate of the Western Himalayan Syntaxis is high, the second invariant strain rate of the Trans Alai Range between the South Tian Shan and the Pamir is the highest of the Western Himalayan Syntaxis. The strain rate of the Tajikistan earthquake epicenter zone is not high before the earthquake, the second invariant strain rate is only 12 nano-strain/yr. Many devastating earthquakes occurred along the junction zone between the South Tian Shan and the Pamir plateau [7], however, the seismicity of the Sarez–Murghab fault zone is less active. The 2015 Mw 7.2 Sarez earthquake is the only strong earthquake that happened during the last 100 years on the Sarez–Karakul fault, the distance between the 2023 earthquake and the 2015 earthquake is only approximately 40 km [8]. Based on the above, there are two issues still not proved yet, whether the 2023 Tajikistan earthquake is linked to the 2015 Mw 7.2 Sarez earthquake, and which fault system the 2023 Tajikistan seismogenic fault belongs to.

Interferometric synthetic aperture radar (InSAR) has been widely used in mapping coseismic deformation to reflect seismic source characteristics [9,10,11,12,13]. It can provide near-field constraints on the determination of the source model. In this study, to understand the coseismic deformation and the seismogenic characteristics associated with the 2023 Tajikistan earthquake and its implications for the earthquake hazard, we derive the SAR line-of-sight (LOS) coseismic deformation using the traditional two-pass differential InSAR technique with ascending and descending orbital Sentinel-1A SAR data. We carry out non-linear and linear inversions to determine the source slip model based on the coseismic InSAR deformation. We further analyze the CFS transfer due to the 2015 and 2023 earthquakes and the tectonic implications.

2. Data and Methods

2.1. SAR Data and Processing Strategy

The Sentinel-1A satellite, which is operated by the European Space Agency images the Earth’s Surface for the purpose of disaster mitigation all day long and its repeated interval is 12 days. The synthetic aperture radar sensor carried on the Sentinel-1A satellite is the C-band and the SAR data imaged by the mode of the interference wide mode (IW) is used for the InSAR processing to obtain surface deformation. Before and after the Tajikistan earthquake, the Sentinel-1A satellite imaged the epicenter area with ascending and descending orbital directions. To obtain InSAR coseismic deformation with high-quality coherence of SAR imagery pairs, shorter temporal and spatial baselines are preferred [14,15]. We selected two SAR imagery pairs in this study and list their information in

Table 2. Detailed information of SAR data.

Sensor	Acquisition Time (D-M-Y)	Orbital Direction	Path Number	Frame Number	Heading Angle (°)
Sentinel-1A	22 January 2023	Ascending	100	117–122	−13.17
	27 February 2023
	21 February 2023	Descending	5	461–466	−166.60
	5 March 2023

, the polarization mode is VV.

The GMTSAR [16,17] is open-source software and is very friendly to scientific researchers. The InSAR processing can be implemented step by step or batch at once. In this study, we absorbed the traditional two-pass differential InSAR strategy and obtained the coseismic deformation based on it. During the InSAR processing, we employed a bigger multi-look factor of 20:4 to increase the quality of the InSAR results. To eliminate the topographic phase and the flat phase, we used the radar topography mission digital elevation model (SRTM DEM) [18] with 90 m resolution as the external data. We also referred to the SRTM DEM when we co-registered the secondary image to the reference one and geocoded the InSAR results. A 400 m wavelength spatial Gaussian filter was applied before unwrapping the interferograms. We masked out the area where the coherence was less than 0.05 in the interferograms and unwrapped them with the statistical-cost network-flow algorithm (SNAPHU) [19] integrated in the GMTSAR. The ascending and descending orbital deformation interferograms are shown in Figure 2. Each fringe in the interferograms corresponds to 10 cm ground displacement along the LOS.

Further, we downsampled the ascending and descending orbital deformation fields with the quadtree algorithm [20] for the source modeling inversion. The variance threshold for quadtree division is 3.14, and the minimum and maximum sizes for each quadtree block are 8 and 128, respectively. The observation numbers are 1368 and 1098 for the ascending and descending orbital deformation fields after down-sampling.

2.2. Source Model Inversion

According to the relationship between the surface displacements and the seismogenic fault motion in a homogeneous elastic half-space [21], the fault model parameters are non-linear to the observations. To obtain the source model, the traditional method converts the above non-linear problem into two steps [20,22], one is the non-linear inversion for the fault geometric parameters, and the second is the linear inversion for the slip distribution. Therefore, we employed this strategy to carry out the source model inversion. Considering the observations were from the same type of data and the same processing flow, we applied equal weight for the observation constraint during inversions.

In recent years, the Bayesian statistical model has been widely used in many geophysical fields to evaluate source models. In the first step, to perform the Bayesian estimation of the fault geometric parameters including the fault top center position, top center depth, strike, dip angle, length, width, and slip-along strike and dip directions, we imported an open-source, stable package emcee [23], which is the Python implementation of the affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) [24]. The likelihood function is constructed based on the homogeneous elastic half-space dislocation model [21]. Assuming a uniform prior probability distribution for each fault parameter, we set a wide range for each fault parameter considering the notable discrepancy between previous source solutions (Table 1). In the definition of fault geometry, the global ranges of the strike direction and the dip angle should be 0.0° ≤ strike ≤ 360.0° and 0° ≤ dip ≤ 90.0°. According to the coseismic deformation, the seismogenic fault of the Tajikistan earthquake is steep. Therefore, we set the range of dip angle was set as [60°, 90°]. For programming convenience, we set the dip range as [60.0°, 120.0°] and the strike angle as [0°, 180.0°]. If the dip is bigger than 90°, we subtract the dip angle from 180°and add 180° to the strike angle to obtain the new dip angle and the new strike angle. This processing can guarantee the strike direction range of [0°, 360°]. Furthermore, we chose the prior ranges for other geometric parameters of the seismogenic fault in the following manner: length ∈ [15.0 km, 50.0 km], width ∈ [4.0 km, 20.0 km], depth ∈ [0.0 km, 10.0 km], strike-slip ∈ [−3.0 m, 3.0 m], dip-slip ∈ [−3.0 m, 3.0m], X ∈ [320.0 km, 362.0 km], and Y ∈ [4200.0 km, 4236.0 km]. We plotted the two-dimensional posterior distribution and one-dimensional histogram of all the parameters in Figure 3. For each fault geometric parameter, we obtained posterior probability density functions and determined the mean of the posterior solution as each best-fitting parameter.

We used an indicator called goodness of fit to judge the inverted parameters and the source model, which is defined as

$$R=1-\frac{\sum {\left(d_{obs}-d_{sim}\right)}^2}{\sum d_{obs}{}^2}$$

(1)

in which d_obs and d_sim represent the vectors of observations and model simulations, respectively. Generally, the model can be reasonable when the goodness of fit is greater than 0.8.

In the second step, to detail the slip on the fault plane, we fixed the fault spatial position, strike direction, and dip angle. Additionally, the mechanism of the seismogenic fault inverted in the first step was used as a priori constraint in the inversion. Based on the fault geometric parameters, we construct a bigger-size fault model with a length of 60 km and a width of 30 km. Then, we discretized the single fault into 450 sub-faults each size 2 km × 2 km along the strike and dip directions. The relationship between the InSAR deformation observations and the sub-faults slip as follows:

$$d_{obs}=G·s$$

(2)

where d_obs are the InSAR deformation observations, G is the Green’s function, and s is the sub-faults slip vector. Green’s function matrix for predicting the ground line-of-sight (LOS) displacement was calculated by assuming the unit slip on each sub-fault based on the elastic dislocation model [21]. To avoid the influence of the ill-conditioned issue in the equation solving of (2), we applied the constrained least-squares method in the obtainment of the source slip model with the smoothing constraint on the slips between neighbor sub-faults, the zero constraints for boundary sub-faults, and positive constrain for the strike-slip for each sub-fault.

3. Results

3.1. Coseismic Surface Deformation

The coseismic deformation of the ascending and descending orbital LOS directions is shown in Figure 4. From Figure 2 and Figure 4, we found that the signal-to-noise ratio is worse than other InSAR results in typical areas, the main reason is that the coherence is low in most parts of the epicentral area due to the ground vegetation coverage. Nevertheless, the fringes in Figure 2 for each interferogram are visible, and it is apparent from Figure 2 and Figure 4 that there are three main deformation zones in the epicentral area of the Tajikistan earthquake. The directions of surface displacements in ascending and descending are opposite each other. According to the geometry of SAR imaging [25], the surface rise or subsidence appears in the same directions in ascending and descending orbital interferograms. On the contrary, the horizontal deformation shows opposite directions. As we know that the strike-slip motion of a fault caused horizontal deformation, we can conclude that the seismogenic fault of the Tajikistan earthquake is strike-slip faulting. From Figure 4, for the ascending orbital deformation field, both the maximum surface displacements of LOS away from and toward the satellite are 0.12 m. For the descending orbital deformation field, the maxima of the surface displacement are 0.15 m and 0.14 m in the LOS away from and toward the satellite.

3.2. Source Model

With the non-linear and linear inversions, we obtained the source model which can recover the downsampled observations well. For the first step, the residual root-mean-square (RMS) between the observation and the simulated displacement is 2.1 cm, the goodness of fit is 0.79 slightly less than, but very close to, the threshold of 0.8. Therefore, we think the geometric parameters of the fault are reasonable. In the second inversion step, the fit degree of the inversion to the observations was improved, the goodness of fit had increased to 0.84, and the residual RMS decreased to 1.9 cm. To check the reasonability of the source model to the ascending and descending orbital deformation, we simulated the deformation fields in two orbital LOS directions with the preferred source model and calculated the residuals between the observed and simulated results. The results, as seen in Figure 4, indicate that our preferred source model can explain the observations well.

The preferred source slip distribution model is shown in Figure 5. The non-linear inversion result is consistent with the first nodal plane result of IPGP and the second nodal plane result of GFZ (Table 1), which indicates that the seismogenic fault of the Tajikistan earthquake is an NNE strike fault and dipping to the SE direction steeply. The source model result (Figure 5) also reveals that the seismogenic fault of the Tajikistan earthquake is blind, located at the southern end of the Sarez–Karakul Fault system, and almost pure left lateral. It is apparent from the source model that the Tajikistan ruptured one asperity with a maximum slip of 1.53 m, and the main slip zone is located within the depths of 5–24 km. Assuming the shear modulus of 33 GPa, the geodetic moment summation of our preferred slip model is 1.99 × 10¹⁹ N·m, corresponding to Mw 6.83, which is very close to most of the previously mentioned moment magnitudes (Table 1).

4. Discussion

4.1. Source Model Reliability and Its Relation to the Sarez-Karakul Fault

The checkerboard test is commonly used to verify the reliability of the inverted model in geophysical research [15]. Therefore, the checkerboard test was carried out to evaluate the preferred source model and the resolution of the observations. We constructed a fault with the same geometric parameters and discretized the fault with the same size sub-faults as the preferred source model. Meanwhile, we pre-defined an asperity that ruptured 10 × 6 sub-faults with 1 m pure left-lateral slip in the constructed fault. Then, the coseismic LOS deformation on the two orbital observation points was simulated with the constructed input source model. In the procedure of the checkerboard test, the same linear inversion method and smoothing constraint as the above inversion was implemented. The recovered model is shown in Figure 6. It is apparent from Figure 6 that most of the slip zone was recovered, and the central part is the most reliable. Even though the boundary of the slip asperity was not inverted well, the simulated observations are enough to constrain the source model and identify the ruptured zone of the asperity. The checkerboard test result indicates that the preferred source model of the Tajikistan earthquake is reliable.

For the geological issue of the secondary fault, Chinnery [26] analyzed the possible reasons for secondary faults using theoretical stress calculation and proposed that the shear stress on most parts of the fault will decrease except at the ends of the fault after the fault movement. The additional stress loaded on the ends of the main fault is large enough to cause the secondary fault rupture. At the northern and southern ends of the Sarez–Karakul Fault, there are some known secondary faults. On the contrary, there are no secondary faults in the central part area, especially in the cross area with the Tanymas Thrust fault. According to the source model of the Tajikistan earthquake, the seismogenic fault is approximately paralleled to the Sarez–Karakul Fault and only about 20 km away from Sarez–Karakul Fault. Additionally, these two faults have the same left-lateral strike-slip characteristics. The spatial relation pattern of these two faults is very similar to Type A, which is a typical mode of the main fault and its secondary faults proposed by Chinnery [26]. Therefore, we infer that the seismogenic fault of the Tajikistan earthquake is a secondary fault of the Sarez–Karakul Fault.

4.2. CFS Change Associated with the 2023 Tajikistan Earthquake

Previous studies have demonstrated that the static CFS change caused by an earthquake on the nearby faults plays an important role in promoting future strong earthquakes or delaying possible earthquakes [27,28]. The static CFS changes estimation has been widely used in the evaluation of the subsequent strong earthquake possibility over the seismic zones. Before the 2023 Tajikistan earthquake, an Mw 7.2 earthquake happened in 2015 within 40 km. Parsons [29] suggested that the distance of the triggering effect can reach up to 240 km. Therefore, we calculated the CFS change (Figure 7) caused by the 2015 earthquake to test the triggering effects of the 2023 Tajikistan earthquake by the open-source software Coulomb v3.4 [30]. The employed finite fault slip model of the 2015 earthquake is from the USGS. The results in Figure 7 show that the increased CFS changes in the depth of 10 km on the southern part of the seismogenic fault of the 2023 earthquake caused by the 2015 earthquake are larger than 0.1 bar, which is the critical threshold of triggering subsequent strong earthquakes pointed out by Hardebeck et al. [31]. The result of the CFS change demonstrates that the 2015 Mw 7.2 earthquake may have contributed to the triggering of the 2023 Tajikistan earthquake.

To evaluate the potential of future seismicity on the nearby faults after the 2023 Tajikistan earthquake, we calculated the CFS changes with the source model inverted in this study as the input model. According to the possible mechanisms of different faults in Figure 1, we calculated the CFS changes on receiver faults with three types of mechanisms including strike-slip, thrust, and normal faults in three depths of 5 km, 10 km, and 15 km. The results of CFS changes are shown in Figure 8. It is apparent from Figure 8a–c that the CFS change in the southern end of the Sarez–Karakul Fault is positive and exceeded the threshold of 0.1 bar. This implies that the earthquake hazard on the southern end of the Sarez–Karakul Fault is enhanced. For the thrust faults, the CFS only increased on the Sarez–Murghab fault notably and was bigger than 0.1 bar too, which reminds us of the earthquake hazard enhancement. From Figure 8g–i, there is no CFS increase on nearby normal faults.

5. Conclusions

In this study, we present a source model for the 2023 Tajikistan earthquake, which is inverted from the Sentinel-1A ascending and descending orbital InSAR deformation. The model demonstrates that the Tajikistan earthquake ruptured one asperity on a blind and almost pure sinistral strike-slip fault, which strike is 28.1° NNE. The released total moment is up to 1.99 × 10¹⁹ N·m (Mw 6.83). The Seismogenic fault of the Tajikistan earthquake is a secondary fault of the Sarez–Karakul Fault System. The 2015 Mw 7.2 Sarez Earthquake plays a triggering role in the occurrence of the 2023 Tajikistan earthquake. We should pay more attention to the earthquake hazard on Sarez–Karakul Fault System and Sarez–Murghab Thrust System due to their CFS increases.