Determination of the Probabilities of Landslide Events—A Case Study of Bhutan

Landslides have been and are prominent and devastating natural disasters in Bhutan due to its orography and intense monsoonal rainfall. The damage caused by landslides is huge, causing significant loss of lives, damage to infrastructure and loss of agricultural land. Several methods have been developed to understand the relationship between rainfall and landslide incidences. The most common method to understand the relationship is by defining thresholds using empirical methods which are expressed in either intensity-duration or event rainfall-duration terms. However, such thresholds determine the results in a binary form which may not be useful for landslide cases. Apart from defining thresholds, it is significant to validate the results. The article attempts to address both these issues by adopting a probabilistic approach and validating the results. The region of interest is the Chukha region located along the Phuentsholing-Thimphu Highway, which is a significant trade route between neighbouring countries and the national capital Thimphu. In the present study, probabilities are determined by Bayes’ theorem considering rainfall and landslide data from 2004 to 2014. Singular (rainfall intensity, rainfall duration and event rainfall) along with a combination (rainfall intensity and rainfall duration) of precipitation parameters were considered to determine the probabilities for landslide events. A sensitivity analysis was performed to verify the determined probabilities. The results depict that a combination of rainfall parameters is a better indicator to forecast landslides as compared to single rainfall parameter. Finally, the probabilities are validated using landslide records for 2015 using a threat score. The validation signifies that the probabilities can be used as the first line of action for an operational landslide warning system.


Introduction
Landslides are the most common and one of the prominent natural hazards in steep terrain affecting human lives and property, and often leading to loss of lives. In [1] the study of a global database of landslide occurrences between 2004-2016 showed that 75% of landslides occurred in Asia, with significant occurrences in the Himalayan arc. The study also showed that the majority of the landslides are shallow landslides and are triggered by rainfall. The high or low intensity reinforced with continuous occurrence of rainfall leads to an increase in pore water pressure due to the constant saturation of soil mass, which leads to landslides [2]. The incidence of landslides is likely to increase in the future due to a growing population, increased construction activities, and exploitation of natural resources. To understand the phenomena of shallow landslides, several models were developed and used in several parts of the world. These models can be broadly categorised as empirical or physical models. Physical models require the use of several parameters involving slope-stability modelling and acquiring the data required to carry out such an analysis are usually challenging [3]. Therefore, the most used technique is the empirical model which uses the observed rainfall and landslide data which are analyzed statistically without considering the physical processes governing the landslide [4][5][6].
The key in utilising such a model is the availability and quality of data which reflects in the results [7]. The performance of the threshold-based models can decrease significantly with a small change in landslide data [8]. References [9][10][11] highlighted that definition of rainfall event and information on landslide data had a significant impact on threshold estimation which often leads to a high number of false alarms for application in an early warning system. Recent developments have led to reproducible conditions for landslide occurrences in the region of interest [7,12,13]. However, one major drawback for using empirical models is the outcome of an absolute value either in terms of rainfall intensity or cumulative event rainfall [14]. Such a result may not always help in setting up an early warning system since a small variation in calculation can lead to catastrophes. To overcome the use of empirical models, probabilistic model using Bayes' theorem was developed by [14] is an efficient choice and the same was applied to the Emilia-Romagna region in Italy. Thereafter, similar work was carried out in several places (e.g., Ha Giang region, Vietnam [15], Sierra Norte De Puebla, Mexico [16], Kalimpong, India [17,18], Chibo, India [19]. The use of probabilistic techniques has also been carried out for severe precipitation conditions [20,21].
The key to understand the thresholds or any analysis is its validation, [6] reported that out of 115 rainfall thresholds determined between 2008-2016, only 69 provided validation of their work. For validation work, only 33% of the work was carried out using an independent dataset. The most common technique for validation is skill score [8,19] or by comparing two different threshold models [22].
The Bhutan Himalayas region has been adversely affected by landslides which are caused by heavy monsoon precipitation. In the case of Bhutan, the literature available on rainfall-induced landslides is very sparse, with only one work showing the determination of event rainfall-rainfall duration (ED) thresholds [13]. The country also suffers from the absence of an operational early warning system and therefore understanding the relation using different models is very critical in developing such an early warning system [13]. This paper determines and validates the probabilities of landslide incidents using the Bayes' theorem by utilizing the available landslide and rainfall data for the Chukha region. The validation of the thresholds has been determined using a threat score which analyses the percentage of correctly predicted landslide events. The probabilities were calculated using the dataset for 11 years (2004-2014) and validated using landslide data of 2015.

