Landslide Susceptibility Mapping Using Statistical Methods along the Asian Highway, Bhutan

Abstract

In areas prone to frequent landslides, the use of landslide susceptibility maps can greatly aid in the decision-making process of the socio-economic development plans of the area. Landslide susceptibility maps are generally developed using statistical methods and geographic information systems. In the present study, landslide susceptibility along road corridors was considered, since the anthropogenic impacts along a road in a mountainous country remain uniform and are mainly due to road construction. Therefore, we generated landslide susceptibility maps along 80.9 km of the Asian Highway (AH48) in Bhutan using the information value, weight of evidence, and logistic regression methods. These methods have been used independently by some researchers to produce landslide susceptibility maps, but no comparative analysis of these methods with a focus on road corridors is available. The factors contributing to landslides considered in the study are land cover, lithology, elevation, proximity to roads, drainage, and fault lines, aspect, and slope angle. The validation of the method performance was carried out by using the area under the curve of the receiver operating characteristic on training and control samples. The area under the curve values of the control samples were 0.883, 0.882, and 0.88 for the information value, weight of evidence, and logistic regression models, respectively, which indicates that all models were capable of producing reliable landslide susceptibility maps. In addition, when overlaid on the generated landslide susceptibility maps, 89.3%, 85.6%, and 72.2% of the control landslide samples were found to be in higher-susceptibility areas for the information value, weight of evidence, and logistic regression methods, respectively. From these findings, we conclude that the information value method has a better predictive performance than the other methods used in the present study. The landslide susceptibility maps produced in the study could be useful to road engineers in planning landslide prevention and mitigation works along the highway.

1. Introduction

A landslide is defined as “the movement of a mass of rock, earth, or debris down a slope”, and they are classified according to the type of slope movement (fall, topple, spread, flow, slide), type of material involved (rock, earth, debris), and the speed of movement [1,2]. Even though gravity is the essential contributor to the occurrence of a landslide event, other triggering factors, such as an earthquake, rainfall, flood, or human intervention [3], significantly increase the likelihood of landslide occurrence. Landslides constitute a major geological hazard and pose considerable risks to the livelihood and lives of the population living in and around the affected area [4]. Prolonged disruption of transportation networks, loss of fertile land, collapse and submergence of buildings, loss of life, etc. are some of the risks associated with a landslide event that can translate into major social impacts and economic loss. Expansion of human settlement into geologically sensitive areas, infrastructure development, and increased agricultural practices result in land use changes that further aggravate the problem of landslides and associated risks [5,6,7]. Cutting slopes for infrastructure development, particularly during road construction, is a major triggering factor for most landslides. Conversely, landslides impede socio-economic activities, such as the development of efficient transportation networks, reservoirs, settlement areas, and agricultural fields, especially in mountainous regions [8]. Hence, to support sound decision-making in building up a national socio-economic development plan that addresses the impact of landslides, information based on risk analysis and landslide assessment concerning the likelihood of landslide occurrences are useful. Risk analysis and landslide assessment are crucial in the development of mitigation and disaster preparedness plans [8]. Landslide assessment is generally considered in terms of landslide susceptibility (“the potential for a given slope to fail compared to others”), landslide hazard (“the potential posed by a landslide to cause damage or loss”), and landslide risk (“the actual or potential damage or loss that may occur as a result of a landslide”) [9,10]. To effectively mitigate landslide risk and prevent landslide hazards, a dependable and detailed landslide susceptibility map (LSM) must be developed for the region [11,12,13] to provide key information to a range of end users. An LSM may be developed for specific applications, such as landfill zoning [14], road corridors [15,16,17,18,19], land use planning [20,21], and reservoir basins [22].

Landslide susceptibility mapping considers many causal factors, which are usually presented as thematic layers in a geographic information system (GIS) platform. Various models and methodologies have been utilized to decide the impact of causal factors on landslide occurrence, of which the multi-criteria decision analysis (MCDA) method based on fuzzy logic [23], analytical hierarchy process (AHP) [24,25,26], and weighted linear combination [27] are most commonly used. In terms of statistical methods, the frequency ratio method is the simplest and easiest to perform, while the information value (IV) [28,29] and weight of evidence (WOE) [30,31] methods are useful in determining the impact of causal factor class on landslide occurrence. Logistic regression (LR) is also used by many researchers for determining the weight of the causal factors [17,26,32,33]. The application of machine learning, such as artificial neural networks (ANN) [34,35,36] and support vector machines (SVM) [37], has also been considered a promising technique for LSMs. Various researchers have studied the comparative performance of these methods. The LR model was found to be more suitable than the frequency ratio method in some studies [12,38], whereas some studies found out that the frequency ratio method is better than the certainty factor [39] and LR [40]. WOE was also established as comparatively less accurate than the fuzzy logic technique [41]. Some methods have been modified to enhance performance in predicting future landslides. Ba et al. [42] improved the IV model by integrating it with an AHP and Gray clustering algorithms. The AHP technique was modified by adopting an interval matrix in deriving the optimal decision matrix [43]. However, only a few studies have developed LSMs for areas along road corridors, and no studies have carried out comparative analysis using the information value (IV), weight of evidence (WOE), or logistic regression (LR) model methods.

