Spatial Landslide Risk Assessment at Phuentsholing, Bhutan

Landslides are one of the most destructive and most recurring natural calamities in the Himalayan region. Their occurrence leads to immense damage to infrastructure and loss of land, human lives, and livestock. One of the most affected regions is the Bhutan Himalayas, where the majority of the landslides are rainfall-induced. The present study aims to determine the hazard and risk associated with rainfall-induced landslides for the Phuentsholing region located in the southwestern part of the Bhutan Himalayas. The work involves developing a landslide risk map using hazard and vulnerability maps utilizing landslide records from 2004 to 2014. The landslide hazard map was generated by determining spatial and temporal probabilities for the study region. The spatial probability was computed by analyzing the landslide contributing factors like geology, slope, elevation, rainfall, and vegetation based on comprehensive field study and expertise about the area. The contributing factors were divided into various classes and the percentage of landslide occurrence under each class was calculated to understand its contributing significance. Thereafter, a weighted linear combination approach was used in a GIS environment to develop the spatial probability map which was multiplied with temporal probabilities based on regional rainfall thresholds already determined for the region. Consequently, vulnerability assessment was conducted using key elements at risk (population, land use/land cover, proximity to road, proximity to stream) and the weights were provided based on expert judgment and comprehensive field study. Finally, risk was determined and the various regions in the study area were categorized as high, medium, and low risk. Such a study is necessary for low-economic countries like Bhutan which suffers from unavailability of extensive data and research. The study is conducted for a specific region but can be extended to other areas around the investigated area. The tool can serve as an indicator for the civil authorities to analyze the risk posed by landslides due to the rapid infrastructure development in the region.


