A Fine Velocity and Strain Rate Field of Present-Day Crustal Motion of the Northeastern Tibetan Plateau Inverted Jointly by InSAR and GPS

Interferometric synthetic aperture radar (InSAR) data from 6 Envisat ASAR descending tracks; spanning the 2003–2010 period; was used to measure interseismic strain accumulation across the Northeastern Tibetan Plateau. Mean line-of-sight (LOS) ratemaps are computed by stacking atmospheric-corrected and orbital-corrected interferograms. The ratemaps from one track with different atmospheric-corrected results or two parallel; partially overlapping tracks; show a consistent pattern of left-lateral motion across the fault; which demonstrates the MERIS and ECMWF atmospheric correction works satisfactorily for small stain measurement of this region; even with a limited number of interferograms. By combining the measurements of InSAR and GPS; a fine crustal deformation velocity and strain rate field was estimated on discrete points with irregular density depending on the fault location; which revealed that the present-day slip rate on the Haiyuan fault system varies little from west to east. A change (2–3 mm/year) in line-of-sight (LOS) deformation rate across the fault is observed from the Jinqianghe segment to its eastern end. Inversion from the cross-fault InSAR profiles gave a shallow locking depth of 3–6 km on the main rupture of the 1920 earthquake. We therefore infer that the middle-lower part of the seismogenic layer on the 1920 rupture is not yet fully locked since the 1920 large earthquake. Benefit from high spatial resolution InSAR data; a low strain accumulation zone with high strain rates on its two ends was detected; which corresponds to the creeping segment; i.e., the Laohushan fault segment. Contrary to the previous knowledge of squeezing structure; an abnormal tension zone is disclosed from the direction map of principal stress; which is consistent with the recent geological study. The distribution of principal stress also showed that the expanding frontier of the northeastern plateau has crossed the Liupan Shan fault zone; even arrived at the northeast area of the Xiaoguan Shan. This result agrees with the deep seismic reflection profile.


Introduction
Surface velocity and strain fields provide an important constraint on geodynamic models of tectonic deformation as well as in the assessment of earthquake hazard [1]. Such velocity fields may be determined using geodetic measurements such as GPS, but GPS observations are too sparse to describe deformation spatial variations in the near-fault area. Increasing Synthetic Aperture Radar (SAR) coverage allows SAR interferometry (InSAR) to be used as a complementary or alternative means of measuring deformation of the Earth's surface. InSAR has been successfully used to determine coseismic deformation of many earthquakes (e.g., the 2003 Bam earthquake [2], 2011 Yushu earthquake [3], and 2015 Nepal earthquakes [4]) and post-seismic deformation (e.g., Atzori et al. 2008 [5] and Gonzalez-Ortega et al. 2014 [6]), as well as deformation related to earthquake swarms (e.g., Kyriakopoulos et al. 2013 [7] and Wicks et al. 2011 [8]). In recent years, InSAR has begun to perform well in measuring the rate and distribution of interseismic strain accumulation with the increasing satellite data volume (e.g., in Northwestern Tibet [9,10], the San Andreas fault system [11][12][13], the Anatolian Fault system [1,14,15], and the Haiyuan fault system in the northeastern margin of the Tibetan Plateau [16][17][18][19]). Satellite InSAR data provide us a chance to give an insight into the near-fault crustal displacement caused by the fault motion, particularly for the shallow aseismic fault slip. It has provided important information on imaging the spatial variation of interseismic coupling on the fault plane in the upper crust. As well as benefitting from its ability to characterize both the long-and short-wavelength deformation with a high accuracy, InSAR also provides new knowledge for improving our understanding the processes of crustal deformation caused by interaction of faults or blocks in some fault transferring zones and the block boundary belts.
In this paper, we focus on the northeastern margin of the Tibet Plateau, where the crustal shortening and lateral extrusion caused by the northeastward growth of the Tibet Plateau is accommodated by the left-lateral strike-slip Haiyuan fault and a series of arcuate structures in its north, and the thrusting Liupan Shan fault and a few folds and thrust faults in its east ( Figure 1). Two large earthquakes occurred on the Haiyuan fault system in the last century. One (Mw > 8.0) with a 220 km long surface rupture occurred on the eastern part of the Haiyuan fault in 1920, resulting in 220,000 deaths [20]. The other one is the 1927 Mw 8.0 Gulang earthquake ruptured south-dipping thrusts situated at the south-eastern end of the Qilian Shan [21]. Between these two ruptured segments, a seismic gap with the 260 km length is identified as the Tianzhu gap with a high seismic potential [17,21]. Here, we combined the InSAR and GPS to give refined estimates of the amplitude and distribution of elastic strain accumulation on those boundary faults with a high spatial resolution. We show how deformation is distributed locally along the faults in this region, what their current strain accumulation rate is. Finally, we discussed the present-day fault kinematic behavior of the 1920 rupture, and the tectonic response of the Tibet-Ordos transition zone to the plateau expansion.

