Landslide susceptibility assessment in a lesser Himalayan road corridor (India) applying fuzzy AHP technique and earth-observation data

Abstract

The Kalsi-Chakrata road corridor, located in the Lesser Himalayas, experiences several landslides every year, resulting in a considerable amount of damage to roads, assets, and other infrastructure and even loss of lives. During the monsoon season (June–August), the disruption of routes due to landslide resulted in economic losses and barred villagers from accessing critical and essential facilities, thus impacting the livelihood of the communities residing along the road corridor. Hence immediate requirement was to systematically assess the landslide susceptibility for the study area that would support in the preparation of the planning and mitigation goal, both short and long term. The present study adopted the fuzzy analytic hierarchy process (fuzzy AHP) method integrated with geospatial technology that may be highly effective for landslide susceptibility assessment in the landslide-prone Lesser Himalayas. The use of validated landslide inventory data and high-resolution remote sensing images for selection and mapping of landslide conditioning factors gratifies the geospatial aspect of local variations across the Kalsi-Chakrata road corridor. For fuzzy AHP model setup, prominent landslide contributing factors viz., slope, aspect, altitude, lithology, proximity to road, fault & drainage, Stream Power Index (SPI), Topographical Wetness Index (TWI), rainfall, land use/land cover (LULC), soil, and seismicity, were mapped and classified into significant classes. The resultant landslide susceptibility map (LSM) shows that about 55% (45.23 km²) of the study area was categorized as a very high and high landslide susceptibility zone. Of this, about 21% (17.7 km²) was within very high LSM zone, while 33% (27.5 km²) fell under high landslide susceptibility categories. Approximately 17.6% (14.5 km²) areas fall within the moderately susceptible zone, where chances of future landslides may be amplified without periodic observation and prospective study. At the village level, it was observed that Jhutaya village, located nearer foothill of Lesser Himalayas, is most susceptible to landslide followed by Dhaira, Chapanu, and Sairi villages. The fuzzy AHP model shows 86.52% accuracy in landslide prediction evaluated through the ROC curve (at 95% confidence level). Hence, the output LSI may be referred by the planners and engineers for monitoring, prevention, and mitigation of landslide hazards and development of infrastructure in the Kalsi-Chakrata road corridor, and the fuzzy AHP methodology adopted may be applied in other areas of the Lesser Himalayas.

Keywords:

• Landslide susceptibility index (LSI)
• fuzzy analytic hierarchy process (fuzzy AHP)
• lesser Himalayas
• road corridor
• GIS

1. Introduction

In the Himalayan terrain, landslides are among the most common natural hazards that carry a significant threat to the socio-economic development in an area and the life and livelihood of the communities residing in those mountainous regions. The statistics published by the Geological Survey of India (GSI) shows that about 0.14 million sq km area in North West Himalayas (Uttarakhand, Himachal Pradesh, and Jammu & Kashmir) are prone to landslide hazard. It is estimated that approximately 13% of Indian land has witnessed landslide events that had caused a loss of over 4.5 million USD (Thai Pham et al. 2018). During the monsoon season (June–August), large numbers of landslides are triggered in the Himalayan mountain chain due to its unique geological setup and tectonic movements. The steep slopes, deep weathering of substrata, and intense southwest Monsoon rainfall in the Lesser and the Middle Himalayas in conjunction develop complex geological conditions that are favorable for the frequent occurrence of landslides (Gabet et al. 2004; Ambrosi et al. 2018). In Uttarakhand, a northern state of India, the landslides are among the most critical natural hazards that cause considerable damages every year. Most often triggered either by heavy precipitation and earthquakes, or anthropogenically due to heavy vehicles' movement or, usually, human activities on unstable slopes. In addition to rugged terrain, high seismicity, and climatic phenomena that are the main triggering factor for landslide occurrence (Kwan et al. 2014; Li et al. 2014), the hydro-geological conditions, and anthropogenic activities are among the key elements that too strongly influence the landslide occurrences in the Himalayan region of India (Rai et al. 2014; Singh et al. 2020).

The Himalayas are an extensive range of landmass shadowing the subcontinent of India that differs from one region to another regarding topography, climate, soils, and rock types. The seismically active Lesser Himalaya, a comparatively recent formation, is bounded by two tectonic planes, i.e., the Main Central Thrust (MCT) in the north and the Main Boundary Thrust (MBT) in the south. Due to intense folding and faulting, this rugged and hilly terrain is highly fractured, sheared, less consolidated, and highly weathered (Pachauri 2010). Moreover, compared to other Himalayas regions, the Lesser Himalayas experience a higher degree of human interference through construction activities, cultivation, deforestation, construction of roads cutting the slopes, more vehicular movements, mining, and other economic activities. These factors make the Lesser Himalayan terrain fragile and weak and often aggravate landslides phenomena, particularly during the Monsoon (Khanduri 2017).

Located in the Lesser Himalayas on the right bank of Amlana River, the Chakrata-Kalsi road corridor is 42 km long corridor situated in the western part of the Dehradun district of Uttarakhand state bordering with Himachal Pradesh. This road corridor well connects the famous hill station of Chakrata with the neighboring densely populated northern plains. During the southwest monsoon season (July–September) every year, the Kalsi-Chakrata road corridor experiences several landslides that incur a considerable amount of damage to roads, buildings, and other assets, even loss of lives. Some of these landslides are even recurring every year (e.g., Amroha landslide). The incidents reported in the national and local newspapers show that this road corridor is frequently blocked by landslide debris in July and August every year. The most common types of landslides in this area include debris slide, debris flow, rock slide, and rockfall. Besides direct damages in the form of property loss and infrastructure disruption, the indirect damages due to landslide include disruption of road transport, interruption, and losses in economic activities (rotten fruits and vegetables) and livelihood of the communities residing in those areas. Significantly, blocked, and sometimes washed out village roads left the mountain communities in a challenging situation barred from access to critical and essential facilities (Figure 1).

Figure 1. Landslide in 2018 razed the village access road and temple near Amroha.

In the last two decades, the Lesser Himalayas has witnessed a substantial increase in developmental activities, which was proven to be disastrous, coupled with extreme climate events (Allen et al. 2016). Therefore, the requirement was obvious to conduct a thorough investigation at a finer resolution applying a reliable statistical technique. It was observed that the area, particularly between Kalsi and Sahiya town, experienced recurring landslides every year, disrupting the road connectivity and often bringing the life to a complete stand-still. The literature survey available on landslides in the Kalsi-Chakrata road corridor reveals that so far, no in-depth and holistic scientific research has been conducted for this area except one published paper on terrain classification and land hazard mapping (Choubey and Litoria 1990). The present investigation was carried out at a sub-regional scale by generating base maps for land hazard evaluation using only five conditioning factors. Hence, to fill in the voids, it is necessary to assess detailed landslide susceptibility applying advanced and predictive contemporary statistical models using advanced geospatial techniques to detect the landslide susceptible areas that require immediate attention.

The predictions associated with landslide studies are always subjected to certain degrees of uncertainties. Hence, investigations concerning ‘when’ ‘where’ and ‘what magnitude’ is considered crucial aspects of prediction (Singh and Sharma 2015). In this regard, the landslide susceptibility mapping (LSM) may be a significant step towards landslide management and mitigation. It effectually helps to delineate the spatial distribution of probable landslide occurrences. The investigation of LSM has been carried out using different prediction methods worldwide, however, their accuracy and predictability vary widely (Chen et al. 2016). Therefore, it is notable to select a present-day approach that locally addresses the landslide issues considering the above speculation. There has been quite a lot of proper research works conceded in the prior two decades to assess landslide susceptibility based on both qualitative and quantitative methods (Varnes 1984; Guzzetti et al. 1999; Pradhan and Lee 2010; Tien Bui et al. 2012, 2015; Thai Pham et al. 2017; Marinos et al. 2019). In recent years, statistical and heuristic landslide prediction models integrated with Geographic Information System (GIS) have been instrumental for improved mapping of hazard footprints. More often, it includes modeling approaches, for example, multi-criteria evaluation, index of entropy, weights of evidence, statistical information value, linear and logistic regression, discriminant analysis, frequency ratio, fuzzy systems, and neuro-fuzzy models, neural networks, random forests (Guzzetti et al. 1999; Lee and Pradhan 2007; Tien Bui et al. 2015; Chen et al. 2017; Thai Pham et al. 2017; Ambrosi et al. 2018; Noorollahi et al. 2018; Taalab et al. 2018; Sur and Singh 2019; Sharma et al. 2020; Singh et al. 2020). It should be kept in consideration that the reliability and accuracy of LSM of a particular area depend on the spatial and temporal frequency of landslides and the modeling methodologies adopted (Glade 2001, Pellicani et al. 2014, Wang et al. 2020).

