MT-InSAR Unveils Dynamic Permafrost Disturbances in Hoh Xil (Kekexili) on the Tibetan Plateau Hinterland

Hoh Xil is an uninhabited extremity secluded on the Tibetan Plateau hinterland. A complete mapping of ground motion variation in Hoh Xil is essential for an in-depth understanding of the terrain’s responses to climate change on the Tibetan Plateau. However, the inaccessibility and extremely harsh environment impeded extensive field investigations on landform alteration and its formative process. Such difficulty can be resolved by interferometric synthetic aperture radar (InSAR), which enables a broad detection of subtle permafrost motions at millimeter precision. This study, for the first time, accomplished a multitemporal InSAR (MT-InSAR) deformation mapping from 2015 to 2020 in Hoh Xil, with a wide coverage of about 200 000 km2. 1592 Sentinel-1 images were processed based on the small baseline subset (SBAS) technique. The results show that Hoh Xil was experiencing dynamic permafrost disturbances. Thawing permafrost with both a linear subsidence rate higher than 2 mm/year and a periodic amplitude over 2 mm was primarily detected in areas of flat or gentle slopes. The InSAR cumulative deformation is highly correlated with permafrost thawing depth. Significant lag times were identified between seasonal oscillation of InSAR deformation and land surface temperature (LST). Thermokarst landforms of retrogressive thaw slumps and thermokarst lakes broadly formed and dynamically evolved as a consequence of permafrost degradation. Particularly, widespread thawing permafrost characterized by the spatial clustering of thermokarst lakes appeared to occur in areas adjacent to large lakes. The discovered dynamic permafrost disturbances in Hoh Xil manifested even the secluded Tibetan Plateau hinterland was facing the threat of climate change.

Abstract-Hoh Xil is an uninhabited extremity secluded on the Tibetan Plateau hinterland. A complete mapping of ground motion variation in Hoh Xil is essential for an in-depth understanding of the terrain's responses to climate change on the Tibetan Plateau. However, the inaccessibility and extremely harsh environment impeded extensive field investigations on landform alteration and its formative process. Such difficulty can be resolved by interferometric synthetic aperture radar (InSAR), which enables a broad detection of subtle permafrost motions at millimeter precision. This study, for the first time, accomplished a multitemporal InSAR (MT-InSAR) deformation mapping from 2015 to 2020 in Hoh Xil, with a wide coverage of about 200 000 km 2 . 1592 Sentinel-1 images were processed based on the small baseline subset (SBAS) technique. The results show that Hoh Xil was experiencing dynamic permafrost disturbances. Thawing permafrost with both a linear subsidence rate higher than 2 mm/year and a periodic amplitude over 2 mm was primarily detected in areas of flat or gentle slopes. The InSAR cumulative deformation is highly correlated with permafrost thawing depth. Significant lag times were identified between seasonal oscillation of InSAR deformation and land surface temperature (LST). Thermokarst landforms of retrogressive thaw slumps and thermokarst lakes broadly formed and dynamically evolved as a consequence of permafrost degradation. Particularly, widespread thawing permafrost characterized by the spatial clustering of thermokarst lakes appeared to occur in areas adjacent to large lakes. The discovered dynamic permafrost disturbances in Hoh Xil manifested even the secluded Tibetan Plateau hinterland was facing the threat of climate change.

