Adopting the I3–R24 rainfall index and landslide susceptibility for the establishment of an early warning model for rainfall-induced shallow landslides

Abstract

Rainfall-induced landslides number among the most devastating natural hazards in the world and early warning models are urgently needed to reduce losses and fatalities. Most landslide early warning systems are based on rainfall thresholds defined on the regional scale, regardless of the different landslide susceptibilities of various slopes. Here we divided slope units in southern Taiwan into three categories (high, moderate and low) according to their susceptibility. For each category, we established separate rainfall thresholds so as to provide differentiated thresholds for different degrees of susceptibility. Logistic regression (LR) analysis was performed to evaluate landslide susceptibility by using event-based landslide inventories and predisposing factors. Analysis of rainfall patterns of 941 landslide cases gathered from field investigation led to the recognition that 3 h mean rainfall intensity ($I_{\mathrm{3}})$ is a key rainfall index for rainfall of short duration but high intensity; in contrast, 24 h accumulated rainfall ($R_{\mathrm{24}})$ was recognized as a key rainfall index for rainfall of long duration but low intensity. Thus, the $I_{\mathrm{3}}$–$R_{\mathrm{24}}$ rainfall index was used to establish rainfall thresholds in this study. Finally, an early warning model is proposed by setting alert levels including yellow (advisory), orange (watch) and red (warning) according to a hazard matrix. These differentiated thresholds and alert levels can provide essential information for local governments to use in deciding whether to evacuate residents.

Introduction

Rainfall-induced landslides number among the most perilous natural hazards, causing severe casualties and economic losses worldwide (Ayalew, 1999; Evans et al., 2007; Tsou et al., 2011; Petley, 2012; Wang et al., 2015; Iverson et al., 2015; Sassa et al., 2015; Fan et al., 2017). Therefore, much efforts has been made to evaluate landslide susceptibility and thereby set criteria for issuing alerts that can save lives and property.

Landslides can be triggered by either rainfall or earthquakes (Dadson et al., 2003; Lee et al., 2004; Lin et al., 2008; Chen et al., 2011). In Taiwan, monsoons and typhoons bring great amounts of rainfall, up to 3000 mm year${}^{-\mathrm{1}}$, and numerous landslides cause casualties every year. Therefore, recognizing the areas where rainfall-induced landslides might occur is an urgent issue. Here we adopted statistical method for the analysis of landslide susceptibility based on the assumption that the predisposing factors that cause landslides in a region are similar and can be used for predicting the locations of landslides in the future (Guzzetti et al., 1999). In previous research, several statistical models have been proposed, as well as widely utilized in landslide susceptibility analysis, and logistic regression was one of the most used methods (Guzzetti et al., 1999; Lee et al., 2004, 2008a, b; Lee, 2014).

Rainfall thresholds for landslides can be categorized as either statistical approaches or deterministic approaches. In the former method, thresholds are decided by collecting historical landslide cases and analyzing their rainfall parameters and the probability lines of rainfall conditions (Caine, 1980; Guzzetti et al., 2008). In the latter method, thresholds are decided by calculating the safety factors of each slope or grid with geomaterial and rainfall parameters (Terlien, 1998; Kim et al., 2010). Statistical rainfall thresholds for shallow landslides have been well discussed (Guzzetti et al., 2007). They can be classified mainly into five categories: intensity–duration (Brunetti et al., 2010; Zhou et al., 2014; Pradhan et al., 2017), accumulated rainfall–duration (Martelloni, 2011; Rosi et al., 2012; Vessia et al., 2014; Gariano et al., 2015; Rossi et al., 2017), accumulated rainfall (Corominas and Moya, 1999; Bell and Maud, 2000), intensity–accumulated rainfall (Hong et al., 2005) and accumulated rainfall–accumulated rainfall (Osanai et al., 2010; Turkington et al., 2014).

Most of the studies set only one threshold for their study areas despite differences in the physical settings (geology, geomorphology and meteorological conditions) of the regions. Recently, some studies have subdivided their study areas into several homogeneous sub-zones to discuss the influence of physical settings on thresholds (Hong and Adler, 2008; Segoni et al., 2014, 2015; Lee et al., 2015; Rosi et al., 2015, 2016; Peruccacci et al., 2017). However, for a smaller area such as slope units, differences in susceptibility may lead to different warning thresholds (Yang and Adler, 2008; Segoni et al., 2015; Lee et al., 2015). For example, the warning threshold of a high-susceptibility slope is likely to be lower than that of a low-susceptibility slope. To reduce this gap in knowledge, we focused on shallow landslides of the debris fall, debris topple, debris slide, earth fall, earth topple and earth slide types proposed by Varnes (1978) and divided slope units according to three different landslide susceptibility levels (high, moderate and low). After that, we established their rainfall thresholds separately. Furthermore, we set alert levels by adopting a hazard matrix and examined whether differentiated warning thresholds for different degrees of susceptibility existed. Moreover, given the importance of validating the performance of a landslide early warning model, especially the false alarms and missed alarms, to make it feasible for further practical application (Calvello et al., 2015; Devoli et al., 2015; Piciullo et al., 2017; Segoni et al., 2018), we also adopted skill scores to verify our results.

Study area

Taiwan is located in the western Pacific Ocean, on the convergent plate boundary zone of the Philippine Sea Plate and the Eurasian Plate. The orogenic uplift rate is 5–7 mm year${}^{-\mathrm{1}}$ (Willett et al., 2003); however, the exhumation rate is also as high as 3–6 mm year${}^{-\mathrm{1}}$ (Dadson et al., 2003) due to the fractured geological materials and the high mean annual precipitation of 2500–3000 mm brought by typhoons and monsoons every year (Hsu, 2013). The frequent natural disasters and high population density (23 million people in 36 000 km${}^{\mathrm{2}})$ of Taiwan make it one of the most exposed countries to multiple hazards (Dilley et al., 2005).

