*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 ◦C. The thawing index TI is the cumulative number of degree days above 0 ◦C for a given thaw season, which can be calculated by [54]:

$$TI = \sum\_{i=1}^{Nr} T\_{i\nu} \qquad T\_i > 0 \tag{3}$$

where *Ti* is the daily temperature on day *i* and *NT* is the number of days in a year with a temperature greater than 0 ◦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 Vd controls the dam height Hd, and is considered as the main stabilizing factor. The watershed area Ab 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]:

$$\text{DBI} = \log(\frac{\text{A}\_{\text{b}} \times \text{H}\_{\text{d}}}{\text{V}\_{\text{d}}}) \tag{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-andfill 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 (Figures 2A and 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.

**Figure 2.** Temporal variations of landslide QLDT01: (**A**,**B**) Google satellite images of landslide boundary changes in 2009 and 2010 and (**C**) landslide evolution in 2010, 2017, and 2020.

The total area of the QLDT01 slide is about 76.5 × 103 <sup>m</sup><sup>2</sup> 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 × <sup>10</sup><sup>3</sup> m2 during 2011–2018 (Figure 2C). The total area for slope failure of QLDT02 is about 131 × 103 m2, 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 × 103 m2 during 2016–2018, which has slowed to 5.5 × <sup>10</sup><sup>3</sup> m2 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.

**Figure 3.** Temporal variations of landslide QLDT02: (**A**,**B**) landslide in Google satellite and Gaofen-2 satellite images from 2009 and 2015; the yellow arrow is the direction of movement of the landslide. (**C**) Landslide characteristics recorded by Gaofen-2 satellite images in 2015, 2018, and 2020. A dammed lake is formed at the toe of the slope, as shown in (**D**).

**Figure 4.** River bank changes of Datong River beneath the twin landslides in 2009, 2010, and 2017. In 2009, the thinnest section of the Datong River was about 66 m. In 2010, the southern river bank expanded northward due to the collapse of Landslide QLDT01, and the width was about 16 m. In 2017, the north bank continued to expand northward, and the thinnest section of the Datong River was 48 m.

#### *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).

**Figure 5.** The downslope velocity is derived from line-of-sight (LOS) deformation using Equation (2). The background map is the shaded relief map derived from UAV DEM. The twin landslides (QLDT01 and QLDT02) and the potential risk zone (PRZ) are marked by red polygons in the bottom right-hand corner of the subfigures. The positive values refer to the movement in the downslope direction. The red triangle denotes the location of the reference point.

#### **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.

**Figure 6.** Variations of precipitation during 2000–2019: (**A**) variations in annual precipitation, (**B**) precipitation intensity, (**C**) number of consecutive drought days in a year.

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.

**Figure 7.** Daily and cumulative precipitation during August–December, 2019. A maximum daily precipitation of 35.8 mm was recorded in August 2009.

#### 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 ◦C, respectively, which is about 0.3 ◦C above the 20-year average (Figure 8A). The MAAT values in the preceding years (2008 and 2014) are about 0.7 ◦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 ◦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 ◦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.

**Figure 8.** Variations in air temperature during 2000–2019: (**A**) mean annual air temperature (MAAT), (**B**) warming days with an annual average temperature greater than 10 ◦C, (**C**) thawing index, (**D**) average temperature of the coldest month.

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.

**Figure 9.** Temporal variations for 2000–2019 in mean daily air temperature (MDAT), mean daily precipitation (MDP), and average downslope velocities of QLDT01, QLDT02, and potential risk zone (PRZ) (Figure 5L). The squares present the average downslope velocities, whereas the corresponding lines show the start and end date of the SAR image pair. The error bars denote the standard deviation of downslope velocities within the red polygons. The bottom panel plots the MDAT and MDP.

The seasonality of significant downslope movement during both the pre- and postfailure 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 × 103 m3 and 4.4 × <sup>10</sup><sup>3</sup> <sup>m</sup>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 (Hd < 3.5 m) and top (Hd > 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 < Hd < 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 Figures 5 and 9). The volume of the PRZ region is about <sup>12</sup> × 105 m3, 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.

**Figure 10.** The dimensionless blockage index (DBI) diagram of landslide QLDT01. When the landslide dam is at the height of 3.5 to 29 m, the mean value of DBI is 2.87, which is in the uncertain domain (UD). When the height of landslide dam is lower than 3.5 m and higher than 29 m, the average DBI is 2.66, which is in the stable domain (SD).

#### **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:


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.


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.

**Author Contributions:** Conceptualization, T.W., J.C. and J.Z.; methodology, J.C. and J.Z.; validation, J.C., J.Z., J.H., X.Z., P.L. and L.Z.; investigation, T.W., J.C., X.W. and X.Z.; field work, J.C., J.Z., J.H., X.Z., X.M. and P.L.; writing—original draft preparation, J.C., J.Z. and T.W.; writing—review and editing, J.C., J.Z., T.W., J.H., X.W., X.M., X.Z.; visualization, J.Z. and J.C.; supervision, T.W.; funding acquisition, J.C. and T.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the National Natural Science Foundation of China (Grant No. 42001072 and 41771076), the fellowship of China Postdoctoral Science Foundation (2021T140702 and 2021M693374), the State Key Laboratory of Cryospheric Science (Grant No. sklcs-zz-2022), the National Cryosphere Desert Data Center Program (Grant No. E01Z790208), the Natural Science Foundation of Gansu Province of China (Grant No. 21JR7RA242), and the CAS "Special Research Assistant program" (Jie Chen).

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflicts of interest.