I. INTRODUCTION
T IBETAN Plateau, recognized as the world's "Third Pole," is the highest and largest elevated plateau on the Earth. The plateau is sensitive to climate change, which has been Manuscript  constantly occurring since the 1950 s [1], [2], [3]. Hoh Xil, also known as Kekexili, is a secluded extremity situated on the Tibetan Plateau hinterland. The average elevation of the region is about 4800 m above sea level, with scarce precipitation and low temperature as its climate features. Hoh Xil is among one of the most vulnerable areas on the Tibetan Plateau to climate change, thus making it a critical area for an in-depth understanding of the terrain's responses to climate change. Tibetan Plateau has a wide permafrost coverage of about one million km 2 [4], [5]. Permafrost disturbances, which refer to the state of permafrost being disturbed, chiefly include the formation of thermokarst lakes, retrogressive thaw slumps, deepened thawing depth, wildfire, thermoerosion, and anthropogenic activities. Climate warming may considerably decrease ground ice stability and accordingly induce permafrost disturbances, thereby dramatically changing the hydrological cycle and hydro-geomorphology on the Tibetan Plateau [6]. Wu et al. [7] argued that cold permafrost (<−1.5 • C) appears to be more sensitive to climate change and human-caused disturbances than warm permafrost (>−1.5 • C). Permafrost degradation characterized by the deepened thawing depth and accelerated ground deformation is prevalent on the Tibetan Plateau, thus resulting in accelerated desertification and environmental deterioration [3], [8], [9], [10], [11]. Liu et al. [12] revealed an increasing trend of active layer thickness (ALT) and seasonal thaw depth over extensive areas on the Tibetan Plateau since 2003. Permafrost disturbances also have a strong control over soil water regimes on the Tibetan Plateau [13]. In the Hoh Xil region, climate warming has raised permafrost temperature and reformed freeze-thaw processes, although such responses of permafrost warming to climate change may be slow and lagged [7], [14]. Ground ice content within permafrost appears to be higher in Hoh Xil than in other areas on the Tibetan Plateau [14]. Precipitation events could be the principal source of water within the active layer of permafrost in Hoh Xil [15]. Thermokarst lakes and thaw slumps, which are the primary evidences of disturbances on unstable and warming permafrost, have been gradually noticed on the Tibetan Plateau [16], [17], [18], [19], [20], [21].
Interferometric synthetic aperture radar (InSAR) is a compelling microwave remote sensing approach to measuring wide-ranging ground motions potentially at millimeter precision. In recent years, there is a growing awareness of the usefulness of differential InSAR (DInSAR) and multitemporal This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ InSAR (MT-InSAR) in extensive permafrost investigation [22]. Liu et al. [23] tentatively justified that the MT-InSAR is capable of capturing seasonal and secular subsidence induced by permafrost thawing, and accordingly provides a complementary tool to conventional long-term in situ measurements of ALT with high spatial resolution over wide areas. Such potential can be further improved by means of multisensor and multigeometry of satellites to capture variable changes in permafrost landscapes with decomposed multidimensional patterns [24]. Short et al. [25] demonstrated the effectiveness of DInSAR stacking in estimating seasonal movements in permafrost landscapes and extracting thaw settlement properties of surficial geology. Provided that core samples are available, it is feasible to establish correlations between InSAR-derived maximum seasonal thaw subsidence and soil water content of active layers [26]. Land surface temperature (LST) can be integrated into InSAR deformation models to further refine seasonal and long-term products [27]. In specific conditions, the longer wavelength L-band InSAR, which potentially provides better coherence and penetration capabilities, may further improve an in-depth understanding of permafrost evolution and degradation as well as thermokarst dynamics [26], [28], [29], [30], [31]. InSAR is also essential for monitoring surface displacements over diverse landscape types in discontinuous permafrost terrain [32]. Thermokarst subsidence after ice-rich permafrost thawing may be captured by L-band interferograms [33]. Furthermore, a number of studies corroborated the usefulness of long-term InSAR in detecting ground motion caused by seasonal permafrost freezing and thawing cycles in difficult terrain conditions on the Tibetan Plateau [27], [34], [35], [36], [37], [38], [39]. In particular, Zhao et al. [35] refined the InSAR deformation model with a consideration of temperature and precipitation to better fit permafrost motions. It turns out that slopes exhibit a consequential influence on the amplitude of permafrost's seasonal deformation on the Tibetan Plateau [27]. As seen from the InSAR-derived ground motions corrected by the ERA-5 atmospheric model, seasonal thawing can exert a progressive influence on ice-rich layers at the permafrost table while rapid mass movements may promote the flattening of the Tibetan Plateau [40].
Notwithstanding the productive InSAR applications in permafrost investigation on the Tibetan Plateau, the entire Hoh Xil region has never been completely mapped. Thus, a thorough understanding of ground motion variation with permafrost disturbances is still lacking. Chen et al. [29] and [30] combined the C-band ENVISAT ASAR images and L-band ALOS PALSAR images to investigate surface displacements restricted in the Beiluhe area, and their studies revealed the apparent anthropogenic influences on permafrost environment. Similarly, in Beiluhe, Wang et al. [41] and [42] processed the TerraSAR-X images for up to three years (2014-2016) to trace seasonal deformation across the Qinghai-Tibet railway and subsequently made effort to retrieve ALT. Lu et al. [43] deployed a pointwise approach of StaMPS-InSAR to extract seasonal deformations caused by permafrost in the Wudaoliang Basin, which only represents a small portion of Hoh Xil. Li et al. [44] monitored periodic and seasonal deformation based on a joint analysis of InSAR and global navigation satellite system (GNSS) leveling series, however, in a small portion of the Wudaoliang Basin. Focusing on a specific natural hazard event, the phenomenon that lake outburst may accelerate surrounding permafrost degradation in Hoh Xil was revealed with the aid of InSAR [45], [46]. Considering the extremely harsh environment in Hoh Xil, it is extremely difficult to perform extensive in situ observations on permafrost thawing-freezing cycles throughout the whole region. Therefore, currently, InSAR remains the only approach to widely map ground deformation and investigate permafrost dynamics in Hoh Xil.
This study, for the first time, releases an MT-InSAR deformation map spanning from 2015 to 2020 across the entire Hoh Xil region with a wide coverage of about 200 000 km 2 . All the InSAR processing based on the 1592 Sentinel-1 images was carefully reviewed for the quality of the derived result of ground deformation. Such MT-InSAR deformation product was further used to investigate spatial variations and seasonal dynamics of permafrost disturbances throughout the whole territory. The novelties of this study chiefly include the following two aspects: 1) by integrating the in situ measurements of multilayer soil temperature, the correlation between InSAR ground deformation and permafrost thawing depth was in-depth investigated and discussed and 2) by a joint analysis of optical imageries and MT-InSAR, the dynamic evolution of typical landform indicators associated with permafrost thawing and degradation, particularly the thermokarst lakes and retrogressive thaw slumps, were broadly revealed. Considering that permafrost in Hoh Xil accounts for approximately 20% of that on the Tibetan Plateau, the derived MT-InSAR products may be useful for promoting a complete estimation of ALT on the Third Pole. To help achieve the goal for better public sharing of scientific data on the Tibetan Plateau and in China [47], the produced MT-InSAR deformation map of Hoh Xil will be shared with National Tibetan Plateau/Third Pole Environment Data Center (TPDC) with the public access.

II. STUDY AREA
Hoh Xil (also known as Kekexili, Fig. 1) is a highland above 4800 m above sea level, connecting Tanggula Mountain in the south and Kunlun Mountain in the north. It is a wide unpopulated land, thereby extremely impervious to modern human impacts. The only transportation connecting the region is the quasi-parallel Qinghai-Tibet railway and highway, stretching along the eastern margin of Hoh Xil (see Fig. 1).
With regard to geology, the Late-Cenozoic volcanoes are widely found in Hoh Xil, and the pyroxenite xenoliths in these volcanic rocks comprise clinopyroxenes and orthopyroxenes [48]. In particular, the shoshonitic to ultrapotassic volcanic rock series are the predominant categories across the entire region [49]. The age of igneous rocks in Hoh Xil was estimated spanning from ca. 217 to 202 Ma [50], [51]. At a regional scale, Hoh Xil is located in the largest Cenozoic sedimentary basin (with a fill thickness greater than 6000 m) [52], [53], whose evolution may show pivotal implications in resolving the formation mechanism and deformation history of the plateau. The Jinsha Suture trends in the East-West direction throughout the Hoh Xil basin which stretches beyond the Qiangtang and Songpan-Ganzi terranes [54]. Zhang et al. [55] argued that the Hoh Xil arc is split from the Kunlun arc, a long-term active magmatic belt in northern TP. Liu and Xia [56] doubted such geological correlation based on the fact that the Hoh Xil arc is primarily overlaid by the Triassic turbidites and no older magmatic rock was recorded.
In terms of hydraulic connection, at the regional scale, the Hoh Xil region is situated at the junction of the northern source the Yangtze River drainage system and Qiangtang Inland Lake region [57]. Lakes are abundantly distributed across the entire region. A majority of non-glacier-fed lakes in Hoh Xil were facing exceptional expansion chiefly as a consequence of increased precipitation [13]. As an extreme event jointly affected by heavy precipitation and glacier melting, on September 14, 2011, the outburst of Zonag Lake (headwater lake) occurred and consequently change the hydraulic connection of Salt Lake (tailwater lake) [45].
The whole Hoh Xil region is prevailed by the semi-arid continental climate, featured by low temperature (with the mean annual air temperature varying from −10.0 • C to −4.1 • C) and high evaporation demand, as well as low precipitation (annually 70.5-291.4 mm) progressively reducing from southeast to northwest, thus generating major natural land covers as alpine meadow, alpine steppe, and alpine desert [6], [57]. Periglacial landforms are prevalently formed across the entire region [58].
The Hoh Xil region is predominantly underlain by continuous permafrost [4], [59]. Higher ground ice content of permafrost was measured in Hoh Xil than in other zones on the Tibetan Plateau [14]. Under the circumstances of regional climate change, the thawing island may individually form within such warm and ice-rich continuous permafrost in which geothermal heat regularly circulates [14], [29]. Zhao et al. [14] observed that from 2004 to 2017, the ALT in Hoh Xil increased from 1.60 to 1.76 m, as a result of increasing melt water. In addition, Zhu et al. [15] showed that an 8.1 mm precipitation in Hoh Xil contributed 49.0% and 30.8% water at depths of 0-10 cm and 10-20 cm within the active layer, respectively. Another source of disturbance to continuous permafrost in Hoh Xil is due to the warm lake bottoms resulting from the deeper water depth than winter ice, thus accelerating permafrost degradation and further forming and enlarging thermokarst lakes; the spatial distribution of these thermokarst lakes is tightly associated with ground temperature and more preferably ice content within permafrost [16], [60]. Moreover, retrogressive thaw slumps, which are indicators of unstable slopes of ice-rich permafrost, were found clustered especially in the Beiluhe region [17], [19], [20]. Lu et al. [45] revealed that a lake outburst event may significantly accelerate regional permafrost degradation, as a consequence of the thermal alteration of permafrost thawing-freezing cycle in addition to changes in hydrological connectivity and soil permeability. After this event, new permafrost had been forming on the Zonag Lake bottom and permafrost degradation even accelerated under the Salt Lake bottom [61], [62].

