Activity and Kinematics of Two Adjacent Freeze–Thaw-Related Landslides Revealed by Multisource Remote Sensing of Qilian Mountain

Abstract

The increase in temperatures and changing precipitation patterns resulting from climate change are accelerating the occurrence and development of landslides in cold regions, especially in permafrost environments. Although the boundary regions between permafrost and seasonally frozen ground are very sensitive to climate warming, slope failures and their kinematics remain barely characterized or understood in these regions. Here, we apply multisource remote sensing and field investigation to study the activity and kinematics of two adjacent landslides (hereafter referred to as “twin landslides”) along the Datong River in the Qilian Mountains of the Qinghai-Tibet Plateau. After failure, there is no obvious change in the area corresponding to the twin landslides. Based on InSAR measurements derived from ALOS PALSAR-1 and -2, we observe significant downslope movements of up to 15 mm/day within the twin landslides and up to 5 mm/day in their surrounding slopes. We show that the downslope movements exhibit distinct seasonality; during the late thaw and early freeze season, a mean velocity of about 4 mm/day is observed, while during the late freeze and early thaw season the downslope velocity is nearly inactive. The pronounced seasonality of downslope movements during both pre- and post-failure stages suggest that the occurrence and development of the twin landslide are strongly influenced by freeze–thaw processes. Based on meteorological data, we infer that the occurrence of twin landslides are related to extensive precipitation and warm winters. Based on risk assessment, InSAR measurements, and field investigation, we infer that new slope failure or collapse may occur in the near future, which will probably block the Datong River and cause catastrophic disasters. Our study provides new insight into the failure mechanisms of slopes at the boundaries of permafrost and seasonally frozen ground.

1. Introduction

A landslide is the downslope movement of soil, rock, and debris under the action of gravity and the landform that results from such movement [1]. The factors triggering slope failure mainly include rainfall, earthquakes, fluvial erosion, excavation, and construction activities [1]. Slope failure occurs frequently in rainy and mountainous areas, often resulting in casualties and property damage [1,2]. Slope failure at high latitudes and high altitudes in cold regions has attracted significant attention for decades. Climate warming and the consequent varied freeze–thaw states in cold regions may accelerate the occurrence and development of landslides [3,4]. Thaw-related landslide activities (such as retrogressive thaw slumps and active layer detachment) are increasing extensively in the permafrost regions of both the Arctic [5,6,7,8] and the Qinghai-Tibet Plateau (QTP) [9,10,11], a phenomenon that is mainly related to extremely warm summers and extensive summer precipitation. Even in very cold permafrost regions, extreme warming events can cause slope failure due to the thawing of ice-rich permafrost [12]. However, slope failures on the boundaries of permafrost and seasonally frozen ground remain little characterized or understood.

Optical and radar remote sensing techniques are becoming one of the most important tools for investigating the activity and kinematics of landslides, especially in rural and mountainous regions [13,14]. The occurrence and development of landslides can be obtained from optical remotely sensed images [15,16,17,18,19]. Rapid and abrupt slope movements can be quantified according to the differences in DEM from two successive orthorectified images [20,21,22] or from single-pass InSAR observations [23,24]. However, optical images are prone to the effects of atmospheric conditions and revisit time, which impede the quantification of landslide activities. In addition, optical-based slope movement is insensitive to gradual slope movement or creep, which is often a sign of slope stability. Interferometric synthetic aperture radar (InSAR) has been increasingly explored and successfully applied to identify and monitor gradual slope deformation with centimeter-to-millimeter accuracy [13,25,26,27,28]. The advantage of InSAR is that it is less affected by cloudy weather and works regardless of whether day or night conditions. Recently, the development of unmanned aerial vehicle (UAV) cameras allows the generation of very-high-resolution optical images and digital surface models (DSMs). The combination of InSAR and UAV data has advanced the understanding of crack development, landslide evolution, and susceptibility [29,30,31,32,33,34].

