Regional Landslide Hazard Assessment Using Extreme Value Analysis and a Probabilistic Physically Based Approach

Abstract

The accurate assessment of landslide hazards is important in order to reduce the casualties and damage caused by landslides. Landslide hazard assessment combines the evaluation of spatial and temporal probabilities. Although various statistical approaches have been used to estimate spatial probability, these methods only evaluate the statistical relationships between factors that have triggered landslides in the past rather than the slope failure process. Therefore, a physically based approach with probabilistic analysis was adopted here to estimate the spatial distribution of landslide probability. Meanwhile, few studies have addressed temporal probability because historical records of landslides are not available for most areas of the world. Therefore, an indirect approach based on rainfall frequency and using extreme value analysis and the Gumbel distribution is proposed and used in this study. In addition, to incorporate the nonstationary characteristics of rainfall data, an expanding window approach was used to evaluate changes in the mean annual maximum rainfall and the location and scale parameters of the Gumbel distribution. Using this approach, the temporal probabilities of future landslides were estimated and integrated with spatial probabilities to assess and map landslide hazards.

1. Introduction

Landslides occur frequently, not only in Korea but also around the world, and cause severe damage to human lives and property. To prevent or reduce damage, injuries, and death caused by landslides, there is a need for landslide hazard analysis, which estimates the probability of a potential landslide occurrence within a given period of time and over a specific area [1,2]. That is, the spatial and temporal probabilities of landslide occurrence should be analyzed to determine landslide hazards. The spatial probability of landslide occurrence, which is also known as landslide susceptibility, predicts where a landslide may occur. A large number of studies have been conducted on landslide susceptibility using a variety of approaches. Landslide susceptibility analyses are generally either statistically or physically based [3,4,5,6]. Statistical approaches acquire knowledge of susceptibility obtained through the statistical analysis of the relationship between landslide occurrences and various conditioning factors [5,7,8,9,10,11,12,13]. However, when applied to a wide area, statistical analysis requires considerable data on both landslide distribution and conditioning factors. In addition, statistical analysis considers the statistical relationship between landslides and conditioning factors exclusively without the consideration of slope failure mechanisms [14]. Therefore, in recent years, physically based analysis, which incorporates the physical processes and mechanisms of landscape occurrence, has been used with a physical slope model to estimate the spatial probability of landslide occurrence independent of its occurrence history [14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32]. This is, therefore, a very promising approach to shallow landslide susceptibility analysis [33]. Moreover, a grid-based analytical structure provides a convenient framework for wide-area coverage in a geographic information system (GIS) environment [34].

In contrast, temporal landslide probability predicts when landslides may occur. Relative to landslide susceptibility analysis, few temporal probability studies have been conducted. In general, temporal probability has been evaluated using the statistical analysis of landslide frequency through long-period multitemporal landslide inventory [35,36,37,38,39,40,41]. However, considerable time and effort are required to construct the necessary multitemporal landslide inventory and, accordingly, these data are not available in most areas. An indirect approach based on the frequency of landslide-triggering events (i.e., earthquake or rainfall events) is, therefore, proposed. In this approach, recognizing that landslides are mainly caused by rainfall events, the temporal probability of such rainfall events is adopted as the temporal probability of landslide occurrence [42,43,44,45,46,47,48,49]. A rainfall threshold for landslide occurrence is determined and then historical rainfall data are analyzed to derive the probability that the rainfall threshold will be exceeded by a certain rainfall event (the exceedance probability). The rainfall exceedance probability is observed as an effective surrogate for temporal landslide probability [37,42,43,45,46,48]. The advantages of this approach are that a complete multitemporal inventory is not required and that temporal probability can be estimated wherever historical rainfall records, which can be easily obtained from rainfall gauges, are available.

In this research, temporal landslide probability was estimated using historical rainfall records. Then, landslide hazard was calculated by combining temporal probability with the spatial probability obtained by conducting a physically based analysis. This approach was adopted in the Jinbu region of Gangwon-do, Korea, where many landslides occurred in July 2006 as a result of extreme rainfall. This approach can also help to solve a global social issue. Achour [50] noted that landslide hazard analysis is one of the significant measures necessary for land use planning and disaster risk reduction, supporting target 15.3 (“By 2030, combat desertification, restore degraded land and soil, including land affected by desertification, drought and floods, and strive to achieve a land degradation-neutral world”) of the United Nations 2030 Sustainable Development Goals (SDGs). Our approach will, therefore, contribute to the achievement of the UN’s SDGs, especially goals 13 (“Take urgent action to combat climate change and its impacts”) and 15 (“Sustainably manage forests, combat desertification, halt and reverse land degradation, halt biodiversity loss”).