III. DATASETS AND METHODOLOGY
Firstly, the MT-InSAR processing procedures are clarified in detail. Secondly, the approach to estimating the InSAR linear velocity of ground deformation and periodic amplitude of seasonal deformation is explained. Thirdly, the acquisition of in situ multilayer soil temperature and the estimation of thawing depth and seasonal oscillation of LST are expressed. Fourthly, the prepared inventories of thermokarst lakes and retrogressive thaw slumps are presented.

A. MT-InSAR Processing
A total of 1592 Sentinel-1 synthetic aperture radar (SAR) images covering seven ascending tracks (Tracks 12, 41, 70, 85, 114, 143, and 158), spanning from May 2015 to August 2020, were acquired in the interferometric wide-swath (IWS) mode. The method of small baseline subset (SBAS) [63], [64], which can minimize the effects of spatial and temporal baseline by deploying multiple master image interferograms, was used for the MT-InSAR processing. Most spatially-correlated atmospheric phases were manually eliminated by the multiple processing of regression and low-pass spatial filtering.
The first step was to generate unwrapped differential interferometric sets. The entire Hoh Xil region was divided into 23 regions for the processing (given in the Appendix, Fig. 14). Hence, a total of 23 InSAR time series were, respectively, processed. The Sentinel-1 SAR images in the same area were processed with precise orbit correction and registration. The interferograms with short temporal baselines ranging from 12 to 72 days between adjacent images were generated and they were multilooked to reduce noises (eight-looks in range and two-looks in azimuth). A 30 m digital elevation model (DEM) of the Shuttle Radar Topographic Mission (SRTM) was then used to remove contributions of terrain phases. About 2000 low-quality interferograms affected by significant decorrelation and strong atmospheric interference were rejected through manual checking. After removing the low-quality interferograms, the baseline network can be still fully connected by the remaining ones rather than consisting of several sets with a temporal gap, and accordingly, it is not difficult to retain a sufficiently constrained graph. After applying the Goldstein-Werner adaptive filter on the remaining differential interferograms, the method of minimum cost flow (MCF) was exploited to unwrap the differential interferometric phases.
Regarding the selection of 23 reference points (given in the Appendix, Fig. 14), due to the lack of broad field surveys and remote sensing imageries meeting sufficient spatial resolution requirements, it is difficult to precisely identify the exposed bedrocks which are usually chosen as stable points on the Tibetan Plateau. As a result, we followed the two rules to manually select the possibly stable reference points, as suggested by Reinosch et al. [38]: 1) the preference is given to selecting pixels with high coherence. The coherence in the Hoh Xil region usually decreases significantly between May and October each year. The pixels with sufficient coherence (>0.8) during this period were selected as the candidate reference points and 2) the reference points were further selected in areas far away from alluvial fans, river/lake bench, valley, and typical thermokarst landforms such as thaw slumps, both of which are considered unstable. This was accomplished by manually checking if there is any change from Plan-etScope CubeSat images acquired in 2016 and 2020. A total of 153 PlanetScope CubeSat images covering all 23 reference points were manually checked. All the reference points were thus selected with no significant geomorphic change observed on the bi-temporal PlanetScope CubeSat images. Following these two rules, may maximally reduce the possibilities that the selected reference points are unstable.
Differential interferograms with unwrapping failures were manually checked and reprocessed by the parameter tuning of the Goldstein-Werner adaptive filter. Those differential interferograms requiring additional reprocessing were mainly during May-July and September-November. Two parameters of the Goldstein-Werner adaptive filter were tuned. The first parameter is the exponent for the non-linear filtering, which was attempted as 0.4 and 0.6. The second parameter is the filtering FFT window size, which was set as 32, 64, and 128, respectively, for tuning until the filtered interferograms show no significant phase jump while unwrapping. We expect such parameter tuning may obtain a filtered interferogram with a relatively continuous phase in the periglacial regions, which usually demonstrate sharp boundaries in seasonal amplitude and can be thus difficult to unwrap.
The unwrapped differential interferometric phase consists of four components and the phase model in each interferogram can be represented as where ϕ i, j represents the phase of the ith pixel in the jth unwrapped differential interferogram, generated by the SAR images acquisitions in t 2 and t 1 . t j equals the time interval of t 2 and t 1 . v j is the linear deformation rate of t j and λ is the wavelength. B ⊥ j is the perpendicular baseline. The DEM height error is z i . The slant range and look angle are ρ i and θ i , respectively. φ atm i (t 2,1 ) accounts for the atmospheric phase component and n i (t 2,1 ) expresses the phase contribution due to decorrelation and thermal noises [63].
The next step was to mitigate the residual phases due to DEM height error and long-wavelength atmospheric delay by performing the time series analysis. According to (1), the slope of the phase ϕ i, j to baseline B ⊥ j dependence can be exploited to calculate the correction value for the height error z i . The slope of the phase ϕ i, j to time t j dependence can be used to estimate the linear deformation rate v j . The regression model of the time series analysis can thus be simplified as where a 0 and a 1 are the coefficients related to the annual linear deformation rates and the height correction values, respectively [45]. In the first regression, a 1 was estimated to simulate the phase components due to height error. After removing the residual terrain phases and updating the time series, a 0 and a 1 were re-estimated in the second regression. Then those pixels with slow deformation rates a 0 were selected to calculate the spatially-correlated atmospheric phases. The atmospheric phases were estimated and removed by the spatial low-pass filter and spatial interpolation. The window size of the spatial filter is fixed rather than the adaptive one. Local turbulent atmospheric phases might be incompletely estimated by the spatial filter with a fixed size. Therefore, the time series could be potentially affected by the residual turbulent atmospheric phases. Those low-velocity points were selected again in the third regression. The residual atmospheric phases were then estimated by the spatial filter and interpolation on these low-velocity points.
The last step was to estimate the line of sight (LOS) deformation. The time series of phases were still a multireference set after removing the residual terrain and atmospheric phases. Hence, the singular value decomposition was carried out to convert the multireference set to a single-reference time series, which was referenced to the time of the earliest scene in the series [64]. The surface deformation along LOS can be consequently derived from the phase signals.

