An integrated quantitative and qualitative approach for landslide susceptibility mapping in West Sikkim district, Indian Himalaya

Abstract

Landslides rank as the third most common natural disaster globally, and the Indian Himalaya Region is no exception, experiencing severe impacts during the rainy season. This study focuses on creating a comparative landslide susceptibility map for the West Sikkim district in India using probabilistic and heuristic approaches. The frequency ratio (FR) and information value (IV) methods are employed for the probabilistic approach, while the analytic hierarchy process (AHP) is used for the heuristic approach. Eleven factors are considered in the analysis. The resulting landslide susceptibility (LS) map demonstrates accuracies of 77% for FR, 74% for IV, and 57% for AHP methods. Preliminary qualitative risk assessment is conducted, incorporating building and population density, as population and buildings are the most vulnerable elements in the society. The LS map with the highest accuracy (from FR) serves as the landslide potential factor, combined with building and population density as the risk damage potential factors for risk zonation. The resulting risk zonation map classifies the study area into high-risk (3%), medium-risk (14%), and low-risk (83%) zones. This study primarily addresses the 3% high-risk area where landslides pose a significant threat to population and infrastructure, aiming to inform policy implementation and mitigation measures.

Keywords:

• Risk assessment
• landslide susceptibility
• probabilistic approach
• heuristic approach

1. Introduction

A mountainous terrain may pose a serious threat to communities living in it owing to landslide disasters that are generated both by unsustainable anthropogenic activities such as mining, road cutting or urbanization, and by natural factors such as earthquakes and excessive rainfall (Martha 2011; Dai et al. 2022). The Himalayan region receives a lot of rainfall and accounts for more than 70% of all fatal landslides in India (Froude and Petley2018). The majority of Himalayan region including Indian Himalayan region is classified as a major risk zone due to the incidence of a variety of natural disasters, including landslides. The Sikkim Himalaya in this context experiences similar fate. Landslides are defined as the movement of debris, rock mass, or earth down a slope as a result of gravity (Cruden 1991; Yalcin et al. 2011). A landslide can inflict significant economic impact, especially in densely populated areas where dislodged rocks hinder highways and sub-roads. As a result of global and micro climate change, population growth, and ever-expanding transportation networks, the mortality toll has risen substantially in recent years (Althuwaynee and Pradhan 2017). Landslide risk has increased in recent decades as a result of population growth, which has resulted in the extension of residential areas; the construction of new settlements on steeper slopes and/or in regions with inadequate foundation (bedrock, soil) in terms of structure and stability (Guzzetti et al. 1999; Yalcin et al. 2011; Zhong et al. 2020). The expansion of infrastructure as a result of population increase in the mountainous terrains of the Himalayas is the greatest source of concern for the incidence of landslides (Jaiswal et al. 2011). Landslides are often discrete events that are not particularly enormous in scale, they can occur with a high frequency in a given region and do significant damage to both human and physical infrastructure (Van Westen et al. 2006). This is due to a higher susceptibility of surface soil to instability, as well as a higher vulnerability of the exposed population (Nadim et al. 2006; Hoyois et al. 2007). As a result of urbanization, more hill slopes are excavated ‘indiscriminately’, leaving them exposed (Jaiswal 2011). Excavated slopes along highways can cause landslides, affecting not only the alignments, automobiles, and people, but also the economy (Van Westen et al. 2006). Landslides are considered to be the third most catastrophic type of natural disaster in the world (Zillman1999; Castellanos 2008). In the Indian context too, landslides are ranked third in terms of devastation caused by natural disasters, trailing only floods and hurricanes (Atlas Magazine 2020). According to Kirschbaum et al. (2010), developing countries account for roughly 80% of all landslide fatalities reported worldwide in the recent decade, with India accounting for nearly 8% of all landslide fatalities (Ghosh 2011). To comprehend a landslide disaster and mitigate its influence on society, it is vital to make a realistic assessment of the projected disaster losses. The evaluation entails quantifying both direct and indirect risks quantitatively, which requires a large number of input factors to measure each component of hazard, vulnerability, and element at risk (Fell et al. 2008; Jaiswal 2011). Most research focus on landslide identification, dynamic modelling, and susceptibility analysis, while few focus on landslide hazard and risk quantification.As a result, defining Landslide susceptibility (LS) is the first step in assessing landslide risk (Kouli et al. 2010; Rozos et al. 2011; Kayastha et al. 2013; Foumelis et al. 2018).Landslide risk assessment, which includes landslide hazard analysis, vulnerability analysis of disaster-bearing bodies, and risk zonation, has steadily become one of the most important subjects in the field of landslide study (Sui et al. 2020).The LS term is frequently used to refer to the location of probable landslides in a certain area, determined by a set of topographical criteria. Future landslides could be induced by the same conditioning parameters that triggered current and previous landslides, according to various research conducted in the past (Carrara et al. 1999). Because qualitative approaches rely on expert opinion, they are very subjective, with the weight of the elements that lead to vulnerability being assessed based on expert knowledge and experience (Clerici et al. 2010; Chalkias et al. 2014; Tsangaratos et al. 2017). Semi-quantitative techniques are qualitative methodologies that work in weighting and rating procedures (Ayalew et al. 2005; Achour et al. 2017).As of now, only a little amount of research has been done to determine the vulnerability and risk of buildings and infrastructures that are vulnerable to landslide disasters (Singh et al. 2021).

Considering data availability, codal provisions, and the geographical conditions of the Indian subcontinent, this work aims to establish a preliminary qualitative risk assessmentwith respect to population and building densityfor hilly terrains of the study area. For the purposes of this definition, risk is a measure of the possibility and severity of a devastating impact on one’s health, property, or environment (AGS (Australian Geomechanics Society Sub-Committee on Landslide Risk Management) 2000). Varnes (1984) defined landslide risk as the likelihood of lives being lost, people being injured, property being destroyed, and economic activity being hampered as a result of a specific landslide event occurring in a specific area and reference time.Landslide risk assessment requires several parameters, making quantification challenging (Van Westen et al. 2006; Li et al. 2014). Many research has been done on risk zonation, however the methods used to quantify landslide risk differ according on the scale, mechanism, and data availability. Also, due to the considerable variation in terrain, geology, climate, etc., no common framework for assessing landslide risk exists (Singh et al. 2017). The spatial examination of susceptibility to landslide hazards and the assessment of exposure to this particular hazard make up landslide risk assessment (Mishra and Sarkar 2020; Nor Diana et al. 2021). Remote sensing and GIS have recently given a new dimension to environmental hazard and risk mapping. This technique can manipulate, acquire, store, retrieve, and model large Spatio-temporal data (Panchal et al. 2015), and can be used to augment other susceptibility mapping methods and models (Akbari et al. 2014). The literature reveals that different strategies were used depending on the objectives, scale, and location of the problem. All of these models have flaws, however, there is no broad consensus on whether method or model is best for assessing landslide risk (Mishra and Sarkar 2020). Apart from illicit hill cutting and mining, increased landslide activity in the Himalayan region has been linked to the construction of road networks (Singh et al. 2014; Froude and Petley2018; Hosenuzzaman et al. 2022).According to the research conducted by Dikshit et al. (2020), about 15% landslides were occurred in the Indian Himalayan Region from 2004-2017 with respect to the global context. The Himalayan landslides are basically triggered by seimic activity, rainfall and growing anthropogenic activities. Lots of studies in the landslide study in the Himalayan region deals in the susceptibility mapping but there are very few or no studies regarding the early warning and risk assesement especially in the eastern Himayalan region (Dikshit et al. 2020). There has also been a lot of studies in climate change and global warming but the exact link is yet to be established in relation to the landslides especially in the context of Indian Himalaya (Dikshit et al. 2020). Some of the past landslides in the Indian Himalaya and their loss are shown in the Table 1.

Table 1. Past landslides in the Indian Himalaya and their losses.

Sl No.	Place	Year	Effect
1	Chamoli landslide	1868	70 people died
2	Nainital landslide	1880	150 people died
3	Chamoli landslide	1970	55 people including many livestocks
4	Assam landslide	1991	300 people died
5	Itanagar landslide	1993	25 people buried alive
6	Kalimpong landslide	1993	40 people died
7	Nagaland landslide	1993	500 people died
8	Luggarbhati landslide	1995	65 people died
9	Deorali, Gangtok landslide	1995	32 people died
10	Okhimath and Malpa slides	1998	380 people died
11	Gyalshing Bazar slide	2000	5 people died
12	Manzing landslide	2005	7 people died
13	Dharla landslide	2007	62 people died
14	Pithoragarh landslide	2009	43 people died
15	Sikkim landslide	2011	More than 111 people died
16	Chungthang landslide	2012	27 people died
17	Kedarnath disaster	2013	More than 5700 casualties
18	Chamoli disaster	2021	More than 200 deaths

