A Large Old Landslide in Sichuan Province, China: Surface Displacement Monitoring and Potential Instability Assessment

: Using advanced Differential Interferometric Synthetic Aperture Radar (InSAR) with small baseline subsets (SBAS) and Permanent Scatter Interferometry (PSI) techniques and C-band Sentinel-1A data, this research monitored the surface displacement of a large old landslide at Xuecheng town, Lixian County, Sichuan Province, China. Based on the MassMov2D model, the effect of the dynamic process and deposit thickness of the potentially unstable rock mass (deformation rate < − 70 mm/year) on this landslide body were numerically simulated. Combined with terrain data and images generated by an Unmanned Aerial Vehicle (UAV), the driving factors of this old landslide were analyzed. The InSAR results show that the motion rate in the middle part of the landslide body is the largest, with a range of − 55 to − 80 mm/year on average, whereas those of the upper part and toe area were small, with a range of − 5 to − 20 mm/year. Our research suggests that there is a correlation between the LOS (line of sight) deformation rate and rainfall. In rainy seasons, particularly from May to July, the deformation rate is relatively high. In addition, the analysis suggests that SBAS can provide smoother displacement time series, even in areas with vegetation and the steepest sectors of the landslide. The simulation results show that the unstable rock mass may collapse and form a barrier dam with a maximum thickness of about 16 m at the Zagunao river in the future. This study demon-strates that combining temporal UAV measurements and InSAR techniques from Sentinel-1A SAR data allows early recognition and deformation monitoring of old landslide reactivation in complex mountainous areas. In addition, the information provided by InSAR can increase understanding of the deformation process of old landslides in this area, which would enhance urban safety and assist in disaster mitigation.


Introduction
There is a large old landslide with an estimated volume of about 3.5 × 10 7 m 3 at Xuecheng town, Lixian County, Sichuan Province, China (103.3206 • E, 31.5673 • N) (hereafter called the XC landslide). It is located on the left bank of the Zagunao river. A ground survey and remote sensing interpretation show that the back part of this landslide has an obvious armchair shape and its deposit mass is well preserved. It can be inferred that this landslide has experienced two slides based on its morphology: the first is the slope failure, which was presumably triggered by a historical earthquake; and the second occurred on the existing landslide mass. Subsequently, as a result of river erosion, small-scale collapses have taken place on the front. Thus, there is significant concern about whether this old landslide will be reactivated, thus leading to a geologic hazard, and the influences on its potential instability.
In recent years, new remote sensing techniques, such as Interferometric Synthetic Aperture Radar (InSAR) and Unmanned Aerial Vehicle (UAV) photogrammetry have been widely used in investigations of landslides [1][2][3][4][5][6][7]. InSAR technique has the advantages of being available at all times and in all weather conditions [8]. In particular, due to the application of SAR satellites of higher resolution and shorter revisit times, and InSAR data processing techniques, multiple temporal InSAR (MT-InSAR) methods have been developed, such as Permanent Scatter Interferometry (PSI) [9][10][11], small baseline subsets (SBAS) [12,13], and Distributed Scatter (DS) approaches [14,15]. Many case studies have demonstrated that PSI and SBAS technologies can facilitate large-scale landslide monitoring, and are suitable for measurement of mass movement at regional and slope scales [15][16][17][18].
Numerical methods can help understand the movement process of large-scale landslides and predict the runout of the potentially unstable mass [19,20]. The commonly used numerical simulation software packages include DAN3D, Massflow, and Mass-Mov2D [21,22]. Among these, MassMov2D is a new open-source landslide model with high operation efficiency, which can be employed to reveal the failure process and kinematic characteristics of landsliding, and aid in monitoring the motion of an unstable rock mass [22].
InSAR has becomes an effective technique for deformation monitoring and identification of landslides, and allows measurement of millimetric ground displacement [4][5][6]. In addition, UAVs, also termed remotely piloted aircraft (RPA) technology, can obtain very high-resolution images, thus increasing the vertical resolution and accuracy of such final products as the digital surface model (DSM) and orthophotos, which are important supplements to conventional investigations of landslides [1,7]. The current study made a further investigation of the XC landslide using UAV and InSAR techniques. Based on PSI and SBAS methods, the deformation rate of this old landslide was measured based on the C-band Sentinel-1A data, and two time-series results were compared. In terms of the InSAR results, we then used the MassMov2D model to simulate the possible motion process and accumulation characteristics of the unstable rock of this old landslide body. Finally, the driving factors of the instability of this landslide and hazard prediction of potential failure were investigated. This study allows us to better understand the relationship between the deformation and potential triggering factors. Furthermore, it also has guiding significance for the assessment of the reactivation of other old landslides in the eastern margin of the Tibetan plateau.

