*4.1. Basic Characteristics*

This rainfall event triggered 667 landslides over an area of 108 km2, and the majority of them (552 landslides) occurred in the Mibei and Yanhua villages (Figure 11a,b). The types of landslides were various, including shallow landslides combined with a small number of flowslides, rockfalls, and debris flows with a total landslide area of 0.75 km2. The largest landslide area was approximately 20,000 m2, the smallest area was 50 m2, and the average landslide area was about 1100 m2. According to the statistics, there were 288 landslides with an area of less than 500 m2, accounting for approximately 43% of all landslides. There were 291 landslides with an area of 500~2000 m2, accounting for approximately 44%. The number of landslides with an area of 2000~5000 m2 and > 5000 m<sup>2</sup> was 71 and 17, respectively (Figure 10).

**Figure 11.** (**a**) Rainfall-induced landslide inventory of the 2019 heavy rainfall event in Longchuan County, the landslide abundance area of this event; (**b**) map showing the zooming of the landslide abundance area.

We calculated the landslide number density (LND) and landslide area density (LAD) within a 1.5 km-radius moving window using a Gaussian density kernel function. The LND and LAD maps indicated that the maximum LAD and LND of the study area were 9.5% and 78/km2, respectively (Figure 11). Landslides had obvious cluster distribution characteristics, and a large number of landslides were concentrated within 2 km of the Mibei Village (Figure 12).

#### *4.2. Factor Analysis*

In order to analyze the relationship between different influencing factors and the occurrence of landslides, we calculated the frequency distribution of landslides and landscape (i.e., non-landslide area) and the LAD of different influencing factors. Figure 13 shows the frequency density distribution of landslide and non-landslide areas, and Figure 14 shows the LAD in different intervals of six influencing factors (the higher the LAD, the more likely the landsliding will occur). For elevation, the frequency density distribution of landslide area and non-landslide area was basically the same (Figure 13a), the peak LAD was situated at elevations from 300 to 450 m, indicating that landslides were more likely to occur within this elevation range (Figure 14a). For the slope angle, the landscape area was clustered between 5 and 20◦, while most of the landslides occurred on slopes with the inclination of 15–35◦ (Figure 13b). Overall, the LAD increased with the increase in the slope angle and was described by an exponential relationship of y = *e*(0.21+0.08*x*), (where *x* is the slope angle and y is the LAD, Figure 14b), suggesting that with the increase in the slope angle, the landslide occurrence possibility also increased. In terms of topographic relief, the relief of non-landslide area was primarily concentrated in the 200~250 m range, whereas the landslide area was primarily clustered in the range of 250~350 m (Figure 13c). Overall, there was a negative logarithmic relationship between the LAD and relief, indicating that the LAD decreased with the increase in relief (Figure 14c). On the part of TWI, landslides were most commonly seen in the range with TWI values between

4 and 6, and there was a positive exponential relationship between the LAD and TWI of y = *e*(−12.16+1.18*x*) (where *x* is the TWI and y is the LAD), and the LAD increased as the TWI increased (Figure 14d). For road density, landslides were primarily distributed in the road density interval between 2 and 4 (Figure 13d). In general, LAD and road density had a linear relationship of y = 0.14*x* + 0.45 (where *x* is the road density and y is the LAD), which shows that landslides were more likely to occur in areas with a high road density (Figure 14e). For the distance to river, landslides were more likely to occur in the range of 100~400 m, and there was no obvious correlation between the LAD and the river distance (Figure 14f).

**Figure 12.** Map showing spatial density of landslides triggered by this rainfall event. (**a**) landslide number density (LND); (**b**); landslide areal density (LAD).

**Figure 13.** Frequency density estimates of landslides and landscape area for different influencing factors; (**a**) slope angle; (**b**) elevation; (**c**) relief; (**d**) road density; (**e**) distance to rivers; (**f**) TWI.

**Figure 14.** Map showing the relationship between the six influencing factors and the landslide areal density (LAD); (**a**) elevation; (**b**) slope angle; (**c**) topographic relief; (**d**) TWI; (**e**) road density; (**f**) distance to rivers.

Figure 15 shows the statistical results of the landslides and the slope aspect. Figure 15a shows the frequency density of the landslides and landscape (i.e., non-landslide area) on different slope aspects. The result demonstrates that the non-landslide area was evenly distributed in all aspects, but most of the landslide area was concentrated in the aspect of 110◦~180◦ (SE to S). The statistical results of LAD show that the peak LAD of 1.4% was present at the aspects from SE to S for the landslides.