Tectonic Setting
In the northeastern margin of the Tibetan Plateau, the Haiyuan fault system is a major active tectonic feature, connecting the seismically active Qilian Shan in the west and the tectonically active Liupan Shan in the east, the latter abutting the relatively stable Ordos block ( Figure 1). It is dominated by left-lateral strike slip faulting, which probably began in the late Pliocene or early Pleistocene, followed by late-stage folding and thrust faulting [23]. At its eastern end, left lateral slip on the N65 • W striking Haiyuan fault zone is transferred into shortening on the generally north trending structures in the Liupan Shan. Three arcuate zones of both strike-slip and thrust faults with associated ramp anticlines, lie about 40-170 km north and northeast of the Haiyuan fault zone. From south to north, the individual structures that comprise this arcuate system are the Tianjin Shan-Mibo Shan, Yanton Shan, and Niushou Shan-Daluo Shan fault zones. Tectonic activity within these areas is generally mild in comparison with that in the Haiyuan structural zone [24].
The Holocene slip rate of the Haiyuan fault system, determined by measuring the offsets and ages of morphological markers, seems to vary from west to east. The left-lateral strike-slip rates in the west have been estimated as 19 ± 5 mm/year on the Lenglongling fault segment, and 12 ± 4 mm/year on the Laohu Shan fault further east [25,26]. However, He et al. [27,28] and Zheng et al. [29] suggest much slower slip rates of 4-5 mm/year on the Lenglongling fault and 4.4 ± 0.4 mm/year on the Laohu Shan fault. The Holocene slip rate estimates along the length of the 1920 rupture (including the Xihuan-Nanhua segment) in the east are 8 ± 2 mm/year [20] and 4.5 ± 2 mm/year [30], which are consistent with an average Quaternary slip rate of 5-10 mm/year [22]. The rates of deformation in the Tianjing Shan and Mibo Shan are about 1.5-2.7 mm/year, notably less than that in the south [24]. GPS profiles across the Qilian-Haiyuan fault show that the left-lateral strike-slip rates and shortening rates from the west to the east, are 4.0 ± 1.0 mm/year and 5.6 ± 1.5 mm/year along the Lenglongling segment, 3.6 ± 1.4 mm/year and 2.3 ± 1.2 mm/year along the Maomao Shan and Laohu Shan segment, and 4.2 ± 1.5 mm/year and 2.6 ± 1.5 mm/year along the 1920 rupture respectively [29]. InSAR results from ERS and Envisat ASAR data give apparent slip rates of 4.2-8 mm/year [16] and 5.3 ± 1.0 mm/year [17] near the junction between the 1920 rupture and Tianzhu gap (the Maomao Shan-Laohu Shan segment). Additionally, they highlight a strong strain concentration in the Laohu Shan fault zone and suggest the presence of creep at shallow depths. Furthermore, the smoothed InSAR time series reveals that the creeping rate on the Laohu Shan segment accelerated in 2007 [18].