Saaty introduced the analytic hierarchy process (AHP) as a valued tool for multi-criteria decision-making (MCDM) adopted as a thriving natural resource management technique and environmental and ecological impact assessment. Literature shows that AHP can be successfully applied in natural resource management, ecological studies, susceptibility mapping, planning, and development (Saaty 1990; Ouma and Tateishi 2014; Hanine et al. 2016; Mallick 2016). However, in practice, while dealing with complex natural environments, the crisp data are often found inadequate to represent the system's relationships and uncertainties (Mallick et al. 2018). In the fuzzy AHP method, the pairwise comparisons or two-dimensional priority matrix through a linguistic variable are performed for both criteria and the alternatives. Such pairwise comparisons are represented by triangular numbers (Ayhan 2013; Mallick et al. 2018). Hence, in many cases, the experts prefer to choose intermediate pairwise comparison ratings for factors in AHP instead of definite ratings because, in complex natural systems, the evaluation method of fuzzy set theory shows a better representation of human experiences, judgments, and thinking styles addressing real-world complications (Gorsevski and Jankowski 2010; Hanine et al. 2016). Zadeh (1965) introduced the fuzzy set theory as a tool to address complex problems. The fuzzy set theory is based on a membership degree for specific attributes of interest known as membership function (Shahabi et al. 2012). Because of this, the application of the theory can be found in many disciplines such as data processing, artificial intelligence, decision support, engineering, medicines, expert systems, operational research, and others (Zimmermann 1992). Over the years, fuzzy AHP has evolved as an imperative MCDM technique in a fuzzy environment where the best choices can be selected. The weight of criteria can be determined according to pairwise comparisons and hierarchical relationships between criteria (Mallick et al. 2018). In this study, the fuzzy set theory was adopted to determine the fuzzy membership function (MF) based on the conditioning factors' spatial relationships.

Several procedures are discussed in various literature to compute the factor weights using the fuzzy AHP technique (Bozbura et al. 2007). Laarhoven and Pedrycz (1983) were among the researchers who applied fuzzy AHP for the first time. Among the successors, Buckley (1985), and later Chang (1996) were phenomenal. Chang introduced a new method based on the degree of probabilities of each criterion using triangular numbers in a pairwise comparison. The next step is to develop a pairwise construction matrix. This task is carried out for a particular level on the hierarchy by placing triangular fuzzy values for the linguistic variables. By comparing each couple, membership functions are then constructed, and the intersections' are determined. The membership value here represents the degree of possibility of the value. The minimum degree of possibility of the situation defines a particular criterion's weight. Subsequently, weights for each criterion are normalized to get the final weights (also known as final importance degrees) for the hierarchy level (Chang 1996). The literature review suggests that LSM generated through the AHP technique was adopted quite frequently by the researchers due to easier understanding and computation. The author observed much literature on applying the AHP technique on LSM in India; however, using the fuzzy AHP technique in this context was less, incredibly scarce for the Himalayas. Moreover, although multiple LSM studies are conducted for Uttarakhand's main highways connecting the "four abodes," scientific research lacks smaller road corridors such as Kalsi-Chakrata are experiencing frequent landslides every year. Hence, the effort has been made to conduct a detailed investigation to delineate landslide susceptible areas applying a GIS-based fuzzy AHP approach following Chang's techniques. The outcome of this study provides a substantial prediction of future landslide susceptible areas in the Kalsi-Chakrata road corridor that can help decision-makers and planners emphasize on appropriate mitigation measures that can lead to the prospective future economic growth of the area.

2. Study area

The Kalsi-Chakrata road corridor lies in the Lesser Himalayan ranges at the western end of the Central Himalayas. Chakrata (2,118 m above mean sea level) is a small hill station and a famous tourist destination, and Kalsi is primarily a tiny Cantonment Town located in the bank of river Yamuna in the foothill of the Lesser Himalayas. The Kalsi-Chakrata road corridor is about 42 km long and located in the Dehradun district of Uttarakhand state (Figure 2). The state highway starts ascending from Kalsi towards Chakrata. The road corridor is linked with Dehradun (Capital of Uttarakhand state) and Saharanpur (a market hub in Uttar Pradesh state) through Vikasnagar-Kalsi road.

Figure 2. Location of the Study Area.

The study area's whole terrain has a very rugged topography with altitude varying from 433 m near Yamuna River and 2456 m near Chakrata town (mean = 1445 m, std. dev. = 585 m). The temperature ranges from 1.4° C in winter to 34° C in summer, with an average annual rainfall varies between 1,339 cm and 1,479 cm. The study area typically experiences three seasons. The winter season starts in October and ends in February, and the summer months extend from March to June. The rainy season begins in July, and it ends in September. Often this area gets heavy rainstorms, sometimes cloudbursts, in July and August, which cause lots of damages resulting in erosion, floods, and landslides. The road corridor is situated in the Siwalik and Middle Himalayas, one of the most seismically active regions. The geology, structures, and lithology of the area have been mapped using geological maps acquired from GSI. In the study area, the outer sedimentary belt consists of the Krol belt, while the inner belts constitute thick sequences of unmetamorphosed sediments. There are seven geological formations observed in the study area, of which the Jaunsar group is the predominant formation that covers more than 90% of the geographical area. During the field investigation, the Mandhali-Chakrata Formation (Proterozoic II) and Chandpur Formation (Proterozoic III lower), where most of the landslides occur, were studied concerning other landslide conditioning factors. The Main Boundary Thrust (MBT), passes through the road corridor, can be distinguished after crossing Kalsi town moving towards Chakrata.

In this study, 34 Census villages and one municipal corporation (Chakrata) intersect with the Kalsi-Chakrata road corridor, have been considered for susceptibility mapping. According to the Census of India (2011) village boundary, the study area spreads over more than 80 sq. Km. with more than 19,679 persons residing in 3,961 households. Agriculture and social forestry play an essential part of livelihood among the villagers. The fruits and vegetables produced here have been transported to the nearby wholesale market in Saharanpur, further connected with the large vegetable market of Delhi.

3. Data, materials and methods

In any LSI study, the first important step is to detect the past and the present landslide locations in an area, followed by the determination of conditioning factors such as terrain, geology, hydrology, climate, and anthropogenic activities that are liable to contribute to the occurrence of landslides (Guzzetti et al. 2012; Deng et al. 2017; Ambrosi et al. 2018). Next, based on literature review and field data analysis, the critical conditioning factors (geo-environmental factors) that might have been responsible for triggering a past landslide have been identified for further consideration in predictive models (Marinos et al. 2019). The geo-environmental factors commonly considered for LSM include terrain, geology, hydrology, climate, and anthropogenic activities (Wang et al. 2019). More emphasis was put to select the high resolution and latest vintage landslide conditional data layers from authentic sources judiciously. These data were duly validated before inclusion in the chosen statistical model through available literature and published open-source information. The following sub-sections briefs the data and materials used and the methodology adopted for fuzzy AHP based LSM:

3.1. Data types, resolution and sources

The data types for various conditioning factors and their resolution, vintage and sources used in this study are the following (Table 1).

Table 1. Data and materials used with their resolution and sources.

