InSAR-Constrained Interseismic Deformation and Potential Seismogenic Asperities on the Altyn Tagh Fault at 91.5–95E, Northern Tibetan Plateau

The present-day kinematic features of the different segments of the Altyn Tagh Fault (ATF) have been well-studied using geodetic data. However, on the eastern segment of the ATF at 91.5–95◦E, high spatial resolution deformation has not been previously reported. Here, we processed 185 interferometric synthetic aperture radar (InSAR) images from three descending tracks of the C band ERS-1/2 and Envisat satellites spanning 1995–2011 and obtained the average deformation velocity field. Results show a left-lateral motion of ~4 mm/year along the fault-parallel direction across the ATF at 91.5–95◦E, which is consistent with Global Positioning System (GPS) observations. The slip deficit rate distribution at shallow depths was resolved through the InSAR deformation velocity using a discretized fault plane. The slip deficit is capable of an Mw 7.9 earthquake, considering the elapsed time of the latest M 7.0 event. Two potential asperities that could be nucleation sites or rupture areas of future earthquakes were delineated based on the coupling coefficient and seismicity distributions along the fault plane. The larger asperity extends more than 100 km along the ATF at depths of 8–12 km. Our InSAR observations support the undeformed blocks model of the Indo-Eurasian collisional mechanism at the northern margin of the Tibetan plateau.