The study area, located in southern Taiwan (red box in Fig. 1), includes a mosaic of 47 1 : 25 000 scale maps (about 7258.5 km${}^{\mathrm{2}})$ and covers densely inhabited and landslide-threatening hillslopes. The elevation ranges from 3243 m in mountain areas to 0 m in plain areas, while the gradient ranges from 87 to 0${}^{\circ }$. The lithological units are mainly sedimentary rocks composed of sandstone, shale, mudstone and conglomerate in the Western Foothills, as well as metamorphic rocks composed of slate, argillite and metasandstone in the Central Range.

Geomorphological and geological settings of the study area. The elevation ranges from 3243 m in the eastern mountain area to sea level in the western plains area. Lithological units are mainly metamorphic rocks in the Central Range and sedimentary rocks in the Western Foothills. Rainfall data of 423 stations in Taiwan (96 of which are located in the study area) were collected for the interpolation and analysis of the triggering rainfall of landslides.

[Figure omitted. See PDF]

Available data

Landslide inventory

Landslide inventories are essential for the assessment of landslide susceptibility or spatiotemporal land changes (Van Westen et al., 2003; Guzzetti et al., 2012; Samia et al., 2017; Valenzuela et al., 2017). Four procedures for the generation of rainfall-induced landslide inventories were followed in this study. First, barren lands were interpreted manually from SPOT 5 images by drawing polygons in ESRI ArcGIS software. Second, aerial photographs and satellite images from Google Earth software were applied to identify whether the barren lands were landslides or agricultural land. In addition, the polygons situated on slopes having gradients greater than 55${}^{\circ }$ were marked as rockfalls according to the classification rules proposed by the Central Geological Survey, Taiwan (Central Geological Survey, 2008). Polygons marked as agricultural land or rockfalls were deleted from the inventories to ensure that only shallow landslides would be analyzed in the study. Third, we randomly selected landslides from inventories and verified the correctness of the locations and boundaries by fieldwork. Finally, event-based landslide inventories, including newly generated landslides and expanded landslides, were identified by comparing inventories before and after each event. In the end, six heavy rainfall events that triggered landslides were chosen and a total of seven landslide inventories were generated in this study (Table 1).

List of landslide inventories generated in this study.

Year	Event
2004	Before Typhoon Mindulle
2004	After Typhoon Mindulle
2005	After Typhoon Haitang
2006	After 0609 torrential rainfall
2007	After Typhoon Sepat
2008	After Typhoon Sinlaku
2009	After Typhoon Morakot

Landslide occurrence time and field investigation

Rainfall conditions such as intensity, duration and accumulated rainfall that induced landslides are key data in application of statistical methods to establish the rainfall thresholds for landslides (Guzzetti et al., 2007, 2008; Brunetti et al., 2010; Peruccacci et al., 2017). To analyze the rainfall conditions for each landslide case used in this study, a flowchart was proposed (Fig. 2). During field investigations, we verified the correctness of the landslide inventories and interviewed local residents to try to acquire the landslide occurrence times, since this information is rarely included in landslide inventories. The accuracy of landslide occurrence times is hard to evaluate due to the lack of video records; however, we focused on interviewing as many residents whose relatives were injured or houses were damaged or destroyed by the landslides as possible. Given the deep impressions left by such memories, we believe that the quality of landslide occurrence times might be improved. In contrast, detailed characteristics of the landslides, such as lithology, geological structure, joint, strength, area, depth and mechanism, were also recorded during the field work. Finally, 941 landslide cases, including their occurrence times (date and hour) and the characteristics of the landslides, were gathered for further analysis.

Flowchart of landslide occurrence time gathering during field investigation (left), locations of landslide cases with the occurrence times used in this study (middle) and the pictures of interviewing residents (right). To improve the quality of this key information, landslide occurrence times were obtained from local residents, especially those whose relatives were injured or houses were damaged or destroyed by the landslides.

[Figure omitted. See PDF]

Slope units

Slope units were used for the analysis of landslide susceptibility in this study (Carrara, 1988; Carrara et al., 1991; Carrara and Guzzetti, 1995; Guzzetti et al., 1999; Schlögel et al., 2017; Yang, 2017). To delineate the boundaries of slope units, $\mathrm{5}\times \mathrm{5}$ m digital elevation models (5 m DEMs) were acquired from the Ministry of the Interior, Taiwan. However, to reduce noise, we developed a Fortran program to obtain smoothed and resolution-reduced $\mathrm{10}\times \mathrm{10}$ m DEMs (10 m DEMS) by calculating the average value of each 2 by 2 grid in the 5 m DEMs. The resolution-reduced 10 m DEMs could generate some differences in the morphological analysis, but the expected scale of the landslide susceptibility in this study was set to 1 : 25 000, so differences smaller than 12.5 m could be ignored according to the relationship between mapping scale and 5 % acceptable error.

This study followed the method proposed by Xie et al. (2004) in delineating slope units according to gullies and ridges. First, gullies and watersheds were analyzed by successively using spatial analysis tools in ArcGIS: fill, flow direction, flow accumulation, stream link (with 2000 used as the threshold) and watershed. Second, reverse DEMs were generated by multiplying DEMs by $-\mathrm{1}$. In the reverse DEMs, ridges became gullies and could be analyzed by the same methods used in the first step. Third, the watersheds of the DEMs and reverse DEMs were transformed from rasters to polygons for further editing by using the “Raster to Polygon” tool in ArcGIS and then cut by each other to delineate the slope units. Finally, slope units were modified manually according to aspect and gradient. It is suggested that the aspect in a slope unit should be within three adjacent directions; e.g., northwest, north and northeast. In contrast, the difference in gradient should not be over 30${}^{\circ }$ in a slope unit, and slope units situated on flat areas, including alluvial deposits and terraces, were deleted. In addition, the area of each slope unit was set to around 5 ha; therefore, slope units smaller than 5 ha were combined with adjacent slope units and those larger than 5 ha were split into several smaller ones. Moreover, slope units delineated by parallel drainage on a dip slope were combined into one slope unit. After the editing, each slope unit is used for disaster prevention.