Data
We collected the Envisat SAR stripmap images for 6 descending tracks, which cover the most part of the northeastern margin of the Tibetan Plateau ( Figure 1). Interferometric data for each track spans the 2003-2010 period with 30 more acquisitions except for the track 247 (12 acquisitions). Interferograms with baselines of less than 200 m were produced for each track by using JPL/Caltech ROI_PAC software. A filled 3 arc sec (90 m) resolution Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) [31] was used to remove the topographic contribution to the interferometric phase changes. As we are only interested in long-wavelength features, the data were multilooked, reducing the sampling interval to~320 m in order to improve coherence and retain the tectonic signal near the fault. The coherence threshold of 0.15 was used to mask out the decoherenced area for each interferograms. Then interferograms were unwrapped using the branch-cut method [32]. The resulting unwrapped interferograms will be selected to construct a chain for stacking in the InSAR ratemap estimation (Section 3.1). GPS data we used is from Gan et al. [22], which is with respect to the Eurasia-fixed Reference Frame.

Stacking InSAR with Atmospheric-Corrected Interferograms
The detection of slow interseismic motion using InSAR requires a combination of multiple radar images over an extended time period to reduce the uncorrelated atmospheric effects and improve the signal-to-noise ratio. The total tropospheric delay is often divided into a stratified delay and a turbulent delay for InSAR applications [33]. Due to the high variability of atmospheric turbulence, the spatial pattern of the turbulent delay is mostly random on each acquisition date and can be removed efficiently by stacking interferograms [34] or through InSAR time series analysis [35]. Variations in atmospheric stratification have a first order effect, with a single path LOS delay reaching ∼10 cm/km in some areas [36]. Furthermore, unlike turbulent effects, stratified delays are not random in time but represent seasonal fluctuations and are therefore difficult to remove by stacking or temporal smoothing compared to turbulence, even when working with a large InSAR dataset. When we retrieve subtle ground displacements (e.g., millimeters per year on the Haiyuan fault), it is necessary to correct interferograms for the stratified tropospheric signals.
Given strong coupling between the linear orbital ramp, long-wavelength component of atmospheric error and long-wavelength tectonic deformation, external independent atmospheric data (ECMWF and MERIS) were used to correct the stratified tropospheric signals firstly, which is verified to be useful for correcting InSAR atmospheric effect [14,37]. Then a linear orbital trend was fitted and removed from the atmospheric-corrected interferograms, using only phase measurements 30 km or further from the fault in the north to avoid the effect of the near-fault gradient in ground deformation. After the atmospheric and orbital corrections, the residual interferograms from independent chains of small-baseline dataset ( Figure S1) were stacked with a total timespan to create a ratemap. Only those pixels which are coherent in all interferograms were averaged with equal weight. The flowchart for retrieving an InSAR ratemap is showed in Figure 2. verified to be useful for correcting InSAR atmospheric effect [14,37]. Then a linear orbital trend was fitted and removed from the atmospheric-corrected interferograms, using only phase measurements 30 km or further from the fault in the north to avoid the effect of the near-fault gradient in ground deformation. After the atmospheric and orbital corrections, the residual interferograms from independent chains of small-baseline dataset ( Figure S1) were stacked with a total timespan to create a ratemap. Only those pixels which are coherent in all interferograms were averaged with equal weight. The flowchart for retrieving an InSAR ratemap is showed in Figure 2.