Figure 16 shows the distribution of the landslide and non-landslide area, and the average landslide area in each land use unit. The result shows that the predominant land type was forest, which accounts for 80% of the study area, followed by cropland land, which accounts for more than 10%. The area of urban area and bare land was less than 1%. Among all land types, shrubland was the most prone to landslides, with roughly 10% of landslides occurring in the 5% area. Landslides were the least developed in cropland, maybe due to the relatively gentle slope of this unit. Furthermore, statistics on the average landslide area of different units suggest that bare land had the largest average landslide area, with more than 1600 m2, followed by forest land, which had an average landslide area of 1200 m2, and cropland had a relatively small average landslide area, only 600 m2.

**Figure 15.** (**a**) The distribution of aspect within landslides and landscape; (**b**) correlations between aspect and landslide area density (LAD).

**Figure 16.** Areal coverage (%) of different land use types for both landslide and landscape overlaid by average landslide area calculated per each unit.

#### **5. Physically Based Landslide Susceptibility Assessment**

#### *5.1. Brief Description of MAT.TRIGRS(V1.0)*

To address the issues of the manual modification of plentiful model parameters and complex data processing in the traditional TRIGRS model, Ma, et al. [72] proposed a new TRIGRS model using Matlab®programming. It can directly read the grid data of TIF format as the input, and then directly exports the prediction results of grid files, which greatly simplifies data preparation and parameter setting. It includes the script files INPUT DATA.m and TRIGRS.m. The INPUT DATA.m file is used to read the TIF input files, and TRIGRS.m is the executable program that can be used to calculate the pressure head and Fs. The minimum Fs and the corresponding pressure head are generated in the TIF format by calculating the pressure head and Fs at various soil depths. More description can be obtained in [72].

In the physically based model, in order to obtain accurate landslide prediction results, sufficient and accurate input data are required [68,73–75]. For the soil thickness distribution, the Z-model developed by Saulnier, et al. [76] was used to evaluate the soil thickness. We assumed that the maximum thickness of the soil in the study area was 5 m and the minimum thickness was 0.5 m based on previous studies [17,77]. Soil thickness can be estimated and calculated by Equation (8). The bedrock in the study area is monzogranite (O3-S1), and the landslide occurred primarily in the weathered soil layer on the bedrock's surface. The soil type of the weathered soil layer is sandy clay loam. Therefore, combined with previous studies [17,52,78], we assigned the corresponding values to mechanical and hydrological parameters including cohesion, internal friction angle, and soil weight of this

soil type (Figure 17 and Table 1). Based on previous experience [57,79], saturated hydraulic diffusivity D0 was set to *D*<sup>0</sup> = 200*Ks* and the initial surface flux (*IZLT*) was generally less than the Ks to one power or more and was set to *IZLT* = 0.01*Ks*.

$$h\_i = h\_{\max} - \left(\frac{Z\_i - Z\_{\min}}{Z\_{\max} - Z\_{\min}}\right) (h\_{\max} - h\_{\min}) \tag{8}$$

**Figure 17.** Maps showing the distribution of slope angle (**a**), soil thickness (**b**) and flow direction (**c**).

**Table 1.** Mechanical properties of the soil.


Simultaneously, in order to account for the uncertainties in the physical process that lead to slope failure, the Monte Carlo simulation, which is a robust and well-known approach in applications concerning probability analyses and reliability studies, was used in this study [56,80]. We considered the uncertainties of two main parameters (cohesion and internal friction angle) that primarily influence the slope failure. To characterize the probability density function (PDF) of the two random variables, the normal PDF was adopted. We assumed that the average and standard deviation of the cohesion were 29 kPa and 9 kPa and those of the internal friction angle were 20◦ and 6◦. Based on the Monte Carlo simulation, the input data were calculated by the TRIGRS model, yielding 1000 predicted pictures of potential landslides in the study area. Finally, the slope failure probability (Pf) of the study area was obtained.

#### *5.2. Landslide Susceptibility Assessment*

Figure 18 shows the distribution of the average value of 1000 predicted pictures calculated by rainfall data over different time periods. From the calculation results, we can observe that the Fs of all raster cells were greater than that before the rainfall event, indicating that all slopes were stable (Figure 17a). In addition, after 12 h of the rainfall (at 8:00 on 10 June 2019) 12-h rainfall reached 86 mm), the Fs of some grid cells in the study area decreased. Particularly, some grid cells with a large slope angle began to fail (Figure 17b). Then, although continuous rainfall occurred in the subsequent stage (after 11 June 2019), the change of Fs in the study area was relatively small, and few new grid units became unstable.

