Deformation and Source Parameters of the 2015 Mw 6.5 Earthquake in Pishan, Western China, from Sentinel-1A and ALOS-2 Data

Abstract

In this study, Interferometric Synthetic Aperture Radar (InSAR) was used to determine the seismogenic fault and slip distribution of the 3 July 2015 Pishan earthquake in the Tarim Basin, western China. We obtained a coseismic deformation map from the ascending and descending Sentinel-1A satellite Terrain Observation with Progressive Scans (TOPS) mode and the ascending Advanced Land Observation Satellite-2 (ALOS-2) satellite Fine mode InSAR data. The maximum ground uplift and subsidence were approximately 13.6 cm and 3.2 cm, respectively. Our InSAR observations associated with focal mechanics indicate that the source fault dips to southwest (SW). Further nonlinear inversions show that the dip angle of the seimogenic fault is approximate 24°, with a strike of 114°, which is similar with the strike of the southeastern Pishan fault. However, this fault segment responsible for the Pishan event has not been mapped before. Our finite fault model reveals that the peak slip of 0.89 m occurred at a depth of 11.6 km, with substantial slip at a depth of 9–14 km and a near-uniform slip of 0.2 m at a depth of 0–7 km. The estimated moment magnitude was approximately Mw 6.5, consistent with seismological results.

1. Introduction

On 3 July 2015 (09:07:46 GMT+8), a Mw 6.5 earthquake struck Pishan in Xinjiang, western China (Figure 1), causing at least four deaths, 48 injuries and hundreds of building collapses. The epicenter was located in the boundary between the southwestern Tarim Basin and the northwestern Tibetan Plateau. The collision between them, due to the northward push from the Indian plate and the subduction of the Tarim Basin, led to the development of the Pamir Plateau, the Karakoram and the west Kunlun orogens, which are associated with many well-known active tectonic features such as the Altyn Tagh fault system [1], where the 2008 Mw 7.1 Yutian earthquake occurred (Figure 1). However, there are no active Holocene faults within the approximately 50 km region surrounding the 2015 Pishan earthquake, according to the map of active tectonics in China [2]. Earthquakes with well-determined source processes would help us better understand the tectonic behaviors of the Tarim and surrounding regions. However, none of the historical earthquakes in the region of Figure 1 have received adequate geodetic coverage. Fortunately, complete Synthetic Aperture Radar Interferometry (InSAR) coverage of the 2015 Pishan event provides a unique opportunity to determine the location, geometry, and focal mechanism of the seismogenic fault, which may play an important role in elucidating the postcollisional convergence of India and Eurasia.

The focal mechanism solutions from the United States Geological Survey (USGS), Global Centroid Moment Tensor (GCMT) and China Earthquake Networks Center (CENC) Centroid Moment Tensor (CMT) indicate that the 2015 Pishan earthquake was primarily thrust faulting with a minor left-lateral strike-slip component. However, there are remarkable differences in focal mechanisms, locations and seismic moments among these results because of the different seismic inversion methods and data sources used [3]. Near-field geodetic measurements with good coverage can provide independent constraints on the source parameters with high accuracy, especially when no surface break is observed to constrain the fault geometry [4,5,6,7]. Such high-precision source parameters can enhance our understanding of the tectonics in the southwestern Tarim Basin.

Figure 1. Tectonic setting of the 2015 Pishan earthquake, western China. The focal mechanism of the 2015 Pishan earthquake is from the GCMT catalog. Red circles are the aftershocks with M ≥ 2.0 from the International Seismological Centre [8] between 3 July 2015 and 22 November 2015. Blue boxes outline the spatial coverage of the Sentinel-1A SAR images (ascending track T056A and descendisng track T136D) and the ALOS-2 Phased Array type L-band Synthetic Aperture Radar-2 (PALSAR-2) images (ascending track P160A). Black lines are the published faults from the map of active tectonics in China [2]. Blue vectors are the interseismic GPS velocities relative to a Eurasian reference frame with 95% confidence [9].

Figure 1. Tectonic setting of the 2015 Pishan earthquake, western China. The focal mechanism of the 2015 Pishan earthquake is from the GCMT catalog. Red circles are the aftershocks with M ≥ 2.0 from the International Seismological Centre [8] between 3 July 2015 and 22 November 2015. Blue boxes outline the spatial coverage of the Sentinel-1A SAR images (ascending track T056A and descendisng track T136D) and the ALOS-2 Phased Array type L-band Synthetic Aperture Radar-2 (PALSAR-2) images (ascending track P160A). Black lines are the published faults from the map of active tectonics in China [2]. Blue vectors are the interseismic GPS velocities relative to a Eurasian reference frame with 95% confidence [9].