Landslide susceptibility factors

Several predisposing factors that might lead to landslides were selected initially in this study to construct a landslide susceptibility model for slope units. These factors included rock mass strength-size classification (RMSSC I–VII), dip slope, average slope, variance of slope, ratio of steep slope, total slope height, average elevation, average curvature, variance of curvature, fault density, fold density, average wetness, rainfall intensity, total rainfall, 3 h mean rainfall intensity ($I_{\mathrm{3}})$ and 24 h accumulated rainfall ($R_{\mathrm{24}})$. The relationships of these factors to landslides could be analyzed through graphic discrimination, including success rate curve, probability of failure curve and difference between landslide and non-landslide groups (Lee, 2014). After that, factor correlation analyses were applied to delete highly related factors to keep the factors used in the landslide susceptibility model as independent as possible (Table 2).

In terms of geological factors, the lithology of a location is essential for the analysis of landslide susceptibility. However, our study area had more than 50 detailed types of lithology, which was unfavorable for the analysis. Therefore, we adopted the 1 : 25 000 rock mass strength–size classification maps from Central Geological Survey, Taiwan, to replace the use of lithology (Franklin, 1975; Central Geological Survey, 2008). In addition, the dip slope inventory interpreted manually from 1 : 5000 aerial photographs by the Central Geological Survey, Taiwan, was also adopted (Central Geological Survey, 2008). The fold density was also calculated by dividing the total length of all the folds by the total area in each slope unit.

Predisposing factors and their coefficient in logistic regression analysis.

Code	Factor item	Coefficient
$L_{\mathrm{01}}$	RMSSC I	–
$L_{\mathrm{02}}$	RMSSC II	–
$L_{\mathrm{03}}$	RMSSC III	$-\mathrm{0.874}$
$L_{\mathrm{04}}$	RMSSC IV	$-\mathrm{0.099}$
$L_{\mathrm{05}}$	RMSSC V	0.314
$L_{\mathrm{06}}$	RMSSC VI	$-\mathrm{0.384}$
$L_{\mathrm{07}}$	RMSSC VII	–
$F_{\mathrm{01}}$	dip slope	0.207
$F_{\mathrm{02}}$	average slope	0.265
$F_{\mathrm{03}}$	variance of slope	0.098
$F_{\mathrm{04}}$	ratio of steep slope	0.344
$F_{\mathrm{05}}$	average curvature	0.016
$F_{\mathrm{06}}$	variance of curvature	0.161
$F_{\mathrm{07}}$	fold density	0.013
$F_{\mathrm{08}}$	average wetness	0.061
$F_{\mathrm{09}}$	3 h mean rainfall intensity ($I_{\mathrm{3}})$	$-\mathrm{0.817}$
$F_{\mathrm{10}}$	24 h accumulated rainfall ($R_{\mathrm{24}})$	0.665
$C$	Constant	0.057

For morphological factors, the average slope and the variance of slope were obtained by averaging and calculating the standard deviations of all the grid cells in the slope unit separately. In addition, shallow landslides are prone to occur on steep slopes; therefore, we also used the “ratio of steep slopes” to present how many steep slopes existed in a slope unit. It was found after trial and error that a threshold of gradient higher than 30${}^{\circ }$ had a higher relationship with landslide susceptibility. Thus, we calculated the area where the gradient was greater than 30${}^{\circ }$ ($A_{>\mathrm{30}})$ as well as the total area ($A_{\mathrm{total}})$ of each slope unit. Therefore, the ratio of steep slope could be calculated by dividing $A_{>\mathrm{30}}$ by $A_{\mathrm{total}}$. The average curvature and variance of curvature could also be calculated by the same method as slope in the ArcGIS software. The average wetness was calculated by averaging the wetness index of grid cells in a slope unit. This factor represents the effect of morphology on soil wetness. When the drainage area is larger and the slope is gentler, the water content in the soil will also be higher and therefore make a slope more prone to failure. The wetness index can be calculated according to the method proposed by Wilson and Gallant (2000) as follows: $$\mathit{\omega }=\ln \left({\displaystyle \frac{A_{\mathrm{s}}}{\tan \mathit{\theta }}}\right),$$ where $\mathit{\omega }$ is wetness index, $A_{\mathrm{s}}$ is the drainage area of a specific grid cell expressed as square meters and $\mathit{\theta }$ is the slope of the grid expressed in radians.

For triggering factors, we collected hourly rainfall data of 423 rain stations provided by Central Weather Bureau, Taiwan (96 of which were located in our study area, shown in Fig. 1) and analyzed both the 3 h mean rainfall intensity ($I_{\mathrm{3}})$ and the 24 h accumulated rainfall ($R_{\mathrm{24}})$ of each station for each rainfall events in Table 2. After that, we used the linear mode of ordinary kriging and applied the default setting in Surfer software to obtain the rainfall distribution of the whole study area. Whenever Taiwan has a typhoon event, the Central Weather Bureau issues disaster prevention alerts. Therefore, we counted the time that the first alert was issued as the beginning of the rainfall event and the time that the alert was canceled as the end of the rainfall event to calculate rainfall amounts. Our reasons for choosing 3 h mean rainfall intensity ($I_{\mathrm{3}})$ and 24 h accumulated rainfall ($R_{\mathrm{24}})$ as factors will be explained in detail in Sect. 4.2.