Velocity and Strain-Rate Field Inversion from InSAR and GPS
A variety of algorithms have been developed to invert geodetic data for strains. The traditional subnetwork methods (e.g. [38,39]) usually generate discontinuities of the strain estimates at the subnetwork boundaries. The geostatistical methods can tackle the discontinuity problem. However, they assume that the deformation field is isotropic and homogeneous, and the accuracy and reliability often depends on the spatial distribution of the data (e.g. [1,[40][41][42][43]). The physics-based models were developed to estimate the elastic and inelastic strain fields from geodetic data (e.g., [44,45]

Velocity and Strain-Rate Field Inversion from InSAR and GPS
A variety of algorithms have been developed to invert geodetic data for strains. The traditional subnetwork methods (e.g., [38,39]) usually generate discontinuities of the strain estimates at the subnetwork boundaries. The geostatistical methods can tackle the discontinuity problem. However, they assume that the deformation field is isotropic and homogeneous, and the accuracy and reliability often depends on the spatial distribution of the data (e.g., [1,[40][41][42][43]). The physics-based models were developed to estimate the elastic and inelastic strain fields from geodetic data (e.g., [44,45]). Shen et al. [46][47][48] developed a series of algorithms to model strains with reasonable weighting and smoothing function optimally determined from observation data, which works well with unevenly distributed data, and can deal with the nonelastic strain accumulation caused by fault creeping in the upper crust. The surface deformation field caused by the strike-slip Haiyuan fault system, is not homogeneous and isotropic, and would changes cross the fault boundary. Previous InSAR results showed a creeping segment (the Laohu Shan segment) on the Haiyuan fault [17], thus the method in Shen et al. [48] was used in our study. We realized Shen's method in a spherical coordinate system due to our big study area (500 × 500 km). Considering the high spatial coverage of InSAR velocity map, an irregular grid that discretized the study area based on the location of fault, is used in InSAR and GPS inversion, which allow us capture the strain variations near the fault.
For each node of the grid, its unknown vector can be related to the 3D velocity of its surrounding GPS observations by the following model: where, u GPS is the vector of GPS observation, l is the unknown vector of the displacement, strain, and rotation on each node of the grid, and A GPS is the design matrix. This observation equation can be expanded as: where, u φ u θ u r is the vector of GPS observation; U φ U θ U r , ω r ω φ ω θ and ε φφ ε φθ ε θθ is the 3D deformation velocity, rotation and strain components on each discrete point; φ 0 and θ 0 is Latitude and longitude coordinates of the node; and r 0 is the earth radius.
For InSAR observations, the observation equation can be established as: where, u los is the vector of InSAR observation, A los is the design matrix. InSAR LOS velocity can be written by using 3D deformation components ( Figure S2) as follow: Then the design matrix A los can be calculated by A GPS as: where α is the azimuth angle of InSAR, and θ is the incidence angle. Finally, the joint observation equations for GPS and InSAR can be formed as: A least-squares solution can be resolved in the form of where, C is the weighting factor, which is constructed by the covariance matrix C multiplying a weighting function G. G accounts for the distance and spatial coverage dependency between the unknown point and observation points (see Shen et al. [48] for details)

Construction of the InSAR Rate Map
For each Envisat ASAR track we collected, we constructed an interferometric stack composed of several chains of small-baseline interferograms, then applied the stacking InSAR method with atmospheric and orbital correction described in Section 3.1, to construct a deformation ratemap. Track 18 was chosen as an example to show the accuracy and reliability assessment of our results. Based on the baseline distribution of the track 18 ( Figure S1), we constructed 3 chains of small-baseline interferograms (including 5 interferograms), which are completely composed of largely cloud-free SAR acquisitions, in order to make sure that the atmospheric contribution in each interferogram can be removed correctly using MERIS data (300 m spatial resolution). Following the steps mentioned above, A MERIS-corrected ratemap is achieved (Figure 3a). A rapid increase in the line-of-sight rate of deformation (∼2 mm/year) in a relatively narrow region across the Haiyuan fault can be seen, qualitatively consistent with left-lateral slip. There is no clear gradient in displacement rate related to tectonic signal on other faults north of the Haiyuan fault. GPS velocities were also projected into the local InSAR LOS direction and agree with the InSAR observations within the errors of the two measurements (see Figure 4). The RMS between the InSAR result and GPS is 0.6 mm/year.
In order to further validate our results using the MERIS corrections for atmospheric path delays, we also examined the ECMWF numerical weather model results for the region. Applying the method of Walters et al. [14], we produced atmospheric delay maps from ECMWF model results for each SAR acquisition. We compared the wet path delay maps derived from MERIS and ECMWF for the times of all the SAR acquisitions. A strong correlation was found between the ECMWF and MERIS-derived wet delay maps for 10 cloud-free acquisitions ( Figure S3) in this area, with a mean correlation coefficient of 0.79 (ranging from 0.52 to 0.92), and a mean RMS misfit of 1 cm. There is also a good correlation for partially and largely cloudy acquisitions, although with a more limited number of data points. This suggests that the ECMWF atmospheric delay is potentially an alternative to MERIS data for correcting interferograms for atmospheric effects in this region.
In order to test the reliability and usefulness of the ECMWF data, we selected 13 interferograms (covering a total cumulative time of 40 years, see Figure S1) not included in the MERIS chains, and stacked 7 chains of these interferograms using the ECMWF correction for atmospheric effects.  Following the steps mentioned above, we processed 12 Envisat ASAR acquisitions from descending track 247, which partially overlaps with track 18. A total of 36 interferograms with baselines of less than 200 m were produced, but only 3 of these showed good enough coherence. Subsequent analysis was based on these 3 interferograms alone. No cloud-free MERIS data is available for the acquisitions on track 247. But because the atmospheric correction based on the ECMWF model works well for track 18, we applied it to the interferograms for track 247 as well. After removing the atmospheric delays in this way and also removing an orbital ramp, we stacked the interferograms to create a ratemap for this track (Figure 5a). Even though only 3 interferograms from track 247 were used, the resulting ratemap agrees very well with the results from track 18. We compared the rates of range change in the overlapping area between the two swaths ( Figure 5b). The two independent results show a very similar pattern with an offset of 2 mm/year in the LOS direction. This suggests that it is the quality of the interferograms that matters, and it is possible to measure slow ground motion with a very limited number of interferograms in this area, provided that appropriate corrections can be made for atmospheric effects. In order to further validate our results using the MERIS corrections for atmospheric path delays, we also examined the ECMWF numerical weather model results for the region. Applying the method of Walters et al. [14], we produced atmospheric delay maps from ECMWF model results for each SAR acquisition. We compared the wet path delay maps derived from MERIS and ECMWF for the times of all the SAR acquisitions. A strong correlation was found between the ECMWF and MERIS-derived wet delay maps for 10 cloud-free acquisitions ( Figure S3) in this area, with a mean correlation coefficient of 0.79 (ranging from 0.52 to 0.92), and a mean RMS misfit of 1 cm. There is also a good correlation for partially and largely cloudy acquisitions, although with a more limited number of data points. This suggests that the ECMWF atmospheric delay is potentially an alternative to MERIS data for correcting interferograms for atmospheric effects in this region.
In order to test the reliability and usefulness of the ECMWF data, we selected 13 interferograms (covering a total cumulative time of 40 years, see Figure S1) not included in the MERIS chains, and stacked 7 chains of these interferograms using the ECMWF correction for atmospheric effects. The ECMWF-corrected result (Figure 3b) agrees well with the MERIS-corrected result.
Following the steps mentioned above, we processed 12 Envisat ASAR acquisitions from descending track 247, which partially overlaps with track 18. A total of 36 interferograms with baselines of less than 200 m were produced, but only 3 of these showed good enough coherence. Subsequent analysis was based on these 3 interferograms alone. No cloud-free MERIS data is available for the acquisitions on track 247. But because the atmospheric correction based on the ECMWF model works well for track 18, we applied it to the interferograms for track 247 as well. After removing the atmospheric delays in this way and also removing an orbital ramp, we stacked the interferograms to create a ratemap for this track (Figure 5a). Even though only 3 interferograms from track 247 were used, the resulting ratemap agrees very well with the results from track 18. We compared the rates of range change in the overlapping area between the two swaths ( Figure 5b). The two independent results show a very similar pattern with an offset of 2 mm/yr in the LOS direction. This suggests that it is the quality of the interferograms that matters, and it is possible to measure slow ground motion with a very limited number of interferograms in this area, provided that appropriate corrections can be made for atmospheric effects. From the experiments above, we can see that the stacking method with atmospheric and orbital correction proposed in this study works well in the Haiyuan fault zone, which is shown by: 1) the consistency between the MERIS-corrected and ECMWF-corrected ratemap (Figure 3), 2) consistency between the results from two adjacent tracks (Figure 5), and 3) consistency between InSAR and GPS ( Figure 4). Therefore, applying this method to the other 5 tracks, we retrieved an InSAR ratemap for the most part of the Northeastern Tibetan area (Figure 6). From the experiments above, we can see that the stacking method with atmospheric and orbital correction proposed in this study works well in the Haiyuan fault zone, which is shown by: (1) the consistency between the MERIS-corrected and ECMWF-corrected ratemap (Figure 3), (2) consistency between the results from two adjacent tracks ( Figure 5), and (3) consistency between InSAR and GPS ( Figure 4). Therefore, applying this method to the other 5 tracks, we retrieved an InSAR ratemap for the most part of the Northeastern Tibetan area (Figure 6).