**Figure 18.** Slope stability conditions, expressed in terms of Safety of Factor (FS) in different time periods of this rainfall event; (**a**) 20:00 on 9 June 2019 (UTC + 8, before rainfall event); (**b**) 8:00 on 10 June 2019 (UTC + 8); (**c**) 20:00 on 11 June 2019(UTC + 8); (**d**) 8:00 on 13 June 2019(UTC + 8).

We calculated the Fs results in the various slope interval over different time periods (Figure 19). The result shows that the Fs of the grids with slope angles between 30 and 40◦ was mostly distributed between 1.3 and 2.5, with an average value of around 1.6. After the onset of heavy rain on 10 June 2019, the Fs of raster cells rapidly decreased, and the Fs of most grids ranged between 0.9 and 1.7, with an average value of about 1.2. From 8:00 on 10 June 2019, although there was rainfall every day at a subsequent stage, the average rainfall was less than 2 mm/h. The low rainfall intensity had a little impact on the slope stability. Rainfall increased to some extent on 12 June 2019, reaching 45 mm in 12 h, and the Fs decreased slightly. For grids with a slope larger than 40◦, we also found the same trend that the Fs of most grid units decreased rapidly after heavy rainfall, and then basically remained unchanged. Overall, the Fs of grids with a slope greater than 40 degrees was much smaller than grids with a slope between 30 and 40 degrees.

**Figure 19.** The statistical results of Factor of safety (Fs) in the various slope interval at different rainfall times; (**a**) slope angle: 30~40◦; (**b**) slope angle: >40◦.

Figure 20 shows the probability distribution of slope failure in different time periods. Obviously, the prediction results of Pf were roughly consistent with the actual landslide distribution. Most areas with a high probability (blue areas) were located on both sides of the river valley, that is, the areas with relatively steep slopes. Before rainfall, almost all the grids in the study area were less than 0.1, indicating that the slope before rainfall was stable. After 12 h of rainfall (at 8:00 on 10 June 2019), the area with steep slopes began to show the instability phenomenon, and the Pf of some grids reached more than 0.6. In the following continuous rainfall, with the decrease in rainfall intensity, there was a slight increase in the area with a high probability of failure.

**Figure 20.** Probability of slope failure (Pf) in different time periods of this rainfall event; (**a**) before rainfall event; (**b**) 8:00 on 10 June 2019; (**c**) 20:00 on 11 June 2019; (**d**) 8:00 on 13 June 2019.

To quantitatively analyze the susceptibility results, we counted the class area, landslide area, and the corresponding LAD of different susceptibility classes before and after rainfall (Figure 21). Based on the natural breaks, the susceptibility level was divided into four classes (i.e., very low, low, moderate, and high). The result shows that before the occurrence of rainfall, most areas belonged to the low susceptibility area, and the majority of landslides were concentrated in very low and low susceptibility areas. With the occurrence of rainfall, the area of low susceptibility areas decreased, while the area of high susceptibility areas increased. The statistical result reveals that 12.1% of the total landslides occurred in the 25.0% of the area which were classified as moderate and high. Meanwhile, the LAD increased with the increase in the susceptibility level, which also shows that the model can effectively predict the potential landslide-prone zone.

**Figure 21.** Susceptibility class distribution and the occurrence of landslides within the study area; (**a**) before rainfall; (**b**) after rainfall.

#### **6. Discussion**

China's southeast area is situated in a subtropical monsoon climate zone with frequent typhoons and rainstorms. The most common types of geological hazards in this area are landslides and debris flows caused by rainfall, which have the characteristics of a small scale of individual hazard point, a large number of groups, and a wide distribution range [18,19]. In mountainous areas, the effect of the orographic amplification of rainfall and the projection of rainfall-vector on hillslopes [81,82] might result in the windward hill-slope receiving more rainfall, leading to more landslides on the hillslope scale [83]. Due to the influence of the monsoon depression and tropical cyclone, the southeast monsoon prevails in the Longchuan area during the summer (June and July). The landslide distribution of this rainfall event indicates that the southeast and south aspect hillslopes are more prone to collapse than the northwest-north aspect ones (Figure 14). The main reason for this phenomenon is that the south slope is mostly windward, which causes more rainfall and splash erosion in the area. Otherwise, the bedrock weathering degree of the south slope will also be high due to the influence of environmental factors such as soil moisture content, surface temperature, light time, and so forth, leading to relatively weak mechanical parameters of rock and soil mass. Therefore, under the condition of heavy rainfall, the south slope is more prone to landsliding.