In Bhutan, the process of road building in the mountainous Bhutanese terrain has aggravated the occurrence of landslides, affected economic growth, and caused loss of lives [44]. Furthermore, the road corridors in Bhutan present a unique condition in that the anthropogenic impacts on the environment are largely due to road construction and very rarely due to sparse human settlement. The objective of this study, therefore, was to investigate the suitability of various methods for developing a landslide susceptibility map in Bhutan along the highways. We derived the relative weights of different classes of landslide causal factors using the IV, WOE, and LR methods, and compared the predictive performance of the methods in generating an LSM.

2. Study Area

For this paper, a road corridor of 80.9 km, averaging 5 km in width, along the Asian Highway-48 (AH-48) stretch from Phuentsholing to Chukha Thegchhen zam (bridge), covering an area of 237.28 km${}^2$, was considered (Figure 1a).

The AH-48 is a major trade route connecting the capital city of Thimphu and rest of Bhutan to India, its main economic partner. The study area is located between 26${}^{\circ }$50${}^{\prime }$34${}^{″}$ N and 27${}^{\circ }$5${}^{\prime }$20${}^{″}$ N latitude and 89${}^{\circ }$23${}^{\prime }$29${}^{″}$ E and 89${}^{\circ }$33${}^{\prime }$12${}^{″}$ E longitude. The highway ascends in altitude from 216 to 2159 m above Mean Sea Level (MSL) over a distance of 40 km from Phuentsholing (Figure 1b) and runs through a geological formation consisting of moderate to highly weathered phyllites [44,45]. During the monsoon season with an average rainfall of 1663.4 mm per year [46], this section of road is subjected to frequent landslides of varying magnitude at a number of locations (Figure 2).

Often, the landslides at a site along the road corridor are shallow in nature and occur mainly during the monsoon season [47]. These landslides are caused due to deforestation and toe cutting on the fragile slope during the road construction process. The landslides at Sorchen and Jumja [48] are most critical, necessitating realignment of the road length as an avoidance strategy in 1992 [44,49]. In 2017, the highest 24-h rainfall of 285.4 mm in the study area [50] triggered a major landslide at three sites, leading to road closure for many days. This affected the transport of goods and people between Phuentsholing and other parts of Bhutan, causing major economic losses and disrupting everyday life.

3. Data

3.1. Landslide Inventory

Landslide distribution—or inventory mapping—is the fundamental information required in determining the size and features of a landslide [3,51]. Since a landslide inventory in Bhutan is not available, the technical report of the Bhutan Land Cover Assessment 2010, National Soil Services Center (NSSC), and Policy and Planning Division (PPD) [52] was used as the primary guide for field visits, as well as digitized satellite imagery from the publicly available data source Google Earth. From this report, in which a class of “degraded land” with subclass “landslide” showed landslide areas, and by interpreting Google Earth images, 120 landslides totaling an area of 2.77 km${}^2$ were identified. This was obtained by digitizing landslides ranging from 122 to 304,625 m${}^2$ area in the Google Earth environment and verified with news reports and field visit. Sixty percent of landslides by area were aligned with the Ministry of Agriculture and Forests, Royal Government of Bhutan (MoAF, RGoB) [52] report, while the rest were obtained from Google Earth images, which were likely more recent and active landslides areas. Figure 3 shows the distribution of landslides in the study area.

3.2. Causative Factor Selection

In the study, the slope failure susceptibility due to underlying causative factors [8] was considered. Land cover, lithology, elevation, proximity to roads, drainage, fault lines, slope aspect, and slope angle were considered as causal factors based on literature [53] and the availability of data for the study area. Each causal factor was mapped and divided into the several equal interval classes described in the legend in Figure 4.

The percentage of landslides against each class of causal factor is illustrated in Figure 5. The geomorphic factors of elevation, slope aspect, and slope angle were derived from corrected DEM of 30 m resolution obtained from the Ministry of Work and Human Settlement, Bhutan. The highway altitude climbs from 216 to 2159 m above MSL over a distance of 40 km (Figure 1b) with diverse climatic conditions, hydrology, and geology. Hence, an elevation map with eight classes of 300 m intervals was produced. The slope aspect indicates the saturation of the slope with water and heat, which affects soil, rock, and vegetation types [11]. The south-facing aspect has a higher frequency (30%) of landslides in the nine proposed classes. The angle of slope in degrees was divided into six classes of equal intervals of 10${}^{\circ }$, with 26.59% of landslides falling in the slope angle class of 30${}^{\circ }$–40${}^{\circ }$. The lithological and fault details were digitized from a 1:50,000 geological map of Bhutan from Long et al. [45].