The main goal of this study is to quantify landslide hazard probability and preliminary qualitative risk zonation in the study area with limited available data. Despite these restrictions, an attempt has been made in this study for risk zonation in West Sikkim district of India. The first section discusses the use of statistical and knowledge-driventechniques to evaluate landslide hazard using data from a landslide inventory. To quantify risk, detailed information on the vulnerability and number of elements at risk is required, as well as the chance of a specific hazard scenario occurring in a given area over a given return period (Vranken et al. 2015). Lack of temporal data on landslide incidence is one reason why landslide risk assessment and mapping remain challenging (Pellicani et al. 2014). Current study tries to analyze the qualitative risk zonation using the landslide potential factors and the risk damage potential.

1.1. Study area

The study area falls in the north eastern part of India in the Lesser and Greater Himalaya of the Eastern Himalaya covering West Sikkim district of Sikkim (Figure 1). Thestudy area West Sikkim which covers an area of 1162 sq. km is bounded by 27.614° N − 27.112° N (latitude) and 88.017° E − 88.353°E (longitude), it ranges56km from North to South 35 km East to West. The elevation in the study area ranges from 269 m to 7047 m above mean sea level and slope ranges from 0 to 86°.

Figure 1. Map of the study area.

Forested ecosystems account for more than 50% of totalstudy area. Sikkim, which makes up 0.20% of India’s total land area, is regarded one of the ecological hotspots to the east of the Himalayas (Govt. of Sikkim, 2022).The lush green vegetation in the area is because of the prolonged and intense rainfall that occurs in both the winter and summer seasons. The mean daily temperature ranges from 4° C (January) to 17° C (July). Kangchenjunga range, which holds the world’s third-highest mountain peak covers northern part of the West Sikkim district (Sharma et al. 2023).Because the structural slope runs from north to south, all of the state’s major drainage flow in a southerly direction. Rangit river is the major river in the study area (Hazra and Krishna 2019). The interaction between denudational and depositional processes responsible for generating the land surfaces of the Sikkim-Himalayan area can be understood in general terms. Landslides and other mass movements are common in the state. This is because of the fragile lithology, rainfall, ancient debris flow, palaeo-landslides, and complex sloping conditions which render the situation sensitive to slope instability (Rawat 2015). During the monsoon season, the study region receives 300-900 mm of rainfall each month, which causes the most landslides. The relative humidity still remains very high during the monsoon seasons, providing adequate moisture for the rock and soil, particularly with clay-rich soils and rocks that allow swelling, causing unpredictable landslides in places. Increasing landslide/soil erosion problems in various locations under different lithology/physiographic zones, relief, angles of slope, and altitudes may pose a threat to rare and endangered species. Sometimes thick piles of debris block the river’s flow paths and when these breaches of natural blocks create flash floods in the downstream regions (Rawat 2015). Many small and large landslides studied in West Sikkim district during the present study and attempts were made to assess the likelihood of susceptibility of these slides. Further these susceptibility factors coupled with risk damage potential are used for the preparation of risk zonation map. Some of the significant landslides surveyed during the field visits in the study area are shown in the figure (Figure 2).

Figure 2. Some major landslides in the study area surveyed during the field visit. A) Som landslide; B) Bojeck landslide; C) Pipaley landslide; D) Kamling landslide; E) Ridang landslide; F) Gerethang landslide.

1.2. Geology of study area

The tectonic heterogeneity, polyphase metamorphism, and unfossiliferous nature of the litho units have led to scholarly debate, appraisal, and re-evaluation (Basu 2013). The geological formations of the state are arranged in reverse tectonic order from the south to the north of the state (Auden 1935). The study area falls in the Lesser and Greater Himalayan division. Low to medium grade Precambrium metapelites and the Baxa carbonate-quartzite association (Daling Group) dominate the landscape of Sikkim Himalaya, which is also covered by high grade gneiss (Darjeeling Group) and deformed granite gneiss (Lingtse Granite Gneiss) (Bose 1891, Acharyya 1989).Gondwana-equivalent Rishi Group (Namchi formation, Rangit Pebble Slate) and Tethyan rocks are among the Palaeozoic Mesozoic rocks (Everest Pelite, Everest limestone, Lachi, and ChhoLahmo Formations) (Ray 1976).The rocks can be seen highly deformed due to the different stages of orogeny (SSDMA 2022). Microfractures, fissures, joints, foliation surfaces, faults (shear), and several generations of folds are all examples of deformation. Each of these structures, which range in age from Proterozoic to Mesozoic, has its own set of tectono-sedimentary, magmatic, metamorphic, and structural features. The principal lithotypes of the Sikkim Darjeeling Himalaya, which have complicated fold-thrust tectonostratigraphic sequences that generate nappes, windows, and klippen, among other structures. The main central thrust (MCT) is the major geological fault and the most critical tectonic zone passes through the study area. The Tista Lineament, the biggest lineament in the district passes almost diagonally through east to west of the district (Figure 3). The planar and the non-planar discontinuities are the chief source of instability in the area.

Figure 3. Geological map of the study area.

2. Materials and methods

The first objective of the current work is to prepare a landslide susceptibility map (LSM) based on the comparative study of different methods for the West Sikkim district of the Sikkim Himalaya. As we have already known by now that there is no fixed guideline or rule for preparing a landslide susceptibility map of any study area, the selection of causative factors for landslides in statistical and data-driven methods generally depends on the nature of the study area and the availability of the data (Ayalew and Yamagishi 2005; Magliulo et al. 2008; Yalcin 2008). The different methods used here are frequency ratio method (FR), information value method (IV) and analytic hierarchy process (AHP). The first two methods are statistical driven probabilistic methods while AHP is the knowledge driven heuristic approach. Probabilistic methods are used to calculate the weights of various conditioning factors and the weights are applied to generate the landslide susceptibility map. On the other hand, the heuristic method used the weights based on the expert judgement.The comparative method helps us to choose the best fit method. As there is no universal guideline for the selection of the landslide conditioning factors, therefore depending on the terrain of the study area and the data availability, eleven factors are considered for the preparation of the landslide susceptibility map. These factors are slope aspect, slope angle, slope curvature, drainage density, elevation, lithology, land use and land cover (LULC), normalized difference vegetation index (NDVI), geomorphology, lineament density, and soil type. The principal objective of the study is preparing a landslide risk zonation map for the study area, which can be used for planning for the policy makers and the administrators. Because landslides wreak havoc on settlement areas, landslide risk assessment must be factored into urban planning strategies.The best fit susceptibility map was used to calculate the interrelationship between each landslide causative factor and landslide distribution in order to assess landslides. The detailed methodology adopted in the study is given below (Figure 4). The data sets and the software used for the study are shown in Tables 2 and 3.

Figure 4. Methodology flowchart for the present study.

Table 2. Data sets used for the present study.

Data type	Source	Year	Spatial Resolution	Output extracted
Satellite imagery	Resourcesat − 1 (LISS IV - Mx), Sentinel 2	2019	5.8 m; 10 m respectively	landslide map, slope aspect, slope angle, slope curvature, drainage density, elevation, land use and land cover (LULC), normalized difference vegetation index (NDVI), geomorphology, building density map
	ALOS PALSAR	2009	10 m
Google Earth image	Google Earth Application	2019	----	landslide map, geomorphology, LULC
Geological map	Published by GSI	2014	----	lithology
Field survey		2017-2021	----	field survey for landslide locations, the ground truth of LULC
SRM data	NBSS&LUP, India	2000	----	soil type
Structural data	Bhuvan portal, NRSC, India	2019	----	lineament density map
Population data	The Gridded Population of World Version 4	2022	----	Population density map

Table 3. Software used for the present study.

Softwares	Variables Extracted
ArcGIS	Integration of thematic layers
Erdas Imagine	Preparation of thematic layers
MS Excel	Plotting graphs

2.1. Landslide mapping

Identification of landslide incidents was made possible through the use of satellite image (Resourcesat, LISS 4) visual interpretation in conjunction with field investigations and previous records from local administration. The satellite image has a spatial resolution of 5 m after correction. The processes of detection, identification, description, and evaluation of the significance of an item and pattern imaged, are referred to as interpretation. The basic elements that are considered while doing the visual interpretation are shape, size, pattern, tone, texture, shadows, location, association, etc. The locations of these landslides, as well as their coordinates and borders, were marked in the map in the form of polygons.In this study, a total of 685 landslides were mapped and these were further divided in the 70:30 ratio using the random sampling in the GIS tool, as most of the researchers use this proportion for their study (Martha et al. 2013). 70% landslides were taken as the training data set, and the remaining 30% landslides were used for the testing of the model. The landslide in the study area has a maximum area of 0.33 square kilometres and a minimum area of 0.000123 square kilometres, with a total area of about 3 square kilometres. The mapped landslides were also validated in during the field survey and in Google Earth image.

2.2. Frequency ratio method