In this study, we combine multisource optical and radar remote sensing images to investigate the activity and kinematics of two adjacent landslides (hereafter referred to as “twin landslides”) located in the boundary regions of permafrost and seasonally frozen ground on Qilian Mountain in the QTP. High-resolution optical satellite and UAV-based images are used to investigate the activity of the twin landslides. InSAR measurements are used to map and quantify the gradual slope movement before and after the failure of the twin landslides. The triggering mechanisms of the twin landslides and their potential risks are analyzed.

2. Study Area

In this study, we investigate two adjacent landslides occurring on the slope of the middle and lower reaches of Datong River on the southeast slope of Qilian Mountain. For convenience, we named the two adjacent landslides as twin landslides and refer to them as “QLDT01” and “QLDT02” throughout the paper (Figure 1C). The study area is located between Tuolai Mountain in the north and Datong Mountain in the south, both of which belong to branches of the Qilian Mountains. The Wari Gaqu River rises from the Tuolai Mountains and ultimately joins the Datong River. A number of brooks from the Datong Mountains = discharge into the Datong River (Figure 1A). Therefore, the study area is a place where many rivers converge. The Tuolaishan and Datongshan faults are distributed in the northwest and southwest of the study area, respectively. The altitude in the study area ranges from 3400 to 3600, which is in the lower boundary of permafrost and seasonal frozen ground on the QTP. The permafrost distribution map is shown in Figure 1B [35].

Qilian Mountain is dominated by a continental alpine semi-humid mountain climate, which is characterized by long, cold, and dry winters and short, cool, and moist summers. Annual average air and ground surface temperatures are approximately −2.4 and 2.9 ${}^{\circ }$C, respectively [36]. The coldest January averages below −11 ${}^{\circ }$C, and the warmest July averages below 25 ${}^{\circ }$C; most of the Qilian Mountains are below 0 ${}^{\circ }$C from December to March, and the highest temperatures range from 4 to 15 ${}^{\circ }$C from April to October. Precipitation ranges from 300 to 500 mm annually and is mainly concentrated during May to September in the form of rainfall [37,38].

3. Data and Methods

3.1. Data Sources

To investigate the activity and kinematics of the twin landslides, multisource high-resolution optical and radar data were acquired during 2009–2020. Google images with a spatial resolution of 0.65 m were acquired using Map Tile Downloader (version number: release 2.3, developed by Centmap Co., LTD., located in Hefei, Anhui Province, China). The details of the remotely sensed images are listed in

Table 1. Summary of the remotely sensed dataset used in this study. The acquisition dates of ALOS PALSAR-1/2 can be found in Table 2.

Data	Resolution (m)	Date (YYYYMM)	Number of Scenes
Google	0.65	200910, 201006, 201712, 202007	4
Gaofen-2	0.8	201512, 201811, 202008	3
ALOS PALSAR-1	7	2008–2010	6
ALOS PALSAR-2	7	2015–2020	13
UAV	0.05	202104	1

