*3.1. Data Sources*

To investigate the activity and kinematics of the twin landslides, multisource highresolution 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. 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◦, heading angle: −10.1◦, range resolution: 4.66 m, azimuth resolution: 3.16 m) and thirteen ALOS-2 (Path: 147, Frame: 750, incidence angle: 36.3◦, heading angle: −10.4◦, 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).

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


#### *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 Δ*A*area [6]:

$$
\Delta A\_{\text{rate}} = \frac{A\_2 - A\_1}{t\_2 - t\_1} \tag{1}
$$

where *A*<sup>1</sup> and *A*<sup>2</sup> are the total area (m2) of landslide in different timeframes and *t*<sup>1</sup> and *t*<sup>2</sup> (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–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. 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.


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

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*los) can be projected into the downslope velocities (*V*ds) with the following equation [47,53]:

$$V\_{\rm ds} = \frac{V\_{\rm los}}{\sin(a\_{\rm aspect} - \beta)\sin\theta\_{\rm inc}\cos a\_{\rm slope} + \cos\theta\_{\rm inc}\sin a\_{\rm slope}}\tag{2}$$

where *α*aspect and *α*slope are the aspect and slope angles, respectively, which can be calculated using the SRTM DEM data; *β* is the flight direction of the SAR satellite; and *θ*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.