According to the GCMT catalog, only two moderate earthquakes (5 < M < 6) were recorded in the earthquake region in the past 40 years, implying that the seismicity is low. There were more (~1012) aftershocks (M ≥ 2) [8] in the region before 22 November 2015 (Figure 1). Among these events, the largest aftershock is an ML 5.0 event, which followed approximately two hours later. These aftershocks are spatially confined to a ~50 by ~40 km nearly NW-SE-trending block (Figure 1). The block almost corresponds to the surface projection of the ruptured thrust fault of the 2015 Pishan earthquake, the distribution pattern of which was similar to that of the 2008 Mw 6.8 Zemmouri earthquake in Algeria [10] and the 2013 Mw 6.6 Lushan earthquake in China [11].

In this study, the coseismic deformation caused by the 2015 Pishan earthquake is measured with the InSAR technique using both ascending and descending Sentinel-1A and ascending ALOS-2 SAR data. The quasi-eastward and quasi-upward surface displacements are determined based on the ascending and descending observations. The fault geometry and slip distribution are retrieved from a combination of nonlinear and linear inversions together with the Helmert variance component estimation (HVCE) weighting method [12]. The stress state in the surrounding area is also investigated through Coulomb failure stress analysis. Finally, we discuss the relationship among the rupture geometry, the fault slip patterns and the potential seismic hazards in this area.

2. Geological Background

The 3 July 2015 Pishan earthquake occurred in the southwest boundary of the Tarim basin, which is the largest inland basin in China. The Tarim basin is bounded by the Tianshan Mountains to the north, the western Kunlun Mountains and the Altyn Tagh Mountains to the south and southeast, and the Pamir Plateau to the west (Figure 1). The basin is covered by a thick foreland sequence of Mesozoic to Cenozoic continental sedimentary rocks with a maximum thickness of 17 km [13]. In the region of this event, there is a famous fault system, the Altyn Tagh fault system, which is located in the Altyn Tagh Mountain range in northern Tibet, extends for ~1200 km along the strike, which consists of the ENE-trending Altyn Tagh fault in the middle, the WNW-trending western Kunlun fault, and the Nanshan thrust belts at the two ends [14,15].

The western part of the Altyn Tagh fault surrounding the Pishan earthquake can be divided into three main segments. (1) One is the western Kunlun thrust system, which has undergone ~100 km of north-south shortening since the early Oligocene [15]. With sparse GPS observations, Li [9] presented an upper bound on present-day slip rates of less than 2 mm/year. Approximately 5 mm/year convergence was inferred for structures within the western Kunlun fault [16]; (2) Another segment is the left-slip Karakax fault parallel to the western Kunlun thrust system, which is a re-activating Jurassic–Cretaceous left-lateral shear zone along the northern boundary of the western Kunlun batholith [17]. Using GPS data between 1993 and 1998, Shen et al. [18] estimated a slip rate of ~9 mm/year, which is very different from the geological result of ~20 mm/year from Landsat images [19]. Similarly, Wright et al. [20] suggested that the Karakax fault had a 5 ± 5 mm/year left-lateral slip rate based on five 500-km-long interferograms from European Remote Sensing-1 (ERS-1) and -2 satellites. Recent geological studies also report lower slip rates that are comparable to the geodetic-derived slip rate [21,22]. The Karakax fault is very seismically active, with a series of Mw ~7 historic earthquakes. The most recent major earthquake is the 2014 Yutian Mw 6.9 event that occurred on its eastern end [23]; (3) The third segment consists of several left-lateral faults and small pull-apart basins distributed along the southwest of the main trace of the Altyn Tagh fault, which has suffered two major normal-faulting ruptures recently (the 2008 Mw 7.1 Yutian event and the 2012 Mw 6.2 Yutian event).

3. InSAR Observations

3.1. SAR Data Processing

During the 2015 Pishan earthquake, several InSAR satellites collected measurements of surface displacement from space. The Sentinel-1A satellite from the European Space Agency collects the C-band data, which are provided as free and open-access data through the Scientific Data Hub. The Japan Aerospace Exploration Agency (JAXA) ALOS-2 satellite collects L-band data, which are delivered from the JAXA ALOS-2/ALOS User Interface Gateway (AUIG2). Because data with different imaging geometries can provide more information about surface deformation [24,25], a more reliable coseismic slip model can be retrieved from interferograms from both ascending and descending orbits of different satellites.

