Constructing physical-based rainfall landslides prediction model: insights from rainfall threshold curves database of slope units

The commonly used rainfall threshold warning method relies heavily on historical rainfall and landslide inventory data, which limits its applicability in regions that lack these data. While physical methods do not rely on landslide inventories to establish warning criteria, the calculation of the safety factor typically requires considerable time. To address these issues, this study integrates physical methods, rainfall threshold warning methods, and slope units to develop a rapid forecasting model for rainfall landslides at a regional scale. A hydrological analysis technique for slope units based on grid cells was developed to calculate the instability probability of slope units. Then, each slope unit was analyzed under 20 levels of antecedent effective precipitation and nearly 200 combinations of rainfall intensity (I) and duration (D) to derive the key fitting parameters α and β of the I–D curves under various rainfall scenarios. The application results from Fengjie County indicate that the model runs in less than 12 min, with missing alarm and false alarm rates of 11.8 % and 21.1 %, respectively, highlighting its excellent potential for practical application. This study is expected to provide insights for the rapid forecasting of rainfall landslides in the impoverished mountainous regions of developing countries.

1 Introduction

Rainfall-induced landslides at a regional scale are among the most common types of natural hazards worldwide. Reports indicate that in the United States, rainfall-triggered landslides and secondary hazards result in 25–50 fatalities and economic losses of approximately USD 2 billion annually (Spiker and Gori, 2003). This loss is even more severe in developing countries in the Third World (Wang and Zhang, 2021; Wang et al., 2023, 2024). In recent years, numerous studies have indicated that regional landslide forecasting is highly effective for hazard prevention and mitigation. Researchers have developed various rainfall landslide forecasting models based on statistical and physical methods (Aristizábal et al., 2016; Baum et al., 2008; Bezak et al., 2016; Bogaard and Greco, 2018; Cuomo et al., 2021; Liang and Uchida, 2021; Medina et al., 2021; Pinho and Augusto Filho, 2022; Tufano et al., 2021; Wang et al., 2013; Zhang et al., 2021, 2019; Moeineddin et al., 2023; Li et al., 2025). However, there are still several unresolved issues in regional landslide forecasting, making accurate and efficient warnings a significant global challenge.

The first major issue is the selection of forecasting methods. The presented statistical approaches generally depend on historical precipitation and landslide inventory data to construct the rainfall threshold curves. Recently, researchers proposed different types of rainfall threshold curves, including Intensity–Duration (I–D), Rainfall event–Duration (E–D), Rainfall event–Intensity (E–I), Intraday rainfall and antecedent effective rainfall (IR–AER), Intensity–Probability (I–P), and Intensity–Duration–Mean areal rainfall( I–D–MEAR) (Brunetti et al., 2010; Hong et al., 2005; Rosi et al., 2020; Zhuang et al., 2014). The I–D curve is the most extensively used among these types. The I–D curve is typically fitted in either Cartesian coordinates or a double-logarithmic coordinate system, and the equation of the curve is governed by two key fitting parameters, α and β, expressed as follows:

$$\begin{array}{}\text{(1)}& I=\mathit{\alpha }{D}^{\mathit{\beta }}\end{array}$$

where α and β are derived from the statistical analysis of historical rainfall and landslide data.

Studies indicate that statistical methods are applicable in regions with abundant historical records of rainfall landslides because these areas can provide sufficient samples for fitting the I–D curve (Bezak et al., 2016; Hong et al., 2017; Kanungo and Sharma, 2014; Kim et al., 2020; Ma et al., 2015; Marra, 2018; Pradhan et al., 2018). However, in the poor mountainous regions of the Third World, many areas that are severely affected by landslides lack professional monitoring devices and rain gauges, potentially limiting the application of statistical approaches (Zhang et al., 2021, 2019). In contrast, physical methods typically rely on hydrological and mechanical analyses to calculate the safety factors of landslides under different rainfall scenarios, thereby reducing the reliance on historical rainfall and landslide observation data. In regions where landslide inventory data are scarce, physical methods could serve as promising alternatives (Zhang et al., 2021, 2019). However, physical methods require historical landslide data to validate the accuracy of the forecasting results, and the safety factor calculation process typically requires a considerable amount of time. This computational burden increases substantially when considering the stability analysis of thousands of slopes at the regional scale, making it difficult to ensure the efficiency of real-time warnings (Zhang et al., 2021).

The second issue pertains to the selection of prediction unit. Clearly defined prediction units enable residents to identify the specific locations where landslides are likely to occur while also providing guidance for local governments to develop emergency schemes. However, the I–D warning curves derived from statistical methods can only provide general trends of hazards within the region but cannot pinpoint the specific locations of landslide occurrences. Grid cells improve the clarity of the prediction results to some extent, as the specific locations of each grid within the area are well defined (Zhang et al., 2021). Researchers have employed grid cells to establish multiple physical forecasting models such as Shallow Landslide Stability model (Montgomery and Dietrich, 1994), Stability Index Mapping (SINMAP) (Tarboton and Goodwin, 1998). The Three-dimensional Fully Distributed Hydrological model-Safety factor (GEOtop-FS) (Rigon et al., 2006), Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability Analysis (TRIGRS) (Baum et al., 2008), High Resolution Slope Stability Simulator (HIRESSS) (Rossi et al., 2013), Hillslope-scale Shallow Landslide Induced Debris Flow Risk Evaluation (H-slider) (Liang and Uchida, 2021), Open and Distributed Hydrological Simulation and Landslides (SHIA_Landslide) (Aristizábal et al., 2016), Shallow Landslides Instability Prediction (SLIP) (Montrasio and Valentino, 2016), and Fast Shallow Landslide Assessment Model (FSLAM) (Guo et al., 2022). However, the morphology of grid cells does not accurately characterize the topographical features of natural hillslopes (Domènech et al., 2019; Zhang et al., 2021), resulting in a lack of clear geomorphological significance. In practical applications, a natural slope can be segmented into a series of grid cells, in which each grid is assigned a different alert level. This indicates that a high warning level in a grid cell does not mean that the entire slope will experience a slide.