Categories	Thematic data	Data & Sensor (resolution/scale)	Data Source (vintage)
Landslide inventory	Landslide inventory	Cartosat (2.5m)	National Remote Sensing Centre (NRSC), India (2017)
		Resourcesat-2/ Linear Imaging Self-Scanning System IV (LISS-IV) (5.8m)
		Google Earth, Published reports	Public Works Department (PWD), GSI portal (http://bhukosh.gsi.gov.in/Bhukosh/Public), Published reports (2001–2017)
Digital terrain model (DTM)	DTM	Cartosat (2.5m), Resourcesat-2/ LISS4 (5.8m)	National Remote Sensing Centre (NRSC), India (2017)
Topographical	Slope	Digital terrain model (10m)	National Remote Sensing Centre (NRSC), India (2017)
	Aspect
	Elevation
	Topographical wetness index (TWI)
Geological	Lithology	Geological map (1:25,000)	Geological Survey of India (GSI) (2015)
	Distance to faults
Hydrological	Distance to drainages	DTM (10m), Toposheet (1:50,000)	National Remote Sensing Centre (NRSC), India (2017), Survey of India (SOI) (2010)
	Stream power index (SPI)	DTM (10m)	National Remote Sensing Centre (NRSC), India (2017)
Meteorological	Rainfall	Rainfall Data (60 years daily records)	Indian Meteorological Department (IMD) (1947–2017)
Pedology	Soil	Soil map (1:25,000)	National Bureau of Soil Survey (NBSS), India, (2010)
Seismic	Ground acceleration	Seismic Zone map, Average shear wave velocity at 30m depth (Vs30)	Bureau of Indian Standards (BIS)(2002) U.S. Geological Survey, 2017
Vegetation	Normalized difference vegetation index (NDVI)	Resourcesat-2/ LISS4 (5.8 m)	National Remote Sensing Centre (NRSC), India (2017)
Anthropogenic	Land-use/ Land-cover (LULC)
	Distance to roads

3.2. Landslide inventory mapping

The landslide inventory dataset helps to understand the factors or conditions that might be accountable for triggering a past landslide or signify the incidence of an active landslide (Marinos et al. 2019). In this study, an extensive landslide inventory dataset was prepared using two primary sources: (1) the landslide records collected from the GSI portal (http://bhukosh.gsi.gov.in/Bhukosh/Public) and the Public Works Department (PWD), Uttarakhand; (2) identification of possible past landslide events from higher resolution periodic earth observation data (LISS-IV satellite images captured in March 2017 and Google Earth satellite images between the year 2001 and 2017). The landslide extent was delineated as precisely as possible, along with metadata for future use. The records collected in this process were examined through extensive field surveys mostly along the road corridor, including the connecting village roads using the Global Positioning System (GPS) instrument and geo-tagged camera during the Monsoon and post-Monsoon periods (Figure 3). It was observed that the most common types of rain-induced landslides in this area include debris slide, debris flow, topple failure, rock slide, and rockfall. The landslide inventory comprises the type of movement, type of material, the slide's dimension, the damages/losses caused by past events. The most widespread and threatening landslide in the area was observed in Amroha, a rotational slide with steep slopes and road undercut areas. The analysis of satellite images and ancillary reports from PWD shows that the Amroha landslide has increased by more than ten times in size over the past two decades and presently covering an area of about 0.05224 sq km. For both modeling and validation, the landslide point locations (demarcated at the centroid of the delineated landslide areas) were used. Though there are no standard procedures for the choosing the training dataset (for model preparation) and validation dataset (for verification of model’s performance) (Guzzetti et al. 2012), however, the most commonly practiced ratio is 70/30 for training and validation samples (Tsangaratos and Ilia 2016, Deng et al. 2017, Nachappa et al. 2020). The compiled landslide inventory finally comprises 55 records as of June 2019, and these were divided randomly into 38 training data (70% of the inventory) and 17 validation data (30% of the inventory).

Figure 3. Field photographs showing (a) rock fall on Kalsi-Chakrata Highway in 2018 (b) slope instabiliy at a landslide site near Sahiya (2017) (c) warning board (rarely found) for frequent rock fall (2017) (d) major landslide area in Amroha (2017).

3.3. Selection and mapping of landslide conditioning factors

The landslide conditioning factors (also called geo-environmental factors), without any predefined criteria for selection, are picked up based on a thorough literature review, and the numbers may vary from a few factors (Pradhan and Lee 2010; Akgun 2012) to a large number of factors (Catani et al. 2013; Tien Bui et al. 2015; Xu et al. 2016; Zhang et al. 2018; Nachappa et al. 2020). In the current investigation, 14 landslide conditioning factors were carefully selected based on hydro-geomorphological characteristics, literature, and expert consultation. These are slope, aspect, rainfall, lithology, normalized differential vegetation index (NDVI), soil, altitude, proximity to road/drainage/fault, seismicity, land-use, and land-cover (LULC), topographic wetness index (TWI), and stream power index (SPI). The conditioning factor thematic layers were arranged from authentic sources using remote sensing and geographic information system (GIS) software with a pixel size of 10 m and categorized applying Jenks natural breaks. Among these factors, the digital terrain model (DTM) was used to produce the thematic layers for topographical parameters such as altitude, slope, aspect, TWI, SPI. The DTM was generated from Cartosat digital elevation model (DEM) (10 m resolution) and LISS-IV satellite images (5 m resolution) in Global Mapper software at a 10 m resolution. The DTM was validated through the Survey of India (SOI) spot heights, GPS reading collected during field visits, and differential global positioning system (DGPS) points already used in other published reports in the Chakrata-Kalsi area. The remaining conditioning factors were generated using ArcGIS 10.8 software from secondary sources selected carefully based on their authenticity, scale, and vintage (Table 1). See Figure 4 for all the conditioning factors.

Figure 4. Thematic maps of the conditioning factors used in this study: (a) slope, (b) aspect, (c) Altitude, (d) rainfall, (e) lithology, (f) distance to roads, (g) distance to drainage, (h) distance to fault, (i) NDVI, (j) LULC, (k) seismicity, (l) soil, (m) TWI, (n) SPI.

The slope, directly influences shear strength, is one of the most critical factors accountable to the slope failure (Zhang et al. 2018; Thai Pham et al. 2019; Wang et al. 2019). The topographies with steeper slopes are usually more susceptible to failure. The slope angle was produced from the DTM in raster format and categorized into five distinct slope classes of <10°, 10°–20°, 20°–35°, 35°–50°, >50° applying Jenk’s Natural Breaks method.

The aspect has a significant impact on vegetation growth. The vegetation cover retains the soil's moisture capacity, which can instigate landslide (Das et al. 2011; Basharat et al. 2016). Aspects represent the slope direction. In this study, the aspect map generated from the DTM was classified into nine distinct sub-classes as flat (-1), north (337.5°–360°, 0°–22.5°), northwest(22.5°–67.5°), east (67.5°–112.5°), southeast (112.5°–157.5°), south (157.5°–202.5°), southwest (202.5°–247.5°), west (247.5°–292.5°) and northwest (292.5°–337.5°).

The variation in altitude influences the geomorphological and biological processes and is often attributed to spatial landscape variability (Raja et al. 2018; Sur and Singh 2019). The elevation may be correlated to different environmental settings, such as fluctuations in rainfall and vegetation types variation. Elevation was categorized into five classes: 39–500 m, 500–1000 m, 1000–1500 m, 1500–2000 m, and 2000–2660 above mean sea level.

For landslide initiation, rainfall is often considered the main triggering factor as it affects slope stability through surface overflow and pore water pressure (Wu et al. 2016; Wang et al. 2019). The spatiotemporal differences in rainfall quantity and intensity can instigate landslides both over large and small areas (Nachappa et al. 2020). The daily rainfall records for eight sub-districts headquarter stations surrounding the road corridor were obtained from the Indian Meteorological Department (IMD) for the past 70 years (from 1947 to 2017). Using the inverse distance weighted (IDW) method in ArcGIS 10.8 software, the mean annual rainfall thematic layer was produced having five classes: 1339–1371 mm, 1371–1398 mm, 1398–1426 mm, 1426–1453 mm, and 1453–1479 mm. It was observed that about 75% of the landslides occurred in the high and very high rainfall classes.

Lithology is another widely accepted and significant conditioning factor in LSM studies (Chen et al. 2017). Lithological units determine the strength and permeability of rocks. From the literature, it is evident that landslides often originate along with a rock stratum with less shear strength and less permeability. We have eight lithological units captured from the geological map. One significant observation from the map is the Mandhali-Chakrata Formation and Chandpur Formation that covers large areas along the road corridor containing highly weathered and less amalgamated greywacke rock materials. The investigation also depicted that the majority of past landslides occurred in the Mandhali-Chakrata Formation and Chandpur Formation.

The anthropogenic factors like distance to roads have been considered a vital conditioning factor for landslide occurrence (Meena et al. 2019). It has been observed that a considerable number of landslides occur across the road cuts mainly due to slope excavation, additional hydrologic load change, deforestation that affects the stability of the slope (Wang et al. 2019). The distance to roads was categorized into six classes: 0–100 m, 100–200 m, 200–300 m, 300–400 m, 400–500 m, and >500 m. The interval was chosen based on the fact that most of the landslides occurred within a distance of 500 m due to road construction activity. In this regard, it may be kept in consideration that there is no standard procedure to determine the distance to roads, drainage, faults as landslide conditioning factor in the literature. By changing the distance to any of these parameters, the output will be different. Hence, based on expert opinion, past studies, and observations during the field investigation, the distance parameters are usually determined for the LSM studies as carried out in this study too.

The distance to drainage controls the slope's stability due to runoff and leads to landslides by eroding slopes (Wang et al. 2019). The presence of drainage near landslides has comparatively more impact on landslide initiation than areas located far-off (Bera et al. 2019). The effect of streams on the slope failure was correlated using five drainage buffer zones were generated, each 100 m wide, using the Euclidean distance, namely, 0–100 m, 100–200 m, 200–300 m, 300–400 m, and >400 m.

The distance to faults is another influential factor for landslide activities as the tectonic breaks usually reduce the adjacent rock strength (Meena et al. 2019). With an increase in distance from faults, the number of landslides decreases sharply (Xu et al. 2016; Wang et al. 2019). The faults were delineated from the geological map obtained from GSI. The tectonically active Himalayan Main Boundary Thrust (MBT) passes through the road corridor, and the same can be detected passing near to Kalsi across the road corridor. The distance to faults was categorized into an interval of 500 m, namely <500 m, 500–1000 m, 1000–1500 m, 1500–2000 m, and >2000 m.

NDVI is another important conditioning factor that strongly influences the likelihood of land sliding. Bare earth without vegetation coupled with steeper slope and higher rainfall is more susceptible to frequent landslides (Singh and Sharma 2015; Sur and Singh 2019). The NDVI values (near to −1) specifies that the bare earth surface devoid of vegetation, while value (near to +1) indicates a higher and healthier vegetation cover. The NDVI data was generated from LISS-IV satellite imagery and classified into five classes: −0.24–0.19, 0.19–0.31, 0.31–0.39, 0.39–0.48 and 0.48–0.77.

LULC directly or indirectly reflects the natural process and human-made activities that influence the slope stability in an area. Usually, the land cover like forested areas has fewer landslides than that of barren lands. In this study, The LULC layer was generated from the LISS-IV satellite image applying the pixel-wise classification technique following the classification scheme followed by the National Remote Sensing Centre (NRSC), which was further validated with available SOI toposheets and NRSC maps. There are 14 distinct LULC classes derived from the study area. It has been noted that 67% of the landslides incidents occurred in LULC classes of sparse vegetation, scrubland, and barren rocky areas where weathering of rocks and erosion is more intense.

The effect of seismicity in the Lesser Himalayas is another critical factor for earthquake-induced landslides (Pachauri 2010). In this study, the Bureau of Indian Standards (BIS) Earthquake map was referred to produce the seismicity layer in GIS, which was further modified after USGS Vs30 data (the time average shear wave velocity at 30 m depth). The final thematic layer had five classes with the ground acceleration varying between 3.5 and 3.84 m/s².

Soil characteristics such as clay content in the mountainous terrain become potential slip zone that can lead to slope failure resulted in landslides . The fundamental soil properties such as texture, organic matter content, structure, and permeability greatly influence the soil erodibility. Hence, the soil layer, generated from the National Bureau of Soil Survey (NBSS) soil map, was categorized into four distinct categories based on the depth of soil cover, soil texture, morphology related to slope, and drainage, and used in this study. These are- (i) moderately shallow, excessively drained, thermic, loamy skeletal soils on moderately steep slopes with loamy surface, moderate erosion and moderate stoniness (found mainly in summits and ridge tops); (ii) moderately deep, excessively drained, thermic, loamy skeletal soils on moderately steep slopes with loamy surface and moderate erosion (found in side slopes); (iii) moderately deep, well drained, thermic coarse loamy soils on moderate slopes with loamy surface and moderate erosion (found in fluvial valleys); and (iv) moderately shallow, excessively drained, thermic, coarse loamy soils on steep slopes with loamy surface, moderate erosion, and strong stoniness (found in side slopes).

Topographic wetness index (TWI) is a ‘steady-state wetness index; that describes topography's influence. It identifies the water saturation zones in an area when there is a high catchment area with a low slope gradient. TWI, after Moore et al. (1988), is defined as: (1) $$\mathit{TWI}=ln\left(\frac{a}{\hspace{0.17em}\text{tan}\hspace{0.17em}\beta }\right)$$(1) where, $ln$= Natural logarithm; $a$= Flow accumulation/specific catchment area of each cell; $\beta$ = Slope gradient in degrees.

In this study, TWI was generated from DTM using ArcGIS software and divided into four categories, namely, 0–4, 4–8, 8–12, and 12–14.

The stream power index (SPI) represents the erosive power of the streams. The slope gradient in the catchment area and the streams' erosive power show a positive relationship. According to Moore et al. 1988 and Wang et al. 2019- (2) $$\mathit{SPI}=ln\left(a\ast \hspace{0.17em}\text{tan}\hspace{0.17em}\beta \right)$$(2) where, $ln$= Natural logarithm; $a$= Flow accumulation/specific catchment area of each cell; $\beta$ = Slope gradient in degreesgradient in degrees

The SPI values generated were further classified into two distinct classes ranging from 0–10 to 10–15.

3.4. Modeling approach for fuzzy AHP

As discussed in the introduction section of this paper, the fuzzy AHP modeling approach in the present study adopted the Chang’s (1996) extent method using triangular fuzzy numbers (TFN) to calculate the fuzzy weights for the selected landslide conditioning factors. The fuzzy AHP modeling was conducted in open source (LPGL 3) package for ‘R’ programming language using the ‘Fuzzy AHP’ package developed by Jan Kaha (version 0.9.5, 2019) (http://github.com/JanCaha/FuzzyAHP/) after necessary adjustment in the code. This package follows the fuzzy comparison matrix approach described by Siebert et al. (2016). While Figure 5 presents the workflow that was followed by the authors for generating the LSM for the study area. The following section briefly details the fuzzy AHP modeling steps followed.

Figure 5. Fuzzy AHP methodology adopted in this study.

3.4.1. Step-1: Development of hierarchical structure

In this step, the AHP method is adopted based on the literature survey and expert opinion specific to the study area. Next, a judgment matrix is developed using pairwise comparison, and decision-makers/experts are invited with their judgments. The matrices were examined for consistency, and eigenvalues and eigenvectors were calculated to compare the priority of the elements (Li et al. 2009). This step is followed by checking of consistency of the matrix judgments in the pairwise evaluation. If the consistency test fails, the pairwise comparison matrix's original values have to be revised (Table 2).

Table 2. Development of hierarchical structure.

code	Linguistic variables	Fuzzy Triangular Scale
		L	M	U
1	Equally important	1	1	1
2	intermediate value between 1 and 3	1	2	3
3	Slightly important	2	3	4
4	intermediate value between 3 and 5	3	4	5
5	Important	4	5	6
6	intermediate value between 5 and 7	5	6	7
7	Strongly important	6	7	8
8	intermediate value between 7and 9	7	8	9
9	Extremely important	9	9	9

3.4.2. Step-2: Determination of degree of membership and computation of fuzzy evaluation matrix

The fuzzy logic in decision-making uses the fuzzy set theory to determine the fuzzy membership function (MF) based on the conditioning factors' spatial relationships. In the fuzzy set theory, the objects can be assigned a membership value ranging between 0 and 1 that shows MF's degree. The scores of pairwise comparisons in Step-1 above are transformed into linguistic variables to assess alternatives under the fuzzy environment following Chang’s fuzzy AHP theory. Therefore, all the landslide conditioning factors were converted in the range of 0 (low susceptible) to 1 (high susceptible). In Chang’s fuzzy AHP method, the fuzzy MF and Si (synthetic extent value) were computed for the ith criteria through the following equation: (3) $$Si=\sum _{j=1}^m{M}_{g_i}^j\times {\left[\sum _{i=1}^n\sum _{j=1}^m{M}_{g_i}^j\right]}^{-1}$$(3) where, $g_i$ is the goal set for ($i$= $1,2,3,4,\dots \dots \dots ..n$) and all the $M_{g_i}^j$($j=1,2,3,4,\dots \dots \dots ..m$) are Triangular Fuzzy Numbers (TFNs).

Here, to obtain ${\sum }_{j=1}^m{M}_{g_i}^j,$ i.e., to perform ‘fuzzy addition’ operation of $\prime m\prime$ extent, sub-totals are calculated for each row of the matrix and new set ($l,m,u$) is obtained, where l is the lower limit value, m is the most promising value and u is the upper limit value. The equation can be presented as- (4) $$\sum _{j=1}^m{M}_{g_i}^j=\left[\sum _{j=1}^m{l}_j,\sum _{j=1}^m{m}_{j,}\sum _{j=1}^m{u}_j\right]$$(4)

Similarly, to get ${\left[{\sum }_{i=1}^n{\sum }_{j=1}^m{M}_{g_i}^j\right]}^{-1}$^, the inverse of the vector is calculated (Table 3).

Table 3. Mean pairwise comparison with respect to GOAL.

GOAL	Slope	Rainfal	Lithology	Aspect	NDVI	LULC	Road buffer	Drainage buffer	Fault buffer	Altitude	Soil	Seismicity	TWI	SPI
Slope	(1.000,1.000,1.000)	(1.000,1.414,3.000)	(1.000,1.189,3.000)	(1.000,2.913,7.000)	(1.000,2.450,4.000)	(2.000,3.223,5.000)	(2.000,3.223,5.000)	(2.000,4.160,6.000)	(2.000,4.820,7.000)	(2.000,3.873,6.000)	(2.000,4.054,7.000)	(2.000,4.605,7.000)	(3.000,4.866,8.000)	(4.000,5.624,9.000)
Rainfal	(0.333,0.707,1.000)	(1.000,1.000,1.000)	(1.000,1.189,3.000)	(1.000,2.214,4.000)	(1.000,2.213,4.000)	(1.000,3.080,6.000)	(1.000,2.449,4.000)	(2.000,3.463,5.000)	(2.000,4.354,7.000)	(1.000,2.632,5.000)	(2.000,3.661,6.000)	(2.000,4.160,6.000)	(2.000,4.738,8.000)	(2.000,5.208,8.000)
Lithology	(0.333,0.841,1.000)	(0.333,0.841,1.000)	(1.000,1.000,1.000)	(1.000,2.515,6.000)	(1.000,1.682,3.000)	(1.000,3.080,6.000)	(1.000,1.861,4.000)	(1.000,3.309,6.000)	(2.000,4.119,7.000)	(1.000,2.449,5.000)	(1.000,3.162,6.000)	(4.000,5.000,6.000)	(4.000,5.232,7.000)	(4.000,6.191,8.000)
Aspect	(0.143,0.343,1.000)	(0.250,0.452,1.000)	(0.167,0.398,1.000)	(1.000,1.000,1.000)	(1.000,1.000,1.000)	(1.000,1.682,3.000)	(1.000,1.414,3.000)	(1.000,2.213,4.000)	(1.000,2.711,4.000)	(2.000,3.223,5.000)	(2.000,3.722,5.000)	(2.000,3.934,6.000)	(2.000,4.527,8.000)	(2.000,4.739,8.000)
NDVI	(0.250,0.408,1.000)	(0.250,0.452,1.000)	(0.333,0.595,1.000)	(1.000,1.000,1.000)	(1.000,1.000,1.000)	(1.000,2.114,6.000)	(1.000,1.189,3.000)	(1.000,2.783,6.000)	(2.000,3.408,6.000)	(2.000,2.999,4.000)	(1.000,2.632,5.000)	(2.000,3.409,6.000)	(2.000,4.787,8.000)	(4.000,5.691,8.000)
LULC	(0.200,0.310,0.500)	(0.167,0.325,1.000)	(0.167,0.325,1.000)	(0.333,0.595,1.000)	(0.167,0.473,1.000)	(1.000,1.000,1.000)	(1.000,1.189,3.000)	(1.000,2.213,4.000)	(1.000,2.213,4.000)	(1.000,3.309,6.000)	(2.000,3.223,5.000)	(2.000,4.160,6.000)	(1.000,3.499,6.000)	(1.000,3.499,6.000)
Road buffer	(0.200,0.310,0.500)	(0.250,0.408,1.000)	(0.250,0.537,1.000)	(0.333,0.707,1.000)	(0.333,0.841,1.000)	(0.333,0.841,1.000)	(1.000,1.000,1.000)	(1.000,2.590,6.000)	(1.000,2.114,6.000)	(1.000,3.080,6.000)	(2.000,4.162,6.000)	(2.000,4.162,6.000)	(1.000,2.829,5.000)	(2.000,3.936,6.000)
Drainage buffer	(0.167,0.240,0.500)	(0.200,0.289,0.500)	(0.167,0.302,1.000)	(0.250,0.452,1.000)	(0.167,0.359,1.000)	(0.250,0.452,1.000)	(0.167,0.386,1.000)	(1.000,1.000,1.000)	(0.333,1.495,6.000)	(1.000,2.783,6.000)	(1.000,2.711,4.000)	(2.000,3.873,6.000)	(1.000,2.514,6.000)	(1.000,2.659,6.000)
Fault buffer	(0.143,0.207,0.500)	(0.143,0.230,0.500)	(0.143,0.243,0.500)	(0.250,0.369,1.000)	(0.167,0.293,0.500)	(0.250,0.452,1.000)	(0.167,0.473,1.000)	(0.167,0.669,3.003)	(1.000,1.000,1.000)	(1.000,2.000,3.000)	(1.000,2.213,4.000)	(1.000,2.783,6.000)	(1.000,2.783,6.000)	(1.000,2.990,6.000)
Altitude	(0.167,0.258,0.500)	(0.200,0.380,1.000)	(0.200,0.408,1.000)	(0.200,0.310,0.500)	(0.250,0.333,0.500)	(0.167,0.302,1.000)	(0.167,0.325,1.000)	(0.167,0.359,1.000)	(0.333,0.500,1.000)	(1.000,1.000,1.000)	(1.000,2.060,4.000)	(1.000,3.080,6.000)	(1.000,2.279,4.000)	(1.000,2.449,5.000)
Soil	(0.143,0.247,0.500)	(0.167,0.273,0.500)	(0.167,0.316,1.000)	(0.200,0.269,0.500)	(0.200,0.380,1.000)	(0.200,0.310,0.500)	(0.167,0.240,0.500)	(0.250,0.369,1.000)	(0.250,0.452,1.000)	(0.250,0.485,1.000)	(1.000,1.000,1.000)	(1.000,2.213,4.000)	(1.000,2.711,4.000)	(2.000,3.722,5.000)
Seismicity	(0.143,0.217,0.500)	(0.167,0.240,0.500)	(0.167,0.200,0.250)	(0.167,0.254,0.500)	(0.167,0.293,0.500)	(0.167,0.240,0.500)	(0.167,0.240,0.500)	(0.167,0.258,0.500)	(0.167,0.359,1.000)	(0.167,0.325,1.000)	(0.250,0.452,1.000)	(1.000,1.000,1.000)	(1.000,2.000,3.000)	(1.000,2.632,5.000)
TWI	(0.125,0.206,0.333)	(0.125,0.211,0.500)	(0.143,0.191,0.250)	(0.125,0.221,0.500)	(0.125,0.209,0.500)	(0.167,0.286,1.000)	(0.200,0.353,1.000)	(0.167,0.398,1.000)	(0.167,0.359,1.000)	(0.250,0.439,1.000)	(0.250,0.369,1.000)	(0.333,0.500,1.000)	(1.000,1.000,1.000)	(1.000,1.861,5.000)
SPI	(0.111,0.178,0.250)	(0.125,0.192,0.500)	(0.125,0.162,0.250)	(0.125,0.211,0.500)	(0.125,0.176,0.250)	(0.167,0.286,1.000)	(0.167,0.254,0.500)	(0.167,0.376,1.000)	(0.167,0.334,1.000)	(0.200,0.408,1.000)	(0.200,0.269,0.500)	(0.200,0.380,1.000)	(0.200,0.537,1.000)	(1.000,1.000,1.000)

3.4.3. Step-3: Evaluation of fuzzy index weights (degree of possibility of value)

In this step, the comparison of fuzzy values (degree of possibility) were computed. The degree of possibility of is determined using the following equations-

The degree of possibility of $M_2=(l_2,m_2,u_2)\ge M_1=(l_1,m_1,u_1)$ is defined by- (5) $$V\left(M_2\ge M_1\right)=\mathit{sub} y\ge x\left[\min \left(u{M}_1\left(x\right),u{M}_2\left(y\right)\right)\right]$$(5) (6) $$=\left\{{}^{\begin{array}{c}1,\\ 0,\end{array}}\frac{l_1-u_2}{\left(m_2-u_2\right)-\left(m_1-l_1\right)}\right\}\begin{array}{c}if m_2\ge m_1\\ if l_1\ge u_2\end{array}$$(6)

Otherwise, (7) $$d^{\prime }\left(A_i\right)=\min V\left(S_i\ge S_k\right)$$(7) where, $i=1,2,3,4,\dots \dots .,k$ and $d$ is the ordinate of the highest intersection point.

For $k=1,2,3,4,\dots \dots ,n$; $k\ne i,$ the weight vector is calculated by- (8) $$W^{\prime }={\left(d^{\prime }\left(A_1\right),d^{\prime }\left(A_2\right),d^{\prime }\left(A_3\right),d^{\prime }\left(A_4\right),\dots \dots \dots \dots ,d^{\prime }\left(A_n\right)\right)}^T$$(8) where, $A_i(i=1,2,3,4,\dots \dots .,n)$ are $n$ elements.

3.4.4. Step-4: Generation of normalized fuzzy decision matrix

Finally, through normalization of fuzzy decision matrix using maximum likelihood function, the normalized fuzzy weight vectors are calculated. (9) $$W={\left[\frac{d\left(A_1\right)}{\sum _{i=1}^{n}d\left(A_i\right)},\frac{d\left(A_2\right)}{\sum _{i=1}^{n}d\left(A_i\right)},\dots \dots ..,\frac{d\left(A_n\right)}{\sum _{i=1}^{n}d\left(A_n\right)}\right]}^T$$(9) where $W$ is a non-fuzzy number.

The weights calculated from the model were used for weighted sum overlay analysis in GIS. The final landslide susceptibility score was obtained for the study area by integrating various landslide factor weights through the following-

3.5. Landslide susceptibility mapping

The weights derived from the fuzzy AHP model were used for weighted sum overlay analysis in ArcGIS software (Version 10.8). The final LSI score was obtained for the study area by integrating various landslide factor weights through- (10) $$\mathit{Landslide}\mathrm{ }\mathit{Susceptibility}\mathrm{ }\mathit{Index}\mathrm{ }\left(\mathit{LSI}\right)=\sum _{j=1}^n{W}_j{w}_{ij}$$(10) where $W_j$ is the normalized fuzzy AHP weights for conditioning factor $j,$ $w_{ij}$ is the weighted value of sub-class $i$ in conditioning factor $j,$ and $n$ is the number of parameters.

i.e., LSI for each grid cell = (0.11 * Slope sub-class weights of the cell) + (0.104 * Lithology sub-class weights of the cell) + … … … … … …+ (0.01 * SPI sub-class weights of the cell)

The LSI values were further classified into five distinct LSM categories (very high, high, moderate, low, very low) using Jenk’s Natural Breaks. The results attained from the fuzzy AHP model were validated using the validation dataset using Receiver Operating Characteristics (ROC) curve and relative landslide density index (R-index).

4. Results

As presented in Table 4, concerning the goal defined in Chang’s fuzzy AHP, the slope is the top priority among the conditioning factors. The next priorities are rainfall, lithology, NDVI, aspect, road buffer, LULC, drainage buffer, fault buffer, altitude, soil, seismicity, TWI, and SPI obtained weights. More precisely, considering each factor weights derived applying fuzzy AHP for all 14 landslide conditioning factors, slope >50° is the most critical landslide conditioning class. Next in the priority classes are rainfall >1450 mm, NDVI < 0.2, slope 35°–50°, rainfall 1427–1450 mm, road buffer < 100 m, drainage buffer < 100 m, fault buffer < 100 m, NDVI 0.2–0.3, lithology class L1 having greywacke content, slope class 20°–35°, road buffer 100 m–200 m, drainage buffer 100 m–200 m, aspect class S and SW. Conversely, slope < 10°, LULC class of agriculture, soil class 13 (moderately shallow, excessively drained, thermic, loamy skeletal soils on moderately steep slopes), lithology class L-7 (micaceous sandstone, purple clay, mudstone), and aspect class flat are the conditioning factors having the least weightages.

Table 4. Normalized fuzzy AHP weights derived for the conditional factors.

Conditioning Factor	Class	Pixel per class	% of pixels	LSD Pixels by Class	% of Landslide pixels	Fuzzy AHP Weights	Fuzzy AHP factor weights	Consistency Ratio
Slope (°)	0–10	50604	6.09	38	0.86	0.00	0.110	0.0067
	10–20	123991	14.93	296	6.67	0.01
	20–35	453913	54.65	106	2.39	0.03
	35–50	196679	23.68	2054	46.25	0.03
	50–89	5406	0.65	1947	43.84	0.04
Aspect	Flat	5542	0.67	0	0.00	0.003	0.093	0.003
	North	33647	4.05	2	0.05	0.01
	NE and NW	171292	20.62	29	0.65	0.014
	S and SW	268986	32.38	2390	53.82	0.025
	E and SE	211828	25.50	389	8.76	0.018
	SE	139298	16.77	1631	36.73	0.023
Altitude (m)	39–500	13271	1.60	5	0.11	0.006	0.055	0.002
	500–1000	144276	17.37	1944	43.77	0.008
	1000–1500	366712	44.15	2480	55.84	0.013
	1500–2000	241042	29.02	12	0.27	0.013
	2000–2660	65292	7.86	0	0.00	0.015
Distance to roads (m)	<100	213244	25.67	2644	59.54	0.03	0.090	0.004
	100–200	147573	17.77	935	21.05	0.03
	200–300	111197	13.39	718	16.17	0.02
	300–400	80720	9.72	129	2.90	0.01
	400–500	56817	6.84	0	0.00	0.00
	>500	221042	26.61	15	0.34	0.00
Distance to drainage (m)	<100	169159	20.37	751	16.91	0.031	0.077	0.003
	100–200	148305	17.86	1075	24.21	0.026
	200–300	134641	16.21	878	19.77	0.016
	300–400	117972	14.20	761	17.14	0.004
	>400	260516	31.37	976	21.98	0.000
Seismicity (m/s2)	3.5–3.58	154076	18.55	725	16.33	0.00	0.030	0.002
	3.59–3.64	122010	14.69	2336	52.60	0.00
	3.65–3.71	179168	21.57	503	11.33	0.01
	3.72–3.77	153579	18.49	877	19.75	0.01
	3.78–3.84	221760	26.70	0	0.00	0.01
SPI	10–15	725142	87.30	4112	92.59	0.007	0.010	0.002
	0–10	105451	12.70	329	0.07	0.003
Distance to faults(m)	<500	206624	24.88	1225	27.58	0.030	0.069	0.003
	500–1000	177003	21.31	1115	25.11	0.024
	1000–1500	150568	18.13	1475	33.21	0.013
	1500–2000	113259	13.64	345	7.77	0.002
	>2000	183139	22.05	281	6.33	0.000
Rainfall (mm)	1297–1325	138104	16.63	0	0.00	0.002	0.104	0.003
	1326–1350	170283	20.50	282	6.35	0.011
	1351–1374	189931	22.87	835	18.80	0.022
	1375–1397	163582	19.69	631	14.21	0.032
	1398–1419	168693	20.31	2693	60.64	0.037
TWI	0–4	702059	84.53	3673	82.71	0.019	0.029	0.004
	4–8	112032	13.49	727	16.37	0.010
	8–12	11466	1.38	30	0.68	0.000
	12–16	5036	0.61	11	0.25	0.000
NDVI	0.50–0.76	110006	13.24	9	0.20	0.000	0.099	0.005
	0.40–0.50	191391	23.04	19	0.43	0.011
	0.3–0.40	319138	38.42	418	9.41	0.021
	0.2–0.3	196906	23.71	2514	56.61	0.030
	<0.20	13152	1.58	1481	33.35	0.037
Soil	Moderately shallow, excessively drained, thermic, loamy skeletal soils on moderately steep slopes	114589	13.80	14	0.32	0.001	0.047	0.002
	Moderately deep, excessively drained, thermic, loamy skeletal soils on moderately steep slopes	153160	18.44	0	0.00	0.009
	Moderately deep, well drained, thermic coarse loamy soils on moderate slopes	25287	3.04	121	2.72	0.016
	Moderately shallow, excessively drained, thermic, coarse loamy soils on steep slopes	537557	64.72	4306	96.96	0.021
Lithology	Shale, Phylite, Siltstone, Quartzite, Limestone, greywacke Conglomerate	441945	53.21	1776	39.99	0.027	0.104	0.004
	Greywacke, Shale, Quartzite, Dolomite, tuff with Dolerite	291387	35.08	2665	60.01	0.022
	Sandstone, Claystone, Siltstone, shale	17539	2.11	0	0.00	0.020
	Multicyclic sequence of brown to grey Silt clay with kankar and micaceous sand with pebbles	17040	2.05	0	0.00	0.016
	Brownish grey, clay, sand and gravel with boulders	8394	1.01	0	0.00	0.012
	Carbonaceous shale, slate, greywacke	10906	1.31	0	0.00	0.005
	Micaceous sandstone, purple clay, mudstone	13737	1.65	0	0.00	0.002
	Carbonaceous slates, quartzite, dolomite with stromatolite and cherry limestone	29564	3.56	0	0.00	0.000
LULC	River	6554	0.79	0	0.00	0	0.083	0.007
	Sandy area	1629	0.20	0	0.00	0.001
	Settlement (urban, rural, commercial)	4587	0.55	0	0.00	0
	Dense Veg	304205	36.63	22	0.50	0.007
	Forest	6853	0.83	0	0.00	0.008
	Agriculture	61753	7.43	0	0.00	0.003
	Sparse Veg	132072	15.90	576	12.97	0.013
	Rocky & Barren land	1801	0.22	200	4.50	0.019
	Open area	286311	34.47	2724	61.34	0.02
	Scrub land	24828	2.99	919	20.69	0.012

The geospatial technology integrated with the fuzzy AHP based MCDM technique focuses more on the correlation of landslide conditioning factors and the spatial distribution of landslide susceptibility. The resultant LSI data overlay on the village boundary depicts that as high as about 55% (45.23 sq km) of the study area fall within the very high and high landslide susceptibility zones (Figure 6). About 21% (17.7 km²) was under very high LSI, while 34% (27.5 km²) fell under high LSI categories (Table 4). Approximately 17.6% (14.5 km²) areas fall within the moderately susceptible zone. These are the areas where future landslides possibilities are high, which may be amplified without periodic observation and prospective study.

Figure 6. (a) Landslide susceptibility map applying fuzzy AHP model (b) Distribution of LSI zones.

This landslide susceptibility assessment technique provides a strong basis to comprehend the nature of the threat for the decision-makers, architects, and engineers to minimize the damages and losses and help mitigate risks caused by existing and future landslides at the village level. The analysis shows that Nevi village, located near Sahiya town, is most vulnerable to landslides. It has 61.6% of its areas fall within a very high LSI zone. The other villages that also show the higher distribution of very high landslide susceptibility are Chapanu (59.8%), Kakari (52.3%), Kanbua (48.9%), and Pajiya (46.7%). Conversely, the villages having very low LSI (less than 5% of the village area are highly susceptible to potential landslides) are Chakrata CB (1.1%), Bantar (2.3%), Vyas Bhood (3.6%), Rani (4.4%), and Mangrauli (4.8%).

Since the likelihood of landslide is notably higher in both very high and high LSM zones, the overlay analysis was performed to identify the villages that have more areas fall in those categories. The result specifies that Nevi (93.4%), Dadhau (92.7%), Chapanu (91.5%), and Kanbua (91.2%) are the villages that are likely to be more susceptible to future landslides. Table 5 presents the villages put in order as per the % of geographical areas under very high and combined very high & high LSI classes. The top ten villages having considerably high landslide susceptibility in descending order are- (i) Kakari, (ii) Kanbua, (iii) Pajiya, (iv) Koti (Bamrad), (v) Malaitha, (vi) Udapalta, (vii) Sainsa, (viii) Samalta Dadauli, (ix) Dadhau and (x) Sainda. In continuation of the present investigation, we recommend conducting a detailed landslide risk reduction assessment in these villages for appropriate preparedness and mitigation in future.

Table 5. Distribution of very high and high LSI areas by village.

S.N.	Village Name	Total Village Area ('000 sqm)	Very High LSI Area ('000 Sqm)	% Very High LSI Area	High LSI Area ('000 Sqm)	% High LSI Area	Combined Very High & High LSI Area ('000 Sqm)	% Combined Very High & High LSI Area
1	Nevi	1,234	760.10	61.6%	392.40	32%	1,152.50	93.4%
2	Chapanu	770	460.60	59.8%	244.00	32%	704.60	91.5%
3	Kakari	1,775	928.70	52.3%	653.00	37%	1,581.70	89%
4	Kanbua	1,620	791.70	48.9%	685.40	42%	1,477.10	91.2%
5	Pajiya	2,059	961.60	46.7%	798.30	39%	1,759.90	85%
6	Koti (Bamrad)	1,189	521.10	43.8%	480.20	40%	1,001.30	84%
7	Malaitha	3,359	1,452.90	43.3%	1,533.00	46%	2,985.90	89%
8	Udapalta	3,527	1,404.00	39.8%	1,469.90	42%	2,873.90	81%
9	Sainsa	1,583	625.70	39.5%	672.90	43%	1,298.60	82%
10	Samalta Dadauli	2,861	997.50	35%	1,294.50	45%	2,292.00	80%
11	Dadhau	536	185.70	35%	311.10	58%	496.80	92.7%
12	Sainda	1,386	459.80	33%	694.00	50%	1,153.80	83%
13	Rikhar	951	314.90	33%	345.60	36%	660.50	69%
14	Kwarana	2,633	868.50	33%	1,017.40	39%	1,885.90	72%
15	Sakani	745	240.00	32%	361.50	49%	601.50	81%
16	Nithala	759	243.70	32%	292.40	39%	536.10	71%
17	Jhutaya	803	256.70	32%	406.80	51%	663.50	83%
18	Thana	3,208	931.50	29%	751.80	23%	1,683.30	52%
19	Amraha	551	156.40	28%	192.10	35%	348.50	63%
20	River Range	9,799	2,424.50	25%	4,752.00	48%	7,176.50	73%
21	Tikri Khera	449	97.00	22%	251.50	56%	348.50	78%
22	Dhaira	1,591	340.30	21%	813.20	51%	1,153.50	72%
23	Sairi	894	179.20	20%	430.40	48%	609.60	68%
24	Koruwa	4,118	807.00	20%	1,632.60	40%	2,439.60	59%
25	Thurau	606	113.30	19%	424.30	70%	537.60	89%
26	Bisoi	1,188	220.20	19%	595.90	50%	816.10	69%
27	Bhugtari	1,602	258.40	16%	893.10	56%	1,151.50	72%
28	Alsi	1,334	160.40	12%	909.60	68%	1,070.00	80%
29	Kalasi	1,770	144.60	8%	505.00	29%	649.60	37%
30	Mangrauli	1,871	89.10	4.8%	783.10	42%	872.20	47%
31	Rani	883	38.70	4.4%	183.30	21%	222.00	25%
32	Vyas Bhood	1,028	37.30	3.6%	86.20	8%	123.50	12%
33	Bantar	1,027	23.70	2.3%	328.20	32%	351.90	34%
34	Chakrata CB	23,079	249.80	1.1%	2,303.20	10%	2,553.00	11%
	Total Area	82,786	17,745	21%	21,496	26%	45,233	55%

5. Model validation

The validation of the LSM results from the fuzzy AHP model was performed through the Area Under the Curve (AUC) approach and relative landslide density (R-index) in terms of prediction probability.

5.1. Receiver operating characteristics (ROC)

The ROC curve represents the true positive rates known as sensitivity versus the false-positive rates (also known as specificity) consistent with the cut-off value or the threshold (Mohammady et al. 2012). In this study, the training sample or the validation data sets were used to validate the LSM model output. The ROC curve was created by plotting the cumulative distribution of the false alarm probability on the y-axis and the variable's detected probability on the x-axis (Qianqian et al. 2017). Next, the LSM was spatially correlated with the landslide validation datasets, and rate curves were produced, and finally, the AUC was estimated (Figure 7).

Figure 7. Validation of landslide susceptibility index through ROC curve and R-index.

The validation result indicates a higher accuracy of the adopted fuzzy AHP model in classifying existing landslides. The AUC value obtained for the Kalsi-Chakrata road corridor, using 30% of the total observed landslides, was 86.52% at a 95% confidence level. This figure depicts that the current model coupled with the use of high resolution latest vintage data has an excellent prediction capability, and the LSM generated has a close resemblance for future landslide scenarios at specified locations. Hence, the fuzzy AHP model adopted here is statistically significant for landslide susceptibility assessment and is endorsed for landslide susceptibility mapping.

5.2. Relative landslide density (R-index)

In addition to the ROC validation above, the LSM output accuracy was tested using the relative landslide density approach (R-index). R-index assesses the relationship among landslide validation datasets and the LSM result (Meena et al. 2019). The R-index can be computed using the following equation- (11) $$R=\left(n_i/N_i\right)/\sum \left(n_i/N_i\right)\times 100$$(11) where $N_i$ is the percentage of landslide susceptible area in each landslide susceptibility class, and $n_i$ is the percentage of landslide in each landslide susceptibility class. The R-index for each vulnerability class is graphically represented in Figure 7 that illustrates the landslide distribution detected in the classes, indicating the consistency of vulnerability levels.

6. Discussion

The study outcomes reveal that the GIS-based fuzzy AHP technique adopted in the present study is a powerful and essential MCDM tool for landslide susceptibility assessment in the Kalsi-Chakrata road corridor. As presented in the introduction section of this paper, the two-dimensional priority matrix (criterion vs. criterion) following Chang’s AHP method is instrumental and the output LSM has closer resemblance with the prevailing landslide sites observed during the field investigation.

Fuzzy AHP deals substantially well with the qualitative conditioning factors (e.g., slope, aspect, soil, geology, LULC) in LSM studies that are difficult to correlate using crisp values in AHP method. Hence, the decision making procedure becomes more comprehensive and rational as both priorities and data are calculated by fuzzy sets (Chen et al. 2011; Mallick et al. 2018). Though a few authors pointed out the lack of matrix consistency in fuzzy AHp studies (Duru et al. 2012), several authors established that the integrated fuzzy sets with AHP give high adaptability in decision making (Mallick et al. 2018; Noorollahi et al. 2018). The fuzzy logic adopted with AHP technique is much easier to comprehend and implement; hence, the weighting can be defined by the experts and data for any measurement scale can be contained within the system (Chen et al. 2011; Feizizadeh et al. 2013; Mijani and Samani 2017).

To better understand the fuzzy AHP model results, the LSI estimated through the fuzzy AHP model was correlated with the landslide conditioning factors adopted in the study. The areas under the very high LSI classes were overlaid on the slope, rainfall, lithology, NDVI, road proximity, and LULC thematic layers. Figure 8 presents the % distribution of the LSI area on these conditioning factors. While correlating the slope classes with the very high LSI areas, it was observed that steeper slope reflects higher landslide susceptibility. The areas with a greater slope than 35° have accounted for about 83% of the very high LSI areas. About 88% of the very high LSI areas are covered by average annual rainfall more than 1,400 mm for rainfall categories. Such higher rainfall areas are located in the study area's central and southern parts between Kasli and Sahiya. The NDVI classes that represent sparse vegetation cover shows more landslide susceptibility. The lithology classes of Chandpur formation and Mandhali-Chakrata formation contains about 92% of the very high landslide susceptible class. Therefore, it is evident that the higher presence of greywacke and unconsolidated rocks are more vulnerable to landslides. About 48% of the area under 200 m proximity from roads shows more landslide susceptibility due to mining activity, cutting slope, and heavy vehicles movement in the study area. Among the LULC classes, the open areas account for more than 53% of the total very high LSI followed by sparse vegetation and forest having 24%, and 18% of their areas fall under very high LSI, respectively. Hence, more attention is needed. It has been observed that landslides in the study area are mainly affected by slope, rainfall, lithology, NDVI, aspect, construction of roads by cutting slopes, and land-use practices. This fact is well supported by the field observations by the author that supports the incidence of past landslides along with the correlation of various conditioning factors in those areas. Overall, the areas located in the southern and south-central part of the road corridor are characterized by a lack of vegetation cover (rocky and barren land class in LULC with NDVI value of less than 0.2) coupled with the steeper slope of >50°, rainfall >1420 mm and greywacke type lithology along the Kalsi-Chakrata highway are more prone to potential landslide occurrences.

Figure 8. Correlation of fuzzy AHP model output and the landslide conditioning factors a) Slope classes, b) Rainfall classes, c) NDVI classes d) Lithology/ Geology Classes, e) Land use classes, f) Road proximity classes.