Introduction
The approximately 1800 km long SW-NE sinistral strike-slip Altyn Tagh Fault (ATF), which is the boundary between the northern edge of the Tibetan Plateau and the Tarim block (Figure 1), has been the focus of an ongoing debate on the role of undeformed blocks (e.g., [1]) and continuum deformation (e.g., [2]) in the Indo-Eurasian collisional models. Based on geomorphological characteristics, the ATF is divided into three segments: the western segment (between 78-84 • E), the central segment (between 84-94 • E), and the eastern segment (between 94-97 • E) ( Figure 1) [3]. Previous studies have suggested that the ATF was formed during the Indosinian period. After the Paleo-Tethys Ocean closed during the Triassic period, the Altyn Tagh ductile shear zone was created by the NE oblique subduction and collision of the Qiangtang block with Gondwana [4,5]. During the early Cretaceous, the Middle Tethyan Ocean closed and the Altyn Tagh ductile shear zone was activated by the subduction of the Gangdise block from Gondwana to beneath the northern plateau. Due to the resistance of the Tarim-Alxa block to the north, the Qilianshan thrust fault and the Northern Tibetan Plateau prototype were formed during the 110-80 Ma period [6]. The closure of the New Tethys and collision of the Indian plate and Eurasian plate during the 60-40 Ma period reactivated the Altyn Tagh ductile shear zone. Consequently, the nascent ATF was formed at that time. During the Oligocene-Miocene period, Figure 1. Seismotectonic map showing the Altyn Tagh Fault (ATF) and its location in the India-Asia collision zone. Fault traces adopted from [8] are superimposed on Shuttle Radar Topography Mission (SRTM) digital elevation models (DEMs). Black triangles represent cities or towns. (a) The ATF is shown as a bold black line. Red circles represent M 7.0+ earthquakes during the period 1900-2017, with circle size proportional to magnitude. The two blue polygons delineate the regions threatened by M 7.0+ earthquakes as mapped by [9]. The black box covers the area shown in (b). (b) Red rectangles show the coverage of the analyzed ERS-1/2 and Envisat synthetic aperture radar (SAR) data, with track numbers indicated. Arrows are Global Positioning System (GPS) velocities with respect to Eurasia [10]. Black circles represent relocated earthquakes with M 3.0+ that occurred during 1980-2012, with circle size proportional to magnitude. and its location in the India-Asia collision zone. Fault traces adopted from [8] are superimposed on Shuttle Radar Topography Mission (SRTM) digital elevation models (DEMs). Black triangles represent cities or towns. (a) The ATF is shown as a bold black line. Red circles represent M 7.0+ earthquakes during the period 1900-2017, with circle size proportional to magnitude. The two blue polygons delineate the regions threatened by M 7.0+ earthquakes as mapped by [9]. The black box covers the area shown in (b). (b) Red rectangles show the coverage of the analyzed ERS-1/2 and Envisat synthetic aperture radar (SAR) data, with track numbers indicated. Arrows are Global Positioning System (GPS) velocities with respect to Eurasia [10]. Black circles represent relocated earthquakes with M 3.0+ that occurred during 1980-2012, with circle size proportional to magnitude. Due to logistical barriers and inclement weather, it is very difficult to carry out fieldwork on the ATF. Thus, paleo-earthquakes are not well documented. According to previous studies, there were at least five large paleo-earthquakes on the ATF from Late Pleistocene to Holocene, and the average recurrence period was approximately 800 years [9,11] The latest paleo-earthquake, near the town of Akesai, was confirmed to occur~665 years BP, with a recurrence of~600 years [12]. Since 1900, several destructive earthquakes with magnitude 7 or above are known from historical records, including the 1924 Minfeng M 7.3 earthquake on the central segment and the 1932 Changma M 7.6 earthquake on the easternmost end [13] (Figure 1). During the last 10 years, two Yutian M 7.3 events occurred, one in 2008 and another in 2012, on the southwestern segment of the ATF, with the two epicenters 100 km apart [14] (Figure 1). In general, the Akesai-Subei-Shibaocheng section of the ATF is believed to have remained unbroken for~700 years and is identified as a seismic gap [9].
The present-day slip rate of the ATF, which is an important kinematic description of the Himalayan collision process and a critical constraint for regional seismic hazard estimation, has been quantified by Global Positioning System (GPS) and Interferometric Synthetic Aperture Radar (InSAR) at different segments. At the western segment, i.e., 78-80 • E, the slip rate is as low as 2-5 mm/year [10,15,16]. Wang and Wright [16] estimated a left-lateral strike-slip rate of 1 ± 2 mm/year at 78 • E of the western ATF based on ERS-1/2 and Envisat SAR images. At 79.3 • E, Zheng et al. [10] obtained a sinistral rate of 1.4 ± 1.3 mm/year from a 200 km wide GPS profile. At 80 • E of the ATF, Wright et al. [15] and Wang and Wright [16] gave a slip rate of 5 ± 5 mm/year and 2 ± 2 mm/year, respectively. At 85-94 • E, the slip rate gradually decreases eastward from~11 mm/year [17][18][19] to~7 mm/year [10,[20][21][22][23][24][25][26][27][28]. Elliott et al. [18] measured the interseismic deformation of the ATF at 85 • E from ERS-1/2 images for the period 1993-2000, and they estimated the left-lateral strike-slip rate to be 11 ± 5 mm/year. Zhu et al. [19] used the same track Envisat data (2003-2010) as Elliott et al. [18], and their slip rate estimate was 8 ± 0.7 mm/year. Using Envisat SAR images, Daout et al. [17] derived a left-lateral strike-slip movement of 10.5 mm/year at 84.5-86.5 • E, below a 17 km locking depth. At 84.8-87.5 • E, Vernant [25] estimated a slip rate of 8.6 mm/year using GPS measurements, with a locking depth of 14.5 km. Shen et al. [23] estimated slip rate between 85-90 • E using GPS measurements, and the results show a left-lateral rate of 9 ± 2 mm/year, with a convergence rate of 0 ± 2 mm/year. At 86.0 • E, Zheng et al. [10] obtained a left-lateral rate of 8.1 ± 0.7 mm/year. At~86.2 • E of the ATF, He et al. [22] reported a slip rate of 9.0 mm/year based on a GPS array. Shen et al. [24] derived the interseismic deformation field over the ATF at 86.8 • E using Sentinel-1 interferograms spanning the period between late 2014 and 2016, with a slip rate and locking depth of 9 ± 2 mm/year and 14 ± 4 km, respectively. At 88.7-91.4 • E, Vernant [25] estimated a slip rate of 9 mm/year using GPS measurements, with locking depth of 25 km. For the segment from 89-91 • E of the ATF, Bendick et al. [20] reported GPS and leveling data that indicated a left-lateral shear rate of 9 ± 5 mm/year and a contraction rate of 3 ± 1 mm/year, while Wallace et al. [26] and Zhang et al. [28] showed a left-lateral slip rate of 9 ± 4 mm/year and 11.9 ± 3.3 mm/year, respectively.  [28]. Therefore, the left-lateral motion along the ATF cannot be described by a single slip rate. Although GPS and InSAR give slip rates at most of the segments of the ATF, GPS station density is not sufficiently high to determine the pattern of strain accumulation at a kilometer scale, while conventional InSAR measurements are often heavily hampered by atmospheric effects and signal decorrelation [29]. Moreover, our target segment (between 91.5 • E and 95 • E) has not been studied using InSAR, and no high spatial resolution measurements have been obtained in that portion of the fault until now. Since the strike of the said segment is nearly perpendicular to the ERS-1/2 and Envisat satellites' flying direction and the region has relatively low vegetation due to its arid to semiarid climate, the use of InSAR is ideal to determine interseismic strain accumulation. Here we present the results of our analysis of 185 SAR images from the C band ERS-1/2 and Envisat satellites spanning 1995-2011 using the InSAR technique.
Mapping the spatial extent of locked areas is key to understanding the segmentation of active faults, and identifying strongly coupled areas helps to constrain the timing, extent, and magnitude of future events (e.g., [30,31]). Here, the InSAR observations provided coverage at high spatial resolution of the ATF at 91.5-95 • E and allowed us to identify spatial slip deficit rates and coupling coefficient variations on the fault plane. Further, potential asperities that may trigger large events in the future were detected, and their implications for seismic hazard assessment are discussed.

InSAR Data Processing Strategy
Processed radar data were acquired by the ERS-1/2 and Envisat satellites on three contiguous descending tracks ( Figure 1, Table 1) operating at a wavelength of 56.3 mm. Between October 1995 and January 2011, 185 radar images covering the eastern ATF region were obtained, providing 16 years of continuous observations. Overlapping with each other, the three tracks together cover about a 250 km long continuous section of the ATF between longitudes 91.5 • E and 95.0 • E ( Figure 1). The basic InSAR processing techniques employed herein were similar to those used in our previous studies [32][33][34]. The interferometric processing was performed with GAMMA software [35]. The two-pass InSAR approach (e.g., [36,37]) was utilized to form interferograms. The effects of topography were removed from the interferograms using a filled 3 arc-sec (~90 m) resolution Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) [38] obtained from the Consultative Group on International Agricultural Research Consortium for Spatial Information (CGIAR-CSI, http://srtm.csi.cgiar.org). To improve the signal-to-noise ratio, interferograms were downsampled to 12 looks in range and 60 looks in azimuth (180 m × 180 m) and were filtered twice (with filters having window dimensions of 128 × 128 and 32 × 32 pixels) using an adaptive filter function based on the local fringe spectrum [39]. This filtering strategy efficiently removes high-frequency noise and makes the phase unwrapping considerably easier. Interferogram phase unwrapping was carried out using the minimum cost flow algorithm [40]. The reference point of phase unwrapping was chosen close to the fault trace, at the center of the interferogram. The final spatial resolution of the InSAR deformation field is 180 m in both azimuth and range directions.
The main limitation of low-amplitude tectonic signal detection using InSAR is that the orbital error and atmospheric delay projected into the line-of-sight direction can be at the same level as the tectonic signals, or even larger, reaching tens of centimeters [29]. To remove residual orbit errors, we applied the following procedure. First, an a priori velocity model by the interferometric time spans was established from a screw dislocation model [41] constrained by the GPS horizontal-component velocity [10]. Second, a residual interferogram was obtained by subtracting the a priori velocity model from the original interferogram. Third, a planar trend surface was estimated as orbital error and was removed from each residual interferogram. Fourth, the a priori velocity model was added back. Atmospheric errors consist of propagation delay through tropospheric and ionospheric layers. Here, we considered only tropospheric noise, because the ionospheric propagation delay is comparatively weak for the ERS-1/2 and Envisat C-band data [42]. The European Centre for Medium-Range Weather Forecasts (ERA-Interim) atmospheric reanalysis was used to correct the atmospheric delay [43]. Specifically, the global ERA-Interim weather model was applied to determine the atmospheric delay close to the time (i.e., within a few hours) of each SAR acquisition. Then, the phase-elevation relation at each grid node of the weather model was estimated and a delay map for each SAR acquisition date was obtained through horizontal and vertical interpolation. Finally, the atmospheric delay-corrected interferogram was achieved by removing the said weather model.
Interferograms spanning a short time interval tend to have higher coherence. However, they make small tectonic signals, such as those across the ATF, more difficult to detect. To minimize geometric decorrelation and errors due to the topography [44] and maximize the deformation signal in the data, we selected pairs of images by constraining the perpendicular and temporal baselines simultaneously (Table 2; Figure 2): (1) If the perpendicular baselines are smaller than 200 m, the temporal baselines can be longer than 1000 d; (2) If the perpendicular baselines are between 200 m and 400 m, the temporal baselines must be less than 1000 day and longer than 500 day; (3) If the perpendicular baselines are between 400 m and 600 m, the temporal baselines must be less than 500 d and longer than 350 day. Stacking the interferograms is one way to enhance the tectonic signal [29,[45][46][47]. By averaging many interferograms over the same area, random noise, such as atmospheric signals, can be subdued, small-scale persistent motion (such as interseismic deformation) can be highlighted, and cumulative trends can be extracted. Stacking N short-baseline and long-timespan interferograms can reduce the variance of atmospheric errors by a factor of N1/2 [29,46]. In this study, we used the Poly-Interferogram Rate And Time-series Estimator (π-RATE) software package developed by [18,46,47] to calculate the stacked velocity. Pixels that are coherent in 18 or more of the component interferograms were included in the stack. Redundant interferograms were used to recover the non-coherent pixels. Forecasts (ERA-Interim) atmospheric reanalysis was used to correct the atmospheric delay [43]. Specifically, the global ERA-Interim weather model was applied to determine the atmospheric delay close to the time (i.e., within a few hours) of each SAR acquisition. Then, the phase-elevation relation at each grid node of the weather model was estimated and a delay map for each SAR acquisition date was obtained through horizontal and vertical interpolation. Finally, the atmospheric delay-corrected interferogram was achieved by removing the said weather model. Interferograms spanning a short time interval tend to have higher coherence. However, they make small tectonic signals, such as those across the ATF, more difficult to detect. To minimize geometric decorrelation and errors due to the topography [44] and maximize the deformation signal in the data, we selected pairs of images by constraining the perpendicular and temporal baselines simultaneously (Table 2; Figure 2): (1) If the perpendicular baselines are smaller than 200 m, the temporal baselines can be longer than 1000 d; (2) If the perpendicular baselines are between 200 m and 400 m, the temporal baselines must be less than 1000 day and longer than 500 day; (3) If the perpendicular baselines are between 400 m and 600 m, the temporal baselines must be less than 500 d and longer than 350 day. Stacking the interferograms is one way to enhance the tectonic signal [29,[45][46][47]. By averaging many interferograms over the same area, random noise, such as atmospheric signals, can be subdued, small-scale persistent motion (such as interseismic deformation) can be highlighted, and cumulative trends can be extracted. Stacking N short-baseline and long-timespan interferograms can reduce the variance of atmospheric errors by a factor of N1/2 [29,46]. In this study, we used the Poly-Interferogram Rate And Time-series Estimator (π-RATE) software package developed by [18,46,47] to calculate the stacked velocity. Pixels that are coherent in 18 or more of the component interferograms were included in the stack. Redundant interferograms were used to recover the non-coherent pixels.

InSAR Velocity Field
Regional GPS measurements show that the ATF is dominated by horizontal left-lateral movement, and the vertical component is negligible relative to the horizontal component [23,28]. Therefore, InSAR-derived deformation in the line-of-sight (LOS) direction is mainly contributed by horizontal displacement. Figure 3 shows a map of the interseismic strain accumulation velocity along

InSAR Velocity Field
Regional GPS measurements show that the ATF is dominated by horizontal left-lateral movement, and the vertical component is negligible relative to the horizontal component [23,28]. Therefore, InSAR-derived deformation in the line-of-sight (LOS) direction is mainly contributed by horizontal displacement. Figure 3 shows a map of the interseismic strain accumulation velocity along the section of the ATF between 91.5 • E and 95 • E. As shown in the figure, the resulting mean LOS InSAR velocity field reveals a clear picture of the interseismic deformation along this segment of the ATF. It shows a gradual change in velocity across the fault and is consistent with the left-lateral sense of motion between the Tibetan plateau and the Tarim plate. Note that cool colors indicate movement away from (i.e., westerly), and warm colors indicate movement toward (i.e., easterly) the satellites, which are flying in descending orbits (i.e., heading southward). The interferograms we used are fairly evenly distributed in time, so the displayed stacks of interferograms approximated the average rate over the total time interval of the stacks. It appears that the velocity difference decreases from~2 mm/year in the west to~1 mm/year in the east along the radar's LOS direction.
Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 20 the section of the ATF between 91.5°E and 95°E. As shown in the figure, the resulting mean LOS InSAR velocity field reveals a clear picture of the interseismic deformation along this segment of the ATF. It shows a gradual change in velocity across the fault and is consistent with the left-lateral sense of motion between the Tibetan plateau and the Tarim plate. Note that cool colors indicate movement away from (i.e., westerly), and warm colors indicate movement toward (i.e., easterly) the satellites, which are flying in descending orbits (i.e., heading southward). The interferograms we used are fairly evenly distributed in time, so the displayed stacks of interferograms approximated the average rate over the total time interval of the stacks. It appears that the velocity difference decreases from ~2 mm/year in the west to ~1 mm/year in the east along the radar's LOS direction.  Figure 6. Arrows are GPS velocities with respect to Eurasia, as adopted from [10]. Stars denote the fault location as derived by screw dislocation inversions. Black lines represent the fault traces mapped by [8]. Dashed ellipses labeled A1, A2, A3, and A4 represent signals dominated by nontectonic components, which are zoomed in on the four right panels. The satellite flight direction and radar look direction are labeled as a solid arrow and an open arrow, respectively.

Validation of Overlapping Regions from Adjacent Tracks
Assuming that the LOS signal is due to surface motion and that this motion is predominantly horizontal and parallel to the local strike of the ATF (N73°E), we converted the mean LOS velocities to fault-parallel (FP) InSAR velocities, taking into account the viewing geometry, including the local radar look angle (e.g., [48]). Equation (1) explains the projecting calculation. In Equation (1), θ is the incidence angle for the pixel measured between the nadir and satellite looking vector for each pixel, α is the horizontal azimuth of the looking vector, and β is the strike of the fault. Dfault is the faultparallel displacement. DLOS is the InSAR LOS displacement: To estimate the uncertainties in the LOS data, we chose two profiles (see Figure 3 for locations) and calculated the differences and standard deviations in fault-parallel velocities in the overlapping  Figure 6. Arrows are GPS velocities with respect to Eurasia, as adopted from [10]. Stars denote the fault location as derived by screw dislocation inversions. Black lines represent the fault traces mapped by [8]. Dashed ellipses labeled A1, A2, A3, and A4 represent signals dominated by nontectonic components, which are zoomed in on the four right panels. The satellite flight direction and radar look direction are labeled as a solid arrow and an open arrow, respectively.

Validation of Overlapping Regions from Adjacent Tracks
Assuming that the LOS signal is due to surface motion and that this motion is predominantly horizontal and parallel to the local strike of the ATF (N73 • E), we converted the mean LOS velocities to fault-parallel (FP) InSAR velocities, taking into account the viewing geometry, including the local radar look angle (e.g., [48]). Equation (1) explains the projecting calculation. In Equation (1), θ is the incidence angle for the pixel measured between the nadir and satellite looking vector for each pixel, α is the horizontal azimuth of the looking vector, and β is the strike of the fault. D fault is the fault-parallel displacement. D LOS is the InSAR LOS displacement: To estimate the uncertainties in the LOS data, we chose two profiles (see Figure 3 for locations) and calculated the differences and standard deviations in fault-parallel velocities in the overlapping areas between neighboring tracks. The differences between neighboring tracks are approximately Gaussian, with mean values close to zero (Figure 4). The standard deviations between these independent estimates of fault-parallel velocities are 1.   Figure 5 shows profiles of fault parallel velocities from InSAR and GPS (from [10]) measurements. Visual comparison of the InSAR LOS velocities and the projected GPS velocities show good agreement between these two independent datasets. The difference between the two velocities has an average of 0.01 mm/year and a standard deviation of 1.05 mm/year. The profiles are consistent overall, with a classic arctangent shape predicted by elastic models across a strike-slip fault [41]. The topographic elevation along the profiles (shaded areas) is shown in order to reveal any potential correlation between the topography and the InSAR-derived velocities. As can be seen in the two profiles, correlation between the deformation signal and the local topography is absent, suggesting that tropospheric effects in the velocity field are insignificant. Note that at a distance of more than a few tens of kilometers from the fault zone, the profiles are satisfactorily flattened, suggesting that the residual orbital errors are negligible. The hypothesis of minimal vertical deformation is supported by the fact that the independent horizontal GPS measurements are in good agreement with the InSAR measurement; if a significant vertical velocity signal was present, disagreement should be found when comparing horizontal GPS displacement with InSAR data.  Figure 5 shows profiles of fault parallel velocities from InSAR and GPS (from [10]) measurements. Visual comparison of the InSAR LOS velocities and the projected GPS velocities show good agreement between these two independent datasets. The difference between the two velocities has an average of 0.01 mm/year and a standard deviation of 1.05 mm/year. The profiles are consistent overall, with a classic arctangent shape predicted by elastic models across a strike-slip fault [41]. The topographic elevation along the profiles (shaded areas) is shown in order to reveal any potential correlation between the topography and the InSAR-derived velocities. As can be seen in the two profiles, correlation between the deformation signal and the local topography is absent, suggesting that tropospheric effects in the velocity field are insignificant. Note that at a distance of more than a few tens of kilometers from the fault zone, the profiles are satisfactorily flattened, suggesting that the residual orbital errors are negligible. The hypothesis of minimal vertical deformation is supported by the fact that the independent horizontal GPS measurements are in good agreement with the InSAR measurement; if a significant vertical velocity signal was present, disagreement should be found when comparing horizontal GPS displacement with InSAR data. Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 20

Fault Geometry
We determined the fault geometry under the following considerations: First, although the target segment of the ATF is composed of two branches (Figure 3), it is difficult to distinguish which branch could be the most related to the velocity change. Moreover, the large-scale deformation signal appears to be controlled more by the fault system as a whole rather than by a single branch. Therefore, we assumed that the two branches merge as a single fault plane at depth. To simplify the fault geometry, we set only one fault plane instead of two branch fault planes. Second, the surface trace of the fault was determined by a screw dislocation model [41], constrained by five fault-normal InSAR profiles (see Figure 3 for their locations). More specifically, we modeled five profiles using screw dislocations on an infinitely long and vertical fault in an elastic half-space ( Figure 6). Three parameters were estimated in the model: slip rate at depth (mm/year), locking depth (km), and location of the fault (km) (i.e., horizontal shift). The modeled surface locations of the fault from the five profiles are shown as stars in Figure 3. The fault surface location is roughly located in the middle of the two branches. Then, we connected the five modeled locations of the fault as a polyline, which is assumed as the surface trace of the fault. In other words, the fault trace was approximated by six linear segments along the strike. Third, resistivity results from two magnetotelluric profiles that traverse the ATF at Akesai and Shibaocheng indicate that the ATF has a nearly vertical extension at depths less than 40 km [49]. Thus, in the modeling, the fault plane was set up to be vertical everywhere. Fourth, no microseismicity was detected at depths of more than 40 km, implying a transition to a ductile zone of the fault [50]. Therefore, we set up the down-dip width of the fault as 40 km.

Fault Geometry
We determined the fault geometry under the following considerations: First, although the target segment of the ATF is composed of two branches (Figure 3), it is difficult to distinguish which branch could be the most related to the velocity change. Moreover, the large-scale deformation signal appears to be controlled more by the fault system as a whole rather than by a single branch. Therefore, we assumed that the two branches merge as a single fault plane at depth. To simplify the fault geometry, we set only one fault plane instead of two branch fault planes. Second, the surface trace of the fault was determined by a screw dislocation model [41], constrained by five fault-normal InSAR profiles (see Figure 3 for their locations). More specifically, we modeled five profiles using screw dislocations on an infinitely long and vertical fault in an elastic half-space ( Figure 6). Three parameters were estimated in the model: slip rate at depth (mm/year), locking depth (km), and location of the fault (km) (i.e., horizontal shift). The modeled surface locations of the fault from the five profiles are shown as stars in Figure 3. The fault surface location is roughly located in the middle of the two branches. Then, we connected the five modeled locations of the fault as a polyline, which is assumed as the surface trace of the fault. In other words, the fault trace was approximated by six linear segments along the strike. Third, resistivity results from two magnetotelluric profiles that traverse the ATF at Akesai and Shibaocheng indicate that the ATF has a nearly vertical extension at depths less than 40 km [49]. Thus, in the modeling, the fault plane was set up to be vertical everywhere. Fourth, no microseismicity was detected at depths of more than 40 km, implying a transition to a ductile zone of the fault [50]. Therefore, we set up the down-dip width of the fault as 40 km.

Preprocessing InSAR Deformation Results
Before modeling, we processed the InSAR-derived deformation field considering the following aspects: First, to reduce the number of points, the InSAR result was downsampled through a regular grid with spacing of 3.6 km. Second, some nontectonic signals are included in the original InSAR deformation field. Considering the large-magnitude nontectonic signals, we only kept the points whose deformation velocities were within ±2 mm/year and discarded the other points. Third, according to the back-slip theory, the geodetic deformation includes two components, a global component accounting for block movement and a local component caused by fault movement. In this study, to remove the global component related to block movement, we split the InSAR data into two parts based on the fault trace that was determined in Section 3.1. and subtracted the corresponding mean value, leaving the residual component that is believed to be related to fault movement. As such, the inversed result is the slip deficit of the fault during the interseismic period. It is worth noting that the GPS surface deformation is not included in the modeling due to its very low spatial resolution ( Figure 1). Moreover, introducing the GPS data caused an additional level of complexity and required dealing with spatially heterogeneous data weight.

Modeling Approach
To detect variations in slip deficit along the fault plane, we subdivided the rectangular fault plane into 4 km (along the strike direction) × 4 km (down-dip direction) patches. The slip deficit distribution model was performed using the steepest descent method (SDM) code [51], which has

Preprocessing InSAR Deformation Results
Before modeling, we processed the InSAR-derived deformation field considering the following aspects: First, to reduce the number of points, the InSAR result was downsampled through a regular grid with spacing of 3.6 km. Second, some nontectonic signals are included in the original InSAR deformation field. Considering the large-magnitude nontectonic signals, we only kept the points whose deformation velocities were within ±2 mm/year and discarded the other points. Third, according to the back-slip theory, the geodetic deformation includes two components, a global component accounting for block movement and a local component caused by fault movement. In this study, to remove the global component related to block movement, we split the InSAR data into two parts based on the fault trace that was determined in Section 3.1. and subtracted the corresponding mean value, leaving the residual component that is believed to be related to fault movement. As such, the inversed result is the slip deficit of the fault during the interseismic period. It is worth noting that the GPS surface deformation is not included in the modeling due to its very low spatial resolution ( Figure 1). Moreover, introducing the GPS data caused an additional level of complexity and required dealing with spatially heterogeneous data weight.

Modeling Approach
To detect variations in slip deficit along the fault plane, we subdivided the rectangular fault plane into 4 km (along the strike direction) × 4 km (down-dip direction) patches. The slip deficit distribution model was performed using the steepest descent method (SDM) code [51], which has been successfully used to invert co-seismic slip [52], afterslip [53], aseismic slip [54,55], and interseismic coupling [56]. The fault movement inversion was implemented based on the constrained least-squares method. During the inversion, a layered earth model from CRUST1.0 [57] was utilized. A smoothing constraint was applied to the stress drop rather than to the slip distribution (e.g., [52]), since a constant stress drop is a more realistic approximation from a physical perspective. An optimal smoothing factor was determined by analyzing the trade-off curve between data misfit and model roughness (e.g., [52]). During the inversion, the slip deficit rates of the six segments of the fault were constrained according to the slip rates derived by the screw dislocation models implemented in Section 3.1. (Figure 7). Since our focus was on the spatial patterns of slip deficit, we did not attempt to analyze the temporal variations in slip deficit rate. Pure strike-slip motion was assumed during modeling.
Remote Sens. 2018, 10, x FOR PEER REVIEW 11 of 20 been successfully used to invert co-seismic slip [52], afterslip [53], aseismic slip [54,55], and interseismic coupling [56]. The fault movement inversion was implemented based on the constrained least-squares method. During the inversion, a layered earth model from CRUST1.0 [57] was utilized. A smoothing constraint was applied to the stress drop rather than to the slip distribution (e.g., [52]), since a constant stress drop is a more realistic approximation from a physical perspective. An optimal smoothing factor was determined by analyzing the trade-off curve between data misfit and model roughness (e.g., [52]). During the inversion, the slip deficit rates of the six segments of the fault were constrained according to the slip rates derived by the screw dislocation models implemented in Section 3.1. (Figure 7). Since our focus was on the spatial patterns of slip deficit, we did not attempt to analyze the temporal variations in slip deficit rate. Pure strike-slip motion was assumed during modeling.   been successfully used to invert co-seismic slip [52], afterslip [53], aseismic slip [54,55], and interseismic coupling [56]. The fault movement inversion was implemented based on the constrained least-squares method. During the inversion, a layered earth model from CRUST1.0 [57] was utilized. A smoothing constraint was applied to the stress drop rather than to the slip distribution (e.g., [52]), since a constant stress drop is a more realistic approximation from a physical perspective. An optimal smoothing factor was determined by analyzing the trade-off curve between data misfit and model roughness (e.g., [52]). During the inversion, the slip deficit rates of the six segments of the fault were constrained according to the slip rates derived by the screw dislocation models implemented in Section 3.1. (Figure 7). Since our focus was on the spatial patterns of slip deficit, we did not attempt to analyze the temporal variations in slip deficit rate. Pure strike-slip motion was assumed during modeling.   Before analyzing the slip deficit rates, a checkerboard test was performed to evaluate the reliability of the inversed slip results. Using the fault geometry mentioned in Section 3.1., we forward-  Before analyzing the slip deficit rates, a checkerboard test was performed to evaluate the reliability of the inversed slip results. Using the fault geometry mentioned in Section 3.1., we forward- Before analyzing the slip deficit rates, a checkerboard test was performed to evaluate the reliability of the inversed slip results. Using the fault geometry mentioned in Section 3.1., we forward-modeled known checkerboard slip distributions with subfault patch areas of 20 × 20 km 2 and 12 × 12 km 2 , with alternating slip values of 0 and 10 mm (Figure 10). The modeled displacements were then downsampled using the same approach as described in Section 3.2.1. to give the same input data distribution and number of displacement points. This is particularly useful in demonstrating where large gradients of slip down to zero could be resolved. The decreased resolution with depth is apparent in the two cases, and the slip of the shallow part of the fault plane (<20 km) is well resolved in the inversion (Figure 10).

Modeling Results
Remote Sens. 2018, 10, x FOR PEER REVIEW 13 of 20 modeled known checkerboard slip distributions with subfault patch areas of 20 × 20 km 2 and 12 × 12 km 2 , with alternating slip values of 0 and 10 mm ( Figure 10). The modeled displacements were then downsampled using the same approach as described in Section 3.2.1. to give the same input data distribution and number of displacement points. This is particularly useful in demonstrating where large gradients of slip down to zero could be resolved. The decreased resolution with depth is apparent in the two cases, and the slip of the shallow part of the fault plane (<20 km) is well resolved in the inversion (Figure 10). The deficit slip rates on the fault are highly variable (Figure 9a), based on which three sections can be identified. First, the western section, S2, has a high deficit rate of 7-9 mm/year near surface (0-8 km), extending for about 30 km along the strike direction. Second, the central section, S3 and S4, has relatively moderate deficit slip rates of 4-6 mm/year down to a 30 km depth. Third, the eastern section, S5 and S6, shows low deficit rates of 0-3 mm/year. Then, the fault coupling coefficient (i.e., the ratio between the deficit rate and the seismogenic loading rate) was computed to evaluate the locking degree of the fault. The coupling coefficient varies from 0 (i.e., fault slip rate matches tectonic loading rate) to 1 (fault is locked). In our case, the tectonic loading rate (also named the fault deep slip rate) was determined from the screw dislocation model mentioned in Section 3.1 (Figure 7). We observed a general pattern of eastward decreasing tectonic loading rates from 8.69 mm/year on the western segment S1 to 5.26 mm/year on the eastern segment S6 (Figures 6 and 7). Figure 9b shows the fault coupling coefficient along the ATF at 91.5-95°E. Similar to the slip deficit rate distribution, the coupling coefficient is also variable at different segments. Segments S3 and S4 are nearly fully locked, while the other segments are partially locked. Figure 3 shows a large-scale deformation pattern, which is supposed to be related to ATF movement. However, we also note four areas that are likely dominated by nontectonic signals (dashed ellipses in Figure 3). The area on the northern end of the central track (A1) shows complex fringes, which is supposed to be the effect of unwrapping error and/or residual atmospheric perturbations in the interferograms. A2 shows a subsidence signal where an alluvial fan is located. We surmise that this subsidence is caused by the consolidation of sediment deposits. A3 covers a salt lake; therefore, we speculate that salt compaction/dilatation accounts for the deformation. A4 coincides with the Lenghu oil fields, showing local subsidence that is likely related to oil production. The deficit slip rates on the fault are highly variable (Figure 9a), based on which three sections can be identified. First, the western section, S2, has a high deficit rate of 7-9 mm/year near surface (0-8 km), extending for about 30 km along the strike direction. Second, the central section, S3 and S4, has relatively moderate deficit slip rates of 4-6 mm/year down to a 30 km depth. Third, the eastern section, S5 and S6, shows low deficit rates of 0-3 mm/year. Then, the fault coupling coefficient (i.e., the ratio between the deficit rate and the seismogenic loading rate) was computed to evaluate the locking degree of the fault. The coupling coefficient varies from 0 (i.e., fault slip rate matches tectonic loading rate) to 1 (fault is locked). In our case, the tectonic loading rate (also named the fault deep slip rate) was determined from the screw dislocation model mentioned in Section 3.1 (Figure 7). We observed a general pattern of eastward decreasing tectonic loading rates from 8.69 mm/year on the western segment S1 to 5.26 mm/year on the eastern segment S6 (Figures 6 and 7). Figure 9b shows the fault coupling coefficient along the ATF at 91.5-95 • E. Similar to the slip deficit rate distribution, the coupling coefficient is also variable at different segments. Segments S3 and S4 are nearly fully locked, while the other segments are partially locked. Figure 3 shows a large-scale deformation pattern, which is supposed to be related to ATF movement. However, we also note four areas that are likely dominated by nontectonic signals (dashed ellipses in Figure 3). The area on the northern end of the central track (A1) shows complex fringes, which is supposed to be the effect of unwrapping error and/or residual atmospheric perturbations in the interferograms. A2 shows a subsidence signal where an alluvial fan is located. We surmise that this subsidence is caused by the consolidation of sediment deposits. A3 covers a salt lake; therefore, we speculate that salt compaction/dilatation accounts for the deformation. A4 coincides with the Lenghu oil fields, showing local subsidence that is likely related to oil production. Figures 3 and 4 show that the deformation of the ATF at 91.5-95 • E is constrained to occur within a width of 30 km across the fault zone. Moreover, at a distance of more than 30 km from the fault zone, the InSAR deformation profiles are satisfactorily flattened ( Figure 5), indicating that the Tarim block to the north and the Qaidam block to the south behave mostly as rigid blocks. These deformation characteristics seem to support one of the two end-member views of the Indo-Eurasian collisional models, i.e., the undeformed blocks model (e.g., [1]). Jolivet and colleagues [21] also detected the wider deformed zone of the ATF at 94 • E using ERS-1/2 and Envisat radar data covering the 1995-2006 period. They also noted that, at a distance of more than a few tens of kilometers from the fault zone, the Tarim and Qaidam blocks do not show obvious deformation. Similarly, other InSAR studies on the other segments of the ATF also show intensive deformation close to the fault zone and, conversely, no obvious deformation on the two sides of the ATF [15][16][17][18][19] except for local deformations related to small-scale faults. Therefore, based on our results and these other available InSAR results, we are inclined to favor the undeformed blocks model in the Northern Tibetan Plateau.

Comparison of Coupling and Seismicity
Seismicity records are valuable in terms of providing information that can both complement and potentially validate InSAR observations. We superimposed relocated earthquakes that occurred between 1980 and 2012 on the coupling ratio map (Figure 9b). We chose those events that were within 10 km distance of either side of the fault trace, which are believed to be related to ATF activities. Considering the magnitude completeness threshold of the seismic catalogue, we only plotted earthquakes larger than M 3.0.
A comparison between the spatial distribution of the earthquakes and the fault coupling coefficient revealed the following: First, most of the earthquakes occurred on segments S1, S4, S5, and S6. No seismicity was detected at depths of more than 40 km, implying a transition to the ductile zone of the fault. Second, considering only the coupling results at the shallow part (<20 km) of the fault plane, which are relatively reliable according to the checkerboard tests, the coupling distribution shows variations along the strike. Segments S3 and S4 are nearly fully locked, with the eastern section of segment S4 having some earthquakes at depths shallower than 10 km. In general, a high coupling value indicates regions where the fault is not slipping and is identified as an asperity [58]. Thus, there may be a large asperity on the fault plane of segments S3 and S4, extending more than 100 km at depths between 8 and 12 km. This asperity may represent the nucleation site or rupture area of future large earthquakes (the white line-enclosed area with a question mark shown in Figure 9b). Meanwhile, a smaller potential asperity extending approximately 20 km at depths shallower than 10 km exists on segment S2. In contrast, coupling values on segments S1, S5, and S6 are not high. A large number of small earthquakes occurred along these segments, releasing some accumulated energy and leading to a relatively low risk of strong earthquakes. However, the occurrence of moderate-sized earthquakes cannot be ruled out. Usually, a geometrical or compositional perturbation and/or heterogeneity in physical properties existing on a fault plane at depth generate an asperity, controlling the variable slip distribution and seismic activity (e.g., [59]). Such spatial variability of interseismic coupling is visible on other major continental strike-slip faults, e.g., the San Andreas Fault [60] and the North Anatolian Fault [61].

Seismic Potential of the ATF at 91.5-95 • E
On the basis of the slip deficit rate distribution shown in Figure 9a, we performed a first-order calculation to compute moment accumulation (M0) using the relation: where µ is the shear modulus and assumed to be 30 GPa, U is the slip-deficit rate from the inversion, and A is the area of the fault patch. According to the paleo-earthquake record, the latest earthquake event near the town of Akesai occurred in 665 ± 40 years BP [10]. Over a 665-year period, the slip deficit would be equivalent to an Mw 7.91 ± 0.01 earthquake. In other words, the estimated moment deficit is equivalent to an Mw 7.91 ± 0.01 earthquake if the entire segment of the ATF studied in this paper ruptures in a single earthquake. Although a large earthquake in this segment is a possibility, the deficit may also be made up by several moderately-sized events or by accelerated afterslip subsequent to large earthquakes on neighboring locked patches.

Present-Day Slip Rate and Locking Depth Variations along the Entire ATF from GPS and InSAR
The present-day slip rate and locking depth along different segments of the ATF have been studied by various researchers using geodetic observations, especially by GPS and InSAR techniques ( Figure 11; Table 3). The slip rate of the western segment at 78-81 • E is constrained between 2 and 5 mm/year [10,15,16]. However, at 84 • E, the slip rate increases to~11 mm/year, and decreases linearly eastward to~4 mm/year at 96 • E [10,[17][18][19][20][21][22][23][24][25][26][27][28]. The large difference of slip rates between 78-81 • E and 84 • E east is expected, considering the bending of the ATF in this segment (Figure 1a). The variation of the slip rate along the ATF may be caused by the heterogeneity in physical properties of the fault plane. In contrast, the locking depth keeps nearly constant at~15 km through the entire ATF. The large locking depth indicates that the ATF is capable of strong events, which is consistent with the paleoseismic and historical earthquake records [10,12].
Remote Sens. 2018, 10, x FOR PEER REVIEW 15 of 20 deficit would be equivalent to an Mw 7.91 ± 0.01 earthquake. In other words, the estimated moment deficit is equivalent to an Mw 7.91 ± 0.01 earthquake if the entire segment of the ATF studied in this paper ruptures in a single earthquake. Although a large earthquake in this segment is a possibility, the deficit may also be made up by several moderately-sized events or by accelerated afterslip subsequent to large earthquakes on neighboring locked patches.

Present-Day Slip Rate and Locking Depth Variations along the Entire ATF from GPS and InSAR
The present-day slip rate and locking depth along different segments of the ATF have been studied by various researchers using geodetic observations, especially by GPS and InSAR techniques ( Figure 11; Table 3). The slip rate of the western segment at 78-81°E is constrained between 2 and 5 mm/year [10,15,16]. However, at 84°E, the slip rate increases to ~11 mm/year, and decreases linearly eastward to ~4 mm/year at 96°E [10,[17][18][19][20][21][22][23][24][25][26][27][28]. The large difference of slip rates between 78-81°E and 84°E east is expected, considering the bending of the ATF in this segment (Figure 1a). The variation of the slip rate along the ATF may be caused by the heterogeneity in physical properties of the fault plane. In contrast, the locking depth keeps nearly constant at ~15 km through the entire ATF. The large locking depth indicates that the ATF is capable of strong events, which is consistent with the paleoseismic and historical earthquake records [10,12]. Figure 11. Variations of (a) slip rate and (b) locking depth along the ATF estimated from GPS (red diamonds) and InSAR (blue circles) as reported in previous studies. The uncertainties are shown at the 1σ level. The estimates in this study are shown as black stars. (A) and (D) are from [16]; (B) is from [15]; (C), (H), (R), and (U) are from [10]; (E) is from [18]; (F) is from [19]; (G) is from [17]; (J) is from [23]; (I) and (Q) are from [25]; (K) is from [22]; (L) is from [24]; (M) is from [27]; (N) is from [20]; (O) is from [26]; (P), (T), and (V) are from [28]; and (S) is from [21].
In this study, we achieved the present-day deformation velocity of the ATF at 91.5-95 • E using the InSAR technique (Figure 3). Constrained by five InSAR deformation profiles, slip rates and locking depths were estimated using a screw dislocation model ( Figure 6). The resultant averaged slip rate and locking depth of the ATF at 91.5-95 • E are 6.4 ± 1.6 mm/year and 15 ± 4 km, respectively (stars in Figure 11). We note that our slip rate is comparable with the 8-10 mm/year slip rate estimated at 94 • E [21] and the 4.5 ± 0.8 mm/year slip rate estimated at 94.6 • E [28] (Figure 11a). Moreover, our locking depth shows good agreement with published results [10,21,28] (Figure 11b). Table 3. Slip rate and locking depth along the ATF estimated from GPS and InSAR measurements as reported in previous studies. The uncertainties are shown at 1σ level. The number listed in the first column is the same as the list in Figure 11. NA stands for not available.

Conclusions
Estimating the fault locking depth and interseismic slip rate is essential for assessing the seismic hazard of an area. With the development of the SAR satellite systems, it is possible to detect fault movement deformation at different phases of a seismic cycle with high spatiotemporal resolution, and standardized methods have been established in recent years (e.g., [62,63]). High spatial resolution deformation fields can also help to detect asperities that may trigger future earthquake events (e.g., [55]). In this study, we obtained the interseismic deformation field of the eastern ATF at 91.5-95 • E using the InSAR technique. The InSAR-derived deformation shows good agreement with GPS observations. Interseismic slip deficit rate distribution was resolved by the InSAR deformation. The slip deficit is capable of an upper bound of Mw 7.9 earthquake considering the elapsed time of the latest M 7 event, providing insight into the seismogenic potential of the ATF. Based on coupling coefficient and seismicity distribution, we delineated two potential asperities on the fault plane where future earthquakes may nucleate or rupture. This study provides an example of seismic potential assessments based on high spatial resolution geodetic observations.
The existing literature shows that the present-day slip rate of the ATF cannot be constrained as a single value. The slip rate decreases eastward from~11 mm/year at 84 • E to~4 mm/year at 96 • E.
For the segment between 91.5 • E and 95 • E, our InSAR-derived slip rate is 6.4 mm/year, which supports the existing eastward decreasing trend. Although different slip rates are shown in different segments of the ATF, the locking depths stay the same, i.e.,~15 km at the entire ATF. Our InSAR-derived locking depth at 91.5-95 • E is 15 km and agrees well with the existing results in the literature. All the InSAR deformation fields of the ATF show that the Tarim and Qaidam blocks behave mostly as rigid blocks, supporting the undeformed blocks model in the northern Tibetan plateau, i.e., one of the two end-member views of the Indo-Eurasian collisional models.