Study Area
Bhutan lies in the eastern part of Himalaya and is surrounded by India to the east, west and south and by the Tibetan plateau to the north. The elevation of the country varies from 100 m to 7500 m [23] with elevations in Chukha region ranging from 1000 m to 4200 m. Figure 1 shows the present study area-the Chukha region. A part of the Phuentsholing-Thimphu highway (also called Asian Highway), which connects Thimphu, the capital of Bhutan, with neighbouring countries, lies in the study region. Said highway is a major trade route with neighbouring countries and is frequently blocked by landslide debris, causing severe economic damage every year.
The landslides in the region are primarily caused due to heavy monsoons and are shallow in nature. The increase in landslide events is frequently activated by toe cutting of the slope and blasting activities for infrastructure development. The increase in anthropogenic activities has escalated deforestation leading to slope instability. The slope failures along the roads usually block the road leading to huge logistic issues and economic loss.
The region comprises of active sedimentary and metasedimentary rocks, gneiss, schist, quartzite, and limestone and is underlain mostly with schistose rocks. The soils in the area consist of weaker phyllites which make the soil texture very fine and the slopes unstable. Geologically, the region is surrounded by major tectonostratigraphic units and structures are Shivalik ranges, Main Boundary Thrust (MBT), the Lesser Himalayan Sequence (LHS), MCT, HHC, and STD [24]. The major parts of the area are underlain by Baxa formation group where the thrust fault, strike and dip of bedding and foliation were more concentrated.  The slope of the region varies from 15°-75° with a majority of the area (35%) lying in the moderate steep slope (30°-45°). The presence of weak geology and weathered rocks makes the region highly susceptible to landslides. Figure 2    The slope of the region varies from 15 • -75 • with a majority of the area (35%) lying in the moderate steep slope (30 • -45 • ). The presence of weak geology and weathered rocks makes the region highly susceptible to landslides. Figure 2 shows typical damages caused by landslides in the area.

Methodology
The methodology involves the use of Bayes' theorem which uses either single or multiple rainfall characteristics to determine one-dimensional and two-dimensional probabilities respectively. The technique can be applied for degrees of variables depending on the quality and quantity of rainfall and landslide data. A detailed explanation of the method along with a solved example has been presented in [14].

One-Dimensional Bayesian Probability
This type of probability approach involves the use of a single rainfall parameter to determine the landslide probability and is defined as: where P(Y|X) = precipitation probability of magnitude Y for landslide occurrence, P(X) = landslide probability regardless of precipitation event of magnitude X, P(Y) = precipitation probability of magnitude X, regardless of landslide occurrence and P(X|Y) = landslide probability for precipitation event of magnitude Y

Methodology
The methodology involves the use of Bayes' theorem which uses either single or multiple rainfall characteristics to determine one-dimensional and two-dimensional probabilities respectively. The technique can be applied for degrees of variables depending on the quality and quantity of rainfall and landslide data. A detailed explanation of the method along with a solved example has been presented in [14].

One-Dimensional Bayesian Probability
This type of probability approach involves the use of a single rainfall parameter to determine the landslide probability and is defined as: where P(Y|X) = precipitation probability of magnitude Y for landslide occurrence, P(X) = landslide probability regardless of precipitation event of magnitude X, P(Y) = precipitation probability of magnitude X, regardless of landslide occurrence and P(X|Y) = landslide probability for precipitation event of magnitude Y.
If the number of rainfall events is M R, number of landslide events being M A , number of precipitation events of degree B be M B and the number of precipitation events leading to slide initiation be M (B|A) , then Equation (1) reduces to the computation of the following frequencies: The parameters being used for the analysis would depend on the key reasons for slide initiation for a region. The variables used in the present are cumulative event rainfall, rainfall duration, and rainfall intensity.

