Landslide susceptibility mapping of Phewa Watershed, Kaski, Nepal

1. Introduction

Landslides are natural phenomena that result in the falling and lateral dislocation of slope materials made of natural rock layers, soil, synthetic materials fill, or mixes of these elements, due to gravity and water [1, 2]. They are triggered by a couple of factors, namely dynamic triggering mechanisms (excessive rainfall and local soil moisture levels) and persisting geological features (steep slopes, weak rock layers, and tectonic activity) [3]. These factors increase the likelihood of landslides occurring in certain areas. Additionally, anthropogenic activities, such as deforestation, mining, and construction can also contribute to the occurrence of landslides by altering the natural landscape and destabilizing the soil. Overall, landslides are the result of a combination of natural and human-induced factors that can lead to the sudden and rapid movement of rock, soil, and debris down the slope. After accounting for landslides brought on by seismic events, 4,862 different landslide incidents were reported globally between 2004 and 2016, which resulted in 55,997 fatalities and economic losses estimated to be in the range of USD 20 billion annually, and the costs of damage may exceed all other costs associated with multi-hazard disasters [4]. Highly developed nations account for 5 % of all deaths caused by natural disasters globally. The remaining 95 % of overall deaths occur in middle and developing countries [5]. Infrastructure is severely damaged by landslides, and during the period from 1950 to 2009, China experienced the highest rate of deadly landslides, followed by Indonesia, India, the Philippines, Japan, Pakistan, and Nepal [6].

Nepal is prone to landslides due to its steep terrain, and unstable geological structure, in addition to high-intensity rainfall during the rainy season, forest destruction, vegetation service charge loads, and frequent earthquakes [7]. Geologically, Nepal being part of the Himalayan regions, comprises soft, fragile rocks, and unstable geological landscapes resulted due to the collision of the Indian and Eurasian plates. These areas are also among those with a high risk of earthquakes because they are tectonically active. As a developing country, Nepal undergoes a variety of construction and development operations. However, due to unplanned and unscientific development, landslides and other natural disasters have occurred. In Nepal Landslides are among the most common catastrophic events, claiming almost 200 lives each year and causing significant damage to infrastructure. Landslides killed 5,813 people in Nepal between 1971 and 2023, 941 people went missing, 2,708 people were injured, 22,319 dwellings were damaged or destroyed, 586,673 families were affected, and 5,654 events occurred [8]. The Table 1 provides additional comparison data, showing that landslides are among Nepal’s severe natural disasters. Landslides in the Kaski district, which includes the Phewa watershed, killed 234 people and left 35 persons missing between 1974 and 2013. Moreover, 152 dwellings were damaged, 5,791 families were affected, and 80 people were injured [9].

Mapping the current landslides is crucial for understanding the connections between the spread of landslides and their causal factors. The probability that landslides will occur in a specific locality is referred to as landslide susceptibility [10]. It is also an effective instrument that assists land planners in developing regions near areas prone to slope instability [11]. The distribution of potentially unstable slopes based on a specific set of Geo-environmental variables is shown on a map of the susceptibility to landslides [12], which divides the land surface into zones according to susceptibility levels [13]. In the form of landslide susceptibility mapping, a significant amount of effort has been made globally to evaluate the landslide-sensitive zones. Among these, Van Westen [14] has thoroughly documented the use of GIS for the zonation of landslide hazards. Geographical Information Systems (GIS) is a tool for data storage, analysis, modeling, and mapping [15], is valuable for analyzing landslide risk, and can play a crucial role in managing and mitigating the impact of landslides on infrastructure and communities. To determine the possibility of a future landslide occurring under specific conditions, landslide susceptibility mapping often evaluates the previous landslide occurrence with numerous triggering elements, including geological, land use/ cover, topographical, and hydrological [16].

In the early 1990s, the manual mapping of landslide risk in Nepal involved field surveys and the collection of data on slope stability, geological conditions, land use, and other relevant factors. This information was then input into a computer system to create maps that identified areas at high risk of landslides. This manual mapping process was time-consuming. Such a mapping of landslide susceptibility started along road corridors and gradually expanded to watersheds and specific areas. Wagner, et. al. [17] designed the computer program named SHIVA to map the soil and rock hazards in Nepal [18]. However, with advancements in technology, such as remote sensing and geographic information systems (GIS), the process of mapping landslide risk has become more efficient and accurate. Today, Nepal continues to use a combination of manual field surveys and computer-based mapping techniques to assess and manage landslide risk.

Different approaches such as heuristic (traditional approach and based on expert knowledge), deterministic (based on data from laboratory tests and field surveys), statistical (based on the mathematical analysis of the field data and subject-matter knowledge), machine learning (based on the correlation between the training datasets previous landslides and the triggering parameters that forecast the high potential area to landslide), are used for LSM. In recent years, several geoscientists and geotechnical professionals have developed and are implementing hybrids of semi-quantitative/semi-qualitative decision-making, statistical/geostatistical approaches, inventory-based probabilistic methods, and so on [19-23]. Several statistical methods have been used widely, including the weights of evidence [24], frequency ratio (FR) [25], and evidential belief function [26]. Due to its reliable calculation and precise results, the statistical approach is becoming more and more popular. An evaluation of several statistical approaches revealed that the FR model outperformed the others [27]. Bourename, et. al. [28] also contrasted five statistical techniques, namely WoE, AHP, Weighting Factor (WF), Logistic Regression (LR), and FR, and concluded that FR had a greater prediction rate of 86.59 percent. For instance, Regmi, et. al. [29] evaluated different statistical techniques, such as the statistical index (SI) and weights of evidence (WoE) approaches and discovered that the frequency ratio (FR) fared better.

Mapping landslide susceptibility typically involves utilizing statistical hazard, landslide distribution, landslide frequency, qualitative hazard, deterministic hazard, and other approaches as standard features [30]. The frequency model is often considered more effective than other approaches for landslide susceptibility mapping due to its ability to directly quantify the frequency of landslides in a given area. This model provides a clear and straightforward assessment of the likelihood of landslides occurring in specific locations, which can be valuable for risk assessment and mitigation efforts. Additionally, the frequency model is often more easily interpretable and can provide a more tangible understanding of landslide susceptibility compared to other approaches. The frequency ratio (FR) model from the statistical method for landslide susceptibility mapping was selected for this study due to its distinct advantages over various other models. To generate Landslide Susceptibility Mapping in the study region and calculate the spatial relationship between two variables i.e., landslide events and causative factors, the Frequency Ratio model is an effective advanced method [31]. The FR model is regarded as a reliable and experimental approach that may be utilized to determine the relationship between landslide events and causal elements and to provide useful Landslide Susceptibility Mapping for the research area [32]. The FR model has a significant advantage in that it can determine the rank of the causative factors concerning a landslide occurrence as well as whether or not a specific range of causative factor values will be dangerous in the event of landslides [33].

Numerous studies have been conducted at international, national, and local levels, yet none have focused on landslide susceptibility modeling to mitigate landslides in the Phewa watershed, despite existing threats from landslides. The primary goal of our research is to assess the relationship between landslide events (dependent variable) and causative factors (independent variables) to create an accurate, updated landslide susceptibility mapping of the study area using the FR model. The study aims to highlight the different zones of landslides in the Phewa Watershed of Nepal.

2. Materials and methods

2.1. Study area

The Phewa watershed extends from 83°48'1.93" east to 83°58'13.18" east, as well as from 28°11'37.83" north to 28°17'27.30" north, with an elevation ranging from 793 meters 2508.81 meters (Fig. 1). In Nepal's mid-hill zone, the watershed is found. Following a Cabinet meeting, the Gandaki province government declared in 2022 that the 123 square kilometers of the Pokhara Lake cluster listed on the Ramsar site will become the Phewa Lake protected watershed area. The watershed area is made up of major parts of Pokhara Metropolitan City and some parts of Annapurna Rural Municipality. The watershed region is comprised of 44.95 percent forest, 38.52 percent agricultural land, 10.76 percent built-up area, 4.29 percent water bodies and barren land, and 1.48 percent flooded area or barren land.

Fig. 1Map of the study area, Phewa Watershed

2.2. Data sources

This study used both primary and secondary sources. Primary data were gathered directly from the field, while secondary data were gathered from published and unpublished records from both private and public organizations. The mapping of the landslide inventory was achieved by analyzing Google Earth images, historical landslide data, and fieldwork. The Phewa watershed’s landslide susceptibility map was created utilizing a variety of maps, including those for land cover/use, DEM, slope, aspect, plan curvature, road distance, stream proximity, profile curvature, and precipitation. The Department of Survey, Government of Nepal, generated topographic maps that were then used to digitize contour layers to construct the DEM. From the DEM, slope maps and aspect maps were generated. The drainage map was used to determine river distances, which were then buffered based on the rivers’ chronological sequence. The topographic map was used to compute the length of the roads. The Department of Hydrology and Meteorology, Government of Nepal, served as the source of rainfall data for the risk assessment (Table 2).

Table 2Types and sources of data that used for mapping

Topographical variable	Value	Resolution	Source
DEM	729-2434 m	30 m	ALOS PALSAR RTC, openly available at https://asf.alaska.edu/
Slope	0-70.96°	30 m	ALOS PALSAR RTC, openly available at https://asf.alaska.edu/
Aspect	–1-360	30 m	ALOS PALSAR RTC, openly available at https://asf.alaska.edu/
Plan curvature	<-0.576->0.371	30 m	ALOS PALSAR RTC, openly available at https://asf.alaska.edu/
Profile curvature	<-0.545->0.386	30 m	ALOS PALSAR RTC, openly available at https://asf.alaska.edu/
LULC (Water bodies, built-up, forest, agriculture, flooded area)	–	30 m	Sentinel Image, openly available in [ESA] at https://www.esa.int.
Stream	Proximity to <100 m->500 m	30 m	Survey department
Precipitation	2,055-4,003 mm	30 m	Department of hydrology and meteorology
Road	Proximity to <100 m->500 m	30 m	Survey department
Landslide Inventory (Training and Testing dataset)	–	30 m	Google Earth Image and Field Survey

2.2.1. Landslide Inventory

Scientists and technicians can determine comparable circumstances in other locations and evaluate the possible risk of future landslides by examining the features and causes of previous landslides. Inside the study area, a total of 46 landslides were mapped. Using the Subset Features Tools in ArcGIS, the landslides in the landslide inventory were randomly divided into a training set area (80 %) with 37 points and a test set area (20 %) with 9 points. These were the dependent variable (landslide polygon) data that were utilized to train the model. Landslides vary in size according to their geography and location (Fig. 2). The approach of using a landslide inventory map based on historic landslide events, field site documentation, and satellite imagery interpretation is a common and effective method for assessing landslide risk. The division of the landslide inventory into a training set and a test set is a standard practice in machine learning and geospatial analysis, allowing for the development and validation of predictive models. The choice of a 80 % training set and a 20 % test set is a reasonable approach, as it provides enough data for both model development and validation.

Fig. 2Landslide inventory map

2.2.2. Landslide triggering factors

The landslide triggering factors are independent variables used to model the landslide susceptibility mapping. Future landslides can be predicted using the relationship between the environment and previous landslides. In this study, eight factors were utilized to identify the landslide-prone areas, which were divided into four categories: topographic, hydrological, geological, and other variables.

Google Earth images and DEM were used to create the landslide inventory map for this study, which was then verified using information from a field survey A high-resolution Google Earth Engine image displaying the landslide scar was used to examine the inventory map. We verified the satellite-based landslide database via a field study and looked at how different landslide types affected vegetation, water supply, physical structures, and the lives of people. The Phewa watershed's forest and ecology have suffered significantly as a result of the landslide hazard. In comparison to other conditioning elements, the slope angle is a factor that is equally significant but has a bigger influence. Due to their steepness and potential for increased environmental stress, areas with greater slopes typically produce more landslides. ALOS DEM was used to calculate the slope of the study area. There are six different classifications of slope [34].

In ArcGIS 10.7.1, the computed slope was then classed into five categories: 0-7.792°, 7.793-16.14°, 16.15-23.1°, 23.11-30.33°, 30.34-39.79°, and 39-8-70.96°. The research area has a slope that ranges from 0 degrees to 70.96 degrees. Most of the research area was highly rocky and steep (Fig. 3(a)). Aspect is also another crucial LSM conditioning component. It influences stratum instability by regulating moisture in the rocks and soil as a result of wind and varying the amount of solar exposure in the study area. Nine classes (Flat, North, Northeast, East, Southeast, South, Southwest, West, and Northwest) were formed using angular values based on the aspect’s facing orientation using ArcGIS 10.7.1 after receiving input through the ALOS DEM [27] (Fig. 3(b)). One of the most crucial initiating elements for a landslide is precipitation [35]. Most of the landslides occurred in areas with the most rainfall. The departments of hydrology and meteorology, Government of Nepal provided 32-year data (1990-2021), which were interpolated in GIS to generate a rainfall map (Fig. 3(c)). A sentinel image available from ESA was used to classify Land use/Cover area. Five classifications, mainly agriculture, water bodies, forest, built-up, and flooded areas, were used as the training dataset. The Land Use/Cover factor was scaled to a spatial resolution of 30 m to ensure smooth model operation. According to classification data, 44.95 % of the area was covered by forest, 38.52 by agricultural, 10.76 % by built-up areas, 4.28 % by water bodies, and 1.47 % by flooded land (Fig. 3(d)).

Another crucial factor for Landslide Susceptibility Mapping is the closeness to streams because these areas are more likely to experience landslides and soil erosion due to the erosive action of water. Based on the Euclidean distance, six buffer zones have been defined around the stream (Fig. 4(a)). Proximity to roads, similar to stream proximity, is a significant factor that directly contributes to landslides. Furthermore, a recently constructed road could contribute to slope instability. For the road map to function properly, the spatial resolution was rescaled to match the resolution of the ALOS DEM. Due to this, six categories were created based on the Euclidean distance from the centerline of roadways: below 100 meters, 100-200 meters, 200-300 meters, 400-500 meters, and above 500 meters [27] (Fig. 4(b)). The surface that is perpendicular to the slope direction is depicted by the plan curvature. It typically comes in three types: laterally convex, laterally concave, and linear, each represented by a positive, negative, or zero number. It was divided into three classes based on Wang et al. [27]: below –0.576, –0.576-0.371, and above 0.371 (Fig. 4(c)). The surface that is parallel with the direction of the steep slope is shown by the profile curvature. It typically comes in three varieties: upwardly convex, upwardly concave, and linear, each represented by a negative, positive, or zero value. Three classifications were developed: below –0.545, –0.545-0.386, and above 0.386, following Wang, et al. [27] (Fig. 4(d)).

Fig. 3Map of a) slope, b) aspect, c) rainfall/precipitation, d) land use/cover

Fig. 4Map of a) stream proximity, b) road proximity, c) plan curvature and d) profile curvature

3. Result

3.1. Landslide susceptibility map

The statistical approach used frequency ratio methods to analyze the quantitative relation between landslides and conditioning factors. The calculation of the landslide susceptibility model (LSM) was performed using Eq. (3) and Eq. (4) on the ArcGIS platform. The landslide susceptibility map was reclassified into five classes using the natural breaks method, with the largest area in the low class (43.226 %), followed by very low (35.234 %), moderate (19.46 %), very high (0.846 %), and high (0.63 %), (Fig. 5 and 6).

Fig. 5Landslide susceptibility map

Fig. 6Classes of landslide susceptibility map

Fig. 7Feature Importance for LSM

4. Model performance and validation

The validity and accuracy were assessed using the area under the curve (AUC) of receiver operating characteristic (ROC) curves, a quantitative method for evaluating the landslide susceptibility model based on prediction accuracy by using ArcGIS 10.7.1 tool. The true positive rate ($y$-axis) and false positive rate ($x$-axis) of the training and validation landslides were used to generate the AUC rate curves. In this work, 80 % of the landslides were utilized to train landslide susceptibility models, while the remaining 20 % were used for model validation. Fig. 8 demonstrates that the frequency ratio model has an AUC value of 0.717, indicating good performance. The AUC value lies between 0 and 1, and the higher its value, the greater the prediction accuracy of the model.

Fig. 8AUC Curve for validation

As the AUC value above 0.5 is generally considered to indicate a satisfactory outcome in Receiver Operating Characteristics (ROC) analysis, the values above 0.7 in this study suggest a good model. Oh, et al. [32] evaluated the effectiveness of the landslide susceptibility analysis using frequency ratio (FR) with iterative random sampling and observed that landslide susceptibility maps were validated with each validation dataset. They conclude that FR enables knowledge-driven analyses of the causal factors to be integrated into the landslide susceptibility analysis while also being widely used in other areas. Acharya, et al. [37] updated an inventory of Bhotang, Nepal using ten conditioning factors and a Google Earth Pro image. They found that the Area Under Curve (AUC) value of FR for validation was 0.713. Budha, et. al. [38], described the FR of the Panchase area, with a success rate of 0.76 AUC, indicating that the model has 76 % predictability. Onagh, et. al. [39] found an AUC of 0.762 during their model validation in the Uttarakhand district of northwestern India. In the Kulekhani watershed of Nepal the AUC for the statistical index method was 0.7578 [40]. Consequently, the AUC found in this study area is 0.717, which is quite comparable to that observed in previous studies in the Himalayas.

5. Discussion

The PR determined the weight of the conditioning factors in the study. Land use/cover was found to be the most important factor, followed by proximity to a stream, proximity to a road, precipitation, and slope. Aspect, plan curvature, and profile curvature were found to be less significant (Fig. 7). This supports previous findings of Basnet, et al. [41] that identified land use, road presence, precipitation, streams, aspect, and slope as primary triggers for landslides in their study area. The soil loss from the Phewa watershed because of the influence of LULC change was estimated by Bista and Basnet [42] using GIS, Remote sensing, and the RUSLE tool. They found that soil loss was 29.3 t/ha/yr in 2007 and 25.4 t/ha/yr in 2017. According to a study conducted by Rowbotham and Dudycha [43], dip slopes, particularly those under cultivation, were identified as the primary sources of slope collapse. The study area was significantly influenced by the closeness of roads, which were frequently built quickly lacking proper technical planning, drainage system, or slope maintenance.

In contrast, a study by Vuillez, et al. [44] found that 84 % of newly generated landslides in the Phewa watershed occur within 40 meters of a road, and more than 40 % of them cross a road. Surface water concentration and slope decline caused by these roads also contribute to slope instability. The presence of steeper slopes increases the likelihood of deformation and failure. Solar radiation and other aspects affect vegetation coverage, surface weathering, and evaporation, all of which influence landslide occurrence. Cooperation between communities upstream and downstream is crucial in reducing landslide risk in the area, where people reside at both levels. Media reports have highlighted landslide events in the Phewa watershed, with heavy rainfall and sloppy terrain being contributing factors. Sharma [45] reported that on August 25, 2014, a landslide triggered by continuous rainfall resulted in three fatalities and one person being reported missing. The incident also led to the destruction of the Typical Restaurant & Lodge located on the shores of Phewa Lake at Anadu, Pumdibhumdi in Pokhara-22. According to the District Police Office in Kaski on July 29, 2015, five people, including women, were killed and 25 others were injured when 12 houses in Dundakhet, Pokhara-23 were buried by landslides [46-48]. On July 10, 2020, around 3:00 a.m., landslide debris buried a house in Gothadi, Sarangkot, killing five people [49]. In an experiment by Budha, et al. [50], the Panchase area adjacent to this study area revealed a larger amount of landslide area in forests and cultivated lands, which are 0.762 km² and 0.358 km² respectively, accounting for 75 % of the total landslide area. Basnet, et al. [41] discovered similar results in a section of the Phewa watershed, indicating that forests have a higher frequency of landslides. Previous research has shown that forests have a higher frequency of landslides. However, this study found that landslides are less concentrated in forested areas because the shadows obstruct landslide images in satellite imagery. It is crucial to carefully choose plant species because there is a considerable chance that landslides will occur in both farmed and forested areas. Landslides can be caused by forests, so to reduce the risk of landslides in the protected Phewa watershed, planners should take into account the importance of deeply rooted vegetation and minor drainage improvements. In determining landslide risk and developing risk mitigation strategies, it is crucial to study the landslide susceptibility mapping because the Phewa watershed is a protected watershed and a part of the Ramsar site.

6. Conclusions

A landslide susceptibility map of the Phewa watershed was generated during this study using a variety of maps and data gathered through fieldwork, past landslide data, satellite data, and Google Earth data. It recognized 46 landslides across the study area and used eight variables, including slope, aspect, land use/coverage, precipitation, proximity to roads and streams, profile curvature, and plan curvature. The frequency ratio (FR) model was applied in the Phewa watershed to produce LSMs, which led to a map that was categorized into five levels: very low, low, moderate, high, and very high. According to the study, the most important variables influencing landslide risk in this region were slope, precipitation, land use/cover, and proximity to streams and roads. It also looked at the relationship between rainfall and landslides, finding that a higher frequency of landslides was linked to increased rainfall. The Area under the Curve (AUC) score of 0.717 for the susceptibility model indicated good performance. The LSMs offer important information to researchers, communities, government authorities, and planners to deal with the risks of landslides in the area, even though they cannot anticipate the frequency or timing. The results of this study will be crucial for comprehending the risk of landslides and developing strategies to mitigate their impact in the protected watershed of the Phewa Lake area.