Slope angle is an important topographic factor affecting the occurrence of landslides. From the spatial distribution of the landslides, we can observe that the landslides were mainly distributed in low mountainous areas, with the sections at elevations within 300~450 m and slopes ranging from 15 to 35◦ (Figure 12). The LAD increased with the increase in slope angle and was described by an exponential relationship, indicating that the landslides of this event more easily occurred in areas with steep slopes (Figure 13b). TWI reflects how surface morphology affects soil groundwater level and moisture content, which is represented by a theoretical measure of the accumulation of flow [84,85]. According to the statistical results, there was an exponential relationship between the LAD and TWI, and the LAD increased as the TWI value increased. Especially when the TWI was greater than 10, the LAD increased rapidly (Figure 14f). Higher soil moisture causes higher pore water pressure and reduces the strength of rock and soil mass. As a result, when it rains, the pore water pressure in these areas rises rapidly, resulting in slope failure.

Anthropogenic factors (such as land-use change, deforestation, hill cutting, etc.) play a significant role in the initiation of landslides in active mountain ranges [86–88]. The construction of roads has significantly altered the slope stability of mountainous areas, making them prone to landslides. When a road is built, the toe of the slope is excavated or the weight of the slope is increased, and the overall stability above the slope is reduced, resulting in the occurrence of new landslides or the reactivation of old landslides[89]. From Figure 13d, we can observe that the landslides of this event were more likely to occur in areas with high road density, illustrating that anthropogenic factors have accelerated the instability of the slope in this area. Furthermore, in the Longchuan area, the majority of local residents have excavated mountains to build houses, leading to a number of nearly vertical artificial slopes. Meanwhile, human activities will fragment surrounding natural slopes and increase the degree of rock weathering, which will also exacerbate slope instability in mountain areas.

The formation lithology of the slope is the material basis of landslides. Granite layers are one of the most common strata in China's southeast coastal regions. Long-term weathering of granite results in widely distributed residual soil layers. For the Longchuan area, the bedrock is monzogranite (O3-S1), and the landslides occurred primarily in the weathered soil layer on the bedrock surface [17,77]. The major influence depth of heavy rainfall was limited to the superficial zone of slopes due to the difference in rainfall intensity and permeability of granite residual soil. This is why the shallow surface zone was severely affected by landslides [17]. A saturated seepage field was formed in the shallow surface zone of slopes as a result of prolonged heavy rainfall. The mechanical strength of saturated soil diminished, and slide failure occurred at the shallow surface saturation zone.

#### **7. Conclusions**

In this work, we established a landslide inventory including all the landslides induced by the 2019 Longchuan heavy rainfall event in Guangdong Province, China. We described the topographical, geological, and hydrological control of landslide hazards. Furthermore, we conducted the physically based susceptibility assessment of shallow landslides based on the MAT.TRIGRS (V1.0) tool. The following conclusions can be drawn: (1) This rainfall event triggered about 670 landslides with a total area of 0.75 km2; the landslides had obvious cluster distribution characteristics, and a large number of landslides were concentrated within 2 km of the Mibei village. (2) The landslide abundance was closely related to slope angle, TWI, and road density but had a low correlation with elevation and distance to rivers. Among them, the LAD increased with the increase in the slope angle and TWI and was described by an exponential relationship. Otherwise, the statistical results of the landslides and the slope aspect showed that most of the landslide area was concentrated in the aspect of 110◦~180◦ (SE to S). (3) The physically based susceptibility assessment results indicated that the prediction results were roughly consistent with the actual landslide distribution, and most areas with a high susceptibility were located on both sides of the river valley. The onset of heavy rain on 10 June 2019 was the main triggering factor of this group-occurring landslides. Our study will be beneficial for understanding the distribution pattern and cause of rainfall-induced shallow landslides in the Longchuan area, and it can provide data and technical support for the prevention of rainfall-induced geological disasters in the southeast mountainous area of China.

**Author Contributions:** Conceptualization, C.X. and S.M.; methodology, S.M.; software, S.M. and X.S.; validation, S.M. and X.S.; formal analysis, S.M. and X.S.; data curation, S.M. and X.S.; writing original draft preparation, S.M.; writing—review and editing, C.X. and S.M.; supervision, C.X.; project administration, C.X.; funding acquisition, C.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the National Institute of Natural Hazards, Ministry of Emergency Management of China (ZDJ2020-10 and ZDJ2021-12).

**Acknowledgments:** We thank Google Earth and Sentinel-2 satellite images for the free access satellite images used in this study. We also express our appreciation of the constructive comments provided by the anonymous reviewers, which is beneficial for the quality of the manuscript a lot.

**Conflicts of Interest:** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence he work reported in this paper.