B. Estimation of Linear Velocity and Periodic Amplitude of Ground Deformation
In order to estimate the linear velocity and periodic amplitude of ground deformation from the derived InSAR time series, a combination of linear and sinusoidal functions was exploited. For the ith pixel of mth images in the time series, the deformation observation derived by MT-InSAR was defined as where ω, v and d represent the averaged period of seasonal variation, annual velocity, and the constant term, respectively. The square root of a 2 plus b 2 equals to the seasonal oscillation amplitude (Amp) t m is the time interval between the SAR image m and the reference image (i.e., the first SAR image). The five parameters of a, b, v, d, and ω were iteratively estimated by the least squares method. The first-order Taylor series expansion was exploited to linearize (3) as where f (P 0 ), f ′ (P 0 ), and P are the function value, the firstorder derivative value at the approximation value P 0 of the five parameters, and the correction of the five parameters estimated by the least squares method, respectively. Initial values of the first four parameters were determined after performing a linear regression in (3), which took 0.0172 as the initial value of ω (the initial period is 365 days for a freezing-thawing cycle). After the linearization, the deformation observation of m images was further expressed as (6), shown at the bottom of the next page. Then, the five parameters were updated by adding the correction values P after each estimation and were brought into the next execution. Once P arrived at a given threshold of 0.001 and converged, the iterated least squares estimation stopped. The periodic amplitude in a freezing-thawing cycle was assumed as the peak-to-peak amplitude [65], [66]. It is the difference between the maximum uplift and subsidence after removing the linear annual contribution, namely the distance between the positive and negative peaks of: The uncertainty of estimated parameters was related to the diagonal entries of the variance-covariance matrix, when the correlation between the five parameters was ignored. The variance-covariance matrix was derived from the residuals and the Jacobi matrix after the last estimation [27], [36], [37], which was calculated by Finally, the estimation results were converted and resampled into raster images with a spatial resolution of 250 m, and each pixel value represents the mean value of all points within it.

C. Estimating Thawing Depth and Seasonal Oscillation of LST
The daily multilayer soil temperature datasets from the in situ measurements at three stations (see their locations in Fig. 1) were collected to interpolate the soil temperature profile and estimate the thawing depth [14], [67]. The thawing depth in this study equals the layer where the temperature is at 0 • in a thawing period. Due to the limited temporal coverage of the in situ observations, the daily 1-km all-weather LST datasets for the Western China (TRIMS LST-TP [68], [69], [70]), which completely covers our InSAR observation period, was exploited to compare the correlation between the LST and the InSAR-derived deformation. The daily LST used in this study was the average of daytime and nighttime products provided by TRIMS LST-TP. Similar to the estimation of seasonal deformation, a combination of linear and sinusoidal functions was used to estimate the seasonal signals of LST from January 1, 2015, to December 31, 2020.

D. Preparing Inventories of Thermokarst Lakes and Retrogressive Thaw Slumps
The thermokarst lake inventory on the Tibetan Plateau produced by Wei et al. [21] was used in this study. The inventory shows that there are totally 54 496 thermokarst lakes in the Hoh Xil region, accounting for about 33.8% of those on the entire Tibetan Plateau. In addition, the inventory of retrogressive thaw slumps was manually interpreted and mapped by comparing the bi-temporal PlanetScope CubeSat images acquired in 2016 and 2020. Those thaw slumps usually have well-defined headwalls due to recently exposed soil layer and poor vegetation coverage [17], which made them clearly visible in the images. The selected PlanetScope CubeSat images were true color composited, and the boundaries of thaw slumps were delineated by manual vectorization to maximally ensure the mapping quality of this inventory.

