Monitoring the Thaw Slump-Derived Thermokarst in the Qinghai-Tibet Plateau Using Satellite SAR Interferometry

Thaw slumps are well-developed within a 10 km wide zone along the Qinghai-Tibet engineering corridor, especially along the Qinghai-Tibet highway and railway. Previous studies have focused on thaw slump instability such as its origin development, headwall retrogression rate, failure scale, and thermal regime, yet the intrinsic dynamic process of surface movement is relatively less known. In this study, we used InSAR based on the L-band ALOS PALSAR images acquired from January 2007 to October 2010 to investigate the distribution of thaw-induced slope failures containing retrogressive thaw slumps and active layer detachment failures along the Qinghai-Tibet highway (QTH). Our InSAR analysis reveals that the maximum annual average sedimentation rates are even up to -35mm·yr in the slope direction to the K3035 thaw slump, and the K3035W active layer detachment failure developed on the west side of K3035. The distribution, failure extent, and stability of the slope failures obtained by our InSAR analysis all agree well with the field investigations. Our study illustrates that InSAR is an effective tool for studying the distribution and processes of the thaw slump-derived thermokarst and provides useful references for evaluating permafrost degradation in response to climate warming and external disturbance on the Qinghai-Tibet plateau.


Introduction
Permafrost degradation is an important indicator of global warming, which has drawn considerable attention in recent decades. Global warming over the last century and further warming in the following long period of time are an indisputable fact according to the Intergovernmental Panel on Climate Change Fourth Assessment Report (IPCC AR4, Acquired from IPCC Report: https://www.ipcc.ch/reports/) [1]. Qinghai-Tibet plateau (QTP) permafrost is characterized by low-middle latitude, high altitude, and high ground temperature compared to polar permafrost which seems to be more susceptible and vulnerable to climate change and external disturbance [2]. Under climatic and humaninduced disturbance, plateau permafrost has been undergoing strong warming and degradation. Thermokarst landforms are important indicators of permafrost degradation, which result from the thawing of ice-rich permafrost or melting of massive ground ice followed by surface subsidence results in irregular, depressed landforms. Generally, thermokarst landforms can be divided into dozens of types such as retrogression thaw slumps, active layer detachment slides, thermal erosion gully, depressions, ice grotto, ice tower, and thermokarst lakes [3]. Thermokarst landforms are common in permafrost regions of the QTP. Temporal dynamics and spatial variability of thermokarst landforms have a strong impact on hydrologic conditions, ecosystems, ground-air carbon exchange, and human infrastructures. However, thermokarst landforms on the QTP, especially nonlake ones, are seldom studied and poorly understood [4].
Thaw slump is one of the most typical and damaging nonlake thermokarst landscape in the QTP. On the one hand, it has a great impact on the cold region environment, resulting in partial melting of permafrost even wear. Meanwhile, the lateral thermal erosion of the thaw slump will lead to deepen the permafrost thaw depth. On the other hand, it has serious damage to infrastructure projects in the QTP, especially for the large Qinghai-Tibet railway (QTR) and QTH, which caused uneven subsidence near the embankment and the embankment toe slip tread. Therefore, in order to avoid potential geohazards, it is urgent to investigate the distributions and dynamics processes of the thaw slump.
Previous investigation on thaw slump mostly reply on in situ observations, including spirit leveling, GPS, and boreholes [5], and only focused on instability such as its origin development, headwall retrogression rate, failure scale, and thermal regime in the local area, yet the long-term intrinsic dynamic process of surface movement and the eventual stabilization in the regional scale are relatively less known. Liu et al. have successfully applied Interferometric Synthetic Aperture Radar (InSAR) to map localized thermokarst subsidence trends at a thermokarst landform from 2006 to 2010, which are induced by constructing a pipeline access road [4]. So InSAR is an effective tool for mapping and studying thermokarst processes at both local and regional spatial scales.
The objectives are as follows: (1) to validate and investigate the distributions and active level of the thaw-induced slope failures with field observations; (2) to quantify spatialtemporal variability in surface deformation at K3035 and K3035W thaw slumps from 2007 to 2010; and (3) to extend time series InSAR application to map and study thaw slumps in the QTP, which have not been well investigated in previous studies.

Study Area.
Considering available ground-based measurements, our study area is located in a 10 km wide zone along the QTH in the Beiluhe basin region; it extends from 92°56 ′ to 93°02 ′ E and from 34°56 ′ to 35°02 ′ N. In a regional setting, this area is mainly underlain by continuous permafrost [6]. The Beiluhe basin region is relatively flat; most elevations distributed between 4500 and 4700 m. The average air temperature is about −5.2°С, where the maximum low and high air temperature reach approximately −38°С and 23°С, respectively, and the frozen period is from September to April. The mean annual ground temperatures are from -2.0 to -0.5°С, and the thickness of the active layer is from 1.5 to 2.0 m [7]. Compared with the scarce precipitation (approximately 290.0 mm·yr −1 ), evaporation is much higher (approximately 1316.9 mm·yr −1 ), resulting in arid conditions characterized by a negative water budget [5,8,9].
Thermokarst landforms, especially thaw-induced slope failures, are well-developed in the Beiluhe basin region. Thaw-induced slope failures are usually initiated by a reduction in shear strength of ice-rich sediments in sloping terrain [10]. Slope failures have two major categories; one is the retrogressive thaw slump and another is the thaw-induced shallow landslides. Retrogressive thaw slumps are usually triggered by the thaw of exposed massive ice, which leads to the development of a steep frozen head scarp that recedes with additional thaw of ice-rich material, while thawinduced landslides commonly occur when high pore water pressures develop at the permafrost table as a result of excess water release from the thaw of ice-rich permafrost [7,11].
According to previous field investigations, there are 25 slope failures within our study area, wherein a part of slope failures was in a stable state [5]. In this study, we focus on the typical K3035 slope failures (location: 34°59′N, 92°58 ′ E), which are located at K3035 mileage of the QTH (see Figure 1). It extends from 92°56 ′ to 93°2 ′ E and from 34°56 ′ to 35°02 ′ N. 2

Journal of Sensors
The formation of the K3035 thaw slump was triggered by digging for pavement material of the QTH in 1992. The digging slope toe blocked the normal thermal equilibrium and resulted in underground ice directly exposed at the surface and permafrost melting. In addition, the underground ice provides a good sliding surface. The saturated clay soil losing freeze strength can easily form slump between the ice surface and the clay soil underwater lubrication. Currently, it has formed a 117 m long and 86 m wide thaw slumping zone. Based on field observation and aerial photographs taken in 2006 and 2010, the thaw slump was still developing [5].

SBAS and Datasets.
In this section, we review the fundamental concepts behind the SBAS algorithm. An exhaustive description of the SBAS algorithm can be referenced in Berardino et al. [12]. This technique is based on the combination of multiple differential interferograms, which is a postprocessing algorithm to generate the date by date displacements and the mean deformation velocity [13].
Let us start by assuming a set of single-look complex SAR images acquired by the same satellite over the same area; the imaging time was t 0 , … , t N−1 ; K interferograms are generated with small temporal and spatial baselines to suppress spatial decorrelation and the residual phase. The differential phase for a generic coherent pixel in the interferogram, generated by combining SAR acquisitions at times, can be expressed in In formula (1), where v is the deformation rate, λ is the transmitted radar wavelength, ϕ t B and ϕ t A are the phases corresponding to the times t B and t A , respectively, and d t B and d t A are the radar line of sight projections of the cumulative deformation at t B and t A with respect to the first scene t 0 assumed as a reference. Δh is the DEM residuals, which is proportional to the perpendicular baseline for each interferogram B ⊥j , θ, is radar satellite look angle. d APS t B and d APS t A are the possible atmospheric delay signal corresponding to the times t B and t A , respectively. Other unknown decorrelation effects and noise sources are included in n j . The following step is to implement the SBAS inversion, retrieving the estimate of the displacement rate and the residual topography. Usually, three main different fitting models can be considered: linear, quadratic, and cubic. The most robust inversion model is the linear model; the other models need high coherent interferograms to provide reliable results. Typically, without any a priori knowledge on the displacement behavior, the best choice is to choose the linear model. In this study, we chose the linear model to inverse the deformation parameters and the height residuals. Then the atmospheric artifacts are mitigated through high-pass filtering in the temporal domain and low-pass filtering in the spatial domain [14,15]. In the end, SBAS uses the singular value decomposition (SVD) approach based on a minimum norm criterion of the deformation rate to derive time series deformation measurements.
A total of 21 L-band ALOS PALSAR images acquired with ascending orbit from January 2007 to October 2010 were used in this study. The SAR data were obtained from the Japan Aerospace Exploration Agency (JAXA), collected with a nominal radar look angle of 34°. Only the HH channel from two fine modes was used for interferometric processing after the FBS data were doubly downsampled in range direction, as listed in Table 1. Compared to other C or X band satellite images, L band ALOS PALSAR penetrates better, which formed high-quality interferograms. Besides the SAR datasets, an accurate reference digital elevation model (DEM), in an ellipsoidal reference system, is needed before going on with the differential processing. We used a DEM to remove the topographic contribution to the interferograms and in geocoding the results (transforming Range-Doppler coordinates into Universal Transverse Mercator map geometry system). If a better DEM is not available, the free shuttle radar topography mission (SRTM) from the United States Geological Survey (USGS) can be used. In this study, we used the 3-arcsecond (~90 m) SRTM, which has a posting 90 m and a vertical accuracy of 16 m. In data processing, we oversampled the DEM to a posting of 25 m to match the interferograms. To improve the signal-to-noise ratio (SNR), all the interferograms were multilooked by the factors of 1 and 5 along the range and azimuth directions, respectively. The ground range pixel dimensions of the resulting products are therefore about 16 × 16 m in the range and azimuth directions.

Results
SBAS analysis using temporal baseline smaller than 10 months and spatial baseline smaller than 3500 m was based on 21 PALSAR images from January 2007 to October 2010. This analysis identified 15 master and 78 candidate SAR image pairs. After manual inspection, 11 low-quality interferograms (large phase discontinuities and low coherence) were identified and discarded, which could significantly decrease the precision in the subsequent parameter estimation (e.g., displacements, DEM residuals, and atmospheric artifacts). 67 interferograms were used for further deformation estimation and time series analysis. The final interferogram network shows that all the pairs are in one subset (see Figure 2).
As no ground control points are available in this area, a high coherent point in a stable region was chosen as a reference and marked by the black triangle in Figure 3(a). Since the thaw slump deformation mainly occurs along the slope direction, the SBAS-derived deformation was transformed from the line of sight (LOS) into the slope direction. The SBAS-derived mean surface displacement rate along the  Journal of Sensors slope over K3035 thaw slump site is illustrated in Figure 3(a). The negative sign of the deformation rate stands for an increasing distance with time away from the satellite (e.g., subsidence), whereas the positive sign of the deformation rate stands for a decreasing distance with time towards the satellite (e.g., uplift). It can be observed in Figure 3 The total time series deformation of the study area is presented in Figure 4. The spatial-temporal changes of the thaw slump in the period from January 2007 to October 2010 can be seen. We also observe that the areal extent and the mean rate of the thaw slump were increasing for the period from January 2007 to October 2010. In addition, the time series of four feature points (A, B, C, and D in Figure 3) was further analyzed and the results are illustrated in Figures 5(a)-5(d), respectively, and simultaneously, four trend lines were added to reflect the deformation tendency of these feature points. From Figure 5(a), the cumulative deformation of point A has exceeded 140 mm in the period from January 2007 to October 2010. That is to say, the thaw slump continued to develop. From Figure 5(b), although there is a seasonal fluctuation (uplifting from October to April and sinking from April to October), it is shown that there is a sinking tendency of point B and the cumulative deformation has exceeded 120 mm by combing the time series curve and trend line. As we know, point B is located in another side of the QTH, which meant that in this location there is instability.
Wang et al. have monitored a thaw slump of the QTP and analyzed the solifluction creep characteristics. The results show that, due to the reduction of the buffering effects of the active layer on heat transfer influenced by thaw slump, the ablation of underlying ice-rich permafrost is increasing and also brings about the increase of the thaw settlement [16]. According to preliminary inference, this phenomenon may be caused by the heat conduction effect of the thaw slump, but the real fact needs further verification. Point C is located in the central area of another thaw slump we found in this study. The location of thaw slump was verified by an optical image from Google Earth image. Obviously, Figure 5(c) shows the same seasonal fluctuation as B, but the cumulative deformation of point C, whose tendency is subsiding, has exceeded 100 mm according to its trend line and time series analysis. In addition, it is clearly found that this newfound thaw slump was developing fast from

Discussion
Slope failures are common in permafrost regions of the QTP, especially during human activities and construction of largescale projects, e.g., QTH, QTR. In this study, we investigated the distribution of thaw-induced slope failures and analyzed the temporal dynamics and spatial variability of the K3035 and K3035W thaw slumps along the Qinghai-Tibet highway (QTH) in the Beiluhe section based on SBAS method. To validate the reliability of the results, we compare the results from the previous study and analyze the causes of thaw slumping.
A total of 33 slope failures were identified by InSAR in the study area ( Figure 6). Of these, twenty-five slope failures whether the location distributions or the level of activity were consistent with the field observations; eight slope failures need further verification. Eleven black circles represent inactive slope failures; fourteen red crosshairs mark the still active landslides; eight red triangles are slope failures which need to be verified; the yellow rectangle represents the position of the K3035 and K3035W thaw slumps. In addition, we focus on analyzing the temporal dynamics and spatial variability of the K3035 thaw-induced slope failures, which is the retrogressive thaw slump developed at Kilometer Post K3035. Previous studies of the K3035 thaw slump focused on its headwall retrogression rate, thermal regime, and stability. The thaw slump was initialed by the digging for pavement material of the QTH in 1992. The thaw of exposed massive ice leads to the development of a steep frozen head scarp that recedes with additional thaw of ice-rich material. Niu et al. have done years of field observation (ground temperature measured and boreholes) on the K3035 thaw slump. They found that during the initial years from 1992 to 2000, the retrogression rate of the thaw slump was fast, with an average retrogression rate of the headwall about 11 mm·yr −1 ; from 2000 to 2002, the headwall retrogression rate was about 4.5 m·yr −1 ; the headwall retrogression rate has slowed down in recent years after 2002 [5,7,17].

Conclusions
This paper has demonstrated the capability of the small baseline InSAR technique for monitoring active thaw slump processes and offers new insights on permafrost degradation on slopes in response to external disturbance. It is suggested that this technique could be also used on a continuous basis to monitor other kinds of thermokarst process (active layer detachment slides, thermal erosion gully, etc.) and hazard along the Qinghai-Tibet engineering corridor.
Our InSAR analysis focused on the dynamics of the thaw slump, which is a complement to previous studies by Niu et al. Overall, we draw the following conclusions: (1) It is proved that InSAR is an effective tool for detecting the location distribution and processes of the thaw slump-derived thermokarst   Journal of Sensors (3) Quantify spatial-temporal changes of the K3035 and K3035W thaw slumps during the period from 1992 to 2000 using SBAS (4) It is found that there was relative large subsidence onto the other side of the QTH at K3035 thaw slump, which may be resulted from a mixture of mudflows or water flows induced from thaw slump through road culvert

Data Availability
The DEM data used to support the findings of this study may be released upon application to the account (http://www.eorc. jaxa.jp/ALOS/en/aw3d30/data/index.htm), who can be contacted at email. The ALOS-1 data used to support the findings of this study were supplied by JAXA under license and so cannot be made freely available. Requests for access to these data should be made to JAXA. Authors obtain ALOS-1 data by applying quota of scientific research data.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.