Post‐Earthquake Fold Growth Imaged in the Qaidam Basin, China, With Interferometric Synthetic Aperture Radar

Questions regarding the development of folds and their interactions with the seismic faults within thrust systems remain unanswered. However, estimating fault slip and earthquake hazards using surface observations and kinematic models of folding requires an understanding of how the shortening is accommodated during the different phases of the earthquake cycle. Here, we construct 16‐years of interferometric synthetic aperture radar time series across the North Qaidam thrust system (NE Tibet), where three Mw 6.3 earthquakes occurred along basement faults underlying shortened folded sediments. The analysis reveals spatio‐temporal changes of post‐earthquake surface displacement rates and patterns, which continue more than 10 years after the seismic events. The decomposition of the Sentinel‐1 ascending and descending line of sight velocities into vertical and shortening post‐earthquake components indicates that long‐term transient uplift and shortening is in agreement with the deformation that might be expected from kinematic models of folding. Long‐term uplift coincides spatially with young anticlines observed in the geomorphology, with steep gradients in the forelimbs, gentle gradients in the back‐limbs, an absence of subsidence in the footwalls, and higher gradients along interpreted bedding planes. Long‐term shortening is also different from the surface displacements expected for typical time‐varying creep on a narrow dislocation interface and shows rates higher than the average convergence across the whole region. These findings provide evidence for anelastic fold buckling during the post‐earthquake phase and highlight the contribution of distributed aseismic deformation to the growth of topography.