2. Materials and Methods

2.1. Study Area

The Jinbu area was selected for testing the proposed approach because it experienced an extreme rainfall event from 14 to 16 July 2006, with numerous landslides being reported (Figure 1). The study area is at latitude 37°33′20″ N–37°39′26″ N and longitude 128°29′49″ E–128°36′36″ E and is mostly mountainous, with an average altitude of about 660 m. The predominant lithological units are Triassic Nokam formation and Imgye granite (Figure 2). These are located on a Precambrian biotite gneiss. The Nokam formation is mostly fine sandstone and sandy shale, and the Jurassic Imgye granite, which occurs as batholiths, consists of granite with syenite and diorite. Ordovician limestone, with a small amount of sandstone and shale, is also distributed across the area. The region was extensively intruded by granitoids during the Daebo Orogeny, which lasted from the Jurassic to the Cretaceous periods [14,51].

To construct a landslide inventory, landslide locations were identified by conducting the comparison of 0.5 m resolution aerial photographs, obtained from the National Geographic Information Institute (Suwon, Korea), taken before and after the 2006 event. A total of 1310 landslides were identified (Figure 1). This is the only landslide record for this study area, given that no other landslide occurrence, either before or after the 2006 event, has been reported. The identified landslides in this area were translational shallow landslides. Their length and width ranged from 30 to 1200 m and from 4 to 20 m, respectively. Their depth to failure plane ranged from 0.5 to 3 m. They started as translational shallow slides and became flow-type landslides as they moved downward. Rainfall data from 1973 to 2017 were obtained from the Sangjinbu rainfall station (latitude 37°39′32″ N and longitude 128°34′41″ E), which is the most reliable rainfall station in the area.

2.2. Evaluation of Temporal Probability

2.2.1. Extreme Value Analysis

In this study, an indirect approach to the evaluation of the temporal landslide probability was adopted. Specifically, the temporal probability of a landslide-triggering rainfall event was evaluated by conducting statistical analysis on historic rainfall data, and the probability of such a rainfall event was adopted as the temporal probability of a landslide occurring. A rainfall threshold, the minimum rainfall required to initiate landslides [48], was first determined. Based on rainfall records and the literature, the rainfall threshold in this area was estimated as 227 mm over a 24 h time period; in July 2006, this threshold was reached and triggered landslides [52]. Once this threshold was determined, its exceedance probability could be calculated. In previous studies, exceedance probabilities have been evaluated using a binomial or Poisson distribution model [37,42,45,46,53,54,55]. However, the use of these models requires an estimate of the mean recurrence interval for the periods between landslide-triggering events. Where no recurrent landslide event has been recorded, such as in this study area, it is impossible to estimate the recurrence interval [56]. Therefore, extreme value analysis, which is able to infer the probabilities of future extremes using past records, was used to evaluate the exceedance probability. This can be applied even in areas where no recurrent landslide-triggering rainfall has been observed. Extreme value analysis is recognized as appropriate for the analysis of the temporal probability of shallow landslides caused by intense rainfall [48]. Therefore, it has been widely used in the context of extreme hydrological events. In extreme value analysis, the maximum rainfall event in a given year, the annual maximum (AM), is considered to follow a generalized extreme value (GEV) distribution [56,57]. Among various GEV distributions, the Gumbel distribution (extreme value type I) has been adopted to estimate the temporal probability of rainfall-induced landslides [39,46,48,49,58,59,60,61,62,63]. In addition, the Gumbel distribution has been applied to determine the frequency of extreme rainfall events in Korea [64]. The cumulative Gumbel distribution is evaluated by the following:

$$F_{GUM}\left(x\right)=exp\left\{-exp\left(-\frac{x-u}{\alpha }\right)\right\},-\infty <x<\infty ,$$

(1)

where $u$ is the location parameter and $\alpha$ is the scale parameter. Subsequently, the exceedance probability of a rainfall event for a specific year is evaluated by the following.

$$p=1-F_{GUM}\left(x\right)=1-exp\left\{-exp\left(-\frac{x-u}{\alpha }\right)\right\}$$

(2)

To evaluate the exceedance probability by adopting the Gumbel distribution, Gumbel parameters (i.e., location and scale parameters) must be estimated. The method of moments, which assumes the equality of population and sample moments, has commonly been applied to estimate the parameters for a given probability model [39,62]. In this study, the method of moments was used to estimate Gumbel parameters.