Introduction
Landslides are considered one of the most prominent and devastative natural disasters leading to immense loss of human lives and infrastructure [1]. The triggering factor for the majority of the landslides is either rainfall or earthquake of which rainfall-triggered landslides can be categorized as very high and intense. The Himalayan region experiences very high rainfall and contributes over 70% of global fatal landslides [2]. Landslides in the Himalayas occur mainly due to monsoonal precipitation leading to huge financial loss and often disrupting the only livelihood of people [3,4]. The global fatal landslide database developed by [2] shows a total of 64 deaths in Bhutan due to nonseismic landslide events from 2004 to 2017, but the actual number of deaths is expected to be very high. The menace is expected to escalate because of the cutting of hills and deforestation emphasizing infrastructure development linked to a rise in population.
The studies on rainfall-induced landslides are quite varied and focus on several aspects like threshold modeling [5,6], monitoring systems [4,7], and computational techniques [8], along with GIS approach. Apart from landslide forecasting, landslide mitigation is important to understand the landslide hazard and the associated risk. Several works have been carried out to develop risk maps using the topography of the area susceptible to landslides [9][10][11]. Others have used ground conditions like the combined effect of rainfall and topography of the region to develop risk maps (South Italy [12], San Francisco Bay area [13]). Ho et al. [14] and Wong [15] used population density along with the natural contributing landslide factors to develop risk maps. Hadmoko et al. [16] utilized a vector-based GIS method to analyze several factors affecting landslides to develop hazard maps of the Central Java region. Pradhan and Abdulwahid [17] developed landslide risk maps using logistic regression and multihazard cases developed by considering precipitation data. Althuwaynee and Pradhan [11] developed a semi-quantitative landslide risk assessment using GIS-based analysis for Malaysia and were able to successfully show the predictive capability of the model. Liu and Miao [18] determined the landslide hazard, vulnerability, and risk in China by using certainty factor and a logistic regression model. Goorchi et al. [19] developed a new technique to analyze the risk linked to a rock slope, including the probability of failure and vulnerability rates. Saldivar-Sali and Einstein [20] proposed a regional landslide risk rating system in the Philippines. Aljeano et al. [21] developed a similar approach to analyze rock fall risk in quarries and the results depicted close correlation to the experimental work conducted for unstable slopes typically used for road construction. The effect of soil properties on slope instability has also been studied by several researchers using reliability-based optimization techniques, wherein the spatial variability related to shear parameters of soil is analyzed [22,23]. Cardarilli et al. [24] conducted an investigation to understand the territorial resilience towards landslide in the Southern Italy region by analyzing the correlations among the causative factors with an aim to reduce loss of lives and economic damage.
All these aforementioned studies have one thing in common, that is, the availability of sufficient data and prior research work. The present study has been performed for the Bhutan Himalayas, which lacks in both aspects. The research conducted in the past is so minimal that a Scopus search using the keywords "landslide*" and "Bhutan" resulted in only 36 documents. The publication of the latest landslide global disaster database has led several researchers to focus on this area and in the recent years, various research works have been conducted to better understand the landslides in the region. The majority of the landslide incidences are triggered due to rainfall, therefore we have focused solely on rainfall-induced landslides. The previous research in Bhutan has been towards forecasting, which led to the development of several thresholds using various statistical [3,25,26] and probabilistic [27] approaches. However, landslide mitigation and understanding the elements at risk are yet to be conducted for any part of Bhutan. Therefore, the present study is a preliminary step towards landslide risk assessment with an aim to perform a complete and comprehensive risk assessment in the future. The study was focused only on the Phuentsholing region, primarily due to two main reasons, the first being its significance as a major trade route and the second being the availability of sufficient data to conduct risk analysis. The work presented reflects a very simple but effective approach to determine landslide risk by computing hazard and vulnerability. Hazard comprises three different components (spatial, temporal, and magnitude) and these are multiplied directly, as an independent relationship is assumed between the components [28,29]. However, the unavailability of a detailed landslide inventory map results in the exclusion of the magnitude component. In the present study, the spatial component was calculated using a weighted sum analysis based on landslide contributing factors which was multiplied with the temporal probabilities determined based on regional rainfall thresholds [26] to develop the landslide hazard map. The predisposing factors for spatial aspect were identified based on a comprehensive field study conducted by the authors and examining the historical landslide locations in the study area. Landslide vulnerability was computed using weighted overlay analysis on key elements at risk by providing weights based on expert opinion. Finally, a landslide risk map was developed by combining both the hazard and vulnerability maps. Such an approach is useful for a country like Bhutan which lacks a proper landslide mitigation taskforce along with an operational early warning landslide system [30]. The present study would give an opportunity to develop risk assessment tools to analyze the slide activity and provide recommendations for mitigation purposes.
Geosciences 2020, 10, x FOR PEER REVIEW 3 of 16 unavailability of a detailed landslide inventory map results in the exclusion of the magnitude component. In the present study, the spatial component was calculated using a weighted sum analysis based on landslide contributing factors which was multiplied with the temporal probabilities determined based on regional rainfall thresholds [26] to develop the landslide hazard map. The predisposing factors for spatial aspect were identified based on a comprehensive field study conducted by the authors and examining the historical landslide locations in the study area. Landslide vulnerability was computed using weighted overlay analysis on key elements at risk by providing weights based on expert opinion. Finally, a landslide risk map was developed by combining both the hazard and vulnerability maps. Such an approach is useful for a country like Bhutan which lacks a proper landslide mitigation taskforce along with an operational early warning landslide system [30]. The present study would give an opportunity to develop risk assessment tools to analyze the slide activity and provide recommendations for mitigation purposes.

