Temporal Probability Assessment and Its Use in Landslide Susceptibility Mapping for Eastern Bhutan

Landslides are one of the major natural disasters that Bhutan faces every year. The monsoon season in Bhutan is usually marked by heavy rainfall, which leads to multiple landslides, especially across the highways, and affects the entire transportation network of the nation. The determinations of rainfall thresholds are often used to predict the possible occurrence of landslides. A rainfall threshold was defined along Samdrup Jongkhar–Trashigang highway in eastern Bhutan using cumulated event rainfall and antecedent rainfall conditions. Threshold values were determined using the available daily rainfall and landslide data from 2014 to 2017, and validated using the 2018 dataset. The threshold determined was used to estimate temporal probability using a Poisson probability model. Finally, a landslide susceptibility map using the analytic hierarchy process was developed for the highway to identify the sections of the highway that are more susceptible to landslides. The accuracy of the model was validated using the area under the receiver operating characteristic curves. The results presented here may be regarded as a first step towards understanding of landslide hazards and development of an early warning system for a region where such studies have not previously been conducted.


Introduction
Rainfall triggered landslides are one of the most devastating naturally occurring disasters across the world [1]. The global dataset of landslide hazards in the 2004-2016 period extracted from Reference [2] showed that almost 75% of the world's fatal landslides occurred in the Himalayan region. Bhutan is no exception to this, and is a part of one of the world's highly landslide-prone regions in the world [3]. The damage caused by landslides in this country has led to casualties and loss of land, affecting people's livelihoods and disrupting the transportation network, which is key to the country's economy. Most of the landslides in the Bhutan Himalayas are triggered by rainfall, especially during the monsoon period [4,5]. Therefore, it is imperative to identify the areas that could be affected by Water 2020, 12, 267 2 of 24 landslides, in order to reduce the probability of damage in the future. The key to achieving this is through a detailed landslide hazard assessment that will help civic authorities to curtail landslide damage through effective land use management.
Landslide hazard may be defined as the probability of a damaging landslide in a spatial ("where") and temporal ("when") context, along with the magnitude ("how large") of the event [6,7]. Landslide susceptibility is defined as the likelihood of landslide occurrence ("where") in an area depending on local terrain conditions [8]. It may be regarded as the first step towards analyzing hazard and risk. Various spatial assessment models for landslide susceptibility have been developed [9][10][11]. However, compared to spatial assessment, there have been fewer attempts to carry out temporal probability assessment (studies have been conducted in Nilgiris, India [12], Hoa Binh, Vietnam [13], and Cameroon [14]). The two main techniques used to assess temporal probability for future landslide occurrences are (i) analysis of potential slope failure and (ii) statistical analysis of past landslide events [13,15]. The first technique involves evaluation of the current slope conditions and determination of the probability for future slope instability, which may be difficult to apply in large study areas [16]. Statistical analysis of past landslide events may be done directly using records of the landslides identified in the study area or, alternatively, it may be performed indirectly by using information related to recurrence of the landslide-triggering events [17]. Direct analysis requires a long time span of historical landslide data which is extremely difficult to obtain, especially in underdeveloped countries. Therefore, an indirect approach analyzing the frequency of occurrence of rainfall was used in this study to determine temporal probability. Even though this approach did not require complete multi-temporal landslide inventory data, it required determination of the relationship between rainfall and landslide incidences. After the calculation of rainfall thresholds, the landslide temporal probability was computed based on the number of times precipitation exceeded the threshold value [16]. As the frequency of rainfall-induced landslides only evaluates how often landslides might occur, it therefore needs to be integrated with spatial probability (susceptibility) and temporal probability to develop a landslide hazard map [17,18].
The prediction of landslide incidences using rainfall thresholds has been successfully carried out for various regions, including Italy [19][20][21], New Zealand [22], Malaysia [23], and the Himalayan arc [24][25][26]. The calculation of rainfall thresholds for landslide triggering can be determined using three main approaches: (i) physically based models [27], (ii) empirical rainfall threshold models [28], and (iii) statistically based models [29]. The physically based models are linked to the physical attributes of the study region and can be difficult to apply in cases of unavailability of an extensive dataset. Empirical models calculate rainfall thresholds based on past rainfall events which led to landslide incidences. The threshold is usually obtained by drawing lower-bound lines to the rainfall conditions that resulted in landslides plotted in Cartesian, semi-logarithmic, or logarithmic coordinates [30]. Statistical models use statistical tools like Bayesian inference or logistic regression to calculate thresholds [31].
In the case of Bhutan, studies to date have focused on the southwest region covering the Phuentsholing-Thimphu Highway (known as the Asian Highway), which connects the national capital Thimphu with neighboring countries. These studies have primarily focused on rainfall estimation and spatial assessment, using various techniques such as a probabilistic approach [5,32], a semi-automatic algorithm approach [26], an empirical approach [4], and machine learning models [33]. The other major highway, Samdrup Jongkhar-Trashigang (S-T), situated in the eastern part of the country, has been neglected, and a landslide study in this region is yet to be conducted. The main aim of this study was to assess landslide susceptibility utilizing temporal rainfall for the S-T highway. The two objectives in this study were (i) to estimate the temporal probability, and (ii) to estimate the landslide susceptibility using a multi-criterion decision-making approach. We addressed three major research themes in the current study: (i) determination of rainfall threshold, probability estimation of the threshold being exceeded, and landslide probability after the threshold has exceeded; (ii) susceptibility of the region with respect to landslides; and (iii) validation of the thresholds and susceptibility map. For this, the rainfall thresholds were determined based on the relation between daily rainfall and past landslide events that occurred between 2014 and 2017. The thresholds were validated using the rainfall records daily rainfall and past landslide events that occurred between 2014 and 2017. The thresholds were validated using the rainfall records of 2018. Thereafter, the exceedance probability of the threshold was calculated and the temporal probability of landslides was determined using a Poisson model. Finally, a landslide susceptibility map was developed using the analytic hierarchy process (AHP), utilizing the determined threshold values. This study was the first attempt to conduct such an elaborate study for the eastern region of Bhutan. The results from the present work can be understood as a preliminary step towards setting up an operational landslide early warning system so that damage to the transportation corridor can be reduced and human lives can be saved.

