**1. Introduction**

In recent years, an increasing number of intense rainfall events have occurred in mountainous areas due to the impact of global climate change, which has dramatically increased the frequency of global rainfall-induced landslides [1,2]. Rainfall-induced landslides are not only widely distributed in the world but also occur frequently and cause significant

**Citation:** Ma, S.; Shao, X.; Xu, C. Characterizing the Distribution Pattern and a Physically Based Susceptibility Assessment of Shallow Landslides Triggered by the 2019 Heavy Rainfall Event in Longchuan County, Guangdong Province, China. *Remote Sens.* **2022**, *14*, 4257. https:// doi.org/10.3390/rs14174257

Academic Editors: Paolo Mazzanti and Saverio Romeo

Received: 29 July 2022 Accepted: 24 August 2022 Published: 29 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

damage to humanity [3–6]. Therefore, a good understanding of the fundamentals of rainfallinduced landslide occurrence, distribution patterns, and susceptibility assessments can provide useful guidance for regional disaster prevention and mitigation, and landscape evolution [7–9].

A new landslide inventory that is generated after a major triggering event (e.g., an earthquake, volcanic eruptions, or heavy rainfall) is referred to as an event-based landslide inventory. Owing to the advancements in earth observation technology, such as multitemporal high-resolution optical satellite remote sensing, more high-quality earthquakeinduced landslide inventories have been developed. In particular, since the 2008 Wenchuan earthquake, the establishment of coseisimic landslide inventories has made great progress. At present, there are roughly 46 detailed coseismic landslide databases mapped as polygons [10–13]. However, unlike earthquake events, the construction speed of landslide inventories triggered by heavy rainfall events is still relatively slow, and currently there are only a few heavy rainfall-induced landslide inventories [14–16]. The main reason is that clouds are often a major obstacle in the affected areas, which may limit the visibility of satellite images and thus affect the visual interpretation of rainfall-induced landslides [15]. At present, there are 16 public landslide inventories triggered by heavy rainfall events around the world, with the majority of these landslide databases being on a small scale. The southeast coastal region in China is economically developed and densely populated. Influenced by monsoon rainfall, this area is also considered a landslide-prone zone [17,18]. Once landslides occur, the social and economic losses in this area will be huge. A comprehensive rainfall-induced landslide database not only contributes to a deeper understanding of the event occurrence but also provides data support for the subsequent in-depth analysis of the formation and evolution of the geological disaster chain[15,19]. However, there are few rainfall-induced landslide inventories for a single event in the southeast coastal region, and thus more analyses are needed for rainfall-induced landslide inventories in this area.

Rainfall-induced landslide susceptibility can provide valuable information for landslide risk assessment. Currently, there are two quantitative methods for assessing the susceptibility for rainfall-induced landslides, which include the data-driven methods based on mathematical methods and physical-based methods that couple the hydrological models and infinite slope stability models. For the data-driven method, the relationship between the influencing factors and the landslide occurrence are analyzed by mathematical models [20–22]. Currently, many models have been widely used in landslide susceptibility mapping, particularly with the development of machine learning technology, such as logistic regression [23,24], random forest [25], artificial neural network [26], convolutional neural network (CNN) [27], support vector machine (SVM) [28], and decision tree [29]. However, the outcomes of landslide susceptibility mapping based on the data-driven method could be subject to considerable uncertainties due to errors and variability in model choice, data selection, system understanding of weighting factors, and human judgment [30,31]. Meanwhile, the data-driven model does not possess the timeliness of emergency assessment for a single triggering event, because it requires sufficient landslide data to establish the susceptibility assessment model. As a consequence, assessment results frequently lag behind practical application and cannot serve the emergency assessment in a short time [32,33]. Otherwise, due to the fact that the majority of these models are trained by regional landslide data and are thus limited by regional geological and geomorphic characteristics [14,34,35], the data-driven model's applicability in different areas is greatly diminished. However, the physically based landslide susceptibility assessment can better solve the above problems.