7. Conclusion

In the preceding two decades, landslides have been considered quite disastrous, coupled with extreme climate events in India's Himalayas. Therefore, to systematically assess the landslide hazard in both the short and long term, this study adopted an integrated GIS-based fuzzy AHP – MCDM approach to understanding the processes that contribute to landslides landslide-prone Lesser Himalayas.

The proposed fuzzy AHP based on Chang’s approach ascertains that this method of GIS-based MCDM can be highly effective for landslide susceptibility assessment. The combined use of validated landslide inventory data and high-resolution satellite images to select and map landslide conditioning factors gratifies the geo-spatial aspect of local variations across the Kalsi-Chakrata road corridor. Regarding choosing the probabilistic LSI model for complex geo-environmental setup, the crisp weights and consistency index method of conventional AHP is not appropriate because of the fuzziness. Following Chang's theory, the fuzzy AHP method is a proven method where the membership values or possibilities are compared during the pairwise comparison. Thus, this method provides more comprehensive, flexible, and substantial results when the decision criteria are qualitative, as in the landslide risk assessment.

The study area's resultant landslide susceptibility zones were identified and mapped into five categories: very low susceptibility, low susceptibility, moderate susceptibility, high susceptibility, and very high susceptibility. The susceptibility map's spatial pattern shows that higher landslide susceptibility areas are mainly concentrated along the south-central and central part of the study area along the Kalsi-Chakrata highway. The analysis presents that about 55% (45.23 km²) of the study area was categorized as a very high and high landslide susceptibility zone (Table 5). About 21% (17.7 km²) was under very high LSI, while 34% (27.5 km²) fell under high LSI categories. At the village level, the result specifies that Nevi (93.4%), Dadhau (92.7%), Chapanu (91.5%), and Kanbua (91.2%) are the villages that are likely to be more susceptible to future landslides.

