Interseismic and Postseismic Shallow Creep of the North Qaidam Thrust Faults Detected with a Multitemporal InSAR Analysis

Understanding the mechanisms by which earthquake cycles produce folding and accommodate shortening is essential to quantify the seismic potential of active faults and integrate aseismic slip within our understanding of the physical mechanisms of the long‐term deformation. However, measuring such small deformation signals in mountainous areas is challenging with current space‐geodesy techniques, due to the low rates of motion relative to the amplitude of the noise. Here we successfully carry out a multitemporal Interferometric Synthetic Aperture Radar analysis over the North Qaidam fold‐thrust system in NE Tibet, where eight Mw> 5.2 earthquakes occurred between 2003 and 2009. We report various cases of aseismic slip uplifting the thickened crust at short wavelengths. We provide a rare example of a steep, shallow, 13‐km‐long and 6‐km‐wide afterslip signal that coincides spatially with an anticline and that continues into 2011 in response to a Mw 6.3 event in 2003. We suggest that a buried seismic slip during the 2003 earthquake has triggered both plastic an‐elastic folding and aseismic slip on the shallow thrusts. We produce a first‐order two‐dimensional model of the postseismic surface displacements due to the 2003 earthquake and highlight a segmented slip on three fault patches that steepen approaching the surface. This study emphasizes the fundamental role of shallow aseismic slip in the long‐term and permanent deformation of thrusts and folds and the potential of Interferometric Synthetic Aperture Radar for detecting and characterizing the spatiotemporal behavior of aseismic slip over large mountainous regions.


Introduction
Recent developments in space geodetic monitoring and in fault rupture inversions have highlighted the variability in fault deformation style in term of where and when slip is accommodated on such structures. On one hand, rate-weakening segments remain locked in the interearthquake (interseismic) period and accumulate elastic strain leading to sudden rupture during earthquakes. On the other hand, rate-strengthening fault regions undergo long-term slow aseismic slip (e.g., Marone, 1998;Zhou et al., 2018). A variety of temporal behavior on aseismic patches has been observed including steady creep (e.g., Murray et al., 2001;Titus et al., 2006;Tong et al., 2013), change in creep rate after a seismic rupture (e.g., Barbot et al., 2013), postearthquake afterslip (e.g., Barbot et al., 2009;Copley, 2014;Floyd et al., 2016;Perfettini & Avouac, 2004;Zhou et al., 2018), and episodic transient creep (e.g., Jolivet et al., 2015;Rousset et al., 2016). As aseismic stress release affects the slip budget and may potentially trigger large earthquakes (Obara & Kato, 2016;Socquet et al., 2017;Vallée et al., 2013) or halt seismic ruptures (Jolivet et al., 2013), the identification of creeping fault regions and the characterization of their temporal evolution, or their relation with earthquakes, are essential for seismic hazard assessments.
Another reason for studying aseismic slip is to constrain the time period within the seismic cycle during which geomorphological structures are formed. As observed for the 2015 Gorkha earthquake in Nepal , as well as the 1998 Fandoqa and 1978 Tabas-e-Golshan earthquakes in Iran (Fielding et al., 2004), many thrust earthquakes occurring on main basal faults remain buried at depth, without causing active range fronts to uplift. The short-wavelength, permanent, Quaternary topographic uplift  (Tadono et al., 2014) and regional active faults. (a) Surface traces of the five overlapping Envisat tracks (descending tracks 90, 319, 047, and 276 and ascending track 455) used in this study with location and focal mechanisms of the eight M w > 5.2 events that occurred across the north Qaidam thrust system since 2003 (Table S1). Purple arrows show the average GPS velocities relative to Eurasia  must therefore occur at other stages of the seismic cycle on adjacent shallow faults or folds above main basal thrusts (e.g., Barnhart et al., 2018;Copley, 2014;Copley & Reynolds, 2014;Zhou et al., 2018).
The North Qaidam thrust (NQT) system is located in the arid, southernmost region of the Qilian Shan-Nan Shan thrust belt in the Northeast part of the Tibetan Plateau, running parallel to the southern margin of the Qilian Shan in a NWW-SEE direction (Fang et al., 2007;Figure 1a). A series of thrusts has been identified, which bound four tectonic units from north to south: the Zongwulong Shan, the Delingha depression, the Olonbulak Shan, and the Aimunike Shan (Fang et al., 2007;Meyer et al., 1998;Métivier et al., 1998;Figure 1b). Between the Qaidam basin and the Nan Shan, 4-6 mm/year of convergence is accommodated  across folded sediments underthrust by low angle basement faults (Burchfiel et al., 1989;Meyer et al., 1998;Tapponnier et al., 2001; Figure 1a). Fault-related foldings are developed away from the main thrusts within the Qaidam basin and can be identified in satellite imagery through light-colored weak deposits exposed in their cores Métivier et al., 1998;Walker, 2006). The region is marked by active seismicity (Figure 1b). In particular, the thrust-fold system experienced three M w 6.3 events in a 150-km along-strike area between 2003 and 2011. First, the 17 April 2003 Delingha earthquake, with a centroid depth of 16 km, ruptured part of the Zongwulong thrust, previously unrecognized as an active fault (Sun et al., 2012;Zha & Dai, 2013;Figure 1b). Following this, the 10 November 2008 and the 28 August 2009 Haixi earthquakes occurred in close proximity to each other to the West of the Delingha rupture (Elliott et al., 2011;Feng, 2015;Guihua et al., 2013;Liu et al., 2015; supporting information Table S1). These M w > 6 events were followed by five M w > 5.2 earthquakes and a period of an increased rate of seismicity in the region with a cumulative moment of 8.77 × 10 17 Nm (according to the solutions of the Global Centroid Moment Tensor Catalog (gCMT) catalog; Ekström et al., 2012). This suggests a high level of earthquake fault segmentation over a relatively short interval of time and a small area (Elliott et al., 2011).