A. Linear Velocity and Periodic Amplitude of Ground Deformation
Typical permafrost disturbances in Hoh Xil (i.e., thermokarst lakes, retrogressive thaw slumps, and permafrost degradation) will cause surface subsidence, which can be captured by InSAR at millimeter precision. Two types of ground deformation maps were produced across the entire Hoh Xil region: 1) the linear velocity of ground deformation [see Fig. 2(a)] and 2) the periodic amplitude of seasonal deformation [see Fig. 2(b)]. The linear velocity of ground deformation is able to represent the rate of permafrost thawing subsidence while the periodic amplitude of seasonal deformation can monitor permafrost freeze-thaw cycles, namely the seasonal thaw subsidence and freeze uplift. The time span of the linear velocity map is from March 2015 to August 2020. Surface displacements were calculated along the LOS. Positive and negative velocities indicate ground motions moving toward and away from the satellite, respectively. As seen from this map [see Fig. 2(a)], it appears that the ground subsidence (i.e., velocity <0, pixels in red) can be prevalently detected across the entire Hoh Xil region. The map statistics indicate that approximately 7.9% of the total area in Hoh Xil was undergoing dynamic ground subsidence with a subsiding rate estimated over 5 mm/year. The uncertainty of model fit variances for the linear velocity is displayed in Fig. 2(c). It turns out that the average uncertainty of linear velocity is about 0.09 mm/year across the entire Hoh Xil region. The mean periodic amplitude of seasonal deformation in the Hoh Xil region is about 2 mm, with an average uncertainty of model fit variances estimated as about 0.90 mm [see Fig. 2(d)]. Approximately 37.6% areas show periodic amplitude of seasonal deformation larger than 2 mm, and 5.8% areas in Hoh Xil exhibit more dynamic seasonal oscillations with the derived peak-to-peak amplitude greater than 5 mm. Those areas with such large seasonal deformation are mainly clustered in the northeastern and mid-northern areas in Hoh Xil. The thawing permafrost, which is defined as with both the linear subsidence rate higher than 2 mm/year and the periodic amplitude over 2 mm in this study, accounts for about 12.1% of the total area in Hoh Xil.

B. Correlations Between InSAR Deformation and Topography
The kernel density plots, as deployed by Daout et al. [36], were drawn to investigate the correlations between the InSARderived ground deformation and the topographic factors in  Hoh Xil (see Fig. 3). The estimated probability density values were normalized to 0-1. Larger (i.e., closer to 1, dark red) and smaller values (i.e., closer to 0, deep blue) represent higher and lower estimated density, respectively. It turns out that the linear velocity of ground deformation is not sensitive to elevation changes [see Fig. 3(a)]. Instead, an obvious correlation can be observed between linear velocity and slope. Namely, when the slope increases, the velocity of freezing or thawing process decreases [see Fig. 3(b)]. Similarly, the periodic amplitude of seasonal deformation was found to be more correlated with variations of the slope. In particular, the seasonal amplitude shows a tendency to decrease when the slope increases. Daout et al. [36] indicated that such spatial variation of seasonal deformation amplitude was determined by the spatial variability of ALT, and the increasing magnitude and variations of seasonal deformation were primarily caused by the ALT thickening. Our observations verified that the thawing permafrost was primarily distributed in areas of flat or gentle slopes. Such a finding is also consistent with the observations on the Tibetan Plateau conducted by Chen et al. [27] and Daout et al. [36].

C. Correlations Between InSAR Deformation and Thawing Depth
Deepened thawing depth is a typical manifestation of permafrost disturbances. The thawing depth outcomes estimated from the in situ measurement of multilayer soil temperature were resampled according to the acquisition date of the Sentinel-1 images during the permafrost thawing period, namely from May to September (see Fig. 4). The estimated maximum thawing depth is 388, 240, and 214 cm for the KKXL, QT08, and Wudaoliang stations, respectively. The results of the InSAR-derived cumulative deformation with reference to the first Sentinel-1 image in May were also calculated for each thawing period (see Fig. 4). The maximum cumulative deformation is 8.1, 5.0, and 8.4 mm at . . .

cos(ω 0 t m )
Higher correlation between cumulative deformation and thawing depth was achieved, which is of 0.99, 0.92, and 0.99 at the KKXL, QT08, and Wudaoliang stations, respectively (see Table I).
We also noticed that there are significant lag times between seasonal oscillation of the InSAR-derived deformation and the LST at both the KKXL, QT08, and Wudaoliang stations (see Fig. 5). The maximum subsidence has a significant delay compared to the maximum LST at these three stations (see Fig. 5). The average lag times between the maximum LST and subsidence are 80, 120, and 86 days for the KKXL, QT08, and Wudaoliang stations, respectively (see Table I). Previous studies have revealed the existence of lag time between seasonal oscillation of air temperature and InSAR deformation over permafrost areas [32], [34], [36], [38], [43]. Our study further confirmed such lag time also may exist between the maximum LST and subsidence. This finding also justified why the InSAR deformation did not equal permafrost thawing depth.

D. Accumulation of Thermokarst Lakes on Thawing Permafrost Close to Large Lakes
Based on the investigation of the thermokarst lake inventory on the Tibetan Plateau, we evaluated that the entire Hoh Xil region has a total of 54 496 thermokarst lakes, with a complete area of 1413 km 2 (on average 0.026 km 2 for each lake). In particular, we discovered that thermokarst lakes on thawing permafrost tend to evolve and accumulate close to large lakes in Hoh Xil. 34 285 thermokarst lakes (approximately 62.9% of the total number of thermokarst lakes) were found distributed within a 20 km distance around the 52 large lakes (with an area exceeding 50 km 2 ) in Hoh Xil. The spatial density of these thermokarst lakes is illustrated in Fig. 6, from which such clustering of thermokarst lakes close to large lakes can be visualized.
Four representative areas in which thermokarst lakes spatially clustered were identified in the vicinity of Xuejing Lake, Ulan Ul Lake, Zonag Lake, and Salt Lake [see Fig. 7(a)]. Subsidence was prevalently identified in these four areas during 2015-2020, as shown by the InSAR-derived linear velocity of ground deformation [see Fig. 7(b)] and the periodic amplitude of seasonal deformation [see Fig. 7(c)]. 4570 thermokarst lakes, accounting for 83.35% of the total number (5483) in these four areas were found with peripheral subsidence, with a maximum rate of 10 mm/year. Substantial subsidence was primarily detected surrounding Salt Lake and Ulan Ul Lake [see Fig. 7(b)]. Such considerable subsidence is associated with the significant expansion of these two lakes. During 2015-2020, Salt Lake and Ulan Ul Lake had expanded 57 and 36 km 2 , respectively. In contrast, Xuejing Lake experienced a relatively small expansion of 8 km 2 during 2015-2020, and its surrounding areas were barely undergoing outspread subsidence as seen from the linear velocity of ground deformation [see Fig. 7(b)]. It appears that the spatial distribution of  thermokarst lakes is densest near Xuejing Lake, thus indicating this area had experienced permafrost degradation at an earlier stage. As a result, lower linear velocity and periodic amplitude were detected near Xuejing Lake. In terms of Zonag Lake, it shrunk 7 km 2 during 2015-2020, however, due to the consequence of an outburst event that occurred in 2011 [45]. We further selected 12 check regions [see Fig. 7(c)] on the periodic amplitude results for further time series investigation (see Fig. 8). The averaged linear velocity and periodic amplitude were subsequently estimated (see Table II). Apparent periodic characteristics of permafrost freezing-thawing process can be observed from the deformation time series for all  12 check regions (see Fig. 8), showing an evident trend of thawing (subsidence) starting in April and May and freezing (uplift) beginning in October and November. The estimated duration of these freezing-thawing cycles for all the 12 check regions is listed in Table II. As seen from the linear velocity, subsidence was observed for all 12 check regions with a rate ranging from 3.1 to 10.4 mm/year, indicating that permafrost thawing has been occurring in those clustered areas of thermokarst lakes. The most dynamic permafrost disturbance was detected close to Salt Lake, where substantial periodic amplitudes of 5.9, 6.7, and 9.5 mm were observed for CR10, CR11, and CR12, respectively. On the whole, our observations demonstrated that the dynamic evolution of thermokarst lakes in Hoh Xil may be formed as a joint consequence of permafrost degradation and the accelerated partial drainage of adjacent large lakes.