Inversion Results from InSAR and GPS
In order to obtain the detailed variations of deformation velocity near the faults, which may be observed by the high-resolution InSAR data, the whole Northeastern Tibetan area was meshed into an irregular grid by using the free software DistMesh [49]. Some points on the Haiyuan and other faults in this region were extracted to be used as the fixed points in discretization. A nominal mesh spacing of 0.05° was used to define the mesh density near the fault. An unstructured mesh was generated with different mesh density varying with the distance from the fault, i.e. dense points near the fault and sparse points far from the fault ( Figure S4). Applying the algorithm presented in the section 3, the horizontal velocities and strain rates on each node of the mesh were inverted from InSAR and GPS data. To avoid the effect from the nonelastic strain accumulation caused by the fault creep, we set the creeping section (103.68° E, 37.11° N to 104.15° E, 37.00° N) of the Haiyuan fault (Laohu Shan segment) identified by Jolivet et al. [17] as the barrier section. Figure 7 shows the resulting horizontal strain rate field and principal compressive /tensile stress direction map, and Figure S5 shows the resulting velocity field. In general, the Northeastern Tibetan Plateau rotates clockwise ( Figure S5), and the velocity decreases from south to north. Compared with the other faults in the north, the Haiyuan fault is the most active one, which demonstrates it is still an important boundary fault for the N-NW expansion of the Tibetan Plateau. This is in agreement with the geological research [24]. The difference between our study and geological studies [25,26] is that the present-day slip rate on the Haiyuan fault system varies little from west to east. A change (2-3 mm/yr) in line-of-sight (LOS) deformation rate across the fault is observed from Jiangqianghe segment to its eastern end, where the northward expansion of the Plateau is transferring into