Unlike the data-driven method, the physically based method does not take into account actual landslide data, but rather simulates the physical process of rainfall-induced landslide occurrence by coupling the hydrological and infinite slope models [36]. The physically based method has been pervasively used because of its high predictive capability and the most acknowledged feasibility for a quantitative assessment of the effects of the individual parameters that contribute to landslide initiation [37] and it is a useful tool for determining the susceptibility zonation of rainfall-induced shallow landslides [38]. In addition, the wide application of GIS technology facilitates the wide application of physical models in large areas [39,40]. Due to its preferable practicability and wide regional applicability, physically based models are popular in the spatial prediction of regional rainfall-induced landslides [41–45]. In recent years, some physically based models for rainfall-induced landslide susceptibility mapping have been developed, such as the TRIGRS model [40], the Slip model [46–48], the GIS-TiVaSS model [45,49], the GIS-TISSA model [50], the CRESTSLIDE model [51,52], and the HIRESSS model [53–55]. Among them, the TRIGRS model, which accounts for transient pore water pressure, can predict the impact of heavy rainfall on groundwater changes in a short period. At present, it is the most widely used physically-based model of slope instability [41,56–58], and has been used in many countries around the world, including Italy, the United States, China, South Korea, and Southeast Asia [59–64]. However, the application of the TRIGRS model in China's southeast area is limited so it is necessary to investigate the applicability of the model in the spatio-temporal prediction of rainfall-induced landslides in the southeast mountainous area.

Longchuan County experienced continuous heavy rainfall from 10 June 2019 to 13 June 2019. Extensive landslides, collapses, and debris flows occurred in the villages of Longchuan county. A total of 352 villages of Longchuan County were devastated to varied degrees, of which Mibei village in Beiling town was the most severely hit with 1571 individuals affected, 120 buildings fully collapsed, and more than 100 houses damaged. The direct economic loss of this event reached CNY 110 million, exerting a significant impact on the normal productivity and lives of local residents. Thus, the objectives of this study are: (1) establishing a landslide inventory including landslides induced by the 2019 Longchuan heavy rainfall event and analyzing the spatial distribution of landslides with topographical, geological, and hydrological factors; (2) conducting the physically based susceptibility assessment based on a new open-source tool named MAT.TRIGRS(V1.0) for predicting the spatiotemporal distribution of rainfall-induced landslides and backanalyzing the response of the rainfall process on the change of landslide stability.

#### **2. Study Area**

Longchuan County is situated in the northeast of Guangdong Province, spanning from 23.8◦N to 24.7◦N of latitude and from 115.0◦ to 115.6◦E of longitude, and covers an area of approximately 3089 km<sup>2</sup> on the surface. The study area is Beiling Town, which is located in the north of Longchuan County and the upper reaches of Dongjiang River. The climate in the region is subtropical monsoon with abundant rainfall and sunshine. The annual rainfall is 1500 mm, and the average temperature is about 22 ◦C. The study area experiences the most rainfall in May, June, and July. The geomorphic unit of the study area is a hilly landform with an elevation range from 100 m to 1100 m (Figure 1). The mountains are steep, and the peaks are conical due to the relatively developed hydrographic nets and strong topographic cutting in this area. As a result, numerous "V" shaped valleys developed in this area, with slope angles ranging from 20 to 50 degrees. The main lithology of the study area is acid intrusive rock of Ordovician and Silurian, mainly monzogranite (O3-S1), which accounts for more than 70% of the rock in the whole study area (Figure 2a). In addition, tuff of Yousheng formation of Middle Cretaceous (K2ys) and quartz mica schist of Daganshan formation of Sinian(Z2djs) also developed in this area (Figure 2a). The main land use type is forest, which accounts for 80% of the whole study area, followed by cropland, accounting for more than 10% (Figure 2b).

**Figure 1.** Mapping shows the location and elevation of the study area; (**a**) Guangdong Province; (**b**) location of Longchuan county; (**c**) the elevation and water net distribution of the study area.

Due to the unique geographical and climatic conditions, Longchuan area experiences several large or small rainstorms every year, making it one of most vulnerable zones to geological disasters. From 10 June 2019 to 13 June 2019, Longchuan County suffered continuous heavy rainfall; this rainfall event triggered a large number of landslides. As far as local people can recollect, since the evening of the 10 June 2019, transportation has been disrupted, communication has been lost, and electricity has been cut off. Meanwhile, the settlement below the mountain was engulfed in mist, and the sound of collapses and landslides was constant.