The study area includes eleven classes of lithology under the three main zones of the greater Himalayan zone, lesser Himalayan zone, and Paro formation, and seven classes of distance to faults. The lithology consists mainly of amphibolite, quartzite, schist, slate, phyllite, marble, paragneiss, limestone, dolostone, and granite dating to the Paleoproterozoic to Ordovician Era [45]. The details of the lithology and the ages of geological formations are given in

Table 1. Lithology and age of geological formations in the study area [45].

	Formation and Unit	Age	Lithology
GHIo	Greater Himalayan Zone: Orthogenesis unit	Cambrian– Ordovician	Granite, Paragneiss, Schist, quartzite
GHImu	Greater Himalayan Zone: Upper metasedimentary unit	Neoproterozoic–Ordovician	Amphibolite, Paragneiss, Quartzite, Schist
Pzpm	Paro Formation: Middle unit	Cambrian–Ordovician	Quartzite, Marble 10 m thick
Pzpm2	Paro Formation: Middle unit (100–200 m thick)	Cambrian–Ordovician	Quartzite, Marble 100–200 m thick
GHIml	Greater Himalayan Zone: Lower metasedimentary unit	Neoproterozoi–Cambrian	Amphibolite, Quartzite, Schist, Paragneiss
Pzph	Lesser Himalayan Zone: Phuentsholing Formation (Baxa Group)	Age range uncertain: Neoproterozoic or younger	Slate, Phyllite, Limestone
pCp	Lesser Himalayan Zone: Pangsari Formation (Baxa Group)	Age range uncertain: Mesoproterozoic–Cambrian	Phyllite, Dolostone, Marble
pCo	Lesser Himalayan Zone: Orthogneiss (Daling-Shumar Group)	Paleoproterozoic	Granite
pCs	Lesser Himalayan Zone: Shumar Formation (Daling- Shumar Group)	Paleoproterozoic	Quartzite, Schist, Phyllite
pCd	Lesser Himalayan Zone: Daling Formation (Daling-Shumar Group)	Paleoproterozoic	Schist, Phyllite
pZj	Lesser Himalayan Zone: Jaishidanda Formation	Neoproterozoic–Ordovician	Schist, Quartzite

. The number of landslides is relatively higher in the Phuentsholing formation (Pzph) region, where the more commonly found lithologies are slate and phyllite. The AH-48 highway is routed over three active main faults: Shumar thrust (ST), Main Central thrust (MCT), and the Southern Tibetan Detachment (STD), plus a few other minor fault lines and folds [45,54]. Landslide density is significantly higher closer to fault lines, with 35% of overall landslides being within 0–500 m of them. Furthermore, one of the major contributors to landslides is land cover [55,56]. Vegetation aids in protecting a slope whereas bare soil or sparse vegetation agitates the occurrence of a landslide [57,58]. Land cover maps from the MoAF/RGoB [52] report were used to derive seven simplified broad land cover categories to meet the aim of the study. Traffic intensity and the cutting of steep slopes during road construction are another factor influencing landslide occurrence [59,60,61]. The proximity of drainage streams also results in saturating the area and causing landslides [62]. Drainage and road maps were acquired from the Bhutan Geospatial portal website [63], and thematic layers were created with six equidistant buffers of 100 m each.

4. Methodology

An inventory of 120 landslide areas and various factors, such as land cover, lithology, elevation, proximity to roads, drainage, and fault lines, slope aspect, and slope angle, was created and converted into 30 × 30 m grid cells in ArcGIS 10.6 to suit the DEM resolution. An 80% training sample equally proportioned between the landslide and non-landslide pixels was generated randomly and was used to create the spatial associations between occurrences of landslides and the causal factors [64,65,66]. The control sample of 20% was used to validate and verify (

Table 2. Number of pixel cell samples for the study area.

	Training Sample (80%)	Validation Sample (20%)	Total
Pixels with Landslide	2159	522	2681
Pixels without Landslide	208,524	52,360	260,884
Total	210,683	52,882	263,565

) the accuracy of the models. The three methodological approaches of IV, WOE, and LR were used to derive the relationship between causal factors and landslide occurrence. The landslide susceptibility indices (LSIs) were generated to produce an LSM. LSIs were divided into the five classes of very low, low, moderate, high, and very high susceptibility using the Jenks natural break classification method [24,53,67]. The Jenks natural break classification method is used to determine the arrangement of values into different classes by minimizing and maximizing each class’s deviation from the class mean and other groups’ means, respectively [12,42,68]. To validate the performance of the methods, the area under the curve (AUC) of the receiver operating characteristic (ROC) was calculated. The control sample was overlaid on the LSM to examine the predictive capability for future landslides. Additionally, a correlation test was performed to assess the level of similarity between LSMs produced using the models. IBM Statistical Product and Service Solutions 25 (SPSS 25) was used for data management and validation. The flowchart in Figure 6 shows the data sources, thematic layers, and methodology applied in the study.

4.1. Information Value Method

The IV model, a simple statistical method for mapping areas vulnerable to landslides by determining the influence of each class of causal factors on landslide occurrence, was found suitable in a number of studies and has been implemented in numerous landslide hazard assessment studies [28,29,42,69]. In this model, the information value $I_i$ for a class i in a thematic layer considering 80% as the training sample is given by:

$$I_i=log\frac{S_i/N_i}{S/N},$$

(1)

where $S_i$ represents the number of pixels of the class containing a landslide, $N_i$ represents the total number of pixels in the class, S is the total number of pixels with a landslide in the layer, and N is the total number of pixels in the layer.

After deriving the information values for each class of causal factor, raster maps were overlaid in a GIS environment. The LSIs were determined by averaging the information value of the causal factors as:

$$LSI=\frac{1}{M}\sum _{n=1}^M{I}_i,$$

(2)

where M is the number of factors considered. The total information value of factors contributing to landslide occurrence is obtained from Equation (2). The influence of these factors on landslide occurrence is lesser if the LSI value is lower, and vice versa.

4.2. Weights of Evidence Model

The WOE method based on Bayes’ theorem was used to find the spatial association between the location of a landslide and a set of contributing factors. Numerous studies have been conducted using WOE, emphasizing its theoretical background and application. A detailed mathematical description and formulation can be found in the literature [61,70,71]. The weights of the causal factor for the presence and absence of a landslide are determined by:

$$W^{+}=ln\frac{P\left(B\right|D)}{P\left(B\right|\overline{D})}$$

(3)

$$W^{-}=ln\frac{P\left(\overline{B}\right|D)}{P\left(\overline{B}\right|\overline{D})},$$

(4)

where $W^{+}$ and $W^{-}$ are the WOE when a binary predictor factor is present or absent, respectively, P is the probability, B is the presence of the desired class of landslide causal factor, $\overline{B}$ is the absence of the desired class of landslide causal factor, D is the presence of a landslide, and $\overline{D}$ is the absence of a landslide.

A simplified form of Equation (3) above can be expressed as [41]:

$$W^{+}=ln\frac{L{S}_{in}\%}{nonL{S}_{in}\%}$$

(5)

$$W^{-}=ln\frac{L{S}_{out}\%}{nonL{S}_{out}\%},$$

(6)

where $L{S}_{in}\%$ and $nonL{S}_{in}\%$ are the proportions of landslide pixels and non-landslide pixels in the class, respectively. $L{S}_{out}\%$ and $nonL{S}_{out}\%$ are the proportions of landslide pixels and non-landslide pixels outside the same class, respectively. $W^{+}$ implies a positive correlation and $W^{-}$ indicates a negative correlation between the causal factor and the occurrence of a landslide [40,41,70]. The weight contrast is given by $C=(W^{+}-W^{-})$, and its magnitude represents the measure of spatial association between the class of the causal factor and the occurrence of a landslide [72,73,74]. The LSI is derived by aggregating the contrasts of all the factors, with higher values of LSI indicating a higher likelihood of landslide occurrence.

4.3. Logistic Regression Model

An LR model is a regression analysis technique used to determine the likelihood of a binary dependent variable from several independent variables. It is used to predict the presence of an outcome based on the values of predictor variables. The method does not require normally distributed data, and the variable can be continuous, discrete, or any combination of both types [75].

In order to perform the LR method on a sample, it is necessary to check for collinearity between the variables [40,76]. If the tolerance (TOL) <0.1 and variance inflation factor (VIF) >10, this indicates multicollinearity between independent variables [77]. The dependent variable (landslide variable) and the independent variable (causal factor) were modeled using the LR application in SPSS 25 to determine the coefficient of each factor. The 80% training sample consisting of an equal proportion of landslide pixels and causal factor pixels was imported into SPSS. The frequency ratio of each class of causal factor was derived as the percentage of landslides against the percentage of the area of the class. The logistic function was applied to causal factors to constrain the values between 0 and 1, where zero indicates the probability of 0% landslide occurrence and one indicates the probability of 100%, according to the equation:

$$P=\frac{1}{(1+ex{p}^{-z})},$$

(7)

where P is the probability of landslide occurrence, and z is the landslide causal factors assumed to be a linear combination of the causal factors $X_i(i=1,2,3,\dots ,n)$, where

$$Z=B_0+B_1{X}_1+B_2{X}_2+B_3{X}_3\dots +B_n{X}_n,$$

(8)

and $B_i$ is the regression coefficient of landslide causal factors.

4.4. Validation of the Models

The AUC-ROC graphically illustrates the performance of a binary classifier for the false positive rate (1-specificity) against the true positive rate (sensitivity). The AUC represents the performance of the success rate and prediction of the model against the occurrence of a landslide [78]. A landslide predicted in an existing landslide area is a true positive outcome, whereas a landslide predicted in a non-existing landslide area is a false positive [79]. Rasyid et al. [12] defined the diagnostic values as “0.5–0.6 (fail), 0.6–0.7 (poor), 0.7–0.8 (fair), 0.8–0.9 (good), and 0.9–1.0 (excellent)” in the AUC curve. The AUC success rate of each model was derived for the LSI computed for each model and the training samples, while the rate of prediction was verified using a control sample. In addition, the training and control samples were overlaid on the LSMs to assess predictive performance. The training sample covered a small area of higher susceptibility classes, whereas the landslide data would lie in higher classes when overlaid on the LSM [12,66]. Correlation coefficients ranging from 0 (no correlation) to 1 (strong relationship) were derived from the correlation analyses between a pair of LSMs developed from the IV, WOE, and LR models.