For landslide susceptibility analysis, the frequency ratio model is a bivariate statistical analysis that examines landslide events and associated factors (Lee and Pradhan 2006). The frequency ratio model is based on the relationship between landslide distribution and each of the factors that affect landslides (Mandal et al. 2018). The frequency ratio value is calculated by dividing the percentage of landslide pixels by the percentage of entire area pixels in each factor’s class (Pradhan and Youssef 2010). A value greater than 1 indicates a higher correlation, while a value less than 1 indicates a lower correlation with respect to the occurrence of landslides. Equation (1)(1) $$FR=\frac{\left(\frac{A}{B}\right)}{\left(\frac{\mathrm{C}}{\mathrm{D}}\right)}=\mathrm{E}/\mathrm{F}$$(1) was used to calculate the frequency ratio (Arabameri et al. 2019). (1) $$FR=\frac{\left(\frac{A}{B}\right)}{\left(\frac{\mathrm{C}}{\mathrm{D}}\right)}=\mathrm{E}/\mathrm{F}$$(1)

Where A is the number of landslide pixels in each factor class, B is the total number of landslide pixels in the study area, C is the total number of pixels in each factor class, D is the total number of pixels in the study area, E is the percentage of landslide events in each factor class of the variables, and F is the relative percentage of the area of the total area for each subclass.

The FR values of the conditioning factors are consolidated together in a GIS environment to produce a landslide susceptibility index (equation 2(2) $$\mathit{LSI}=\sum \left(FR\right)i\mathrm{ }(\mathrm{i}=\mathrm{ }1,\mathrm{ }2,\mathrm{ }\dots \dots \mathrm{n})$$(2) ). (2) $$\mathit{LSI}=\sum \left(FR\right)i\mathrm{ }(\mathrm{i}=\mathrm{ }1,\mathrm{ }2,\mathrm{ }\dots \dots \mathrm{n})$$(2)

Here, i = weights of landslide conditioning factors.

All thematic maps are reclassified according to the FR values for each factor class and then integrated into GIS to produce the landslide susceptibility index map (Sharma and Mahajan 2018).

2.3. Information value model (IVM)

The risk of landslide hazards of any area can be calculated based on the contribution and the relationship between causative factors (Sarkar et al. 2013). The information value model is a statistical method for spatial prediction of an event based on the relationship of the parameter and the event. This is a bivariate statistical approach that is based on the frequency of the landslide event on Landslide Inducing Factor spatial distribution over the study region (Cao et al. 2016; Banerjee et al. 2018). It is an indirect statistical approach and could be used to assess the spatial relationship between the probability of the occurrence of landslides and the classes of the conditioning factor (Du et al. 2019). In this model, the probability of occurrence of landslides is influenced by factor information value (Luo et al. 2019). The summation of the landslide pixel ratios per factor class (Si/Ni) to the total landslide pixel per factor (S/N) gives statistics for the landslide susceptibility index (Sharma 2018). In this process, the spatial prediction of landslides is based on the relationship between the landslide causative factors and landslides. The performance of this model depends on the variables of landslide conditioning being crossed with the landslide distribution (Singh and Kumar 2018). To evaluate the weights of the predictors, it is appropriate to measure the ratio of the density of landslide occurrence in each class of a causative factor to the density of landslides in total area (Chen et al. 2014). The equation (3)(3) $$=\ln \frac{\mathit{Npix}\left(Si\right)/\mathit{Npix}\left(Ni\right)}{\sum \mathit{Npix}\left(Si\right)/\sum \mathit{Npix}\left(Ni\right)}$$(3) for the calculation of weight is given below: $$\mathrm{Wt}=\ln \frac{\text{landslide density within a factor class}}{\text{landslide density within the study area}}$$ (3) $$=\ln \frac{\mathit{Npix}\left(Si\right)/\mathit{Npix}\left(Ni\right)}{\sum \mathit{Npix}\left(Si\right)/\sum \mathit{Npix}\left(Ni\right)}$$(3) where, Wt = Weight of a factor class; ln = natural logarithm; Npix(Si) = Number of pixel of landslide within class i; Npix(Ni) = Number of pixel of class i; ΣNpix(Si) = Number of pixel of landslide within the whole study area; ΣNpix(Ni) = Number of pixel of the whole study area.

The natural logarithm is used to give negative weights when the density of the landslide is lower than normal, and positive weights when it is higher than normal. Positive weight (Wt) indicates the direct correlation between predictor factor and the occurrence of landslides and also the presence of predictors in landslides. The weight value for each class of causative factors was calculated using the aforementioned equation. To prepare a map of landslide susceptibility, the landslide susceptibility index (LSI) was calculated for each pixel by summing the weighted values of each factor class, as shown below. (4) $$\mathit{LSI}=\sum \left(\mathit{IVM}\right)i\mathrm{ }(\mathrm{i}=\mathrm{ }1,\mathrm{ }2,\mathrm{ }\dots \dots \mathrm{n})$$(4)

Here, i = weights of landslide conditioning factors.

2.4. Analytic hierarchy process

The analytical hierarchy process (AHP) is a decision-making tool that assists in the resolution of complex problems using simple criteria (Panchal and Shrivastava 2021). The problem is divided into hierarchical criteria in AHP. These criteria are weighed against one another. The elements’ weights are assigned based on the expert’s assessment, and the weight values range from 1 to 9 (Saaty 1980; Saaty and Vargas 2001, Saaty 2005; Sharma and Mahajan 2018). In this study already published work is taken as the expert assessment (Mirdda et al. 2020). Weights are assigned to landslide factors and different groups within each factor in the landslide analysis. A rating of 9 indicates exceptional success from one factor to the next, whereas a rating of 1 indicates the same value for both. Typically, the decision-maker’s preference defines the distribution of values for similar trigger variables (Mirdda et al. 2020). In this study, ratings based on expert assessment and field observation are both subjective and objective. The Saaty scale determines the priority of each class (Table 4). The rankings are calculated using the eigen vector method, and the consistency of the solution is checked using the consistency ratio (Saaty1987). The equation (Equation 5(5) $$CI=\frac{\left(\mathrm{\lambda }\max \text{-}\mathrm{n}\right)}{(\mathrm{n}-1)}$$(5) & Equation 6(6) $$CR=\frac{\mathrm{CI}}{\mathrm{RI}}$$(6) ) below can be used to check the consistency of the weights for relative importance assigned during the pairwise comparison. (5) $$CI=\frac{\left(\mathrm{\lambda }\max \text{-}\mathrm{n}\right)}{(\mathrm{n}-1)}$$(5) where λ max = largest eigenvalue of the matrix, n = order of the matrix (6) $$CR=\frac{\mathrm{CI}}{\mathrm{RI}}$$(6)

Table 4. Scale for pairwise comparison (Saaty 1995).

1	Equal importance
3	moderate importance of one over another
5	strong or essential importance
7	very strong or demonstrated importance
9	extreme importance
2,4,6,8	intermediate values
Reciprocals	reciprocals for inverse comparison
1.1-1.9	when the elements being compared are nearly equal

where RI = random index (average of consistency index

The random index (RI) varies with the matrices order. Inconsistency is accepted up to 10% (CR 0.1). If the CR is greater than 10%, the factor rating and subjective value judgement must be revised. The revision may involve reassessing factors and their classes. Table 5 shows Saaty’s random index value. From equation 7(7) $$\mathit{LSI}={\sum }_{i=0}^n\left(\text{cw x fw}\right)$$(7) , landslide susceptibility index (LSI) is measured. (7) $$\mathit{LSI}={\sum }_{i=0}^n\left(\text{cw x fw}\right)$$(7) where fw = weight of factor, cw = weight of classes within each factor.

Table 5. Random index values based on matrix size (Saaty 1980).

n	1	2	3	4	5	6	7	8	9	10	11
RI	0.0	0.0	0.52	0.89	1.11	1.25	1.35	1.40	1.45	1.49	1.51

2.5. Validation of the model

Receiver Operating Characteristic (ROC) Curve: The ROC curve shows how the probability of a true positive rate varies when compared to the probability of a false positive rate. This is a type of cross-validation used to ensure the work’s authenticity and validity. The two distinct data sets are used to test the model’s performance in the sample (Martha 2011). The training data (70% sample) were used for the success rate curve and the testing data (30% sample) is used for the prediction rate curve. The success rate curve obtained by plotting the cumulative landslide area against the cumulative total study area. The prediction rate curve is the relationship between the cumulative percentage region of non-modelled landslides and the cumulative percentage region of the study area. These curves provide coherence in landslide susceptibility evaluation (Sarkar et al. 2013).

2.6. Landslide conditioning factors

2.6.1. Slope aspect

The direction of the maximum slope with respect to magnetic north is indicated by the slope aspect. The aspect of the slope has an impact on erosion, evapotranspiration, desertification, solar heating, and soil water retention (Khan et al. 2019). Most researchers use this slope dimension to map landslide vulnerability. The slope map is extracted from ALOS PALSAR DEM of 12.5 m resolution which is further resampled to 5 m resolution to maintain the uniformity with the other layers. The slope aspect is divided into ten classes (Figure 5A).

Figure 5. Landslide conditioning factors A) aspect, B) slope, C) curvature, D) drainage density, E) lithology, F) LULC, G) NDVI, H) elevation, I) lineament density, H) geomorphology and K) soil type.

2.6.2. Slope gradient