Study Area
Lixian County is located in the northwest of Sichuan Province, 200 km from Chengdu ( Figure 1). This area is close to the Longmen Shan fault zone, between the Songpan-Ganzi block in the west and Sichuan Basin in the east. Associated with complex tectonic movement, four earthquakes with magnitude greater than M S 7.0 have occurred around this area during the past 100 years. Due to the combination of rugged topography, earthquakes, heavy rainfall, and human activity, this region is prone to landslides. Previous work shows that this area is among the worst hit by geological disasters caused by the 2008 Wenchuan Mw7.9 earthquake [23]. Following this shock, numerous old landslide deposits and unstable slopes have experienced large continuous deformation [24,25]. The study area has a large topographic relief in which the highest eleva m and the lowest is 1400 m, and an average elevation of 2700 m (Figure 2a). T annual temperature in the valley is between 10 and 13 °C, and the annual ra tween 400 and 700 mm. The rainy season is from May to September, and accou thirds of the annual precipitation. The monthly average rainfall from May t largest, and is more than 70 mm [27]. Most streams pertain to the Zagunao R is a major tributary of the Minjiang River. The slopes on both sides of the Zag are steep with sparse vegetation and strong erosion. As a result, geologic haz debris flows and landslides occur frequently [27]. The lithology of the study includes Silurian gray phyllite with limestone of the Maoqun Formation (Smx Devonian Carbonaceous phyllite or phyllite with quartzite of the Weiguan (Dwg1, Dwg2), Carboniferous and Permian basalt and conglomeratic limestone Upper Triassic feldspathic quartz sandstone of the Zhunwei Formation (T3zh ternary deposits are present in the bottom of the Zagunao River valley, includ proluvial, eluvial, and old landslide masses. The study area has a large topographic relief in which the highest elevation is 5900 m and the lowest is 1400 m, and an average elevation of 2700 m (Figure 2a). The average annual temperature in the valley is between 10 and 13 • C, and the annual rainfall is between 400 and 700 mm. The rainy season is from May to September, and accounts for two-thirds of the annual precipitation. The monthly average rainfall from May to July is the largest, and is more than 70 mm [27]. Most streams pertain to the Zagunao River, which is a major tributary of the Minjiang River. The slopes on both sides of the Zagunao River are steep with sparse vegetation and strong erosion. As a result, geologic hazards such as debris flows and landslides occur frequently [27]. The lithology of the study area mainly includes Silurian gray phyllite with limestone of the Maoqun Formation (Smx 4 and Smx 5 ), Devonian Carbonaceous phyllite or phyllite with quartzite of the Weiguan Formation (Dwg 1 , Dwg 2 ), Carboniferous and Permian basalt and conglomeratic limestone (C + P), and Upper Triassic feldspathic quartz sandstone of the Zhunwei Formation (T 3 zh). Thin Quaternary deposits are present in the bottom of the Zagunao River valley, including alluvial, proluvial, eluvial, and old landslide masses.

Landslide Characteristics
The XC landslide is located on the left bank of the Zagunao River, where the valley is deep and narrow, and exhibits a "V" shape. The landslide has a "U" shape in plane view, in which most of the material was displaced from the upper part and deposited in the lower part ( Figure 3). Based on its topographical features, it is a typical large-scale old landslide with an area of approximately 1.4 km 2 , and is about 1900 m long and 850 m wide. The main sliding direction is 130°. The main scarp has an armchair-like shape with a height of about 260 m. The elevations of the main scarp and the toe are 2750 and 1500 m, respectively, with an elevation difference of about 1250 m. The slope angle averages 35°, and reaches 50~60° at some places. According to field investigations and remote sensing interpretations, two sliding events can be distinguished, and the main scarp and lateral scarps of two slides can be clearly identified ( Figure 3). The first slide is the slope failure, and the second event occurred on the deposit due to the first event. The gullies and surface water flows can be clearly seen in the whole landslide mass. Previous research [28] indicates that the hydrological paths are closely related to the micro-topography evolution and stability of landslides. Based on UAV-based DSM data with 0.1 m resolution, the hydrological paths of the landslide mass were extracted using GRASSGIS software [29]. Figure 3 shows that the areas with the densest hydrological paths are in the middle part of the old landslide.

Landslide Characteristics
The XC landslide is located on the left bank of the Zagunao River, where the valley is deep and narrow, and exhibits a "V" shape. The landslide has a "U" shape in plane view, in which most of the material was displaced from the upper part and deposited in the lower part ( Figure 3). Based on its topographical features, it is a typical large-scale old landslide with an area of approximately 1.4 km 2 , and is about 1900 m long and 850 m wide. The main sliding direction is 130 • . The main scarp has an armchair-like shape with a height of about 260 m. The elevations of the main scarp and the toe are 2750 and 1500 m, respectively, with an elevation difference of about 1250 m. The slope angle averages 35 • , and reaches 50~60 • at some places. According to field investigations and remote sensing interpretations, two sliding events can be distinguished, and the main scarp and lateral scarps of two slides can be clearly identified ( Figure 3). The first slide is the slope failure, and the second event occurred on the deposit due to the first event. The gullies and surface water flows can be clearly seen in the whole landslide mass. Previous research [28] indicates that the hydrological paths are closely related to the micro-topography evolution and stability of landslides. Based on UAV-based DSM data with 0.1 m resolution, the hydrological paths of the landslide mass were extracted using GRASSGIS software [29]. Figure 3 shows that the areas with the densest hydrological paths are in the middle part of the old landslide.

SAR Datasets
Sentinel-1A, a space mission of the Europe Space Agency (ESA), was launched on 3 April 2014 [30]. Compared to other SAR missions, due to its revisit time of 12 days and short baselines (i.e., orbital separations), Sentinel-1A is promising in the application of InSAR [31,32]. The SAR data used in this study were taken from Sentinel-1 satellites, and were acquired with the Interferometric Wide (IW) swath mode. Its pixel resolution is about 3.5 m in the range direction and 14 m in the azimuth direction. A total of 63 ascending scenes from 8 April 2018 to 21 April 2020 were acquired and processed with an incidence angle of 39.45°. POD Precise Orbit Ephemerides published by the Sentinel-1A satellite were used to remove systematic errors caused by orbit errors (https://qc.senti-nel1.eo.esa.int/, accessed on10 April 2021). In addition, an ALOS PALSAR DEM of 12.5 m resolution (https://search.asf.alaska.edu/, accessed on10 April 2021) was used as the reference elevation dataset for removing the topographic phase and geocoding the InSAR results.

UAV Data
The UAV used for photogrammetry was a FEIMA Intelligent Aerial Survey System D2000. This is a small long-endurance quadrotor UAV system, which can carry optical to microwave sensors to obtain multi-source data [33]. UAV data sets used in this work were collected on 23, 24, and 27 September 2020. The ground resolution was set to 5 cm, the images were designed with an overlap of 80%, and the flight altitude was about 160 m. UAV manager software developed by FEIMA Company was used to process images (http://www.feimarobotics.com/zhcn/productDetailManager, accessed on10 April 2021). This software can generate georeferenced spatial data, such as 3D point clouds, DSMs, and digital orthophoto models (DOM) [33]. Finally, we obtained the DSM and DOM data with a resolution of 0.1 m.

SAR Datasets
Sentinel-1A, a space mission of the Europe Space Agency (ESA), was launched on 3 April 2014 [30]. Compared to other SAR missions, due to its revisit time of 12 days and short baselines (i.e., orbital separations), Sentinel-1A is promising in the application of InSAR [31,32]. The SAR data used in this study were taken from Sentinel-1 satellites, and were acquired with the Interferometric Wide (IW) swath mode. Its pixel resolution is about 3.5 m in the range direction and 14 m in the azimuth direction. A total of 63 ascending scenes from 8 April 2018 to 21 April 2020 were acquired and processed with an incidence angle of 39.45 • . POD Precise Orbit Ephemerides published by the Sentinel-1A satellite were used to remove systematic errors caused by orbit errors (https://qc.sentinel1.eo.esa.int/, accessed on 10 April 2021). In addition, an ALOS PALSAR DEM of 12.5 m resolution (https://search.asf.alaska.edu/, accessed on 10 April 2021) was used as the reference elevation dataset for removing the topographic phase and geocoding the InSAR results.

UAV Data
The UAV used for photogrammetry was a FEIMA Intelligent Aerial Survey System D2000. This is a small long-endurance quadrotor UAV system, which can carry optical to microwave sensors to obtain multi-source data [33]. UAV data sets used in this work were collected on 23, 24, and 27 September 2020. The ground resolution was set to 5 cm, the images were designed with an overlap of 80%, and the flight altitude was about 160 m. UAV manager software developed by FEIMA Company was used to process images (http://www.feimarobotics.com/zhcn/productDetailManager, accessed on 10 April 2021). This software can generate georeferenced spatial data, such as 3D point clouds, DSMs, and digital orthophoto models (DOM) [33]. Finally, we obtained the DSM and DOM data with a resolution of 0.1 m.

MT-InSAR Processing
Time-series analysis based on multiple SAR acquisitions can be performed using a set of differential interferograms generated from the co-registered image stack. Multireferencing and master-referencing are the two approaches used to build a stack of interferograms, and these two techniques were developed to carry out DInSAR analysis, including PSI, and SBAS [34].
The PSI method was first proposed by Ferretti et al. [35], and bounces satellite radar signals of naturally occurring permanent scattering points, such as man-made objects (buildings, roads) and bare rocks [36]. This technique is based on a set of pixels that have a stable phase over long periods of time, and are used to estimate the motion rate [37]. The SBAS technique is a different time series method based on multi-master images, and proposed by Berardino et al. [38]. This method overcomes the drawbacks of the poor coherence of some interferograms caused by the use of one master image. In addition, because only interferograms with small baselines are selected for the time series analysis, the deformation results are denser and more reliable. However, in the SBAS algorithm, the interferograms are subjected to multi-look processing to further reduce decorrelation noise. Thus, spatial filtering is not needed in the SBAS method, which avoids increasing the noise contribution of non-dominant scatterers by coarsening of the resolution [34]. In the PSI method, all interferograms can be created with respect to a single master, which allows reduction of the noise contribution of the master image prior to phase unwrapping because it is present in all interferograms. Therefore, the two approaches can be considered to be complementary in the usual case in which a dataset contains pixels with a range of scattering characteristics [34].
In this study, for the PSI method, a master image on 12 May 2019 was chosen based on the minimization of the perpendicular and temporal baselines (Figure 4a). Low and high band pass filters with window sizes of 1200 m and 365 days were then applied to remove the spatial and temporal distributions of the atmospheric variations. Interferogram generation and inversion were performed using the SARscape software [39]. For the SBAS method, the maximum normal baseline and temporal baseline were set to 2% and 60 days, respectively, and interferogram pairs were generated (Figure 4b). All unwrapped interferograms obtained were analyzed individually to inspect large unwrapping errors before the time series generation. To ensure low coherence of the interferograms, we removed these interferogram pairs [40]. The phase unwrapping was conducted using Delaunay Minimum Cost Flow (Delaunay MCF) [41] and the Goldstein filter was used to perform interferometric filtering [42]. After phase unwrapping, 36 ground control points (GCPs) were selected to correct the unwrapped phase. The selection criteria can be found in detail in previous work [43]. Preliminary displacements were estimated by a linear model that is commonly used. The 12.5 m DEM was used to remove the topographic phase and geocoding the SAR results in the SARscape software. The result was geolocated to the WGS84 map geometry.

MassMov2D Model
MassMov2D is a numerical model used to simulate the expansion (runout) and deposition of mass over a complex topography by approximating the sliding mass to a homogeneous one-phase fluid [44]. The model describes the mass movement as a two-dimensional flux taking advantage of the shallow water equations [45]. This formula is derived from the general Navier-Stokes equations under the hypothesis that the vertical components of velocity and pressure are negligible with respect to the horizontal components, and that the vertical pressure profile can be considered to be almost hydrostatic [46]. The code is freely available and implemented in GRASSGIS [22].
The input data includes the sliding surface topography (topographic data without the sliding body) and the initial sliding mass thickness. Other input parameters that characterize the mass material are the Chèzy (turbulent) coefficient (m·s −2 ), internal friction angle of the sliding mass (degree), and the fluid flow rate (m·s −1 ) [22,47]. The turbulent term summarizes all velocity-dependent factors of flow resistance and the density of the

MassMov2D Model
MassMov2D is a numerical model used to simulate the expansion (runout) and deposition of mass over a complex topography by approximating the sliding mass to a homogeneous one-phase fluid [44]. The model describes the mass movement as a twodimensional flux taking advantage of the shallow water equations [45]. This formula is derived from the general Navier-Stokes equations under the hypothesis that the vertical components of velocity and pressure are negligible with respect to the horizontal components, and that the vertical pressure profile can be considered to be almost hydrostatic [46]. The code is freely available and implemented in GRASSGIS [22].
The input data includes the sliding surface topography (topographic data without the sliding body) and the initial sliding mass thickness. Other input parameters that characterize the mass material are the Chèzy (turbulent) coefficient (m·s −2 ), internal friction angle of the sliding mass (degree), and the fluid flow rate (m·s −1 ) [22,47]. The turbulent term summarizes all velocity-dependent factors of flow resistance and the density of the debris. The fluid rate defines the transaction velocity of the sliding mass when it changes from a solid state to a fluid state.
The work flow of this study is shown in Figure 5. The basic data comprised 63 ascending C-band Sentinel-1A data, the 12.5 m DEM, UAV data, and monthly average rainfall and lithology data. First, we analyzed the surface displacement of the landslide with SBAS and PSI techniques. Then, we selected potential unstable rock mass (motion rate < −70 mm/year) to simulate the kinematic and accumulation characteristics based on the MassMov2D software. Finally, by combining the monthly average rainfall and lithology distribution, we investigated the potential driving factors of the landslide.
Remote Sens. 2021, 13, x FOR PEER REVIEW debris. The fluid rate defines the transaction velocity of the sliding mass when it from a solid state to a fluid state.
The work flow of this study is shown in Figure 5. The basic data compris cending C-band Sentinel-1A data, the 12.5 m DEM, UAV data, and monthly aver fall and lithology data. First, we analyzed the surface displacement of the lands SBAS and PSI techniques. Then, we selected potential unstable rock mass (moti −70 mm/year) to simulate the kinematic and accumulation characteristics base MassMov2D software. Finally, by combining the monthly average rainfall and distribution, we investigated the potential driving factors of the landslide.

InSAR Time Series
The PS points refer to naturally occurring permanent scattering points, such made structures. The buildings, roads, and bare rocks in the study area are di along both side of the gully with sparse vegetation, which increase the likelihood ing a persistent scatterer within a resolution cell ( Figure 3).
As shown in Figure 6a, 1164 PS points were identified in the whole slope ar PS points are mainly distributed in the middle and lower parts, whereas those in t part are relatively sparse due to vegetation cover. Based on the InSAR result a

InSAR Time Series
The PS points refer to naturally occurring permanent scattering points, such as manmade structures. The buildings, roads, and bare rocks in the study area are distributed along both side of the gully with sparse vegetation, which increase the likelihood of finding a persistent scatterer within a resolution cell ( Figure 3).
As shown in Figure 6a, 1164 PS points were identified in the whole slope area. These PS points are mainly distributed in the middle and lower parts, whereas those in the upper Remote Sens. 2021, 13, 2552 9 of 21 part are relatively sparse due to vegetation cover. Based on the InSAR result and UAV data, the area of deformation rate <−10 mm/year was delineated as the presumed landslide boundary. Figure 6a and Figure  Overall, the PSI result is consistent with that of SBAS. The deformation rate in the middle part of the slope is the largest, and has an average range of −55 to −80 mm/year, whereas that of the upper part and toe area of the slope is relatively small, having a range of −5 to −20 mm/year.  Figures 6b and 7b show the temporal evolution of the LOS displacement at different points over the observation period. The displacement of point 1 outside the presumed landslide area is almost 0. Point 6 has the largest cumulative displacements of −142 mm and −136 mm using PSI and SBAS methods, respectively. However, point 2 (upper part of the slope) and points 3, 4, and 5 (toe area) display roughly the same displacement between −50 and −70 mm. Comparing the InSAR results of the two techniques (Figures 6 and 7), we can observe that, in the case of PSI, only accumulation areas that are generally located on gentle slopes are detected, whereas SBAS provides the deformation rates of the area with dense vegetation and the steepest sectors of the landslide. Thus, we believe that the SBAS technique is more effective than PSI at detecting landslides, particularly in rugged mountainous areas [48]. Figure 8 shows the time-series deformation along the radar LOS based on the SBAS method. From April 2018 to Sep 2018, there was no obvious cumulative deformation in the slope area. Subsequently, the obvious cumulative deformation signals occurred. In the middle part of the slope, in particular, the cumulative deformation gradually increased, and the deformation rate remained 5.6 mm/month. By July 2019, the total cumulative displacement reached 80 mm. Subsequently, the deformation rate of the region increased to 6.3 mm/month, and the cumulative deformation reached 160 mm by 21 April 2020. In contrast, the deformation rate in the upper part of the slope was relatively small. From July 2019, the average deformation rate was 2.1-2.5 mm/month.   Figures 6b and 7b show the temporal evolution of the LOS displacement at different points over the observation period. The displacement of point 1 outside the presumed landslide area is almost 0. Point 6 has the largest cumulative displacements of −142 mm and −136 mm using PSI and SBAS methods, respectively. However, point 2 (upper part of the slope) and points 3, 4, and 5 (toe area) display roughly the same displacement between −50 and −70 mm. Comparing the InSAR results of the two techniques (Figures 6 and 7), we can observe that, in the case of PSI, only accumulation areas that are generally located on gentle slopes are detected, whereas SBAS provides the deformation rates of the area with dense vegetation and the steepest sectors of the landslide. Thus, we believe that the SBAS technique is more effective than PSI at detecting landslides, particularly in rugged mountainous areas [48]. Figure 8 shows the time-series deformation along the radar LOS based on the SBAS method. From April 2018 to Sep 2018, there was no obvious cumulative deformation in the slope area. Subsequently, the obvious cumulative deformation signals occurred. In the middle part of the slope, in particular, the cumulative deformation gradually increased, and the deformation rate remained 5.6 mm/month. By July 2019, the total cumulative displacement reached 80 mm. Subsequently, the deformation rate of the region increased to 6.3 mm/month, and the cumulative deformation reached 160 mm by 21 April 2020. In contrast, the deformation rate in the upper part of the slope was relatively small. From July 2019, the average deformation rate was 2.1-2.5 mm/month. Figures 9 and 10 display the LOS velocity and cumulative displacement along the longitudinal profile of the slope. It can be observed that the lithology of the upper part of the slope is Carboniferous and Permian basalt and limestone (C + P). The deformation rate in this area is relatively small, with a range of about −30 to −10 mm/year, and the cumulative displacement is between −20 and −60 mm (Figure 10). The lithology of the middle part of the slope body is Devonian phyllite with quartzite (Dwg 1 , Dwg 2 ). In the middle part of the slope, the deformation rate increases rapidly, with an average of about −60 mm/year. The maximum deformation rate is −86 mm/year, and the cumulative displacement is between −120 and −160 mm ( Figure 10). With the decrease in slope height, the deformation rate gradually decreases. At the toe of the slope, the cumulative displacement is almost 0 ( Figure 10). Figure 11 shows the deformation rates on three latitudinal profiles of the slope. Evidently, the deformation rate has a decreasing trend from the middle to each side. In the upper and foot parts of the slope, the rate is relatively small, with an average ranging from −20 mm/year to −10 mm/year. By comparison, the middle part with the largest displacement is about 200 m wide, in which the deformation rate is about 70 mm/year, implying that this area is unstable where the lithology is Devonian phyllite with quartzite (Dwg 1 , Dwg 2 ).   of the slope body is Devonian phyllite with quartzite (Dwg1, Dwg2). In the middle part of the slope, the deformation rate increases rapidly, with an average of about −60 mm/year. The maximum deformation rate is −86 mm/year, and the cumulative displacement is between −120 and −160 mm (Figure 10). With the decrease in slope height, the deformation rate gradually decreases. At the toe of the slope, the cumulative displacement is almost 0 ( Figure 10).    Figure 7a. the slope, the deformation rate increases rapidly, with an average of about −60 mm/year. The maximum deformation rate is −86 mm/year, and the cumulative displacement is between −120 and −160 mm (Figure 10). With the decrease in slope height, the deformation rate gradually decreases. At the toe of the slope, the cumulative displacement is almost 0 ( Figure 10).   upper and foot parts of the slope, the rate is relatively small, with an average ranging from −20 mm/year to −10 mm/year. By comparison, the middle part with the largest displacement is about 200 m wide, in which the deformation rate is about 70mm/year, implying that this area is unstable where the lithology is Devonian phyllite with quartzite (Dwg1, Dwg2).

Prediction of Potential Unstable Rock Mass
Combined with UAV images, we selected the above-mentioned potential unstable rock mass (motion rate < −70 mm/year) to simulate the kinematic and accumulation characteristics. The area of this rock mass is 46,000 m 2 (Figure 12a). Based on the "volumearea" power law [49], the volume was estimated to be 550,000 m 3 with an average thickness about 12 m. To evaluate the influence of different material properties on the movement and deposit characteristics, we conducted a sensitivity analysis using different simulations sets; in each, a single parameter was varied while others took the fixed values of the defined variation range. Table 1 lists the ranges of input parameters. The setting of rheological parameters was estimated by following previous studies [22,50,51]. The simulation results showed that Chèzy coefficients and fluid rate influence the model output, including the sliding velocity and movement process of unstable rock mass, but the predicted landslide boundaries and depositing characteristics of the landslide are similar. Compared with the Chèzy efficiency, fluid rate has less effect on the accumulation distribution and movement characteristics [47]. Overall, for the Chèzy coefficient and fluid rate, reasonable parameters based on previous studies [22,47,50,51] ensure the accuracy of the simulation results. Thus, in the text, we discuss one of the simulation results with the Chèzy coefficient and fluid rate of 200 m·s −2 and 10 m·s −1 , respectively. More details about the sensitivity analysis can be found in the Supplementary Materials.

Prediction of Potential Unstable Rock Mass
Combined with UAV images, we selected the above-mentioned potential unstable rock mass (motion rate < −70 mm/year) to simulate the kinematic and accumulation characteristics. The area of this rock mass is 46,000 m 2 (Figure 12a). Based on the "volumearea" power law [49], the volume was estimated to be 550,000 m 3 with an average thickness about 12 m. To evaluate the influence of different material properties on the movement and deposit characteristics, we conducted a sensitivity analysis using different simulations sets; in each, a single parameter was varied while others took the fixed values of the defined variation range. Table 1 lists the ranges of input parameters. The setting of rheological parameters was estimated by following previous studies [22,50,51]. The simulation results showed that Chèzy coefficients and fluid rate influence the model output, including the sliding velocity and movement process of unstable rock mass, but the predicted landslide boundaries and depositing characteristics of the landslide are similar. Compared with the Chèzy efficiency, fluid rate has less effect on the accumulation distribution and movement characteristics [47]. Overall, for the Chèzy coefficient and fluid rate, reasonable parameters based on previous studies [22,47,50,51] ensure the accuracy of the simulation results. Thus, in the text, we discuss one of the simulation results with the Chèzy coefficient and fluid rate of 200 m·s −2 and 10 m·s −1 , respectively. More details about the sensitivity analysis can be found in the Supplementary Materials. Figure 12 illustrates the simulated velocity of the potentially unstable rock mass. The whole sliding process lasts about 100 s. At 10 s, the material in the source area starts to slide down from the shear outlet. At 30 s, the mass rushes down streams with a high velocity of more than 30 m·s −1 . At 40 s, a small amount of material begins to reach the foot of the slope, while the remainder of the material extends downward at a high speed. At 60 s, the landslide materials gradually accumulate on the riverbed, and the sliding velocity decreases rapidly. At 80 s, most mass is deposited at the foot of the slope, blocking the river channel and forming a barrier dam. Little material remains on the slope surface. At 100 s, the accumulation process ends. the slope, while the remainder of the material extends downward at a high sp s, the landslide materials gradually accumulate on the riverbed, and the slidin decreases rapidly. At 80 s, most mass is deposited at the foot of the slope, bl river channel and forming a barrier dam. Little material remains on the slope 100 s, the accumulation process ends.

Possible Triggering Factors of Old Landslides
In the study area, many old landslides are present along the banks of th River, of which some deposits and unstable slopes had large continuous dis after the 2008 Wenchuan Mw7.9 earthquake [24]. The most frequent landslides occurred in the Devonian and Silurian strata, particularly Devonian phyllite w ite (Dwg1, Dwg2) [27]. The InSAR data show that most areas with high motio located in this stratum (Dwg1, Dwg2) ( Figure 13). From the meteorological sta study area, we collected the data of average monthly rainfall in the Xuecheng and analyzed the relationship between the LOS deformation rate and precipita randomly selected PS points near point 6 (locations are shown in Figure 6a) ( The results show that there is an obvious correlation between the LOS deform and rainfall (Figure 14). In rainy seasons, particularly from May to July, the d  In the study area, many old landslides are present along the banks of the Zagunao River, of which some deposits and unstable slopes had large continuous displacement after the 2008 Wenchuan Mw7.9 earthquake [24]. The most frequent landslides in this area occurred in the Devonian and Silurian strata, particularly Devonian phyllite with quartzite (Dwg 1 , Dwg 2 ) [27]. The InSAR data show that most areas with high motion rates are located in this stratum (Dwg 1 , Dwg 2 ) ( Figure 13). From the meteorological station of the study area, we collected the data of average monthly rainfall in the Xuecheng area [27], and analyzed the relationship between the LOS deformation rate and precipitation at five randomly selected PS points near point 6 (locations are shown in Figure 6a) (Figure 14). The results show that there is an obvious correlation between the LOS deformation rate and rainfall (Figure 14). In rainy seasons, particularly from May to July, the deformation rate was relatively high. From May to July in 2018, the deformation rate in some periods was as large as 1 mm/day ( Figure 14). Considering the geological and rainfall conditions of the area, our research suggests that the possible driving factors of reactivation of the XC landslide include: (1) Weathering of hard basalt and conglomeratic limestone (C + P) is slower than that of soft phyllite (Dwg 1 , Dwg 2 ). Thus, the slope of the hard rock area is steeper than that of the soft rock area, and a greater surface water flow is concentrated in the gentle slope area. Figure 2 shows that the densest hydrological paths are in the middle part of the old landslide. This implies that the flow regime may be confined to the gentle slope area, resulting in a significant increase in groundwater level in response to rainfall.
(2) The water can soften phyllite, of which the mechanical strength decreases greatly, and is thus prone to slope instability under gravity. The rainfall is abundant and concentrated in the study area, which can greatly reduce the strength of this stratum (Dwg 1 , Dwg 2 ), and then induce new slope instability. (3) Finally, river erosion can cause vacant surfaces on both sides of the river, the lowering the landslide's stability. Figure 3 shows the eroded landslide toe and some small-scale collapses developed at the toe of the old landslide, which indicates that the landslide mass has experienced erosion on the concave bank of the river over a long period.

PEER REVIEW 14 of 19
XC landslide include: (1) Weathering of hard basalt and conglomeratic limestone (C + P) is slower than that of soft phyllite (Dwg1, Dwg2). Thus, the slope of the hard rock area is steeper than that of the soft rock area, and a greater surface water flow is concentrated in the gentle slope area. Figure 2 shows that the densest hydrological paths are in the middle part of the old landslide. This implies that the flow regime may be confined to the gentle slope area, resulting in a significant increase in groundwater level in response to rainfall.
(2) The water can soften phyllite, of which the mechanical strength decreases greatly, and is thus prone to slope instability under gravity. The rainfall is abundant and concentrated in the study area, which can greatly reduce the strength of this stratum (Dwg1, Dwg2), and then induce new slope instability. (3) Finally, river erosion can cause vacant surfaces on both sides of the river, the lowering the landslide's stability. Figure 3 shows the eroded landslide toe and some small-scale collapses developed at the toe of the old landslide, which indicates that the landslide mass has experienced erosion on the concave bank of the river over a long period.      Figure 6a).

Impacts of Potential Failure
Both of the banks of the Zagunao River are steep slopes, where many old landslides may be reactivated by heavy rainfall or seismic shaking [25]. The simulation of this work suggests that the unstable rock mass could collapse and form a barrier dam at the Zagunao river with an area of about 0.02 km 2 ( Figure 15). The maximum thickness of the dam is 16 m, and it has an accumulated length of about 800 m along the river. Based on the accumulation characteristics of the barrier dam, the average thickness of the main dam body was estimated to be about 13 m. Then we delineated the scope of the upstream dammed lake using 0.5 m DEM data. Based on ArcGIS, the possible upstream lake area was measured to be 2.5 × 10 4 m 2 , and the volume of the impoundment lake was calculated to be 1.5 m × 10 5 m 3 . Constrained by the narrow valley, some high-speed sliding mass may cross the river and extend to the opposite bank. As a result, the slide may destroy roads and buildings on both sides of the Zagunao River. In addition, because the barrier dam may block the river, the precipitation and surface runoff could flow into the gully, resulting in the rise of the upstream water level, which may threaten the residents. A natural disaster chain of "landslide-barrier lake-flooding" may easily occur in this area [52,53]. After the occurrence of landslides, their deposits may block rivers and form landslide dams. Due to flooding caused by rainstorms, the upstream water level rises rapidly, thus forming a barrier lake. In the western mountainous areas of China, the mountain disaster chain events linking landslides, barrier lakes, and outburst floods have occurred frequently in recent years, such as the 2018 Baige landslide [54] and 2000 Yigong landslide [55], which resulted in high casualties or significant economic loss.
flooding caused by rainstorms, the upstream water level rises rapidly, th rier lake. In the western mountainous areas of China, the mountain disa linking landslides, barrier lakes, and outburst floods have occurred fre years, such as the 2018 Baige landslide [54] and 2000 Yigong landslide [55 in high casualties or significant economic loss.

Application Prospects for InSAR and UAV Techniques
A large number of monitoring results show that the deformation of slide is often characterized by creep slip, and the InSAR technique can e the slow creep deformation signals before the occurrence of the landslid

Application Prospects for InSAR and UAV Techniques
A large number of monitoring results show that the deformation of the ancient landslide is often characterized by creep slip, and the InSAR technique can effectively capture the slow creep deformation signals before the occurrence of the landslides [8,32,56]. As a result, numerous scholars have identified ancient landslides using the InSAR technique and high-resolution aerial photos over large areas [15,57,58]. Based on the high-quality long-time displacement information, the cumulative displacement-time curve of the typical monitoring points can be obtained, which can be used to help discriminate the evolution stage or predict the stability of ancient landslides [57]. This work is of significance to reactivation monitoring and early identification of ancient landslides. In addition, highprecision remote sensing images and terrain data based on UAV flight surveys can be conducive to the identification of the landslide boundary, scarp, and fissures [58]. The surface deformation information obtained using InSAR and UAV techniques can significantly improve the early identification and evolution analysis of ancient landslides at a large scale [57]. In addition, InSAR and UAV techniques can also be used to trace the historical deformation evolution and analyze the potential inducing factors based on actual landslide events [59,60]. This is of significant value for the analysis of landslide triggering factors and post-disaster stability assessment.

Conclusions
This work used two InSAR techniques, PS and SBAS, to monitor the reactivation of a large-scale old landslide in the eastern margin of the Tibetan plateau. Its deformation rate, potential sliding boundary, and driving factors were studied using the Sentinel-1A dataset and UAV data. In addition, a hypothetical collapse body was numerically simulated to predict the potential hazard of this old slope failure. Furthermore, relationships between the LOS deformation rate and the precipitation were analyzed. The main conclusions are as follows.
(1) Both InSAR techniques can provide useful information about the landslide, and SBAS can produce smoother and more abundant displacement time series data. From April 2018 to April 2020, the deformation in the middle of the slope was the largest, and the cumulative displacement was between −120 and −160 mm, with an average LOS deformation rate of −70 mm/year. (2) A strong correlation exists between the LOS deformation rate and rainfall. In rainy seasons, particularly from May to July, the deformation rate was relatively high. From May to July in 2018, the deformation rate reached 1 mm/day. (3) The whole sliding process of the hypothetical collapse body lasts about 100 s, and leads to the collapse and formation of a barrier dam at the Zagunao river. Constrained by the narrow valley, the occurrence of a landslide in this area may cause the natural disaster chain of "landslide-barrier lake-flooding", thus posing a significant threat to upstream and downstream areas. As a result, this possible hazard should receive greater attention.
Our result shows that the application of Sentinel-1 images and InSAR techniques to landslide monitoring is an effective and low-cost method. The combination of the UAV and InSAR techniques can be beneficial for the investigation of inaccessible landslide areas. In future research, more monitoring information from different sources, such as Global Navigation Satellite System (GNSS) data, can be added as a supplement to the InSAR results to verify the long-term deformation characteristics of the old landslide. In addition, further field investigations should be carried out to improve the analysis of the occurrence time and failure mechanism of the old landslide.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/rs13132552/s1, Figure S1: Distribution of collapse body thickness in different Chèzy coefficients, Figure S2:Deposit thickness along longitudinal profile A-A' in different Chèzy coefficients, Figure S3: Predicted landslide boundaries in different Chèzy coefficients, Figure S4: Max sliding velocity in different Chèzy coefficients, Figure S5: Simulated sliding velocity of the unstable rock mass with Chèzy coefficient and fluid rate being 800 m/s and 10 m/s, Figure S6: Maximum sliding velocity and maximum thickness versus fluid rates.
Author Contributions: C.X. proposed the research concept and participated in data curation and analysis. S.M. designed the framework of this research, processed the InSAR data, and wrote the manuscript. X.S. participated in data analysis and revised the manuscript. X.X. and A.L. reviewed and edited the paper. All authors have read and agreed to the published version of the manuscript.