In contrast, slope units can represent the topographical features of landslides more accurately, and their boundaries are easily discernible in field environments. Currently, there are various methods for extracting slope units, including the DEM-based hydrological process analysis method (Turel and Frost, 2011), r.slopeunits method (Alvioli et al., 2020), curvature watershed methods (Yan et al., 2021), MIA-HSU methods (Wang et al., 2019), and multi-scale image segmentation methods (Huang et al., 2021). In recent years, researchers have developed forecasting models utilizing slope units, validating their promising application potential in predicting rainfall-induced landslides (Wang et al., 2023; Zhang et al., 2021). Addressing the issues outlined in regional landslide forecasting, this study focuses on the integration of slope units, physical methods, and rainfall parameterized warning techniques to develop a rapid forecasting model applicable to large areas on a scale of thousands of square kilometers. Within this model, we no longer pay attention to the positional relationship between the rainfall data of a landslide and the I–D curve, but concentrate on the key fitting parameters α and β of the I–D curve for each slope unit. To facilitate this, we developed a rainfall infiltration simulation technique rooted in grid cells within slope units and subsequently utilized physical methods to analyze the instability probability for slope units under different rainfall scenarios. For each slope unit, we designed rainfall scenarios comprising various antecedent rainfall levels combined with hundreds of rainfall intensity and duration combinations. This allowed us to obtain the key parameters α and β of the I–D curves for different rainfall scenarios, thereby constructing a database of parameters α and β under various antecedent precipitation levels. A case study in Fengjie County, in the Three Gorges Reservoir area, was conducted to validate the reliability of the proposed method. This research is expected to provide valuable insights for regional landslide forecasting in impoverished mountainous areas in the developing world.

2 Methodology

2.1 The slope unit extraction method MIA-HSU

Currently, there are various methods for extracting slope units, each requiring distinct key parameters and yielding different extraction results (Wang et al., 2023). Among these, the MIA-HSU method generates slope units with homogeneous slope gradient and aspect within each unit. In recent years, it has been utilized in rainfall landslide early warning systems in the mountainous regions of southwestern China (Wang and Zhang, 2021; Wang et al., 2025). In this study, we employed the MIA-HSU method to extract slope units (Wang and Zhang, 2021; Wang et al., 2019, 2023). In the MIA-HSU method, each HSU(homogeneous slope unit) is defined as a continuous and homogeneous geomorphological entity. This definition implies that terrain features related to slope and aspect are uniform within each HSU, with boundaries indicating transitions in topographical features. The MIA-HSU method consists of two steps. The first step involves partitioning the Digital Elevation Model (DEM) into small regions with homogeneous terrain characteristics. In this step, the MIA-HSU method utilizes terrain curvature analysis to identify ridge and valley regions (Fig. 1a) and then extracts the morphological skeleton lines of ridge and valley areas to characterize topographic relief. Morphological algorithms (such as dilation and erosion) were used to extract the morphological skeletons of ridges, valleys, and flat areas from the DEM (Fig. 1b), ultimately connecting these skeletons into a closed network (Fig. 1c). Thus, each small region within the network exhibits uniform geomorphological characteristics. The second step involves merging small adjacent regions. The MIA-HSU method employs the principal component analysis (PCA) method to derive fitted planes from localized terrain regions, followed by the implementation of vector similarity criteria to merge adjacent small regions, thereby generating HSUs (Fig. 1d).

Figure 1 HSU extraction process: (a) the identification of ridge and valley areas; (b) the morphological skeleton line extraction for ridge and valley areas; (c) the morphological skeleton closed network; (d) HSU extraction result.

2.2 The HSU hydrologic simulation technique based on grid cells

2.2.1 The identification for row and column information of grid cell within HSUs

From a geometric perspective, an HSU can be regarded as a spatial polygon that signifies a landform entity with homogeneous terrain features in the field environment. At the regional scale, there is obvious heterogeneity in the topography and boundary characteristics among different HSUs (Wang and Zhang, 2021; Wang et al., 2019, 2023), resulting in the immaturity of hydrological analysis methods based on slope units. In contrast, hydrological analysis methods based on grid cells are well-developed. Some researchers have employed grid cells integrated with an infinite slope model or the limit equilibrium method to conduct regional landslide assessment or prediction (Gu et al., 2014; Wang et al., 2023; Zhang et al., 2021; Zhuang et al., 2016). In this study, each HSU was conceptualized as a composition of grid cells with similar microtopographic features, as illustrated in Fig. 2a. For each HSU, we utilized GIS spatial analysis tools to quantify the number of grid cells contained within it and their corresponding row and column positional information, thus establishing a comprehensive database that includes the position information of the grid cells within each HSU.

Figure 2 The diagram for HSU-grid cell hydrological connection: (a) grid cells contained within HSU (r_i, c_j represent the row and column of grid cells contained within HSU, respectively).

2.2.2 Initial water content assignment of HSUs

