Landslide Hazard Prediction Based on UAV Remote Sensing and Discrete Element Model Simulation—Case from the Zhuangguoyu Landslide in Northern China

Abstract

Rainfall-triggered landslides generally pose a high risk due to their sudden initiation, massive impact force, and energy. It is, therefore, necessary to perform accurate and timely hazard prediction for these landslides. Most studies have focused on the hazard assessment and verification of landslides that have occurred, which were essentially back-analyses rather than predictions. To overcome this drawback, a framework aimed at forecasting landslide hazards by combining UAV remote sensing and numerical simulation was proposed in this study. A slow-moving landslide identified by SBAS-InSAR in Tianjin city of northern China was taken as a case study to clarify its application. A UAV with laser scanning techniques was utilized to obtain high-resolution topography data. Then, extreme rainfall with a given return period was determined based on the Gumbel distribution. The Particle Flow Code (PFC), a discrete element model, was also applied to simulate the runout process after slope failure under rainfall and earthquake scenarios. The results showed that the extreme rainfall for three continuous days in the study area was 151.5 mm (P = 5%), 184.6 mm (P = 2%), and 209.3 mm (P = 1%), respectively. Both extreme rainfall and earthquake scenarios could induce slope failure, and the failure probabilities revealed by a seepage–mechanic interaction simulation in Geostudio reached 82.9% (earthquake scenario) and 92.5% (extreme rainfall). The landslide hazard under a given scenario was assessed by kinetic indicators during the PFC simulation. The landslide runout analysis indicated that the landslide had a velocity of max 23.4 m/s under rainfall scenarios, whereas this reached 19.8 m/s under earthquake scenarios. In addition, a comparison regarding particle displacement also showed that the landslide hazard under rainfall scenarios was worse than that under earthquake scenarios. The modeling strategy incorporated spatial and temporal probabilities and runout hazard analyses, even though landslide hazard mapping was not actually achieved. The present framework can predict the areas threatened by landslides under specific scenarios, and holds substantial scientific reference value for effective landslide prevention and control strategies.

1. Introduction

Rainfall-induced landslides are one of the most frequent and deadly natural hazards around the world, and have been reported to cause considerable damage on a yearly basis [1,2,3]. China is the country most severely affected by landslides in Asia, due to its large-scale mountainous terrain, complicated climate conditions, and massive engineering activities [4,5]. The recent EMDAT database indicates that the number of events and people impacted by climate-related disasters in China rank at the forefront among countries [6]. Hence, it is important to develop and implement a landslide hazard assessment framework for local authorities and researchers. This is widely considered to be a key step in landslide investigations, risk management, and mitigation, since it can determine the spatial and temporal probability of the occurrence of landslides in target areas, along with their mode of propagation, size, and intensity [7].

Generally, the standard procedure for landslide hazard assessment involves three aspects, namely the susceptibility [7], temporal, and magnitude probabilities of landslide events [8,9]. However, a complete landslide hazard assessment is an operational challenge commonly faced by the scientific community in practice. This is mainly associated with data availability, especially the severe scarcity of data on landslide date and magnitude [10]. Hence, most quantitative methods within this topic are related to the spatial probability assessment of landslides, which is also expressed as susceptibility in the literature, conducted at a various scales, ranging from the regional to slope scale. For large areas, many methods have been developed, which can be roughly divided into three categories, namely heuristic models [11], data-driven models [12,13], and physical-based models [14,15]. For local or small scales, the slope failure probability is often determined by using numerical analysis, such as the Monte Carlo method [16,17], and various simulation software [18,19,20], where the input parameter uncertainties of soil materials and random rainfall patterns are often combined.

Regarding the temporal probability of landslides, this can be determined by analyzing the historical activities of a single hazard, where the exceedance probability or occurrence frequency are estimated as the real probability of the landslide [21,22]. However, this is difficult for regions without a complete landslide inventory or historical records. Hence, some studies utilize the temporal probability of triggering events as a proxy, mainly including rainfall and earthquakes [23,24]. In addition, the probability of whether a landslide can reach the elements at risk should also be also clarified, namely by so-called runout analysis, especially for spatially explicit cases [7]. Numerical simulation is normally employed as the method of choice in this topic, and many numerical models based on various mathematical foundations have been proposed and developed, including but not limited to PFC [25,26], SPH [27], DDA [28], UDEC [29], and their coupling methods. Recently, rapid mass movement simulation [30] and the depth-integrated continuum method [31] were also reported in the literature. During the analysis procedure, the deposit depth, velocity, impact force, and energy were the most important indicators for characterizing the landslide intensity.