5. Results

Three LSMs were produced using the IV, WOE, and LR models with 80% training samples. The detailed calculations and values of each class of causal factors derived using the three methods are given in

Table A1. Information value, weight of evidence, and logistic regression indices for landslide causal factors.

		No. of Landslide Pixels in Class $\left({\mathit{S}}_{\mathit{i}}\right)$	Total No. of Pixels in Class $\left({\mathit{N}}_{\mathit{i}}\right)$	Landslide Ratio Area	Total Ratio Area	Information Value	Frequency Ratio (FR)	Standar- Dized FR	${\mathit{W}}^{+}$	${\mathit{W}}^{-}$	Weight of Evidence (C)
Slope angle (degree)	0–10	96	10,750	4.447	5.102	−0.060	0.871	0.071	−0.138	0.007	−0.144
	10–20	405	44,423	18.759	21.085	−0.051	0.890	0.079	−0.117	0.029	−0.146
	20–30	556	76,822	25.753	36.463	−0.151	0.706	0.000	−0.348	0.156	−0.504
	30–40	574	54,488	26.586	25.863	0.012	1.028	0.138	0.028	−0.010	0.037
	40–50	374	19,257	17.323	9.140	0.278	1.895	0.509	0.639	−0.094	0.734
	>50	154	4943	7.133	2.346	0.483	3.040	1.000	1.112	−0.050	1.162
Aspect	Flat	2	294	0.093	0.140	0.000	0.664	0.355	−0.410	0.000	−0.410
	North	2	17,654	0.093	8.379	−1.956	0.011	0.000	−4.505	0.087	−4.591
	North East	88	27,963	4.076	13.273	−0.513	0.307	0.161	−1.181	0.101	−1.281
	East	168	30,236	7.781	14.351	−0.266	0.542	0.289	−0.612	0.074	−0.686
	South East	450	31,550	20.843	14.975	0.144	1.392	0.751	0.331	−0.072	0.402
	South	655	34,561	30.338	16.404	0.267	1.849	1.000	0.615	−0.182	0.797
	South West	494	30,132	22.881	14.302	0.204	1.600	0.864	0.470	−0.105	0.575
	West	238	21,278	11.024	10.100	0.038	1.091	0.588	0.088	−0.010	0.098
	North West	62	17,015	2.872	8.076	−0.449	0.356	0.187	−1.034	0.055	−1.089
Elevation (m)	0–300	108	4796	5.002	2.276	0.342	2.197	0.610	0.787	−0.028	0.816
	300–600	554	15,863	25.660	7.529	0.533	3.408	1.000	1.226	−0.218	1.444
	600–900	254	12,265	11.765	5.822	0.306	2.021	0.553	0.704	−0.065	0.769
	900–1200	166	17,450	7.689	8.283	−0.032	0.928	0.201	−0.074	0.006	−0.081
	1200–1500	260	33,634	12.043	15.964	−0.122	0.754	0.145	−0.282	0.046	−0.328
	1500–1800	316	39,074	14.636	18.546	−0.103	0.789	0.156	−0.237	0.047	−0.284
	1800–2100	338	35,574	15.655	16.885	−0.033	0.927	0.200	−0.076	0.015	−0.090
	>2100	163	52,027	7.550	24.694	−0.515	0.306	0.000	−1.185	0.205	−1.390
Distance to road (m)	0–100	368	31,726	17.045	15.059	0.054	1.132	0.536	0.124	−0.024	0.148
	100–200	300	23,159	13.895	10.992	0.102	1.264	0.780	0.234	−0.033	0.268
	200–300	265	18,694	12.274	8.873	0.141	1.383	1.000	0.324	−0.038	0.363
	300–400	158	16,087	7.318	7.636	−0.018	0.958	0.216	−0.042	0.003	−0.046
	400–500	147	14,199	6.809	6.740	0.004	1.010	0.312	0.010	−0.001	0.011
	>500	921	106,818	42.659	50.701	−0.075	0.841	0.000	−0.173	0.151	−0.324
Distance to drainage (m)	0–100	1107	91,695	51.274	43.523	0.071	1.178	1.000	0.164	−0.148	0.312
	100–200	502	60,914	23.252	28.913	−0.095	0.804	0.027	−0.218	0.077	−0.295
	200–300	281	29,545	13.015	14.023	−0.032	0.928	0.350	−0.075	0.012	−0.086
	300–400	135	14,204	6.253	6.742	−0.033	0.927	0.348	−0.075	0.005	−0.081
	400–500	78	7440	3.613	3.531	0.010	1.023	0.597	0.023	−0.001	0.024
	>500	56	6885	2.594	3.268	−0.100	0.794	0.000	−0.231	0.007	−0.238