After obtaining the grid cell information contained within each HSU, conducting a rainfall infiltration analysis for these grid cells represents a complex and important task. One issue that cannot be overlooked is initial moisture content. For landslides in the Three Gorges Reservoir area of China, the soil typically experiences a prolonged dry winter before the rainy season (May to September). Previous research indicates that the residual moisture content of slopes before the rainy season averages approximately 7 % to 8 % (Wang et al., 2023). Accordingly, this study categorizes the initial water content into two components: the residual moisture content (θ_res) and the moisture content increment caused by antecedent precipitation (θ_ante). Here, θ_res reflects the average moisture level of the soil prior to the rainy season, while θ_ante indicate the increase in moisture content due to antecedent effective precipitation prior to landslide occurrence.

In this study, each grid cell is stratified into ten soil layers, each with a thickness of 0.2 m (Fig. 2b). For the Three Gorges Reservoir area, the regional landslides triggered by rainfall were mainly shallow (with thicknesses of 2–3 m). Therefore, variations in residual moisture content within the soil depth were disregarded, and the same residual moisture content value was assigned to each soil layer. Following this, we applied steady-state infiltration theory to simulate the distribution of moisture content across the soil layers influenced by antecedent precipitation, thereby allocating the antecedent rainfall to each soil layer. The calculation for θ_ini of each soil layer within the grid cell is as follows:

$$\begin{array}{}\text{(2)}& {\mathit{\theta }}_{\mathrm{ini}}\left(k\right)={\mathit{\theta }}_{\mathrm{ante}}\left(k\right)+{\mathit{\theta }}_{\mathrm{res}}(k=\mathrm{1},\mathrm{2},\mathrm{3}\mathrm{\dots },n)\end{array}$$

Where n represents the number of soil layers, and here n=10; θ_ini(k) indicates the initial moisture content of each soil layer; θ_ante(k) refers to the moisture change in each soil layer due to previous precipitation; θ_res stands for the residual moisture content in the grid cell.

2.2.3 Rainfall infiltration process simulation of grid cell

After obtaining the initial moisture content distribution, the 1-dimensional Richards infiltration equation was used to solve the moisture content distribution in the grid cell during the rainfall infiltration process.

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left[D\left(\mathit{\theta }\right)\cdot {\displaystyle \frac{\partial \mathit{\theta }}{\partial z}}\right]-{\displaystyle \frac{\partial K\left(\mathit{\theta }\right)}{\partial z}}\end{array}$$

Where D(θ) represents the value of soil water diffusivity under unsaturated conditions and has $D\left(\mathit{\theta }\right)=K\left(\mathit{\theta }\right)/\frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}{\mathit{\psi }}_m}$.

The finite difference scheme outlined above was formulated for numerical simulation of hydrological processes. The lower boundary, identified as impermeable, is based on the maximum soil depth of the grid cell. The upper boundary of each grid cell was designated as an infiltration boundary. When the rainfall intensity I(t) is less than the infiltration capacity of the topsoil, all precipitation infiltrates into the soil, and no runoff is generated. In this scenario, the infiltration boundary of precipitation was governed by the following differential equation:

$$\begin{array}{}\text{(4)}& -D\left(\mathit{\theta }\right){\displaystyle \frac{\partial \mathit{\theta }}{\partial z}}+K\left(\mathit{\theta }\right)=I\left(t\right)\end{array}$$

When the rainfall intensity exceeded the soil infiltration capacity, the excess portion was transformed into overland flow. At this point, the rainfall infiltration boundary was governed by the following equation:

$$\begin{array}{}\text{(5)}& \mathit{\theta }={\mathit{\theta }}_{\mathrm{s}}\end{array}$$

Where θ_s is the saturated water content of the grid cell.

2.2.4 Soil water content generation of HSU

Following the calculation of the soil moisture for individual grid cells, the soil water distribution of the HSU was computed as follows:

$$\begin{array}{}\text{(6)}& {\mathit{\theta }}_{\mathrm{HSU}}\left(k\right)={\displaystyle \frac{\sum _{k=\mathrm{1}}^n\mathit{\theta }\left(k\right)}{n}}\end{array}$$

Where θ_HSU(k) represents the moisture content of the kth layer of the HSU, θ(k) denotes the moisture content of the kth layer in the grid cell. n is the number of soil layers (n=10).

2.3 HSU_prob: the calculation of instability probability of HSUs

2.3.1 Profile extraction

After calculating the soil water content within each HSU, analyzing the stability of HSUs during the rainfall infiltration process is another important task. At present, the time required to carry out 3D analysis for each HSU on a large regional scale is too large, so extracting the calculation profile of the HSU becomes a reasonable selection. Currently, there is no uniform method for extracting the calculation profile of HSUs. Some reasonable assumptions are summarized as follows: the position of the profile line should reflect the elevation difference between the front and back edges of the slope, and the centroid point of the HSUs should be on the calculated profile to ensure that the soil weight on both sides of the calculated section is relatively uniform, and the areas of the two sections should be close to each other.

Based on these considerations, we developed a fast extraction algorithm HSU-profile (Wang and Zhang, 2021; Wang et al., 2023) for HSU profiles at large regional scales, which can be divided into three steps.

Figure 3 HSU instability probability calculation diagram. (a) Extraction of HSU boundary points; (b) profile line extraction of HSU polygon; (c) random search for critical slip surface; (d) enlarged view of the sliding mass for detailed visualization.

First, the highest elevation point H of the HSU polygon is connected to centroid point C to obtain line segment L₁, which intersects the HSU polygon at point J₁ (Fig. 3b). Line segment L₁ divides the HSU polygon into two parts, and the areas of the two parts, S₁ and S₂ are calculated to obtain the area ratio $A=S_{\mathrm{1}}/S_{\mathrm{2}}$.