Methodology

Landslide susceptibility analysis

The main purpose of landslide susceptibility analysis is to determine the effectiveness of each predisposing factor and the relative possibility of landslide occurrence in a specific area. Several methods can be used to analyze landslide susceptibility. The deterministic method uses a physical model and geotechnical material properties to determine the safety factor of slopes; however, precise parameters of materials are difficult to obtain, especially on a regional scale (Montgomery and Dietrich, 1994; Van Westen and Terlien, 1996). The qualitative and semi-quantitative methods rely on the experience and knowledge of the experts who carried out the analysis; however, these results might vary from one expert to another. The machine learning method uses multiple samples to build a model by trial and error; however, it is time consuming (Gorsevski and Jankowski, 2010; Yeon et al., 2010; Yilmaz, 2010; Marjanovic et al., 2011; Lee et al., 2012; Song et al., 2012). The statistical method also requires numerous samples for the training; however, it is more efficient, especially when dealing with regional-scale analyses, and can avoid the uncertainty of material parameters as well as differences in expert experience. Recently, nonlinear analysis, a statistical method, has been used for the analysis of complex landslide phenomena. Methods such as logistic regression (Yilmaz, 2010; Lee et al., 2012, 2014, 2015; Schlögel et al., 2017) and discriminant analysis (Lee et al., 2004, 2008a, b) are often used to analyze landslide susceptibility. In this study, we applied logistic regression (LR) to evaluate the susceptibility of each slope unit (Guzzetti et al., 1999; Ayalew and Yamagishi, 2005). The LR function can be expressed as follows: $$P={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{-z}}},$$ $$z=\sum _{i=\mathrm{1}}^m{L}_i{w}_i+\sum _{j=\mathrm{1}}^n{F}_j{w}_{m+j}+C,$$ where $P$ is landslide susceptibility, $L_i$ is RMSSC factor ($L_{\mathrm{01}}$ to $L_{\mathrm{07}}$ in Table 2), $F_j$ is other factors ($F_{\mathrm{01}}$ to $F_{\mathrm{10}}$ in Table 2), $w_i$ and $w_{m+j}$ are regression coefficients and $C$ is a constant. Six event-based landslide inventories in this study were used to label whether or not landslides occurred in the slope units. After that, all the slope units were divided randomly into two groups, one for training the model and the other for validation. The index indicating landslide or non-landslide was set as the dependent variable, and all the landslide susceptibility factors were set as covariates in SPSS for training of the model. After iterative training, the regression coefficients of each landslide susceptibility factor, as well as the success rate curve (SRC), the prediction rate curve (PRC) and the area under the curve (AUC), were reported in SPSS. The AUC can be used to examine if the model predicts landslides well, and the regression coefficients can be used for the prediction of landslide susceptibility. During the process of training, several details required attention. Because the non-landslide samples outnumbered the landslide samples, we randomly selected equal numbers of non-landslide and landslide samples for the training so as to avoid the effect of difference in quantity. In addition, different samples could have led to different results when selecting non-landslide samples randomly. To reduce this effect, we prepared several sets of randomly selected samples for the analysis of landslide susceptibility to test if the model was stable enough; i.e., the AUC would not vary severely when validating models with different sets of samples. Finally, the individual landslide susceptibilities of slope units were calculated with this model and classified into high, moderate and low susceptibility levels.

The 3 h mean rainfall intensity ($I_{\mathrm{3}})$ and 24 h accumulated rainfall ($R_{\mathrm{24}})$ were used as short-term and long-term rainfall indexes for the establishment of rainfall thresholds.

[Figure omitted. See PDF]

$I_{\mathrm{3}}$–$R_{\mathrm{24}}$ rainfall index and thresholds

Rainfall-induced landslides are triggered by either high intensity rainfall or high accumulated rainfall (Larsen and Simon, 1993; Corominas and Moya; 1999; Yu et al., 2006). To identify rainfall indexes responsible for landslides, the triggering rainfall, including the rainfall intensity ($I_{\mathrm{1}}$, $I_{\mathrm{2}}$, $I_{\mathrm{3}}$, $I_{\mathrm{4}}$, $I_{\mathrm{5}}$ and $I_{\mathrm{6}})$ and accumulated rainfall ($R_{\mathrm{6}}$, $R_{\mathrm{12}}$, $R_{\mathrm{24}}$, $R_{\mathrm{48}}$ and $R_{\mathrm{72}})$ of different time windows of each landslide case were analyzed according to the landslide occurrence time. The results revealed that 218 landslides occurred within the 3 h following the highest rainfall intensity, and 242 occurred within the 3 h following the 2nd or 3rd highest rainfall intensity (i.e., induced by high rainfall intensity), accounting for nearly 49 % of the landslide cases gathered in this study (Table 3). From these results, it became clear that in Taiwan, $I_{\mathrm{3}}$ is the most important index for landslides induced by rainfall of short duration but high intensity. In contrast, 481 landslides occurred close to the end of the rainfall events (i.e., induced by high accumulated rainfall), accounting for about 51 % of the total cases. Furthermore, analysis of the different accumulated rainfall indexes showed that 24 h accumulated rainfall had the lowest coefficient of variation (Table 4), indicating that this index was less dispersive than others and might be more suitable for serving as an accumulated rainfall index for establishing rainfall thresholds. The coefficient of variation can be calculated as follows: $$C_{\mathrm{v}}={\displaystyle \frac{\mathit{\sigma }}{\mathit{\mu }}},$$ where $C_{\mathrm{v}}$ is the coefficient of variation, and $\mathit{\sigma }$ and $\mathit{\mu }$ are the standard deviation and average of accumulated rainfall of all the cases used in this study respectively.