The present study contributes significantly to providing a useful prediction for the decision-makers and authorities in adopting appropriate approaches to minimize the potential damages and losses from existing and future landslides in the Kalsi-Chakrata road corridor. It has been observed that landslides along the road corridor are predominantly more affected by slope, rainfall, lithology, NDVI, aspect, construction of roads, and land use practices. Due to the uncertainties inherent in the landslide conditioning factors, certain uncertainties will always persist in landslide susceptibility assessment. The subjectivity of evaluation by experts also may have some limitations. Therefore, the fuzzy AHP approach in the future may consider other significant conditioning factors such as spatiotemporal change of rainfall distribution and frequency under the climate change scenarios. The fuzzy AHP model applied in this study shows 86.52% accuracy in predicting landslides evaluated through the ROC curve (at 95% confidence level). Thus, the current model has an excellent prediction capability, and the LSM generated for this study has a close resemblance for future landslide scenarios at specified locations. Hence, the fuzzy AHP method applied herein may be encouraged for landslide susceptibility assessment in the other Lesser Himalayan road corridors.

Acknowledgements

Authors like to thank the anonimous reviewers for their constructive inputs on the manuscript.

Data availability statement

The data that support the findings of this study are available from the corresponding authors, upon reasonable request.

Disclosure statement

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