E. Dynamic Evolution of Retrogressive Thaw Slumps
Based on the derived InSAR deformation map, we identified two representative areas of accumulated retrogressive thaw slumps in Hoh Xil, namely the Chumar River area and the Beilu River area (see Fig. 9). These two areas are characterized by both the increased values of linear velocity [see Fig. 9(a)] and periodic amplitude [see Fig. 9(b)]. In particular, a max-  III   AVERAGED LINEAR VELOCITY AND PERIODIC AMPLITUDE CALCULATED  FROM THE INSAR-DERIVED LOS DEFORMATION FOR RETROGRES-SIVE THAW SLUMPS MAPPED WITHIN THREE CHECK REGIONS  (CR13-CR15) imum value of accumulated subsidence of 141 and 108 mm during 2015-2020 was detected for the Chumar River area and the Beilu River area, respectively. In addition, a maximum periodic amplitude of 15.3 and 14.8 mm was discovered for the Chumar River area and the Beilu River area, respectively. We also prepared an inventory of thaw slumps for these two areas through the manual inspection on the bi-temporal PlanetScope imageries acquired in October 2016 and October 2020. The spatial distribution of these manually mapped thaw slumps is demonstrated in Fig. 9(b). We further selected three check regions (CR13-CR15 in Fig. 9) on the linear velocity maps for further time series investigation (see Fig. 10). CR13 is located in the Chumar River area while CR14 and CR15 are situated in the Beilu River area (see Fig. 9). All these three check regions were identified with many thaw slumps as seen from the manually prepared thaw slump inventory. The boundaries of these thaw slumps were mapped in October 2016 (blue in Fig. 10) and October 2020 (red in Fig. 10), respectively. Compared to 2016, these thaw slumps significantly expanded and increased for all three check regions in 2020. In particular, the area of thaw slumps had expanded 137%, 125%, and 94% for CR13, CR14, and CR15, respectively. Moreover, the number of thaw slumps had increased by 77%, 74%, and 20% for CR13, CR14, and CR15, respectively. We calculated the LOS deformation averaged for all the mapped thaw slumps, and then plotted the time series and fit curves for the three check regions (see Fig. 10). It turns out that all thaw slumps exhibit an evident trend of seasonal fluctuation, namely uplift and subsidence in the freezing and thawing season, respectively. In addition, all three check regions show an apparent trend of subsidence, and the Chumar River area shows relatively large subsidence velocity (12.0 mm/year for CR13) and periodic amplitude (7.6 mm for CR13) than the Beilu River area (CR14 and CR15, Table III). Such observation from InSAR is in accordance with the abovementioned increase in thaw slump area and number, which indicates that the Hoh Xil region was experiencing the dynamic evolution of retrogressive thaw slumps.