As part of the Copernicus mission, Sentinel-1 is composed of a constellation of two satellites (A and B) carrying a C-band SAR instrument to assure data continuity provided by the ERS-1/2 and Envisat satellites. Sentinel-1A, launched on 3 April 2014, operates in four image modes with different resolution and coverage: Strip Map (SM), Interferometric Wide Swath (IWS), Extra Wide Swath (EWS), and Wave (WV) modes. During these modes, both the IWS and EWS modes employ the novel Terrain Observation with Progressive Scans (TOPS) acquisition mode [26]. With the TOPSAR technique, in addition to steering the beam in range as in ScanSAR, the beam is electronically steered from backward to forward in the azimuth direction for each burst, avoiding scalloping and resulting in homogeneous image quality throughout the swath [27]. In this study, both ascending track T056A and descending track T136D Sentinel-1A IWS mode TOPS data covering the 2015 Pishan event rupture were acquired between 24 June 2015 and 24 July 2015 (Table 1). The original TOPS data are composed of three sub-swaths, with each sub-swath consisting of a series of bursts, which can provide a large coverage of 250 km with resolutions of 5 m and 20 m along the range and azimuth directions, respectively. We only used sub-swath 1 (bursts 1–4) and sub-swath 2–3 (bursts 8–11) from tracks T056A and T136D, respectively, to construct the TOPS interferograms, which is enough to capture main deformation signal of the event.

Table 1. Details of satellite synthetic aperture radar (SAR) images used in this study.

Satellite	Track	Master	Slave	Perp. B	Inc. Angle	Azi. Angle	σ †	α ǂ
		YYYYMMDD	YYYYMMDD	m	°	°	mm	km
Sentinel-1A	T056A	20150630	20150724	−53	33.7	−169.3	6.0	15.0
	T136D	20150624	20150718	26	41.4	−9.8	2.5	9.9
ALOS-2	P160A	20150222	20150726	−245	36.4	−169.7	6.9	7.9

^† Standard deviation calculated with all points in the non-deforming area; ^ǂ The e-folding correlation length scale of 1-D covariance function.

Table 1. Details of satellite synthetic aperture radar (SAR) images used in this study.

Satellite	Track	Master	Slave	Perp. B	Inc. Angle	Azi. Angle	σ †	α ǂ
		YYYYMMDD	YYYYMMDD	m	°	°	mm	km
Sentinel-1A	T056A	20150630	20150724	−53	33.7	−169.3	6.0	15.0
	T136D	20150624	20150718	26	41.4	−9.8	2.5	9.9
ALOS-2	P160A	20150222	20150726	−245	36.4	−169.7	6.9	7.9

^† Standard deviation calculated with all points in the non-deforming area; ^ǂ The e-folding correlation length scale of 1-D covariance function.

ALOS-2, the successor of the ALOS satellite, was launched on 24 May 2014 carrying an L-band PALSAR-2 instrument to continue the L-band SAR observation. The ALOS-2 satellite also supports several image modes, including the Spotlight mode with a width of 25 km, the Fine mode (Strip Map) with a swath width of 70 km, and ScanSAR (Wide Swath), with a width of 350 km. Compared with the C-band SAR, the PALSAR-2 preserves high coherence in interferometry over regions of rugged topography with high rainfall and dense vegetation because of its longer wavelength [28]. In this study, ascending track P160A ALOS-2 Fine mode data spanning from 22 February to 26 July, 2015 (Table 1) were collected for standard interferograms.

The Sentinel-1A SAR data and ALOS-2 PALSAR-2 data were processed using the Switzerland GAMMA software [29]. All interferograms were generated from the Single Look Complex (SLC) products. The multi-look ratio between the range and azimuth direction was 10:2 for Sentinel-1A and 8:8 for ALOS-2 data. ALOS-2 Fine mode interferometry is the same as the traditional interferometry. For TOPS interferometry, however, an accuracy of a few thousandths of a pixel co-registration in the azimuth direction is required [30]. Otherwise, phase jumps between subsequent bursts would be observed. To ensure very high co-registration accuracy, a method considering the effects of the scene topography and a spectral diversity method [31] considering the interferometric phase of the burst overlap region were used. After the high-quality co-registration between the TOPS SLC data, the effects of topography were removed from the interferograms using agency-provided orbits together with the 90-m-resolution Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) [32]; slight height-related signals still remained in the interferograms. The interferograms were then filtered using a power spectrum filter [33] to reduce the effects of phase noise and unwrapped using the branch cut method [34]. The interferograms were finally geocoded to the WGS84 geographic coordinates with 90-m resolution.

3.2. InSAR Coseismic Deformation