The rate of change in elevation with respect to a horizontal surface is known as slope gradient. Slopes fail when downward shear stresses exceed the shear strength of soil or rock mass. Thus, slope stability is directly affected by slope angle (Regmi et al. 2014). The probability of the slope failure increases with increasing slope gradient (Panchal and Shrivastava 2021). The slope map is extracted from ALOS PALSAR DEM of 12.5 m resolution which is further resampled to 5 m resolution. The slope gradient map is divided into 5 class (Figure 5B).

2.6.3. Curvature

Slope curvature is a widely used conditioning factor for the landslide susceptibility mapping (Pourghasemi et al. 2013). Curvature is another DEM-based derivative that may impact slope stress and groundwater distribution. The slope curvature regulates water inflow and outflow and plays an important role in slope failures (Ramesh and Anbazhagan 2015). Curvature is prepared from DEM and is divided into three classes, viz. concave, flat and convex (Figure 5C).

2.6.4. Drainage density

The drainage density is defined as the total stream length of all orders in stream eroded topography (Demir et al. 2013). The inclusion of drainage density is required because it causes mass wasting and slope failure. In the study area, the drainage density is divided into five classes (Figure 5D).

2.6.5. Lithology

Because lithological units have varying landslide susceptibility values, they are vital to be considered. The resistance of rocks to weathering and erosion is dependent on their lithology, it is an important factor when assessing landslide susceptibility (Singh and Kumar2018). The map is used at a scale of 1:25000. There are six classes of lithology based on published work and field survey (Rawat 2015; GSI 2016; Biswakarma et al. 2020). These are Permafrost Area, Central Crystalline, Everest Pelite, Tourmaline Granite, Lingtse Gneiss and Daling Group (Figure 5E).

2.6.6. Land use and land cover (LULC)

LULC is one of important landslide conditioning factors used by almost all the researchers in landslide susceptibility mapping. Landslides are favoured by changes in land-use patterns, such as agricultural expansion and construction, as well as deforestation, which disrupts and changes the environment (Jaiswal et al. 2010). Resourcesat-2A LISS-4 multispectral (2019) bands 2, 3, and 4 were used along with ground-truth testing to plan land use/land cover classification. For land use/land cover class in the study area, supervised classification techniques were used to make interpretations. 50 training samples for each class were used for the supervised classification. A hybrid classification method (both supervised and unsupervised) was further used to address the cumulative errors. To categorise the land use/land cover class, a Google Earth map and an NRSC-prepared LULC map were used as a reference. Seven distinct land use/land cover classes were identified in the study area. Snow/Glaciers, Barren Land/Rock, Waterbodies, Fallow Land, Built-up, Agriculture/Plantation, and Forest Cover are the mapped classes (Figure 5F).

2.6.7. Normalized difference vegetation index (NDVI)

The NDVI is calculated using the equation 8(8) $$\text{NDVI}=\left(\text{NIR}\text{-}\mathrm{R}\right)/\left(\text{NIR}+\mathrm{R}\right)$$(8) . (8) $$\text{NDVI}=\left(\text{NIR}\text{-}\mathrm{R}\right)/\left(\text{NIR}+\mathrm{R}\right)$$(8) where, NIR is the near infrared band and R is the red band of the electromagnetic spectrum (Fiorucci et al. 2019). The value of NDVI ranges between −1 to +1, where −1 indicates portion with no vegetation and +1 with the dense vegetation. ERDAS IMAGINE^® software was used for this stage, which generated a matrix of real (floating point) numbers in the range between −0.73 and + 1.0 for the image (Figure 5G). The vegetation index uses variations in plant, soil, and water reflective properties that reflect different wavelengths of visible and non-visible light to carry out the calculation. This study utilized this variation in reflective properties to effectively differentiate between green plants and non-plant objects (Ni et al. 2019).

2.6.8. Elevation

Land elevation is one of the parameters which affect landslide risk and the researchers use it frequently for landslide susceptibility mapping (Yilmaz2010; Park 2011; Yalcin et al. 2011; Du et al. 2017; Kose and Turk 2019). The elevation is related to other variables, such as rainfall, types of soil, types of vegetation, and coverage of vegetation, which affect stability on the slope. In this study, the elevation is divided into ten classes, with 269 m as the lowest elevation and 7047 m as the highest elevation (Figure 5H). Elevation is prepared from ALOS PALSAR DEM and is divided into ten classes.

2.6.9. Lineament density

Lineament features in hilly terrain indicate joints, fractures, faults, and other structural features of weakness and thus directly relate to the occurrence of landslides. Geological structures such as fractures, faults, and cleavages of rock have a significant role to play in building up pore pressure (Ramakrishnan et al. 2013).There are a number of active landslides in the fault zones and major lineaments in the area. The lineament map was prepared using the Bhuvan portal of National Remote Sensing Centre, India (https://bhuvanapp1.nrsc.gov.in/thematic/thematic/index.php) and Resourcesat-2A LISS 4 satellite data (Figure 5I). For the present study, ten classes are prepared in lineament density.

2.6.10. Geomorphology

To create a geomorphic map, visual analysis of Resourcesat-2A LISS 4 (2019) data was conducted at a scale of 1:25,000. The reference map for preparing the geomorphology map was adopted from published work (Anusuya 2017). The area of study is divided into ten geomorphological units such as, Very High Mountain; Narrow Valley; Ridge; Perpetual Snow/Glacier; Escarpment; Low Mountain; High Mountain; Rocky Cliff; Less High Mountain and Mid Mountain (Figure 5J).

2.6.11. Soil texture

Both topography and rock types govern the soil texture. The soils, with minor erosion in the fluvial valleys, are moderately thick, well-drained, fine loamy soils with the loamy surface. The soils on the cliffside are very soft, extremely drained, while the loamy surface soils on the tops and ridges are relatively shallow, overly drained, coarsely loamy (Rawat 2015).Five different classes of soil are demarcated in the study area, which is Fine loamy; Coarse loamy; Loamy skeletal; Fine silty and Snow-covered (Figure 5K). The soil data has been prepared from the Soil Resource map obtained from NBSS & LUP, India (2000).

2.7. Preliminary qualitative risk assessment w.r.t population and building

For as long as there have been slopes and landslides, there has also been a need for some kind of risk assessment and management. The more conventional adaptations of risk assessment and management concepts, carried out in a qualitative manner, have been put into effect for landslide risk zoning for the purpose of urban planning and highway slope management for several decades (Hungr et al. 2005). The capacity of resources to sustain damage as a result of the landslide potential has been taken into consideration as the landslide risk, and the mapping of these risk potentials has been assigned the title of landslide risk map (Banshtu et al. 2020). In order to perform an initial assessment of the landslide hazard zonation in the study region, the landslide susceptibility map is essential. In order to determine the landslide hazard zone, it is necessary to identify the places that are susceptible to landslides and to perform a calculation that estimates the likelihood of a landslide occurring within a specific time frame. Due to the fact that it is usually difficult to demarcate this time span, the landslide hazard is frequently represented by landslide susceptibility (Brabb1984; Banshtu et al. 2020). Landslide risk is calculated with the following equation 9(9) $$\mathit{NDBI}=\mathrm{ }\frac{\mathit{SWIR}-\mathit{NIR}}{\mathit{SWIR}+\mathit{NIR}}$$(9) as per Hungr et al. 2005: $$\mathit{Landslide}\mathit{Risk}=\mathit{Hazard}x\mathit{Vulnerability}x\mathit{Exposure}$$

Where hazard is the spatial and temporal probability, vulnerability is the degree of loss and exposure is the elements at risk like population, buildings, infrastructures, etc. However, for the current study the different version of the landslide risk is used as per Banshtu et al. 2020, which is given as follows. This is done as per the available data to assess the preliminary qualitative risk w.r.t. to population and building $$\mathit{Preliminary}\mathit{Landslide}\mathit{Risk}=\mathit{Landslide}\mathit{potential}x\mathit{Risk}\mathit{Damage}\mathit{potential}$$

Where landslide potential is represented by the landslide susceptibility map prepared using eleven conditioning factors mentioned previously, whereas the risk damage potential is represented by population density map and the building density map. The current study is more of a qualitative study rather than the quantitative study; therefore, the quantitative estimation of risk is not mentioned.

2.8. Population density

The population is the most exposed and affected element in any kind of disaster. Therefore, to assess the population density and the preparation of population density map is essential for the disaster risk study. The population data was collected using the census data from the Government of Sikkim (2022). The population of the West Sikkim district according to the 2011 census is 1,36,435. The data was also cross validated with the population density data of “The Gridded Population of World Version 4” (NASA). The data collected were further plotted in the GIS platform to produce the population density map (Figure 6A).

Figure 6. Risk damage potential factors A) population density map B) building density map.

2.9. Building density

The building data was prepared through the indexing method using normalized difference building index (NDBI). The data was further validated using the building data collected from the Bhukosh, GSI (Govt. of India). The satellite data which was used to analyse the NDBI is the Sentinel 2 data which has the spectral resolution of 10 m further enhanced to 5 m after orthorectification and resampling. NDBI uses the NIR and SWIR bands of the electromagnetic spectrum to differentiate the built-up structure from the other structures. (9) $$\mathit{NDBI}=\mathrm{ }\frac{\mathit{SWIR}-\mathit{NIR}}{\mathit{SWIR}+\mathit{NIR}}$$(9)