10.1029/2019JB017692
Another striking feature is the nucleation of a cluster of M w ∼ 5 events in 2004, including two high-angle ruptures on 4 May 2004 with a high strike-slip component, and the 10 May 2004 M w 5.5 event, which shows a very similar location and mechanism to the 17 April 2003 Delingha event (Figure 1b). Whether this release of seismic strain 1 year after the 17 April 2003 M w 6.3 event belongs to its postearthquake (postseismic) phase, or implies a change of coupling of the fault system, remains unclear and requires further investigation. The 2003-2011 Qaidam sequence of earthquakes thus provides a rare opportunity to capture the strain accumulation before and after three main events within a thrust-fold system with an Interferometric Synthetic Aperture Radar (InSAR) time series approach. This allows to study the link between the release of the elastic stress along active thrust faults during earthquakes and the change of aseismic coupling along folds and faults.
In this paper, we present a Multitemporal InSAR analysis in this challenging mountainous area using Envisat C-band data archives (Figure 1). We reconstruct the relative spatiotemporal evolution of the surface displacements from 2003 to 2011, to quantify the strain distribution of the thrust-fold system at different stages of its earthquake cycle and to characterize the relationship between aseismic slip and seismic strain release during sequences of earthquakes.

InSAR Time Series Processing
Producing spatially continuous InSAR time series displacement maps over a large-scale, tectonically active region is challenging, in part due to the long-wavelength change of atmospheric conditions that can obscure the small magnitude tectonic signal (Daout et al., 2018;Grandin et al., 2012). In the case of an earthquake, another obstacle is that we need to simultaneously capture the slow interseismic and postseismic rates of movement, with their small displacement gradients masked by tropospheric phase delays, as well as instantaneous surface displacements of earthquakes (coseismic phase) that cause high displacement gradients and phase unwrapping challenges.

Differential InSAR Processing
We process the complete Envisat data archive along four ∼300-km-long and 100-km-wide overlapping descending tracks (90, 319, 047, and 276) and one ascending track (455) between 2003 and 2011. We define optimal small temporal (B t < 2 years) and perpendicular (B p < 300 m) baselines for an interferometric network of 303 interferograms that connects all 102 acquisitions with redundancy ( Figure S1). As several images are incomplete in the along-track direction, we also add long baselines interferograms to provide connections for full coverage. Differential InSAR interferograms are generated in 2 and 10 looks in the range and azimutal directions with the New Small Baselines Subset processing chain (Doin et al., 2011) based on the ROI_PAC software (Rosen et al., 2004) and using the 30-m resolution ALOS PRISM DEM (Tadono et al., 2014). We notice an increase of phase coherence when using the ALOS PRISM compared to interferograms constructed with the 90-m resolution SRTM (Farr et al., 2007) with a 5-10% decrease in the interferometric phase variance ( Figure S2).

