Estimating Rainfall Thresholds for Landslide Occurrence in the Bhutan Himalayas

Abstract

Consistently over the years, particularly during monsoon seasons, landslides and related geohazards in Bhutan are causing enormous damage to human lives, property, and road networks. The determination of thresholds for rainfall triggered landslides is one of the most effective methods to develop an early warning system. Such thresholds are determined using a variety of rainfall parameters and have been successfully calculated for various regions of the world at different scales. Such thresholds can be used to forecast landslide events which could help in issuing an alert to civic authorities. A comprehensive study on the determination of rainfall thresholds characterizing landslide events for Bhutan is lacking. This paper focuses on defining event rainfall–duration thresholds for Chukha Dzongkhag, situated in south-west Bhutan. The study area is chosen due to the increase in frequency of landslides during monsoon along Phuentsholing-Thimphu highway, which passes through it and this highway is a major trade route of the country with the rest of the world. The present threshold method revolves around the use of a power law equation to determine event rainfall–duration thresholds. The thresholds have been established using available rainfall and landslide data for 2004–2014. The calculated threshold relationship is fitted to the lower boundary of the rainfall conditions leading to landslides and plotted in logarithmic coordinates. The results show that a rainfall event of 24 h with a cumulated rainfall of 53 mm can cause landslides. Later on, the outcome of antecedent rainfall varying from 3–30 days was also analysed to understand its effect on landslide incidences based on cumulative event rainfall. It is also observed that a minimum 10-day antecedent rainfall of 88 mm and a 20-day antecedent rainfall of 142 mm is required for landslide occurrence in the area. The thresholds presented can be improved with the availability of hourly rainfall data and the addition of more landslide data. These can also be used as an early warning system especially along the Phuentsholing–Thimphu Highway to prevent any disruptions of trade.

1. Introduction

Rainfall or earthquake triggered landslides are common in some parts of the world, causing loss of human life and property [1]. A global dataset of landslide disasters [2] showed that three-quarters of all landslide events occurred in the Himalayan arc between 2004 and 2016. Bhutan is one of the highly susceptible landslide zones in the Himalayan region [3]. The majority of landslides in Bhutan are initiated due to heavy monsoon precipitation and aggravated due to the increase in human activities. The increasing number of landslide events in Bhutan can be attributed to complex geological conditions, steep slopes, climate change, type of soil, and tectonic activity. Landslides in Bhutan mostly occur during the monsoons during which the torrential rainfall leads to several flash floods and landslide triggering, cutting off many parts of south Bhutan from the rest.

The relationship between the amount of rainfall associated with landslide occurrences is generally studied using either an empirical or physical based approach [4,5,6,7]. Physical process models are based on numerical models which study the relationship between rainfall, pore water pressure, soil type, and volumetric water content that can lead to slope instability. Such a study is usually site specific due to variation in soil properties. It is a challenge to extend this approach to large areas, as the extensive data that is required are usually not available, and their use for an early warning system is either experimental or prototype based [8,9,10]. On the other hand, empirical methods study the landslides that are caused by rainfall events—both massive downpour that triggers instantaneous landslides and the low but continuous antecedent rain that destabilizes the slope and triggers the landslide. Although, there are many factors like rainfall, earthquake, geology, soil type etc. involved for landslide triggering, in the present study, precipitation rates have been considered as this is the primary cause of several changes in soil properties, pressure variations, etc. The minimum quantity of precipitation requisite for landslide occurrences known as thresholds can be determined using empirical models. The limit is defined by lower-bound lines to the precipitation conditions causing landslides and plotted in Cartesian, semi-logarithmic, or logarithmic coordinates [11]. Contingent upon the kind of available rainfall data, empirical thresholds can be summarized as follows: (1) thresholds which combine rainfall data obtained from specific rainfall events [12], (2) thresholds involving antecedent parameters [6], and (3) alternating thresholds, like hydrological thresholds [13]. Therefore, several works can be found depicting rainfall thresholds based on empirical techniques [4,7,14,15,16,17,18]. The present study highlights the importance of antecedent rainfall along with the determination of cumulative event-rainfall–duration thresholds for an operational early warning system. A recent review on rainfall thresholds [19] showed that the thresholds could be used to predict landslide events at various geographical extents and also in a broad spectrum of physical settings [10]. The study also found that, for setting up an early warning system using empirical rainfall thresholds, various factors needs to be taken care of: (i) collection of reliable and large rainfall and landslide datasets, (ii) selecting threshold parameters depending on landslide characteristics and precipitation data, (iii) defining the events and using an objective and standardized methodology, (iv) validation of the thresholds determined. The recent development is on defining objective thresholds using semi-automated algorithms [20,21,22,23].

In the context of Bhutan, it is difficult to study the relationship using the physical based approach as the data required is not available. Therefore, this study is lying on an empirically based approach, which defines thresholds using available records of daily rainfall and landslide data in the period 2004–2014. The available scientific literature for landslides in Bhutan is very sparse with no operating landslide early warning system [24,25]. Only one study has been carried out to determine site-specific rainfall thresholds [25], which determined thresholds using an algorithm-based approach, CTRL-T (Calculation of Thresholds for Rainfall-Induced Landslides Tool). However, more extensive work can be found for the Indian Himalayas [11,26,27,28] and Nepal Himalayas [29]. The study also discusses the significance of antecedent rainfall and the possibility to use the output as a decision tool for landslide determination. The effect of antecedent rainfall on landslides was analysed for various time durations: 3, 7, 10, 20, and 30 days. Further, the analysis has been validated using a statistical indicator, i.e., threat score for 2015 rainfall data. The paper is presented in five sections; the first two sections are listed as ‘Study Area’ and ‘Rainfall and Landslide data’. The third section is listed as ‘Methodology’ where explanation on threshold calculation and the uncertainty associated with it are described. The fourth heading is marked as ‘Results and Discussions’ summarizing threshold determination along with understanding the importance of antecedent rainfall. The findings of the present study are summarized in section five ‘Conclusions’.

2. Study Area

Bhutan is situated in the eastern section of the Himalaya with an area of 38,394 km² (Figure 1). The nation is enclosed with the Tibetan Plateau in the north and Indian states Arunachal Pradesh in the east, Bengal and Assam in the south, and the Darjeeling–Sikkim Himalayas in the west. Bhutan is divided into 20 Dzongkhags and has elevations varying from 150 m to 7570 m [30]. The elevation of the present study area Chukha Dzongkhag ranges from 1000 m to 4200 m. Most of the people in the region depend on agriculture and livestock for their livelihood.

Geologically, the Chukha region belongs to the Lesser Himalayan formation which includes a wide variety of sedimentary and low-grade metamorphic rocks. The study area is mainly comprised of tectonically active metasedimentary rocks such as phyllite, schist, quartzite, and limestone. The northern part belongs to the Higher Himalayan crystalline rocks comprising mainly of garnetiferous mica-shist, quartzite, and gneiss. The soils in the region are mostly moderate to high weathered and are comprised of weaker fracture phyllites, which make the soil texture very fine, and slopes are very unstable [25]. Based on the practical experiments, the soil can be classified as poorly graded sand (SP) mixed with pebbles. The geology of the region described in [31]. The tectonic setting of the region is very similar to Nepal and Indian Himalayas [25]. The morphology of the unstable slopes leading to failures are complex and managed by several factors including lithology and rock type.

The region receives yearly rainfall of 4000–6000 mm with heavy bursts reaching up to 800 mm/day [18]. More than three-quarters of the annual rainfall occurs during the monsoon period of June to September [25]. Figure 1b depicts the landslide occurrences between 2004–2014 and the rain gauge coordinates used for the analysis. The monsoon rainfall leads to erosion of rocks causing widespread slope instability and mass movements. The increase in anthropogenic activities has escalated deforestation leading to slope instability. The landslide typology in the region can be described as rock fall, rockslide, debris flow, debris slide, and earth slide [25]. The major effects due to this monsoonal rainfall of this region are roadblocks. As mentioned earlier, the Phuentsholing-Thimphu highway (also known as Asian Highway) situated in Chukha Dzongkhag is one of the major road connections of Bhutan, which is key to the trade, passes through these landslides affected areas. The landslide events along the highway pose serious logistic problems as it affects the economy of the region. Figure 2a,b illustrates the damage caused by landslides along the highway.

3. Rainfall and Landslide Data

The daily rainfall data were collected from the three rain gauges from Malbase, Gedu, and Chukha. These three rain gauges are maintained by the National Center for Hydrology and Meteorology, Royal Government of Bhutan. The presence of a properly distributed rain gauge network is essential for determining rainfall thresholds. Figure 3 depicts the cumulative rainfall distribution for the study period (2004–2014). The rainfall data collected for 2004–2014 shows that about 78% of the annual rainfall falls during the monsoon (June–September) with pre-monsoon (March–May) and post-monsoon (October–December) contributing 18% and 4%, respectively.

The landslide data were provided by project Dantak of Border Road Organization (BRO), of the Govt. of India. The landslide locations are confined along the Phuentsholing–Thimphu highway, most prominent trade route, which connects the capital city Thimphu with neighbouring countries by road. The total number of landslides during the study period was 248, out of which 63, 81, and 105 landslides occurred in Malbase, Gedu, and Chukha regions respectively. The landslides in the region are mostly shallow with depths ranging from few decimetres to metres [25]. The catalogue of landslides included coordinates and date of failures. Figure 4 depicts the year-wise landslide occurrences.

4. Method

This section explains the procedure to determine thresholds. The determination of thresholds using empirical methods can be divided into two primary steps. The first step is the collection of rainfall and landslide data, and the next step is to apply empirical models to determine the thresholds. An event rainfall is determined by the number of consecutive days of rainfall before the landslide incident. This helps in calculating the number of days of rainfall before landslide and the total rainfall. Thereafter, the rainfall events leading to landslides are plotted in a log(E) vs log(D) graph and the distribution fitted to the power law equation. The equation of the threshold is $E=\left(\alpha \pm \Delta \alpha \right).D^{\left(\gamma \pm \Delta \gamma \right)}$, where E is cumulated event rainfall (mm), D is duration (h), α is intercept, and γ is the slope of threshold curve [4]. The uncertainties Δα and Δγ are determined using a bootstrap nonparametric statistical technique [32]. This uncertainties measure the variation of the threshold around a central tendency line, which depends on multiple factors, but primarily on the number and the distribution of the empirical data points representing different rainfall conditions that have resulted in landslides [33]. The distribution of rainfall conditions which have resulted to triggering of landslides is fitted the power law equation in a log-log graph. The thresholds were determined using the methodology proposed by [4] and further modified by [32] for various exceedance probabilities ranging from 1% to 50%.

The equation involves the discrete and continuous maximum likelihood function estimation to fit the data in agreement with the equation. There are two assumptions involved with the use of power law equation. (1) With the increase in the cumulated rainfall, there is a nonlinear increase in the probability of landslides. It asserts that the possibility of landslide decreases when threshold reduces and vice versa. (2) With the increase in rainfall duration, the occurrences of slope failure reduce [27].

The study of antecedent rainfall for landslide incidences is important as it may lead to an increase in soil moisture content leading to slope instability. The impact of antecedent rainfall should be a site-specific study and may not always hold good for other regions with similar geological and rainfall conditions [34]. The variation in soil moisture content across an area is difficult to accurately determine as it depends on various factors like the variation in soil type, depth, climatic variation, etc. [7]. Various authors have used different periods to determine the correlation between antecedent rainfall and number of days for landslide triggering. [35,36,37,38] examined for 3, 4, 18, and 180 days respectively. [14] used 7, 10, and 15 days, whereas [39] assessed 2, 5, 15, and 25 days based on a trial and error basis. In this study, we considered 3, 7, 10, 15, 20, and 30 days.

5. Results and discussion

In this study, an event rainfall–duration threshold has been determined using available rainfall and landslide data [32,40]. The definition of rainfall and landslide event to determine any kind of threshold is very critical [18]. In the present study, to determine the thresholds, the landslides after the initial failure were not considered, i.e., if, on a particular day, five landslides occur, they will be recognized as one landslide event. This approach is similar to other works [40] and reduces the number of rainfall events with landslides. After that, landslide events which initiated due to very low rainfall values (lower than 25th percentile) were discarded to determine more accurate thresholds [25] as such landslide incidences may not be solely initiated by rainfall. Thus, the number of landslides reduced to 51 landslide events. With the method explained in the previous section, the threshold came out to be $E=\left(5.68\pm 1.80\right).D^{\left(0.70\pm 0.04\right)}$. Figure 5a depicts the 51 rainfall conditions, which led to landslides, threshold at 5% (T_5,B) exceedance probability and the associated uncertainty (Figure 5b,c). The relative uncertainties for the $\gamma$ parameter are less than the acceptable value of 10% [32]. However, the uncertainty for α is slightly higher at 31.7%, which can be attributed to limited number of empirical data used to determine thresholds [25]. The power-law function and maximum likelihood estimations are used to get the best fitted curve by utilization the present set of data. The lower boundary points are determined using regression analysis and a standard deviation is used to understand the distribution of normally distributed data from the mean value [27]. The graph is drawn on a logarithmic scale. [25] also defined ED thresholds for Chukha region using a semi-automatic empirical approach and defined threshold as E = (7.3 ± 2.0). D^{(0.71 ± 0.06)} using 43 landslide data points for 2004–2014. The thresholds defined for Kalimpong region, India by [11] depicted ED relationship as E = 3.52. D^0.41.

To evaluate the importance of antecedent rainfall the rainfall events during the monsoon of year 2012 were analyzed, a total of 50 rainfall events occurred which led to eight landslide events. As there were more than 1000 rainfall events during 2004–2014, the landslide events of 2012 were selected to understand the significance of antecedent rainfall. Figure 6 depicts the daily (R_day) and 30-day (R_30-day) antecedent precipitation values for 2012. The daily rainfall values were comparatively higher on the day of landslide occurrence. Adopting the calculated threshold for triggering landslide, the number of daily rainfall events exceeding the reference period is high which depicts that antecedent precipitation plays an important role in landslide initiation [41]. The landslide events during this period were characterized by an R_30-day value of 350 mm. Considering the R_30-day value of 350 mm, the number of rainfall events exceeding it is 14. This analysis resulted in fewer data points, which depict the significance of antecedent rainfall, and thus its effect on the distribution of soil moisture during the initiation of triggering rainfall [41].

To determine the number of antecedent rainfall days significant for landslide triggering, a trial and error approach considering various days was used [11,27]. The daily rainfall on the day of rainfall for landslide events was plotted against antecedent rainfall for various time durations of 3, 7, 10, 20, and 30 days (Figure 7). The diagonal line divides the graph to differentiate the scattering bias of daily rainfall (ordinate) and antecedent rainfall (abscissa). The diagonal depicts that the daily precipitation data on the day of failure and the antecedent precipitation prior to failure are same [27]. Figure 7 shows graphs of daily rainfall corresponding to antecedent rainfall of various periods. As observed in Figure 7, the majority of the landslide events are biased towards antecedent rainfall as compared to daily rainfall. In the case of 3-day antecedent rainfall, 19.6% of the landslides are biased towards daily rainfall and the remaining 80.4% (41 of the 51 landslides) are biased towards 3-day antecedent rainfall prior to failure. Similarly, for 7-day antecedent conditions, the biasness of daily rainfall towards landslide initiation decreases to 13.7%. In the remaining cases (10-day, 20-day, and 30-day) the effect decreases to 1.9%. Such a comparison between the daily and antecedent rainfall data would give more clarity with the availability of hourly data. The observation of plots for similar bias antecedent conditions shows that a greater number of points in the case of 10 days is concentrated along the abscissa in comparison to other cases where a scattering of points along the abscissa is prevalent [27]. Therefore, it can be concluded that an antecedent rainfall of 88.35 mm for a minimum of 10 days provides the best correlation for triggering of landslides in the region.

The use of thresholds for an operational early warning system can be justified by validating it with an independent dataset which is lacking in various studies conducted, as mentioned in the review article by [19]. The present study validates the determined thresholds using the rainfall data of 2015 by determining the threat score [3,42]. Threat score (TS) is defined as the number of true positive cases (TP) divided by the summation of true positive (TP), false negative (FN), and false positive (FP) cases [43].

$$TS=\raisebox{1ex}{$TP$}\!\left/ \!\raisebox{-1ex}{$TP+FN+FP$}\right.$$

(1)

During this period, eight landslide events occurred due to a total of 46 rainfall events. As determined earlier, a R_30-day precipitation value of 350 mm could be used as a pre-filter, and it was found that only 11 of the 46 rainfall events exceeded the value. Thereafter, the biasness of the landslide events with respect to daily and antecedent rainfall was carried out and the results show that only one case was slightly biased towards the daily rainfall whereas the rest were biased towards antecedent rainfall. This analysis shows that a R_30-day value could be used as a pre-filter for determining thresholds and antecedent rainfall plays a significant part for landslide initiation in the region. Finally, the threat score for 2015 landslide event was calculated using the ED threshold value of 53.3 mm. The results determined were TP equal to six, FN and FP were two and three respectively, and so TS was found to be 0.54. The result shows that the rainfall thresholds can be used the first step and eventually the threshold effectiveness will be improved in time when additional data will be collected, as shown in other long-term projects [12]. However, when using as an early warning system the effect of daily as well as antecedent rainfall needs to be considered.

6. Conclusions

Rainfall-induced landslides are one of the major destructive natural disasters in Bhutan. However, minimal study has been conducted to develop a landslide early warning system, either regionally or locally. This paper attempts to determine rainfall thresholds in terms of event rainfall–rainfall duration and antecedent rainfall for landslide incidences using daily rainfall data for Chukha region located in the south-western part of Bhutan Himalayas. The analysis was carried out from the available rainfall and landslide occurrences between 2004 to 2014. The majority of the landslide occurrences are along the Phuentsholing–Thimphu highway, which is an important road connection for the country as it connected the capital with India and used for trade purposes. The thresholds determined for the region depicts a minimum event rainfall of 55 mm for short duration of 24 h can trigger landslides. Thereafter, the significance of antecedent rainfall was carried out using the rainfall events of 2012 monsoon. The analysis revealed the importance of antecedent rainfall and therefore it was necessary to determine the antecedent time window necessary for landslide occurrence. The biasness of the rainfall events which resulted in landslides revealed that a 10-day antecedent rainfall would provide the best correlation for landslide occurrences in the study region. Further, it can be stated that R_30-day value of 350 mm could be used as a pre-filter before the use of ED thresholds. The derived rainfall thresholds can be improved with the availability of hourly rainfall data to use for an effective warning system.