Surface data on landslides, which can directly reflect geomorphological features such as terrain, plants, and infrastructure, are necessary for hazard assessment strategies [32,33]. Compared with ground-based instruments, an increasing number of studies are using remote sensing methods with a high spatial resolution to obtain these data, for instance, manned aircraft and satellite remote sensing [34,35]. However, atmospheric conditions and timely data acquisition are common issues when it comes to these methods [36]. Hence, researchers are exploring the use of manned-control UAVs for sensing, the photogrammetry of landslide surfaces, and obtaining beyond-surface data, taking advantage of its flexible, effective, and low-cost characteristics. The most recent research regarding on-site landslide studies indicates that UAVs have almost become a necessary tool to support investigations [37,38]. On the one hand, UAVs can get closer to slopes in a contactless way, thus offering a rapid response option in the aftermath of natural hazards. On the other hand, UAVs can be employed for ground displacement monitoring and geological structures to characterize both the dynamic and static facets of landslides [39,40,41]. Additionally, it should be noted that UAVs should be used for mapping/inventorying recent landslides, not for old landslides, where they will be lost.

The main aim of the present study was to develop a framework regarding the prediction and assessment of landslide hazards under extreme conditions (rainfall and earthquakes), where remote sensing and numerical simulation were the most relied upon tools. The ZhuangGuoYu landslide (ZGYL) in Tianjin city of northern China was selected as the study area. Specifically, our objectives included the following: (i) obtaining high-resolution topography data by employing a UAV, (ii) determining the failure probability of the landslide under a specific scenario based on a numerical simulation, and (iii) forecasting the landslide kinematic process after sliding and comparing the landslide hazard under various scenarios.

2. Study Area

2.1. General Settings

ZhuangGuoYu landslide (ZGYL) is located in the Jizhou District of Tianjin city in northern China (Figure 1), with center coordinates of 117°18′36″E, 40°07′23″N. The elevation of the landslide ranges from 110 ASL to 205 ASL, with a length of 220 m and a width of ~100 m. Its average slope angle is 30°, and the total area of the landslide is 2.1 × 10⁴ m², with an approximately 11 × 10⁴ m³ volume. Regarding geological data, the sliding materials mainly contain artificial filling soil, Quaternary Holocene sediments, and Quaternary Holocene alluvial sediments. The bedrock is the Wumishan Formation of Jurassic (J_xw), and the outcrop lithology mainly involves dolomites and dolomitic limestone, where a brownish grey flint band can be found occasionally. According to a field survey, the rock–soil mass has been strongly to extremely weathered, thus, the sliding surface is the boundary between the upper layer and the bedrock (Figure 2). In addition, considering that the landslide has been undergoing continuous deformation in recent years, the local civil administration conducted a drilling investigation in the field. The drilling data revealed that the average depth of the slope was approximately 4.5 m.

The study area has a temperate monsoon climate, characterized by distinct seasons and obvious seasonal fluctuations regarding dry and wet conditions. The topography is higher in the north and lower in the south, with an average temperature of 11.6 °C for the whole year and 24 °C for summer. The long-term annual precipitation ranges from 351 mm to 1213 mm, with a mean value of 678.6 mm. An extremely uneven distribution regime for precipitation is observed, where more than 80% of the total annual precipitation happens during June and September.

The area is highly urbanized, with a village at the toe of the landslide, which has a population of approximately 1100 inhabitants. This village is the main area affected by the landslide, especially considering its potential for failure. Human activities in the area have developed rapidly over recent decades. On the one hand, the settlement quickly expanded, and important transportation infrastructure was also constructed (e.g., Highway G1N and G102). On the other hand, the region has alternating hills and plains, thus, a large amount farmlands have been developed.

2.2. Landslide Details

The area was once a small-scale landslide at the toe of the slope after heavy rainfall in the summer of 2015, with a volume of more than 1000 m³. The plants have been damaged here, but other sections of the slope are still covered by trees. The local civil administration set a warning point for geohazards here, but no professional monitoring instruments were implemented, besides a simple protective net (Figure 3).