Next, the lowest elevation point L and centroid C are connected to form line segment L₂, as shown in Fig. 3b. Determine the intersection point J₂ between L₂ and the polygon of the slope unit is determined. At this point, the HSU was divided into two parts by line segment L₂, and the areas of the two parts, S₃ and S₄, were calculated to obtain the area ratio $B=S_{\mathrm{3}}/S_{\mathrm{4}}$.

Finally, $\left|A\right|$ and $\left|B\right|$ are compared. A smaller absolute value of A indicates that line segment L₁ divides the areas on both sides of the HSU polygon more evenly. In this case, L₁ is selected as the profile line. Otherwise, the line segment L₂ was chosen as the profile line.

2.3.2 Calculation of safety factor F_s calculation

For each HSU, the Monte Carlo method was used to generate a large number of potential polyline-type slip surfaces (Fig. 3c), and the random walk method (Greco, 1996; Zhang and Zhang, 2006) was employed to search for the critical slip surface. The random search times for the sliding surface of each HSU are set to 500 times, then the infinite slope model was used to calculate the safety factor F_s of each potential slip surface as follows:

$$\begin{array}{}\text{(7)}& F_{\mathrm{s}}={\displaystyle \frac{\tan \mathit{\phi }}{\tan \mathit{\alpha }}}+{\displaystyle \frac{c+u_{\mathrm{s}}\tan \left({\mathit{\phi }}^b\right)}{{\mathit{\gamma }}_{\mathrm{s}}{D}_{\mathrm{s}}\cos \mathit{\alpha }\sin \mathit{\alpha }}}\end{array}$$

Where c is the effective cohesion of the soil, φ is the effective internal friction angle of the soil, and r_s is the average unit weight of soil above the slip surface. φ^b is related to the matric suction; when the matric suction is low, it is close to the internal friction angle (Zhang et al., 2018). D_s is the thickness of the soil layer above the slip surface. u_s represent the matric suction, which can be described by the Van Genuchten model(Van Genuchten, 1980):

$$\begin{array}{}\text{(8)}& S_{\mathrm{e}}={\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}={\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\left({\mathit{\alpha }}_{\mathrm{w}}\times u_{\mathrm{s}}\right)}^n}}\right]}^m\end{array}$$

Where S_e represents the saturation degree, θ denotes the soil water content of the HSU, θ_s and θ_r are the saturated and residual water content, respectively. The parameters α_w, n and m characterize the shape of the soil–water characteristic curve, with the relationship $n=\mathrm{1}-\mathrm{1}/m$.

2.3.3 HSU_prob calculation

According to the saturated-unsaturated rainfall infiltration theory, the mechanical parameters of the soil (such as cohesion force c (kPa) and internal friction angle φ (°)) are significantly affected by soil moisture content fluctuations. The variation in soil mechanical parameters during the process of rainfall infiltration is very complex, and it is generally acknowledged that dry soil prior to rainfall infiltration exhibits higher mechanical strength (characterized by elevated parameter values). As rainwater continues to infiltrate, the soil water content gradually increases, leading to a decreasing trend in mechanical parameters, such as cohesion and internal friction angle. Consequently, the mechanical parameters of the soil within each HSU are not fixed, but spatial uncertainty exists to some extent. In this context, employing probabilistic analysis methods to calculate the instability probability of an HSU is a more reasonable choice. Probability density functions (such as normal or uniform distributions) are commonly used to describe the uncertainty of the geotechnical parameters. The normal distribution is considered suitable for small areas or watersheds where hydrogeological parameters can be collected in detail, whereas a uniform distribution is more applicable for larger areas, where it is difficult to acquire detailed hydromechanical parameters (Wang and Zhang, 2021; Wang et al., 2023).

In this study, we utilized a uniform distribution to simulate the uncertainty of the mechanical parameters within the HSUs. The soil mechanical parameters in the unsaturated state before rainfall were taken as the upper bound, while those in the fully saturated state were considered the lower bound, thereby establishing the upper and lower value boundaries for the mechanical parameters within the HSU, as indicated in Eqs. (9) and (10):

$$\begin{array}{}\text{(9)}& {\displaystyle }c\in \left[c_{\mathrm{lower}},c_{\mathrm{upper}}\right]\text{(10)}& {\displaystyle }\mathit{\phi }\in \left[{\mathit{\phi }}_{\mathrm{lower}},{\mathit{\phi }}_{\mathrm{upper}}\right]\end{array}$$

where c_upper and c_lower represent the upper and lower bounds of c (kPa), respectively, φ_upper and φ_lower represent the upper and lower bounds of φ (°), respectively. The Monte Carlo method was employed to randomly select the values within these bounds. The instability probability of the HSU was calculated using Eq. (11).

$$\begin{array}{}\text{(11)}& {\mathrm{HSU}}_{\mathrm{prob}}={\displaystyle \frac{{\mathrm{Sum}}_{F_{\mathrm{s}}<\mathrm{1}}}{m}}\end{array}$$

where m represents the number of random selections for the mechanical parameters and m is set to 500.

2.4 The obtainment of key fitting parameters α and β for I–D curves of HSUs

In this study, an HSU is regarded as unstable when the value of HSU_prob exceeds 50 %. Then, the rainfall intensity and duration data with HSU instability under different rainfall scenarios were recorded to obtain the key fitting parameters α and β for the I–D curves of each HSU, thereby establishing a database of parameters α and β. The detailed steps are as follows.