For the geocoded interferograms, especially the ALOS-2 interferogram, strong long-wavelength signals or ramps remained, which may be attributed to the inaccurate orbits. Fortunately, there is no significant ionospheric disturbance found in the ALOS-2 interferogram. In addition, the Pishan earthquake occurred in the southwestern corner of the Tarim Basin, which has a dry climate, but there are still atmospheric phase delays in the interferograms, especially in the southwest part. The correlation analysis between phase and topography indicates that topography-correlated atmospheric effect correction is recommended. To remove both the residual orbital errors and topography-correlated atmospheric delays, a linear function among location (x, y), elevation (h) and error phase [35] were estimated with observations away from the deformed areas. Figure 2 shows the coseismic interferograms with topography-correlated atmospheric delay and long-wavelength orbital errors corrected.

Figure 2. Coseismic interferograms obtained from the ascending Sentinel-1A track T056A (a), the descending Sentinel-1A track T136D (b) and the ascending ALOS-2 track P160A (c). The interferograms are rewrapped with an interval of 6 cm.

Figure 2. Coseismic interferograms obtained from the ascending Sentinel-1A track T056A (a), the descending Sentinel-1A track T136D (b) and the ascending ALOS-2 track P160A (c). The interferograms are rewrapped with an interval of 6 cm.

The Sentinel-1A and ALOS-2 interferometry provided excellent coverage and coherence for the Pishan earthquake. The surface deformation is clearly visible from both ascending and descending interferograms because of a high-quality spatio-temporal correlation between images acquired for a relatively flat and dry desert with little vegetation (Figure 2). The bulls-eye pattern of negative range change in the three interferograms represents motion toward the satellite over an area of approximately 42 by 28 km. The relative maximum LOS displacements are 14.6, 14.0 and 14.0 cm for the T056A, T136D and P160A interferograms, respectively. The slight discrepancy between the two ascending interferograms can be attributed to the different look angles of Sentinel-1A (~34°) and ALOS-2 (~37°).

Although the main errors including orbital errors and topography-correlated atmospheric delay were removed from the interferograms, there are still some un-modeled errors remained in the corrected interferograms. Because the residual errors of interferograms are highly spatially correlated, a 1D covariance function [36] in the non-deforming areas was used to describe the uncertainty characteristics of the interferograms (Figure 2). The standard deviations were 6.0 mm, 2.5 mm and 6.9 mm for the T056A, T136D and P160A interferograms, respectively, indicating that the interferograms had a low noise level.

3.3. 2.5D Surface Deformation

Theoretically, if LOS displacements with two different incidence angles on both ascending and descending tracks are available, it is possible to recover a three-dimensional surface deformation field from these measurements using a least-squares method. However, the error of the north component of the displacement is approximately 10 cm when applying this method, which is much larger than that of the east and vertical components (~1 cm) [20]. The difficulty in the north component estimation is mostly attributed to the poor view of the north component afforded by the near-polar orbits of the satellites. Instead, a technique of 2.5-dimensional (2.5-D) analysis [37] is adopted to obtain quasi-eastward and quasi-upward components of the 2015 Pishan earthquake from ascending and descending LOS data. Because of the near-polar orbital geometry (about 10° deviation), the LOS contribution of north component (sine of azimuth angle) is small. Then, the LOS displacements U associated with orbits can be written as

$$U\approx U_e{s}_e+U_u{s}_u,$$

(1)

where U_e and U_u are quasi-eastward and quasi-upward displacements, respectively. In addition, s_e and s_u are the direction cosines in east and vertical direction, respectively.

Figure 3 shows the 2.5-D surface deformation map (2-D spatial distribution of the 2-D displacement) of the Pishan earthquake, including a quasi-eastward displacement map and a quasi-upward displacement map. The greatest deformation occurred in southwestern Pishan, approximately 30 km away. The ground uplift region (red region in Figure 3b) with a maximum value of 13.6 cm located in the northeast compared with the major ground subsidence region (blue region in Figure 3b) with a maximum value of 3.2 cm. The horizontal displacements were far less than the vertical displacements caused by the Pishan earthquake. The reported focal mechanisms from USGS, GCMT and CENC (Table 2) indicate two candidate seismogenic faults, southwest dipping and northeast dipping. However, our InSAR measured surface displacements only stand for a southwest dipping fault [38]. The potential location of the seismogenic fault extending to the ground is very close to the blind southeastern Pishan fault, which is an Early-Middle Pleistocene fault in the active tectonics map of China [2] about which little is known. In addition, the 5.6 cm eastward ground displacements suggest the presence of dextral strike-slip components during this event.

Figure 3. Surface displacement of quasi-eastward (a) and quasi-upward (b) components derived from a combination of the Sentinel-1A ascending and descending tracks and the ALOS-2 ascending track.

Figure 3. Surface displacement of quasi-eastward (a) and quasi-upward (b) components derived from a combination of the Sentinel-1A ascending and descending tracks and the ALOS-2 ascending track.

Table 2. Source parameters of the 2015 Pishan earthquake.