To illustrate the ground deformation of the slope, SBAS-InSAR techniques were utilized to obtain the displacement of the ZGYL over the past decade (see Section 3.2 for technical details). It was found that the pixel displacements within the slope were relatively larger than those beyond the slope (Figure 4). During the past decade, the general trend of deformation has gradually increased. In particular, since the year of 2020, the cumulative displacement has exceeded 30 mm, thus indicating that the ZGYL is under the creeping stage.

Regarding the landslide mechanism, the investigation report from local authority mentioned that the deformation was mainly triggered by rainfall during the past years. This was determined from the three following main aspects: (i) No extreme events were recorded in the study area during the last decade, including earthquakes, typhoons, tropical cyclones, blasting, etc. (ii) It was not possible that potential engineering activity (slope cutting, etc.) disturbed the slope stability, given that the site has been protected. (iii) The continuous deformation by SBAS-InSAR leads us to believe that the triggering factor for activating the landslide should have long-term effect, rather than short-term effects. However, it is difficult to analyze the quantitative relationship between landslide deformation and rainfall, since detailed rainfall data on a daily/monthly basis on-site is not available. In addition, it should be noted that Tianjin city is located along the North China Plain seismic belt, which is a seismically active region in the eastern Chinese mainland [42]. Several strong earthquakes have been reported in this region, including the 1966 Ningjin M7.2 earthquake (~300 km from Tianjin) and the 1976 Tangshan M7.8 earthquake (~100 km from Tianjin), causing huge casualties and economic losses [43]. Studies have confirmed that the region poses a high risk of earthquake hazards, with the potential for strong earthquakes in the future when stress accumulation reaches a certain level [44,45]. This leads us to believe that an earthquake may also trigger the catastrophic failure of the ZGYL, although its current creeping deformation is not associated with this. Hence, both extreme rainfall and earthquake scenarios should be taken into account in stability analysis and hazard assessment.

3. Materials and Methodologies

3.1. Modeling Framework

Overall, the modeling framework proposed in this study includes several steps, as follows (Figure 5): (i) the Gumbel distribution is used to determine extreme rainfall scenarios based on historical rainfall data; (ii) the UAV remote sensing technique is applied in the field to generate high-resolution topography data; (iii) the Geostudio software v2024 is used to compute the factor of safety and the probability failure of the slope under different scenarios; (iv) the results regarding the stability analysis are used to determine the input parameters for the PFC; (v) a runout analysis of the landslide is conducted by PFC modeling, and kinetic indicators are analyzed; and (vi) the landslide hazard under different scenarios is assessed and compared. The detailed method for each step will be introduced in the following sections.

3.2. Deformation from SBAS-InSAR

The small baseline subset (SBAS) interferometric synthetic aperture radar (InSAR) method is a technique used for measuring ground deformation over large areas, which was proposed in 2002 [46]. It has been widely applied for landslide deformation monitoring due to its high accuracy and temporal continuity, especially for high mountain regions, where physical accessibility is seldom achieved [35,47].