Study Area
The study area is Phuentsholing, situated in southwestern Bhutan, and is located from 26°50' N to 26°59' N latitude and 89°15' E to 89°28' E longitude covering an area of around 131.7 km 2 . The area is further divided into six subdivisions (known as chiwog), which are Lingdaen, Pachhu, Dzongkhag Tromde, Dophugchen Wangduegatshel, Deling Marpji, andChong Geykha Dophulakha ( Figure 1). Landslides in the study region are mainly due to intense monsoon rainfall, weak geology, and anthropogenic activities. The landslide typology in the study area can be characterized as rock fall, rockslide, debris flow, debris slide, and earth slide [24,25]. Gansser [31]  Landslides in the study region are mainly due to intense monsoon rainfall, weak geology, and anthropogenic activities. The landslide typology in the study area can be characterized as rock fall, rockslide, debris flow, debris slide, and earth slide [24,25]. Gansser [31] described the regional geological setting of the Bhutan Himalayas and characterized the majority of the country as part of the Himalayan Crystalline Complex. The areas that are most prone to landslides are heavily fractured and weathered rocks of phyllite, slate, and schist that contain high amounts of clay [26,27]. The study region has been identified as geologically weak with the presence of several shear zones, including weathered and folded phyllite rocks turning to clay [25].
The Phuentsholing-Thimphu highway, which is a part of the study region, is a vital lifeline as the majority of the import and export goods of the country are carried out through this route. The landslides lead to blockage of roads at several sections leading to a considerable amount of economic losses in repair and maintenance every year. The highway was the first motorable road of the country, constructed in 1961, and it crosses through areas susceptible to landslide, subsidence, erosion, and rockfalls. The instability in certain sections of the road has been as long as 20-50 years old. A 2005 report on the landslide investigations along the highway showed that a landslide at a 360 m section of the highway between the 17 km and 18 km point, known as Sorchen landslide, was so severe that the road alignment, which had four bends, eventually was left with only one bend over time. Figure 2 depicts two damage pictures taken during landslide events in Chukha Dzongkhag along the Phuentsholing-Thimphu Highway after the 2017 monsoon.
Geosciences 2020, 10, x FOR PEER REVIEW 4 of 16 geological setting of the Bhutan Himalayas and characterized the majority of the country as part of the Himalayan Crystalline Complex. The areas that are most prone to landslides are heavily fractured and weathered rocks of phyllite, slate, and schist that contain high amounts of clay [26,27]. The study region has been identified as geologically weak with the presence of several shear zones, including weathered and folded phyllite rocks turning to clay [25]. The Phuentsholing-Thimphu highway, which is a part of the study region, is a vital lifeline as the majority of the import and export goods of the country are carried out through this route. The landslides lead to blockage of roads at several sections leading to a considerable amount of economic losses in repair and maintenance every year. The highway was the first motorable road of the country, constructed in 1961, and it crosses through areas susceptible to landslide, subsidence, erosion, and rockfalls. The instability in certain sections of the road has been as long as 20-50 years old. A 2005 report on the landslide investigations along the highway showed that a landslide at a 360 m section of the highway between the 17 km and 18 km point, known as Sorchen landslide, was so severe that the road alignment, which had four bends, eventually was left with only one bend over time.

Data Used
The understanding of various causative factors and preparation of thematic layers are necessary to develop hazard and vulnerability maps. In the present study, a risk map was developed by analyzing the hazard contributory factors (geology, slope, elevation, vegetation, rainfall) with vulnerability factors (land use, population, proximity to road and stream). Hazard has been defined as a probability of a specific danger likely to occur in a given time period whereas risk refers to hazard multiplied by the potential loss of life [20].
Geology is one of the most significant contributory factors to landslide initiation [32]. The region has four significant geological groups, which can be subcategorized into a total of seven types, which are: (i) Baxa Group (Phuentsholing Formation, Pangsari Formation), (ii) Daling-Shumar Group (Daling Formation, Orthogneiss unit, Shumar Formation), (iii) Lesser Himalayan Zone (Jaishidanda Formation), and iv) structurally lower Greater Himalayan section (lower metasedimentary unit) ( Figure 3). The Baxa Group can be characterized as heavily fractured and weathered phyllites, slates, and schists [25]. The Daling-Shumar Group comprises tan weathering, cliff-forming quartzite and muscovite-biotite schist and phyllite [33]. The majority of the settlement area is dominated by the Baxa formation group wherein the thrust fault, strike, and dip of bedding and foliation are more concentrated. Similarly, such features were also delineated in the structurally lower Greater Himalayan Zone.