A. Discussion on the Methodology
This study processed 1592 Sentinel-1 SAR images using the SBAS algorithm. The biases in retrieving deformation may result from the following sources. The first source is the interferogram network configuration. Using interferograms with a relatively short temporal baseline may cause potential 'fading effect' [71]. Also, a small amount of low coherence interferograms were still reserved in order to retain a sufficiently constrained baseline network. The second source is the fixed single-phase model. The phase model used in the article did not consider the phase contribution of dry snow depth and soil moisture. Thinning of dry snow depth could counteract the actual surface uplift due to permafrost freezing in winter. As for soil moisture, the increasingly drying of thawed layers on the top of permafrost may induce spurious uplift signals to InSAR observations. These signals could also mask the actual surface subsidence during the permafrost thawing period [72]. The third source is the selection of reference pixels. It is difficult to identify completely stable reference pixels. In addition, the whole Hoh Xil region was divided into 23 processed blocks, whose reference pixels are totally different. If the selected pixels have potential uplift or subsidence, the retrieved deformation results would be biased, and the discontinuities at the block boundaries could be observed. The fourth source is from the limitation of atmospheric phase removal by spatial filtering. The fixed filter window size is not conductive to remove the local turbulent atmospheric. The low pass spatial filter may not estimate the seasonal component of the long-wavelength atmospheric signal [73]. Hence, the residual atmospheric phase would result in biases to deformation results.
In terms of the fading effect, the SBAS algorithm deploying the short temporal baseline (12-72 days in this study) interferograms assumes closed phases, however, phases may be in fact not closed due to soil and vegetation variations, thus resulting in biased results of deformation time series brought by the misclosure phase [71]. To quantify the potential influence of this fading effect on the derived velocity, we processed the datasets in a test region in Hoh Xil (34.  length) with the same reference points, but four different temporal baselines: 12-60 days, 24-72 days, 12-120 days, and 12-240 days. In order to exclude the potential bias brought by the data processing, we estimated the derived velocities before and after the atmospheric phase removal and singular value decomposition, respectively. Both results show that the velocities estimated by the shorter temporal baselines (12-60 days and 24-72 days) are larger than those by the longer temporal baselines (12-120 days and 12-240 days) (see Fig. 11). For both shorter temporal baselines (12-60 days and 24-72 days), the estimated velocities after the atmospheric phase removal and singular value decomposition are 1.5 and 2.2 times larger than the longer temporal baselines of 12-120 days and 12-240 days, respectively (see Table IV). Thus, the potential fading effect that the velocity decreases as the temporal baseline increases exists in this study. The cause of this effect is still unclear, possibly brought by soil, snow, and vegetation variations on the Tibetan Plateau. More experiments are needed to deeply investigate the physical causes of this effect.
Phase contributions from dry snow and soil moisture may bias the observed deformation time series [38], [72], [74]. We estimated such potential contribution of dry snow as [75] where ϕ snow is the snow phase term, and λ is the radar wavelength of 5.55 cm. d snow is the snow depth which was estimated from the product of daily snow depth in China based on the passive microwave brightness temperature [76], [77], [78], [79]. θ i is the radar incidence angle of 33.80 • . ϵ snow is the dielectric constant of snow, which was set to 1.298 according to Liu et al. [80]. Here, the ground surface with snow cover days longer than zero and the ground surface temperature under 0 • was considered to be covered by dry snow. We masked the snow phase result with snow cover days shorter than 12 days because it exceeds the shortest interval of the interferometric pairs. The estimated phase contribution from dry snow is shown in Fig. 12. Dry snow may only affect the periphery of the Hoh Xil region, and the estimated phases were significantly higher in the northeast and west (>1.4 rad) than those in the southeast (<0.4 rad), and its phase contribution mainly accumulates in the freezing period from November to April. Regarding the soil moisture, according to Reinosch et al. [38] and Zwieback et al. [72], its change may induce a spurious deformation, which approximately equals to λ/8 and around 10%-20% of the radar wavelength, namely 5.55-11.10 mm for the C-band Sentinel-1A interferometric pairs. The snow cover in Hoh Xil is usually shallow, patchy, and frequently with short duration (primarily less than 12 days). Most snow in Hoh Xil may have completely melted before summer. Thus, we consider that the 12-72 days interferograms could be less affected by the presence of seasonal snow cover, Fig. 13. Examples of dynamic evolution of (a) thermokarst lakes and (b) thaw slumps, as seen from the bi-temporal PlanetScope CubeSat images.
which may have already included the signal of seasonal snow accumulation and melting. Even though, it is true that some interferograms may suffer from decorrelation when spanning snow onset and melt, which were thus identified as low-quality interferograms and manually removed.
Regarding the deformation model fitting, we did not fix the seasonal period of the freeze-thaw cycle because the temperature and precipitation conditions may have interannual differences in Hoh Xil, which could influence the varying freeze-thaw state of permafrost [81]. As a result, we set the seasonal period as a free parameter although it may complicate the model fitting by making the iterative Taylor expansion necessary.
Discontinuities at the track boundaries may be produced when stitching the individual track results in the multistack analysis. Even for the LOS deformation, such discontinuities still exist due to the track-to-track variations in local incident angle. To solve this problem, we exploited the distance-weighted average to the two adjacent tracks at the overlapping regions and accordingly such track-to-track variations can be to some extent reduced.
The MT-InSAR outcome is ground deformation. It alone may insufficiently support the determination of the dynamic evolution of observed landforms. As a result, the additional interpretation of optical imageries is fundamentally needed. This is also the reason why we prepared the inventory of thermokarst lakes and thaw slumps from the optical images. For instance, as seen in Fig. 13, the dynamic evolution of thermokarst lakes and thaw slumps can be detected from the bi-temporal PlanetScope CubeSat images acquired in 2016 and 2020. Multiple small lakes were found to merge into large lake [see Fig. 13(a)]. The enlargement of thaw slumps can be also noticed, indicating that they are in the active stage of rapid development [see Fig. 13(b)]. These observations from optical imageries are essential to help determine the dynamic evolution of those thermokarst landforms.

B. Discussion on the Results
By virtue of the MT-InSAR map generated across the entire Hoh Xil region, we discovered dynamic permafrost disturbances were broadly occurring. In particular, we noticed that thermokarst landforms, namely thermokarst lakes and retrogressive thaw slumps, widely formed and dynamically evolved in Hoh Xil.
Thermokarst lakes are the typical manifestation of thawing and warming permafrost undergoing degradation, and they have been found occurring on the Tibetan Plateau [16], [21], [60]. Our result shows that thermokarst lakes appear to evolve and cluster close to large lakes in Hoh Xil. This might be due to higher water and ground ice content, which may facilitate the development of an intermediate layer in upper permafrost in those areas. The study conducted by Niu et al. [16] along the Qinghai-Tibetan engineering corridor may to some extent justify this hypothesis, who revealed that the ice content plays a more important role than the ground temperature in determining the spatial distribution of thermokarst lakes. Moreover, the warmer lake bottoms with water depth deeper than winter ice may further accelerate permafrost degradation, and subsequently from and enlarge thermokarst lakes [16]. Future field work of borehole sampling and ground-penetrating radar (GPR) is needed for a deeper understanding of this phenomenon. In addition, our study shows that dynamic permafrost degradation, which was characterized by the significantly deepened thawing depth and accelerated ground deformation, was detected close to Salt Lake. Such permafrost degradation in Salt Lake area was in accordance with the observations by Lu et al. [45], potentially resulting from the thermal alteration of the permafrost thawing-freezing cycle due to the changes in the hydrological connectivity and soil permeability.
On the Tibetan plateau, retrogressive thaw slumps are a typical type of thermokarst landform across ice-rich permafrost on slope terrain [17], [19], [20], [82]. Retrogressive thaw slumps on the Tibetan Plateau may evolve from fluvial-thermal erosion on ice-rich permafrost, and may subsequently lead to slope failures of active layer detachment at the thawing interface [83]. This study has identified two areas with the development and accumulation of thaw slumps in Hoh Xil, namely the Chumar River area and the Beilu River area. Our finding of concentrated thaw slumps in the Beilu River area is in agreement with the previous studies [17], [19], [20], [84], both of which declared a rapid evolution of retrogressive thaw slumps in the past decade. Thaw slumps developed in the Chumar River area are mainly distributed along the Tianchi Valley and the Haiding Nor River, and they have not been noticed by previous studies. Our study justifies that thaw slumps can be extensively triggered with large number in an ice-rich permafrost terrain undergoing degradation. In addition to the Tibetan Plateau, such a high concentration of thaw slumps is not rare in areas vulnerable to climate warming, especially in high-arctic environment [85], [86], [87], [88].
InSAR only measures the surface displacements and may not directly estimate the thaw depth, namely the instantaneous level at which the permafrost soil warms to 0 • C. Encouragingly, this study shows that the InSAR-derived cumulative deformation is highly correlated with the deepened thawing depth, which is a typical manifestation of permafrost disturbances. Furthermore, as seen in Fig. 4 Table V). This justified that the vegetation type at the KKXL station is alpine desert while that at the Wudaoliang station is alpine meadow (see Table I, [92]). Alpine meadow regions usually have higher soil water content than alpine desert areas [42], [93]. Moreover, the KKXL station has lower SoC content (7.59, 8.82, 7.02, 4.63, 3.03, and 2.84 g/kg) compared to the Wudaoliang station (14.29, 14.58, 12.88, 8.08, 4.82, and 3.56 g/kg) at both six depths, respectively (see Table V). This is also reasonable because higher SoC content is usually found in areas with more vegetation cover and higher water content. It turns out that slight changes in SoC content can have a large impact on the soil properties of mineral soils. The maximum seasonal amplitude of thaw subsidence is proportional to the soil water content in the active layer and accordingly vegetation type may determine the subsidence magnitude [26]. More soil water is thus needed for the alpine desert at the KKXL station to approach equivalent cumulative deformation at the Wudaoliang station. As a result, a thicker permafrost thawing depth has been reached at the KKXL station (see Fig. 4, Table I). Such a phenomenon has also been revealed by Wang et al. [32], who showed that in sub-Arctic subsidence rate and deformation pattern may vary due to the different soil water content in an active layer, which regulates the seasonal freezing-thawing cycle of permafrost at a local scale.
Previous studies have shown that there may exist lag time between air temperature and InSAR deformation over permafrost areas on the Tibetan Plateau [34], [35], [36], [38], [43]. This study further revealed that such lag time may also exist between the maximum LST and subsidence. The lag times observed in this study are relatively long, remaining 80, 120,  SOC IN THE KKXL AND WUDAOLIANG STA-TIONS AS RECORDED BYTHE NATIONAL SOIL INFORMATION GRIDS OF  CHINA [89], [90], [91] and 86 days at the KKXL, QT08, and Wudaoliang stations, respectively. Such delay of maximum subsidence is due to the fact that the downward heat diffusion from the ground surface to the bottom of the active layer fundamentally needs transferring time [34], [36]. In particular, a longer lag time between temperature and deformation indicated deeper ALT in terms of heat transfer within the soil that used in a diffusive model [34]. However, it should be noted that all these three stations are located in a relatively flat valley instead of mountain areas on the Tibetan Plateau. Cold mountain climate may expedite the freezing process of permafrost, thus shortening the thawing period and generating thinner ALT [38], [94]. As a result, the observed longer lag times in this study may also be attributed to the flat terrain in which the three stations are located. In addition, the QT08 station is situated in an area of human settlement and accordingly the subsidence may be additionally affected by anthropogenic impacts.