Source	Lon.	Lat.	Strike	Dip	Rake	Depth	Length	Width	Slip	Mw
	°	°	°	°	°	km	km	km	m
USGS-BW	78.154	37.459	98	34	72	23.0	-	-	-	6.2
GCMT	78.14	37.58	109	22	85	15.6	-	-	-	6.4
CENC	78.1	37.5	115	23	72	22	-	-	-	6.5
Uniform slip model	78.057	37.571	114.0	23.6	92.6	8.8	22.1	10.1	0.59	6.4
	±0.4 km	±0.3 km	±1.6	±1.5	±3.2	±0.4	±0.5	±1.0	±0.06
Distributed slip model	78.057	37.571	114.0	25.0	97.0	-	36	40	0.34	6.5

Table 2. Source parameters of the 2015 Pishan earthquake.

Source	Lon.	Lat.	Strike	Dip	Rake	Depth	Length	Width	Slip	Mw
	°	°	°	°	°	km	km	km	m
USGS-BW	78.154	37.459	98	34	72	23.0	-	-	-	6.2
GCMT	78.14	37.58	109	22	85	15.6	-	-	-	6.4
CENC	78.1	37.5	115	23	72	22	-	-	-	6.5
Uniform slip model	78.057	37.571	114.0	23.6	92.6	8.8	22.1	10.1	0.59	6.4
	±0.4 km	±0.3 km	±1.6	±1.5	±3.2	±0.4	±0.5	±1.0	±0.06
Distributed slip model	78.057	37.571	114.0	25.0	97.0	-	36	40	0.34	6.5

4. Source Parameter Modeling

4.1. Data Reduction and Weighting

The InSAR observations shown in Figure 2 consist of millions of pixels, and it is impractical to use all of them to invert the model. To reduce the number of points and to improve the computational efficiency, several methods can be used to downsample the original interferograms, including uniform sampling [39], quadtree sampling [40], resolution-based quadtree sampling [41] and equation-based sampling [42]. We adopted the equation-based sampling approach, which can distinctively sample near-field points according to a deformation gradient and prevent the far-field noises based on the difference in Green’s function [42]. The location-dependent look vectors of each observation were then calculated considering the differences of incidence and azimuth angle based on the precise orbit data and local topography. Finally, 502 and 568 observations were obtained from the Sentinel-1A tracks T056A and T136D interferograms, respectively, and 416 observations were obtained from the ALOS-2 track P160A interferogram.

Because three different datasets were used for modeling, an appropriate determination of the weighting of each dataset was required during inversions. Although the empirical errors of the three interferograms derived from the 1-D covariance function (Section 3.2) can be used to weight the three datasets directly, Xu et al. [12] suggested that the HVCE method has various advantages over other methods, especially when the number of redundant observations is relatively large in a joint inversion. In this study, the HVCE method was employed to weight the three different observations to derive the coseismic slip distribution of the Pishan earthquake. In the HVCE method, the relative weight ratio is obtained according to the posterior misfit information of each dataset. During inversion, the starting value for the relative weight ratio among the datasets was given as 0.131:1:0.174 based on the variances (standard deviations) of the three interferograms. After 5 iterations, the variances of unit weight became almost uniform for all three datasets, and the relative weight ratio converged to 0.344:1.000:0.495.

4.2. Finite Fault Slip Model

The source parameters for the 2015 Pishan event were modeled through joint inversions of the ascending and descending Sentinel-1A and the ascending ALOS-2 data using the analytical solutions of a rectangular dislocation in a homogeneous, elastic half-space [43]. During inversion, a two-step strategy was employed to solve the fault parameters of the one-segment model. In step one, a nonlinear inversion was used to determine a set of model parameters (fault geometry) by minimizing the misfit under the assumption of a uniform slip on a rectangular fault; in step two, a linear inversion was used to estimate the slip distribution on the modeled fault plane.

In this study, a hybrid minimization algorithm basing on multi-peak particle swarm optimization (M-PSO) [44] was used to invert nine fault parameters, including location, strike, dip, length, depth, width and slip, by minimizing the misfits between the observed and the model predictions, assuming a Poisson ratio of 0.25 and a shear modulus of 3.32 × 10¹⁰ N/m² [45]. A Monte Carlo bootstrap simulation technique [46] was used to estimate the uncertainties of the fault parameters. Model solutions from 100 simulations perturbed with noise from the statistical properties based on previous 1-D covariance functions were used to estimate the standard deviation from their distributions. Table 2 shows all fault parameters and their errors. In general, the errors are relatively small. The inversion results indicate that the rupture fault is a blind thrust fault dipping approximately 24° to the southwest. The inverted geodetic moment magnitude is Mw 6.4, which is consistent with the GCMT result and slightly lower than that of CENC, but larger than the value from USGS.