2.4.1 Setting the antecedent effective rainfall levels AER_i (i=1, 2, 3 …, n)

The antecedent effective rainfall(AER) has a significant impact on landslide occurrence. Previous research indicates that in the Three Gorges Reservoir area, the minimum value of AER before landslide occurrence is 0 mm, whereas the maximum value of AER can exceed 170 mm (Wang and Zhang, 2021). Therefore, 20 different levels of AER ranging from 0 to 200 mm were established with intervals of 10 mm.

2.4.2 Design of the combination of rainfall intensity (I) and duration (D)

For each antecedent rainfall level, we categorized rainfall intensity (I) into eight levels to represent the variation from light to heavy rainstorms: 2, 5, 10, 20, 30, 40, 50, and 60 mm h⁻¹. The rainfall duration (D) ranged from 1 to 24 h, with intervals of one hour. Consequently, 192 combinations of I and D were generated for each AER level.

2.4.3 Generation of fitting parameters α, β of the I–D curves

For each combination of rainfall intensity and duration data, the method outlined in Sect. 2.2 is used to determine the soil water distribution within each HSU, and the corresponding value of HSU_prob was computed using the method described in Sect. 2.3. If the HSU is unstable, the corresponding intensity and duration data can serve as data points for fitting the I–D curves. Subsequently, a power function was utilized to fit these data points to obtain the key fitting parameters α and β of the I–D curve. As presented above, the fitting parameters α and β for the I–D curve of each HSU can be generated, thereby establishing a database for α and β at different AER levels.

2.5 Warning Mode

In practical applications, the antecedent effective rainfall(AER), rainfall intensity (I), and duration (D) for each HSU can be computed using Quantitative Precipitation Estimation (QPE) and Quantitative Precipitation Forecasting (QPF) products provided by the meteorological department (Wang and Zhang, 2021). Next, we analyzed the relationship between the actual value of AER and the 20 levels of AER documented in the database, thereby determining the level that is closest to the antecedent effective rainfall data of the HSU. The α and β values corresponding to this level were retrieved from the database for the following assessments.

1. If I≥αD^β, the data point (I, D) is above the warning curve; thus, the warning should be released.

2. Conversely, if I<αD^β, it signifies that the data point (I, D) is below the warning curve; therefore, no warning should be issued.

Figure 4 The flow chart of the fast warning mode based on parameter α and β database.

The programming languages Fortran 95 and Python 3.1 were employed to compile the algorithms outlined in Sect. 2.1–2.5, and the overall flowchart of the warning mode is depicted in Fig. 4.

3 Study area and data

3.1 Study area and slope unit data

Fengjie County is situated in the eastern region of the Three Gorges Reservoir area, with geographical coordinates ranging from 109°1^′17^′′ to 109°45^′58^′′ E and 30°29^′19^′′ to 31°22^′33^′′ N, covering a total area of 4087 km². The region has a subtropical humid monsoon climate with an annual average rainfall of 1500–2000 mm. The rainy season occurs from May to September, accounting for 70 % of the annual precipitation. The terrain is primarily mountainous and the Yangtze River flows across the region from west to east. Geological hazards, such as landslides, debris flows, and collapses, are widely distributed in Fengjie County, with rainfall landslides posing the most significant threat. Based on the 7 m DEM of Fengjie (Fig. 5a), the MIA-HSU method was employed to extract the slope units, resulting in the identification of 17 547 HSUs (Fig. 5b and c). Histograms of the slope gradient and area distribution of the HSUs are presented in Fig. 5d–e. As shown in Fig. 5d, the slope gradients of the HSUs follow a normal distribution, with 85.4 % of the slopes falling within the range of 10 to 30°. Figure 5e illustrates that the average area of the HSUs is 0.23 km², with 53.9 % of the slope units having an area less than 0.25 km². Because the sliding depth of shallow landslides typically ranges from 2 to 3 m, the majority of HSUs can be classified as small-to medium-scale landslides (with volumes under 500 000 m³).

Figure 5 Division of HSUs in Fengjie County. (a) Fengjie DEM; (b, c) extraction results for selected regions: enlarged view; (d) histogram of slope distribution of HSUs; (e) histogram of area distribution of HSUs.

3.2 Soil mechanical parameter c (kPa) and φ (°) data of HSUs

The rainfall-triggered shallow landslides within the study area are mainly composed of quaternary clay and silt, which are classified as fine-grained soils (Wang and Zhang, 2021; Wang et al., 2023). Field investigations indicate that the sliding soil is fully or even oversaturated, with some soil mass transitioning into mudflow during the sliding process. The laboratory moisture content tests revealed that the soil water content under these conditions approached or exceeded the liquid limit. To obtain detailed soil mechanical parameters under different moisture states, we conducted extensive field sampling across Fengjie County, resulting in 312 sampling points, as depicted in Fig. 6f. For each sampling point, direct laboratory shear tests were performed to derive the soil mechanical parameters c (kPa) and φ (°) at the liquid and plastic limits, respectively. Based on geological survey data provided by the Fengjie County Land Bureau, the dry density of soil within a 10 m thickness ranges from 1.7 to 1.8 g cm⁻³ (Wang et al.,2021). Therefore, the dry density of each soil samples is randomly selected within this range. In accordance with the ASTM-d6528 (ASTM D6528-17, 2017) standard, 312 groups of liquid-plastic limit tests and 624 groups of undrained direct shear tests were performed to obtain the mechanical parameters of each sample at both the liquid and plastic limit water contents. Subsequently, ArcGIS spatial analysis tools were utilized to generate distribution maps of c (kPa) and φ (°) under plastic and liquid limit moisture conditions, as shown in Fig. 6g–j.