2.2.2. Nonstationary Approach

Previously, extreme value analysis assumed processes to be stationary, which means that historical rainfall data are invariant over extended time periods. However, increases in extreme rainfall frequency and intensity caused by climate change have been reported by many recent studies. The stationary assumption, with unchanging AM values and Gumbel parameters, is therefore not valid. A nonstationary model should respond to changes in AM rainfall and consequent changes in Gumbel parameters. Zeng et al. [65] proposed an expanding window approach to analyze nonstationarity in AM rainfall. The expanding window begins with a given minimum size at a fixed starting point, but as the analysis progresses into the time series, the window expands to include each new data value rather than only a finite and constant widow size [66]. Initially, a 20-year window of historical data, as suggested previously [67], was used to evaluate the mean AM and Gumbel parameters. Next, the period window was expanded to include 21 years. This procedure was repeated until all the data years were included. In this study, rainfall data from 1973 to 2017 were available, and the mean AM rainfall and Gumbel parameters were initially calculated from 1973 to 1992. Using the expanding window, the means and parameters were then obtained for 1973 to 1993, etc. This approach is able to reveal any nonstationary trend in mean and statistical parameters.

2.3. Evaluation of Spatial Probability

Spatial probability was calculated using a physically based approach. As recent landslides triggered by heavy rainfall are mainly shallow, the infinite slope model, which is widely used for shallow-depth slope failure, was used as the physical model. Previously, the infinite slope model has mainly been used to assess the stability of individual slopes. With the rapid development of GIS-based analysis, the infinite slope model can now be used to analyze shallow landslide susceptibility over broad areas.

The infinite slope model evaluates a factor of safety (FS) based on the concept of limit equilibrium. It analyzes the equilibrium of a potentially unstable mass by comparing the driving force, the force tending to slide along the failure plane, with the resisting force, which is the shear strength. If the groundwater height is assumed to be $z_w$ above the sliding plane (Figure 3), the equation used for calculating FS using the infinite slope model is as follows:

$$\mathrm{FS}=\frac{c+\left(\gamma D-{\gamma }_w{z}_w\right){\cos }^2\alpha \tan \varphi }{\gamma D\sin \alpha \cos \alpha } ,$$

(3)

where c and $\varphi$ are the cohesion and effective friction of the slope materials; $\gamma$ and ${\gamma }_w$ are the unit weights of the slope materials and water, respectively; D is the vertical depth to the failure surface; $z_w$ is the groundwater level; and $\alpha$ is the angle of the slope.

A transient hydrological model was used to estimate $z_w$. The transient infiltration model is a process used for estimating pore pressure changes during rainfall infiltration. This was used in conjunction with a grid-based regional slope stability model (Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability Model, TRIGRS) [68], which estimates shallow landslide occurrence by combining the transient pressure increases caused by rainfall and infiltration [69]. The infiltration model was based on Iverson’s [70] solution, which provides a theoretical background to the influence of hydrologic processes on landslide locations and occurrence times derived using the Richards’ equation [18,27]. This evaluates transient infiltration by modeling pore water pressure. TRIGRS was coupled with Monte Carlo simulation (MCS) for the spatial probability assessment carried out in this study. In MCS, the values of predictive variables are randomly generated according to their probability density functions (PDFs). This technique is widely used for probabilistic analysis because of its robustness and conceptual simplicity. Here, FS values were calculated from sets of these randomly generated input values. After numerous iterations, the failure probability was evaluated from the repeated FS values. This calculation was carried out for all pixels throughout the study area, and the results were mapped as the spatial distribution of landslide probability.

2.3.1. Construction of a Spatial Database of Input Parameters

The physical slope model requires input parameters such as geometric characteristics and strength parameters for the slope materials. Geometric input parameters can be obtained from topographic data and strength parameters can be acquired from field investigation and laboratory tests. A digital elevation model (DEM) was constructed to obtain the geomorphic attributes of the study area using digital topographic maps acquired from the National Geographic Information Institute. Then, a thematic map for the slope angle (Figure 4) was created using DEM. In addition, soil thickness, which is the depth to the failure surface, was acquired from digital soil maps obtained from the National Institute of Agricultural Science (NIAS, Jeollabuk-do, Korea). The obtained map was converted into a grid-based (raster) layer, providing a thematic map of soil thickness with a 10 m resolution (Figure 5).