Unwrapping and Time Series Analysis
We apply a series of corrections to the wrapped phase of the differential interferograms in order to reduce the variability of phase and thus enable unwrapping over large areas (e.g., Grandin et al., 2012; Figure 2). An example of a long-baseline wrapped interferogram of the ascending track 455 is displayed in Figure 2a. The corrections include (1) a linear term in range direction to remove residual orbital errors and effects of the satellite clock drift (Fattahi & Amelung, 2014;Zhang et al., 2014;Figure 2b), (2) linear empirical estimation of the phase delays due to differences in the stratified troposphere with a mask on deforming areas Figure 2b), and (3) the separation of the coseismic signal ( Figure 2d). All those corrections are estimated by searching on small moving subwindows the proportionality between (1) phase and range, (2) phase and elevation, or (3) phase and a coseismic template that maximizes the local phase coherence within the subwindows , respectively. The 2008 and 2009 earthquake templates (Figure 1c) are produced by stacking a pile of well-unwrapped coseismic interferograms containing the earthquakes signals that we obtained using a filter that preserves the high phase gradients (Figure 2c). We check the unwrapping reliability of this set of interferograms by computing the phase misclosure of the interferometric network. The scaled templates of coseismic displacements are then subtracted from all wrapped interferograms (Figure 2d). The remaining phase, corrected from orbital phase errors, Figure 2. Processing strategy adapted for this large mountainous area illustrated with the long-baseline preseismic, coseismic, and postseismic interferograms computed between the two acquisition dates 18 June 2004 and 08 January 2010 of the ascending track 455. Before the critical step of phase unwrapping, we apply an empirical correction of the stratified-atmospheric phase delays and the residual phase ramps (b), and we remove the coseismic signal computed from a stack of interferograms (c). After phase unwrapping (d), the previous corrections are reintroduced (e) and a final cross estimation of the stratified and ramp delays is performed and adjusted to be consistent within the interferometric network. The wrapped phase is between − and , while the unwrapped phase is shown between −8 and 8 rad. Interferograms are displayed in radar geometry and superimposed on the hill-shaded Digital Elevation Model in radar geometry. Enlargements around the earthquake area (black boxes) are also displayed at the top of each panel.
elevation-correlated phase delays, and coseismic displacements, is afterward more easily unwrapped, across both high mountain ranges and rupture areas, even for long time frame interferograms (Figure 2e).
In addition to this signal-separation procedure, the unwrapping is performed using a specific scheme (Daout et al., 2017). First, we multilook by a factor of 8 in range and 40 in azimuth and replace the amplitude of the interferograms by the phase collinearity (Pinel-Puyssegur et al., 2012) computed on the reduced phase ( Figure 2d). This allows an improvement in the apparent phase coherence of the earthquake areas, which is usually very low due to the high fringe rates (Figure 2a). We then filter with a weighted average of the interferometric phase computed on sliding windows of eight pixels and impose an unwrapping path that goes from the high to low collinear areas defined by the filter (Grandin et al., 2012). Afterward, we compute the phase difference between this smoothed unwrapped phase and the set of less-strongly filtered interferograms previously obtained from the adaptative filter step in Figure 1c. The residual high-frequency signal is assumed to be between − and and added to the smooth unwrapped phase. We then reintroduce the template as well as the previous corrections (orbital and atmospheric) previously removed to each interferogram ( Figure 1b) to reconstruct the full unwrapped phase signal (Figure 2f).
In order to capture the north-south variability of the differences of delays in the stratified troposphere, we then cross estimate, on the continuous unwrapped phase, the Line-of-Sight (LOS)-elevation relationship with the azimuthal component. We also added orbital ramps estimations, such as where LOS(x, y) are the LOS displacements for each interferogram, x, y, z are, respectively, the range direction, the azimuth direction, and the elevation, and a, b, c, d, e, and f the inverted coefficients. Inconsistencies of parameters a, b, c, d, e, and f are detected by a least-square inversion of all independent coefficients before applying the correction ( Figure 2f). The time series analysis is then performed with the New Small Baselines Subset method López-Quiroz et al., 2009). Inconsistencies in the interferometric network are mapped to detect unwrapping errors, and high-priority unwrapping paths are imposed, if necessary, to correct them in new iterations of unwrapping and time series analysis.

Time Series Results
The In addition, track 276 images capture only partially the feature in the far range of the track with therefore a larger radar incidence angle than track 047. The 10 May 2004 event is also unfortunately only captured by the descending track 047 due to the absence of acquisitions before 10 May 2004 for tracks 455 and 276. In the following, we focus on the extraction of these various phases of the seismic cycle with a successive parametric time series decomposition approach.