The ratio generally rectifies the terrain illumination differences and the atmospheric effects to get the desired values. The values of NDBI generally ranges from −1 to 1. The negative values generally indicate the non-built class whereas the values close to 1 are confirmative of built class. Keeping this relation in mind, the result obtained from NDBI were classified into 5 classes using the natural breaks in the GIS software. The highest values close to 1 were further extracted as the building class (Figure 6B).

3. Results

3.1. Frequency ratio

The weights of the sub-factors of the landslide conditioning factors based on the frequency ratio method are given in the Table 6. The south facing slopes are gaining the maximum FR weights in the study area. South-east, south and south-west facing slopes have obtained the weights 1.76, 1.61 and 1.15 respectively. South-east slope covers 10% among the subclass which is around 18% landslide area among the ten sub-class. On the contrast north facing slope are having the lowest FR weights, where north, north-west and north-east have FR weights of 0.24, 0.33 and 0.45 respectively. Slope gradient above 40° have obtained the higher weights. The highest value is obtained for sub-class 45°-86° followed by 40-45° sub-class and the values are 1.92 and 1.17 respectively. The sub-class 45°-86° covers an area of 14% but landslide area in the class is 28%, which shows that this class has maximum landslide probability. The lowest value for slope gradient is obtained for 15-25° sub-classes. For the curvature the highest weight of 1.01 is obtained for convex curvature, followed by concave curvature (FR = 1). The flat curvature has the lowest weight. In the drainage density the highest weight is obtained for the sub-class 0-1 sq. km followed by 1-2 sq. km and 4-8 sq. km. The weights are 1.8, 0.8 and 0.68 respectively. The lowest weight is obtained for sub-class 2-3 sq. km. 3500-4000 m elevation is showing the highest FR weight of 3.77 followed by 3000-3500 m sub-class with a weight of 3.74. The lowest weight for elevation is given by class above 4500 m. For the LULC class the highest weight is given by built-up sub-class with a weight of 5.17 followed by fallow land with a weight of 1.42. The lowest weight in this class is given by sub-class barren rock/land. NDVI class is showing highest weight for sub-class 0.2−0.4 with a weight value of 3.47. The lowest weight of 0.06 is obtained for sub-class −0.73-0. For the lithology, Everest pelite is showing the maximum value of 3.48 followed by Daling group with 1.08. Permafrost area has the lowest value. In the geomorphology section High Mountain class has the maximum value of 5.27 followed by Narrow valley and Rocky cliff with values 2.01 and 1.77 respectively. The lowest value is obtained for Perpetual snow/glacier and Very High Mountain. Coarse loamy soil type has the maximum value of 1.57 followed by Loamy skeletal with value of 0.997. Snow covered sub-class has the minimum value. For the lineament density class, the highest value is obtained for 0–150 subclass and the lowest value is obtained for 1400–1600 sub-class. The graphical representation in the Figure 7 shows the weights obtained by different conditioning factors and their classes.

Figure 7. Graphical representation of weights (FR) from different conditioning factors A) aspect, B) slope, C) curvature, D) drainage density, E) lithology, F) LULC, G) NDVI, H) elevation, I) lineament density, J) geomorphology, K) soil type.

Table 6. Weights obtained from the probabilistic methods (FR & IV).

Conditioning factors	Classes	Area%	Landslide%	FR-Wt.	IV-Wt.
Aspect	Flat (-1)	6.62	6.11	0.92	−0.08
	North (0-22.5)	10.61	4.40	0.41	−0.88
	Northeast (22.5-67.5)	8.27	3.77	0.46	−0.79
	East (67.5-112.5)	17.45	18.21	1.04	0.04
	Southeast (112.5-157.5)	10.04	17.70	1.76	0.57
	South (157.5-202.5)	18.12	29.13	1.61	0.47
	Southwest (202.5-247.5)	8.98	10.34	1.15	0.14
	West (247.5-292.5)	12.28	8.06	0.66	−0.42
	Northwest (292.5-337.5)	4.84	1.61	0.33	−1.10
	North (337.5-360)	2.80	0.67	0.24	−1.43
Slope	0–5	6.90	6.33	0.92	−0.09
	5–10	4.96	3.66	0.74	−0.30
	10–15	6.24	4.76	0.76	−0.27
	15–20	8.29	5.83	0.70	−0.35
	20–25	9.43	6.12	0.65	−0.43
	25–30	11.25	8.23	0.73	−0.31
	30–35	14.85	12.25	0.83	−0.19
	35–40	12.14	11.30	0.93	−0.07
	40–45	11.19	13.09	1.17	0.16
	45-85.6	14.75	28.43	1.93	0.66
Curvature	Concave	45.63	45.69	1.00	0.00
	Flat	8.79	8.03	0.91	−0.09
	Convex	45.58	46.28	1.02	0.02
Drainage Density	0–1	26.75	50.26	1.88	0.63
	1–2	23.93	19.30	0.81	−0.22
	2–3	20.23	11.18	0.55	−0.59
	3–4	13.93	9.01	0.65	−0.44
	4–8	15.16	10.25	0.68	−0.39
Lithology	Permafrost Area	11.66	0.00	0.00	0.00
	Central Crystalline	64.45	61.71	0.96	−0.04
	Everest Pelite	6.14	21.34	3.48	1.25
	Tourmaline Granite	1.15	0.00	0.00	0.00
	Lingtse Gneiss	3.45	2.69	0.78	−0.25
	Daling Group	13.15	14.26	1.08	0.08
LULC	Snow/Glaciers	4.71	0.25	0.05	−2.93
	Barren Land/Rock	15.94	0.00	0.00	0.00
	Waterbodies	0.34	0.00	0.01	−4.92
	Fallow land	6.93	9.84	1.42	0.35
	Built_Up	0.54	2.78	5.17	1.64
	Agricullture/Plantation	9.00	2.22	0.25	−1.40
	Forest cover	62.55	84.90	1.36	0.31
NDVI	−0.73-0	4.64	0.26	0.06	−2.89
	0–0.2	10.64	11.97	1.12	0.12
	0.2–0.4	16.29	56.55	3.47	1.24
	0.4–0.6	29.71	28.80	0.97	−0.03
	0.6–1	38.73	2.42	0.06	−2.77
Elevation	269–1000	7.39	6.57	0.89	−0.12
	1000–1500	14.17	9.17	0.65	−0.44
	1500–2000	15.25	10.56	0.69	−0.37
	2000–2500	14.12	7.36	0.52	−0.65
	2500–3000	12.65	9.60	0.76	−0.28
	3000–3500	7.73	28.88	3.74	1.32
	3500–4000	7.09	26.73	3.77	1.33
	4000-4500	8.79	1.13	0.13	−2.05
	4500–5500	10.77	0.00	0.00	0.00
	5500–7047	2.04	0.00	0.00	0.00
Lineament Density	0–150	60.78	77.96	1.28	0.25
	150–400	4.04	2.68	0.66	−0.41
	400–600	5.51	2.41	0.44	−0.83
	600–800	6.76	4.97	0.74	−0.31
	800–1000	6.86	3.95	0.58	−0.55
	1000–1200	6.33	3.83	0.60	−0.50
	1200–1400	4.92	3.19	0.65	−0.43
	1400–1600	1.73	0.35	0.20	−1.59
	1600–1800	2.13	0.45	0.21	−1.55
	1800–2004	0.95	0.22	0.23	−1.47
Geomorphology	Very High Mountain	12.36	0.00	0.00	0.00
	Narrow Valley	2.34	4.70	2.01	0.70
	Ridge	10.66	8.17	0.77	−0.27
	Perpetual Snow/Glacier	5.36	0.00	0.00	0.00
	Escarpment	5.54	7.06	1.27	0.24
	Low Mountain	4.47	2.02	0.45	−0.80
	High Mountain	7.31	38.52	5.27	1.66
	Rocky Cliff	4.99	8.85	1.77	0.57
	Less High Mountain	19.61	17.02	0.87	−0.14
	Mid Mountain	27.37	13.66	0.50	−0.70
Soil type	Fine_loamy	22.91	10.00	0.44	−0.83
	Coarse_loamy	34.68	54.62	1.57	0.45
	Loamy_skeletal	35.14	35.04	1.00	0.00
	Fine_silty	0.37	0.34	0.91	−0.09
	Snow_covered	6.91	0.00	0.00	0.00

3.2. Information value method