Distance to fault lines (m)	0–500	773	39,183	35.804	18.639	0.285	1.921	1.000	0.653	−0.237	0.890
	500–1000	443	37,506	20.519	17.842	0.062	1.150	0.547	0.140	−0.033	0.173
	1000–1500	160	26,458	7.411	12.586	−0.228	0.589	0.217	−0.530	0.058	−0.587
	1500–2000	289	20,786	13.386	9.888	0.133	1.354	0.667	0.303	−0.040	0.342
	2000–2500	257	18,081	11.904	8.601	0.143	1.384	0.684	0.325	−0.037	0.362
	2500–3000	103	12,317	4.771	5.859	−0.088	0.814	0.350	−0.206	0.011	−0.217
	>3000	126	55,885	5.836	26.585	−0.657	0.220	0.000	−1.516	0.249	−1.765
Lithology	GHIo	43	16,863	1.996	8.012	−0.604	0.249	0.000	−1.390	0.063	−1.453
	GHImu	231	13,688	10.724	6.503	0.217	1.649	0.462	0.500	−0.046	0.546
	Pzpm	172	12,505	7.985	5.941	0.128	1.344	0.362	0.296	−0.022	0.318
	Pzpm2	46	4079	2.136	1.938	0.042	1.102	0.282	0.097	−0.002	0.099
	GHIml	429	93,368	19.916	44.360	−0.348	0.449	0.066	−0.801	0.364	−1.165
	Pzph	471	14,049	21.866	6.675	0.515	3.276	1.000	1.187	−0.178	1.364
	pCp	277	13,127	12.860	6.237	0.314	2.062	0.599	0.724	−0.073	0.797
	pCo	36	3247	1.671	1.543	0.035	1.083	0.276	0.080	−0.001	0.081
	pCs	215	8727	9.981	4.146	0.382	2.407	0.713	0.879	−0.063	0.941
	pCd	150	22,637	6.964	10.755	−0.189	0.647	0.132	−0.435	0.042	−0.476
	pZj	84	8189	3.900	3.891	0.001	1.002	0.249	0.002	0.000	0.002
Land cover	Forests	895	172,299	41.454	81.781	−0.295	0.507	0.000	−0.679	1.167	−1.847
	Meadows	33	4416	1.528	2.096	−0.137	0.729	0.006	−0.316	0.006	−0.322
	Cultivated agricultural land	68	10,104	3.150	4.796	−0.183	0.657	0.004	−0.420	0.017	−0.438
	Bare areas	711	1900	32.932	0.902	1.562	36.517	1.000	3.598	−0.390	3.988
	Shrubs	365	14,685	16.906	6.970	0.385	2.425	0.053	0.886	−0.113	0.999
	Water bodies	57	2366	2.640	1.123	0.371	2.351	0.051	0.855	−0.015	0.870
	Built-up areas	30	4913	1.390	2.332	−0.225	0.596	0.002	−0.518	0.010	−0.527

* S/N = 2159/210, 683 = 0.0102.

. In each of the models, the LSIs of the classes under different causal factors corresponding to each pixel in the map were divided into five susceptibility levels and mapped as shown in Figure 7.

5.1. Information Value Method

The information values of each class under causal factors were determined using Equation (1) and by overlaying the causal factor map on the landslide distribution map. As indicated in Table A1, the IV is directly proportional to the slope angle. A slope angle of >50 degrees had the highest IV value of 0.483 in the class. All southern exposure aspects had higher IV in classes, with a maximum IV of 0.267 on the southern aspect. The 300–600 m elevation class had the highest IV, followed by the 0–300 m and 600–900 m elevation classes. The least IV value was at a higher altitude of >2100 m. The highest IV for the distance to road factor was 0.141 in the 200–300 m class, with classes over 300 m having lower IV values. In contrast to distance to road, the distance to drainage factor showed no specific relationship of its proximity with the occurrence of a landslide. However, the IV value for the closest distance class, 0–100 m, was 0.071, which indicates a higher likelihood of occurrence of landslides within the class, but all other classes had insignificant IVs. As shown in Table A1, the proximity to a fault line suggests a greater likelihood of landslide occurrence. The 0–500 m class had the highest IV of 0.285, whereas the >3000 m class had the lowest value of −0.657. For lithology factors, Pzph had the highest IV of 0.515, followed by PCs with 0.382 and pCp with 0.314, suggesting a greater probability of landslide occurrence in these classes. Among all factor classes, the land cover bare area class had the highest IV of 1.562, and, as expected, the forest and built-up area class had the minimum IV. To create the LSI, the IV of the classes under different causal factors corresponding to each pixel in the map was derived using Equation (2).

5.2. Weights of Evidence Model