Two-Dimensional Bayesian Probability
The two-dimensional Bayesian probability determines the conditional probability of landslide event due to the common occurrence of two rainfall parameters: where Y, Z denotes the combined probability of having a certain or range of value of any two variables. If, Y equals rainfall intensity and Z is rainfall duration, the likelihood of landslide occurrence due to a rainfall event of given duration and intensity is expressed using Equation (5). The model allows any set of rainfall parameters to be used to calculate thresholds and its significance can be compared with prior probability [14]. Also, it can also be used for n-variables by modifying the equation accordingly like understanding the effect of a combination of rainfall intensity, duration, event rainfall and antecedent rainfall on landslide incidence.

Data
For the present study, the daily rainfall data was collected from rain gauge installed at Chuka and maintained by the National Center for Hydrology and Meteorology of the Royal Government of Bhutan. The daily rainfall data were collected for 1 January 2004 to 31 December 2014. The variation of the rainfall during the entire duration has been depicted in a box and whisker plot ( Figure 3). The average rainfall during 2004-2014 was 1663.4 mm with maximum rainfall of 2926.6 mm occurring in 2012. The monsoonal rain (June-August) contributed about 80% of annual rainfall.
The landslide data for the study was provided by the Border Roads Organization, (Project DANTAK), Government of India, which maintains the Phuentsholing-Thimphu highway. The temporal distribution of landslides indicates that majority of the landslides initiate during the monsoonal period, and the spatial distribution depicts that landslide events occur along the highway. The spatial distribution of landslides is depicted in Figure 1b and involved materials are quartzite, phyllite, gneiss. Landslide event has been identified as only those landslides mapped as single points which were triggered due to rainfall. Also, if multiple events occurred on a day it was counted as a single event. The effect of the spatial distribution of landslides was considered by drawing a buffer radius of 5 km around the rain gauge [7]. A large search radius is usually required when the study area suffers from low rain gauge density [13]. A total of 123 landslide events occurred during 2004-2014. Considering the above mentioned factors, the number of landslide events was reduced to 66.
For the present study, the daily rainfall data was collected from rain gauge installed at Chuka and maintained by the National Center for Hydrology and Meteorology of the Royal Government of Bhutan. The daily rainfall data were collected for 1 January 2004 to 31 December 2014. The variation of the rainfall during the entire duration has been depicted in a box and whisker plot (Figure 3). The average rainfall during 2004-2014 was 1663.4 mm with maximum rainfall of 2926.6 mm occurring in 2012. The monsoonal rain (June-August) contributed about 80% of annual rainfall.

Results
The key to determine any thresholds is to define rainfall and landslide event. In this study, rainfall event is defined as the total number of consecutive days of rainfall. The rainfall intensity (mm/day) is determined by dividing the corresponding total rainfall (mm) by the rainfall duration (days). The landslide event was described as the single rainfall triggered slide activity. The landslide event for analysis includes the spatiotemporal distribution, i.e., dates and coordinates of each landslide. The total number of rainfall and landslide events determined was 480 and 66, respectively. The results of one-dimensional Bayesian probability are depicted in Figure 4a-f.
The results show that probability reaches the highest value of 0.5 for cumulative event rainfall of 300-350 mm and a probability of 0.33 for rainfall duration of 9 days. In case of rainfall intensity, the likelihood is 0.5 for the intensity of 80-90 mm/day and 0.227 for the intensity of 30-40 mm/day. This shows that heavy rainfall as well as slow and continuous rainfall can lead to landslide incidence in the region. The landslide data for the study was provided by the Border Roads Organization, (Project DANTAK), Government of India, which maintains the Phuentsholing-Thimphu highway. The temporal distribution of landslides indicates that majority of the landslides initiate during the monsoonal period, and the spatial distribution depicts that landslide events occur along the highway. The spatial distribution of landslides is depicted in Figure 1b and involved materials are quartzite, phyllite, gneiss. Landslide event has been identified as only those landslides mapped as single points which were triggered due to rainfall. Also, if multiple events occurred on a day it was counted as a single event. The effect of the spatial distribution of landslides was considered by drawing a buffer radius of 5 km around the rain gauge [7]. A large search radius is usually required when the study area suffers from low rain gauge density [13]. A total of 123 landslide events occurred during 2004-2014. Considering the above mentioned factors, the number of landslide events was reduced to 66.

Results
The key to determine any thresholds is to define rainfall and landslide event. In this study, rainfall event is defined as the total number of consecutive days of rainfall. The rainfall intensity (mm/day) is determined by dividing the corresponding total rainfall (mm) by the rainfall duration (days). The landslide event was described as the single rainfall triggered slide activity. The landslide event for analysis includes the spatiotemporal distribution, i.e., dates and coordinates of each landslide. The total number of rainfall and landslide events determined was 480 and 66, respectively. The results of one-dimensional Bayesian probability are depicted in Figure 4a-f.
The results show that probability reaches the highest value of 0.5 for cumulative event rainfall of 300-350 mm and a probability of 0.33 for rainfall duration of 9 days. In case of rainfall intensity, the likelihood is 0.5 for the intensity of 80-90 mm/day and 0.227 for the intensity of 30-40 mm/day. This shows that heavy rainfall as well as slow and continuous rainfall can lead to landslide incidence in the region.   The probability values are also hugely affected by the number of landslide events in each bin. There can be a huge difference in the values due to slight variation in the data points. For the present study, the number of landslide events for rainfall intensity in the bin 10-20 mm/day is 26, however, the number of points for 80-90 mm/day is only 5. This suggests that validation of the probabilities needs to be carried out and the results need to be enhanced with the addition of more data.
The result for two-dimensional Bayesian probability for landslide occurrences in Chukha has been depicted in Figure 5. The results indicate that maximum probability reaches 1 for the intensity of 30-40 mm/day with a duration of 7 days as well as the intensity of 40-50 mm/day with a duration of 6 days. However, the value of probability reaching an extreme amount can be attributed to low sample size and a slight variation in sample size can immensely affect the data. [18] determined probabilities for Kalimpong region which is roughly 100 km from the study area depicted that a precipitation event with intensity higher than 30 mm/day for 3 days represents the highest probability of 0.67. A sensitivity analysis of the calculated probabilities was performed to determine the variation of the results with a change in rainfall event parameters [19]. The input values of the parameters were varied and the changes in values of probabilities were ascertained. The variations include the definition of rainfall event definition and the period of analysis. The rainfall event for the probability determination was defined as the number of consecutive days of rainfall. For sensitivity  The probability values are also hugely affected by the number of landslide events in each bin. There can be a huge difference in the values due to slight variation in the data points. For the present study, the number of landslide events for rainfall intensity in the bin 10-20 mm/day is 26, however, the number of points for 80-90 mm/day is only 5. This suggests that validation of the probabilities needs to be carried out and the results need to be enhanced with the addition of more data.
The result for two-dimensional Bayesian probability for landslide occurrences in Chukha has been depicted in Figure 5. The results indicate that maximum probability reaches 1 for the intensity of 30-40 mm/day with a duration of 7 days as well as the intensity of 40-50 mm/day with a duration of 6 days. However, the value of probability reaching an extreme amount can be attributed to low sample size and a slight variation in sample size can immensely affect the data. [18] determined probabilities for Kalimpong region which is roughly 100 km from the study area depicted that a precipitation event with intensity higher than 30 mm/day for 3 days represents the highest probability of 0.67. A sensitivity analysis of the calculated probabilities was performed to determine the variation of the results with a change in rainfall event parameters [19]. The input values of the parameters were varied and the changes in values of probabilities were ascertained. The variations include the definition of rainfall event definition and the period of analysis. The rainfall event for the probability determination was defined as the number of consecutive days of rainfall. For sensitivity analysis, we varied the minimum daily rainfall in intervals of 1mm and the corresponding probabilities were calculated which depicted similar trend as shown in Figures 4 and 5. The probabilities were recalculated by decreasing the analysis period and the changes in probabilities were in accordance to the number of events being considered. The results showed the increase in probability values for certain bins due to few numbers of landslide and rainfall events. The results show that the probabilities varied between 12%-22% for various bins of different rainfall parameters, however the peak values of landslide probability were in the same rainfall class. The analysis also suggested that one-dimensional results are more sensitive compared to two-dimensional Bayesian analysis. analysis, we varied the minimum daily rainfall in intervals of 1mm and the corresponding probabilities were calculated which depicted similar trend as shown in Figures 4 and 5. The probabilities were recalculated by decreasing the analysis period and the changes in probabilities were in accordance to the number of events being considered. The results showed the increase in probability values for certain bins due to few numbers of landslide and rainfall events. The results show that the probabilities varied between 12%-22% for various bins of different rainfall parameters, however the peak values of landslide probability were in the same rainfall class. The analysis also suggested that one-dimensional results are more sensitive compared to two-dimensional Bayesian analysis. The validation of the probabilities was determined using the landslide record for 2015. During this period, there were 38 rainfall events and five landslide events with cumulative rainfall being 2698.83 mm ( Figure 6). The validation of the probabilities was determined using the landslide record for 2015. During this period, there were 38 rainfall events and five landslide events with cumulative rainfall being 2698.83 mm ( Figure 6). Reference [13] defined regional based thresholds using an algorithm-based approach and determined that 69.7 mm of total rainfall for 24 h can cause landslides. Considering this value, the probability in terms of event rainfall is 0.14. To evaluate the performance of any model, threat score is considered as a useful metric [21]. Threat score is defined as: where TS is threat score, TP is true positive, FN and FP being false negatives and false positive respectively [25]. For the present study the various values of the Equation (6) Reference [13] defined regional based thresholds using an algorithm-based approach and determined that 69.7 mm of total rainfall for 24 h can cause landslides. Considering this value, the probability in terms of event rainfall is 0.14. To evaluate the performance of any model, threat score is considered as a useful metric [21]. Threat score is defined as: where TS is threat score, TP is true positive, FN and FP being false negatives and false positive respectively [25]. For the present study the various values of the Equation (6), TP = 4, FP = 1, FN = 2, which makes TS = 0.57. The result shows that the model has a good potential to be used for early warning system and can be improved with the further addition of data.