The weights of the different classes of the landslide conditioning factors based on the information value method are given in the Table 6. The weights obtained from the information value method are following the similar trend as the frequency ratio method. The southing facing slopes have higher values as compared to the other classes in the slope aspect. The highest value is obtained for the south-east, south and south-west slope with values 0.57, 0.47 and 0.14 respectively. The lowest values are obtained for northerly facing slopes. The slope gradient with 40° have shown the positive value for the information value weight. This shows the positive correlation towards the occurrence of landslide. The highest valueclasses are 45–86° and 40–45° with values 0.66 and 0.16 respectively. The lowest value for the slope gradient is obtained for the 20–25° class with value of −0.43. Convex curvature is more vulnerable to landslides with highest information value weight of 0.015, followed by concave curvature with weight of 0.0011 and the lowest value is obtained for flat curvature. For the drainage density, the highest value is obtained for the class 0–1 sq. km and the lowest value is obtained for 2–3 sq. km sub-class. In the elevation factor the highest value is obtained for the sub-class 3500–4000 m followed by 3000–3500 m sub-class. The lowest value in this category is given by the sub-class 4000–4500 m. Built-up, Fallow land and Forest cover are showing the positive correlation IV weight of 1.64, 0.35 and 0.31 respectively. Waterbody class in this factor has attained a minimum value of −4.9, which indicates the minimal probability of landslide in the class. The vegetation index is showing the maximum value for the sub-class 0.2−0.4 followed by 0−0.2 with values 1.24 and 0.118 respectively. The sub-classes −0.73-0 and 0.6−1 for NDVI have shown the minimum and maximum values respectively for IV method. Lithology is an important factor for consideration while doing the landslide susceptibility analysis. Everest Pelite and Daling Group have shown more positive correlation the towards the occurrence of landslide having IV weights 1.25 and 0.081 respectively. Lingtse Gneiss is this factor class have the lowest value with weight of −0.25. High Mountain sub-class in the geomorphology factor have the highest value followed by Narrow Valley, Rocky Cliffs and Escarpment having weights 0.70, 0.57 and 0.24 respectively. The lowest value for the geomorphology is obtained for the classes Low Mountain and Mid Mountain respectively. According to IV weight the highest correlation for the occurrence of landslide in soil type class is obtained for the Coarse loamy soil with weight of 0.45. The minimum weight in this category is obtained for Fine loamy soil with a weight of −0.82. In the lineament density parameter, the only sub-class with positive correlation for the occurrence of landslides is 0-150. All other sub-classes in the lineament density are showing negative correlation with a minimum value of −1.59 for the sub-class 1400-1600. The graphical representation in the Figure 8 shows the weights obtained by different conditioning factors and their classes.

Figure 8. Graphical representation of weights (IV) from different conditioning factors A) aspect, B) slope, C) curvature, D) drainage density, E) lithology, F) LULC, G) NDVI, H) elevation, I) lineament density, J) geomorphology, K) soil type.

3.3. Analytic hierarchy process (AHP)

The weights of the classes of the landslide conditioning factors based on the AHP method are given in the Table 7. The weights of overall conditioning factors using AHP study is shown in the Table 8. The weights obtained in this method are mainly based on heuristic approach, there results obtained by this method may vary with probabilistic method. In the slope aspect the highest value is obtained for the north-west and the west sub-class followed by the southerly sub-class. The lowest value of 0.018 is obtained for the flat aspect. The highest weight in the slope gradient is obtained for the sub-class 45–86°, followed by 40–45° sub-class with AHP weight of 0.31 and 0.2 respectively. The class has the lowest weight with 0.015 weightage. In the curvature the highest weightage is obtained for convex curvature, followed by the flat curvature and concave curvature having the weightage 0.63, 0.3 and 0.06 respectively. For the drainage density the highest value is obtained for 4–8 sq. km followed by 3-4 sq. km having the weightage 0.57 and 0.22 respectively. The lowest weight for drainage density is obtained for the sub-class 0–1 sq. km with a weight of 0.03. For the elevation the highest weight is obtained for 5500–7047 m class, followed by 4500–5500 m with weights 0.28 and 0.2 respectively. The lowest weight of 0.017 is obtained for sub-class 269–1000 m. The LULC is important parameter and the highest value is obtained for built-up class followed by fallow land with weight value of 0.31 and 0.2 respectively. The lowest value in the LULC class is obtained for the Snow/Glaciers class having the weight 0.028. The highest weight in the vegetation index is obtained for 0.6−1 class, followed by 0.4−0.6 with the weightages 0.445 and 0.279 respectively. The lowest value for vegetation index is given by sub-class −0.73-0 with value of 0.035. For the lithology the highest value is obtained for the Daling Group followed by LingtseGniess and the lowest value is obtained for the Permafrost area. The weights of Daling group, LingtseGniess and Permafrost area 0.39, 0.26 and 0.024 respectively. In the geomorphology section the maximum value is obtained for the High Mountain, Narrow Valley and Rocky cliff with values 0.198, 0.155 and 0.144 respectively. The lowest value in this section is given by weights 0.017 for Perpetual Snow/Glacier sub-class. The coarse loamy soil type has the maximum impact for the landslide occurrence followed by loamy skeletal soil type. The snow-covered area in the soil type class has the lowest value. The values of coarse loamy, loamy skeletal and snow-covered area are 0.557, 0.2 and 0.03 respectively. The highest weight for the lineament density is obtained for the 1800-2004 sub-class having the weight value of 0.29. The lowest value in this section is given by 0-150 sub-class with a value of 0.015. The graphs show the different weights obtained from the AHP method (Figure 9).

Figure 9. Graphical representation of weights (AHP) from different conditioning factors A) aspect, B) slope, C) curvature, D) drainage density, E) lithology, F) LULC, G) NDVI, H) elevation, I) lineament density, J) geomorphology, K) soil type.

Table 7. Weights obtained for different class using AHP method.

Conditioning Factors	Class	1	2	3	4	5	6	7	8	9	10	Criteria Weight
Aspect	Flat (-1)	1										0.018
	North (0-22.5)	2	1									0.03
	Northeast (22.5-67.5)	3	2	1								0.031
	East (67.5-112.5)	4	2	3	1							0.042
	Southeast (112.5-157.5)	7	6	6	4	1						0.091
	South (157.5-202.5)	9	7	8	3	1	1					0.107
	Southwest (202.5-247.5)	8	6	6	4	2	1	1				0.109
	West (247.5-292.5)	7	5	5	6	3	3	2	1			0.15
	Northwest (292.5-337.5)	6	4	6	6	5	4	3	2	1		0.216
	North (337.5-360)	5	3	4	4	4	3	4	3	1	1	0.207
CI = 0.145; RI = 1.49; CR = 0.097
Slope	0–5	1										0.015
	5--10	2	1									0.019
	10--15	3	2	1								0.025
	15--20	4	3	3	1							0.037
	20--25	5	4	4	3	1						0.053
	25--30	6	5	5	4	3	1					0.074
	30--35	7	6	6	5	4	3	1				0.104
	35--40	8	7	7	6	5	5	3	1			0.158
	40--45	9	8	8	7	6	6	4	2	1		0.205
	45-85.613	9	9	8	8	8	7	6	3	3	1	0.310
CI = 0.139; RI = 1.49; CR = 0.093
Curvature	Negative	1										0.061
	Flat	7	1									0.302
	Positive	8	3	1								0.637
CI = 0.054; RI = 0.52; CR = 0.093
Drainage Density	0--1	1										0.031
	1--2	3	1									0.055
	2--3	5	4	1								0.121
	3--4	7	5	3	1							0.228
	4--8	9	8	7	4	1						0.565
CI = 0.098; RI = 1.11; CR = 0.088
Lithology	Permafrost Area	1										0.024
	Central Crystalline	4	1									0.05
	Everest Pelite	5	5	1								0.101
	Tourmaline Granite	8	3	5	1							0.18
	Lingtse Gneiss	9	7	4	2	1						0.259
	Daling Group	9	9	4	3.5	2	1					0.386
CI = 0.123; RI = 1.25; CR = 0.099
LULC	Snow/Glaciers	1										0.028
	Waterbodies	2	1									0.046
	ForestCover	2	2	1								0.060
	Agricullture/Plantation	8	7	5	1							0.167
	Fallow land	7	6	4	2	1						0.205
	BuiltUp	8	7	6	4	2	1					0.311
	Barren Land/Rock	5	2	2	2	1	1	1				0.182
CI = 0.088; RI = 1.35; CR = 0.067
NDVI	−0.73-0	1										0.035
	0−0.2	2	1									0.053
	0.2−0.4	8	6	1								0.187
	0.4−0.6	9	6	2	1							0.279
	0.6−1	6	5	5	2	1						0.445
CI = 0.107; RI = 1.11; CR = 0.096
Elevation	269-1000	1										0.017
	1000-1500	2	1									0.02
	1500-2000	3	2	1								0.025
	2000-2500	4	3	2	1							0.036
	2500-3000	5	4	3	3	1						0.052
	3000-3500	9	8	7	5	3	1					0.096
	3500-4000	9	9	8	7	4	2	1				0.126
	4000-4500	8	7	7	5	5	3	3	1			0.153
	4500-5500	6	6	8	6	6	6	3	2	1		0.201
	5500-7047	7	8	9	7	6	7	4	3	2	1	0.276