Signal Separation With a Parametric Approach
In track 047, the surface displacement time series relative to the acquisition date 25 May 2004 of points very close to the gCMT centroid of the 10 May 2004 M w 5.5 earthquake show an abrupt and clearly coseismic ∼2 cm step between the 16 April and the 25 June acquisitions (Figures 4a and 4b).
South of the epicentral area of the 17 April 2003 Delingha earthquake, we observe an increase of surface displacements that are localized at the southern side of the anticline of the north Qaidam range. The presence of such steep and narrow surface displacement gradients at the surface leads us to associate the deformation to transient aseismic afterslip (Barbot et al., 2009) and discard any viscoelastic relaxation in the shallow crust overlying a blind thrust (Lyzenga et al., 2000). From 2004 to 2011, no aftershocks with M w > 4.5 occurred in this area, while surface displacements gradually increase to finally reach about 6 cm toward the satellite in both ascending and descending geometries (Figures 3 and 4). A common approximation for afterslip following a stress change (Perfettini & Avouac, 2007) is a logarithmic function of the surface displacements, u(t), equal to log(1 + t−t 0 ), with the characteristic time and t 0 the time of the stress-changing event based on the rate-and-state-friction parametrization (Marone, 1998). This approximation provides a reasonably good fit with a characteristic time, 2003 , of 4 months (red lines in Figure 4b). In order to separate the coseismic signal due to the 10 May 2004 earthquake and the postseismic signal from the 17 April 2003 Delingha earthquake and to minimize residual atmospheric and topographic noise, we estimate from descending track 047 the best-fit offset and logarithmic postseismic surface displacements for all pixels (Figure 4c). Each cumulative LOS displacement for each pixel can thus be described as where t is the SAR acquisition time, A 2004 and B 2003 correspond, respectively, to amplitude of the step and the logarithmic function, and 2003 is the characteristic time. We solve the inverse problem for each pixel independently with a singular-value decomposition approach neglecting small eigenvectors inferior to 0.1 and associating an error bar to each date equal to the phase misclosure from the time series analysis (blue vertical bars in Figure 4b). In the same way, for ascending track 455 and descending track 276, with no acquisitions before the 10 May 2004, we map the lateral variability of the afterslip by estimating the best-fit logarithmic postseismic surface displacements for all pixels with a characteristic time of 4 months ( Figure 5), such as A small-amplitude displacement signal that amounts to ∼2.0 cm is observed for the 10 May 2004 event, consistent with slip remaining buried at depth (Figure 4c, left). Postseismic slip is, on the contrary, measured south of this 2004 lobe of surface displacement, with mainly a steeper signal along a 13-km-long and 6-km-wide surface patch (Figure 4c, right, and Figure 5). Despite the poor sampling of tracks 276 and 455, the observed signal is consistent within the three geometries (Figure 5a), which suggests a mainly vertical uplift. The three surface displacements highlight a sharp, high increase of the uplift signal at the southern side of the Zongwulong Shan that coincides with the topographic front of the Delingha anticline (Fang et al., 2007;. To the North, the exposed geology of the Zongwulong Shan, which rises 4.5 km high, is predominantly composed of Paleozoic units, while south of this mountain region, the anticline cuts across Pleistocene alluvial and lacustrine deposits (Figure 5b).

Two-Dimensional Fault Model
To get insights about the fault-system geometry along the main crustal shortening direction, we use a two-dimensional Bayesian inversion tool developed by . We model one or more connected edge dislocations embedded in an elastic half-space (Savage, 1998;Segall, 2010) with variable dips, depths, down-dip widths, and fault slip. The synthetic surface displacements are compared to the observed LOS displacements of postseismic afterslip along a 70-km-long and 10-km-wide profile centered on and perpendicular to the Delingha anticline (as defined in Figure 5a but 20 km longer). For the northernmost segment, we invert for the depth, LD 1 , the horizontal position in reference to the center of the profile, D 1 , the width of the segment, W 1 , and a dip angle, 1 (Figure 6a). For the other dislocations, i, we invert for their horizontal distance, D i , and vertical distance, H i , to the previous dislocation ( Figure 6a). We also add an unknown azimuthal linear trend to tie the InSAR data to the model. The degree of freedom for multisegment models is also reduced by forcing a uniform horizontal slip component, called shortening, across all segments, such that mainly the segment dips are controlling the synthetic displacement. To do so, we only invert for the shortening, S h , and project this shortening on a dip-slip component, DS, on each segment, i, We find that a single dislocation underestimates the amount of displacements close to the Delingha anticline ridge and overestimates the amount of displacement, 5 to 10 km north of it, because the synthetic displacement spreads with a long wavelength ( Figure S8). A two-segments model improves the fit to the observations but still does not capture well all its complexity. The steep gradient close to the surface trace of the fault linked to the steepening of the fault as it approaches the surface is not well predicted, and the long wavelength signal 7 to 25 km away from the fault trace is not well matched for both tracks ( Figure S8).
The most convincing results were obtained using three segments ( Figure 6a and Table S2). The model shows that the LOS data are in agreement with a ramp-flat-ramp geometry that extends over a range of ∼20 km.   Table S2.
From this model, as well as the active faults mapping and work of Yin, Dang, Wang, et al. (2008) and Fang et al. (2007), we therefore propose in Figure 5c the structural schematic cross section of the region. We infer that the 17 April 2003 Delingha earthquake has occurred on the ∼30 • north-dipping fault plane at ∼15 km depth (Sun et al., 2012; Table S1), which corresponds to the modeled down-dip extent of the northernmost deeper slipping segment that roots below the Zongwulong high ranges. We also observe a long wavelength of surface displacements similar to the 2003 aseismic signal for the 10 May 2004 earthquake, ∼10 km north of the Delingha topographic front (orange line in Figure 5c). According to the model, this more distributed signal over a length of ∼15 km is linked to slip on the northernmost deeper segment. As the focal mechanism of the May 2004 earthquake is similar to that of the 2003 earthquake (Table S1), we also suggest that it may have occurred on the ∼30 • north-dipping deeper segment or on a ∼30 • north-dipping splay segment.

LOS Velocity
Unfortunately, due to the poor sampling of track 455, the interseimic linear velocity from the ascending geometry presents more noise than signal. As for the postseismic following the 2003 event, we explore several characteristic times for the postseismic surface motion of the 2008 and 2009 ruptures, 2008 and 2009 , and choose a characteristic time of half a month. We also impose that displacements from the logarithmic functions to be in the same sense of direction as their associated earthquake steps and solve the constrained inverse problem with a sequential least-square algorithm ( Figure S9). Resulting coseismic and postseismic ground motion maps ( Figure 7b) show a very smooth signal, which validates this pixel-by-pixel approach. We also compare the surface displacements for the four events along a 50-km-long and 10-km-wide profile in Figure 7a, where the lines correspond to a sliding median over the LOS displacements.
A smooth displacement signal reaching ∼6.8 cm and a shorter wavelength displacement signal with a higher maximum of ∼35 cm toward the satellite and a minimum of ∼8 cm away from the satellite are observed A broad, low-amplitude surface displacement signal (∼2 and ∼1.5 cm for descending and ascending geometries, respectively) is then measured during the postseismic phase of the 2008 earthquake. The maximum of this signal is further north than the maximum of the coseismic surface displacements (Figure 7b), indicating a down-dip afterslip for the 2008 north-dipping scenario (Guihua et al., 2013) or an up-dip afterslip for the conjugate south-dipping scenario (Elliott et al., 2011). The 2009 earthquake caused stronger postseismic surface displacements (∼4.5 and ∼2.7 cm for descending and ascending geometries, respectively) than the 2008 earthquake. The pattern of the 2009 afterslip surface displacement is similar to the 2009 coseismic surface displacement pattern but extends also slightly further north (Figure 7b).

Steady-State Aseismic Slip
The 8 Figure S9). Due to the near-polar orbit of Envisat with an almost north-south flight direction, the LOS displacements are almost insensitive to the overall north-south 4-6 mm/year of shortening. They therefore reflect mostly the local vertical movements (Figure 9). In mountainous areas, the aliasing of oddly sampled seasonal time-varying stratified atmospheric signal pollutes the decadal linear trend of the time series Daout et al., 2018). To highlight the benefits of previously applied atmospheric corrections and to show that they did not remove any major gradient of displacements, we compare our results (in Figure 8) with a velocity map obtained without empirical signal-elevation corrections ( Figure S11 obtained without applying equation (1) on the unwrapped interferograms). The resulting map shows similar results but with stronger velocity gradients across mountain ranges for time series without atmospheric correction, which result in discrepancies in the overlapping areas of the four tracks.
At first order, the relative velocity field, referred with respect to points within the Qaidam basins (pixels with no relative surface displacements in light colors or white in Figure 8), shows displacements away (blue tones in Figure 8) from the satellite from the Qilian Shan to the southern side of the NQT and displacements toward the satellite (red tones in Figure 8) within the Qaidam basin, which might be related to the overall north-south shortening. However, no large major gradient of displacements between the Kunlun and Qilian Shan is visible, and the velocity map shows both ground velocity away and toward the satellite. Steeper velocity gradients are observed within Qaidam Basin and illuminate slow aseismic slip at many segmented thrust faults a few tens of kilometers long of more than ±1 mm/year. More particularly, several thrusts, in the Qaidam basin, indicate a movement toward the satellite, at longer wavelength than the short wavelength of elevation changes, likely linked to the shallow steady-state slip at high-angle thrusts running beneath young anticlines and uplifting them (Figure 9a). For example, such uplift is visible, from west to east, south of the NQT, across the Luliang thrusts, the Xitieshan thrusts (XT), the Aimunik thrust, or the OLT (red arrows Figure 8a). The part of NQT that broke during the 2008-2009 earthquake sequence does not show major velocity toward the satellite. More particularly, no creep is observed across the shallow north-dipping portion of the OLT fault system that did not break during the 2008-2009 sequence.
Small localized surface displacements (red tones) also appear at several fault jogs. For example, east of the 2003 postseismic slip, the Delingha anticline shows a geometrical discontinuity along strike (Figure 9b). The uplift signal (red tones within the jog) is of a similar pattern than the one observed for the 2003 afterslip across the Delingha anticline with a stronger signal at the fault jog and a broader and lower amplitude signal further north. It might be interpreted as volumetric deformations or by an aseismic slip in this geometrical irregular area of the thrust system.
In general, nontectonic deformations, such as anthropogenic activities, landslides, and sediment compaction, are abundant on our velocity map and are superimposed to slower-moving tectonic displacements (Figure 8). implied by the geomorphology would suggest that the hanging wall is going down. Instead, the XT fault may here be acting as a barrier to groundwater flow and sedimentation, and the observed sharp signal may be due to a long-term increase of water content or sedimentation infill processes in this endorheic basin.
Over the Xitieshan mine in the central part of our map, a 4-km-long subsiding area can be seen with a clear linear trend in the decadal time series of surface displacement (Figure 9a). South of the XT, the apparent large velocity gradients are associated with the intense anthropogenic surface-changing activities such as crop farming with reworking of fields and groundwater flows. Furthermore, several small-scale landslide on low-angle slopes can be observed in the InSAR time series. For example, at an inter-segment of the XT, in the west of our study area, a clear scarp in soil is visible in the optical images in a small area with relatively large velocity gradients of up to a few centimeters per year (Figure 10d). The velocity map shows areas of decorrelation in the closest proximity of the scarp, likely indicating the fastest displacement rates. These are followed downslope by areas characterized by displacements away from the satellite and toward the satellite, respectively, indicating the zone of depletion and accumulation of the landslide. This small landslide at a XT jog may furthermore attest the activity of the fault as it may have been triggered by shallow slip on the fault.  Figure 10e shows alluvial fans observed in the north of the study area. The mega fan in the center of the image covers smaller and older fans. Textural changes between the older fans and the more recent fan are also observed in the optical images, these possibly indicating different hydrological regimes and water content in the sediments. On the surface of the most recent mega fan, relatively small signs of activity are detected with InSAR only in association with active debris flow channels, where ongoing erosional processes are related to displacements away from the satellite. Elsewhere on the surface of the mega fan, no other clear surface displacements are detected. However, on all other older fans, we observe significant and widespread signals away from the satellite with a clear linear trends ranging from −3 to −5 mm/year (time series of Figure 10e). As these frequent features indicate no sign of soil slip in the optical images, possible explanations in the long-term displacement rate is the combined effects of slow sediment compaction and water volume changes within older fans.
There are more regions with widespread motion away from the satellite that are very difficult to interpret, such as signals across a major back-thrust structure within the northern boundary of the NQT at latitude 38 • 75 (Figure 10f). The negative velocity gradient across alluvial fans north of the back-thrust might again be related to sediment compaction. However, the abrupt change of the displacement signal further north also corresponds to sharp changes in the sediment colors in the optical images and could instead be controlled by a fault. Though no clear fault scarp is visible, a tectonic structure could be covered by loose sediments that get reworked quickly by erosional processes. Assuming a tectonic origin of the deformation, the sharp velocity gradient could be associated to a left-lateral fault motion north of the back-thrust or to local subsidence linked to a flat portion of the thrust system.