Data Used
The understanding of various causative factors and preparation of thematic layers are necessary to develop hazard and vulnerability maps. In the present study, a risk map was developed by analyzing the hazard contributory factors (geology, slope, elevation, vegetation, rainfall) with vulnerability factors (land use, population, proximity to road and stream). Hazard has been defined as a probability of a specific danger likely to occur in a given time period whereas risk refers to hazard multiplied by the potential loss of life [20].
Geology is one of the most significant contributory factors to landslide initiation [32]. The region has four significant geological groups, which can be subcategorized into a total of seven types, which are: i) Baxa Group (Phuentsholing Formation, Pangsari Formation), ii) Daling-Shumar Group (Daling Formation, Orthogneiss unit, Shumar Formation), iii) Lesser Himalayan Zone (Jaishidanda Formation), and iv) structurally lower Greater Himalayan section (lower metasedimentary unit) ( Figure 3). The Baxa Group can be characterized as heavily fractured and weathered phyllites, slates, and schists [25]. The Daling-Shumar Group comprises tan weathering, cliff-forming quartzite and muscovite-biotite schist and phyllite [33]. The majority of the settlement area is dominated by the Baxa formation group wherein the thrust fault, strike, and dip of bedding and foliation are more concentrated. Similarly, such features were also delineated in the structurally lower Greater Himalayan Zone. The other important factors are slope and elevation, which vary from 0°-75° and 183-2500 m, respectively. The majority of the study region (>80%) falls under 45° suggesting that the area is situated on a relatively flat plateau. More than half of the region has an elevation of less than 600 m and only 5% of the region is under high elevation zone (1500-2500 m).
The vegetation map was derived using the Normalized Difference Vegetation Index (NDVI) from the Landsat image (2017) and was classified into vegetated (86.76%) and nonvegetated areas (13.24%) (Figure 4c). NDVI always generates a value between −1 and +1 with high positive NDVI values (green) meaning high vegetation, while water usually has negative NDVI values (yellow and red), with urban features usually having values near zero. The maps were developed using the image analysis toolbar in ArcGIS 10.4. The other important factors are slope and elevation, which vary from 0 • -75 • and 183-2500 m, respectively. The majority of the study region (>80%) falls under 45 • suggesting that the area is situated on a relatively flat plateau. More than half of the region has an elevation of less than 600 m and only 5% of the region is under high elevation zone (1500-2500 m).
The vegetation map was derived using the Normalized Difference Vegetation Index (NDVI) from the Landsat image (2017) and was classified into vegetated (86.76%) and nonvegetated areas (13.24%) (Figure 4c). NDVI always generates a value between −1 and +1 with high positive NDVI values (green) meaning high vegetation, while water usually has negative NDVI values (yellow and red), with urban features usually having values near zero. The maps were developed using the image analysis toolbar in ArcGIS 10.4.  The study area is located in the subtropical and temperate climatic zone, which receives very high annual rainfall up to 6000 mm. As per the Köppen climate classification [34], the climate is subtropical monsoonal receiving 10 times more rainfall during the wettest summer month as compared to dry winter months. The contribution of monsoonal rainfall (June-September) is 80% of the annual rainfall. However, pre-and post-monsoonal months can also suffer from heavy rainfall due to local convection and tropical cyclone activity over the Bay of Bengal [35], with January and February being relatively drier [25]. Figure 4d shows the rainfall map based on the annual rainfall measurements from 2004 to 2014 collected from the National Center for Hydrology and Meteorology of the Royal Government of Bhutan.
A total of 64 landslides occurred from 2004 to 2014, which was collected from Project Dantak, the Border Roads Organization, Government of India. The data included the coordinates of the landslide events which were collected from the scarp along with the dates of occurrences. The number of landslide events was counted as a single landslide event if multiple events occurred on the same day and coordinates were selected where maximum damage or loss of land prevailed. Table  1 shows the number of landslides in each class of landslide contributing factors.  The study area is located in the subtropical and temperate climatic zone, which receives very high annual rainfall up to 6000 mm. As per the Köppen climate classification [34], the climate is subtropical monsoonal receiving 10 times more rainfall during the wettest summer month as compared to dry winter months. The contribution of monsoonal rainfall (June-September) is 80% of the annual rainfall. However, pre-and post-monsoonal months can also suffer from heavy rainfall due to local convection and tropical cyclone activity over the Bay of Bengal [35], with January and February being relatively drier [25]. Figure 4d shows the rainfall map based on the annual rainfall measurements from 2004 to 2014 collected from the National Center for Hydrology and Meteorology of the Royal Government of Bhutan.
A total of 64 landslides occurred from 2004 to 2014, which was collected from Project Dantak, the Border Roads Organization, Government of India. The data included the coordinates of the landslide events which were collected from the scarp along with the dates of occurrences. The number of landslide events was counted as a single landslide event if multiple events occurred on the same day and coordinates were selected where maximum damage or loss of land prevailed. Table 1 shows the number of landslides in each class of landslide contributing factors.