Type and the proportion of landslide occurrence times.

Type of landslide occurrence time	Amount (percentage)
Type A: within the 3 h following the highest rainfall intensity	218 (23 %)
(landslide induced by high rainfall intensity)
Type B: within the 3 h following the 2nd or 3rd highest rainfall intensity	242 (26 %)
(landslide induced by high rainfall intensity)
Type C: near the end of the rainfall event	481 (51 %)
(landslide induced by high accumulated rainfall)
Total	941 (100 %)

Based on these data and previous studies (Cheung et al., 2006; Liao et al., 2010), 3 h mean rainfall intensity ($I_{\mathrm{3}})$ and 24 h accumulated rainfall ($R_{\mathrm{24}})$ were respectively chosen as the short-term and long-term rainfall indexes for the establishment of the rainfall threshold (Fig. 3). We chose 3 h mean rainfall intensity here instead of 3 h accumulated rainfall to focus on rainfall of short duration but high intensity. Similarly, we chose 24 h accumulated rainfall to focus on rainfall of long duration but low intensity. Finally, rainfall thresholds were decided by plotting the $I_{\mathrm{3}}$ and $R_{\mathrm{24}}$ rainfall index of historical landslides in the $I_{\mathrm{3}}$–$R_{\mathrm{24}}$ diagram (Fig. 4). Here we used the ellipse as the threshold line, and the parameters $a$ (semi-major axis) and $b$ (semi-minor axis) of the ellipse were set according to the slope of best fit line obtained from the least square method. Thresholds such as 90, 60, 30, 15 % were determined according to the percentage of historical cases that could be enveloped under the threshold line; e.g., the 90 % threshold ($T_{\mathrm{90}\mathit{\%}})$ included 90 % of the historical cases. A higher threshold indicates a more dangerous condition for the occurrence of landslides. The original warning values of $I_{\mathrm{3}}$ and $R_{\mathrm{24}}$ of the 90, 60, 30, 15 % thresholds were equal to the semi-minor axis and semi-major axis of each threshold respectively. After that, $I_{\mathrm{3}}$ was rounded to the nearest 5 mm h${}^{-\mathrm{1}}$ and $R_{\mathrm{24}}$ was rounded to the nearest 50 mm for operational purposes, such as the evacuation of residents.

Establishment of $I_{\mathrm{3}}$–$R_{\mathrm{24}}$ rainfall thresholds for shallow landslides. The best fit line was derived by least square method, and the ratio of a and b was used as the ratio of the semi-major axis and semi-minor axis in the ellipse threshold line.

[Figure omitted. See PDF]

Coefficient of variation of different accumulated rainfall indexes.

Accumulated rainfall indexes	Coefficient of variation
6 h accumulated rainfall ($R_{\mathrm{6}})$	0.68
12 h accumulated rainfall ($R_{\mathrm{12}})$	0.47
24 h accumulated rainfall ($R_{\mathrm{24}})$	0.38
48 h accumulated rainfall ($R_{\mathrm{48}})$	0.41
72 h accumulated rainfall ($R_{\mathrm{72}})$	0.45

Landslide early warning model and validation

Landslide early warning model and alert considering both landslide susceptibility and rainfall thresholds.

[Figure omitted. See PDF]

The landslide early warning model in this study considered both landslide susceptibility and rainfall thresholds and alerts were determined by using a hazard matrix. As mentioned above, the LR method was applied to analyze the susceptibility of each slope unit. After that, all the slope units were categorized into high, moderate and low susceptibility levels. We consequently established rainfall thresholds for each susceptibility level separately and then set alerts of red, orange, yellow and green according to the level of danger.

High-susceptibility slopes might be more susceptible to rainfall. Hence, the alerts were set as red (extreme danger level) for rainfall conditions exceeding the 60 % threshold line, orange (high danger level) for those between the 60 and 30 % threshold lines, yellow (medium danger level) for those between the 30 and 15 % threshold lines and green (low danger level) for rainfall conditions lower than the 15 % threshold line (Table 5). For moderate-susceptibility slopes, the alerts were set as red for rainfall conditions exceeding the 90 % threshold line, orange for those between the 90 and 60 % threshold lines, yellow for rainfall conditions between the 60 and 30 % threshold lines and green for rainfall conditions lower than the 30 % threshold line. Low-susceptibility slopes should be less susceptible to rainfall. Hence, the alerts were set as orange for rainfall conditions exceeding the 90 % threshold line, yellow for those were between the 90 and 60 % threshold lines and green for those lower than the 60 % threshold line.

Several methods can be used for the validation of a landslide early warning model (Segoni et al., 2014, 2018; Gariano, 2015; Rosi et al, 2015; Piciullo et al., 2017; Krøgli et al., 2018). According to the analysis of Segoni et al. (2018), compiling a contingency matrix and calculating skill scores is the most commonly used method in recent years. Therefore, we applied this method and quantitatively validated our model with probability of detection (POD, also known as hit rate), probability of false detection (POFD, also known as false alarm rate) and probability of false alarm (POFA, also known as false alarm ratio). The contingency matrix is presented as Table 6. Comparing the observed events and the forecasted events produces four outcomes: true positive (TP), true negative (TN), false positive (FP) and false negative (FN). Because warnings against natural hazards are always issued by taking a village as a unit in Taiwan, here we validated our model on the village scale. TP indicated the number of villages for which warnings of slope units were issued and landslides did occur, while TN indicated the number of villages for which no warning was issued and no landslide occurred. In contrast, FP indicated the number of villages for which warnings were issued but no landslides occurred, also known as false alarms, while FN indicated the number of villages for which no warnings were issued but landslides did occur; i.e., missed alarms. In addition, in our study, we defined red (extreme danger level) and orange (high danger level) alerts as warnings issued for evacuation; i.e., the alarm zone in our model. In contrast, yellow (medium danger level) and green (low danger level) alerts were considered to indicate no need for evacuation; i.e., the no alarm zone in our model. POD, POFD and POFA can be calculated by the following equations: $$\mathrm{POD}={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}},$$ $$\mathrm{POFD}={\displaystyle \frac{\mathrm{FP}}{\mathrm{TP}+\mathrm{FP}}},$$ $$\mathrm{POFA}={\displaystyle \frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{TN}}}.$$ The ranges of POD, POFD and POFA are all between 0 to 1, and their optimal values are 1, 0 and 0, respectively.