In this study, the raw images used were Sentinel-1A data, which were downloaded from the website of the Alaska Satellite Facility (https://asf.alaska.edu/, accessed on 10 September 2024). The data spanned ascending and descending track Sentinel-1A SLC images from November 2014 to November 2023. Since the calculated deformation is along the LOS direction of the satellite, we converted the line-of-sight deformation to the deformation along the slope direction of the landslide. To achieve this goal, the ascending and descending tracks were first combined to obtain the two-dimensional deformation. During this procedure, Precise Orbit Ephemerides (POD) data from the European Space Agency were selected for orbital position correction.

3.3. UAV Remote Sensing

In this study, the DJI Mavic3E UAV with an RTK module was utilized to conduct high-resolution remote sensing of the landslide. The important parameters for the UAV and camera are shown in

Table 1. The important parameters of the DJI Mavic3E UAV.

Weight	Duration for Single Fight	Cruising Speed	Max Climbing Speed	Pixel of Camera	Resolution of Images
0.915 kg	0.5 h	15 m/s	8 m/s	2 × 107	5280 × 3956

. The technical procedure was mainly as follows: (i) We determined the working region in Google Earth and saved it as a KMZ file, which was subsequently imported into the UAV. (ii) We set the controlling points, which were evenly distributed in the region, as well as the UAV routes. It should be noted that two overlapping routes were employed to achieve the oblique photography. Since the study area had a relatively undulating terrain and the features of the ground objects were not obvious, the heading overlap was set to 83% and the lateral overlap was 75% when planning the route to ensure smooth data processing. (iii) The Pix4Dmapper software v4.5.6 was used to process the image data and obtain dense point cloud data to generate a digital surface model (DSM), orthophotos, and oblique models of the survey area (Figure 6). These data can not only be used to present the overall terrain conditions of the landslide, but also to conduct accurate numerical analysis due to their high resolution (0.08 m resolution in this study).

3.4. Extreme Rainfall Analysis

Landslide hazard assessments should be based on the determination of scenarios [7,48]. For the ZGYL, rainfall is the most important factor triggering instability, hence, extreme rainfall scenarios should be analyzed to predict the potential landslide stability and hazard in advance. Generally, there are two main types of situations when considering the temporal probability of triggering events. The first one is designed (specific) events according to historical data. However, this may be subjective, because the rainfall or cumulative rainfall thresholds for activating landslides are sometimes determined by expert knowledge or statistical results. The second deals with return periods, the concept of which has been widely used to illustrate the degree of cumulative rainfall or rainfall intensity in the literature [24,49,50]. Considering that the main objective of the present study is to predict the landslide hazard under extreme scenarios, we adopted a return period for determining rainfall scenarios. The main idea in calculating the return period of rainfall is to fit the frequencies of historical rainfall and then predict the extremes within a given probability. In this study, the long-term historical rainfall during 1980 and 2017 was applied to generate the extreme rainfall of the study area (Figure 7). An analysis regarding the number of days with rainfall for every single month showed that three continuous days of rainfall was the most frequent event during the rainy season in the study area. The largest continuous 3-day rainfall for each month was taken into account and then used to calculate the rainfall versus the return period with the Matlab program based on the Gumbel distribution.

The probability density function f(x) and cumulative density function F(X) for the Gumbel distribution are as follows [51,52]:

$$f(x)=({\displaystyle \frac{1}{\alpha }})exp[-({\displaystyle \frac{x-\beta }{\alpha }})]-exp[-({\displaystyle \frac{x-\beta }{\alpha }})]$$

(1)

$$F(x)=exp[-exp[-({\displaystyle \frac{x-\beta }{\alpha }})]$$

(2)

where α is the scale parameter and β is the location parameter. Their relationship can be denoted as follows:

$$\mu =\beta +0.5772\alpha$$

(3)

where μ is the mean value of the distribution. The parameters α and β can be obtained through maximum likelihood estimation.

3.5. Landslide Stability Evaluation and Runout Analysis

3.5.1. Failure Probability of the Landslide

As mentioned above, landslide hazard involves the occurrence probability, which was determined by numerical simulation in Geostudio software v2024 in this study. It is based on the finite element method, and can compute both the seepage (ground water level, pore water pressure, etc.) and stability of a slope. In particular, the software can obtain the failure probability by considering the spatial variability of soil properties. For the ZGYL, the modeling and analysis procedure were conducted as follows:

(i)

A geological model of the landslide was first established with grid shapes of quadrilaterals and triangles. Considering the geometry and volume of the slope, the size of the grids was set to 0.5 m, thus determining a total of 6027 nodes and 5866 elements (Figure 8).

(ii)

The whole slope was divided into a sliding body and bedrock according to the cross-section of the landslide, and the hydraulic boundary conditions were added according to the results from the extreme rainfall analysis. Because the on-site engineering investigation did not reveal the ground water level, the eater level was not set. The seepage field of the slope under rainfall scenarios was computed in the SEEP/W module of the software, and the results from this step were taken to conduct the stability analysis using the SLOPE/W module. During this procedure, when the deterministic parameters were used as the inputs, the factor of safety (FS) was obtained. The failure probability was computed when the distribution functions of the shear strength parameters were defined (i.e., stochastic parameters). In this study, the deterministic values for the parameters were set according to

Table 2. The parameters for the slope stability analysis in Geostudio.

Section	Unit Weight (KN/m3)	Cohesion (KPa)	Friction Angle (°)	Elastic Modulus (KPa)	Poisson’s Ratio	Volumetric Moisture Content
Sliding body	20	30	22	1200	0.3	0.4
Bedrock	21	50	30	20,000	0.33	0.3

, which were mainly provided by the archived landslide report from the local authority of the geo-environmental department. The cohesion and friction angle followed a normal distribution based on the existing literature [18], where the max and min values were determined by laboratory tests.

(iii)

As analyzed in Section 2.2, the landslide is under the slow-moving stage, and future earthquake scenarios will possibly cause catastrophic failure. In China, the seismic hazard assessment and seismic design of buildings are commonly based on seismic intensity zonation maps [53]. Tianjin city is classified as the VII level in the latest version of the Seismic Ground Motion Parameters Zonation Map of China, representing a peak ground acceleration of 0.15 g [54]. Regarding the seismic duration, the assessment in this study was not for a specific event, but for a general scenario. We referred to the most famous earthquake within the study area, namely the Tangshan M7.8 earthquake, which lasted 14–16 s (https://www.britannica.com/event/Tangshan-earthquake-of-1976, accessed on 10 September 2024). Hence, the duration of the earthquake scenario for the ZGYL was set as 20 s. For a comparison, the duration of the Luding M6.8 earthquake in China in 2022 was ~20 s [55].

3.5.2. Landslide Runout Analysis Based on Discrete Element Model

For a slow-moving landslide, a high risk is mainly posed by rapid movement after catastrophic failure, rather than the creeping stage. Hence, the PFC was used to conduct the runout analysis in this study, which determined the spread distance of the slope. The PFC is a numerical simulation tool based on the discrete element model, capable of modeling both the macroscopic and microscopic features of rocks or soils. It can simulate the entire process from separation to final failure in rock–soil masses, fitting well with the landslide initiation and movement in the real-world scenarios. Hence, the PFC is a powerful tool for landslide runout analyses and kinematic characteristics [56].

Geological models of the ZGYL were established in the PFC at the 2D and 3D scales, respectively (Figure 9). However, we were not aiming to confirm the PFC2D conclusion by the results of PFC3D, especially given that it was impossible for the parameters of two calculations to be totally the same. These calculations were performed to compare the relative trends of kinematic indicators under different scenarios. The reliability of the results can be guaranteed if the trends fit well with each other between two scaling operations. The porosity was set to 0.2 to randomly generate particles, and the sliding body and bedrock layers were contacted in the ball–ball manner for 2D and the ball–wall manner for the 3D modeling. The particle sizes of the sliding body ranged from 0.4 to 0.7 m, whereas those for the bedrock were 0.5–1 m. The total numbers of generated particles for the ZGYL were 8195 (2D) and 8589 (3D), respectively.

It should be noted that the initiation moment and movement process of landslides are influenced by rainfall magnitude and intensity, so they normally exhibit different kinetic characteristics under different scenarios. However, the PFC struggles to analyze rainfall infiltration and seepage, thus making it unable to reflect the effects of rainfall or earthquakes on mechanical parameters. Therefore, the strength reduction method was applied to the parameters in this study when simulating the landslide kinetics to distinguish different scenarios. Given that the Geostudio software v2024 can simulate the slope stability, the strength reduction method integrates the advantages of Geostudio and the PFC. In detail, the factor of safety (Fs) of the slope was first calculated by using Geostudio under different scenarios. The Fs value was considered as the ratio of the stability between the natural scenario (no rainfall or earthquake) and the triggering one. Next, the coefficient was utilized to reduce the natural parameters from the natural scenario (

Table 3. The important original parameters for the runout analysis in PFC.

Section	Density (kg/m3)	Cohesion (KPa)	Friction Angle (°)	Bulk Modulus (MPa)	Shear Modulus (MPa)	Poisson’s Ratio
Sliding body	21	30	22	60	50	0.35
Bedrock	21	50	30	300	100	0.35

). For example, if the FS obtained from Geostudio under the rainfall scenario with 50-year return period was 0.5, the input parameters for the PFC simulation should be half of that for the natural scenario.

As seen in Figure 9, a total of four particles were selected in the 2D model to monitor the dynamics of the velocity and displacement. They cover different parts from the top (#1) to the toe (#4) of the profile, thus, they can reflect the detailed kinetic characteristics at different locations. We did not analyze the average values of all particles, since calculation is time consuming when it comes to a large amount of particles. In addition, the analysis regarding a specific particle is considered as more meaningful, given that a particle can be affected by the impacts or push forces from nearby particles.

4. Results

4.1. Extreme Rainfall and Landslide Stability Analysis

The results of the extreme rainfall versus the return period obtained from the Gumbel distribution are shown in Figure 10. It can be seen that the rainfall with a return period of 20 years, 50 years, and 100 years was 151.5 mm, 184.6 mm, and 209.3 mm, respectively. Given that most rainfall events occurring within one single day were concentrated in several hours, it can be inferred that the rainfall intensity in the study area could even reach 100 mm/h. Subsequently, rainfall events with 50-year and 100-year return periods were set as the hydraulic boundary conditions, and the safety of factor and failure probability under each scenario were determined by using Geostudio. The results showed that the Fs values were 0.897 (rainfall with a 50-year return period), 0.819 (rainfall with a 100-year return period), and 0.956 (earthquake), respectively. The failure probabilities were 84.3% (rainfall with a 50-year return period), 92.5% (rainfall with a 100-year return period), and 82.9% (earthquake), respectively (Figure 11). The results indicated that both heavy rainfall and earthquake scenarios posed a high risk of inducing instability in the ZGYL, and the stability of the landslide was worse under heavy rainfall than under earthquake scenarios. Regarding the displacement, the max values under rainfall scenarios were 0.73 m and 1.22 m, respectively. It should be noted that these were probably not the actual displacements in practice, because the finite element method is not suitable for simulating large deformation. In fact, when the slope had a deformation of more than 1 m, the slope failure may have already happened. Under the earthquake scenario, the displacement was only 0.09 m, which was mainly because some of the forces provided by horizontal seismic acceleration were towards the slope, which improved the slope stability.

4.2. Landslide Runout and Kinetic Analysis

The PFC 2D simulation results of the ZGYL showed that the duration of the landslide movement lasted approximately 60 s. The movements of the particles at different moments for the scenarios of rainfall with a 100-year return period and earthquake are shown in Figure 12 and Figure 13 for comparison. It can be seen that the movement was fast during the first 10 s after the landslide failure, and many particles reached the flat area at the toe of the slope. From 20 s to 60 s, slow movement and the accumulation stage were mainly observed. In comparison to the rainfall scenario, the particles seemed to accumulate faster under the earthquake scenario, and, in particular, the accumulation pattern was more compact at 20 s. From 20 s to 60 s, the dynamic changes of the particles under the earthquake scenario were smaller than those under the rainfall scenario, which also indicated that the kinetics process under the earthquake scenario was faster. However, the runout distance under the earthquake was smaller than that when it came to the rainfall scenario, thus meaning that the overall speed of the particles under the earthquake scenario was slower.

The results regarding the PFC3D also confirmed the conclusions from the 2D simulation. The kinetic processes under different scenarios are shown in Figure 14. It should be mentioned first that the color bars in the subfigures are different because the magnitudes of the kinetic indicators varied with time. It can be seen that the runout distance was larger when the return period of rainfall was larger, with more particles accumulating at the toe of the slope rather than stopping on the slope. Next, the velocity and total displacement of the four monitoring particles were obtained and compared. It was found that the initiation moments for different particles varied, which was associated with the positions of the particles, as well as the movement characteristics of the surrounding particles. When the velocity of a specific particle fluctuated frequently or changed suddenly, it may have been affected by the impacts or push forces from other neighboring particles. Hence, the dynamics of the velocities for the monitoring particles are difficult to summarize quantitatively, but at least one peak velocity was observed for each particle (Figure 15). In addition, the max velocity of all the particles under the rainfall scenario with a 100-year return period (P = 1%) was the largest compared with other scenarios, thus indicating the highest hazard under this scenario. Regarding the runout distance (Figure 16), the max displacement of the particles also had the peak value when it came to the rainfall scenario with P = 1%. The max displacements for the rainfall with P = 2% for the four particles were 122 m (#1), 123 m (#2), 144 m (#3), and 146 m (#4), respectively; the max displacements were 132 m, 142 m, 149 m, and 162 m, respectively for the rainfall with P = 1%. Under the earthquake scenario, the peak displacements were 121 m, 119 m, 133 m, and 143 m, respectively. Considering that a settlement is located at the toe of the slope, it can be highlighted that the landslide hazard under the rainfall scenario was greater than the earthquake scenario.

5. Discussion

Three major topics are discussed in this section. Firstly, the relationship between the UAV remote sensing and the discrete element model for landslide disaster prediction is presented. Secondly, the modeling and analysis uncertainties within the current landslide hazard prediction framework are addressed. Finally, the results associated with the landslide hazard assessment and zonation will be debated.

UAVs have been one of the most important tools for landslide investigation and monitoring. Particularly, they can be applied both before and after a landslide to obtain information. As it stands, the current methods/models for landslide assessment invariably require multi-temporal landslide information. On the one hand, although it is also a common practice to utilize remote sensing imagery to characterize original slopes prior to an event, this approach faces several technical challenges, such as resolution limitations and availability in terms of time. On the other hand, UAVs can get closer to landslides to generate high-resolution images, and users can use UAVs at any time. On the other hand, UAVs offer a rapid response option in the aftermath of hazard events, particularly when ground investigations are time-consuming. In this study, the most important role of UAV remote sensing was to obtain a high-resolution DSM of the slope, which can be used as critical input data for discrete element simulations. Considering that the ZGYL is currently in the creeping deformation stage without catastrophic motion, it is difficult to determine whether the application of a UAV in our case belongs to the pre- or post-event time period. Certainly, open-source DEM/DSM data are a common option for most studies. For example, the 12.5 m ALOS, 30 m ASTER, and 90 m SRTM terrain data have been widely used worldwide [57,58]. However, the impact of data resolution on outcomes always exists, and some research has investigated the impact of data resolution and the scaling effect [5,59]. In fact, data with a resolution of more than 10 m are normally suitable for cases with a large area, such as regional-scale landslide susceptibility and hazard assessments. Regarding the local or slope scale, the sensitivity of the model to spatially varied parameters and inputs limits its accuracy for application. In this study, UAV remote sensing derived 0.08 m resolution DSM data, which also contained vegetation information. It is apparent that the rock–soil mass will also move along with the plants on the slope, so a DSM is more accurate for landslide hazard simulations than a DEM. In addition, images from UAV remote sensing will help us to recognize the boundaries of elements at risk (roads and houses, etc.) for vulnerability assessment, which needs subsequent work in hazard zonation.

The uncertainties of the results in the present study are related to several important aspects, among which, the that first should be mentioned is the determination of parameters. Two categories of numerical simulation were included in the proposed framework for landslide hazard assessment, therefore, the accurate parameters for hydrological and mechanical properties were the solid foundation for reliable outcomes. However, the parameter values were determined only according to limited survey and laboratory results in our procedure, incorporating the probabilistic distribution of the shear strength parameters in Geostudio rather than their spatial distribution. Indeed, the characterization of spatial heterogeneity for soil properties has been always a technical issue faced by the scientific community [60]. Hence, as reported by Hürlimann et al. [2], techniques that consider the spatial uncertainties of soil properties may be of a high interest and importance. In fact, some methods have been developed to address this issue [61,62], which can improve the accuracy and reliability of the proposed framework. Moreover, the mechanical parameters were determined according to the strength reduction by comparing the factors of safety under different scenarios, which also caused uncertainties in the final results. This is mainly associated with the inherent limitations of numerical simulation methods. Geostudio can compute the failure probability of the slope by considering the probability distribution of the parameters, but it struggles to simulate macro large deformation and movement [18,63]. The PFC is suitable for analyzing landslide kinetics, but it seldom incorporates rainfall infiltration. Therefore, the present framework combined these two tools together to achieve all three components of the landslide hazard concept, namely spatial probability, temporal probability, and landslide kinetics under a given scenario. It should also be mentioned that the incorporation of strength reduction into numerical simulation modeling is a common method for conducting slope stability quantification in the literature [64,65]. Therefore, using the safety factor obtained by Geostudio as a strength reduction factor was considered to be rational in this study. Moreover, given the previous limitations regarding quantitative landslide hazard assessments at the slope scale, the uncertainty during the determination of parameters was acceptable.

Furthermore, one uncertainty with regard to the stability analysis originated from the setting of the extreme rainfall scenario. The rainfall was determined by the statistical data from continuous 3-day rainfall events during rainy season in the region, where the number of days with rainfall fit with the most frequent single rainfall event, and the rainfall intensity was evenly distributed. This was actually a hydrological simplification, without an adequate analysis of the rainfall infiltration and effective recharge. Given the diversity of the evolution of rainfall-induced landslides under random rainfall patterns in the real world, e.g., the evolution time and the cumulative landslide area [66], the uncertainty associated with the volume of the ZGYL may be quite large. In particular, the landslide initiation and failure may be highly variable with respect to the various distributions of rainfall peaks, which was not considered in this case. Certainly, it is also important to admit these limitations rather than evading them, which can help users to improve modeling accuracy in the future. Hence, the present modeling strategy is more likely a guideline for predicting landslide hazards quantitatively under extreme conditions. The outcomes regarding the comparison of the landslide hazards between the rainfall and earthquake scenarios can be utilized to provide early warnings for threatened areas and take measures for risk mitigation. Last but not least, the extreme rainfall under a given return period may also change in the future due to global climate change, especially regarding the aspects related to the frequency and abundance of extreme rainfall events.

Typically, landslide hazard assessment is expressed as zonation or mapping based on GIS, ranked by the probability values in various sections. However, the proposed framework enables revealing the potential threatened regions impacted by mass movement under a given scenario. The spatial and temporal probabilities of the landslide occurrence were computed, but the spatial distribution of the landslide intensity (i.e., the probability of an area being impacted by a given landslide intensity) was not provided in our case. In fact, initiation and runout hazard analyses are the most difficult challenges within landside hazard assessment. The first critical issue is landslide volume characterization under different scenarios [7]. It is evident that the landslide volume varies with the rainfall intensity or earthquake magnitude, but most available codes seldom consider the failure process [67,68]. Users commonly set the sliding surface according to the investigation results, which are consistent under all scenarios. Moreover, every landslide event will produce a different distribution of intensity and different probability of impact based on its dynamics. Hence, iterative calculations for various scenarios may be necessary to determine the probabilities of the occurrence of landslides with different magnitudes at particular locations below the slope, which is highly time-consuming. Given the limitations mentioned above, the present study only simulates landslide kinetics under a specific scenario without hazard zonation. However, the analysis results still revealed the regions that were most likely to be affected, as well as the most important indicators (velocity and displacement) for characterizing the mass movement. In the future, analysis regarding the magnitude–frequency relationship may be a potential solution, in order to create landslide hazard maps for specific return periods [69,70].

6. Conclusions

There is no doubt that landslide hazard assessment is critically important for stakeholders responsible for civil protection and land use planning. In the present work, a framework for landslide hazard assessment at the slope scale was proposed. Given that a complete hazard modeling procedure is commonly ignored, all the three different aspects associated with this topic were incorporated into the framework, namely spatial probabilities, temporal probabilities, and landslide runout analyses. UAV remote sensing was utilized to obtain high-resolution terrain data, and scenarios were determined according to the extreme analysis regarding the rainfall time series. The application results with the ZhuangGuoYu landslide in Tianjin city of China confirmed the framework’s applicability. The results showed that the spatial probabilities of the landslide under extreme rainfall scenarios were larger than those under earthquake scenarios. The failure probabilities computed by Geostudio were 84.3% (rainfall scenario with P = 2%), 92.5% (rainfall scenario with P = 2%), and 82.9% (earthquake scenario), respectively. The landslide kinetics and the regions affected by mass movement were revealed by PFC modeling at a 3D scale. Two important kinetic indicators were analyzed by four monitoring particles. The max velocities ranged from 19.8 m/s to 23.4 m/s, whereas the displacement varied from 119 m to 162 m, and the rainfall with a 100-year return period was still the most dangerous scenario.

Overall, the present model strategy combined multiple methods and tools to achieve landslide hazard assessment. The outcomes indicated that UAV remote sensing, numerical simulation, and extreme value analysis can be integrated to improve the quantitative hazard assessment, even though the “real” landslide hazard zonation was not created. Meanwhile, the existence of multiple sources of uncertainties within this work was also stated, which need to be transferred to end-users.

The proposed framework can be easily replicated in the areas with similar geological contexts. In particular, the prediction of the locations which may be affected by future slope failures under specific scenarios has twofold practical applications. It can be incorporated into the development program of urbanization and spatial planning to avoid the most dangerous locations. Moreover, the outcomes from the hazard assessment can be combined with the spatial distribution of people and assets to identify risk exposure. Finally, we recommend undertaking efforts to develop techniques regarding runout analyses or magnitude–frequency relationships of landslides to further improve quantitative hazard mapping.