Method
The method revolves around the determination of some important landslide risk factors (i.e., hazard, vulnerability, elements at risk, and potential damage). Landslide risk is defined as the maximum potential loss due to a landslide event for a definite area and during a certain period [18]. Hence, landslide risk can be calculated as: A complete quantitative risk analysis includes explicitly the economical exposure of elements at risk like buildings, infrastructures, and assets. However, due to the unavailability of such factors, population is used as a proxy to obtain a relative estimation of exposure and the final output is a map of qualitative risk classes.
The first factor is "hazard" which has been often erroneously used synonymously with "susceptibility" [36], but there are major differences between them. Landslide hazard is the probability of a landslide event of a certain magnitude within spatial and temporal limits [11], whereas landslide susceptibility is the likelihood of landslide occurrence depending on local conditions [37].
The present study determines the landslide hazard by considering spatial (P s ) and temporal (P T ) components [33,34]. The spatial aspect was extracted by analyzing all the landslide contributing factors and applying a weighted linear combination (WLC) method. WLC is quite often used by landslide researchers, but an accurate method of determining weights is not defined and is based on expert opinion [38]. It is based on determining weighted averages of landslide causative parameters pertinent to the study area and available data. The parameters are classified and multiplied with the provided weight in a GIS overlay environment, and the final output is determined by adding weighted averages.
The temporal probability was calculated using the rainfall threshold determined by Dikshit et al. [26] using the Poisson probability model. Readers are referred to [3,39,40] for the detailed explanation on determining the temporal landslide probability.
Landslide vulnerability is the degree of loss of a particular element at risk due to the occurrence of a certain magnitude of the event [41,42]. Glade and Crozier [43] highlighted that there is no fixed technique for quantitative vulnerability assessment of elements at risk with respect to different landslide magnitudes and types. Therefore, vulnerability is usually determined by either assigning a fixed value to elements at risk (0 for no damage and 1 for complete damage), or assigning values based on expert opinion and knowledge about the area [44]. Thereafter, the potential damage is calculated by combining the vulnerability and elements at risk. In the present study, the potential damage map was developed using the WLC method by providing the exposure of the weights to causative factors expressing both elements of vulnerability and exposure. The factors considered for vulnerability were land use/land cover (LULC), proximity to road, and proximity to stream, and for exposure was population.
The classification of derived thematic layers can be conducted using several techniques like equal intervals, natural breaks, and quantile and statistical techniques which are default GIS processes [45]. In the present study, the natural breaks method was used to classify values as it is quite popular among the researchers [45][46][47].

Landslide Hazard Assesment
As described in the methodology, the spatial component (susceptibility) was based on the WLC method and the temporal probability was determined based on the thresholds obtained using an empirical-based approach. The susceptibility map was based on landslide causative factors and the weights were provided as per Equation (3).
where LS is landslide susceptibility. The weights are assigned using GIS overlay capability, based on its influence on landslide in comparison to other factors and performing a weighted sum to combine the factors [38]. The fault data were not added as the majority of the landslides are rainfall-triggered and can be attributed as a low-seismic region [48]. The temporal probability is calculated by multiplying annual exceedance probability and the probability of landslide occurrence. The annual exceedance probability is the probability of threshold being exceeded in a specific year, which is determined using the Poisson probability model. The maps derived using both the components resulted in the hazard map. Figure 5 represents the landslide hazard map, which has been divided into five classes: very low The results show that the central-eastern parts of the area are more highly hazardous than the 4 km radius of central Phuentsholing, which comes under high hazard zone. The eastern and western parts of the region fall under very low to medium hazard zone. This could prove to be very helpful in landslide prevention mapping and could be used for land use planning and in inhibitory works by governmental agencies. The results show that the central-eastern parts of the area are more highly hazardous than the 4 km radius of central Phuentsholing, which comes under high hazard zone. The eastern and western parts of the region fall under very low to medium hazard zone. This could prove to be very helpful in landslide prevention mapping and could be used for land use planning and in inhibitory works by governmental agencies.

Landslide Vulnerability Assessment
The vulnerability assessment is defined as the expected degree of loss for elements at risk due to a certain event, with values ranging between 0 (no damage) and 1 (complete damage). Its assessment is based on estimation of vulnerability of life and property. A detailed risk assessment would involve considering several elements at risk, like population density, residential and commercial buildings, and their relationship with the type of landslides. However, due to the unavailability of such detailed data, the component of elements at risk could not be considered. The potential damage map was developed using four layers of elements related to exposure and vulnerability, which were LULC, population, proximity to road, and proximity to stream Figure 6ad and a similar approach to susceptibility mapping was used. Apart from expert opinion, the weights (Equation (4)) provided indicate a logical progression. For example, the built-up area would have more loss due to a higher population, thereby providing them with a higher multiplier. Such regions are specifically used for either commercial or residential purposes or both. These areas can be dangerous on the hills as most of the houses are built on slopes and the roads are constructed by hill cutting. Agriculture covers a significant portion of the study area as it is a means of livelihood for the people staying there. Miscellaneous region includes water bodies and mining sites, and any effect on such sites due to landslides can lead to economic and ecological impact. Similarly, the road network is a lifeline to the people and to the economy of the country, therefore a higher multiplier was provided. This makes the proximity to road factor critical as road cuts are quite common in this region and the risk decreases as distance increases away from the road network. The relevance to the proximity to stream factor is due to the presence of stream channels which could suffer from erosion and hill cutting especially during heavy monsoons. The cutting of hills increases slope instability, making the area more vulnerable to landslides, and it usually increases in regions of higher drainage

Landslide Vulnerability Assessment
The vulnerability assessment is defined as the expected degree of loss for elements at risk due to a certain event, with values ranging between 0 (no damage) and 1 (complete damage). Its assessment is based on estimation of vulnerability of life and property. A detailed risk assessment would involve considering several elements at risk, like population density, residential and commercial buildings, and their relationship with the type of landslides. However, due to the unavailability of such detailed data, the component of elements at risk could not be considered. The potential damage map was developed using four layers of elements related to exposure and vulnerability, which were LULC, population, proximity to road, and proximity to stream Figure 6a-d and a similar approach to susceptibility mapping was used. Apart from expert opinion, the weights (Equation (4)) provided indicate a logical progression. For example, the built-up area would have more loss due to a higher population, thereby providing them with a higher multiplier. Such regions are specifically used for either commercial or residential purposes or both. These areas can be dangerous on the hills as most of the houses are built on slopes and the roads are constructed by hill cutting. Agriculture covers a significant portion of the study area as it is a means of livelihood for the people staying there. Miscellaneous region includes water bodies and mining sites, and any effect on such sites due to landslides can lead to economic and ecological impact. Similarly, the road network is a lifeline to the people and to the economy of the country, therefore a higher multiplier was provided. This makes the proximity to road factor critical as road cuts are quite common in this region and the risk decreases as distance increases away from the road network. The relevance to the proximity to stream factor is due to the presence of stream channels which could suffer from erosion and hill cutting especially during heavy monsoons. The cutting of hills increases slope instability, making the area more vulnerable to landslides, and it usually increases in regions of higher drainage density and higher elevation. The total population in the study area is 33,444 with 18,043 males and 15,401 females as per census in Bhutan, 2017 [49].

Landslide Risk Assessment
The type of damage linked to landslide risk in the present study has been confined to loss of life and/or property along with a road network. The risk map was developed by multiplying the landslide hazard and vulnerability map. Such a study would also help in estimating the vulnerable area and the impact from any landslide event in the foreseeable future. In order to understand the vulnerability to human lives, the risk to population was estimated and the high-risk zones were identified ( Figure 8). The obtained map was classified into five classes (very low (0-0.12), low (0.13-0.21), moderate (0.22-0.37), high (0.38-0.55), and very high (0.56-1)) based on the natural break classification technique [47]. Table 2 depicts the area and affected population under various risk  zones whereas Table 3 determines the various landslide risk zones for all six subdivisions (chiwogs). It can be inferred that 20.2% of the population with an area of 15.13 km 2 lives in high to high-risk landslide zones.

Landslide Risk Assessment
The type of damage linked to landslide risk in the present study has been confined to loss of life and/or property along with a road network. The risk map was developed by multiplying the landslide hazard and vulnerability map. Such a study would also help in estimating the vulnerable area and the impact from any landslide event in the foreseeable future. In order to understand the vulnerability to human lives, the risk to population was estimated and the high-risk zones were identified ( Figure 8). The obtained map was classified into five classes (very low (0-0.12), low (0.13-0.21), moderate (0.22-0.37), high (0.38-0.55), and very high (0.56-1)) based on the natural break classification technique [47]. Table 2 depicts the area and affected population under various risk zones whereas  Table 3 determines the various landslide risk zones for all six subdivisions (chiwogs). It can be inferred that 20.2% of the population with an area of 15.13 km 2 lives in high to high-risk landslide zones.

Conclusions
The occurrence of landslides and their growing incidence is a major threat to population and developments to the mountainous region, especially in developing countries. This has led to a need for comprehensive study on various aspects of landslides in such regions. However, such studies in the Bhutan Himalayas is very scarce and studies need to be conducted to understand the landslide effects in the region. Therefore, a landslide risk map has been developed for the Phuentsholing region of southwestern Bhutan. The study area has been selected since the major trade route, Phuentsholing-Thimphu Highway, which connects Bhutan with neighboring countries, passes through it and a landslide in this area can hamper people's lives and lead to economic losses. The method involved assessing the risk by multiplying the hazard and vulnerability maps. The landslide hazard map was developed by determining the spatial and temporal probabilities of landslides. The spatial probability was determined by performing a weighted sum analysis in an overlay environment for landslide causative factors (geology, slope, rainfall, vegetation, elevation), whereas temporal probability was determined using Poisson's probability model based on the already calculated regional thresholds. Similarly, landslide vulnerability was determined by providing weights based on expert knowledge gathered after carrying out extensive field study to various elements at risk involved. Finally, the landslide risk map was developed by combining the landslide hazard and vulnerability maps and was divided into five categories depending on various levels of risk associated. The analysis showed that 11.7% of the area comes under the high to very high category of risk, whereas 15.85% falls under moderate risk. One of the major problems for determining landslide risk in such a region is the unavailability of longer and precise rainfall and landslide records, with often under-reporting of the landslide damages. However, the present study can be considered as a preliminary step towards achieving a deep understanding of landslide risk at a local scale. The maps can be further modified and improved with the availability of more detailed data like the distribution of population in terms of gender, and the landslide type. The areas along the highway are expected to be at considerable risk in the near future as more hills are cut for residential and economic development purposes. Such an analysis may be used as a tool to develop a building code based on landslide risk, which Bhutan is presently lacking. This analysis would help the concerned authorities to understand the potential risk in certain sections and could be used for landslide mitigation purposes as there is a rise of infrastructure development in the region.
Author Contributions: Conceptualization, A.D. and R.S.; methodology, A.D. and S.A.; data curation, S.A.; writing-original draft preparation, A.D. and R.S.; writing-review and editing, B.P. and A.M.A.; supervision, B.P. and R.S. All authors have read and agreed to the published version of the manuscript.