Contingency matrix for the validation of landslide early warning model. True positive (TP), true negative (TN), false positive (FP) and false negative (FN).

	Observed events
		Yes	No
Forecasted events	Yes	True positive	False positive
		(TP)	(FP)
	No	False negative	True negative
		(FN)	(TN)

Results and discussions

Landslide susceptibility analysis

After several calibrations, the resultant model was obtained. The coefficients for each factor of LR are given in Table 2, and the landslide susceptibility of each slope unit was also calculated. To evaluate the quality of a predicted model, the success rate curve (SRC) and prediction rate curve (PRC) (Chung and Fabbri, 1999) were mapped, and then the area under the curve (AUC) was used to describe the model's ability to distinguish landslide and non-landslide (Yesilnacar and Topal, 2005). A higher AUC value indicated a better model for the prediction of landslides. If the AUC value was 0.5, it meant that the model did not predict the occurrence of the landslide better than a random approach. If the AUC value was 1.0, the capability of the model for predicting a landslide was perfect.

In our study, the AUCs were 0.745 and 0.691 in training and validation respectively, indicating that our LR model could identify 60 % of the landslides in the top 25 and 30 % of the highest susceptibility areas during training and validation (Fig. 5). These results showed that the LR model was acceptable in both the training and the validation. For a statistical landslide susceptibility analysis, it is essential to use as many samples as possible. However, we used slope units instead of grid units in this study for application to disaster prevention. This led to the reduction of samples, since one slope unit might equal hundreds of grids. Therefore, our AUC might not be considered high in comparison to a grid-based landslide susceptibility model.

To avoid over-training, it was necessary to validate the capability of the model. One common method is to divide the study area into sub-regions such as left and right, one for training and the other for validation (Chung and Fabbri, 2008). But this method might cause the loss of a training pattern in a small or particular geological region if the study area is extensive. To overcome this problem, we used multi-event data from the same area for training and testing. The data used in this study were therefore randomly divided into two portions, and several sets of data were tested. This approach would also solve the problems mentioned above.

Area under the curve (AUC) of training and validation of landslide susceptibility analysis.

[Figure omitted. See PDF]

$I_{\mathrm{3}}$–$R_{\mathrm{24}}$ rainfall threshold

Rainfall thresholds for southern Taiwan. The values were calculated as 90, 60, 30 and 15 % of the original threshold. After that, $I_{\mathrm{3}}$ was rounded to the nearest 5 mm h${}^{-\mathrm{1}}$ and $R_{\mathrm{24}}$ was rounded to the nearest 50 mm for operational purpose (e.g., evacuation). The original values are shown in parentheses.

[Figure omitted. See PDF]

We gathered a total of 941 landslide cases in this study and used 240 cases located in southern Taiwan, consisting of 110 high-susceptibility cases, 84 moderate-susceptibility cases and 46 low-susceptibility cases, to establish a susceptibility-based regional landslide early warning model. The ellipse-shaped $I_{\mathrm{3}}$–$R_{\mathrm{24}}$ rainfall thresholds for three different landslide susceptibility slopes are presented in Fig. 6. For practical use, the original threshold values of $I_{\mathrm{3}}$ and $R_{\mathrm{24}}$ (as shown in the parentheses in Table 7) were separately rounded to the nearest 5 mm h${}^{-\mathrm{1}}$ and the nearest 50 mm. It was found that the threshold values gradually increased as the susceptibility of slope units decreased for the same alert level, indicating that greater rainfall amounts would be needed when issuing alters on less susceptible slope units. These results showed that establishing rainfall thresholds according to different landslide susceptibilities and then setting alert levels by adopting a hazard matrix not only provided differentiated thresholds but also avoided the over- or underestimation of the thresholds for slopes.

Alerts and the corresponding danger levels, as well as suggested actions.

Alert	Danger level	Suggested action
Green	Low	–
Yellow	Medium	Notice announcements
Orange	High	Evacuation
Red	Extreme	Forced evacuation

After the establishment of the landslide early warning model, we converted the model to an early warning system (EWS) connected to the QPESUMS, which provides nearly real-time radar rainfall data, for disaster prevention. Based on the alerts present in the system, the corresponding danger levels and suggested actions for residents around the warning slope are shown in Table 8. During a yellow alert, residents should listen for further announcements and prepare for evacuation if the alert is raised to orange. When an orange alert is issued, residents should evacuate as quickly as possible because landslides are likely to occur, according to the validations shown in the next section. Finally, when a red alert is issued, evacuation may need to be enforced to protect residents from injury.

$I_{\mathrm{3}}$–$R_{\mathrm{24}}$ rainfall thresholds and alert of (a) high-susceptibility slope units (b) moderate susceptibility slope units and (c) low-susceptibility slope units for southern Taiwan.

[Figure omitted. See PDF]

Disastrous landslide cases in 2016 and the results of validation. Warnings could have been issued for all landslide cases in advance or at the time of occurrence (all times referred to are in UTC+8).