Discussion
Our time series measurements show the potential for InSAR to detect and characterize the spatial and temporal behavior of slow-slip faulting and its relation to seismic activation.
Postseismic surface displacement observations bring first insights into the geometry of the NQT. There have been several studies that aim to constrain the geometry of the NQT fault system and have been undertaken for the three M w 6.3 events (Elliott et al., 2011;Guihua et al., 2013;Feng, 2015;Liu et al., 2015;Sun et al., 2012;Zha & Dai, 2013), but their conclusions differ on the geometry of the 2003 and 2008 deeper earthquakes, with the dip direction being ambiguous. Structural cross sections of the NQT also have several different interpretations in the literature, with both south-dipping and north-dipping structures suggested to offset large amounts of sediments in the cross section proposed by Fang et al. (2007) and Yin, Dang, Wang, et al. (2008). Sun et al. (2012) inferred a 31 • north-dipping plane for 2003 Delingha earthquake, from the inversion of seismological waveforms from the mainshock and relocation of aftershocks recorded with a local network. On the contrary, the inversion of global body waves and surface waves (Zha & Dai, 2013) suggests a best-fitting high-angle south-dipping solution. In this study, surface displacements associated with the postseismic slip of the 2003 earthquake leave little ambiguity on the activity and the dip direction of this shallow part of the ZT and clearly illuminate a north-dipping geometry, with slip reaching the surface ( Figure 5). For the 2008-2009 earthquake sequence, Elliott et al. (2011) andFeng (2015) constrained nearly coplanar high-angle south-dipping planes using surface displacement maps based on InSAR measurements for the 2008 and the 2009 Haixi M w 6.3 earthquakes. However, the analysis of a GPS time series, field data, and aftershock relocations presented by Guihua et al. (2013) suggests that the 2008 earthquake might have occurred on a north-dipping low-angle fault with left-lateral oblique thrust motion. As we observe in this study surface displacements for the 2008 and 2009 postseismic that roughly coincide in space, it may suggest that both earthquakes have triggered the same velocity-strengthening asperities. This supports a model in which the 2008 and the 2009 events rupture an identical south-dipping structure (Elliott et al., 2011). However, addressing this ambiguity robustly would require a three-dimensional modeling of both seismological and InSAR time series data, which is beyond the scope of this study.
We report a large post-earthquake slip variability, with ∼6 cm of uplift 3 years after the 2003 earthquake ( Figure 4) and, on the contrary, small and large, but short-duration, afterslip following the 2008 and 2009 earthquakes, respectively (Figure 7). Shallow creep is in agreement with a rate-strengthening fault behavior of the uppermost part of the crust. However, a notable result of this study is the duration and the spatial extent of the observed surface displacement following the 2003 earthquake. Afterslip larger than coseismic displacements has already been observed for other earthquakes, such as the 2004 M w 6.0 Parkfield earthquake in California. However, in this case, the main aseismic release occurred within several months after the main shock (e.g., Barbot et al., 2009;Johanson et al., 2006), which highlights that large magnitude and long duration afterslip can be triggered by moderate earthquakes. This raises the question of the validity of existing velocity-strengthening mechanical models, based on static stress changes and faults frictional properties only and the possibility of triggering a large amount of aseismic slip on wide fault patches due to fault geometrical controls (Ong et al., 2019;Romanet et al., 2018). In addition, spatial agreement of the 2003 afterslip surface displacements with the Delingha anticline suggests that permanent plastic deformation along the folds has been activated by the release of elastic stress during the 2003 Delingha earthquake on the basal seismic fault. This plastic deformation may play a role in the long duration of the observed postseismic surface displacements.
Another notable observation from the time series of displacements of track 047 is the step in the surface displacements after the 10 May 2004 M w 5.5 earthquake. A change in the surface displacements appears to start after the acquisition on the 25 May on track 047 (acquisition date marked with a black arrow on time series, Figure 4b), several days after the 10 May 2004 M w 5.5 earthquake. This step is consistent in all time series points in the near field to the creeping fault. To further highlight this point, we perform an Independent Component Analysis with four components for descending track 047 ( Figure S13). The first two components, which explain the majority of the signal ( Figure S13c), are clearly associated with the logarithmic postseismic slip after the 2003 earthquake and the May 2004 earthquake offset. The third and fourth components explain small variances of the total signal. We, however, interpret the fourth component as residual tropospheric delays due to its spatial correlation with topography and a seasonal pattern in time.
As the third component shows a change in surface displacements after the acquisition on the 25 May, it is potentially associated with a transient slow-slip behavior beneath the Delingha anticline as previously described on the cumulative displacement time series. It may be attributed to the small stress redistribution following the M w 5.5 earthquake and a north-south migration of the stress transfer. Such behavior could point to the sensitivity of aseismic slip to small perturbations of stresses (Obara & Kato, 2016). Unfortunately, it is here difficult to conclude due to the poor temporal sampling, the small surface displacements, and the high level of noise.
Two dimensional modeling solutions (Savage, 1998;Segall, 2010), such as that applied here, are useful due to their simplicity. They allow to improve our understanding of the deformation by simplifying the slip representation and by reducing the number of free parameters. In this study, the solution applied allows to explore the influence of the geometry of several connected thrust segments on the surface displacements. The model shows the compatibility of the surface displacements with a ∼20-km-wide ramp-flat-ramp structure. The LOS displacements highlight a steep signal that coincides with the Delingha anticline over 13 km along strike. Similarly to other regions such as the Kerman province in Iran (Walker, 2006), surface displacements are explained with a thrust that steepens as it approaches the surface resulting in the folding of sediments on a comparable wavelength. Further north, at ∼10 km distance to the steep displacement gradients, another ramp structure branching from the flat-ramp thrust root is also required to explain the broader signal. The Delingha thrust-fold system is thus stepping out into the basin away from the main Zongwulong ranges as it is also observed in other tectonically comparable regions (Chen et al., 2007;Walker, 2006). Moreover, the modeled shortening of ∼9 cm for 3 years of cumulative displacements implies a horizontal slip rate of 3 cm/year, more than four times larger than the shortening rate of 4-6 mm/year over the whole NQT, from the Kunlun to the Qilian Nan ranges derived from the GPS velocities ; Figure  S12). We can therefore conclude with confidence that the aseismic slip measurements across the Delingha anticline are not steady-state processes. Clearly, this two-dimensional representation could be improved by using more complex Earth models compared to a uniform elastic half-space and/or a three-dimensional solution. A two-dimensional approximation introduces a bias in the dip-slip component because it neglects the finite length of the fault creeping patches. In addition, although a fault plane must evolve beneath the Delingha anticline, part of the signal might also be associated with plastic deformation of the folds in the hanging wall of the reverse fault.
Strike-slip earthquakes have been observed to nucleate or terminate near along-strike geometrical irregularities along faults (azimuth changes, extensional or compressive jogs, and branching of secondary faults), with fault inter-segments usually associated with local slip minima (e.g., Klinger, 2010;Lasserre et al., 2005;Manighetti et al., 2009Manighetti et al., , 2015Wesnousky, 2006). While along-strike segmentation of strike-slip faults is discussed considerably in the literature, less is known about the down-dip segmentation of reverse faults in intracontinental settings (e.g., Copley & Reynolds, 2014;Elliott et al., 2011;Mackenzie et al., 2016). The moderate sizes of the earthquakes in this area might suggest segmentation control on the length of the ruptures. The proposed two-dimensional model shows a change of dip angle with depth that might have played a role of barrier for the 2003 rupture propagation. Finally, as we observe significant differences in the amount and extent of afterslip for the 2003 and 2008 events, even if their magnitudes, depths, and mechanisms are very similar, this supports the view that transient aseismic slip behavior primarily depends on the size and geometry of velocity-strengthening areas, rather than the amount of stress changes from the main shock (Wimpenny et al., 2017;Zhao et al., 2017).
The simple parametric decomposition into functional basis functions does not allow the quantification of heterogeneities in the temporal behavior of the afterslip, which may be activated at different times due to changes in stresses associated with the main shock (Marone, 1998;Perfettini & Avouac, 2007). In addition, due to the limited number of acquisitions after the 2008 and 2009 earthquakes, the relaxation time scales of the two afterslip periods are not robustly estimated. However, this parametric decomposition offers a simple linear fit to the time series data, despite the poor temporal sampling of the observations, in order to map the spatial distribution and amplitude of afterslip for the three M w 6.3 reverse events. It also provides a simple way to reduce undesired signal in the time series, therefore allowing the extraction of the average LOS velocity displacements from 2003 to 2011 over the entire region.
The ∼2,700-m-high Qaidam basin is the largest basin inside the Tibetan Plateau, bounded by the 1,600-km-long NE-striking left-lateral Altyn Tagh Fault to the northwest, the Qilian Shan-Nan Shan fold-thrust belt to the northeast, and the Kunlun transpressional fold-thrust system to the south. This once-giant basin, mainly composed of Cenoizoic sands, gravels, and thick gypsum deposits (Fang et al., 2007;Métivier et al., 1998;, results from the exhumation of the Kunlun and Qilian Shan bounding ranges, which have slowly closed the basin from the surface hydrological flow and atmospheric circulation (Métivier et al., 1998;Tapponnier et al., 2001). In the Qaidam basin, the present deformation is the opposite to what is reported along the major structures of the Tibetan Plateau such as the Altyn Tagh, the Kunlun, or the Haiyuan Fault transpressional system, north and south of the Qaidam basin, where strain is localized along lithospheric weaknesses that dominate the regional deformation (Daout et al., , 2018Meyer et al., 1998;Tapponnier et al., 2001). Instead, the LOS velocity across the NQT (Figure 8) presents multiple segmented thrust faults, few tens of kilometers long and 10-40 km wide, which slip at a rate of a few millimeters per year or less. The slow slip rates of these segmented thrusts suggest that these faults are not major structures that extend into the lithosphere with the geometry and rheology they have in the upper crust but that some deformation is transferred from the main ranges via décollements or shear zones beneath them (Figure 5c). The aseismic behavior of those segmented thrusts and folds, likely linked with gypsum deposits, also demonstrates the importance of rate-strengthening areas for slip-on shallow thrust faults (Copley & Reynolds, 2014;Mackenzie et al., 2016). Spatial agreement of steady-state uplift, or transient post-earthquake slip (as seen here with the 2003 or 2009 afterslip), with narrow geomorphological features, suggests that the Quaternary landscape development in the Qaidam basin is governed in part by the aseismic phase of the seismic cycle. In addition, the large amount of sedimentary and basin infill processes observed in the LOS velocity map (Figure 10) attests to the quantity of sedimentation involved in the upstream catchment filling. This forms an intimate link with tectonic processes that smooth out the relief into a broad plateau of uniform height Métivier et al., 1998;Tapponnier et al., 2001). The ongoing building of the Tibetan Plateau seems therefore to result from the combination of basin infill processes and vertical partitioning of the deep-seated convergence along shallow and high-angle thrusts creating the short-wavelength permanent topography through aseismic processes.

Conclusion
In this paper, we present a multitemporal InSAR processing method to assist the unwrapping of interferograms affected by high-fringe rates, as well as an approach to decompose time series of surface displacement into functional basis functions. Aided by these methods, we successfully derive high-Signal to Noise Ratio, large-scale, and continuous time series maps of surface displacements in the high mountain ranges of the northeastern part of the Tibetan Plateau. We image the surface displacement of the three main earthquakes that occurred in this 150-km-long area in May 2004, October 2008, and August 2009, as well as the slow transient and steady-state aseismic slip before and after those earthquakes. We construct a continuous 300 km × 300 km relative LOS velocity map in this challenging mountainous area. We report a large