Inversion Results from InSAR and GPS
In order to obtain the detailed variations of deformation velocity near the faults, which may be observed by the high-resolution InSAR data, the whole Northeastern Tibetan area was meshed into an irregular grid by using the free software DistMesh [49]. Some points on the Haiyuan and other faults in this region were extracted to be used as the fixed points in discretization. A nominal mesh spacing of 0.05 • was used to define the mesh density near the fault. An unstructured mesh was generated with different mesh density varying with the distance from the fault, i.e., dense points near the fault and sparse points far from the fault ( Figure S4). Applying the algorithm presented in the Section 3, the horizontal velocities and strain rates on each node of the mesh were inverted from InSAR and GPS data. To avoid the effect from the nonelastic strain accumulation caused by the fault creep, we set the creeping section (103.68 • E, 37.11 • N to 104.15 • E, 37.00 • N) of the Haiyuan fault (Laohu Shan segment) identified by Jolivet et al. [17] as the barrier section. Figure 7 shows the resulting horizontal strain rate field and principal compressive /tensile stress direction map, and Figure S5 shows the resulting velocity field. In general, the Northeastern Tibetan Plateau rotates clockwise ( Figure S5), and the velocity decreases from south to north. Compared with the other faults in the north, the Haiyuan fault is the most active one, which demonstrates it is still an important boundary fault for the N-NW expansion of the Tibetan Plateau. This is in agreement with the geological research [24]. The difference between our study and geological studies [25,26] is that the present-day slip rate on the Haiyuan fault system varies little from west to east. A change (2-3 mm/year) in line-of-sight (LOS) deformation rate across the fault is observed from Jiangqianghe segment to its eastern end, where the northward expansion of the Plateau is transferring into shortening along the Liupan Shan fault due to the obstruction of the Alxa block in the north and Ordos block in the east.  From the horizontal strain rate field (Figure 7a), strain accumulation is strongly localized on the left-lateral strike-slip Haiyuan fault zone. This is similar to the discovery of previous studies on the Karakoram fault [10] and the Anatolia fault [1]. Wright et al. [50] suggests that strain accumulation in inter-seismic stage is a major feature of large strike-slip faults. High strain rate regions also are observed to follow the other faults in the north, although the magnitude is smaller than that of strain rate on the Haiyuan fault. Benefit from high spatial resolution InSAR ratemap, some detail strain rate variation near the fault can be detected, which can not be observed by GPS only as shown in Li et al. [51]. In Figure 8, an interesting phenomenon caused by the fault creep is revealed clearly by our fine strain rate field. The creeping segment reveal by our InSAR and Jolivet et al. [17] located at a low strain accumulation zone, instead the areas next to its two ends are holding a high strain rate. From the horizontal strain rate field (Figure 7a), strain accumulation is strongly localized on the left-lateral strike-slip Haiyuan fault zone. This is similar to the discovery of previous studies on the Karakoram fault [10] and the Anatolia fault [1]. Wright et al. [50] suggests that strain accumulation in inter-seismic stage is a major feature of large strike-slip faults. High strain rate regions also are observed to follow the other faults in the north, although the magnitude is smaller than that of strain rate on the Haiyuan fault. Benefit from high spatial resolution InSAR ratemap, some detail strain rate variation near the fault can be detected, which can not be observed by GPS only as shown in Li et al. [51]. In Figure 8, an interesting phenomenon caused by the fault creep is revealed clearly by our fine strain rate field. The creeping segment reveal by our InSAR and Jolivet et al. [17] located at a low strain accumulation zone, instead the areas next to its two ends are holding a high strain rate. observed to follow the other faults in the north, although the magnitude is smaller than that of strain rate on the Haiyuan fault. Benefit from high spatial resolution InSAR ratemap, some detail strain rate variation near the fault can be detected, which can not be observed by GPS only as shown in Li et al. [51]. In Figure 8, an interesting phenomenon caused by the fault creep is revealed clearly by our fine strain rate field. The creeping segment reveal by our InSAR and Jolivet et al. [17] located at a low strain accumulation zone, instead the areas next to its two ends are holding a high strain rate. From the direction of principal stress (Figure 7b), the Haiyuan fault system is dominated by a strike-slip motion with a relatively small thrust component (please note that the pink E-W strike-slip components is covered by the black thrusting components), which is consistent with the focal mechanism of the moderate earthquakes happened on it in recent decades ( Figure S6). An abnormal tension zone can be see along the Daluo Shan fault, which is contrary to the previous knowledge of squeezing structure. Recent geological study [52] confirm that there exists a normal faulting component in the junction between the Niushou Shan fault and Daluo Shan fault. The regional seismicity map also shows that a Mw 5.3 earthquake with normal faulting occurred nearby on 10 August 1987 ( Figure S6). We also note that the southern part of the Liupan Shan fault shows a dominating extension that is out of our expectation. That is because local GPS observation dominates the inversion result in absence of high-coherenced InSAR data in this region, where one of GPS velocity in the west of the fault is small than the ones in the east ( Figure S7), thus caused a tensional stress zone.