Some studies have used geodetic measurements of coseismic surface displacements to document growth of anticlines (e.g., Belabbès et al., 2009;Nissen et al., 2007;Stein & King, 1984;Tizzani et al., 2013). However, recent satellite-based measurements have highlighted the observation that the surface ground deformation associated with thrust earthquakes often does not match the present-day or Quaternary active geomorphology in extent and that it is unlikely that repeated earthquakes of the same type could reproduce the observed topography (e.g., Ainscoe et al., 2017;Barnhart, Brengman, et al., 2018;Copley, 2014;Elliott, Jolivet, et al., 2016;Mackenzie et al., 2016). In addition, although simplified models are valuable in many ways, modeling shortening and the associated uplift across a fold-and-thrust belt system with one or two elastic dislocations is an arguable simplification of the complicated history and geometry of a thrusting morphology, which may involve several interconnected faults and folds that interact with each other, multiple bends, and anelastic deformation as in the upper plate around the bends (e.g., Daout, Barbot, Peltzer, Doin, et al., 2016;Davis et al., 1983;Medwedeff & Suppe, 1997;Sathiakumar et al., 2020;Tapponnier, Meyer, et al., 1990;Whipple et al., 2016). Some other geodetic observations have recently documented or/ and modeled episodic aseismic ground deformations (e.g., Mariniere et al., 2020) or long-periods of after-slip that correlate with the present-day geomorphology (e.g., Barnhart, Lohman, & Mellors, 2013;Copley, 2014;Daout, Sudhaus, Kausch, et al., 2019;Elliott, Bergman, et al., 2015;Fielding et al., 2004;Mackenzie et al., 2016;Wimpenny et al., 2017;Zhou et al., 2018). These measurements suggest that the permanent deformation in fold-and-thrust belts might be sometimes created by distributed off-fault deformation (i.e., anelastic buckling of the medium), or aseismic slip on secondary faults branching from the main earthquake fault or around it, and which occur during stages of the earthquake cycle other than the seismic event. Despite the likely importance of those phenomena, little is known about the link between seismic slip on the main fault and slip on secondary faults or distributed deformation following a seismic event. To understand how topography is generated by repeated earthquake cycles and/or by other processes, and to constrain the period within the earthquake cycle during which geomorphological structures are formed, it is necessary, therefore, to bridge the discrepancy between simple geodetic models of fault slip and geological models of fold-and-thrust belts, by combining measurements of surface displacements over a long-time period and with realistic structural fold-and-fault geometries.
The North Qaidam thrust (NQT) system is located at the southernmost part of the Qilian Shan-Nan Shan thrust belt and exhibits segmented anticlines striking in a northwest-southeast direction ( Figure 1). In this NE part of the Tibetan Plateau, while major strike-slip faults, such as the Kunlun and Haiyuan Faults, have received considerable attention from the tectonic and geodetic communities, data from the thrust systems in the South Qilian Shan are not widely reported, and little is known about the structural geology of the region (e.g., Fang et al., 2007;Guihua et al., 2013;Meyer et al., 1998;Tapponnier, Zhiqin, et al., 2001;Yin, Dang, Wang, et al., 2008;Yin, Dang, Zhang, et al., 2008). According to sparse global navigation satellite system (GNSS) measurements, about 4-6 mm/yr of convergence is accommodated there in an overall N22°E direction (Liang et al., 2013;Wang et al., 2017) ( Figure S1). The region is marked by active seismicity with three M w 6.3 events occurring in a 150 km along-strike area between 2003 and 2009 ( Figure 1). The 17 April 2003 Delingha earthquake ruptured part of the Zongwulong thrust with a centroid depth, derived from local seismic network, of 16 km (Sun et al., 2012), in a section of the North Qaidam Shan which was previously unrecognized as active (Figure 1b). It was followed by extensive and long-lasting logarithmic post-seismic slip that coincided with the Delingha anticline and continued into 2011 along a north-dipping creeping structure (Daout, Sudhaus, Kausch, et al., 2019). As reported by the authors (Daout, Sudhaus, Kausch, et al., 2019), shallow creep is in agreement with a rate-strengthening fault behavior of the uppermost part of the crust, but the duration, the amplitude, and the spatial extent of the observed surface displacements following the 2003 earthquake, during the 2003-2011 Envisat period, is a notable particularity, despite the low strain rate of the area. The 10 November 2008 and the 28 August 2009 Haixi earthquakes occurred in close proximity to each other within the Olongbulak Shan, to the west of the Delingha rupture (Figure 1), and were followed by a period of an increased rate of seismicity in the region and aseismic slip (Daout, Steinberg, Isken, et al., 2020;Daout, Sudhaus, Kausch, et al., 2019;Elliott, Parsons, et al., 2011;Feng, 2015;Guihua et al., 2013;Liu, Xu, Wen, & Fok, 2015;. Daout, Steinberg, Isken, et al. (2020) modeled the 2008 earthquake with a ∼32° north-dipping fault that roots under the Olongbulak pop-up structure at ∼12 km depth (top edge of the fault) (Figure 1b), as well as down-dip aferslip along a coplanar north-dipping plane. They also inferred three south-dipping segmented ∼55-75° high-angle faults for the 2009 earthquake with fault tip depths between ∼2.5 and 4.5 km and three post-seismic faults with similar geometries than the three earthquake patches (Daout, Steinberg, Isken, et al., 2020). They, therefore, inferred a low-angle north-dipping décollement that roots at the bottom of the fold-andthrust belt structures of the NQT transferring the shortening rate from the north, along the Qilian Shan, to the south, in the Qaidam basin, in agreement with the southward expansion of the South Qilian Shan since the Neogene (Pang et al., 2019) (Figure 1b).
In this study, we conduct a Multi-Temporal Interferometric Synthetic Aperture Radar analysis across this northeastern part of the Tibetan Plateau using both C-band Envisat (2003C-band Envisat ( -2011 data derived from Daout, Sudhaus, Kausch, et al. (2019), and C-band Sentinel-1 (2014-2019) data (Figures 1 and S2). We processed the Sentinel-1 data set and corrected both Envisat and Sentinel-1 data sets for tropospheric path delays using the newly ERA-5 atmospheric model provided by the European Center for Medium-Range Weather Forecasts (ECMWF). To capture and quantify the spatio-temporal change of the aseismic strain release that follow the 2003 Delingha ( Figure 2) and the 2008-2009 Haixi (Figure 3) earthquakes, we extract linear ground velocities over different time spans with a parametric decomposition. We then decompose both ascending and descending line of sight (LOS) velocities into an N22°E horizontal velocity, which is the main direction of shortening, perpendicular to the overall orientation of the North Qaidam thrust systems, and a vertical surface velocity (Figures 4 and 5), to infer two structural models of the folds and associated faults for the two areas. Top right inset: regional setting and location of the study area (red box). (b) Two topography profiles marked AA' and BB' in (a) with regional conceptual interpreted fault structures and earthquake locations derived from Daout, Sudhaus, Kausch, et al. (2019), Daout, Steinberg, Isken, et al. (2020), and the available cross-sections of Yin, Dang, Wang, et al. (2008), Fang et al. (2007), and Guihua et al. (2013).   Figure S3). Additionally, we processed three Sentinel-1 tracks (descending tracks 004 and ascending tracks 172 and 099) acquired in interferometric wide-swath mode (width of ∼250 km) from 2014 to 2019 ( Figure S2). The processing is carried out using the New Small Baselines Subset (NSBAS) processing chain (Doin, Lodge, et al., 2011, Doin, Twardzik, et al., 2015, which includes routines from the ROI_PAC software (Rosen et al., 2004) and precise azimuthal coregistration with the spectral diversity DAOUT ET AL.
10.1029/2020JB021241 5 of 15  technique for Sentinel-1 data (Grandin, 2015;Grandin, Klein, et al., 2016;Métois et al., 2020). In order to avoid possible biases linked to the inclusion of short temporal baselines (Ansari et al., 2020), we include interferograms with both short and long temporal baselines in the analysis ( Figure S3) that we successfully unwrap by applying empirical atmospheric corrections before unwrapping and by imposing an unwrapping DAOUT ET AL.
After unwrapping, we re-introduce the empirical corrections that had been applied on the wrapped phase and construct cumulative time series on a pixel-by-pixel basis with the NSBAS method (Doin, Twardzik, et al., 2015;López-Quiroz et al., 2009). Residual unwrapping errors are first detected by visual inspection on interferograms or with high phase misclosure computed from the time series analysis (Doin, Lodge, et al., 2011;López-Quiroz et al., 2009). Errors are then corrected by imposing high-priority paths during a new unwrapping iteration. As part of the time series analysis, remaining unwrapping errors are further automatically identified and resolved using an iterative technique (López-Quiroz et al., 2009). In order to improve the referencing estimation and exclude surface periglacial processes from this analysis, we then, mask all pixels showing in the timeline of ground deformation, seasonal displacements higher than 3.5 mm associated with frost heave and thaw settlement of the permafrost active layer (Daout, Dini, Haeberli, et al., 2020). Cumulative LOS displacements (LOS(t)) from both Envisat and Sentinel-1 data are, subsequently corrected for tropospheric effects using the ERA-5 atmospheric model provided by ECMWF (Doin, Lasserre, et al., 2009;Jolivet et al., 2011), and referred to an 100 × 100 pixels area within the Qaidam basin of low elevation, low deformation gradient and high coherence. Afterward, in order to map and quantify the spatio-temporal strain evolution and associated uncertainties, we proceed to a simple linear parametric decomposition of the co-seismic and post-seismic ground deformations for the two earthquake areas.

Parametric Decomposition
Sentinel-1 ascending and descending tracks 004, 172, and 099 are first decomposed, assuming zero displacement during the first acquisition, into a single linear trend from 2014 to 2019 (V 2014−2019 in mm/yr) in both rupture areas ( Figure S4): where t is the SAR acquisition time. In order to quantify the spatio-temporal where t 09 , the timing of the 2009 earthquake; H, the Heaviside function; and A 09 and V 09−11 correspond to the amplitude of the August 2009 coseismic surface displacements (in mm) and the mean velocity from 2009 to 2011 (in mm/yr), respectively (Figure 3). The 10 November 2008 earthquake occurred south of the Olongbulak Shan with a depth of ∼12 km and produced short-term and low-amplitude down-dip afterslip (Daout, Steinberg, Isken, et al., 2020). Conversely, the 28 August 2009 shallow event occurred at a depth of ∼5 km and produced afterslip around the co-seismic area (Daout, Steinberg, Isken, et al., 2020 it is important to note that the short-term post-seismic observed in the Envisat data during the 2009-2011 period, might be a combination of the post-seismic surface displacements from both Haxi earthquakes.

Uncertainty Estimation
Two ascending and one descending Sentinel-1 track, with a total of 248 acquisitions, illuminate at a regional scale the LOS surface velocity from 2014 to 2019 over a 120,000 km 2 area of northeastern Tibet ( Figure S4). Uncertainty maps, derived from the time series decomposition ( Figure S5), are of 0.6 [0.4-0.8] mm/yr on average (with 95% confidence intervals), with a maximum uncertainty of up to ∼2 mm/yr in Qilian and Kunlun high ranges or at the bottom of some localized valleys. Those average low uncertainties, as well as the good agreement of the data within the profile in the overlapping area of the two ascending track 099 and 172 validates the quality of the results ( Figure S6a). Histograms of differences of the overlap are approximatively Gaussian with mean close to zero and standard deviation of less than ∼0.5 mm/yr ( Figure S6b). Descending track 004 also shows a good agreement with the regional Envisat velocity map from 2003 to 2011 obtained in Daout, Sudhaus, Kausch, et al. (2019) (Figures S6c and S7). However, the histogram of differences has an higher standard deviation for Envisat data, which can be explained by the larger noise of this data set due to the poorer temporal sampling in comparison to Sentinel-1 data and to the difference in time span ( Figure S6d). Velocity maps also show some high velocity gradients toward the satellite in many basins bounded by high ranges (Figure S4), where uncertainties are also higher ( Figure S5). Because interferograms have been corrected using the ERA-5 atmospheric model, and because of the good agreement between adjacent ascending tracks and the descending Envisat and Sentine-1 velocity maps, we discard large tropospheric residual signals (>∼1 mm/yr) in the velocity maps that arise from the oddly sampled seasonal stratified delays (Daout, Doin, Peltzer, Lasserre, et al., 2018;Doin, Lasserre, et al., 2009). We associate the apparent correlation with the topography as due to the effect of a long-term increase of water content within the sedimentary basins. Those interferometric synthetic aperture radar (InSAR) observations are in agreement with the water mass increase and gradual expansion of water levels of many lakes of North-East and Central Tibet of the past few decades, amplified by increase in surface air temperature and precipitation, melting of glaciers and thawing of ice-rich permafrost (Bibi et al., 2018;Daout, Dini, Haeberli, et al., 2020).

Forward Modeling
To get insights into deformation mechanisms and compare the observations with first-order dislocation models, we compute two-dimensional surface displacements, for a 20 km-long edge dislocation dipping at 35° to the north and embedded in an elastic half-space, using the formulations provided by Savage (1998) ( Figure S8). Three synthetic surface displacements are plotted for variable locking depths (LD = 0, 5, or 10 km). We choose not to perform any inverse modeling to compare shortening and vertical surface displacement observations with the simplest fault geometry model, as we believe that a distributed fault model could fit almost any kind of surface deformation. We aim at demonstrating in the following sections that at least part of the signal can be explained by anelastic bulk deformation.

Results: Observational Constraints
To first-order, the LOS velocity maps ( Figure S4) show, for both ascending and descending tracks, relative displacements away from the satellite (negative values) from the southern front of the Qilian and Qilian-Nan ranges to the North Qaidam thrust systems (within the North Qaidam ranges), and displacements toward the satellite (positive value) within the Qilian or Qilian-Nan ranges and within the Qaidam Basin, likely associated with the N22°E ∼4-6 mm/yr of shortening across the whole region, derived from GNSS measurements (Liang et al., 2013;Wang et al., 2017). Smaller-wavelength LOS gradients are correlated with Quaternary active fold-and-thrust belts on localized areas of the Qaidam Basin (Yin, Dang, Wang, et al., 2008;Yin, Dang, Zhang, et al., 2008) as across the Lenghu Anticlinorium, the south-dipping Aimunike Thrusts and north-dipping Xitieshan thrusts, or the Olongbulak and Delingha ranges ( Figure S4). Assuming a tectonic origin for the deformation, the velocities are probably related to the slow aseismic shortening and uplift of segmented thrust faults and folds a few tens of kilometers long within the Qaidam basin. However, as the signals here are relatively small, we choose not to interpret them further. In the following sections, we will focus on the two major deformation signals of the velocity map in Figure S4, which are linked to the post-seismic ground deformation following the 2003-2009 Qaidam earthquake sequence (Figure 1).
During the 2014-2019 period, both ascending and descending LOS surface displacements highlight a sharp signal of up to ∼1.5 mm/yr that coincides with the topographic front of the Delingha anticline, in agreement with shallow deformation processes (Figures 2a and 2d). In comparison with the descending view, the peak of the ground deformation in ascending view is shifted a few kilometers southward (Figure 2d). This apparent difference in the peak of velocity gradient is likely due to the ∼NS shortening sensitivity difference between the two geometries, which will be extracted in the following section.
While surface displacements are in good agreement with a logarithmic decay during the 2003-2011 period (as seen from Figure 2c and as modeled by Daout, Sudhaus, Kausch, et al., 2019), Sentinel-1 measurements suggest deformation rates do not follow a 1 t decay, as expected for deformation due to afterslip (e.g., Zhou et al., 2018) ( Figure S9).

Post-seismic Deformation of the 2009 Haixi Earthquake
Decomposition  (Figures 3a and S10). The co-seismic surface displacement pattern does not match the Olongbulak Shan topography ( Figure S10). A spatial pattern similar to the 2009 Haixi earthquake is observed during the 2009-2011 aseismic period with Envisat, but extending slightly further north, as expected from shallow afterslip (Figures 3a and 3c and S10). During this short period after the earthquake, rates are 1 and 6 cm/yr in the hanging wall of the south-dipping fault that broke in 2009, with a sharp and narrow maximum peak in the northwestern part, and − 1-0 cm/yr, in a restricted area of the footwall of the south-dipping thrust (Figure 3). The 2009-2011 short-term post-seismic velocities show an along-strike segmentation of the strain release, which can also be modeled with three south-dipping thrust-fault (Daout, Steinberg, Isken, et al., 2020;Elliott, Parsons, et al., 2011).
The 2014-2019 surface velocity maps in ascending and descending views, show surface displacements patterns that differ from the Envisat co-seismic and post-seismic surface measurements from 2009 to 2011, with lobes toward the satellite matching almost perfectly the recent geological structures of the Olongbulak Shan, observed in satellite imagery (Figure 3a). For the descending geometry, a maximum surface LOS velocity of 4 and 5 mm/yr is measured on the north-western part of the uplifted ranges, at the contact between the Paleozoic bedrock in the north, and the deformed sedimentary sequences of the Olongbulak Shan (Guihua et al., 2013;Yin, Dang, Wang, et al., 2008). The descending western profile there (blue line in profile AA' of Figure 3c), suggests steep and asymmetric surface velocities with 2-3 mm/yr of subsidence in the footwall. Positive surface velocities in the hanging wall of the Sentinel-1 data (blue profile in Figure 3c) are slightly more distributed than in the 2009-2011 Envisat profile (green profile in Figure 3c), with secondary lobes 2.5 km south of the maximum peak. The ascending profile (red profile in Figure 3c), shows an even more distributed pattern, with a slower velocity peak of 2-3 mm/yr, in the northern front, shifted to the south in contrast to the descending geometry, and a ∼8 km-wide LOS pattern, in the NS direction, toward the satellite, which is correlated with the whole Olongbulak pop-up structure. On the two eastern profiles BB' and CC', post-seismic velocities during the 2014-2019 period are wider than the co-seismic surface displacements in both the two views, as for the 2009-2011 post-seismic velocities, but differ from the 2009-2011 period by showing additional short-wavelength peaks of LOS surface velocities toward the satellite in the southern part of the ranges. In other words, in contrast to the Envisat post-seismic measurements that suggest a single lobe of deformation that are in agreement with three shallow high-angle faults, Sentinel-1 observations indicate a more distributed signal, with deformation not only on the shallow northern part of the ranges, which broke during the 2009 event, but also across discrete sections of the southern part of the Olongbulak Shan.

Inferred Structural and Kinematic Models
Both the Olongbulak and Delingha ranges have complicated surface expressions with both north-dipping and south-dipping scarps and with topography modified by erosion. Structural interpretations are, therefore, difficult for both areas, as illustrated by the variety of interpretations found in the literature (Fang et al., 2007;Guihua et al., 2013;Yin, Dang, Wang, et al., 2008). Geodetic measurement of surface displacements represents a good opportunity to show the utility of InSAR in elucidating structural models of fault-and-fold geometries. Decomposition of the LOS velocity signal into a N22°E horizontal and a vertical velocity (Figures 4 and 5, with associated uncertainties in Figures S11 and S12) is possible if we assume no ground deformation in the direction perpendicular to the shortening , which is here a reasonable assumption given the absence of strike-slip motion. From this decomposition, we observe that vertical surface displacements appear to affect the whole Delingha anticline and Olongbulak structure (Figures 4a and 4b and 5a and 5b). Maximum horizontal displacements rates of 0.5-1 cm/yr are higher within these two structures than in the surrounding areas of the North Qaidam thrust systems and higher than the average GNSS shortening rate of ∼4-6 mm/yr, which is distributed across the whole North Qaidam fault system (Liang et al., 2013;Wang et al., 2017) (Figure S1), indicating that measured surface displacement rates represent transient deformation induced by the earthquakes. Given the moderate size of the shallow earthquakes and the short-wavelength of the post-seismic signal, we assume possible viscoelastic relaxation in the lower crust or mantle can be neglected.
The structural cross-section across the Delingha anticline shown in Figure 4c is inferred from the observed signal, our own geomorphological mapping, and available models from Yin, Dang, Wang, et al. (2008), Fang et al. (2007, and Daout, Sudhaus, Kausch, et al. (2019) (Figure S13). The Delingha fold shows sharp contacts between the truncated Pleistocene units on the southern sides of the fold and the Holocene deposits within the Delingha basin, in the footwall (Figure 4b, (Yin, Dang, Wang, et al., 2008)). The exhumed Paleogene soft sediments in the core of the anticline (black arrows in Figure 4b, left) are light-colored in comparison to the harder and more resistant Pleistocene limbs of the anticline. We interpret the folding of the sedimentary beds to be caused by the fault step up over a ramp (Figure 4c) (e.g., Walker, 2006). The geometry of the anticline is controlled by two kink bands, branching where the fault step up (kink 1 in Figure 4c), and at the end of the fault tip (kink 2) (Brandes & Tanner, 2014;Suppe & Medwedeff, 1990). When the fault propagates, the beds roll from a flat position into the steep kink bands. Daout, Sudhaus, Kausch, et al. (2019) modeled the short-term logarithmic after-slip determined using Envisat data from 2003 to 2011 with a 2-dimensional segmented ramp-flat-ramp north-dipping structure that steepens as it approaches the surface (orange line in Figure S13). The absence of the broader long-wavelength signal, ∼10 km north of the high-velocity gradient at the southern topography front, in the Sentinel-1 measurements, suggest that the deformation is essentially concentrated at shallower depth during 2014-2019. Horizontal surface displacements extracted from the 2014-2019 period across the Delingha anticline (blue profile in Figure 4c) are in agreement with movement of the hanging wall toward the south and with a shallow north-dipping thrust fault with steep horizontal gradient localized on the tip of the dislocation (Figure S8). However, an additional localized extension on the southern limb of the anticline is also measured, which is not in agreement with the dislocation models of Figure S8. We suggest that this observation might be the effect of the buckling of the crest, which collapses the southern limb and results in relative shortening across the blind north-dipping fault that is larger than the overall shortening over the region (Figure 4c).
In addition, contrary to what is expected from an elastic dislocation ( Figure S8), no subsidence is observed north and south of the hanging wall (red profile in Figure 4c). Vertical velocity profiles show an asymmetric uplift, with steep velocity gradients in the southern limbs, which coincides with the region with high horizontal motion to the south, a maximum of uplift at the front of the crest, and a gently uplift gradient on the northern dipping limb. Kinematic fault-bend fold models predict that, slip on the dipping ramp produces uplift of the back-limb above the ramp alone, whereas fault-propagation and detachment fold model uplift the crest (Brandes & Tanner, 2014;Suppe, 1983). We propose that the measured vertical velocity with Sentinel-1 data corresponds to the long-term uplift of the fold, which is here consistent with a fold underlain by a north-dipping thrust fault (Brandes & Tanner, 2014;Daëron et al., 2007;Suppe, 1983;Walker, 2006), with a steep but eroded forelimb and a more gently dipping back-limb, as conceptually proposed in Figure 4c. The findings suggest that the 2003 coseismic rupture in the basement induced aseismic fault slip in the sedimentary cover, which has evolved from localized slip shortly after the event (during the 2003-2009 period, approximatively) to a more distributed deformation (during the 2014-2019 period, approximatively) across the anticline. This transient strain release is still ongoing in 2019, more than 16 years after the earthquake.
From our InSAR measurements and the structural cross-section of Yin, Dang, Wang, et al. (2008), Guihua et al. (2013) and Daout, Steinberg, Isken, et al. (2020), we infer for the Olongbulak pop-up, the structural cross-section of Figure 5c. From the inverted shortening rates (Figure 5a), the Olongbulak Shan is mainly moving to the north, with a concentration of the shortening gradient on the northwestern part of the range. However, a movement in the SW direction is also observed on the southwestern part of the system, highlighting the activity of the north-dipping branches, in the southern side of the Olongbulak range and some extension between the northern and southern side of the pop-up structure. Curiously, the N22°E velocity difference across the whole structure (i.e., the shortening) is close to zero (Figure 5c), which is in disagreement with elastic dislocation models ( Figure S8). This absence of shortening motion in the far-field of the structure, in contrast to the observed motion within the structure, may indicate a motion driven from below with the vertical motion. From 2014 to 2019, vertical ground motion is observed across the whole ∼10 kmlarge structure with three maxima that coincided with the three inferred anticlines, numbered (1), (2), and (3) in Figure 5b, from our own geomorphological analysis and bibliographic reviews (Elliott, Parsons, et al., 2011;Guihua et al., 2013;Yin, Dang, Wang, et al., 2008). Differences between ascending track 099 and descending track 004, due to the N22°E shortening (Figure 3a), are mainly observed in the northern part of the pop-up structure, where the LOS surface velocity is also sharp and asymmetric. Surface velocities there may be compatible with shallow aseismic slip on a south-dipping fault ( Figure S8), located on the upward prolongation of the south-dipping and shallow 2009 rupture ( Figure 5c). However, higher rates of surface velocities are located between the basin deformed sedimentary sequences of the Olongbulak range and the Paleozoic bedrock, north of it ( Figure 5b). As we did not identify fault traces or folds in the satellite imagery where the maximum uplift is measured and as two different rock types are identified from the satellite imagery and from the geological mapping of Yin, Dang, Wang, et al. (2008), we interpret the high surface displacements here as due to bedding planes slip along the stratigraphic contacts between a major anticline to the south (marked 1) and the Paleozoic strata to the north. In this complex structure, with imbricated Cenozoic and Mesozoic strata and multiple folds of different ages, old bedding planes might have been reactivated as faults following more recent events, such as the 2009 south-dipping Haixi earthquake (Figure 5c).
In the southern part of the Olongbulak Shan, two secondary uplift lobes are colocated with the two identified southern anticlines, marked (2) and (3) (Figure 5a); small symmetric localized N22°E horizontal components are also observed there. The two folds show truncated geological units with sharp contacts on the southern sides between the uplifted rocks and the more recent deposits on the footwall (Figure 5b). Those two surface expressions are in agreement with two folds associated with shallow north-dipping faults. The exhumed soft sediments in the core of the two anticlines (black arrows in Figure 5b) are light-colored in comparison to the harder and more resistant rocks exposed on the two limbs of the folds (black dashedlines in Figure 5b). Because of differential erosion, the back-limbs look steeper in comparison to the crest and the forelimbs (Figure 5b). Because the surface displacements coincide with the two fold morphologies and because the symmetric shortening signal is not in agreement with what would be expected from slip on thrust faults ( Figure S8), we here again interpret the signal not as the direct result of slip along a north-dipping thrust structure, but as distributed fold-growth deformation of the upper layers induced by slip during the 2008-2009 Haixi events and afterslip following them. The deformation might be accommodated by multiple mechanisms such as off-fault deformation along bedding planes or visco-plastic flexural shear creating mainly vertical uplift in combination with internal shearing with no net shortening across the whole structure. To conclude, the dip-angle between the main south-dipping thrust underlying the Olongbulak pop-up that ruptured during the 2009 Haixi earthquake (conceptual black arrow in Figure 5c, insert) and the horizontal shortening (conceptual blue arrow in Figure 5c, insert) creates a reverse north-dipping component of motion that must be accommodated during the seismic cycle either by seismic or aseismic deformation (conceptual red arrow in Figure 5c, insert). The measurements suggest that this vertical partitioning is, at least in part, compensated by aseismic post-seismic folding across the whole structure and localized along the anticlines observed in the geomorphology.

Aseismic Transient Phenomena
Since the discovery at different crustal depths of aseismic transient phenomena from observations made using geodetic networks, large efforts have been made by the geodetic and earthquake communities to better understand and locate these events occurring at low rupture speed and to improve the understanding of their interaction with earthquakes (e.g., Behr & Bürgmann, 2020;Collettini et al., 2011;Hawthorne et al., 2016;Obara & Hirose, 2006;Rousset et al., 2016;Socquet et al., 2017). The findings of this study illustrate that moderate earthquakes can trigger anelastic transient responses of the shallow crust, which lasts for decades. It shows that visco-plastic deformation can be produced in low strain regimes by moderate stress changes and at shallow depths. These observations may support weak fault zone models, documented both at shallow (e.g., Niemeijer & Spiers, 2005;Collettini et al., 2011) and greater depths (e.g., Behr & Bürgmann, 2020), where velocity-strengthening frictional behavior is enhanced by the presence of weak materials, such as phyllosilicate-rich meta-sediments or sediments, within heterogeneous fault zones. In those models, weak rocks or viscous matrix surround mafic and more competent lenses, and may accommodate the transient deformation by pressure solution, foliation-parallel sliding, viscous deformation and micro-folds (e.g., Collettini et al., 2011;Niemeijer & Spiers, 2005;Lavier et al., 2013;Reber et al., 2014). This study suggests that slow and transient aseismic events, detected both at shallow and greater depths during the inter-or post-seismic periods, might not be always linked with slip on a thin main fault or slab driving elastic deformation of the surrounding medium. On the contrary, they might also be associated with anelastic response of the medium itself, either triggered by an earthquake on the main adjacent fault-segment (post-seismic transients), or by a transient increase in stress concentrations within the more competent rocks of an heterogeneous shear zone (slow-slip transients). The increased number of InSAR satellite missions and the improved temporal resolution of new satellites open wide perspectives for accessing the lag between the co-seismic forcing and the anelastic response of the crust, and therefore better estimating the effective viscosity controlling those mechanisms.

Conclusion
The improved spatial and temporal sampling offered by the new generation of Sentinel-1 radar data in both ascending and descending geometries, has enabled us to resolve the vertical and shortening component across two fold-and-thrust belts in the Qaidam Basin, China. The study reveals a link between the post-seismic phase and the long-term geomorphology of the tectonic structures. We interpret the non-logarithmic long-term ground aseismic surface displacements observed up to 16 years after the seismic events as offfault volumetric and anelastic deformation associated with active folding. In the two examples, co-seismic slip is translated into permanent deformation during the post-seismic phase through localized signals that match young anticlines observed in the morphology, and that cannot be simply explained by slip along faults driving elastic deformation of the medium. The results demonstrate, along with other observations and models (e.g., Ainscoe et al., 2017;Bonanno et al., 2017;Johnson, 2018), that it not always valid to assume that fold growth above blind reverse faults is co-seismic or that it is appropriate to model it by slip on a fault, because a significant part of the finite shortening may occur as distributed off-fault deformation (flexural slip, interbed flow, shearing) during the post-seismic or inter-seismic phases. The study highlights that the contribution of distributed deformation and slip-partitioning between secondary structures needs to be considered in earthquake cycle models, analyses of fault-related folds, and the assessments of earthquake hazard from surface measurements.