CI= 0.145; RI= 1.49; CR= 0.097
Geomorphology	Perpetual Snow/Glacier	1										0.017
	Very High Mountain	2	1									0.027
	Low Mountain	3	2	1								0.034
	Mid Mountain	4	3	3	1							0.054
	Ridge	9	8	8	6	1						0.115
	Less High Mountain	8	7	7	5	2	1					0.128
	Rocky Cliff	7	6	6	4	3	2	1				0.144
	Narrow Valley	6	3	3	2	4	3	2	1			0.155
	High Mountain	7	5	5	3	4	3	2	2	1		0.198
	Escarpment	6	4	3	3	2	1	1	1	1	1	0.129
CI= 0.142; RI= 1.49; CR= 0.095
Soil	Snow_covered	1										0.033
	Fine_loamy	6	1									0.098
	Fine_silty	4	2	1								0.11
	Loamy_skeletal	6	3	3	1							0.201
	Coarse_loamy	9	7	6	5	1						0.557
CI= 0.108; RI= 1.11; CR= 0.096
Lineament Density	0-150	1										0.015
	150-400	2	1									0.02
	400-600	3	2	1								0.026
	600-800	4	3	3	1							0.038
	800-1000	5	4	4	3	1						0.055
	1000-1200	6	5	5	4	3	1					0.077
	1200-1400	7	6	6	5	4	3	1				0.104
	1400-1600	8	7	7	6	5	5	3	1			0.161
	1600-1800	8	8	7	7	6	5	5	2	1		0.212
	1800-2004	9	8	8	7	7	6	7	3	2	1	0.29
CI= 0.14; RI= 1.49; CR= 0.094

Table 8. Weights obtained for conditioning factors using AHP method.

Conditioning Factors	1	2	3	4	5	6	7	8	9	10	11	Criteria Weight
Curvature	1											0.015
Elevation	2	1										0.019
Aspect	3	2	1									0.024
NDVI	3	3	2	1								0.033
Soil	4	4	3	2	1							0.044
Geomorphology	5	4	4	3	2	1						0.059
Lineament Density	6	5	5	4	3	2	1					0.081
Drainage Density	6	6	6	4	4	3	2	1				0.106
Slope	7	6	7	5	5	4	3	2	1			0.144
Lithology	8	7	8	6	6	5	4	3	2	1		0.192
LULC	9	8	9	7	7	6	5	4	3	3	1	0.283
CI = 0.085; RI = 1.51; CR = 0.056

These normalised weights are further used to prepare the landslide susceptibility map. The landslide susceptibility index based on the weighted value of all the three methods (FR, IV and AHP) were determined. The weight of the calculated values has been depicted in Tables 6 and 7 for FR, IV and AHP respectively. There is always a question mark while preparing the landslide susceptibility map, therefore it is recommended to go for comparative methods and choose the best method. As already discussed above, in this paper probabilistic and the heuristic methods are applied to make a comparative study for the production of landslide susceptibility mapping.

3.4. Validation of the model

Validation is required to determine the comparative assessment of the techniques used for landslide susceptibility zonation. Receiver operating characteristic(ROC) curve is used for the current study to determine the validity of the methods used. The receiver operating characteristic (ROC) curve is a useful tool for evaluating the performance of any classification scheme, and it was used to validate all the methods used in this study, including IV, FR, and AHP. As the cut-off probability changes, the ROC curve plots the probability of true positive identified landslides versus false positive identified landslides (Martha 2011). ROC is also showing a relation between the sensitivity and the specificity. The ROC curve for all the methods is shown in the Figure 10. The value of the ROC is 77%, 74% and 57% for FR, IV and AHP methods respectively. These values are well beyond the acceptable limit and hence shows that the methods applied for the current study are having good accuracy for all the methods. The result of the susceptibility map was validated using the success rate curve (SRC) and the prediction rate curve (PRC) for the methods. The 70% training samples were used to prepare the SRC and the remaing 30% testing samples were used to prepare the PRC. All the methods inferred that SRC showing good match with PRC (Figure 11). This suggests that the techniques and factors used in the study are capable of accurately predicting the spatial occurence of landslides.

Figure 10. ROC curve obtained from the different susceptibility methods.

Figure 11. Success and prediction rate curves of the different methods used.

3.5. Preliminary qualitative risk assessment w.r.t population and building

For the landslide risk the best of the three susceptibility models used in the study were considered. In the current study, the highest accuracy was obtained for the frequency ratio model, hence susceptibility map prepared from FR model was used as landslide potential (LP) factor. For the risk damage potential (RDP) the two important exposing factors were chosen as already discussed earlier. The factors are population and building which are most affected by the landslide disaster. The population density per sq. km in the current study varies from 10 to 3293 (Figure 6A). The areas with the higher density are more risk prone which is obviously understood. Similarly, the building density from the current analysisranges up to 190 building per sq. km irrespective of the sizes of the building (Figure 6B). Here in this study, we have preferred the numbers of building rather than their size. The mapped building was also validated with the building data of GSI. 40 building data were validated from the GSI and were obtained a 100% accuracy for these data (Figure 12). The final risk map prepared was classified into three zones as the high-risk zone, medium risk zone and the low-risk zone. From the present study it is obtained that 83% area falls under the low-risk zone, 14% falls under the medium risk zone and the 3% falls under the high-risk zone.

Figure 12. Map showing validation of building data.

4. Discussion

Preparation of landslide susceptibility map in the mountainous terrain is always a challenging task, as the selection of landslide conditioning factors is another challenge due to the thick forest cover leading to moist soil, slope hiding under natural shade, inaccessibility of high elevation and many other factors in the mountainous terrain (Sharma 2018). Considering all these factors, high-resolution remote sensing data (LISS-4) and ancillary data for the application of statistical models (FR and IV) and knowledge driven methods have been used for the generation of LSM for the present study area. Main causative factors were selected and used for the assessment of landslide susceptibility based on landslide polygons prepared and high-resolution remote sensing image interpretation. The main purpose of this study is to compare the probabilistic and the heuristic methods for the preparation of landslide susceptibility mapping.

A relation between the slope aspect and the slope failures were established using the IV and FR model, which infers that southerly facing slopes are more prone to landslides. Previous research done in the Himalayan terrain w.r.t slope aspect also reveals the same interpretation (Sarkar et al. 2013, Martha 2011). This indicates the probability of occurrence of landslides in these slopes is maximum. Because the Himalayan south-facing slopes receive sunlight throughout the day, there are more anthropogenic activities (Selvakumar et al. 2009; Rai et al. 2014; Bhardwaj et al. 2021). As anthropogenic activities increase, the slope becomes more vulnerable. The south west monsoon winds that travel west along the southern slopes of the Himalaya bring more rainfall to the southern slopes, which is also one of the major reasons for these failures (Fu et al. 2018). However, in case of AHP the maximum weights are obtained for north-west facing slopes followed by the south facing slopes, showing coherence with the probabilistic methods. The probable reason for north facing slopes contributing higher value is that the right bank of the main stream is dominant towards north (Bera et al. 2019).

In terms of slope gradient, it is always challenging to establish a perfect relationship between the slope angle and the slope failures. Previous research in Darjeeling Himalayas, however, shows that debris slides in the lower slope and rockslides typically occurred along steeper slopes (Kanungo and Sarkar2003; Ghosh et al. 2009). Nevertheless, in the present study, failures were not categorized, but the overall failures were considered, which shows that slopes with the gradient >40° are very vulnerable as per the ratings obtained from all the three models used here. This correlation supports the hypothesis that steeper slopes have a higher downhill component of force, causing gravitational sheer stress on gradient components and slope failure (Dai et al. 2001; Chapin et al. 2011). It is also found that the slope gradient influences the shear force, so with the increase of the slope gradient, the probability of the slope failures increases (Lee and Min 2001). Apart from this the component responsible for slope failure are also the cohesion of the slope material, angle of repose and the moisture condition (Selby 1993).

The curvature of the slope controls the inflow and outflow of water. Studies also reveal that curvature of a slope influences the moisture content dependent on saturation and erosion (Singh and Kumar 2018). Ridge sites had a positive curvature of the slope, which represented an upwardly convex surface. Places on the shallow and valley floor had a negative curvature, the surface is concave upwards (Fan et al. 2020). In the present study, it is observed that the convex slope are more susceptible to landslides followed by concave surface and flat surface are least susceptible and the same can be compared with the research done by other researchers (Gnyawali et al. 2020). Convex hill slopes, where seepage forces are directed upward from the steep, lower section of the slope, have the highest failure potentials in saturated homogeneous materials (Reid and Iverson 1992).

In general, the area with the high drainage density is more susceptible to slope failures (Singh and Kumar 2018), but in the present study, the negative correlation is seen for drainage density for the probabilistic method.Wang et al. 2003 also found in his research that low drainage density is most likely for landslide. He also mentioned that surface runoff is reduced by the high infiltration capacity, which erodes the ground surface and forms a gully. However, for the heuristic method, the positive corelation is seen showing higher weightage for the higher density value. The area with the lowest drainage is showing the highest weighted value. This means that local terrain influences are more important than drainage density alone. The other reason for this negative correlation is that the lowest class consists of the highest pixel values, i.e. it covers the maximum area, so probability increases. That is why it is always important to consider the multiple factors for these mappings.