Landslide	Lithology	Landslide	Alert &	Occurrence of	Early ($+$),
susceptibility		area	issuing time	landslide	late ($-$)
Shihwen landslide (Shihwen village, Chunri Township, Pingtung County)
Low	Weathered	61 500 m${}^{\mathrm{2}}$	Orange, 11:00,	05:00,	$+\mathrm{18}$ h
	sandstone		14 September 2016	15 September 2016
Zhongmin landslide (Zhongmin Rd., Yanchao District, Kaohsiung City)
High	Mudstone interbedded	3500 m${}^{\mathrm{2}}$	Orange and red, 04:00,	10:00,	$+\mathrm{6}$ h
	with thin sandstone		28 September 2016	28 September 2016
Houcuo landslide (Houcuo Ln., Qishan District, Kaohsiung City)
Low	Conglomerate	4000 m${}^{\mathrm{2}}$	Orange, 03:00–04:00,	03:30,	0 h
			28 September 2016	28 September 2016

Validation of landslide early warning model

We validated our model with two kinds of data: (1) three disastrous shallow landslides in 2016 and the occurrence times provided by witnesses; (2) a landslide inventory of two historical typhoon events and the occurrence times reported by newspapers.

The first set of validation data showed that orange or red alerts could have been issued in advance for all of the disastrous landslides before the landslides occurred, according to the rainfall path in the $I_{\mathrm{3}}$–$R_{\mathrm{24}}$ diagram (Fig. 7; Table 9).

The Shihwen landslide occurred on a low-susceptibility slope. From the rainfall histogram and $I_{\mathrm{3}}$–$R_{\mathrm{24}}$ diagram (Fig. 7a), we knew that the occurrence time was quite close to the end of the rainfall event, and that the $I_{\mathrm{3}}$ was only 2.3 mm h${}^{-\mathrm{1}}$ while the $R_{\mathrm{24}}$ was 507.5 mm, indicating that accumulated rainfall might have been the principal cause of this case. The rainfall path showed that on 14 September, the alert was raised to yellow at 10:00 (GMT+8) and then to orange at 11:00 during the downpour. Then the precipitation fell for several hours, and the alert was lowered to yellow. However, when it rained again, the alert was raised back to orange at 18:00 and at 23:00, and the landslide occurred at 05:00 on 15 September.

The Zhongmin landslide occurred on a high-susceptibility slope. From the rainfall histogram and $I_{\mathrm{3}}$–$R_{\mathrm{24}}$ diagram (Fig. 7b), it was found that the occurrence time was also quite close to the end of the rainfall event and that the $I_{\mathrm{3}}$ was 8.3 mm h${}^{-\mathrm{1}}$ while the $R_{\mathrm{24}}$ was 479 mm. The high rainfall intensity (74 mm h${}^{-\mathrm{1}}$ at 04:00 on 28 September) and accumulated rainfall might both have contributed to this case. The rainfall path showed that on 28 September, the yellow alert was raised to orange and red at 04:00 during the high intensity rainfall mentioned above. After that, although the rainfall soon fell off, the landslide occurred 6 h later, at 10:00 on 28 September, during an orange alert.

Disastrous landslide cases in 2016, their rainfall histograms, and their rainfall paths in the $I_{\mathrm{3}}$–$R_{\mathrm{24}}$ diagram. The results showed that orange or red alerts could have been issued in advance for these disastrous landslides.

[Figure omitted. See PDF]

The Houcuo landslide also occurred on a low-susceptibility slope. From the rainfall histogram and $I_{\mathrm{3}}$–$R_{\mathrm{24}}$ diagram (Fig. 7c), we found that the occurrence time was close to the time that the highest rainfall intensity showed in the rainfall event, and the $I_{\mathrm{3}}$ was 24.3 mm h${}^{-\mathrm{1}}$ while the $R_{\mathrm{24}}$ was 291.3 mm, indicating that high rainfall intensity might have been responsible for this case. Due to this intensity, the rainfall path showed that the alert was raised from green to yellow and then to orange within just 1 h, from 03:00 to 04:00, on 28 September, and the landslide occurred at around 03:30.

For the second set of validation data, we applied the kriging method to interpolate spatial rainfall data and analyzed the alerts for each slope unit hour-by-hour. The results showed that the hit rates in the two historical typhoon events were all sufficiently high, according to the accumulative warning numbers relative to the numbers of landslide slopes (Fig. 8; Table 10).

During Typhoon Mindulle in 2004, landslides occurred in 10 911 slope units, including 5129 high-susceptibility slopes, 2750 moderate-susceptibility slopes and 3032 low-susceptibility slopes. According to newspaper reports, several landslides occurred at 10:00 and between 15:00 and 16:00 on 2 July 2004; however, most of the landslides occurred between 06:00 and 13:00 the next day, 3 July 2004 (blue dashed box in Fig. 8a). From the alert history (Fig. 8a), it was found that the peak number of orange and red alerts matched the reported occurrence times quite well. In addition, orange alerts, indicating the need for evacuation, had been issued for 8283 slope units during the whole event, accounting for 75.9 % of the slope units where landslides occurred in this event.

Validation of warnings issued for slope units and the hit rate during Typhoons Mindulle and Haitang.

Typhoon event	Landslide occurrence time (GMT+8)	Number of landslide slope units	Number of slope units	Hit rate
(year)	reported by newspapers	(number of high-, moderate- and	for which orange alerts
		low-susceptibility slope units)	had been issued
Mindulle (2004)	Mainly 06:00–13:00,	10 911	8283	75.9 %
	3 July 2004	(5129; 2750; 3032)
Haitang (2005)	From 05:00, 19 July	10 804	10 245	94.8 %
	to 06:00, 20 July 2005	(2592; 2355; 5857)