The calculation of spatial probability using a physically based method (Equation (3)) requires strength parameters (cohesion and friction angle), unit weight, and the hydraulic conductivity of slope materials, since these are indispensable input parameters in physically based analyses [71]. These parameters should be obtained from laboratory tests of soil samples collected from the field. Twenty soil samples were collected from the study area. The geotechnical parameters were obtained from direct shear tests conducted on the slope materials. Other laboratory tests (such as permeability tests) were carried out to obtain the hydraulic conductivities and unit weights of the soils. These were considered as deterministic parameters.

Slope material is composed of soils developed from underlying rocks, and soil properties such as geotechnical and hydrological parameters are strongly affected by the rock type involved. The soil samples used for laboratory tests were, therefore, collected from areas with different underlying rock types. The geotechnical and hydrological parameters listed in

Table 1. Input parameters for the physically based model.

Geological Formation	Friction Angle (°)		Cohesion (kPa)		Unit Weight (kN/m3)	Hydraulic Conductivity (m/h)
	Mean	COV (%)	Mean	COV (%)
Felsite	20.8	15	17.5	30	17.8	0.171
Quartzite in Sambangsan	40.6	15	8.0	30	19.3	0.096
Imgye granite	35.2	15	3.8	30	23.2	0.089
Jeongseon limestone	28.4	15	4.4	30	17.9	0.019
Sandstone in Nokam	40.2	15	10.4	30	18.4	0.090
Sandstone in Gobangsan	37.1	15	7.8	30	18.7	0.084

are linked to the underlying rock type. While the mean values of the strength parameters were calculated using test data, there was substantial uncertainty because the quantity of test data was limited compared with the size of the study area. Therefore, high values for their coefficients of variation (COV) (namely, 30% for cohesion and 15% for friction angle) were used in this analysis. To estimate the groundwater height using the hydrological model, hourly rainfall data were acquired from Sangjinbu automatic weather station (Figure 6). The landslides began around 10:00 on 15 July, and 227 mm (in the previous 24 h) was considered as the landslide-triggering rainfall threshold [52].

2.3.2. Monte Carlo Simulation

In probabilistic analyses, using MCS, cohesion and friction angle, the major sources of uncertainty as a consequence of limited sampling and spatial variability, were considered as random variables. Their statistical properties, as required for MCS, were based on a normal PDF, as previous studies have suggested [21,72,73,74,75,76,77,78,79,80]. Their assigned means and COVs are shown in Table 1. In the MCS process, uniformly distributed random numbers between zero and one were generated, then random values for each input variable were calculated using the generated random numbers corresponding to the cumulative normal distribution function for each input variable. The generated parameters were used, in combination with the deterministic input data, to calculate 5000 individual FS values for each pixel. Then, the failure probability, or the proportion of cases in which FS was less than 1, was determined. This process was conducted for all pixels in the study area to produce a landslide probability distribution map.

3. Results

3.1. Temporal Probability of Landslide Occurrence

The Sangjinbu rainfall time series from 1973 to 2017 were used to determine the nonstationary character of the local rainfall data using the expanding window approach. The maximum 24 h rainfall value for each year was calculated from hourly rainfall data and designated as the AM rainfall for that year. Then, the mean AM values and Gumbel parameters were calculated for the first 20 years of data (1973–1992). By adding rainfall data for each year to the initial 20 years of data, new mean AM and Gumbel parameters were then derived.

Table 2. Mean AM rainfall and the parameters of the Gumbel distribution using an expanding window.

No.	Data Period	Location Parameter	Scale Parameter	Mean of AM Rainfall
1	1973–1992	110.42	50.12	140.4
2	1973–1993	109.38	50.91	138.8
3	1973–1994	106.18	50.40	135.6
4	1973–1995	109.10	50.01	138.5
5	1973–1996	111.02	50.46	140.0
6	1973–1997	113.05	49.37	141.7
7	1973–1998	115.45	49.85	144.0
8	1973–1999	118.18	49.97	147.1
9	1973–2000	116.54	49.98	145.2
10	1973–2001	116.88	49.73	145.1
11	1973–2002	119.60	50.49	148.7
12	1973–2003	120.84	49.87	149.6
13	1973–2004	120.57	50.19	148.9
14	1973–2005	121.34	50.31	149.2
15	1973–2006	123.82	50.77	153.1
16	1973–2007	122.81	50.35	151.9
17	1973–2008	124.10	50.94	152.9
18	1973–2009	126.16	50.33	155.3
19	1973–2010	126.28	50.81	155.0
20	1973–2011	128.40	50.59	157.9
21	1973–2012	125.95	51.43	155.7
22	1973–2013	126.25	50.86	155.6
23	1973–2014	125.42	50.59	154.6
24	1973–2015	122.44	51.02	152.1
25	1973–2016	124.28	51.21	154.4
26	1973–2017	122.87	51.13	153.1