It is generally found that the landslides in the Daling Group of rocks of Sikkim Himalaya are comparatively high (Sarkar et al. 2013). This group mainly contains interbedded chlorite-sericite schist/phyllite, feldspathic greywacke, biotite phyllite/mica schist, biotite quartzite, mica schist, variegated cherty phyllite, chert etc.However,in all the methods more weights were obtained for Everest Pelite. Everest Pelite Formation consists of PhyIlite, quartzite, biotite gneiss (Rawat 2015). The phyllitic, quartzite and gneissic lithological units in Everest Pelite has a high landslide density because it weathers quickly and produces less cohesive soil (Sarkar et al. 2013).All these rocks types have weak planes and get oversaturated with water during the continuous rainfall. As a result, these rocks have the highest probability of slope failures.

The human interference has accelerated the landslide incidence in the study area as supported by the results obtained from all the methods. In the Land use and land cover category, the Built-Up area shows the highest weighted value rating, followed by Fallow land and Forest cover. The other classes show a negative correlation. The highest rating in the Built-Up area is because of the unplanned and unscientific construction in vulnerable slopes (Jesiya and Gopinath 2019). The higher value obtained for the fallowlandas obtained in the current study can be attributed to the inherent properties of the LULC class (Anbalagan et al. 2015). On the other hand, the positive correlation on the forest cover is that according to predictive models, the area under forest cover is better correlated with landslides, which highlights the role of high precipitation in the Himalayan region (Sharma 2018).

NDVI suggests the characteristics of vegetation production, and it mainly changes the hydrological and soil shear strength characteristics of slope soil (Hong et al. 2017). The Normalized Difference Vegetation Index (NDVI) is one of the most important indicators used to determine vegetation cover over a wide range of time periods and in specific areas by employing remote sensing technology (Aburas et al. 2015). NDVI is the indexing that gives the idea about the vegetation cover in an area with the values ranging from −1 to +1, where +1 has the highest vegetation cover. The negative value generally indicates the cloud or snow cover. Above 0.3, the vegetation cover increases, and between 0.1to 0.3 generally represents the barren land, and zero value is for water bodies (Gandhi et al. 2015). The probabilistic analysis shows the positive corelation of NDVI with the occurrence of slope failures. So, in this scenario, the weighted value obtained shows a positive correlation. However, the heuristic approach is showing the negative correlation.

The elevation in the study area varies from 269 m to above 7000 m, but the more susceptible class, according to probabilistic method, are the classes where the elevation is ranging from 3500-4000 m and 4000-4500 m. The elevation class, mainly catchment areas with high surface runoff, tends to be more landslide-prone (Du et al. 2019). The spatial relation between landslides and elevation classes suggests that slope failures in the study region occurred mainly in locations with an elevation 3500-4500 m, often catchment areas with high surface runoff, which tend to be more landslide-prone. This is the reason why these two classes showing high weighted value, as the area with the elevation ranging from 3500-4500 m contains mostly catchment areas. Studies also reveal that in the hilly terrains, moderate elevations are more prone and susceptible to slope failures (Kose and Turk 2019). For the AHP method the higher values are shown for the higher elevation. Nakileza and Nedala (2020), in their research found similar results as of AHP method.

Lineaments of any types such as joints, foliation in rocks decrease resistance, material cohesion, and make surface water to percussion and it is the explanation for the high probability of failure in the area of structural discontinuity (Sarkar et al. 2013). A lineament has a minimal impact zone on the rock power relative to a thrust that has a broader region of impact (Martha et al. 2013). In the present study, lineament features include all the linear features, be it a joint or thrust. Both the geomorphic and the structural features taken together for the lineament density. This is the reason why the present study is showing the negative correlation with increasing lineament density for the probabilistic method. However, for the heuristic method the increasing lineament density is showing positive corelation with the occurrence of landslides.

Geomorphologically, the folded ridge and the highly dissected mountain regions have the highest potential for landslide occurrence. The hilly and mountainous regions are significantly more susceptible to landslides than plain and plateau regions (Roy et al. 2019). The same has been observed in the current study for all the methods used, with the highest weighted value in High Mountain followed by Narrow Valley, Rocky Cliff, and Escarpment. The reason for showing the high value in the Narrow Valley is that Narrow Valley contains steep sides, which makes the area susceptible to landslides.

The soil texture is strongly associated with the probability of a landslide incidence (Basu and Pal 2018). All soils fail under the effect of shear, where shear is a mixture of coherence and internal stress (Shooshpasha and Shirvani 2015). Cohesion decreases with an increase in sand content and increases with the clay content. Therefore, any mass of debris or siliceous sandy soil with poor cohesion is likely to fail when saturated on steep slopes (Bell 2013). Research also suggests that coarse-loamy, gravel-loamy, and loamy-skeletal soil texture classes indicate a higher risk of landslide occurrence compared to other soil texture classes (Roy et al. 2019). For the present study, the Coarse-loamy soil type shows the highest weighted value indicating the highly susceptible class in terms of landslide occurrence. Coarse loamy soil also contains a high percentage of sand as compared to clay content. Therefore, in the study area, steep slopes of over 40° occur frequent landslides consisting of soils with predominant sand fractions.

In the study the different comparisons have been made for the preparation of the susceptibility map keeping in the mind that both the statistical method (probabilistic) and the knowledge driven method (heuristic) are taken into consideration. The maps produced by using different methods for preparing landslide susceptibility map are shown in the Figure 13. The high accuracy in the probabilistic method is indicated of the unbiased in the analysis, however there can be some biasness in giving the ranking for heuristic methods. A comparatively less value in the ROC for AHP method might be the result of the biasness in the rating process. Therefore it is suggested to compare our results with the statistical method. Comparative study is very crucial when we have to prepare landslide risk zonation map, because the susceptibility map produced will be the landslide potential factor for the landslide risk zonation map. For the current study as the highest accuracy was obtained for the FR method, therefore the landslide susceptibility map obtained form it is considered for risk zonation. On the other hand, many researchers have used building and population as the risk damage potential for calculation of risk zonation map. Though the quantification is not done in the current study for the risk but it gives a very good qualitative risk zonation of the study area. Susceptibility map should not be the ultimate goal of the study. There can always be risk element associated with the susceptibility study. It is also very clear from the susceptibility map and the risk map that when the risk damage potential is associated, the outcome is different and the same is very obvious in the risk zonation map (Figure 14).

Figure 13. Landslide susceptibility map using different methods. A) FR method, B) IV method, C) AHP method.

Figure 14. Final risk zonation map of the study area.

5. Conclusions

The present study basically dealt with landslide susceptibility map using various methods, including the frequency ratio (FR) method and information value (IV) method for a probabilistic approach, and the analytic hierarchy process (AHP) for a heuristic approach. All methods performed well within acceptable accuracy limits during validation. The ultimate objective was to create a landslide preliminary qualitative risk zonation for the study area. The FR method, which yielded the best fit landslide susceptibility model (LSM), was used for the risk map. The accuracy, as measured by the receiver operating characteristic curve, was 77% for FR, 74% for IV, and 57% for AHP. The susceptibility maps were categorized into three classes—high, medium, and low susceptibility—using natural breaks. All three models indicated that land use/land cover (LULC), lithology, and geomorphology factors had the highest influence on landslide susceptibility. Specifically, the built-up class had the highest weightage among LULC, Everest Pelite among lithology, and high mountain class among geomorphology. Such studies also aid in prioritizing causative factors. After susceptibility mapping, the preparation of the risk zonation was crucial. Consistent with previous research, population and buildings were considered as indicators of risk and damage potential, given their exposure and vulnerability. The building map, created using the normalized difference built-up index (NDBI) and validated against GSI data, showed good accuracy. The resulting risk zonation map was divided into three zones: high-risk, medium-risk, and low-risk. The findings revealed that 3% of the study area fell under the high-risk zone, 14% under the medium-risk zone, and the remaining 83% under the low-risk zone (Figure 15). This type of research is highly valuable, particularly for planners, policymakers, and stakeholders involved in the area’s planning and development.

Figure 15. Pie chart for risk zonation percentage.

Compliance with ethical standards

Disclosure of potential conflicts of interest

Ethical approval

This article does not contain any studies with human participants or animals performed by any of the authors.

Contributions

P.B. and V.J.- Writing Original Draft, Conceptualization, Formal analysis, Investigation, Methodology, Software, Visualization, Review & Editing; I.R.C -Review & Editing, H.G.A.- Funding allocation, Review & Editing H.A., A.A.A., and M.M. - Review & Editing. All authors have read and agreed to the published version of the manuscript

Acknowledgements

PB is thankful to Guru Gobind Singh Indraprastha University, India for financial support through Indraprastha Research Fellowship vide F.No. GGSIPU/DRC/2019/1683. All the authors are thankful to the Dean, USEM, Guru Gobind Singh Indraprastha University, India for providing the facility in carrying out the research work.The authors are thankful to the editors and potential reviewers.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available on request from the corresponding author.

Additional information

Funding

This project was funded by Princess Nourah bint Abdulrahman University Research Supporting Project Number PNURSP2022R241, Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia.