{'fraction_non_alphanumeric': 0.035962155176540035, 'fraction_numerical': 0.022320787529715846, 'mean_word_length': 5.213740458015267, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 10, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Application of Data-Driven Techniques For Landslide Susceptibility Prediction Along An Earthquake Affected Road Section in Kashmir Himalaya


 The interaction of seismic events with geo-environmental conditions and anthropogenic activities may exacerbate the risk of landslide hazard in a mountainous region. As an example of this, 2005 Kashmir earthquake triggered a large number of shallow to deep slope failures, which was further intensified in following years by human activities notably along road networks, posing a long-term hazard. Hence, this study was planned to evaluate the effectiveness of landslide susceptibility prediction along earthquake affected road-section of Neelum Highway using six different data-driven models. We applied analytical hierarchy process as heuristic approach, weight of evidence and index of entropy as statistical models and multi-layer perceptron, support vector machine and binary logistic regression (BLR) as machine learning models. Initially, 224 landslides locations were marked through field surveys to prepare landslide inventory which was further randomly divided into training (70%) and testing (30%) datasets. Then, 13 landslide causative factors (LCFs) were extracted from geo-spatial database and analysed by measuring collinearity among factors and assessing their contribution in landslide occurrence using different feature selection methods for inclusion in susceptibility modelling. Thereafter, six employed models were trained to produced landslide susceptibility maps of investigated road-section. Finally, the area under receiver operating characteristics (AU-ROC) curve and various statistical measures were applied to validate and compare the performance of modeled landslide susceptibility. The results revealed that no collinearity issue exists among all 13 LCFs, and all six models exhibited satisfying performance in predicting landslide susceptibility of study area. However, BLR model have produced most promising and optimum results as compared to other models with AU-ROC (0.881), Matthew’s correlation coefficient (0.609), Kappa coefficient (0.604), accuracy (0.797) and F-score (0.787). The outcomes of this study can be used as pertinent guide for preventing and managing the landslide disaster risk along Neelum Highway and beyond.


Introduction
Landslide is a kind of complex geohazard in the mountainous terrains around the world that imposes an enormous threat to population, sources of livelihood and infrastructure. It occurs as a result of the interplay of geological events, geo-environmental factors and anthropogenic activities (Aleotti and Chowdhury, 1999;Guzzetti, 2005). Globally, the Himalaya is one of the most favourable and dynamic mountainous terrain for destructive landslide hazard where combine impact of periodic seismicity, monsoon rains and human activities along unstable slopes often aggravate the phenomena of landsliding (Gerrard, 1994;Rusk et al., 2022). The region of Kashmir Himalaya is falling in seismically active zone of western Himalayan belt, had faced devastating seismic events in the past with most recent catastrophic earthquake of 2005 (Gardezi et al., 2021;Sana and Nath, 2017). This earthquake event and aftershocks had triggered more than two thousands of landslides over an affected area of around 7500 square kilometres. The devastation had badly affected the major components of Neelum Highway and resulted in severe damaging and blocking of main route even three months after the main event (Basharat et al., 2014a;Rahman et al., 2014). This mega event and post-earthquake shaking highly deformed the rocks and significantly increased vulnerability of slopes for landsliding which leads to the development and reactivation of landslides upon subsequent interaction with other extrinsic variables (Khattak et al., 2010). Therefore, it is crucial to assess the possible zones of slope failures in this region using modern practices to formulate effective landslide mitigation strategies. For this purpose, the landslide susceptibility mapping is considered as a primary component in landslide hazards and risk assessment in past decades (Huang et al., 2020).
Landslide susceptibility represents the spatial probability of landslides in a particular region based on local geo-environmental conditions and causative factors. Over the years, the number of studies were performed by several investigators e.g., (Achour and Pourghasemi, 2020;Aditian et al., 2018;Chen et al., 2019;Guzzetti, 2005;Hong et al., 2020;Huang et al., 2020;Merghadi et al., 2020;Saha and Saha, 2020;Wang et al., 2021;Wang et al., 2016;Zhang et al., 2018) to assess the landslide susceptibility using different modelling approach, which can be generally categorised into deterministic and data-driven models. In general, the deterministic models calculate the quantitative stability factor of slope in certain area to predict landslides events. These models are usually applied to a local area having simple landslide types and homogamous geologic and geomorphic conditions, however the hazard complexity and extensive data acquisition make it difficult to analyse over large scale in heterogeneous conditions (Ciurleo et al., 2017). In contrast, data-driven models exhibit more reliable results with high accuracy in predicting regional landslide susceptibility and can be categorized into heuristic, statistical and machine learning models (Huang et al., 2020). The heuristic approach can execute landslide susceptibility prediction (LSP) by measuring the weight values for each conditioning factor based on expert opinion derived from the field knowledge using rating techniques, mainly include analytical hierarchy process (AHP) and weighted linear combination (WLC) (Shahabi and Hashim, 2015). On the other hand, statistical approach predicts the probability of landsliding through numerical relationship between observed landslides and causative factors (Erener and Düzgün, 2012). The statistical techniques such as certainty factor , evidential belief function (Ding et al., 2017), frequency ratio (Li et al., 2017;Yalcin et al., 2011), generalized linear models (Goetz et al., 2011;Huang et al., 2020), information value method (Chen et al., 2014;Du et al., 2017), index of entropy (IoE) and weight of evidence (WoE) (Liu and Duan, 2018;Regmi et al., 2014) have been successfully applied in various studies to produce landslide susceptibility maps (LSM) with significant model accuracy. Similarly in last few years, machine learning (ML) models are rapidly evolving in predicting landslide susceptibility due to high precision rate and adequate handling of complex data. The commonly used ML methods in LSP include AdaBoost (Hong et al., 2018), binary logistic regression (Huang et al., 2020;Pradhan and Jebur, 2017), decision tree (Dou et al., 2019), gradient boosting Sahin, 2020), K-nearest neighbors (Pradhan and Jebur, 2017), kernal logistic regression (Aditian et al., 2018;Park et al., 2013), multilayer perceptron (MLP) Zare et al., 2013), naïve bayes Tien Bui et al., 2012), neuro fuzzy (Pradhan, 2013), random forest Dou et al., 2019) and support vector machine (SVM) Pradhan, 2013).
Although, a number of approaches have been proposed for producing LSMs in existing literature, but still no consensus developed among scholars around the globe that which model is more reliable and best suited for predicting landslide susceptible areas in different terrains.
However, researchers recommend that these models should be applied and tested in different geo-environmental conditions, where landslide activities are governed by multiple processes to validate the quality and accuracy of model used (Pourghasemi and Rahmati, 2018).
Moreover, few studies suggest the practical application of different models and their comparative analysis in a similar area should be investigated to obtain a robust model for LSP. Therefore, in light of foregoing consideration, it is indispensable to implement and compare popular data-driven techniques in a given area for better understanding of landslide susceptibility prediction.
As the most disastrous and tectonically active region of Kashmir Himalaya, Muzaffarabad-Neelum trunk road was chosen as example for present study. This road-section experienced a devastating 7.6 Mw earthquake that triggered enormous landslides along natural and cut slopes with subsequent blockage of accessed roads Owen et al., 2008). The main seismic event and its aftershocks developed extensive fissures, cracks, lateral and back scarps in surrounding slopes of this road-section, which resulted in frequent and progressive slope failure specially in extreme weather conditions. Hence, it is essential to understand and recognize the persistent patterns of landslide hazard for mitigating future slope failures in this region. Most of studies have been conducted to identify the occurrence, size, spatial distribution and pattern of landslides associated with 2005 Kashmir earthquake e.g., (Ahmed et al., 2021;Basharat et al., 2014a;Basharat et al., 2014b;Basharat et al., 2016;Jadoon et al., 2015;Kamp et al., 2008;Khan et al., 2013;Khattak et al., 2010;Owen et al., 2008;Saba et al., 2010;Shafique, 2020;Shafique et al., 2016). However, only few studies were focused on landslide susceptibility mapping using conventional quantitative techniques such as Basharat et al. (2016) and Kamp et al. (2008) undertook LSP for earthquake affected areas based on semiquantitative approach (AHP). Recently, Ahmed et al. (2021) demarcated the landslide prone areas along a small portion of Neelum road using WOE method and investigated the geotechnical aspect of specific landslides. Though these studies demonstrate landslide hazard maps which were modulated by limited number of landslide conditioning factors, whereas a comprehensive study on spatial prediction of landslide susceptibility using popular machine learning models is still lacking for this region. We therefore undertook an extensive analysis of machine learning techniques (MLP; SVM; BLR) and their comparison with heuristic (AHP) and traditional statistical methods (IoE; WoE) using a broad range of causative variables to evaluate a robust model for accurate LSP along earthquake affected road-section in Kashmir Himalaya.

Study area
The study area lies in the western part of main landmass of the State of Jammu and Kashmir and is approximately situated between the longitudes 73°26' E to 73°56' E and latitudes 34°21' N to 34°40' N. It is positioned along the main drainage of Neelum River and covers an area of ∼287.97 sq.km (Fig. 1c, d). The Athmuqam city is located at eastern, Nauseri in the central while Muzaffarabad city falls in the western part of the study area connected by a major trunk road. This road is the lifeline for more than 0.3 million inhabitants of Neelum valley. Moreover, it has a significant strategic importance as it is located along Pakistan-India Line of Control.
Topographically, the study area dominantly characterized by hilly terrain with elevation ranges from 632 m to 2857 m above sea level (asl). The terrain in northern and central regions of the study area are relatively higher with altitude exceeding ~1100 m asl as compared to western region which have mean elevation of ~800 m asl (Fig. 1d). The settlements in the study area are widely dispersed along gentle slopes, river and streams terraces and valley floors increasing pressure on the land and environment due to high population density and small land holdings.
The climate of the region is humid, humid tropical and subtropical highland with average annual rainfall of 1400 mm. Mean maximum and minimum temperatures in the area during summer are about 20-34°C and in winter 2-18°C, respectively (Bashir and Hanif, 2018).

Geomorphological setting
The study area being a western part of Lesser and Sub-Himalayas is classified into three geomorphological zones i.e., mountainous terrain, the lower hills of valley slopes and foot slopes of fluvial sediments along the major tributaries. The valley floor in the northern part is deeply incised by Neelum river and its feeding tributaries having a relief of about 2200 m. The rapid flow of the river in this region with an average discharge of nearly 230 m 3 /s has resulted in fluvial incision forming a steeper low valley with slopes exceeding 50°. This landscape reflects a high relief energy resulting a V-shaped valley with high erosion rate. The back slopes above these valley slopes are less steep with average slopes of about 15-30° . Talus or slided-mass, fluvial deposits, terraces, alluvial fans and glaciofluvial deposits are major geomorphic accumulations observed in the study area.

Geological and tectonic settings
Geologically, the area is comprised of sedimentary, metamorphic and igneous rocks which ranges from Precambrian to Miocene. The lithostratigraphic units exposed in the northern part of the study area are Precambrian Salkhala Formation comprising of talc-quarzitic schists, graphitic schist with subordinate marble and Cambrian Jura Granite of gneissic and porphyritic type with doloritic intrusions at certain places. The Panjal Formation consist of greenstone, andesitic flows with bedded tuff, and meta-carbonates with agglomeratic slate, phyllite and mica schist exposed in the central part. The Tertiary fluvio-clastic deposits of Murre Formation containing sandstone, mudstone and siltstones formed as a consequence of Himalayan orogeny (Bossart et al., 1988) occupies the central and southern part of the study area. The Paleocene-Eocene carbonate deposits are unconformably overlies along the cherty dolomitic limestone of Cambrian Abbottabad Formation near Muzaffarabad. The quaternary deposits of fluvial origin are mostly exposed along the lower slopes and Terraces (Hussain and Khan, 1996;Shah, 2009).
The study area lies in the core and limbs of prominent structure i.e., Hazara-Kashmir Syntaxis (HKS), which implicates the region as complex tectonic zone (Calkins et al., 1975;Hussain et al., 2020). Tectonically, the area is sub-divided into three units: a) the youngest Tertiary Molasse deposit in the core of HKS surrounded by the Main Boundary Thrust (MBT) and Pre-molasse sequences, b) the tholeiitic-alkaline volcanic and agglomeratic slates, extending as a linear belt between the MBT and the Panjal Thrust, c) the oldest Precambrian meta-sediments and Cambrian granites exposed on the eastern and western limbs of HKS. The MBT in the study area constitute a suture along which remnants of a Paleozoic ocean floor sediment, the Panjal Volcanics and the agglomeratic slates have been thrusted over the Murree Formation (Baig, 2006;Hussain et al., 2009) (Fig. 1b, c). The western limb of the HKS is more  (Kaneda et al., 2008). This event created many cracks on the hillslope resulted in slope failure of bed rock as well as colluvial deposits along major roads, natural slopes and other man-made slope cuts .

Dataset preparation
The landslides susceptibility modelling in a given area involves the collection of data on landslide occurrence, geospatial database construction on landslide conditioning factors, and probability modelling for landslide susceptible areas based on relationship between landslide events and their triggering factors followed by validation of produced models (Guzzetti, 2005;Regmi et al., 2014). For this purpose, the dataset we used to construct different thematic layers includes ALOS PALSAR DEM with 12.5 m resolution, geological maps (GSP), published reports and research articles, Landsat-8 Images (USGS), Google Earth imageries, rainfall data (TRMM; Pakistan Metrological Department) and detailed field mapping.

Landslide inventory mapping
The preparation of accurate landslide inventory is an important step for investigating the possible landslide hazards and susceptibility mapping (Guzzetti et al., 2012). These inventory maps demonstrate information about the landslide locations, its dimensions and types of mass movements that have triggered in the past. In present study, the landslide distribution has been extracted from previous literatures (Basharat et al., 2014a;Kamp et al., 2008;Owen et al., 2008), historical records of GSP, and visual interpretation of Google Earth images (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020). The extensive field survey was also carried out during March-April 2021 for validation, precise demarcation of landslides area and their characteristics. For this purpose, the landslides dataset recorded systematically during field observation using GPS based Field MOVE Clino (Android App) and finally converted into vector database in GIS for preparation of landslide inventory map. A total of 224 landslide events were recorded during this study and are illustrated in (Fig. 2).
The landslide distribution analysed during fieldwork reveals that geo-environmental and anthropogenic factors are controlling mass movements in most part of the study area, whereas co-seismic landslides are mostly concentrated along the proximity of BBF particularly around the north-western portion of hanging wall slopes of fault trace. Furthermore, the majority of landslides are rock falls/slides and debris slides, while the types of movement are translational, rotational and complex (Fig. 3). The clustered failures were mostly shallow involving bedrock to superficial accumulations, whereas size of the materials ranges from few m 3 to thousands of m 3 . Based on landslides inventory mapping, the study area can broadly be categorized into three sections: a) northern section is a narrow dissected valley characterized by rockslide, rockfall and debris slide of steeply jointed rocks of Precambrian metamorphic and Cambrian granites; b) comparatively broader valley in the central section having less steep slopes (<50 ○ ) triggering mixed landslides mostly associated with road cuts and artificial water channels in Tertiary molasses and Quaternary fluvial terraces; c) in the southern section, rockfalls-mixed landslides developed in Cambrian dolomitic limestone and Paleocene-Eocene marine sediments mostly clustered along the proximity of active BBF. In current study, the compiled landslide inventory was converted as landslides point data based on sampling strategy specified by , whereas equal number of nonlandslides points were randomly selected from stable areas. This dataset was then converted into training (70%) and testing (30%) data using ArcGIS "Subset" tool for modelling landslide susceptibly and its validation.

Landslide causative factors
The selection of landslide causative factors (LCFs) generally rely on the available data, scale and nature of study area as there is no set criteria or universal agreement until now (Nefeslioglu et al., 2008). Our study area lies in the Himalayan region where active tectonics, geo-environmental variables and anthropogenic activities are main contributors of mass movements. Therefore, thirteen parameters consisting of geological factors (lithology and proximity to fault), geomorphic factors (slope, aspect, curvature, elevation, topographic wetness index (TWI), proximity to drainage and lineaments density), environmental factors (Normalized Difference Vegetation Index (NDVI), landuse and precipitation) and anthropogenic factor (proximity to roads) were selected as LCFs and were standardized to a common cell size (12.5 x 12.5 m).

Geological factors
The fragile lithology, sheared and jointed rocks affects the stability of materials which often leads to landslide event along steep slopes (Basharat et al., 2014a;Kamp et al., 2008). Therefore, in present study lithology and distance to fault were taken as input variables for landslide susceptibility modelling.

Lithology
Lithology is classified in terms of rock and soil composition, whereas, its mechanical properties and composition controls the permeability and strength of slope forming materials, thus considered as an important parameter in landslide susceptibility mapping (Ayalew and Yamagishi, 2005;Yalcin et al., 2011). The information on the spatial distribution of lithological formations were identified from the geological maps of GSP (Hussain and Khan, 1996;Hussain et al., 2006;Hussain et al., 2009). Eleven geological formations are exposed, which comprised of lithological units derived from igneous, metamorphic and sedimentary origin such as granite, doleritic intrusions, schist, slates, phyllites, marble, basaltic volcanic flows, dolomitic limestone, cherty nodular limestones, marine shales, sandstone, siltstone and mudstone with quaternary fluvial terrace, alluvial fans and colluvial deposits at certain places (Fig. 4h).

Proximity to faults
The faults are categorized as zone weakness and are responsible for rocks fracturing and jointing which controls the failure mechanism in the form of rock sliding and wedge sliding.
The tectonic stresses in the proximity to faults controls the permeability of terrain and affect the slope stability is considered as the important influencing parameter for slope failures in landslide susceptibility Pradhan and Jebur, 2017). The study area being part of North-eastern Himalayas of Pakistan is recognized as the zone of high seismicity with active thrust faulting often contributing seismic hazard accompanied by lethal landsliding (Gardezi et al., 2021). The most recent seismic event of 2005 had produced thousands of landslides with extensive cracks, fissures and bulging slopes in the region (Basharat et al., 2014a;Kamp et al., 2008;Owen et al., 2008). The field observations recorded in the study area revealed that rocks fracturing was common over a proximity of 500 m from fault, whereas the intense zones of fractures were seen in the range of 0-150 m. Therefore, for understanding the relationship between fault presence and landslide occurrence, five buffer zones were created along the major fault lines in ArcGIS "Euclidian distance" tool and classified accordingly: <150 m, 150-300 m, 300-450 and >500 m ( Fig. 4i).

Geomorphic factors
Geomorphological parameters are considered as major influencing terrain factor in landslide occurrence (Pourghasemi and Rahmati, 2018). The geomorphic parameters we used in this study are briefly described below:

Slope angle
Slope angle plays a vital role in contributing landslide activity and is considered an important causative factor in susceptibility modelling of landslide hazards (Huang et al., 2020;Wang et al., 2020a). The common slope failure generally depends on the degree of cohesion in slope material and inclined angle where gravity exerts a sliding force on slope mass . In our case, the slope gradient ranges from 0-80.27° and were reclassified into seven classes based on landslide distribution analysed during field study. The slope angle categorized with10° interval as: <10°, 10-20°, 20-30°, 30-40°, 40-50°, 50-60°and <60° (Fig. 4a).

Aspect
The aspect of slope is generally described as slope facing direction having significant impact on several processes (vegetation distribution, discontinuities patterns, weathering and moisture retention) which can directly or indirectly effects the landslide phenomena (Shahabi and Hashim, 2015) . The environmental conditions such as snow melt precipitation pattern and sunlight exposure also influence the slope stability with subsequent changes along particular slope direction (Ayalew and Yamagishi, 2005). The slope aspect in this study was classified as N, NE, E, SE, S, SW, W, NW and flat (Fig. 4b).

Curvature
Slope curvature is one of the important geomorphological factors widely used as LCFs in landslide susceptibility mapping (Dou et al., 2015). It represents the concavity, flatness and convexity of terrain surfaces which stimulates the flow velocity and divergence of runoff water and controls the occurrences of landslide (Bera et al., 2019). The curvature values for current study were reclassified as flat (zero curvature), concave (+tive curvature) and convex (-tive curvature) (Fig. 4c).

Elevation
The altitude has a strong relationship with mass movement and is considered as primary risk factor that controls the secondary causative variables such as weather conditions, vegetation, human activities and other conditions (Aditian et al., 2018). Various landslide studies have shown that slope failures are scarce at higher altitude .
In current study, the elevation data obtained from DEM was sub-divided into five classes with an average interval of 500m (Fig. 4d).

TWI
TWI indicates the topographic control on local hydrological conditions, evaluated through soil saturated zones of an area with respect to catchment locations (Wang et al., 2020a).
The degree of water distribution in a soil increases the pore pressure which affect the soil, rock and vegetation condition on upslope resulting in slope failure (Sun et al., 2020). According to Beven and Kirkby (1979), it is defined as: where α is expressed as the flow accumulation in certain point and β is measured as the local degree slope (radian). In present study, the TWI values were calculated through DEM using Arc GIS operations and were reclassified as: <5, 5-9 and >9 (Fig. 4l).

Proximity to drainage
The proximity to drainage is considered as an important landslide causative factor, influencing the slope stability in terms of erosion along the toes of slopes and by increasing the water saturation level in lower part of the slope (Bera et al., 2019). This effect become adverse when water laden drainage borders the steep slopes and sliding zones because the fluvial incision and lateral cutting gradually favours the condition for mass movement (Sun et al., 2020). In the study area, the proximity to drainage was calculated through "Euclidean distance" tool in GIS and was classified into 4 classes with an interval of 100 m: 0-100 m, 100-200 m, 200-300 m and > 300 (Fig. 4k).

Lineaments density
The lineaments pattern in a specific area can be expressed as linear topographic signatures controlled by geomorphic and tectonic processes. These features may affect the stability of slopes and increases the probability of slope failures (Ayalew and Yamagishi, 2005). The lineament patterns were prepared through visual interpretation of Landsat 8 OLI imagery by using PCI Geomatica software. The extracted lineaments were processed in ArcGIS by Euclidean density function to generate the lineaments density map with five classes: 0-0.45 km/km 2 , 0.45-0.9 km/km 2 , 0.9-1.35 km/km 2 , 1.35-1.9 km/km 2 and 1.9-3.5 km/km 2 (Fig. 4g).

NDVI
The NDVI is measured as degree of vegetation distributed in a region and is closely associated with water flow and penetration and soil erosion which influence the landslide activity along slope surfaces (Huang et al., 2020). Therefore, it has been effectively utilized as a crucial parameter in landslide prediction modelling. In current study, the values of NDVI were computed using Landsat 8 OLI multi spectral imagery of 15 th Sep, 2019 (https://earthexplorer.usgs.gov) by following formula: where R is computed as the red spectral band and NIR is near infrared band of electromagnetic spectrum (Tucker, 1979). These values range from -0.15 to 0.56 and were categorized as <0.1, 0.1-0.3 and >0.3 (Fig. 4e).

Landuse
The landuse is frequently used contributing factor in susceptibility modelling for landslide especially in the mountainous terrain where natural process effect the slope stability of a region . The changes in land cover due to landscape modification for infrastructure, deforestation and farmland on steep slope directly influence the hydrological conditions which affect the slope failure mechanism (Wang et al., 2020b). The landuse map of study area was prepared by analysing Landsat 8 OLI images through supervised image classification technique in ArcGIS environment. The output of landuse types were categorized into 5 classes i.e., barren land, forest, grass land, settlements and water body (Fig. 4f).

Precipitation
Various studies suggest that extreme events of precipitation affect the slope stability by triggering new landslides or reactivation of old landslide. Thus, the relationship of landslide occurrence with precipitation had been adopted as a critical LCF in predicting the landslide hazard around the globe (Huang et al., 2020). The effect of climate variations in study area during intense precipitation of monsoon season and winter westerlies increases the infiltration and erosion rate along vulnerable slopes which influence the probability of rapid slope failure (Rahman et al., 2014). In this paper, mean annual precipitation data of study area  were interpolated through inverse distance weighted (IDW) function in ArcGIS to generate the precipitation map. This map was further classified into four classes: <1000 mm, 1000-1200 mm, 1200-1400 mm and >1400 mm (Fig. 4m).

Anthropogenic factors
The anthropogenic activities such as construction of buildings, roads, terracing, stone quarries, etc., disturb the equilibrium state of slope material and increase the vulnerability of landslide hazard along steep slopes (Regmi et al., 2014). In this research, the proximity to roads was selected as the main anthropogenic factor influencing slope failures in the study area.

Proximity to roads
Roadcut being a site of anthropogenically disturbed slopes, can act as a catalyst for landslide occurrence in a mountainous region (Ayalew and Yamagishi, 2005). The excavation and obsolete blasting techniques for construction of roads alter the natural stability of slope and accelerate the probability of landslide incidences (Das et al., 2010). In present study, the landslide failure zones due to slope cutting were observed at a distance of approximately 100-200m to major road network. Therefore, we have generated proximity to road map comprising of four buffer zones: 0-100m, 100-200m. 200-300 and > 300m (Fig. 4j).

Methodology
The methodological approach adopted for comparative assessment of different landslide prediction models is illustrated in figure 5 and categorically explained below.

Landslide predictors analysis
In this research, factor analysis of LCFs were adopted to identify the relevant predictors by measuring covariation among multiple variables and assessing the relative importance of causative factors in landslide occurrence.

Collinearity analysis
In landslide studies, collinearity among predictors may occur due to linear association of two or more independent variables which reduces the predictive performance of susceptibility models. Therefore, it is important to diagnose and eliminate the collinearity problem between input variables before data mining . To address this issue, variance inflation factor (VIF) and Pearson correlation were performed as statistical tests in this research.
VIF helps to quantify multicollinearity arise due to linear relationship between one variable and another variable in a model (Di Napoli et al., 2020). Mathematically, it is computed as: here, 2 is the coefficient of determination calculated by regression analysis of i th predictors to other predictors. If VIF = 1, it indicates no correlation between one LCF to the remaining ones, whereas VIF >5 or 10 reflects higher multicollinearity (James et al., 2013).
On the other hand, the Pearson correlation is a statistical test that computes the extent of linear relationship between two predictors and their association with each other . The resultant linear association expresses how strong or weak relation exists between them. The coefficient of Pearson correlation (r) between two LCFs "a" and "b" is calculated as: here, n refers to the size of data points, ai and bi denote the data points of a and b indexed with i, whereas ̅, ̅ and , are the means of observations and standard deviations of data points of variables a and b. The outcomes of the "r" range between -1 (negative correlation) and +1 (positive correlation), while zero value indicate no correlation among two variables. However, according to Nettleton (2014) the coefficient values of factors >0.7 represent high collinearity which can affect the model performance and should be eliminated before data modeling.

Factors contribution analysis
The landslides and their contributing factors vary from area to area depending on local geo-environmental conditions (Hong et al., 2020). Therefore, evaluating the role of different LCFs in terms of their importance for landslide occurrence will help to enhance the predictive capability of susceptibility models. For this purpose, the gain ratio and relief-F methods were applied in this study to estimate the relative importance of each conditioning factor in triggering here, IG (A) and E(A) are the information gain and entropy values with respect to an attribute "A". The function E( ) = ∑ − ( ) 2 ( ); where ( ) indicates probability score of i th class by contributing to overall scores of attribute "A" classes. The IG function of a given attribute "A" with respect to class "X" is calculated in Eq. (6).
where, H(X) represents the class "X" entropy and ( ) is an average entropy of attribute "A" with respect to class "X" and are computed using Eq. (7) and Eq. (8), respectively: in Eq. (7), ( ) is probability score of j th class corresponding to overall scores. Whereas in Eq. (8), m represents total number of class partitions, is the probability score of j th partition (Quinlan, 1993).
Relief-F is primarily used as data pre-processing technique that computes the quality and relevance of conditioning factor in a multiclass dataset based on correlation between predictors and their classes (Urbanowicz et al., 2018). Firstly, the relief-F randomly select an instance of a specific attribute, then look for the k nearest neighbours of same and different classes of selected instance termed as hits and misses, respectively. Finally, the rank of each predictor is evaluated by averaging the probability weights of all hits and misses (Kononenko, 1994). The computed scores of attributes in a given data are generally range between -1 (least significant) and +1 (most significant) (Kira and Rendell, 1992), however in landslide studies, if the factor contribution is equal to 0 or < 0, then it is considered as redundant factor and must be ignored in susceptibility modeling (Ali et al., 2021).

Landslide susceptibility models
In this paper, we have selected six popular models comprising of AHP, IoE, WoE, ANN, SVM and LR for predicting the landslide susceptibility of study area. The working mechanisms of these models are briefly described in the following subsections.

Heuristic model
AHP is the only heuristic approach we employed for the production of landslide susceptibility map. It is a well-known semi-qualitative method that involves expert based judgements to compute the criteria weights of different contributing variables through relative pair-wise comparison (Franek and Kresta, 2014;Saaty, 1990). The method has been successfully applied in various studies related to susceptibility mapping for landslides hazard assessments (Bera et al., 2019;Huang et al., 2020;Park et al., 2013;Shahabi et al., 2014;Yalcin et al., 2011). In our case, the practical framework adopted for computing weights of factor and their classes through AHP comprising of four stages i) hierarchic arrangement of causative factors, ii) pairwise comparative judgement (expert opinion) based on relative importance of factors using rating scale (1 to 9 and its reciprocal), iii) weight estimation of each parameter and its classes by computing the principal eigenvector, and iv) consistency checking to eliminate redundant judgments through consistency ratio (CR) which should be < 0.10 (Saaty, 1990) The CR was calculated as the ratio of consistency index (CI) and random inconsistency index (RI). Here CI is defined as: where, λ max and indicates the largest eigenvalue and total number of factors/ classes in a given matrix, respectively. Whereas, RI values were adopted from simulated average CI values of pairwise matrix specified by Saaty (1990) and Franek and Kresta (2014). Finally, the landslide susceptibility index (LSI) was produced in accordance with weighted summation procedure using following equation: here, denotes the weight value of LCF ; is j t h class weight of factor and is total number of conditioning factors.

Statistical models
The statistical models can assist in predicting the susceptibility of landslide through various statistical operations by interpreting a relationship of landslides distribution with their causative variables in a given area (Huang et al., 2020). To comprehend the statistical approach for landslide susceptibility assessment in the study area, we utilized WoE and IoE models.

Weight of evidence
The WoE is well established and reliable geo-statistical technique evolved from information and Bayesian theory and is considered as an effective tool for assessment of landslide hazard (Saha and Saha, 2020). The core concept of WoE in landslide susceptibility modeling is to predict likelihood of landslide occurrence under same conditions by establishing a relationship between dependent (LCFs) and independent variables (presence and absence of landslide events (Ding et al., 2017). To achieve this target, firstly we calculated positive ( + ) and negative ( − ) weight scores for each class of LCFs using following equations: here, is the probability; and ̅̅̅̅ indicates the presences and absence of respective causative factor; is the presence of landslide events and ̅̅̅ is the absence of landslide events. Afterwards, the weight scores of + were subtracted from − to find the weight contrast ( ) of each class of LCF, which represents the spatial relationship of landslide events with classes of different causative factors. Finally, the LSI was calculated by the summation procedure using Eq. (13).
where, is weight contrast score for th class of factor and representing total LCFs (Regmi et al., 2014).

Index of entropy
The index of entropy is bivariate statistical approach, successfully utilized by various investigators in the assessment of landslide susceptibility by measuring the extent of various causative factors (Wang et al., 2016). In landslide susceptibility modelling, IoE helps to estimate the weight values of each LCF influencing the occurrence of landslide using following equations: = log 2 (17) = − , = (0, 1), = 1, . . . , Here in Eq. (18), is the information coefficient, whereas is the resultant weight values of each contributing factor in Eq. (19) Liu and Duan, 2018). Finally, the using IoE model, was prepared by following equation: where, is total number of LCFs used in mapping; the greatest number of classes denoted by ; represent the number of classes within a specific LCF; is the secondary classified value of each class, whereas the was obtained from Eq. (19).

Machine learning models
The Machine learning methods are well-understood statistical algorithms, efficiently explored in multiple disciplines to solve the complexity of big-data for making accurate predictions in contrast to traditional modeling techniques . In recent years, ML techniques have been extensively used to predict the susceptibility of landslides in various geo-environmental conditions with a significant level of predictive accuracy. In present case study, three popular ML approaches i.e. MLP, SVM and BLR were adopted to assess the landslide susceptible zones.

Multi-layer perceptron
The MLP is well known class of feed-forward neural network that uses a supervised learning algorithm to identify complex non-linear associations between predicting variables and response variables (Taud and Mas, 2018). The typical structure of MLP model consist of input/ visible layer, hidden layer and output layer interconnected with each other through a set of neurons. The neurons present in theses layer are connected by corresponding weights which are trained to generate a desired output using activation function (Murtagh, 1991). In this study, we have implemented a single hidden layer MLP neural network to recognize the likelihood of landslide occurrences with respect to input LCFs. The schematic approach of this hierarchical structure is illustrated in figure 6. The learning parameters (epochs = 821, learning rate = 0.3) and activation function (sigmoidal) were calibrated according to (Tien Bui et al., 2016).

Support vector machine
SVM is widely accepted as a robust predictive model in machine learning algorithms primarily used to solve classification problems. The core concept of SVM is to maximize the margin between classes and to segregate data points into different classes by establishing a suitable decision boundary or optimal hyperplane (Wu et al., 2008). In studies related to landslide susceptibility, it has been observed that the predictive performance of SVM model is better than other conventional techniques, however, the accuracy is associated with applied kernel function and other tuning parameters (Zhang et al., 2018). The radial base function (RBF) is proven to be an effective kernel function in LSP; therefore, RBF kernel was used in this study and can be expressed as: here, is kernel function; and are input vectors of landslide training data and is kernel parameter whose values ranges between 0 and 1. The optimum values of tuning or hyperparameters such as C-regularization, kernel coefficient (γ) and loss function (ɛ) were adopted from . Afterwards, these parameters were employed in SVM algorithm using training dataset of normalized LCFs to model the susceptibility of landslides in this study.

Binary logistic regression
BLR is considered as one of most appropriate predictive ML algorithm used for classification of two-class dataset. It computes the probability of any event by establishing a relation between one target binary-variable and several predictors or independent variables (Lever et al., 2016). In this study, we used the presence (1) or absence (0) of landslides as target or dependent variable ( ) and the LCFS as independent variable to map the probability ( ) of landslide occurrence under the given conditions (Huang et al., 2020). Mathematically, the BLR model can be formulated as: = + 1 1 + 2 2 + ⋯ + where, is the intercept of model; 1 , 2 , … are the regression coefficients and 1 , 2 , .

Model validation and comparison
The evaluation of landslide prediction models is a core part of landslide susceptibility mapping for achieving an effective model performance and to ensure the scientific significance of models used. The validation of these models is commonly assessed through several evaluation metrics like area under receiver operating characteristics (AU-ROC) curve, crossvalidation, confusion matrix, Root Mean Squared Error etc (Frattini et al., 2010). In this study, we use AU-ROC curve, Matthew's correlation coefficient (MCC), Cohen's Kappa (K) coefficient, accuracy (ACC), F-score and confusion matrix fourfold-plot to evaluate model's accuracy and to discriminate among resultant models based on their performance scores. The AU-ROC curve was constructed by plotting sensitivity on horizontal axis and 1-specificity on vertical axis for both training and validation data, whereas the quantitative estimation of AU-ROC values was used to assess model quality; < 0.6 indicates poor results, > 0.8 represents very good results (Huang et al., 2020;Pourghasemi and Rahmati, 2018).
here P and N represents total number of landslides and non-landslide counts, respectively; TP (true positive) and TN (true negative) are correctly classified pixel counts; and FP (false positive) and FN (false negative) are incorrectly classified pixel counts (Achour and Pourghasemi, 2020;Chen et al., 2021) . The performance values for LSP models evaluated by F-Score and accuracy varies from 0 to 1 (values close to 1 indicates as high predictive models and vice versa), whereas the maximum and minimum value of MCC range between +1 (best) and -1 (worst), respectively (Wang et al., 2021).
The Cohen's kappa coefficient measures the degree of agreement between observed (Pobs) and expected (Pexp) outcomes using following expression: where, Pobs and Pexp are computed by following equations: here, T is the total pixels used in training; n is the percentage of correctly classified pixels (Ali et al., 2021). The values range between 0 and 1. The values close to 0 represent a random agreement, 0.21-0.40 indicates fair agreement, 0.41-0.60 imply a moderate agreement, 0.61-0.80 shows substantial agreement, whereas values >0.80 demonstrate a strong level of agreement (Cohen, 1960).

Collinearity analysis
In case of landslide susceptibility modeling, the assessment of collinearity is important to diagnose unnecessary variables which may lead to decrease the predictive performance of resultant models. For this purpose, Pearson's correlation and VIF tests were performed to assess the degree of linear association among thirteen LCFs. The results of Pearson's corelation (all <0.7) indicates that all conditioning factors have low degree of association between them (Fig. 7).Similarly, the highest value of VIF is quite less than the threshold limit i.e., 5, (Wang et al., 2020b) indicating no collinearity issue among all LCFs (Table 1). Therefore, the present analysis suggest that all thirteen predictors are suitable for susceptibility modelling in this study.

Factor contribution analysis
Evaluating relative contribution of causative factors on occurrence of landslide is an important task because it helps to detect most significant predictors for efficient landslide susceptibility modeling. In present work, GR and Relief-F methods were implemented to perform this task. The calculated GR values for all thirteen LCFs shows that lithology (0.124), landuse (0.110) and slope angle (0.084) are the dominant contributory factors, while lowest AM scores of TWI (0.009) and curvature (0.007) are within critical limit, therefore all factors were taken into consideration for landslide susceptibility prediction (Fig. 8a). Likewise, the obtained importance using Relief-F indicate that among thirteen LCFs slope angle (0.0877), landuse (0.0837), precipitation (0.0753) and lithology (0.0456) are the top most important factors, whereas TWI (0.0153), aspect (0.0130) and curvature (0.0055) are the least important factors in landslide probability modelling (Fig. 8b).

Landslide susceptibility maps
The spatial prediction of landslide in the study area was explored by three different approaches i.e., heuristic (AHP model), statistical (WoE and IoE models) and machine learning (MLP, SVM and BLR models) to produce LSMs using thirteen landslide conditioning factors.

LSP using heuristic model
The class wise weights of each LCF in AHP analysis were computed through pairwise relative comparison matrix using knowledge driven approach based on extensive field studies.
The result of this comparison matrix represents a reasonable consistency with values <0.10 for all instances, therefore the calculated class weights of all thirteen parameters seem to be more realistic (Table 2). Similarly, the factor weights of all thirteen LCFs were evaluated by adopting the same technique. The CR for this analysis is 0.07 showing an acceptable level of consistency which implies that the calculated factor weights are reliable for predicting heuristic landslide susceptibility model. The final result of class weights and factor weights of 13 LCFs were used in Eq. (10) to produce LSIAHP for the target area. The predicted LSI was reclassified into five classes using natural breaks classifier in GIS to produce LSM representing very low (22.2%), low (26.3%), moderate (26.7%), high (17.6%) and very high (7.2%) susceptible zones (Fig. 9,   11). On the other hand, quantitative analysis of actual landslides in each susceptibility class indicates that maximum percentage of landslide pixels lies in high to very high susceptible zones (Fig. 11).

LSP using statistical models
The WoE and IoE models were used as statistical approaches to compute class weights of each thematic layer by analysing a spatial relationship between landslide training dataset as dependent variable and all LCFs as independent variables (Table 3). The above analysis suggests that lithology, landuse, elevation, proximity to faults, lineament density and slope are most significant conditioning factors in triggering landslides, whereas the Muzaffarabad Formation, barren lands, moderate elevation, areas within 500m of faults, high lineament density and relatively higher slopes (40°-60°) are major factor classes influencing susceptibility of landslide occurrence in the study area.
The resultant LSIs for both statistical models were constructed by applying prediction scores in Eq. (13) and Eq. (20) using summation procedure in ArcGIS environment. These LSIs were than reclassified into five susceptible classes i.e., very low (VL), low (L), moderate (M), high (H) and very high (VH) using widely used Jenks classification scheme (Chen et al., 2021;Merghadi et al., 2020) to produce final LSMs of the study area (Fig. 10) (Fig. 11). Furthermore, the distribution of existing landslide pixels in each susceptibility classes depicts that most of landslide counts can be observed in H to VH susceptible zones of both statistical LSMs, representing a reliable prediction result (Fig. 11).

LSP using machine learning models
In this paper, three standalone machine learning algorithms (MLP, SVM and BLR) were implemented in Weka 3.8.4 software to compute prediction scores of training dataset, which were than processed in GIS to build pixel-based LSIs of entire study area. The resultant LSI of each model was reclassified into very low, low, medium, high, and very high susceptible zones based on Jenks classifier to generate the map of landslide susceptible areas (Fig. 12).
According to all LSMs, the higher susceptibility can be observed along Neelum highway, proximity to main faults and areas bordering Neelum river and its feeding tributes especially

Model performance and comparison
The performance of each model employed in this study were assessed by examining the spatial concordance of actual data of landslide and non-landslide with the prediction results of each model using application of several validation methods i.e., AU-ROC, ACC, F-score, MCC, Kappa and confusion matrices (Fig. 13, 14). The training dataset was utilized to compute the performance of each model in terms of their success rate, whereas the predictive ability of derived LSP models was evaluated by using testing dataset of observed landslides and nonlandslides.  (Table 4).
On the other hand, the prediction scores obtained through AU-ROC for all six models  (Table 4). Furthermore, the contingency plots of all six models illustrating percentages of correctly and incorrectly classified landslide and non-landslide picture elements validate the reliability of employed models in predicting spatial distribution of landslide in study area (Fig. 14). The overall performance scores comparison of all six models indicates that BLR, MLP and SVM (machine learning models) had significantly higher learning and prediction capabilities than statistical (IoE and WoE) and AHP models, respectively.