Validation of alerts issued for villages during Typhoons Mindulle and Haitang by using contingency matrix and skill scores.

	Observed events
	Typhoon Mindulle (in 2004)			Typhoon Haitang (in 2005)
		Yes	No		Yes	No
Forecasted events	Yes	220	30	Yes	194	24
	No	9	77	No	28	12
POD		0.961			0.874
POFD		0.280			0.667
POFA		0.120			0.110

Warning history of (a) Typhoon Mindulle and (b) Typhoon Haitang showing that the time at which our model issued alerts matched the landslide occurrence times reported by newspapers.

[Figure omitted. See PDF]

Typhoon Haitang in 2005 was another event of concern. Landslides occurred in 10 804 slope units, including 2592 high-susceptibility slopes, 2355 moderate-susceptibility slopes and 5857 low-susceptibility slopes. According to newspaper reports, landslides occurred between 05:00 on 19 July and 06:00 on 20 July 2005 (blue dashed box in Fig. 8b). From the alert history (Fig. 8b), it was found that landslides occurred immediately after the number of orange and red alerts increased sharply, and the peak number of orange and red alerts also matched the reported occurrence times quite well. Orange alerts had been issued for 10 245 slope units during the whole event, accounting for 94.8 % of the slope units where landslides occurred in this event. These results revealed that our model could provide valuable information for evacuation and disaster prevention.

In addition, the second set of validation data was also used to validate the warnings issued for villages during two typhoon events by adopting the contingency matrix and skill scores. According to the event-based landslide inventories, if any landslides were located in a village, the village was classified as “Yes” for observed events. If orange or red warning alerts were issued for slope units in a village, the village was classified as “Yes” for forecasted events. Based on these rules, the numbers of true positive (TP), true negative (TN), false positive (FP) and false negative (FN) were counted and the skill scores were calculated (Table 11). The probabilities of detection (PODs) of the two typhoon events were 0.961 and 0.874 respectively, indicating that most of the villages where landslides occurred could have been warned in advance. The probabilities of false detection (POFDs) of the two typhoon events were 0.280 and 0.667 respectively, suggesting that the model performed well for Typhoon Mindulle but might not be as ideal for Typhoon Haitang. Lastly, the probabilities of false alarm (POFAs) of the two typhoon events were 0.120 and 0.110 respectively, which meant that our model would not issue an excessive number of false alarms and was feasible for disaster prevention.

Conclusions

This study attempted to establish regional rainfall thresholds for shallow landslides according to their landslide susceptibility levels and set alerts with a hazard matrix to provide more detailed results for disaster mitigation.

Logistic regression (LR), a statistical method, was applied in this study to analyze the landslide susceptibilities of slope units. The areas under the curve (AUC) were 0.745 and 0.691 in the training and validation respectively. Due to our use of slope units instead of grid units in this study for application to disaster prevention, the number of our training samples was less, since one slope unit might equal hundreds of grids. Therefore, our AUC might not be considered high as compared to a grid-based landslide susceptibility model, but it was still acceptable for practical use.

This study also examined the relationships between rainfall indexes and the occurrence of landslides. From 941 landslide cases we gathered, it was found that 3 h mean rainfall intensity ($I_{\mathrm{3}})$ and 24 h accumulated rainfall ($R_{\mathrm{24}})$ were the most dominant short-term and long-term parameters responsible for rainfall-induced landslides in Taiwan. There were 460 cases (about 49 %) occurred within the 3 h following the highest, 2nd and 3rd rainfall intensities, while 24 h accumulated rainfall had the lowest coefficient of variation of the long-term rainfall indexes. The $I_{\mathrm{3}}$–$R_{\mathrm{24}}$ rainfall index was therefore used to establish rainfall thresholds.

We categorized the slope units into three landslide susceptibility levels (high, moderate and low) and then separately established a susceptibility-based regional rainfall threshold. We also set three alert levels, including red (extreme danger level), orange (high danger level) and yellow (medium danger level), by adopting a hazard matrix for application to evacuation decisions. It was found that the threshold values gradually increased as the susceptibility of slope units decreased for the same alert level, indicating that greater rainfall amounts would be needed when issuing alters on less susceptible slope units.

Validations using three disastrous shallow landslides in 2016 and two landslide inventories of historical typhoon events showed that, for the landslide cases in 2016, orange or red alerts could have been issued before the landslides occurred and the hit rates of the alerts issued for slope units in the two historical typhoon events were 75.9 and 94.8 % respectively, which are sufficiently high for a landslide early warning model. In addition, the skill scores applied to the validation of alerts issued for villages during two typhoon events showed that the probability of detection (PODs) was 0.961 and 0.874, the probability of false detection (POFDs) was 0.280 and 0.667, and the probability of false alarm (POFAs) was 0.120 and 0.110, respectively, indicating that our model could be used for landslide early warnings.

It can concluded that classifying landslide susceptibility and establishing rainfall thresholds separately not only provides refined thresholds but also avoids overestimation or underestimation of the thresholds for slopes, especially when considering the application to disaster prevention.

The data used in this paper are part of an ongoing project. In this phase, the data are not available online but can be accessed by contacting the corresponding author and Central Geological Survey, MOEA, Taiwan for collaboration purposes.

The authors declare that they have no conflict of interest.

This article is part of the special issue “Landslide early warning systems: monitoring systems, rainfall thresholds, warning models, performance evaluation and risk perception”. It is not associated with a conference.

Acknowledgements

The authors would like to thank the Central Geological Survey, Taiwan, for supporting this research financially and for providing helpful comments on the research. Edited by: Samuele Segoni Reviewed by: Ascanio Rosi and one anonymous referee