. Gaofen-2 is a Chinese high-resolution optical satellite that was launched in August 2014 and carries two panchromatic and multispectral charge-coupled device camera sensors. We used panchromatic images with a spatial resolution of 0.8–1 m, which were obtained from the China Centre for Resources Satellite Data and Application (http://www.cresda.com/CN/, accessed on 10 October 2021). Based on the selection of images free of cloud and snow cover, four Google and three Gaofen-2 images obtained during 2019–2020 were used.

The L-band Phased Array type L-band Synthetic Aperture Radar (PALSAR) sensor is mounted on the Advanced Land Observation Satellite (ALOS)-1/2. Six ALOS-1 (Path: 477, Frame: 750, incidence angle: 38.7${}^{\circ }$, heading angle: −10.1${}^{\circ }$, range resolution: 4.66 m, azimuth resolution: 3.16 m) and thirteen ALOS-2 (Path: 147, Frame: 750, incidence angle: 36.3${}^{\circ }$, heading angle: −10.4${}^{\circ }$, range resolution: 4.29 m, azimuth resolution: 3.77 m) SAR images were chosen to map and quantify the ground movements of the twin landslides. We expected better performance from the L-band PALSAR data with a wavelength of about 24 cm. This is because L-band electromagnetic waves can penetrate deeper into snow and vegetated surfaces [39], leading to higher interferometric coherence [40].

To investigate potential unstable zones, we conducted a UAV survey and used structure-from-motion/multiview stereo photogrammetry to map the twin landslides and their surroundings on 29 April 2021. We used the DJI Phantom 4 RTK flying platform with an altitude of 100 m above the ground surface. The heading and lateral overlap of flying were 85% and 60%, respectively. We obtained the very-high-resolution (VHR) images and a DEM with a resolution of around 5 cm/pixel using Agisoft PhotoScan software. The uncertainty of the relative positions was estimated to be around 2–3 cm.

The temperature and precipitation data in the Qilian Mountains from 2000 to 2019 obtained by the National Meteorological Station of China were used to describe the climatic conditions (http://www.cma.gov.cn/, accessed on 10 November 2021).

3.2. Mapping of Twin Landslides from Optical Remote Sensing

Google, Gaofen-2, and PALSAR-1/2 satellite images were used to determine the occurrence and development of the twin landslides. The Gaofen-2 images were geometrically corrected using ENVI5.3 software. Then, the one arc-second Shuttle Radar Topography Mission (SRTM) DEM product was used for image orthorectification. Due to destruction of the integrity of the original stratum, landslide features such as changed vegetation and soil collapse can be identified from high-resolution optical images [16,41]. The boundaries of the twin landslides and adjacent shorelines were outlined based on visual inspection by three experienced researchers. To further evaluate the change characteristics, we estimated the rate of landslide areal growth $\Delta A_{\mathrm{area}}$ [6]:

$$\Delta A_{\mathrm{rate}}=\frac{A_2-A_1}{t_2-t_1}$$

(1)

where $A_1$ and $A_2$ are the total area (m${}^2$) of landslide in different timeframes and $t_1$ and $t_2$ (year) are the corresponding time points.

VHR optical images have proven useful in identifying landslide features such as small cracks or ground discontinuities [42,43]. In this study, to map the potential unstable zones surrounding the twin landslides, we identified cracks through visual inspection of the VHR UAV optical images. As the spatial resolution of UAV images is 5 cm, cracks with width larger than 5 cm were very likely to be identified.

3.3. InSAR for Ground Deformation Monitoring

The InSAR technique detects ground movements by comparing the phase differences between SAR images acquired from slightly different positions at different times [44]. Differential InSAR (DInSAR) and multi-temporal InSAR (MTInSAR) have frequently been used to measure slope movements in both permafrost and nonpermafrost regions [27,28,45,46,47]. As there are very limited descending PALSAR-1/2 images, all the archived ascending PALSAR-1/2 images available that covered our study area were examined for ground deformation monitoring. The interferometric coherence decreases rapidly in the thaw season (May–October) and is slower in the freeze season (November–next April). To mitigate the decorrelation impact, we selected image pairs with temporal spans of less than 150 days and perpendicular baselines shorter than 500 m. Considering the accuracy of SRTM DEM and the maximum perpendicular baseline (591 m in our case), we estimated that the residual topographic phase would be about 0.8 radians, corresponding to 1.5 cm in the InSAR measurements. Relying on the interferometric coherence and phase quality, only six ALOS PALSAR-1 images and thirteen ALOS PALSAR-2 images, taken during 2009–2010 and 2015–2020, respectively, were selected. This causes severe disconnection between SAR images and does not allow the use of MTInSAR approaches such as small-baseline subset InSAR [48].

We calculated ground movement using the DInSAR technique, which was conducted using the commercial GAMMA software [49]. We constructed three and nine interferograms for PALSAR-1 and -2, respectively. The range and azimuth look numbers were 2 and 5 for PALSAR-1 and 2 and 4 for PALSAR-2, generating ground pixels of approximately 15 m × 15 m. The one arc-second SRTM DEM product was used to remove the topographic phase of each interferogram. The temporal and perpendicular baselines are presented in

Table 2. The interferogram pairs from ALOS PALSAR-1/2 and their temporal and perpendicular spatial baselines.

ALOS PALSAR-1				ALOS PALSAR-2
ID	Master–Slave (YYYYMMDD)	Time Span (Days)	${\mathrm{B}}_{\perp }$(m)	ID	Master–Slave (YYYYMMDD)	Time Span (Days)	${\mathrm{B}}_{\perp }$(m)
1	20090630–20090815	46	−32	1	20151009–20151218	70	200
2	20091231–20100215	46	491	2	20171201–20180209	70	−98
3	20100703–20100818	46	116	3	20180209–20180601	112	−44
				4	20180601–20180727	56	−6
				5	20181116–20190125	70	15
				6	20190125–20190531	126	−60
				7	20190531–20190726	56	78
				8	20190726–20190906	42	−10
				9	20200306–20200529	84	169

. We applied a power spectrum adaptive filter to mitigate the phase noise and mask out decorrelation areas with a coherence threshold of 0.6 [50]. We unwrapped all the interferograms using the minimum cost flow approach [51]. To compare the deformation between PALSAR-1 and -2, a local reference point with high coherence nearby the twin landslides was selected for calibration of the unwrapped phase. Tropospheric artifacts may contaminate the ground deformation in mountainous regions. As our study area was very small, we mitigated tropospheric artifacts by fitting the topographic-related components [52]. Residual atmospheric and orbital errors were mitigated using a linear deramping approach.

We calculated the light-of-sight (LOS) movement from each interferogram. By dividing the time interval between the interferogram pairs, we calculated the deformation velocities along the LOS direction. Assuming the slopes move purely along the downslope direction, the InSAR-estimated LOS velocities ($V_{\mathrm{los}}$) can be projected into the downslope velocities ($V_{\mathrm{ds}}$) with the following equation [47,53]:

$$V_{\mathrm{ds}}=\frac{V_{\mathrm{los}}}{sin({\alpha }_{\mathrm{aspect}}-\beta )sin{\theta }_{\mathrm{inc}}cos{\alpha }_{\mathrm{slope}}+cos{\theta }_{\mathrm{inc}}sin{\alpha }_{\mathrm{slope}}}$$

(2)

where ${\alpha }_{\mathrm{aspect}}$ and ${\alpha }_{\mathrm{slope}}$ are the aspect and slope angles, respectively, which can be calculated using the SRTM DEM data; $\beta$ is the flight direction of the SAR satellite; and ${\theta }_{\mathrm{inc}}$ is the local incidence angle, which can be calculated using the SAR geometry and SRTM DEM data. To reduce the noise in the calculation of slope, aspect, and local incidence angles, we applied a Gaussian filter with a 7 × 7 window (around 200 m) to the SRTM DEM.

3.4. Climatic Factors

Air temperature and precipitation data from 2000 to 2019 were used to analyze their impacts on the evolution of the twin landslides. We calculated four temperature indicators: mean annual air temperature (MAAT), thawing index, warming days, and average temperature in the coldest month of the year. To account for warming days, we calculated the number of days with a daily temperature higher than 10 ${}^{\circ }$C. The thawing index $\mathrm{TI}$ is the cumulative number of degree days above 0 ${}^{\circ }$C for a given thaw season, which can be calculated by [54]:

$$TI=\sum _{i=1}^{N_T}{T}_i,T_i>0$$

(3)

where $T_i$ is the daily temperature on day i and $N_T$ is the number of days in a year with a temperature greater than 0 ${}^{\circ }$C.

We calculated the annual total precipitation, precipitation intensity, extreme precipitation, and the number of consecutive drought days in a year. The precipitation intensity is the ratio between the total precipitation and the duration of precipitation days, and represents the average amount of precipitation in a certain duration. Daily precipitation of between 10 and 25 mm is defined as moderate rainfall by the World Meteorological Organization. However, because the annual precipitation is about 450 mm, we consider a daily precipitation of higher than 15 mm to be extreme precipitation in our study. The consecutive drought days is the number of days without precipitation.

3.5. Risk Assessment

We evaluated the potential risks related to the twin landslides and their surroundings in the same slopes. A landslide dam forms when a landslide reaches the bottom of a valley and causes partial or complete blockage of a river [55]. The sudden collapse of landslide dams and the rapid release of water storage poses a great risk of flooding downstream [56]. The dimensionless blockage index (DBI) has been developed for the prediction of potential risks of a landslide dam by linking the stability of a landslide dam to three geomorphic parameters [57]. The dam volume ${\mathrm{V}}_{\mathrm{d}}$ controls the dam height ${\mathrm{H}}_{\mathrm{d}}$, and is considered as the main stabilizing factor. The watershed area ${\mathrm{A}}_{\mathrm{b}}$ indirectly controls the channel flow and flow power, and is the main factor influencing dam instability. The dam height is an important variable for evaluating the stability of landslide dams against overtopping and pipeline failure. Thus, the DBI can be expressed as [57]:

$$\mathrm{DBI}=log\left(\frac{{\mathrm{A}}_{\mathrm{b}}\times {\mathrm{H}}_{\mathrm{d}}}{{\mathrm{V}}_{\mathrm{d}}}\right)$$

(4)

As only QLDT01 has caused the formation of a landslide dam, we calculated the DBI only for QLDT01. The dam height was obtained from UAV-based DEM data. The volume of the landslide dam was calculated from high-resolution UAV-based DEM using the cut-and-fill volume tool in the Global Mapper software. The hydrological analysis tool was used to calculate the catchment area from UAV-based DEM in ArcGIS (version number: release 10.7, developed by Environmental Systems Research Institute, Inc., located in RedLands, CA, USA).

4. Results

4.1. Spatiotemporal Variations of the Twin Landslides

The occurrence and development of the twin landslides are shown below (QLDT01 in Figure 2 and QLDT02 in Figure 3). Based on visual interpretation of optical and PALSAR-1/2 SAR backscatter images, we infer that QLDT01 occurred between October and December of 2009, whereas QLDT02 occurred sometime between October and November of 2015. The failure of QLDT01 caused the mud and rubble to slide into the river channel and almost blocked the Datong River (Figure 2B). Crack features can be found in the headwall regions for both QLDT01 and QLDT02 as far back as 2009 (Figure 2A and Figure 3A). As there are no high-resolution images in our study before 2009, we cannot precisely determine the exact initialization time of these cracks. Compared to 2009, the cracks in QLDT02 had significantly enlarged during 2009–2015 before its failure. The development of QLDT01 is slow, and its slide into the Datong River almost stopped during 2010–2018. Contrary to QLDT01, the headwall of QLDT02 continued to slowly retreat during 2009–2015. The mass of QLDT02 slid along the northwest side of the headwall region and caused the formation of a dammed lake at the foot of the slope.

The total area of the QLDT01 slide is about 76.5 × 10${}^3$ m${}^2$ following the slope failure in 2009. The landslide body slid into and dammed the Datong River. QLDT01 slowly expanded at an areal growth rate of 0.5 × 10${}^3$ m${}^2$ during 2011–2018 (Figure 2C). The total area for slope failure of QLDT02 is about 131 × 10${}^3$ m${}^2$, which is about double that of QLDT01 (Figure 3C). A small dammed lake has formed at the toe of QLDT02. The areal growth rate of QLDT01 is 10.7 × 10${}^3$ m${}^2$ during 2016–2018, which has slowed to 5.5 × 10${}^3$ m${}^2$ during 2018–2020.

The slope failure of QLDT01 completely dammed the Datong River and rerouted its flow (Figure 4). The width of the Datong River beneath QLDT01 was 66 m before the slope failure in 2009. The landslide body slid into the river and reached to about 4 m beyond the northern bank when the slope collapsed. The river quickly expanded towards the northern bank, whereas the river’s width changed to 16 m in 2010. Under continuous fluvial erosion, the northern bank expanded northward by about 30 m during 2010–2017, whereas the river’s width changed to 48 m in 2017. In other words, the average bank erosion rate was about 4 m/year during 2010–2017.

4.2. InSAR-Derived Downslope Movement of the Twin Landslides

We derive the downslope movement of twin landslides before and after their failure from PALSAR-1/2 InSAR measurements during 2009–2010 and 2015–2020 (Figure 5). Significant downslope movement is observed within the twin landslides, whereas the maximum displacement rate reaches up to 15 mm/day. We observe strong displacement of up to 5 mm/day outside the twin landslides. To evaluate the potential risks related to the twin landslides and their surroundings, we outline one polygon adjacent to the QLDT02 based on the phenomenological features from UAV images (potential risk zone (PRZ) in Figure 5L).

During the summer before the failure of QLDT01, significant downslope velocities up to 15 mm/day are observed in the boundary and central regions of the landslide body, whereas the mean value is about 4 mm/day (Figure 5A). Five years later, after the failure of the slope, the mean downslope velocities are smaller than 0.5 mm/day in both the summer and winter seasons during 2015–2020 (Figure 5D–F).

We observe that for QLDT02 the mean downslope velocities are about 1.6 mm/day with a maximum value of 5 mm/day from July 2009 to August 2010, i.e., the periods just before and after the failure of QLTD01 (Figure 5A–C). A distinct scarp can be observed in the high-resolution optical image at the head of the landslide body (Figure 3A), which may cause severe InSAR decorrelation and result in no measurements in these regions. The failure of QLDT02 occurred during October and December of 2015; however, there are no valid InSAR measurements due to this severe decorrelation. We observe that significant downslope velocities with mean values of about 2.3 mm/day are pronounced in QLDT02 during July–February of 2016–2020 after slope failure (Figure 5E,H,K). On the contrary, QLDT02 is inactive during March–June (Figure 5F,G,I,J,L).

5. Discussion

5.1. Triggering Mechanisms

5.1.1. Precipitation

In general, changing precipitation patterns increase subsurface saturation and pore pressure, which increase the likelihood of slope failure [58]. Extensive or extreme precipitation and rapid snow/ice melt are therefore likely to increase the frequency and magnitude of landslides [59,60]. The precipitation data near the twin landslides show fluctuating upward trends during 2000–2019 (Figure 6A). The annual precipitation is 451 mm in 2009 and 448 mm in 2015, which is significantly higher than the mean annual precipitation since 2000. The annual precipitation in the preceding year is higher than in the year of slope failure.

To evaluate the impacts of precipitation events on the occurrence of the twin landslides, we graphed the daily precipitation before slope failure (Figure 7). An extreme precipitation event of 35.8 mm is recorded in August 2009. The accumulated precipitation was about 193 mm during August–September 2019 (Figure 7). The number of consecutive drought days is among the lowest during 2000–2009 (Figure 6C). Extensive rainwater may increase the pore water pressure and reduce the shear strength in weak soil layers. Thus, we infer that the occurrence of QLDT01 may have been primarily triggered by extensive precipitation. In 2015, there was no extreme precipitation event such as that in 2009, however, the annual precipitation was higher than the 20-year average (Figure 6A). The number of consecutive drought days was below the average (Figure 6C). Therefore, we presume that increased precipitation is likely one of the triggering factors of landslide QLDT02.

5.1.2. Freeze–Thaw Processes

Climate warming and disturbance may have strong impacts on slope stability in cold environments [3,4]. In permafrost areas, rocks are glued together by ice filling their cracks and crevices. Freeze–thaw processes are characterized by variability in subsurface temperature and moisture content, which results in substantial fluctuations of shear strength (cohesion and friction angle) and drives landslide initiation [61]. The transition from perennially frozen to seasonally frozen ground accelerates the effect of freeze–thaw processes on both bedrock and unconsolidated material [62]. As the air temperature increases, the warming and thawing of permafrost may weaken rock faces and the inherent stability of permafrost, leading to slope failure [63]. In mountain permafrost regions, e.g., the European Alps, Canada, and the Tibetan plateau, researchers have recorded an increasing tendency of landslide activities due to the warming climate [9,64,65].

The MAAT shows an obvious warming trend during 2000–2019 (Figure 8A). The MAAT in 2009 and 2015 is 2.32 and 2.36 ${}^{\circ }$C, respectively, which is about 0.3 ${}^{\circ }$C above the 20-year average (Figure 8A). The MAAT values in the preceding years (2008 and 2014) are about 0.7 ${}^{\circ }$C lower than in the failure years. The strong fluctuations in air temperature may amplify freeze–thaw processes and thus affect slope stability. In 2009, the warming days (the number of days with air temperature above 10 ${}^{\circ }$C), the thawing index, and the average temperature in the coldest month were all above their 20-year averages (Figure 8B–D). This suggests that the warming events in 2009 might have been one of the triggers of QLDT01 failure. On the contrary, the warming days and thawing index in 2015 are lower than their averages. However, the mean temperature of the coldest month in 2015 is about 0.7 ${}^{\circ }$C) above the 20-year average, which suggests a warm winter. Warm winters may slow down freezing processes, allow the soil water to remain in an unfrozen state for a longer time, thereby increasing the risk of slope failure.

The activities and kinematic patterns in the twin landslides and their surroundings have been derived from InSAR measurements. We observe that the average downslope velocities in QLDT01 and QLDT02 exhibit distinct seasonality (Figure 9). During the early thawing periods from May to early July, the slopes are in an inactive state. In this stage, soil thawing is shallow and does not reach the sliding surface, resulting in limited downslope movement. During the late thawing and early freezing periods from late July to the next January the slopes are in an active state, with average downslope velocities up to 4 mm/day. In the late thawing stage, the sliding surface is thawed, which results in significant downslope movement. Despite the shallow soil being frozen during the early freezing season, downslope movement remains significant, as the sliding surface is in a thawed state. During the early freezing period from February to April the slopes become inactive, as the sliding surface is in a frozen state.

The seasonality of significant downslope movement during both the pre- and post-failure stages suggest that the occurrence and development of the twin landslides were strongly influenced by freeze–thaw processes. The seasonal pattern is different from the seasonal deformation corresponding to freeze uplift and thaw subsidence due to ice–water phase change constraints in the active layer [66,67]. To put this work in a spatial context, we compare our study with several freeze–thaw-related slope instability studies on the QTP. Meng et al. and Hao et al. [34,45] observed deformation velocity up to 100 mm/year with a linear trend assumption using the multi-temporal InSAR technique on an earthflow in Yushu, QTP. Dini et al. [68] characterized different magnitudes of LOS deformation over different types of slope instability in the eastern Himalayas. Hu et al. [53] found similar seasonal patterns of downslope velocity up to about 3 mm/day during the active stage in several rock glaciers in the East Kunlun Mountains. The less pronounced downslope velocity may be primarily related to the kinematic behaviors of rock glaciers.

5.1.3. Other Triggering Factors

The slope failure of QLDT01 may be partially attributed to fluvial erosion at the slope toe and its geomorphological characteristics. QLDT01 is situated at the confluence of several rivers. The north bank of the Datong River facing landslide QLDT01 is the Wari Gaqu River, which flows into the Datong River. This results in high runoff flow of the Datong River, which is usually accompanied by transverse expansion when it is scoured downward along the river. Under continuous erosion by river flow or streams, the slope toe becomes too steep to hold itself, consequently resulting in slope failure [69]. Moreover, QLDT01 has a slope of about 17.5 degrees on average, with a slope height of about 66 m, which make it prone to slope failure.

Shaking from earthquakes may be a direct triggering factor of the QLDTL02 failure. According to data from the China Earthquake Networks Center (http://www.ceic.ac.cn, accessed on 10 November 2021), a Mw 5.2 earthquake occurred on 23 November 2015, with a focal depth of 10 km and a direct distance of 24 km from QLDT02 (Figure 1A). The Tuolaishan fault is the seismogenic fault of this earthquake according to the Qinghai Earthquake Administration, China (www.qhdzj.gov.cn, accessed on 10 November 2021). Earthquakes increase the occurrence of landslides due to ground shaking, liquefaction of susceptible sediments, and swelling of soil materials caused by shaking, which allows water to seep in rapidly. In addition, earthquakes can alter friction at the base of landslides, thus accelerating their movement over several days or weeks [70,71].

5.2. Hazard Analysis

We evaluate the stability of the landslide dam of QLDT01 based on the DBI calculation. According to the DBI criterion proposed by Ermini and Casagli [57], the state of landslide dam can be categorized as a stable domain (DBI < 2.75), an uncertain domain (2.75 < DBI < 3.08), and an unstable domain (DBI > 3.08). The height of the landslide dam ranges from 0 m at the toe of the landslide to about 40 m at the south bank of the Datong River. The landslide dam volume and catchment area are 69 × 10${}^3$ m${}^3$ and 4.4 × 10${}^3$ m${}^2$, respectively. Relying on different dam heights, the calculated DBI ranges from 2.55 to 2.97, with an average of 2.78 (Figure 10). We find that the toe (${\mathrm{H}}_{\mathrm{d}}$ < 3.5 m) and top (${\mathrm{H}}_{\mathrm{d}}$ > 29 m) of the slope in QLDTL01 can be considered as a stable domain, as their DBI is lower than 2.75. The landslide dam is in the uncertain domain in the middle of the slope (3.5 m < ${\mathrm{H}}_{\mathrm{d}}$ < 29 m), which accounts for 70% of the entire slope. Thus, we infer that the QLDT01 is at risk of further slope collapse.

While only a portion of the slope (QLDT02) has collapsed, we evaluate the stability of the noncollapsed regions of the slope and the potential risks. Two long cracks (about 300–400 m) could be observed as of 2009. While one crack (QLDT01) collapsed in 2009, only a small portion of another crack developed into a landslide (QLDT02) in 2015. Based on the high-resolution UAV DEM (Figure 11B), we find that the slope height varies significantly and the slope gradient is large, providing geomorphological conditions for slope creep. In addition, many new cracks are found in the noncollapsed regions of the slope (Figure 11A), suggesting the occurrence of strong internal movement. Moreover, continuous InSAR-derived downslope movements are observed, further confirming the instability of the noncollapsed slope (PRZ in Figure 5 and Figure 9). The volume of the PRZ region is about 12 × 10${}^5$${\mathrm{m}}^3$, which is 1.6 times larger than that of landslide QLDT01. In addition, there are temporary houses in the area for locals to graze animals. A potential slope failure may completely block the Datong River and cause a catastrophic disaster.

6. Conclusions

We have documented the spatiotemporal evolution of two adjacent landslides on the southeast slope of Qilian Mountain during their pre-failure and post-failure stages from 2008 to 2020 by integrating multisource optical and radar remote sensing techniques. The main conclusions are as follows:

1.

The occurrence and development of the twin landslides and the adjacent river bank have been determined using high-resolution optical and radar images.

2.

Significant downslope movements are observed within the twin landslides (up to 15 mm/day) and their surrounding slopes (up to 5 mm/day). The downslope movement exhibits distinct seasonality. During February-early July, the downslope velocity is nearly inactive; during late July-next January, mean velocity of about 4 mm/day is observed. The seasonality of downslope movement during both pre- and post-failure stages suggests that the occurrence and development of the twin landslide are strongly influenced by freeze–thaw processes.

3.

Combining data on precipitation, air temperature, and InSAR-based deformation history, we infer that the occurrence of the twin landslides is mainly related to extensive precipitation, warm winters, and freeze–thaw processes.

4.

From the UAV-based geomorphological features, InSAR-derived downslope movements, and dimensionless blockage index, we infer that new collapse and slope failure events may occur within the twin landslides and their surroundings, which may completely block the Datong River and cause a catastrophic disaster.

Our study demonstrates the capability of multisource high-resolution remote sensing techniques to monitor landslide activities in cold regions. As the impacts of climate warming becoming more extensive, freeze–thaw-related slope instability in climate-sensitive regions (the boundary regions of permafrost and seasonally frozen ground, in this case) should be afforded greater attention.