Figure 6 State of Landslide Soil Before and After Rainfall. (a) Soil approaching plastic limit moisture content before rainfall; (b) soil approaching plastic limit moisture content before rainfall; (c) soil in a fluid state after rainfall; (d) soil in a fluid state after rainfall; (e) fully saturated and liquefied soil after rainfall; (f) soil sampling locations; (g) c (kPa) at plastic limit moisture content; (h) φ (°) at plastic limit moisture content; (i) c (kPa) at liquid limit moisture content; (j) φ (°) at liquid limit moisture content.

3.3 Rainfall data

Rainfall data sources include Quantitative Precipitation Forecasting (QPF) products and Quantitative Precipitation Estimation (QPE) products. The QPF product obtained from the local government of Fengjie County is typically utilized to forecast future rainfall at a regional scale, which can provide rainfall forecast products for the next hour. QPE data are applied to estimate historical regional rainfall at a regional scale and are essential for determining the antecedent effective rainfall (AER), which can be computed as follows:

$$\begin{array}{}\text{(12)}& \mathrm{AER}=\sum _{i=\mathrm{1}}^n{a}^n{R}_i\end{array}$$

Where AER is the antecedent effective rainfall, a is the attenuation coefficient, which is equal to 0.84, based on the research of the Fengjie count (Wang and Zhang, 2021), n is the number of days before the landslide occurs.

4 Case study: Rainfall-induced landslides of 31 August 2014

From 30–31 August 2014, Fengjie experienced continuous heavy rainfall, triggering a series of landslide hazards that resulted in over 30 fatalities and an economic loss of CNY 580 million. Based on the QPE data with the resolution of 1 km for 15–31 August, the effective precipitation for the 15 d prior to the landslide hazards is shown in Fig. 7. As illustrated in Fig. 7, the maximum precipitation during this period was 179.10 mm, which occurred in the northwestern region of the area. The hourly QPF data for 31 August (with a resolution of 1 km) are presented in Fig. 8a–l. As illustrated in Fig. 8a–d, the rainfall was minimal from 00:00 to 08:00 Beijing Standard Time (UTC+8), with a maximum cumulative rainfall of 12.2 mm 08:00 BST. As shown in Fig. 8e–g, rainfall began to increase rapidly at 10:00 BST, reaching a maximum cumulative precipitation of 92.40 mm by 14:00 BST in the northwestern region of Fengjie County. Figure 8h–l indicate that from 16:00 to 24:00 BST, the cumulative rainfall remained constant, suggesting that the rainfall process had ceased.

Figure 7 Precipitation Data Processing (Effective Precipitation from 15 to 30 August 2014).

Figure 8 Radar forecast precipitation data for 31 August 2014. (a) 02:00 BST; (b) 04:00 BST; (c) 06:00 BST; (d) 08:00 BST; (e) 10:00 BST; (f) 12:00 BST; (g) 14:00 BST; (h) 16:00 BST; (i) 18:00 BST; (j) 20:00 BST; (k) 22:00 BST; (l) 24:00 BST.

The Land and Resources Bureau of Fengjie County provided data on landslide points triggered by rainfall on 31 August. This heavy precipitation triggered 583 landslides, which were mainly distributed in the northwestern region (as indicated by the red and green solid points in Fig. 9). This study utilized the QPE (Fig. 7) and QPF data (Fig. 8a–l) as inputs to forecast landslide hazards for 31 August.

Figure 9 Prediction results at 02:00 to 24:00 BST. (a) 02:00 BST; (b) 04:00 BST; (c) 06:00 BST; (d) 08:00 BST; (e) 10:00 BST; (f) 12:00 BST; (g) 14:00 BST; (h) 16:00 BST; (i) 18:00 BST; (j) 20:00 BST; (k) 22:00 BST; (l) 24:00 BST.

The landslide forecast results from 02:00 to 24:00 BST are shown in Fig. 9a–l. It can be seen from Figs. 8 and 9 that there is a good correlation between the spatial distribution of unstable HSUs and rainfall characteristics. As presented in Fig. 9a–d, at the beginning of the rainfall process (before 08:00 BST), the majority of the HSUs remained stable owing to the minimal rainfall. Unstable HSUs began to emerge in the northwestern region starting at 10:00 BST, coinciding with the rapid increase in rainfall. Additionally, as the rainfall progressed, the number of unstable HSUs increased swiftly and spread towards the central and southern regions (Fig. 9f–g). Notably, many unstable slope units appeared within several hours after heavy rainfall ceased. Figure 9h–l reveal that from 16:00 to 24:00 BST, although the heavy rainfall essentially ended, the number of unstable HSUs continued to rise because of the moisture infiltration of the saturated top soil, reaching a total of 3987 at 24:00 BST.

This study employs the Receiver Operating Characteristic (ROC) method to analyze the predictive performance of the HSU (Fawcett, 2006). For physically model-based slope units, the ROC method describes the following four possible states using a contingency table:

1. True Positive (TP): HSU contain landslide points and exhibit instability;

2. True Negative (TN): HSU does not contain landslide points and does not exhibit instability;

3. False Positive (FP): HSU does not contain landslide points but exhibits instability;

4. False Negative (FN): HSU contains landslide points but does not exhibit instability.

Table 1 Analysis of Forecast Results for the 831 Case Study.