Under the determined fault geometry, the slip on the rupture plane shows a linear relationship with surface displacements based on classic linear-elastic dislocation theory. Fixing the fault geometry for the optimal fault plane determined in the uniform slip modeling enables the fault length and width to be extended to 36 km along the strike and 40 km along the down-dip direction, respectively, and then discretized into 1440 patches with a size of 1 km by 1 km. Then, a constrained least squares method is used to solve the objection function:

$$F\left(s\right)=Gs-d^2+{\kappa }^{2}L{s}^2,$$

(2)

where G is Green’s matrix relating unit slip at the patches to the predicted displacement, $s$ is the slip vector on each patch, d is the observed InSAR LOS displacement, L is the second-order finite difference approximation of the Laplacian operator used to avoid unphysical oscillating slip distribution, and k² is the smoothing factor.

Previous studies (e.g., [44,47,48]) indicate that the fault geometry (especially the dip angle) derived from the uniform slip model is not optimal for the distribution slip model. In this study, the grid search method proposed by Feng et al. [44] is used to determine the optimal dip and smoothing factor simultaneously. In their method, a log function $f\left(\delta ,{\kappa }^2\right)=\log \left(\psi +\xi \right)$ is used, where $\delta$ is the dip angle, ψ is the model roughness and ξ is the residuals. The optimal dip angle and smoothing factor derived from the log function are 25° and 2.5°, respectively (Figure 4). The dip angle is slightly different from the nonlinear inversion results.

Figure 4. Plot of the log-function of the sum of root-mean-square error and model roughness for a range of models with different smoothing factors (k²) and dip angles. Red cross is the global minimum, which is chosen for the inversion presented in Figure 5.

Figure 4. Plot of the log-function of the sum of root-mean-square error and model roughness for a range of models with different smoothing factors (k²) and dip angles. Red cross is the global minimum, which is chosen for the inversion presented in Figure 5.

The preferred coseismic slip distribution and uncertainties from joint inversion of the ascending and descending Sentinel-1A and the ascending ALOS-2 observations are shown in Figure 5. The rupture fault is dominated by thrust motion with a slight strike-slip component. The total released geodetic moment is approximately 6.1 × 10¹⁸ N m (Mw 6.5), which is in agreement with the CENC result but larger than those of GCMT and USGS. The predominant slip occurred on the fault with a peak magnitude of 0.89 m at a depth of 9–14 km (width of 21–33 km). In addition to the main slip, the upper layer 0–7 km is found to have a near-uniform slip with a value of ~0.2 m and with 20% of the total moment, but there is no finding of surface rupture from field investigation organized by Institute of Geology, China Earthquake Administration. The shallower asperity corresponds to an Mw 6.0 sub-event (compared with the maximum aftershock of ML 5.0) with rupture close to the surface, although there is a relatively large uncertainty in the same area. The average errors of slip are ~5 cm, except for the area close to the surface, where the error increases to 7 cm, which can be attributed to reduced measurement constraints.

Figure 5. The 1 km × 1 km finite fault model of the 2015 Pishan earthquake (a), standard deviations in slip from the Monte Carlo estimation with 100 perturbed datasets (b) and the sum of moment release along strike distance (c).

Figure 5. The 1 km × 1 km finite fault model of the 2015 Pishan earthquake (a), standard deviations in slip from the Monte Carlo estimation with 100 perturbed datasets (b) and the sum of moment release along strike distance (c).

Figure 6 shows the simulated interferograms and residuals from our best fitting slip model. It is clear that the general patterns of both Sentinel-1A and ALOS-2 observations can be sufficiently explained by the distributed slip model. There are no notable residual fringes in the surrounding region of the seismogenic fault. The standard deviations between the InSAR observations and the simulated LOS displacements are 5.1 mm, 2.4 mm and 4.2 mm for tracks P056A, T136D and T160A, respectively, which are close to the noise levels of these InSAR observations. The correlation coefficient between the observations and predictions is 99.4%.

Figure 6. Modeled interferograms for the ascending Sentinel-1A track T056A (a), the descending Sentinel-1A track T136D (b) and the ascending ALOS-2 track P160A (c); and their residuals (d–f) with the distributed slip model. The white line is the top boundary of the distributed slip model, and the black dashed lines represent the surface projections of the modeled fault.

Figure 6. Modeled interferograms for the ascending Sentinel-1A track T056A (a), the descending Sentinel-1A track T136D (b) and the ascending ALOS-2 track P160A (c); and their residuals (d–f) with the distributed slip model. The white line is the top boundary of the distributed slip model, and the black dashed lines represent the surface projections of the modeled fault.

5. Stress Readjustment