**Figure 2.** (**a**) Geological map of the study area obtained from 1:200,000 geological maps published by China Geological Survey (http://dcc.cgs.gov.cn/, accessed on 1 July 2022); (**b**) the land use type map of the study area derived from the 10-m resolution global land cover results [65].

#### **3. Data and Method**

#### *3.1. Landslide Mapping*

The availability of high-resolution satellite images on the Google Earth (GE) platform allowed us to conduct a detailed visual interpretation of landslides [66,67]. The remote sensing images used for landslide interpretation in this study are based on the GE platform. It was important that the high-resolution satellite image covered the entire study area, and the dates of images before and after the rainfall event were mainly in January 2019 and January 2021. Meanwhile, given the relatively long interval between the images before and after the rainfall event, we obtained the Sentinel-2 images with 10 m resolution as a supplementary (the pre- and post-events images were 17 April 2019 and 24 September 2019, respectively) (Figure 3). The landslide inventory was checked by Sentinel-2 images to ensure that the interpreted landslides were caused by the 2019 rainfall event. The reason for selecting these two images was that they had the closest interval between rainfall events without cloud cover in the study area. Landslides were identified by visual interpretation and mapped as polygons. Since the study area has high vegetation coverage, landslides can be better delineated by satellite images before and after this event. Figure 3 depicts the Sentinel-2 satellite images before and after the rainfall. According to remote sensing images, most landslides triggered by this event were small and medium-scale shallow landslides, and a majority of them were located near Mibei village, showing obvious group-occurring characteristics (Figure 4).

**Figure 3.** Mapping shows the Sentinel-2 images before and after the rainfall event; (**a**) satellite image before rainfall event taken on 17 April 2019; (**b**) satellite image before rainfall event taken on 24 September 2019.

**Figure 4.** (**a**) Aerial photograph of Mibei village after the rainfall event, houses are damaged by rainfall-induced landslides; (**b**) group-occurring shallow landslides; (**c**) the landslide damaged the hillside residences, and the floors on the second floor crashed on the first floor; (**d**) road damage caused by landslides (Picture source: http://www.gdlctv.com/Pc/index/new\_detalis.html?id=3320, accessed on 25 June 2022).

#### *3.2. Rainfall Data*

We collected the precipitation data over the past two decades from 2000 to 2020 in Longchuan County (Figure 5). The results show that the average rainfall remained between 1200 and 2400 mm, with periodic fluctuations. The annual rainfall in 2006 and 2017 was unusually high, reaching almost 2300 mm or more. In comparison, the annual rainfall in 2019 was low with 1500 mm, which was roughly the same as the recent 20-year average (Figure 5a). Comparing the monthly rainfall in 2019 with the average value over the past two decades (Figure 5b), we also found that the rainfall from March to June in 2019 was higher than the monthly average rainfall in the last 20 years. The precipitation in June of 2019 was 300 mm, slightly higher than the monthly average rainfall of 250 mm in previous years.

**Figure 5.** Monthly rainfall data of Longchuan County in the past 20 years from 2000 to 2020; (**a**) monthly and annual average rainfall data over the last 20 years; (**b**) comparing the monthly rainfall in 2019 with the average value over the last two decades.

We obtained the data for the rainfall every 12 h based on the rainfall stations of China Meteorological Administration. Eight national rainfall stations within 50 km of the study area were utilized for interpolation, and the most popular Kriging interpolation algorithm was used to obtain the spatial distribution of rainfall (Figure 6). The results show that this rainfall event occurred primarily from 10 June 2019 to 13 June 2019 (Figure 7). The cumulative rainfall was basically the same, remaining at 210 to 220 mm, with rainfall in the west slightly higher than that in the east (Figure 6). Figure 8 shows the distribution of daily rainfall from 10th to 13th of June during this rainfall event. The rainfall on 10 June 2019 was the heaviest, peaking at around 120 mm, accounting for more than half of this rainfall event. The rainfall for the next three days was expected to be around 20–40 mm. Otherwise, the spatial change of daily rainfall in the study area from June 10 to 13 was relatively small, and the difference of daily rainfall of the study area was essentially maintained within 10 mm.

**Figure 6.** Mapping shows the spatial distribution of total rainfall from 10 June 2019 to 13 June 2019.

**Figure 7.** Data of the rainfall every 12 h from the national rainfall stations in the study area from 1 May to 30 June.

**Figure 8.** Mapping shows the distribution of daily rainfall from the 10th to 13th of June during this rainfall event; (**a**) 10 June 2019; (**b**) 11 June 2019; (**c**) 12 June 2019; (**d**) 13 June 2019.

#### *3.3. Data Related to Other Influencing Factors*

To assess the role of topographic, geologic, and hydrologic factors on the distribution of rainfall-induced landslides, we obtained several terrain metrics (i.e., elevation, hillslope gradient, and topographic relief) and lithologic and hydrological data. The elevation data were derived from ALOS PALSAR DEM with 12.5 m resolution, which were then resampled into a 5 m resolution based on the bilinear algorithm. The hillslope gradient and slope aspect were derived from the DEM data. In addition, we estimated the topographic relief from the elevation range within a 1.0 km radius. TWI was computed using GRASS GIS and the DEM data. Drainages were also derived from DEM by AcrGIS. The road data were downloaded from the OpenStreetMap Data (https://master.apis.dev.openstreetmap.org/export#map= 11/35.2510/103.4308, accessed on 5 June 2022). The lithology data are obtained from 1:200,000 geological maps published by China Geological Survey (http://dcc.cgs.gov.cn/, accessed on 1 July 2022). The land use type data were derived from the 10 m resolution global land cover results [65]. The spatial distribution of the above influencing factors was converted into a raster format with a grid cell size of 5 m. Finally, seven influencing factors were considered for the statistical analysis, including the elevation, hillslope gradient, relief, slope aspect, land use type, road density, and distance to river (Figure 9). Meanwhile, the relationship between different influencing factors and the occurrence of landslides were analyzed by the polygon feature.

**Figure 9.** Mapping showing the distribution of the influencing factors in the study area; (**a**) slope angle; (**b**) topographic relief; (**c**) topographic wetness index; (**d**) aspect; (**e**) road density; (**f**) distance to river.

#### *3.4. TRIGRS Model*

The TRIGRS model (Transient rainfall infiltration and grid-based regional slopestability model) is a widely used and effective evaluation model of rainfall-induced shallow landslide susceptibility [68,69]; the model was developed by the United States Geological Survey (USGS) [40,70] and written by Baum et al in FORTRAN [40], and it needs specific input parameters, mainly including rainfall parameters, soil mechanics parameters, and hydrological parameters of the study area. Following the determination of the parameters, the grid stability caused by the change of transient pore water pressure of each grid during the rainfall period is calculated based on the GIS platform for the purpose of evaluating the slope stability of all grids in a certain rainfall period.

Iverson [36] linearized the solution of the Richards equation and this serves as the foundation for the infiltration models for wet initial conditions. It consists of a steady component and a transient component of seepage. The steady seepage is determined by the initial depth of the water table and steady infiltration rate. Under steady infiltration, the slope is stable. Transient infiltration is the short-term change in pore water pressure caused by rainfall. The infinite slope model is then applied using the computed transient pore water pressure. The generalized solution in TRIGRS is:

$$\begin{split} \psi(Z,t) &= (Z-d)\beta + 2\sum\_{n=1}^{N} \frac{l\_{n}}{\mathcal{X}\_{s}} H(t-t\_{n}) [D\_{1}(t-t\_{n})]^{\frac{1}{2}} \\ \overset{\text{res}}{=} & \left\{ \text{ierfc}\left[\frac{(2m-1)d\_{\mathcal{L}\mathcal{X}} - (d\_{\mathcal{L}\mathcal{X}} - \mathbb{Z})}{2[D\_{1}(t-t\_{n})]^{\frac{1}{2}}}\right] + \text{ierfc}\left[\frac{(2m-1)d\_{\mathcal{L}\mathcal{X}} + (d\_{\mathcal{L}\mathcal{X}} - \mathbb{Z})}{2[D\_{1}(t-t\_{n})]^{\frac{1}{2}}}\right] \right\} - \\ 2\sum\_{n=1}^{N} \frac{l\_{n}}{\mathcal{X}\_{s}} H(t-t\_{n+1}) [D\_{1}(t-t\_{n+1})]^{\frac{1}{2}} \sum\_{n=1}^{\infty} \left\{ \text{ierfc}\left[\frac{(2m-1)d\_{\mathcal{L}\mathcal{X}} - (d\_{\mathcal{L}\mathcal{X}} - \mathbb{Z})}{2[D\_{1}(t-t\_{n+1})]^{\frac{1}{2}}}\right] + \text{ierfc}\left[\frac{(2m-1)d\_{\mathcal{L}\mathcal{X}} + (d\_{\mathcal{L}\mathcal{X}} - \mathbb{Z})}{2[D\_{1}(t-t\_{n+1})]^{\frac{1}{2}}}\right] \right\} \end{split} \tag{1}$$

where *ψ* is the groundwater pressure head; *t* is time; *N* is the total number of time intervals; Z is depth below the ground surface in the vertical coordinate direction; d is the depth of steady-state water table; *dLZ* is the depth of the impermeable basal boundary; *<sup>β</sup>* = cos2 *<sup>δ</sup>* − (*IZLT*/*Ks*), *<sup>δ</sup>* is the slope angle; *IZLT* is the steady surface flux; Ks is the saturated hydraulic conductivity; *InZ* is the the surface flux or rainfall intensity for the nth time interval; *<sup>D</sup>*<sup>1</sup> = *<sup>D</sup>*0/ cos<sup>2</sup> *<sup>δ</sup>*, *<sup>D</sup>*<sup>0</sup> is the saturated hydraulic diffusivity; and *<sup>H</sup>*(*<sup>t</sup>* − *tn*) is the Heaviside step function in which *tn* is the time at the nth time interval in the rainfall infiltration sequence.

$$
\operatorname{ierfc}(\eta) = \frac{1}{\sqrt{\pi}} \exp\left(-\eta^2\right) - \eta \operatorname{erfc}(\eta) \tag{2}
$$

where *er f c*(*η*) is the complementary error function.

The model calculates infiltration (*I*) at each cell as the sum of precipitation (*P*) and any runoff from upslope cells (*Ru*), with the caveat that infiltration cannot exceed the saturated hydraulic conductivity (*Ks*):

$$I = P + R\_{\text{u}}, \text{ if } P + R\_{\text{u}} \le \mathcal{K}\_s \tag{3}$$

$$I = \mathcal{K}\_{\mathfrak{s}\prime} \text{ if } P + \mathcal{R}\_{\mathfrak{u}} > \mathcal{K}\_{\mathfrak{s}} \tag{4}$$

When *P* + *Ru* exceeds *Ks* in a cell, the excess is considered runoff (*Rd*) and is diverted to nearby downslope cells.

$$R\_d = P + R\_\mu - K\_{\mathfrak{s}\_\prime} \text{ if } P + R\_\mu - K\_\mathfrak{s} \ge 0 \tag{5}$$

$$R\_d = 0,\text{ if }P + R\_\text{\textquotedblleft}R - K\_\text{\textquotedblright} < 0\tag{6}$$

The TRIGRS model calculates the slope stability using an infinite-slope stability analysis (Equation (7)), as described in Iverson [36]. The ratio of resistant basal coulomb friction

to gravitationally induced downslope basal driving stress characterizes the instability of an infinite slope in the failure analysis [71]. This ratio *Fs*, is computed at depth *Z* by

$$F\_s(Z,t) = \frac{\tan \phi'}{\tan \delta} + \frac{c' - \psi(Z,t)\gamma\_w \tan \phi'}{\gamma\_s Z \sin \delta \cos \delta} \tag{7}$$

where *c* is the cohesion of the soil, *ϕ* is the friction angle, *γ<sup>s</sup>* is the unit weight, and *γ<sup>w</sup>* is unit weight of groundwater.

The flow chart of this study is shown in Figure 10.

**Figure 10.** Flow chart of this study.