The Present-Day Kinematics for the Rupture of the 1920 Earthquake
To study the present-day fault movement of the 1920 rupture, we constructed fault-perpendicular profiles across the atmospheric-corrected ratemaps from the track 18 covered the 1920 rupture. We projected all LOS velocities onto this profile, then a conventional screw dislocation arctan model [53] was used to invert for the best-fitting slip rate and locking depth. The minimum misfit corresponds to a best-fitting model with slip rate of 5.9 mm/year and locking depth of 3 km shown by the dashed red line in Figure 4a for the MERIS-corrected ratemap. The inversion for the ECMWF-corrected ratemap gives a similar result (Figure 4b: a slip rate of 5 mm/year with a locking depth of 6 km). Our InSAR results are consistent with the GPS velocity of 4.2 ± 1.5 mm/year given by Zheng et al. [29], and Holocene slip rate of 4.5 ± 2 mm/year given by Li et al. [30].
A shallow locking depth is often linked to creep at shallow depths. This kind of creep might be disconnected from the steadily slipping section below the seismogenic depth [16]. Such creeping behavior has been observed on the western end of the 1920 rupture, i.e., the Laohu Shan segment [16][17][18]. Our inversions show a shallow locking depth (3-6 km) for the main rupture zone of the 1920 earthquake. Our InSAR profiles (Figures 4 and 5) also show that the LOS velocities vary rapidly in a relative narrow zone across the fault, but not as sharp as in the Laohu Shan segment (~5 km) showed in Jolivet et al. [17]. We therefore infer that the middle-lower part of the seismogenic layer on the 1920 rupture is not yet fully locked since the 1920 large earthquake, which may explain why no moderate earthquakes have ever occurred here since 1920. The question remains whether it is caused by the creep or by the post-seismic afterslip, which may be answered by more observations in future decades.