According to GIS spatial statistics, 583 landslides triggered by rainfall on 31 August were contained within 425 HSUs. In this study, these HSUs are taken as benchmark values for the calculation of TP, TN, FP, and FN, and the missing alarm rate (MAR) and false alarm rate (FPR) can be calculated as follows:

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathrm{MAR}=\mathrm{100}\mathit{\%}\times \mathrm{FN}/\mathrm{425}\text{(14)}& {\displaystyle }\mathrm{FPR}=\mathrm{100}\mathit{\%}\times \mathrm{FP}/(\mathrm{FP}+\mathrm{TN})\end{array}$$

The detailed forecast results for 2–24 h are shown in Table 1. As shown in columns 7–8 of Table 1, with increasing rainfall duration, the MAR gradually decreases, while the FPR gradually increases. Taking the result of 24 h as an example, the MAR of 24 h is 11.8 % and the FPR is 21.1 %, indicating that the prediction result can satisfy the requirement of early warning practice.

According to the ROC method, the precision and accuracy of the prediction results were calculated as follows:

$$\begin{array}{}\text{(15)}& {\displaystyle }\mathrm{Precision}=\mathrm{TPR}/(\mathrm{TPR}+\mathrm{FPR})\text{(16)}& {\displaystyle }\mathrm{Accurancy}=(\mathrm{TP}+\mathrm{TN})/(\mathrm{TP}+\mathrm{FN}+\mathrm{TN}+\mathrm{FP})\end{array}$$

Table 2 provide the calculation results of precision and accuracy at 24 h. As shown, the precision of the forecasting results is 80.7 %, with an accurancy value of 79.1 %, indicating the proposed warning mode has satisfactory comprehensive forecasting performance.

Table 2 Calculation Results of Precision and Accuracy at the 24th Hour.

5 Discussion

5.1 The discussion on the computational efficiency

For emergency warnings during the rainy season, the swift release of warning information is crucial for local authorities to develop emergency plans and to evacuate residents from landslide-prone areas. Therefore, local governments not only seek satisfactory accuracy in the warning model but also require minimal time. To evaluate the computational efficiency of the proposed model, a standard laptop was utilized to execute the forecast for landslide hazards on 31 August. The device specifications and computation times are presented in Table 3. As shown in Table 3, for the regional scale covering several thousand square kilometers, the prediction model can rapidly complete real-time warnings for the next 24 h within 12 min, indicating that its computational efficiency can satisfy the requirements of emergency warning.

Table 3 Analysis of computational efficiency of the prediction model.

Table 4 The selected HSUs that report false alarms at 24:00 BST.

5.2 Further analysis of prediction performance

Using the 24 h prediction results as an example, we randomly selected seven HSUs with false alarms for further analysis (Table 4). Given that existing research indicates that antecedent effective rainfall (AER) and daily accumulated rainfall are crucial sensitive factors affecting the forecasting performance of HSU-based physical models (Wang et al., 2023), the AER of these HSUs, the AER levels assigned by the database and daily accumulated rainfall are presented in Table 4. As shown in Column 5 of this table, the relative error ranges from 0.7 % to 6.3 %, indicating that the 20 levels of the AER designed in the database can accurately reflect the effective antecedent rainfall characteristics of the HSUs. The average rainfall intensity, duration, and cumulative rainfall data at 24:00 BST are shown in Columns 6–8. As seen in Column 6, the cumulative rainfall for the seven HSUs ranges from 12 to 29.8 mm, with average rainfall intensities range from 0.5 to 1.25 mm h⁻¹, which can be classified as light to moderate rain type. The instability probability (HSU_prob) of these HSUs was calculated to investigate the causes of false positives. As shown in Column 9, among the seven HSUs with false alarms, five had an instability probability of less than 50 %, indicating that these HSUs did not experience instability during the rainfall process. Therefore, we cautiously conclude that although the prediction model exhibits preferable operational efficiency, it may increase the false positive rate to some extent.

To investigate the potential for reducing the false-alarm rate, we selected four HSUs from Table 4 for further analysis. Figure 10a–d presented the I–D curves and cumulative precipitation distribution histograms for these HSUs.

Figure 10 The I–D Curves of HSUs before the adjustment of rainfall process. (a) 8535; (b) 8561; (c) 6066; (d) 13 108.

For each HSU, the QPF data from 00:00 to 24:00 BST were discretized into 12 sets of rainfall intensity and duration data points at 2 h intervals (represented by black and red solid dots). The black solid dots positioned below the I–D curve indicate that the HSU is stable at that moment, whereas the red solid dots located above the curve signify false alarms at the current forecasting hour. As shown in Fig. 10a–d, the red false alarm points for the four HSUs are all situated very close to the I–D curve, nearly tangent to it. This proximity suggests that slight spatial adjustments to these points could alter the forecast results. Another important issue is that some of the black solid dots correspond to a cumulative rainfall of 0 mm, indicating that the rainfall process had not yet begun. Therefore, it is necessary to adjust the spatial positions of data points I and D based on the actual initiation time of the rainfall process, thereby facilitating an in-depth investigation of the causes of the false alarms.

Table 5 Rainfall process adjustment for HSU with number of 8535. Bold values indicate the recalculated rainfall intensity and duration values after adjustment.