A coseismic fault slip induces stress change to the surrounding area that may affect the state of neighboring faults and seismicity rates. To estimate the state of stress, we adopted the static Coulomb failure stress change [49]. From the Coulomb failure law, the static Coulomb failure stress $\Delta CFS$ on a given receiver fault plane following an event is given by

$$\Delta CFS=\Delta \tau +{\mu }^{\prime }\Delta \sigma ,$$

(3)

where Δτ is the fault-parallel shear stress change on the receiver fault plane (positive when sheared in the direction of fault slip), Δσ is the fault-normal stress change (positive when unclamping) and μ′ is the effective coefficient of friction. Increasing Coulomb stress, due to the stress readjustment following the event, tends to promote the failure of that plane. Decreasing Coulomb stress, representing stress release, inhibits the failure. Although the Coulomb stress following an event is relatively smaller than the tectonic stress required for an earthquake, events can be sufficiently triggered by a low Coulomb stress of 0.1 bar [49,50].

In this study, the Matlab-based software Coulomb v3.4 [51] was used to calculate the coseismic static Coulomb stress change triggered by the Pishan earthquake. The effective coefficient of friction and the shear modulus were set to be 0.4 and 3.32 × 10¹⁰ N/m² [45], respectively. The receiver faults in the studied area have a strike of 114°, a dip of 25°, and a rake of 97°, which are the same as the source fault parameters of the 2015 Pishan event.

The calculated Coulomb stress at a depth of 8 km (the top boundary of the main slip asperity) is shown in Figure 7a. Additionally, to illustrate the distribution of Coulomb stress with depth, a cross-section perpendicular to the rupture fault of the event is shown in Figure 7b. It is clear from Figure 7 that the pattern of Coulomb stress is controlled by both the initiation and termination ends of the main rupture fault. The modeled stress readjustment shows that the Pishan earthquake increased the Coulomb stress by approximately 0.5 bar (0.05 MPa) on both ends of the surface trace of the fault and slightly decreased the Coulomb stress on the West-Southwest extended region of the rupture. The results indicate that this event has increased the static Coulomb stress loading on the shallower (S1) and deeper (S2) sections (Figure 7b), which may bring them close to future failure.

Figure 7. Calculated static Coulomb failure stress changes caused by the 2015 Pishan event with the preferred InSAR distributed slip model and a friction coefficient of 0.4. (a) The distribution of Coulomb stress change at a depth of 8 km; (b) Cross-section of Coulomb stress change through profile AB.

Figure 7. Calculated static Coulomb failure stress changes caused by the 2015 Pishan event with the preferred InSAR distributed slip model and a friction coefficient of 0.4. (a) The distribution of Coulomb stress change at a depth of 8 km; (b) Cross-section of Coulomb stress change through profile AB.

6. Discussion

The fault geometry and slip characteristics of the 2015 Pishan earthquake are well constrained using both the ascending and descending Sentinel-1A and the ascending ALOS-2 SAR data (Table 2). This is the first time that a moderate earthquake in this region, the northwestern edge of the Tibetan Plateau, has been investigated with high-quality geodetic observations. We superimposed the earthquake locations from USGS, GCMT and CENC onto the 2.5-D surface displacement map (Figure 3) in order to discuss the hypocenter differences between the seismic CMT solution and the InSAR inversion. The locations from the different agencies appear to have been shifted eastward approximately 10 km to the peak vertical displacement region, which is far larger than the SAR geocoding errors. The comparison between the results from InSAR and seismology (Figure 7) shows a similar pattern. Because the InSAR data provide complete spatial coverage in this study, the location of the rupture fault has small errors (0.4 km). Thus, the hypocenter differences from this study indicate that the location from CENC may be more accurate than locations from USGS and GCMT. These differences are most likely caused by seismic network station distributions in western Tibet and the different methods used by the agencies. The location from the InSAR inversion indicates that the fault segment that ruptured in the 3 July 2015 event is likely to be the western segment of a blind southeastern Pishan thrust fault (Figure 6).

The centroid depth in the uniform slip model and the depth of maximum slip in the distributed slip model are 10.8 km and 11.6 km, respectively, which are shallower than the depths from the USGS and CENC and slightly shallower than the depth from GCMT (Table 2). The predominant slip occurred at a depth of 9–14 km. That the InSAR centroid depth is shallower than the seismic depths is consistent with the median discrepancy of 5 km from previous studies [52]. The first plausible reason is variations in the properties of the local upper crust, as our results are derived from a homogeneous, elastic half-space model, which is the first-order approximation to the true crust. However, in some instances, the depth from the layer crust model is shallower than that from a half-space model [3]. The other possible cause is that the seismic depth is the rupture initiation depth, whereas the depth in the InSAR model is the centroid depth of the rupture fault or the depth of the peak slip. Mai et al. [53] found that the shallow thrust rupture tends to propagate upward from the deeper sections of the fault toward the surface from more than 80 finite-source rupture models of over 50 earthquakes. Another possible explanation for the shallow InSAR centroid depth is the distribution of the measurements. In contrast to the almost far-field data used in the seismic inversion, the InSAR data used in our source modeling is almost near-field (0–100 km) observation, which can provide a strong constraint to the shallow slip on the rupture fault.