In this model, the weights and contrast were determined using Equations (5) and (6) and are shown in Table A1. Similarly to the findings with IV, the highest contributors to landslides are the degree of slope of >50 degree (WOE = 1.162), south exposure aspect (WOE = 0.797), 300–600 m elevation (WOE = 1.444), 200–300 m distance to road (WOE = 0.363), 0–100 m distance to drainage (WOE = 0.312), 0–500 m distance to fault lines (WOE = 0.890), lithology group Pzph (WOE = 1.364), and bare area in land cover (WOE = 3.988). The LSI was computed by totaling the WOE contrasts.

5.3. Logistic Regression Model

The TOL and VIF calculated for this study (

Table 3. Multicollinearity indices for causal factors.

	Collinearity Statistics
	Tolerance	VIF
Slope angle	0.953	1.050
Aspect	0.980	1.021
Elevation	0.414	2.413
Distance to road	0.966	1.035
Distance to drainage	0.990	1.010
Distance to fault lines	0.856	1.168
Lithology	0.396	2.527
Land cover	0.982	1.018

Dependent variable: landslide.

) were more than 0.39 and less than 2.6, respectively, indicating that there was no multicollinearity between variables. We could, therefore, use all the variables for LR analysis. The forward step-wise LR analysis shows that all factors with an estimated logistic coefficient have a significance value of less than 0.05. This indicates that all the independent variables have an influence on the landslide occurrence variable (

Table 4. Coefficients of each factor derived from logistic regression.

	B	S.E.	Wald	df	Sig.	Exp (B)	95% C.I. for EXP(B)
							Lower	Upper
Slope angle	1.419	0.092	239.425	1	0.000	4.131	3.452	4.944
Aspect	1.518	0.083	337.914	1	0.000	4.564	3.882	5.366
Elevation	0.487	0.115	18.005	1	0.000	1.628	1.300	2.039
Distance to road	0.483	0.066	54.081	1	0.000	1.622	1.426	1.845
Distance to drainage	0.467	0.056	69.624	1	0.000	1.596	1.430	1.781
Distance to fault lines	1.137	0.074	237.883	1	0.000	3.118	2.698	3.602
Lithology	1.435	0.112	163.460	1	0.000	4.200	3.371	5.234
Land cover	3.882	0.059	4268.221	1	0.000	48.522	43.188	54.515
Constant	−7.896	0.097	6693.400	1	0.000	0.000

B = logistic coefficient; S.E. = standard error of estimate; Wald = Wald chi-square; df = degree of freedom; Sig. = significance; EXP(B) = exponentiated coefficient; 95% C.I. for EXP(B) = 95% confidence interval for EXP(B).

), and, therefore, all causal factors were considered for analysis.

The highest regression coefficient was computed for land cover at 3.882, followed by aspect, lithology, slope angle, distance to fault lines, elevation, distance to road, and distance to drainage. The Z value or LSI was computed in a GIS environment using Equation (8), as follows:

$$\begin{array}{c}\hfill Z=-7.896+1.419\times slopeangle+1.528\times aspect+0.487\times elevation+0.483\times distancetoroad+\\ \hfill 0.467\times distancetodrainage+1.137\times distancetofault+1.435\times litology+3.882\times landcover.\end{array}$$

(9)

5.4. Validation

For the IV, WOE, and LR models, the AUC-success rates were 0.873, 0.872, and 0.866, while the prediction rates were 0.883, 0.882, and 0.880, respectively. From this, we conclude that all the models are suitable methods for generating LSMs (Figure 8), although the IV model shows a slightly better predictive performance.

Figure 9 shows that the overlaid training sample covered 33.2%, 27%, and 16.2% of the study region for IV, WOE, and LR, respectively, in the higher susceptibility classes. When overlaying the training landslide and control landslide samples on the LSM, the area coverage above high susceptibility was 87.4% and 89.3% in the IV model, 84.1% and 85.6% in the WOE model, and 69.6% and 72.2% in the LR model. The IV and WOE models appear to be more reliable, since their values have lesser disparity. From the correlation analysis, the correlation coefficients at a significance level of 0.01 (two-tailed) were 0.699 or greater, with the highest coefficient of 0.845 between WOE and LR (

Table 5. Correlation between the methods.

Methods	Information Value	Weight of Evidence	Logistic Regression
Information Value	1	0.755 **	0.699 **
Weight of Evidence	0.755 **	1	0.845 **
Logistic Regression	0.699 **	0.845 **	1

** Correlation is significant at the 0.01 level (two-tailed).

).

6. Discussion

LSMs were generated from the IV, WOE, and LR models considering the relationship between causal factors and landslides. The LSMs were compared for predictive performance using AUC-ROC by overlaying the training and control samples on LSMs and performing correlation coefficient tests. Higher susceptibility was observed in the extreme southern and northern parts of the focus area, where human settlements, steeper slopes, and more fault lines can be found. Factors such as roads and landcover were found to have a high correlation with the occurrence of landslides, as found by [80]. From the LR model, it was found that the land cover contributes the most to landslide occurrence out of all the factors, followed by aspect, lithology, slope angle, distance to fault lines, elevation, distance to road, and distance to drainage. As shown in Table A1, the most significant contributors to landslide occurrence from each factor were slope angle of >50, southern exposure aspect, 300–600 m elevation, 200–300 m distance to road, 0–100 m distance to drainage, 0–500 m distance to fault lines, lithology group Pzph (Phuentsholing formation), and bare land cover.