Study Area
The study area was the Samdrup Jongkar-Trashigang (S-T) highway, which is a 180 km stretch of road located in the eastern part of Bhutan, which covers 1880 km 2 ( Figure 1). The region was selected as it connects eight districts (known as "dzongkhags" in Bhutanese) and is a major route for the people residing around the highway. The highway is critical as it is the only transportation network connecting eight dzongkhags in East Bhutan, and it is a lifeline for the people residing in these areas. The transportation corridor has a history of severe slope failures in the form of frequent landslides, rockfalls, and mudflow. The region falls in a sub-tropical zone, where heavy rainfall is frequent. The region receives its maximum rainfall during the monsoon seasons (June-September). The rainfall pattern in this region can be described as low-intensity and long-duration with occasional intermittent outbursts. Figure 1. Location of (a) Bhutan; (b) spatial distribution of landslides in the study area considered in this analysis. The zones were categorized based on spatial coverage of rain gauges, rain gauge location, and elevation difference.
The area includes five different geological groups, which are Baxa Group, Daling Shumar Group, and the Greater, Lesser, and Sub-Himalayan zones. The geology varies from the Baxa Group in the south to the Greater Himalayan Zone in the north. The rocks found in the region are predominantly dark grey to green fine-grained phyllites, slate varying from dark brown to black, and Figure 1. Location of (a) Bhutan; (b) spatial distribution of landslides in the study area considered in this analysis. The zones were categorized based on spatial coverage of rain gauges, rain gauge location, and elevation difference.
The area includes five different geological groups, which are Baxa Group, Daling Shumar Group, and the Greater, Lesser, and Sub-Himalayan zones. The geology varies from the Baxa Group in the south to the Greater Himalayan Zone in the north. The rocks found in the region are predominantly dark grey to green fine-grained phyllites, slate varying from dark brown to black, and fine-grained, medium-Water 2020, 12, 267 4 of 24 to thick-bedded quartzite with thin to very thin grey-black fine-grained phyllite interbeds [34,35]. The average thickness of the quartzite is about 100 m, but the individual bands of quartzite range from 10 cm to 2 m. The orientation of the latter is 48 • NW, with an average dip 40 • towards the slope direction of the slide. The quartzite in the landslide area undergoes brittle deformation with many irregular joints [36]. Figure 2 depicts examples of the damage caused by landslides along the highway.
Water 2020, 12, x FOR PEER REVIEW 4 of 24 fine-grained, medium-to thick-bedded quartzite with thin to very thin grey-black fine-grained phyllite interbeds [34,35]. The average thickness of the quartzite is about 100 m, but the individual bands of quartzite range from 10 cm to 2 m. The orientation of the latter is 48° NW, with an average dip 40° towards the slope direction of the slide. The quartzite in the landslide area undergoes brittle deformation with many irregular joints [36]. Figure 2 depicts examples of the damage caused by landslides along the highway.