Expansion Frontier of the Northeastern Tibet Plateau
Previous studies suggested that the northeastern Tibet is expanding laterally since the late Cenozoic (20-10 Ma), which is partly accommodated by the left-lateral strike-slip along the East Kunlun and Haiyuan faults and crustal shortening of the Tibet-Ordos transition zone [54]. It is still controversial about the route and scale of lateral motion caused by the outward growth of the Northeastern Tibetan Plateau. Some studies suggested that the Plateau is extruding laterally along the EW trending Qinling tectonic zone [55,56], but some results argued that the ENE lateral expansion of the plateau has being penetrating into the southern margin of Ordos Block or even its interior due to spatial differences in the internal lithosphere structure of the Ordos [57][58][59]. From the distribution of principal stress ( Figure 9) we achieved here, the expanding front of the northeastern plateau has crossed the Liupan Shan fault zone, even arrived at the northeast area of the Xiaoguan Shan. This is also revealed by deep seismic reflection profile. The deep crustal "architecture" suggests that the Xiaoguan Shan to the east of the Liupan Shan fault zone has been affected by the northeastward growth of the plateau [60].
Remote Sens. 2019, 11, x FOR PEER REVIEW 14 of 18 to spatial differences in the internal lithosphere structure of the Ordos [57][58][59]. From the distribution of principal stress ( Figure 9) we achieved here, the expanding front of the northeastern plateau has crossed the Liupan Shan fault zone, even arrived at the northeast area of the Xiaoguan Shan. This is also revealed by deep seismic reflection profile. The deep crustal "architecture" suggests that the Xiaoguan Shan to the east of the Liupan Shan fault zone has been affected by the northeastward growth of the plateau [60].

Conclusions
We used InSAR, combined with two different estimates of atmospheric path delays, to estimate the interseismic strain accumulation across the Northeastern Tibetan Plateau. The results show that the MERIS-based and ECMWF-based corrections both work well in removing the atmospheric delays and retrieving the interseismic deformation, even when stacking with a smaller number of interferograms. The combination of InSAR and GPS observations has given refined estimates of the deformation velocity and strain rate field with a high spatial resolution in this region. It revealed the strain variation near the fault in great detail. Such a fine strain field enabled us to gain new insights about the present-day kinematics of the 1920 rupture and the lateral expansion of the northeastern Tibet Plateau. With more and more new satellites being operated, the combination of multi-satellites data with a timespan of a few decades will allow us to investigate time-dependent slip behavior of the fault, which enable us to discriminate different deformation mechanism.

Conclusions
We used InSAR, combined with two different estimates of atmospheric path delays, to estimate the interseismic strain accumulation across the Northeastern Tibetan Plateau. The results show that the MERIS-based and ECMWF-based corrections both work well in removing the atmospheric delays and retrieving the interseismic deformation, even when stacking with a smaller number of interferograms. The combination of InSAR and GPS observations has given refined estimates of the deformation velocity and strain rate field with a high spatial resolution in this region. It revealed the strain variation near the fault in great detail. Such a fine strain field enabled us to gain new insights about the present-day kinematics of the 1920 rupture and the lateral expansion of the northeastern Tibet Plateau. With more and more new satellites being operated, the combination of multi-satellites data with a timespan of a few decades will allow us to investigate time-dependent slip behavior of the fault, which enable us to discriminate different deformation mechanism.