shows the mean AM and location and scale parameters for different time periods under this method.

Linear regression was then used to check for temporal trends in the mean AM and Gumbel parameters to estimate the AM values and Gumbel parameters for future years. Figure 7 shows the results of the linear regression for mean AM. The linear regression equation used was as follows:

$$Mt=\left(0.7677\times N_t\right)-1390.2,$$

(4)

where $Mt$ is the mean AM for the year $N_t$. Using this equation, the mean AM for future target years can be estimated. The Gumbel parameters were estimated, for future years, using linear regressions between the mean AM values and the location and scale parameters (Figure 8 and Figure 9). Equation (5) is the linear regression equation used for estimating location from the mean AM, while Equation (6) is the regression equation used for the scale and mean AM:

$$u_T=0.9863\times Mt-27.077,$$

(5)

$${\alpha }_T=0.0215\times Mt+47.231,$$

(6)

where $u_T$ is the location parameter, $Mt$ is the mean AM, and ${\alpha }_T$ is the scale parameter for the target year. Using these equations, the mean AM and the scale and location parameters of the Gumbel function for each of the next $N_i$ years from the final year of rainfall data, 2017, were estimated for $N_i$ = 10, 50, 100, and 150 (

Table 3. Mean AM and statistical parameters of the Gumbel distribution over the next $N_i$ years.

Period $N_i$ (Years)	Mean AM	Location Parameter	Scale Parameter
10	160.6	131.35	50.69
50	191.3	161.64	51.35
100	229.7	199.50	52.17
150	268.1	237.35	53.00

). From these values, future exceedance probabilities can be calculated using the Gumbel distribution function, and these values can be considered as the temporal probability of landslide occurrence.

Table 4. Temporal probability of landslide occurrence over the next $N_i$ years.

Period $N_i$ (Years)	Temporal Probability
10	0.2197
50	0.3659
100	0.6145
150	0.8575

shows the temporal landslide probabilities over the four analyzed time periods.

3.2. Spatial Probability of Landslide Occurrence

Figure 10 maps the distribution of the spatial landslide probability. A receiver operating characteristics (ROC) graph was used to evaluate model performance. In the analysis, true class (landslide occurrence) is compared with modeled class (landslide prediction) using a confusion matrix [3]. Here, the analyzed grid cells (i.e., modeled class) were classified as unstable or stable for comparison with the landslide occurrence (i.e., true class). In previous studies [27,30,33,77,81,82,83], a landslide probability greater than 10% has been used as the criterion for an unstable area; hence, a grid cell with a probability greater than 10% was classified here as unstable. It was found that 71.3% of the observed landslide locations were classified as unstable. That is, the true positive rate (TPR; the number of correctly predicted landslide pixels over the total number of landslide occurrence pixels) was 0.713. In addition, 27.4% of the nonlandslide pixels were mapped as unstable; that is, the false positive rate (FPR) was 0.274. Model performance was evaluated as 71.9% on the basis of the area under the curve (AUC) value.

3.3. Landslide Hazards

Landslide hazard was evaluated by multiplying the temporal probability of landslide, obtained from extreme value analysis, by spatial probability, which was obtained using a physically based model and MCS. The landslide hazards for four future time periods (10, 50, 100, and 150 years) were calculated (

Table 5. Proportion (%) of landslide hazard values.

Time Period (Years)	Landslide Hazard
	$0–0.1$	$0.1–0.2$	$0.2–0.3$	$0.3–0.4$	$0.4–0.5$	$0.5–0.6$	$0.6–0.7$	$0.7–0.8$	$0.8–0.9$	$0.9–1.0$
10	85.90	9.21	4.89
50	82.02	5.76	5.31	6.91
100	78.02	5.23	3.36	2.98	3.19	5.66	1.56
150	74.47	6.26	3.01	2.37	2.14	2.11	2.30	3.22	4.12

) and are mapped in Figure 11. As expected, the landslide hazard values increased as the time period increased. In the 10-year-period landslide hazard map, no pixels had landslide hazard values of over 0.2. Moreover, for the 10- and 50-year periods, all pixels were less than 0.5, which means that there was a low landslide hazard area. When the landslide hazard maps for 10 and 50 years were compared with the spatial probability map, landslide hazards were substantially lower than the spatial probability of landslide occurrence. This is because the temporal probabilities that were multiplied by the spatial probabilities to obtain the landslide hazards over the 10- and 50-year periods were small: 22.0% and 36.6%, respectively.