C. Future Direction
Future work will give priority to the following two aspects. The first direction will be devoted to retrieving ALT in Hoh Xil based on the derived InSAR outputs and soil moisture data over the region. In particular, the future model should be directed to complement the current ALT retrieval model proposed by Liu et al. [95], with the specific applicability to relatively dry soils on the Tibetan Plateau, and with a particular focus on an estimation of the water content of multiple layers with a consideration of different soil types from both in situ monitoring and remote sensing observations. The second perspective would be directed to an extension of the developed model, targeting on a complete InSAR mapping and ALT inversion over the entire Tibetan Plateau. In addition, it could be interesting to deploy the InSAR approaches of distributed scatterers for the purpose of obtaining higher density of scatterers in low-coherence areas [96], [97], and to involve multitemporal change detection from multisensor platforms [98] for automated and rapid discrimination of land process variations (such as a wider mapping of thaw slumps) in Hoh Xil. Our InSAR outcomes will be openly shared with public access, and additional modification, improvement, and interpretation of these data from users are anticipated.

VI. CONCLUSION
In this study, a total of 1592 Sentinel-1 SAR images spanning from 2015 to 2020 were manually processed based on the SBAS algorithm, for generating the MT-InSAR deformation outcomes across the entire Hoh Xil region on the Tibetan Plateau hinterland. Two types of products were derived particularly for the regional permafrost investigation, namely: 1) the linear velocity of ground deformation and 2) the periodic amplitude of seasonal deformation. The linear velocity of ground deformation is capable to estimate the rate of permafrost thawing subsidence while the periodic amplitude of seasonal deformation can be adapted to determine the seasonality of permafrost thaw subsidence and freeze uplift. The produced MT-InSAR outcomes show that the Hoh Xil region was broadly experiencing dynamic permafrost disturbances, which were primarily recognized as the formation of thermokarst lakes, retrogressive thaw slumps, and deepened thawing depth. The thawing permafrost with both a linear subsidence rate higher than 2 mm/year and a periodic amplitude over 2 mm was mainly identified in areas of flat or gentle slopes. The InSAR cumulative deformation is highly correlated with the permafrost thawing depth, potentially induced by diverse soil compositions and vegetation types that may cause the variance in soil water content. In addition, the significant lag times were identified between the seasonal oscillation of InSAR deformation and the LST, demonstrating the need for the transferring time for the downward heat diffusion that initiated from the ground surface to the bottom of the active layer. Thermokarst landforms, specifically thermokarst lakes and retrogressive thaw slumps, widely formed and dynamically evolved. In particular, widespread thawing permafrost characterized by the accumulation of thermokarst lakes was likely to occur in areas adjacent to large lakes in Hoh Xil. Such dynamic evolution of thermokarst lakes, observed with a maximum linear velocity of 10.4 mm/year and periodic amplitude of 9.5 mm, may be considered as a joint consequence of permafrost degradation and accelerated partial drainage brought by the higher content of water and ground ice. Moreover, massive retrogressive thaw slumps were triggered and spread as a result of permafrost degradation, and they were found clustered on hillslopes in the Chumar River area and the Beilu River area, with a maximum accumulated subsidence of 141 and 108 mm during 2015-2020, respectively. This study demonstrated the usefulness of remote sensing and particularly MT-InSAR in broad investigation of typical permafrost disturbances of thermokarst lakes, retrogressive thaw slumps, and deepened thawing depth in the unpopulated hinterland of the Tibetan Plateau, where the extreme harsh environment prohibited extensive in situ observations and field work throughout the whole region. We anticipate this MT-InSAR outcome may be useful for further retrieving progression of ALT in Hoh Xil, and for an in-depth understanding of the linkage between permafrost processes and climate change on the Tibetan Plateau. APPENDIX See Fig. 14.