The northern region of the study area had the greatest slope angle, and the lowest slope angle was in the south (Figure 4). Landslide occurrence was directly proportional to the slope angle, which agrees with the general conclusion of expecting fewer landslides on gentler slopes because of lower associated shear stresses [11]. This also conforms to the findings of other research [81]. The south, southeast, and southwest had the highest LSIs in aspect classes in the present study. This may be due to the high humidity in the south-facing aspect [11]. Field knowledge also dictates that south-facing slopes are warm, wet, and forested. Most of the study area falls within the elevation class of >2100 m above mean sea level, which is less susceptible to landslide occurrence. In our study, we observed that elevation classes 300–600 m and 0–300 m had a higher influence on landslide occurrence, which gradually decreased as elevation increased. Other studies have also demonstrated similar findings [82]. The stability and resistance to weathering processes of rocky cliffs at higher altitudes may explain this phenomenon [60]. All classes with a distance less than 300 m from the road had a higher influence on landslide occurrence than distances over 300 m. However, the highest LSI was in the 200–300 m class rather than a distance closer to a road. This causal factor contributes to landslides when the slope is eroded from the bottom of the slope [83]. Various studies have shown different patterns depending on the interval of the class. Some studies have shown higher susceptibility closer to roads [81,84,85], while others presented similar findings to the current study [78,86,87]. The results for the occurrence of landslide and distance to drainage suggest that the 0–100 m class is the main contributor. The influencing factor here could be the degree of saturation of material in the area and streams eroding the slopes [25]. The study area has three main fault lines in the south and one fault line in the north. Under the distance to faults factor, the 0–500 m class has the most impact on the probability of landslides, while other classes contribute comparatively less. A landslide occurs closer to a fault because of fractures in the rock masses [36,88], as concluded in other studies [69,89]. The influence of lithology on landslide occurrence is highest where the class consists of phyllite, slate, and schist, with the highest value of LSI [40]. From the results, the bare area land cover class has the highest influence on the occurrence of landslides in the classes considered. Bare areas lacking vegetation are more prone to landslides because vegetation helps prevent erosion through an anchoring effect [90].

The AUC of 0.8 and higher suggests that all three models performed well and produced reliable LSMs. The IV model had the highest AUC, while the WOE and LR models had slightly lower values. An LSM developed for the Zigui–Badong area near the Three Gorges Reservoir in China found that the IV model performs better than the LR model, although the deviation was small [91]. Ozdemir and Altural [40] have also reported that the WOE model has a higher AUC than the LR model. By contrast, a study conducted in the city of Mizunami, Japan found that the LR model has a better predictive performance than the WOE model [84]. The correlation analysis indicated that the LSM generated from the WOE and LR models had similar distribution patterns, whereas the LSM by the IV method was less similar to other methods (Table 5). When landslide pixels were overlaid onto the LSM, it showed that the LSM produced with the IV model was better, since it predicted a higher percentage of landslides of higher susceptibility in both the training and control groups (Figure 9). Figure 10 shows some of the existing landslide photographs overlaid onto the LSM produced with the IV method, with most of the landslides lying in areas with high landslide susceptibility.

7. Conclusions

This paper presents a comparative performance of the IV, WOE, and LR models in developing landslide susceptibility maps along a road corridor in Bhutan. The LSM assessment on the studied road corridor showed that the probability of landslide occurrence is greater in the southern area, with the mid-region being more stable. The influence of each class of causal factor is evident from the IV and WOE models. The highest contributors to landslide occurrence from each factor were slope angle of >50, southern exposure aspect, 300–600 m elevation, 200–300 m distance to road, 0–100 m distance to drainage, 0–500 m distance to fault lines, lithology group Pzph (Phuentsholing formation), and bare area in land cover. As a whole, the LR model has greater advantage in identifying the influence of causal factors on the probability of a landslide. In order from most to least influential, the factors were: land cover, aspect, lithology, slope angle, distance to fault lines, elevation, distance to road, and distance to drainage.

The LSMs generated by the three methods were evaluated for suitability by deriving the AUC and overlaying the actual landslide map on the LSM. The correlation coefficients for all the models were greater than 0.699, indicating strong correlation. The WOE and LR models were determined to be similar with a correlation coefficient of 0.85. The validation showed the success rate and prediction rate of all the models to be suitable. The comparative test of the models showed that the IV model was better than the WOE and LR models. We therefore recommend the IV as the most suitable method to predict future landslides. However, its performance can be improved by using higher-resolution images and large-scale maps.

We conclude that the methods considered in the present study can be suitably performed to study areas in Bhutan to develop an accurate and reliable LSM. The LSM thus developed may be used to support planning and decision-making for landslide prevention and mitigation activities during the construction, operation, and maintenance of the road network in Bhutan.