The preferred slip model shows that the rupture of the 2015 Pishan earthquake is dominated by thrust faulting (4%–5% strike-slip component). The peak slip has a 0.89 m thrust and 0.01 m right-lateral slip, respectively, with a rake angle of 89°. These results are different from the seismic results from the USGS, GCMT and CENC CMT solutions based on broadband seismograms. The geodetic and seismic moments generally agree, indicating that the largest part of the rupture radiated seismically. Additionally, the differences between these findings may suggest that the geodetic estimates are probably affected by early (first ~20 days) post-seismic and aftershock deformation.

Using high-resolution seismic reflection data acquired in the foothills of the western Kunlun range, Jiang et al. [54] revealed that there are several thrusts, fault-propagation folds, and triangular zones rooted in detachments that transfer shortening from the mountain northward into the Tarim Basin. The deeper asperity of the 2015 earthquake indicates a partial rupture of the ramp fault at a 14–18 km depth that cuts through Paleozoic strata upward to the bottom of the Paleocene strata, transferring slip from the Proterozoic decollement beneath the range to the shallower Eocene layer under the Tarim Basin. Such main slip behavior is similar to that of the 2013 Mw 6.6 Lushan earthquake, which also occurred on a blind thrust fault in the boundary zone between the Longmen Mountain and the Sichuan Basin [11]. Both moderate events indicate an outward migration of the deformation front from the interior of the Tibetan Plateau, resulting in its lateral extrusion.

Approximately 0.2 m of slip was found on the fault patches at a 0–7 km depth from the InSAR modeling. The shallower asperity is distributed over a very broad area and releases approximately 20% (a Mw 6.0 event) of the total moment; however, it is accompanied by larger uncertainties of approximately 7 cm (30% of slip). This shallow asperity might be an artifact due to the use of elastic solution during the inversion. The surface deformation in the near field of the rupture fault is partly caused by an essentially inelastic response of the upper brittle crust, which was observed in the 2003 Bam earthquake [55] and the 2011 Burma earthquake [44]. Additionally, from the static Coulomb failure stress changes calculated from both the uniform slip model and the distributed slip model, the area of shallower asperity suffers increased stress changes, which indicates that the slip in the shallow layer may have arisen from postseismic afterslip or aftershocks triggered by the deeper asperity (Figure 7). If the accumulated strain energy in the shallower section is released by the afterslip or aftershocks, there will be no second, moderate earthquake shortly thereafter in the shallower section, such as the case in the 2008 and 2009 Qaidam Mw 6.3 earthquake sequences [56].

Using a balanced cross-section technique, Jiang et al. [54] estimated the rate of thrust slip along a detachment surface of 2.35 mm/year. Assuming that the 2015 event released all of the accumulated strain energy, an earthquake recurrence interval of 200–300 years for a moderate earthquake in this segment can be obtained from a similar slip rate and the average coseismic slip (0.59 m). Such a recurrence interval is comparable to 300–1100 years for major earthquakes along the western Altyn Tagh fault [57], which also indicates that the event released most of the previously accumulated strain energy in the segment.

7. Conclusions

In this study, a 2.5-D coseismic deformation map of the 2015 Pishan earthquake was derived from both the ascending and descending Sentinel-1A TOPS and the ascending ALOS-2 Fine Mode SAR satellite data. The substantial deformation concentrated in an area of approximately 42 by 28 km, and the maximum displacements in uplift and subsidence were approximately 13.6 cm and 3.2 cm, respectively. The vertical deformation pattern associated with non-linear inversion results shows that the southeastern Pishan fault is most likely responsible for the 2015 Pishan earthquake. The optimal dip and rake angles were 25° and 97°, respectively, indicating that this thrusting event had a slight dextral strike-slip component. The main slip occurred at a depth of 9–14 km, with a maximum slip of 0.89 m at a depth of 11.6 km, and a near-uniform slip of 0.2 m occurred at a depth of 0–7 km. The total released moment was approximately 6.1 × 10¹⁸, equivalent to Mw 6.5. The calculated earthquake recurrence interval of 200–300 years is comparable to that of major earthquakes along the western Altyn Tagh fault. Our investigation suggests that the rupture of this event has released most of the previously accumulated strain energy.