Data
Data from a total of 347 landslides which occurred from 2014 to 2018 were collected from Border Road Organization, Government of India under Project DANTAK (Figure 1b). The data included the dates and geographic coordinates of landslide location based on field observations, interactions with locals, and media reports. The types of movement in the region based on Cruden and Varnes' [37] classification are: debris slides, rock falls, earth flows, and rotational landslides. The field visit conducted in November 2017 revealed the landslides to be shallow, with depths ranging up to few meters, and able to be mapped as single points. The yearly distribution (Figure 3a) of landslides shows that the majority of the landslides occurred in 2017 and 2018, whereas the monthly distribution

Data
Data from a total of 347 landslides which occurred from 2014 to 2018 were collected from Border Road Organization, Government of India under Project DANTAK (Figure 1b). The data included the dates and geographic coordinates of landslide location based on field observations, interactions with locals, and media reports. The types of movement in the region based on Cruden and Varnes' [37] classification are: debris slides, rock falls, earth flows, and rotational landslides. The field visit conducted in November 2017 revealed the landslides to be shallow, with depths ranging up to few meters, and able to be mapped as single points. The yearly distribution (Figure 3a) of landslides shows that the majority of the landslides occurred in 2017 and 2018, whereas the monthly distribution ( Figure 3b) shows that 89% of landslides occurred between the months of June and September. It is Water 2020, 12, 267 5 of 24 often difficult to determine the rainfall conditions responsible for failures, due to lack of rain gauge density and high distances between rain gauges and landslide points [26,38]. As multiple landslides can occur during a rainfall event, subsequent landslides for a single rainfall event after the initial failure were not considered for the threshold analysis. Thus, in this study, we defined a landslide event as "single landslide-triggering event", in which the landslide events after the initial failure were not considered for threshold estimation [16,26]. A rainfall event was defined as a period of continuous rainfall separated by dry (without rainfall) period. Using all these criteria, the total numbers of rainfall and landslide events during the threshold determination time frame (2014-2017) were 477 and 104, respectively. The zones were divided based on spatial coverage of rain gauges [4,26]. A buffer radius of 15 km around each rain gauge was selected to divide the region into zones. In terms of landslide events for respective zones, half of the landslide events occurred in Zone 2 (49) followed by Zone 3 (40) and Zone 1 (15).
Water 2020, 12, x FOR PEER REVIEW 5 of 24 ( Figure 3b) shows that 89% of landslides occurred between the months of June and September. It is often difficult to determine the rainfall conditions responsible for failures, due to lack of rain gauge density and high distances between rain gauges and landslide points [26,38]. As multiple landslides can occur during a rainfall event, subsequent landslides for a single rainfall event after the initial failure were not considered for the threshold analysis. Thus, in this study, we defined a landslide event as "single landslide-triggering event", in which the landslide events after the initial failure were not considered for threshold estimation [16,26]. A rainfall event was defined as a period of continuous rainfall separated by dry (without rainfall) period. Using all these criteria, the total numbers of rainfall and landslide events during the threshold determination time frame (2014-2017) were 477 and 104, respectively. The zones were divided based on spatial coverage of rain gauges [4,26]. A buffer radius of 15 km around each rain gauge was selected to divide the region into zones. In terms of landslide events for respective zones, half of the landslide events occurred in Zone 2 (49) followed by Zone 3 (40) and Zone 1 (15). The daily rainfall data used for this study ( Figure 4) was collected from two rain gauges at Deothang and Kanglung, managed by the National Centre for Hydrology and Meteorology, Bhutan The daily rainfall data used for this study ( Figure 4) was collected from two rain gauges at Deothang and Kanglung, managed by the National Centre for Hydrology and Meteorology, Bhutan Water 2020, 12, 267 6 of 24 (http://www.hydromet.gov.bt). The average cumulative yearly rainfalls for Deothang and Kanglung for 2014-2018 were 3495 mm and 1020 mm, respectively, of which 89% occurred during the monsoon season (June-September). The higher precipitation observed in Deothang is due to its location on the windward side of the mountain.
Water 2020, 12, x FOR PEER REVIEW 6 of 24 (http://www.hydromet.gov.bt). The average cumulative yearly rainfalls for Deothang and Kanglung for 2014-2018 were 3495 mm and 1020 mm, respectively, of which 89% occurred during the monsoon season (June-September). The higher precipitation observed in Deothang is due to its location on the windward side of the mountain.

Rainfall Threshold Estimation
A rainfall threshold determines the minimum rainfall conditions necessary for landslide initiation in a specific region [39], and various researchers have attempted to quantify thresholds using several approaches [39][40][41]. A recent review article [40] on rainfall threshold estimation explained in detail the various approaches currently in use, along with their merits and demerits. Of the various techniques, empirically based approaches are widely used because of the simplicity and ease with which they can provide an accurate approximation of minimum precipitation conditions. Various rainfall thresholds using different parameters have been developed, such as ID (intensityduration) [21], ED (cumulative event rainfall duration) [4], and AD (antecedent rainfall duration) [16,42]. The most commonly used rainfall variables for threshold estimation are daily rainfall, antecedent rainfall, and cumulative rainfall [40]. However, the choice of rainfall variable with which to determine thresholds is primarily dependent on the type of landslides in the region [20].
For the S-T highway, monsoonal rainfall occurs with interruptions and can be characterized mostly low-intensity and long-duration events along with occasional extreme events, making the choice of antecedent rainfall appropriate [30]. Antecedent rainfall is a significant factor for landslide triggering, especially in less impermeable soils, as it lowers soil suction and increases pore water pressure [13]. The use of antecedent rainfall was based on analysis of historical landslide pattern, the field visit, and previous studies conducted in other regions of Bhutan. Although estimation of the number of days to be considered to analyze the effect of antecedent rainfall was a challenge, it has been widely accepted that antecedent rainfall over 15 days to 30 days plays a crucial role for landslide initiation in the Himalayas [43]. The calculation of the antecedent period prior to landslide incidence is usually based on a trial-and-error approach, ranging from 3 days to 120 days [30,44,45].
For this study, the correlation between daily and antecedent rainfall conditions was analyzed for six different time periods (3,5,7,10,15, and 30 days) (Figure 5a-f). The analysis of the various antecedent rainfall time periods was conducted using the method proposed by Zezere et al. [42]. Blue dots represent the daily rainfall, whereas the orange points depict the antecedent rainfall values for respective time periods. The best discrimination between daily and antecedent conditions, according

Rainfall Threshold Estimation
A rainfall threshold determines the minimum rainfall conditions necessary for landslide initiation in a specific region [39], and various researchers have attempted to quantify thresholds using several approaches [39][40][41]. A recent review article [40] on rainfall threshold estimation explained in detail the various approaches currently in use, along with their merits and demerits. Of the various techniques, empirically based approaches are widely used because of the simplicity and ease with which they can provide an accurate approximation of minimum precipitation conditions. Various rainfall thresholds using different parameters have been developed, such as ID (intensity-duration) [21], ED (cumulative event rainfall duration) [4], and AD (antecedent rainfall duration) [16,42]. The most commonly used rainfall variables for threshold estimation are daily rainfall, antecedent rainfall, and cumulative rainfall [40]. However, the choice of rainfall variable with which to determine thresholds is primarily dependent on the type of landslides in the region [20].
For the S-T highway, monsoonal rainfall occurs with interruptions and can be characterized mostly low-intensity and long-duration events along with occasional extreme events, making the choice of antecedent rainfall appropriate [30]. Antecedent rainfall is a significant factor for landslide triggering, especially in less impermeable soils, as it lowers soil suction and increases pore water pressure [13]. The use of antecedent rainfall was based on analysis of historical landslide pattern, the field visit, and previous studies conducted in other regions of Bhutan. Although estimation of the number of days to be considered to analyze the effect of antecedent rainfall was a challenge, it has been widely accepted that antecedent rainfall over 15 days to 30 days plays a crucial role for landslide initiation in the Himalayas [43]. The calculation of the antecedent period prior to landslide incidence is usually based on a trial-and-error approach, ranging from 3 days to 120 days [30,44,45].
For this study, the correlation between daily and antecedent rainfall conditions was analyzed for six different time periods (3,5,7,10,15, and 30 days) (Figure 5a-f). The analysis of the various antecedent rainfall time periods was conducted using the method proposed by Zezere et al. [42]. Blue dots represent the daily rainfall, whereas the orange points depict the antecedent rainfall values for respective time periods. The best discrimination between daily and antecedent conditions, according to Water 2020, 12, 267 7 of 24 the method suggested by Reference [42], was observed for 30 day antecedent rainfall and was accepted as the metric for threshold calculation.
Water 2020, 12, x FOR PEER REVIEW 7 of 24 to the method suggested by Reference [42], was observed for 30 day antecedent rainfall and was accepted as the metric for threshold calculation.
(a) 3 Days The threshold determination was performed using a scatter chart for daily (RTH) and 30 day antecedent rainfall (R30ad), and was calculated for all three zones. The graph was generated using the rainfall and landslide data from 2014-2017. The threshold equation was obtained by using the lower  The threshold determination was performed using a scatter chart for daily (R TH ) and 30 day antecedent rainfall (R 30ad ), and was calculated for all three zones. The graph was generated using the rainfall and landslide data from 2014-2017. The threshold equation was obtained by using the lower end points in the scattered graph [13,16]. The threshold equations of various zones are depicted in Figure 6.
Water 2020, 12, x FOR PEER REVIEW 9 of 24 end points in the scattered graph [13,16]. The threshold equations of various zones are depicted in Figure 6.

Validation of Rainfall Threshold
The significance of any landslide study is determined by the validation of the results obtained. One review of rainfall threshold studies [40] emphasized the importance of validation of rainfall thresholds for conducting landslide studies. The rainfall threshold validation was performed using rainfall and landslide data from 1 January to 31 December 2018. During this period, a total of 52 landslide events occurred, out of which 80% of landslides (41) happened during the monsoon season. The threshold equation for Zone 1 was RTH = 150 − 0.24R30ad, and its validation for the monsoon of 2018 is depicted in Figure 7. The threshold exceedance axis depicts the value of RTH with respect to R30AD value for each day, wherein the positive values indicates landslide occurrence. This figure shows that a heavy measure of rainfall occurred, exceeding the threshold. The abrupt increase in the magnitude of daily rainfall or constant rise in 30 day antecedent rainfall is shown by the rise in the threshold curve.
During the validation period, the threshold was exceeded nine times, out of which seven times landslides occurred. No landslides were reported on 2 June and 7 June, even though the threshold was exceeded. This observation can be attributed to the fact that a landslide event does not generally happen with the increase in threshold curve, and sometimes happens a couple of days later as the result of a difference in pore pressure because of changes in the measure of antecedent rainfall. From 1 October to 30 December, there was no threshold exceedance, and no landslides occurred during that period. Similar validations were carried out for Zone 2 and Zone 3 threshold equations using 2018 rainfall and landslide data. These results indicate that the threshold model performed well for landslide forecasting in 2018.

Validation of Rainfall Threshold
The significance of any landslide study is determined by the validation of the results obtained. One review of rainfall threshold studies [40] emphasized the importance of validation of rainfall thresholds for conducting landslide studies. The rainfall threshold validation was performed using rainfall and landslide data from 1 January to 31 December 2018. During this period, a total of 52 landslide events occurred, out of which 80% of landslides (41) happened during the monsoon season. The threshold equation for Zone 1 was R TH = 150 − 0.24R 30ad , and its validation for the monsoon of 2018 is depicted in Figure 7. The threshold exceedance axis depicts the value of R TH with respect to R 30AD value for each day, wherein the positive values indicates landslide occurrence. This figure shows that a heavy measure of rainfall occurred, exceeding the threshold. The abrupt increase in the magnitude of daily rainfall or constant rise in 30 day antecedent rainfall is shown by the rise in the threshold curve.
During the validation period, the threshold was exceeded nine times, out of which seven times landslides occurred. No landslides were reported on 2 June and 7 June, even though the threshold was exceeded. This observation can be attributed to the fact that a landslide event does not generally happen with the increase in threshold curve, and sometimes happens a couple of days later as the result of a difference in pore pressure because of changes in the measure of antecedent rainfall. From 1 October to 30 December, there was no threshold exceedance, and no landslides occurred during that period. Similar validations were carried out for Zone 2 and Zone 3 threshold equations using 2018 rainfall and landslide data. These results indicate that the threshold model performed well for landslide forecasting in 2018.

Temporal Probability of Landslide Initiation
The temporal probability of landslide incident is determined as the product of annual exceedance probability (AEP) and probability of landslide occurrence [16]. AEP is defined as the probability of the threshold being exceeded in a given year [46], and is determined using the Poisson probability model defined as [16] P(N(t) = n) = e (−λt) n!
where N(t) represents the number of landslide incidences during time t and λ is landslide occurrence rate. The exceedance probability for time t is calculated as [13] P[N(t) where µ is the mean recurrence interval between subsequent landslides determined from landslide inventory. The determination of temporal probability was based on the following assumptions: (i) probability of landslide incidence is correlated to the probability of rainfall threshold being reached or exceeded, and (ii) landslides will not or will seldom occur when precipitation value is less than the threshold value [13,16,47].
The annual temporal probabilities for different zones of the study region are depicted in Figure  8, and their distribution along with the threshold equations is presented in Table 1. For Zone 1, the threshold value was exceeded 55 times during the simulation period, and out of these 55 cases, landslides were triggered in 16. The estimated probability P[L|(R > RT)] for Zone 1 was 0.29. Similarly, for Zone 2 and 3, the threshold value was exceeded 76 and 64 times in the 4 year period, leading to 31 and 21 landslides being triggered and contributing temporal probability values of 0.41 and 0.33, respectively. The probability of having one or more rainfall events in any given year varied from 0.29 to 0.41. The highest probability values were obtained for Zone 2, followed by Zone 3 and Zone 1. This variation was also observed in the number of landslide incidences for each zone. These precipitation events were also capable of triggering multiple landslides during the monsoon period.

Temporal Probability of Landslide Initiation
The temporal probability of landslide incident is determined as the product of annual exceedance probability (AEP) and probability of landslide occurrence [16]. AEP is defined as the probability of the threshold being exceeded in a given year [46], and is determined using the Poisson probability model defined as [16] P(N(t) = n) = e −λt (−λt) n n! (1) where N(t) represents the number of landslide incidences during time t and λ is landslide occurrence rate. The exceedance probability for time t is calculated as [13] P[N(t) where µ is the mean recurrence interval between subsequent landslides determined from landslide inventory. The determination of temporal probability was based on the following assumptions: (i) probability of landslide incidence is correlated to the probability of rainfall threshold being reached or exceeded, and (ii) landslides will not or will seldom occur when precipitation value is less than the threshold value [13,16,47].
The annual temporal probabilities for different zones of the study region are depicted in Figure 8, and their distribution along with the threshold equations is presented in Table 1. For Zone 1, the threshold value was exceeded 55 times during the simulation period, and out of these 55 cases, landslides were triggered in 16. The estimated probability P[L|(R > RT)] for Zone 1 was 0.29. Similarly, for Zone 2 and 3, the threshold value was exceeded 76 and 64 times in the 4 year period, leading to 31 and 21 landslides being triggered and contributing temporal probability values of 0.41 and 0.33, respectively. The probability of having one or more rainfall events in any given year varied from 0.29 to 0.41. The highest probability values were obtained for Zone 2, followed by Zone 3 and Zone 1. This variation was also observed in the number of landslide incidences for each zone. These precipitation events were also capable of triggering multiple landslides during the monsoon period.

Landslide Susceptibility Mapping
Landslide susceptibility can be defined as the probability of spatial occurrences of slope failures for a given set of geo-environmental conditions [48], and its determination is one of the crucial steps needed to understand identify potentially landslide-prone sections for any study region. Several studies around the world have been conducted towards the development of landslide susceptibility maps (LSM) using various methods [8]; however, there seems to be no consensus as to the best method for analysis [49]. Aleotti and Chowdhury [50] categorized LSM methods as either quantitative or qualitative. Qualitative models are mostly based on expert opinion, whereas quantitative models are data-driven, which makes them more reliable. The quantitative approaches include several kinds of techniques such, including statistical, deterministic, and other approaches [51][52][53]. In the case of statistical approaches, it is assumed that the parameters affecting landslide events in the past will be the same in future [54], and these analyses can be categorized into bivariate and multivariate [49]. In bivariate analysis, the factors affecting landslides are compared with landslide inventory data by providing weights based on landslide causative factors. The most frequently used methods in bivariate models are overlay, index-based, and weight-of-evidence analyses [8,51,55]. Bui et al. [6] performed a comparison between a bivariate approach (statistical index) and a multivariate approach (logistic regression) for Vietnam, and found equal forecasting capability. However, one of the main issues with the use of a quantitative approach is the assignment of weights to the landslide-affecting factors [56][57][58]. The use of GIS has been proven to be a powerful tool with which to validate the significance of factors, and it has been used for multi-criterion decision

Landslide Susceptibility Mapping
Landslide susceptibility can be defined as the probability of spatial occurrences of slope failures for a given set of geo-environmental conditions [48], and its determination is one of the crucial steps needed to understand identify potentially landslide-prone sections for any study region. Several studies around the world have been conducted towards the development of landslide susceptibility maps (LSM) using various methods [8]; however, there seems to be no consensus as to the best method for analysis [49]. Aleotti and Chowdhury [50] categorized LSM methods as either quantitative or qualitative. Qualitative models are mostly based on expert opinion, whereas quantitative models are data-driven, which makes them more reliable. The quantitative approaches include several kinds of techniques such, including statistical, deterministic, and other approaches [51][52][53]. In the case of statistical approaches, it is assumed that the parameters affecting landslide events in the past will be the same in future [54], and these analyses can be categorized into bivariate and multivariate [49]. In bivariate analysis, the factors affecting landslides are compared with landslide inventory data by providing weights based on landslide causative factors. The most frequently used methods in bivariate models are overlay, index-based, and weight-of-evidence analyses [8,51,55]. Bui et al. [6] performed a comparison between a bivariate approach (statistical index) and a multivariate approach (logistic regression) for Vietnam, and found equal forecasting capability. However, one of the main issues with the use of a quantitative approach is the assignment of weights to the landslide-affecting factors [56][57][58]. The use of GIS has been proven to be a powerful tool with which to validate the significance of factors, and it has been used for multi-criterion decision analysis [49,59]. The decision analysis technique combines primary-and secondary-level weights for every causative factor, where primary weights are similar to the bivariate approach and secondary weights are expert-opinion-based. For secondary weights, the analytic hierarchy process (AHP) has become popular and has been successfully applied for decision-making systems [60,61]. AHP uses a pairwise relative comparison between every landslide-causative factor. Generally, AHP consists of five key steps: (a) simplify the decision process into its component factors, (b) distribute the factors in a hierarchy process, (c) allocate numerical values to analyze the relative significance of each factor, (d) compose a comparison matrix, and (e) provide weights to every factor by calculating normalized principal eigenvectors [62].
To determine susceptibility, a variety of factors responsible for landslides in the study region were considered. Parameter selection depends on various factors, such as landslide type, data availability and reliability, and adopted methods [63]. For the present study, we used eight landslide-conditioning factors based on the characteristics of the area and prepared from various data sources (Table 2).  The above-mentioned thematic layers were combined by using a weight-of-factors approach determined by AHP to develop the landslide susceptibility map. The use of AHP to develop landslide susceptibility maps has been successfully applied in various regions [61,66,67]. The weights required to carry out AHP were calculated by performing pairwise comparisons for each landslide factor and assigning values from 1 to 9 [63,[68][69][70]. Table 3 shows the pairwise comparison and priority calculation, along with rankings of all indicators. These values were based on an expert's opinion and were placed in n × n matrix, where n is the number of factors.   The AHP reduced the inconsistencies formed due to the subjectivity of different experts' opinions by computing a consistency index (CI) and consistency ratio (CR), which were determined by where λ max represents the largest Eigenvector of the matrix and n represents the total causative factors (order of the matrix) used in the generation of the LSM. RI (random index) is the average value of CI for a randomly generated pairwise matrix and can be accepted only when CR values are less than 10% [71]. For the present study, the average consistency index was estimated for a sample size of 500 and its value was 0.039 (3.9%), which was considered acceptable. Several authors have calculated and estimated different RIs based on various simulation methods and the total number of matrices involved in the process (  The values of CR cannot be negative and can attain a maximum value of 0.3. Values of CR less than 0.1 are considered acceptable; if this is not achieved, new attempts are made until the value is acceptable [80]. However, the values of CR are dependent on the analysis type and the number of criteria being considered. In some cases, CR > 0.1 may not be considered critical, and values of CR ranging from 0.15 to 0.3 can also be considered acceptable. A matrix will be considered consistent, according to Saaty [71], if λmax < n + 0.1((λmax) − n) Finally, all the landslide-causative factors and classes were integrated by a method of weighted overlay in ArcGIS to generate the landslide susceptibility index (LSI). LSI = k=1 n W k W jk (6) where W k is the weight of the causative factor, W jk is the rank value for factor class j of causative factor k, and n represents the total causative factors selected. Ranks of criteria were calculated based on priority or weight values. The highest priority considered was Rank 1, while the lowest priority considered was the last rank in the AHP.

Model Development and Validation
The landslide inventory data was randomly categorized into two datasets, i.e., training (70%) and testing (30%). The landslide susceptibility map based on 8 causative factors using AHP is depicted in Figure 10. The map was divided into 4 classes: very low (0-0.24), low (0.25-0.49), moderate (0.5-0.74) and high (0.75-1) according to natural breaks to define the class intervals in the susceptibility map. From total of 242 landslides, more than 78.1% falls under the high zone. Whereas 11.98% falls under the moderate zone and combined 9.92% falls under the low and very low zones (Table 5).
where Wk is the weight of the causative factor, Wjk is the rank value for factor class j of causative factor k, and n represents the total causative factors selected. Ranks of criteria were calculated based on priority or weight values. The highest priority considered was Rank 1, while the lowest priority considered was the last rank in the AHP.

Model Development and Validation
The landslide inventory data was randomly categorized into two datasets, i.e., training (70%) and testing (30%). The landslide susceptibility map based on 8 causative factors using AHP is depicted in Figure 10. The map was divided into 4 classes: very low (0-0.24), low (0.25-0.49), moderate (0.5-0.74) and high (0.75-1) according to natural breaks to define the class intervals in the susceptibility map. From total of 242 landslides, more than 78.1% falls under the high zone. Whereas 11.98% falls under the moderate zone and combined 9.92% falls under the low and very low zones (Table 5).  Validation of the susceptibility maps was performed for randomly selected data from the inventory data using receiver operating characteristics (ROC). This is an effective way to analyze the quality of predictive techniques [81]. The ROC curve is plotted between true positive rate (sensitivity) on the Y axis against false positive rate (specificity) on the X axis. The terms "sensitivity" and "specificity", which are used to plot ROC curves, are defined as follows.
where true positive (TP) is the number of actual landslides predicted correctly, and true negative (TN) is the total number of non-occurring landslides predicted correctly. False positive (FP) is the number of actual landslides inaccurately predicted as non-occurring landslides, and false negative (FN) is the number of non-occurring landslides inaccurately predicted as actual landslides. [82]. The area under the ROC curve (AUC) was also used to determine the quality of the prediction by analyzing the model's ability to forecast the occurrence or nonoccurrence of predefined events [83]. The results of the success rate curve of the AHP model had an AUC of 0.798, corresponding to a prediction accuracy of 79.8% ( Figure 11). Validation of the susceptibility maps was performed for randomly selected data from the inventory data using receiver operating characteristics (ROC). This is an effective way to analyze the quality of predictive techniques [81]. The ROC curve is plotted between true positive rate (sensitivity) on the Y axis against false positive rate (specificity) on the X axis. The terms "sensitivity" and "specificity", which are used to plot ROC curves, are defined as follows.

Sensitivity = TP TP + FN (7)
Specificity = TN FP + TN (8) where true positive (TP) is the number of actual landslides predicted correctly, and true negative (TN) is the total number of non-occurring landslides predicted correctly. False positive (FP) is the number of actual landslides inaccurately predicted as non-occurring landslides, and false negative (FN) is the number of non-occurring landslides inaccurately predicted as actual landslides. [82]. The area under the ROC curve (AUC) was also used to determine the quality of the prediction by analyzing the model's ability to forecast the occurrence or nonoccurrence of predefined events [83]. The results of the success rate curve of the AHP model had an AUC of 0.798, corresponding to a prediction accuracy of 79.8% ( Figure 11). The results of pairwise comparison, priority estimation, and ranking of all the criteria could be applied to other study areas for susceptibility assessment. Several high-resolution satellite image datasets are required to better understand the locations and perform these assessments. In cases of unavailability of high-resolution satellite images, Google Earth images could be useful for preparation of some indicators. The use of Google Earth images could also be helpful for accurate The results of pairwise comparison, priority estimation, and ranking of all the criteria could be applied to other study areas for susceptibility assessment. Several high-resolution satellite image datasets are required to better understand the locations and perform these assessments. In cases of unavailability of high-resolution satellite images, Google Earth images could be useful for preparation of some indicators. The use of Google Earth images could also be helpful for accurate identification of landslide locations and conduction of extensive field studies. All the requirements for more accurate analysis depend on study location tectonics, data availability, and proper methods. Therefore, based on local geology and tectonic conditions, these results could be transferable and applicable in other locations in both small-and large-scale areas.

Conclusions
Landslides are the most frequently occurring natural hazards, especially in the Himalayan regions, which suffer from heavy monsoonal rainfall and subsequent landslides. In this study, the temporal probability of landslide events was determined using rainfall and landslide data from 2014-2017 along Samdrup Jongkhar-Trashigang highway in East Bhutan. The highway was divided into three zones based on land use, topography, and rain gauge coverage for the determination of temporal probability. Thereafter, a landslide susceptibility map was developed using AHP. The results of temporal probability were validated with landslide event dataset of 2018 to understand the applicability of the model. The conclusions from the study can be summarized as follows.
(1) Threshold determination was performed considering the antecedent rainfall duration approach.
The selection of antecedent rainfall was based on previous studies which have highlighted its significance, especially in the Himalayan region and in Bhutan. A trial-and-error approach was adopted to determine the number of antecedent days required for landslide initiation, and a 30 day antecedent rainfall was adopted for threshold determination. (2) The temporal probability of landslide occurrence was determined based on the Poisson model, and the validation results based on 2018 data revealed that the model performed well. (3) The landslide susceptibility map of the region was developed using AHP classified into four categories (Very Low, Low, Moderate, and High), and the results showed that 78% of the region falls under into high-susceptibility zone. The performance of the model was assessed using the area under ROC and an accuracy of 79.8% was achieved. However, due to the dynamic nature of land use and rainfall patterns, the susceptibility map will need to be updated from time to time.
The present study on rainfall threshold estimation and the development of a susceptibility map for eastern Bhutan along Samdrup Jongkhar-Trashigang highway is an important study in the context of the Bhutan Himalayas, for which study on both these aspects is lacking. The current study can be regarded as a preliminary step towards risk management, which could be supported by conducting a hazard and vulnerability assessment of the region. The temporal probabilities determined can be integrated with the susceptibility map to obtain a landslide hazard map. However, the results might be improved by increasing the number of landslide events and using precipitation data with a higher temporal resolution. The results from the present study also could prove helpful to civic authorities in identifying key sections of the road which are most vulnerable to landslides, and undertaking strict measures to prevent slope failures, strengthening the transportation network and saving human lives.