The landslide hazard values for the 100- and 150-year periods were, as expected, larger than those in the 10- and 50-year periods. The proportion of high-hazard pixels (hazard value > 0.5) in the 100- and 150-year periods was greater: for 100 years, 7.22%, and for 150 years, 11.75%. This is because of the greater temporal probabilities found for the 100- and 150-year periods: 61.5% and 85.8%, respectively.

4. Discussion and Conclusions

We proposed a process to estimate temporal landslide probability using extreme value analysis and spatial probability using a physically based model. In previous studies, temporal probability was estimated by using frequency analysis of historical landslide occurrences or, indirectly, of rainfall events that triggered landslides. For this, sufficient data on repetitive landslides or recurrent rainfall events are required. However, in many cases it is practically impossible to obtain sufficient data on either landslide occurrence or recurrent landslide-triggering rainfall events. Therefore, this study adopted extreme value analysis to evaluate temporal probability. This approach can be applied in areas where a multitemporal inventory is not available or where a single landslide event has occurred. Extreme value type I distribution, also known as the Gumbel distribution, was used to analyze time series rainfall data and estimate the triggering threshold exceedance probability. Moreover, to accommodate the nonstationary character of rainfall records in a time of climate change, the expanding window method was adopted, and changes in AM rainfall and the Gumbel parameters were estimated. Rainfall data from 1973 to 2017 were used and AM rainfall and Gumbel parameters for four future periods (10, 50, 100, and 150 years hence) were estimated using linear equations derived by the expanding window method. The temporal probabilities of landslide occurrence for the next 10, 50, 100, and 150 years were calculated, and their values were 0.2197, 0.3659, 0.6145, and 0.8575, respectively. Spatial probabilities were evaluated using a physically based approach and the infinite slope model, in conjunction with probabilistic analysis. Input parameters were obtained from laboratory tests and a DEM, with the strength parameters considered as random variables because of their uncertainty and variability. Subsequently, MCS, known as the most complete probabilistic method, was used to account for the uncertainty and estimate the spatial probabilities. Finally, temporal and spatial probabilities were combined to estimate the landslide hazard for future periods.

The proposed approach overcomes the shortcomings of previous studies that determined temporal probability from the frequency analysis of recurrent events. When a multitemporal inventory of historical landslide data is not available or no recurrent landslide events have occurred in an area, as in our case, the existing approach cannot estimate temporal probability. Therefore, an extreme value model (based on the Gumbel distribution) was used to obtain temporal probabilities from available time series rainfall data for the study area. This approach can be used when the conventional approach is impossible. In addition, our approach estimated nonstationary temporal probability, which previous stationary extreme value analyses could not calculate, by using an expanding analytical window. In this manner, the temporal probabilities of landslide occurrence for several different periods were obtained and combined with spatial probabilities, obtained from the probabilistic physically based approach, to evaluate landslide hazard. In previous work, spatial probabilities were estimated using statistical analysis or machine learning methods. However, it has been argued that translating the results of statistical analysis into spatial probabilities may not be appropriate, since the landslide susceptibility index from statistical analysis could not be replaced directly with spatial probability [84]. Therefore, this study used a more appropriate approach by estimating, in combination with probabilistic analysis, physically based spatial probabilities.

However, our proposed approach has some limitations. In this study, only the strength parameters of slope materials were considered as random variables in the MCS in order to limit computational time and effort. However, uncertainty and variability also could be involved in hydrological parameters and the unit weight of soil, and these could also be considered as random variables in future research. In addition, climate change and its influence on landslide occurrence can vary substantially, even over small distances [85]. In our experience, rainfall in this study area has a nonstationary trend, but other rainfall gauges located some distance from the study area suggest a more stationary character. Therefore, it is critical to carefully scrutinize stationarity in rainfall data. Finally, since information about elements at risk and their vulnerability was not available for our study area, it was not possible to fully assess landslide risks—that is, the accurate assessment of threat to life and property. To reduce or prevent the damages and fatalities caused by landslide occurrences effectively, landslide risk should be evaluated in future studies.