Conclusions
Landslides can be considered as one of the most devastating natural disasters due to their reoccurrence, economic damages and loss of human lives. Most landslides are shallow and caused due to rainfall which can be torrential downpours or slow and continuous. Therefore, it is imperative to understand the relation between rainfall and landslide incidences. Significant numbers of landslides are occurring throughout the Himalaya and as a part of it, Bhutan is not left out. The studies on understanding the relation of landslide incidences with rainfall have been rarely studied in this part of the world, so the present paper calculates the results of a probabilistic approach using Bayes' theorem using data between 2004-2014 to understand the effect of rainfall parameters for landslide initiation. The results were determined using a single along with a combination of rainfall parameters. The results were further validated using the landslide data of 2015 for the event rainfall parameter. The conclusions from the study are: • The use of a probabilistic approach can be a better approach than empirical thresholds as the latter provides a single value of a specific rainfall parameter for landslide incidences.

•
The use of two-dimensional probability for determining probabilities for landslide events is better as compared to one-dimensional as the latter depicts that a single rainfall parameter may not be a significant factor to trigger landslides.

•
The validation of the thresholds for event rainfall parameter depicts that the model has an accuracy of 57%. However, with the addition of more landslide records and temporal rainfall data, the accuracy will improve. The use of such technique would help in setting up an operational early warning system and help in landslide mitigation.
Author Contributions: R.S. analyzed the data and wrote the article, K.D. collected the data and contributed in writing the article.