In this study, an HSU with number 8535 is taken as an example to illustrate the process of adjusting the spatial positions of data points I and D. As shown in Table 5, the rainfall process for this HSU started at 12:00 BST and ended at 24:00 BST with a duration of 12 h. The start time of rainfall was used as the starting point to recalculate the rainfall intensity during the rainfall process, as shown in the bold text in Table 5. The adjusted average rainfall intensity was significantly higher than the values prior to adjustment. This means that the adjustment of the rainfall process led to notable changes in the spatial locations of the data points I and D. As shown in Fig. 11a, after updating the positions of data points I and D, the HSU does not exhibit any false alarms. Figure 11b–d present the updated forecast results for the other three HSUs after the adjustment. As shown in Fig. 11a–d, following the adjustments, three out of these four HSUs were able to release accurate results. Therefore, we advise that practical warning applications should consider the influence of the difference in rainfall processes of HSUs on the prediction results.

Figure 11 The I–D Curves of HSUs after the adjustment of rainfall process. (a) 8535; (b) 8561; (c) 6066; (d) 13 108.

5.3 Comparison with existing approaches in the Fengjie county

Statistical methods have been widely used to establish rainfall thresholds for landslides. Liu et al. (2024) categorized Chongqing into four distinct subregions based on topography, tectonic structures, lithological stratigraphy, and the distribution and characteristics of landslides, and established the rainfall thresholds for landslides in each subregion, as shown in Table 6 (Liu et al., 2024). Fengjie county belongs to Subregion II, with its rainfall threshold is listed in the second row of Table 6. In this table, R represents the cumulative daily rainfall, while R_e signifies the antecedent effective rainfall.

Table 6 Equations of the thresholds defined for each subregion (Liu et al., 2024).

Table 7 Prediction performance parameters of Liu's model.

The prediction results derived from statistical methods are illustrated in Fig. 12, with detailed prediction performance parameters provided in Table 7. As depicted in Fig. 12, there are 5462 unstable HSUs, primarily located in the northwest and southern mountainous regions. According to Table 7, 256 landslide events are successfully predicted by Liu's model. However, the statistical method exhibits a MAR of 40 % and a false alarm rate of 30.4 %, resulting in low overall precision and accuracy. In contrast, the missing and false alarm rate of the HSU_ID model are 11.8 % and 21.1 % at the 24th hour (Table 1), thus demonstrating superior forecasting performance.

Figure 12 Prediction outcomes generated by Liu's model.

5.4 Limitations and future work

The method proposed in this paper reduces dependence on landslide inventory data and avoids the issue of excessive computation time associated with physical models. This offers promising application prospects for emergency early warning in regions of third-world countries where landslide inventory data are scarce. However, it should be noted that the calculation of HSU_prob relies on regional-scale mechanical parameters. Therefore, establishing the distribution of mechanical parameters in Southwest China through field sampling, experiments, and spatial analysis is an important task for the future. Another important task lies in optimizing the model algorithm, because the physical mechanisms of rainfall landslides are extremely complex, and the soil-water coupling process under rainfall involves significant nonlinearity (Apip et al., 2010). Some machine learning methods that incorporate physical frameworks are expected to accurately describe this issue (for example the PINN). Therefore, improving the model algorithm by combining machine learning with physical approaches is another important task for the future.

6 Conclusion

Currently, the operational forecasting of rainfall-induced landslides over regional scales of thousands of square kilometers faces significant challenges. Conventional physical and statistical approaches have shown limitations in terms of achieving satisfactory results. This study utilized HSU as a basis to integrate physical models and rainfall threshold methods for a warning model applicable to large-scale regions. The warning model employs HSU as a prediction unit to improve the clarity of the warning results, physical methods are utilized to develop the warning criteria, thereby reducing the overreliance on historical observational data, and a database of rainfall parameters across different rainfall scenarios is constructed, which enhances the efficiency and applicability of the warning model. The prediction performance was validated through a case study of “8.31” rainfall landslides in Fengjie County. The conclusions are as follows.

1. A rainfall-triggered landslide warning model was established by integrating HSUs, physical approaches, and rainfall parameters.  Initially, a grid-based HSU hydrological analysis technique was established to determine the soil moisture content distribution within the HSUs during different rainfall hours. Subsequently, computer graphics algorithms, random search techniques, and infinite slope models were used to develop a regional-scale HSU stability analysis method. Soil mechanics parameters at the limit of water content and probability density functions were used to describe the spatial uncertainty of the soil mechanical parameters within the HSU during rainfall infiltration, allowing for the calculation of the instability probability of the HSU.  Different rainfall scenarios were simulated to derive rainfall intensity I and duration D data that can trigger HSU instability, thereby constructing early warning curves for the rainfall thresholds of the HSU.

2. A database for the I–D curve fitting parameters α and β across various AER levels was established. This database includes α and β data for 17 547 HSUs across 20 AER levels, amounting to a total of 350 940 records, thus offering substantial data support for rainfall-induced landslide predictions in Fengjie County. In practical applications, it is sufficient to quickly issue warning information by assessing the relationship between the values of I and αD^β, thereby reducing the time required to calculate the safety factors using conventional physical models. The calculation efficiency test indicates that the warning mode can perform forecasts for thousands of kilometers within a runtime of less than 12 min, thereby meeting the operational needs for real-time warnings over large regional scales.

3. The case study indicates that the distribution trends of unstable HSUs align well with rainfall characteristics. As the rainfall duration increased, the MAR gradually decreased, while the FPR continued to increase. Taking the 24 h forecast results as an example, the MAR was 11.8 %, while the false alarm rate was 21.1 %. ROC analysis revealed that the accuracy of the forecast result at this moment was 80.7 %, with a precision of 79.1 %, reflecting satisfactory overall forecasting performance. Further discussion of the false alarm rate suggests that adjusting the spatial locations of rainfall intensity and duration data points based on the rainfall characteristics of each HSU may be conducive to reducing false alarm rates.