Discussion
The spatial prediction of landslide susceptibility at various scales provides a foundation for effective landuse planning and hazard management in a complex mountainous landscape (Aditian et al., 2018;Erener and Düzgün, 2012). Owing to its unique topographic configuration, tectonic fragility and harsh climatic conditions, the Kashmir Himalaya is considered as a place sensitive for landslide disasters. Being a part of this region, the studied area has faced an enhanced rate of landsliding in recent decades due to mega earthquake of 2005 and rapid growth in socio-economic activates, which exacerbate the risk of landslide hazard (Rahman et al., 2014;Sana and Nath, 2017). However, the exploration of LSP models using current practices and their comparison with traditional techniques to evaluate an optimal model for landslide susceptibility assessment is still lacking in this region. To address this issue, we employed and compared widely used data-driven techniques such as machine learning (MLP, SVM and BLR), statistical (WoE and IoE) and heuristic (AHP) to generate LSMs of earthquake affected road-section, a part of western Himalaya in Kashmir.
In general, the LSP modelling is mainly dependent on selection of factors contributing landslides, as it has a substantial impact on predictive performance of models. The criteria for the selection of LCFs is still debatable because the variable that influence landslides may varies under different geo-environmental conditions. However, it is important to employ variable selection techniques in order to maximise the performance of LSP model by removing noisy factor before training phase Reichenbach et al., 2018). To ensure this, Pearson correlation and VIF tests were performed which demonstrates that degree of liner association among all the thirteen LCFs lies within the threshold limits specified by Dou et al.
(2020) ( Fig. 7; Table 1). Additionally, the selected LCFs were analysed using GR and Relief-F methods to find the contribution of each factor in aiding to initiate landslides. Results revealed that the slope angle, lithology, landuse and precipitation were the most significant factors while curvature, TWI and aspect were the least important factors influencing landslides in study area (Fig. 8). Our findings are an agreement with the previous studies (Basharat et al., 2014a;Basharat et al., 2016;Kamp et al., 2008;Rahman et al., 2014), where these factors played a dominant role in occurrence of landslides. Similarly, the distribution of landslides and their causes inferred from field observations also complement our findings; for instance, fragile lithology on terrain slopes aggravated the phenomena of landsliding when it rains and the rapid and unplanned construction activities also triggered shallow landslides. On the basis of these analysis, 13 LCFs i.e., lithology, proximity to faults, slope angle, aspect, curvature, elevation, TWI, proximity to drainage, lineaments density, NDVI, landuse, precipitation and proximity to roads were selected in this article to predict landslides susceptibility.
Numerous studies were performed to generate regional scale LSMs and proposed suitable models for the assessment of susceptibilities. Recently a review highlighted that more than 500 research articles based on 70 different models were published in 2005-2016 , however it is pertinent to mention here that few studies (Ahmed et al., 2021;Basharat et al., 2016;Kamp et al., 2008)   addition to the pre-existing landslides data along targeted road-section, our LSMs predicted the landslide potential zones with a high rate of predictive accuracy as compared to previous studies. The generated LSMs of all six models are well harmonized with one another, although the ratio of LSZs with respect to area percentage distribution is different in each model (Fig.   11). Generally, a uniform pattern of high and very high LSZs in all models observed along an earthquake affected road section of Neelum Highway, however it was noticed in field investigations that landsliding is the most common phenomena with potential hazards along slope angle ranging from 40° to 60° and mainly observed in the fragile lithologies especially along faults, roads, and drainages. On the other side, the low and very low LSZs are mainly dispersed in the areas with slope angle less than 20°, far from roads and drainage networks, away from faults and less affected by anthropogenic activities. In general, the machine learning models i.e., BLR, MLP and SVM predicted well with respect to observed landslide percentages of 75.0%, 73.04% and 72.79% for high and very high LSZs in comparison with statistical models (WoE= 70.11%, IoE= 68.54%) and AHP model (69.06%) of area under study (Fig. 11).
Furthermore, different accuracy assessment methods i.e., AU-ROC, ACC, F-score, MCC,  (Pourghasemi and Rahmati, 2018). As a result of these analysis, it can be deduced from the model's comparison in section 5.4 that all the six models used in this study have shown a considerable predictive capability to produce LSMs, although the machine learning models in general have higher rate of predictive accuracy when compared with other models. In addition to this, the machine learning algorithms have many benefits in the model processing and relatively simple to interpret without bulk of background information (Achour and Pourghasemi, 2020;Wang et al., 2021). Whereas, the other data-driven techniques used in this study need several datasets to train and to test each model. Unfortunately, not a single machine learning technique model was used for landslide susceptibility mapping in the seismically active zone of studied region. Only Kamp et al. (2008) produced a regional level LSM of 2005 earthquake affected areas with a prediction accuracy of 67% by using AHP technique, whereas our findings revealed 74% accuracy of landslide prediction by AHP model which has a considerably better degree-of-fit to the assessing data. Overall, BLR method outperformed the other models and exhibits a highest capability in terms of landslide prediction and LSZs distribution, followed by MLP, SVM, IoE, WoE and AHP model ( Fig. 13; Table 4).
Being an important portion of the discussion section, the limitations of the techniques and dataset employed in this study can be summarised as; i) currently, there are no standards and code of practices established to map and evaluate the susceptible areas for landslide occurrence.
Therefore, it is always a challenge in selection and adoption of an appropriate technique/ method for the evaluation of landslide susceptibility, which resulted in an uncertain condition on account of its reliability and credibility of outcomes obtained by the adopted method. ii) in this study, we were focused on earthquake affected road-section to generate LSMs and to compare their predictive performances. For this, landslide inventory was prepared by compiling field-data and remote sensing data. However, due to limitations of historical images, we were not able to exactly identify the landslides before 2005 Kashmir earthquake particularly of smaller size. iii) another limitation of this study is that it doesn't propose the causative mechanism of slope failure due to a complex relationship with triggering factors, thus it can only be employed as an initial step in assessment of landsliding phenomena. iv) the analysis based on AHP model is one of the limitations, as it strongly relies on opinion of geologist/expert and it refers to their subjectivity/ field knowledge in a specific area. v) an important problem can be seen that we didn't classify landslides by their type of movement in order to generate LSMs and considered each type of landslide occurred by following same processes, whereas the causative factors and triggering factors may vary for different mode of movement. vi) finally, the last and most important point is that the results obtained for the assumption of landslide hazard zones are only viable for similar geo-environmental conditions.
Therefore, the predictions based on present geo-environmental conditions and/or the mechanisms that controls the slope failure are not a guarantee for the future, if these conditions change.
To improve the presented methodology of landslide susceptibility zonation, above mentioned limitations (i) to (v) must be considered for developing new LSP models. Recently, novel machine learning and deep learning approaches as well as their ensembles explored in landslide susceptibility mapping with significant prediction accuracies. So, employing these techniques and comparing with other statistical and multicriteria decision making methods can help to find an optimal model for LSP in studied region of Kashmir Himalayas.

Conclusion
The landslide disasters in earthquake prone region of Kashmir Himalayas are mainly controlled by combined interaction of multiple geo-environmental and anthropogenic drivers which badly affected the human population in terms of infrastructure deterioration and socioeconomic development. Hence, it is indispensable to assess the critical vulnerable zone of slope failures in this region for emergency management /risk analysis and minimizing oncoming landslide hazards. The present study was focused to analyse landslide susceptibility along a road-section of Neelum highway where natural slopes were severely disturbed by postearthquake activities. We developed six different LSP models by applying heuristic (AHP), statistical (WoE and IoE) and machine learning (BLR, MLP and SVM) techniques based on site-specific causative factors with detailed landslide inventory and compared their prediction performance to evaluate robustness of employed models. Our analysis revealed that all the six models show promising results, however machine learning models exhibited remarkable LSP in comparison with statistical and heuristic models. The BLR reflected highest prediction accuracy in a seismogenic zone of Kashmir Himalaya followed by MLP, SVM, WoE, IoE models, while AHP was found to be least effective method in demarcation of landslide susceptible areas along targeted road-section. The performance scores derived through statistical metrics provides a valuable indication that the employed models performed well in predicting spatial representation of landslide susceptibility, however, the machine learning portrays highest predictive performance, comparatively. The majority of predicted landslides distribution pattern was observed in western and central part of study area, particularly in the proximity of road/faults and deeply dissected valley slopes which replicates the actual situation of slope failures evident from field observation. The analysis of factor contribution in landsliding and field observations suggested that unstable and fragile lithology on steep slopes were dominant contributors in landsliding mainly controlled by climatic and anthropogenic factors. In conclusion, we presumed that the established methods have a reliable capability to predict landslide potential zones, which could assist the local government and decision makers in mitigation and prevention of landslide disasters along the highway road-corridors.