Assessment of climate change impact on hydrology of a transboundary river of Bhutan and India

Water is considered as one of the most dominant natural resources and a large cluster of living organisms depend on it. This basic resource is now under threat due to increase in population and urbanization, resulting in impacts of climate change such as increase in global average temperature, rainfall variability and severe extreme events (Shrestha et al. 2021). According to IPCC, climate change effects and adaptation measures are of major global concern (IPCC 2014). Although water is one of the abundant resources in Bhutan, there is a major risk to it due to the climate induced changes and an extensive climate change impact study is urgently needed in the region. Based on the observation, the unpredictability of climate change is likely to affect the Indian subbasins as this region is highly reliant on glacier and snow melt to meet the freshwater demands. The glacier retreat in the Indian high mountains greatly influences the region's population and economy (Singh et al. 2017). A study in Vijila river, which is a transboundary river between Lithuania and Belarus, states that transboundary basins are at higher threat because of the shift in seasons and the two countries' policymakers need to act jointly (Čerkasova et al. 2018). A study in Bhramaputra river, which is a cross border river between India and Bangladesh, states that longer periods of flooding combined with transboundary rivers can affect millions of populations and this can be addressed with a combination of various remote sensing methods with modeling (Dubey et al. 2021). Therefore, to minimize the risks associated with climate change, the climate change must be understood, quantified and incorporated into regional water management schemes.

The potential impact of climate change in the Himalayas is said to be more obvious, as glaciers and melting snow dominate the river (Lutz et al. 2014). According to Oertli et al. (2005), the freshwater ecosystem is closely related to its natural hydrology at multiple geographical and temporal scales. Therefore, river hydrology remains a very important factor in preserving the river ecosystem's structure and function. The elements of flow regime including high, low, and medium flows are critical for biodiversity conservation. The drastic change in different climate parameters will hence affect the river hydrology which will have a tremendous impact on water management in the basin (Budhathoki et al. 2020). Change in timing and magnitude of the stream flow may also change the operating rules of reservoirs like flood control, generation of hydropower, water supply and recreational applications (Eum et al. 2017).

The rivers from Bhutan flow to India, but apart from mutual comprehension, there has been no agreement between the two nations on the management of these transboundary rivers. A relationship between Bhutan and India that helped both the countries to develop is said to be a good example of transboundary water management (Biswas 2011). The major source of revenue for Bhutan is the hydroelectric power which is exported to India and the country is gearing up towards the construction of many more hydro projects in future. Bhutan has agreed to produce electricity and supply to India and India has agreed to help Bhutan in building up those projects. Water accessibility has been a major problem in Bhutan, and it needs construction of storage for proper water management. If Bhutan wants to construct storage over those rivers in future, there may arise a conflict between the two nations (Khan 2017). With the changing climate affecting the river hydrology, the economic development and the relationship of Bhutan with India can both be affected. As reported by Goldenman (1990), one of the major challenges when a country must share a resource with another country is to develop a principle, procedure and institutions that can be agreed by both the countries, while at the same time adapting to the changing climate. Moreover, only a few studies concerning the transboundary river of Bhutan in relation to the changing climate have been conducted. With a limited number of climate projection studies on mountain terrain (Akhtar et al. 2008), especially for Bhutan, the outcome of this study will help policymakers and people from Bhutan and India understand the future climate scenario in the study area. In addition, transboundary rivers require a genuine utilization and regulation of water resources, and this will be more strenuous without considering climate change. In this regard, the Indicators of Hydrologic Alterations (IHA) method has been widely used to precisely distinguish the flow regimes with the help of ecological information and its interconnection with respect to hydrological and ecosystem elements (Richter et al. 1998; Brouziyne et al. 2021) including inter and intra annual variability of the streamflow conditions (Stefanidis et al. 2016).

Due to the harsh weather, rugged terrain, and few observed data, very few studies have been conducted in the Indian Himalayan region (Singh et al. 2017) which limits the knowledge of possible impacts of climate change and increases the uncertainties associated with it. In this regard, the purpose of the study is to assess future climate variables and analyze climate change trends in the Wangchu and Raidak river area. In addition, this study projects future runoff using the hydrological model Soil and Water Assessment Tool (SWAT) due to climate change in the study area, and evaluates possible alteration in the hydrological regime of Wangchu and Raidak river due to climate change. Overall, by understanding climate change and its impact on the hydrology of a river shared between Bhutan and India, policymakers, water managers and individuals from both nations will gain insight into the future choices and adaptation measures needed, and water security can be retained at individual, domestic and global levels.

The Wangchu River in Bhutan and Raidak River in India is a transboundary river that flows from Bhutan to India. After leaving Bhutan this river enters Alipurduar, a district in the western Bengal State of India where it is known as Raidak River (Fakhruddin 2015). The river originates at Mount Akungphu at the Himalayan Range at an altitude of 6,400 masl in Bhutan and in India it confluences with the river Brahmaputra at chainage of 370 km (Figure 1).

In Bhutan, the Wangchu river plays an important economic and social role as it supports two large hydroelectric plants at Tala (1,020 MW) and Chukha (336 MW). The average net monthly outflow to India is approximately 5,200 MCM, which is 7.4 percent of Bhutan's total river flow into India (NEC 2016). Downstream communities in India are mainly using the river for agricultural purposes. The total area covered by the watershed is 4,820 km² and the length of the river is 160 km. Table 1 shows the main physiographical characteristics of the study area.

The basin comprises 7 meteorological stations for rainfall and temperature measurement and 5 hydrological stations for discharge measurement. The geo-spatial data used for hydrological simulations using the SWAT model were topographic, soil and land use data. The topography of the study area was extracted from the digital elevation model (DEM) using the ArcGIS map. Six different types of soil were found in the basin with Dystric Nitosols covering the maximum area in the selected basin. As per the investigation conducted by FAO in the harmonized World Soil Database, 75% of the area is composed of sandy clay loam and remaining 24% by loam (Fakhruddin 2015). The land cover includes forest, wood land, open shrub land, wooded grass land and other land use type. As per the report by MoAF, in the upper part of the study area the pressure on forest is increasing with the rapid growing population and developmental activities and as per the study conducted by Roy & Basak, 2018, in the lower part of the study area, especially near the outlet, between 1999 and 2006, the agricultural and allied lands (fallow, river vacant and barren land) have increased by almost 27.39 km² and water bodies have decreased by 2.61%. In the lower region (India), the development of the tea industry and the recent growth in the tourism sector has altered the land use pattern with different associate infrastructural development.

The projection of climate and assessment of hydrologic alteration depends on meteorological data. The meteorological data used for the study included daily precipitation, maximum and minimum temperature. The data were collected from the National Center for Hydrology and Meteorology in Bhutan. The details of the meteorological and hydrological stations used in the study are illustrated in Table 2. Apart from ground-based data, global gridded data such as APHRODITE daily rainfall of resolution 0.25 × 0.25° and CPC daily maximum and minimum temperature data of resolution 0.50 × 0.50° were used for the climate analysis.

The data for the ground-based stations were available mostly from 1996. For the climate change impact studies, the baseline data required was from 1980 and the data at the Indian counterpart was not available. Due to limited data and inconsistency in data period at all the stations, the observed data could not be used for climate projection. Therefore, the validated temperature and precipitation global datasets were used for the study. The details of global gridded datasets used are given in Table 3.

The future climate prediction generated by climate models plays an essential part in determining the impacts due to changing climate. Based on the research conducted by Babel et al. (2014), the GCMs are not preferred for hydrological modelling in mountain areas. Hence, three RCMs, namely, ACCESS-CSIRO-CCAM (ACCESS), CNRM-CM5-CSIRO-CCAM (CNRM) and MPI-ESM-LR-CSIRO-CCAM (MPI) of 0.5° resolution, were used for the future climate dataset. Similar RCMs were used in other studies with similar topography in South Asia (Budhathoki et al. 2020; Shrestha et al. 2021). The CORDEX RCM data were downloaded from http://cccr.tropmet.res.in/home/index.jsp. There are generally four scenarios for which the climate projection is made, namely RCP 2.6 which is considered to be low emission scenario, RCP 4.5 and RCP 6 which are considered to be scenario, with medium emission and RCP 8.5 which is assumed to be scenario with high emission. In this study, the climate change over the study area was shown by change in daily rainfall, maximum and minimum temperature under two emission scenarios, RCP 4.5 and RCP 8.5 (Gebre & Ludwig 2015; Nikam et al. 2018; Budhathoki et al. 2020; Shrestha et al. 2021).

In this study, climate data was divided into historical and future periods. The historical period is named as baseline (BL) period from 1980 to 2005. The future study has been carried out for two future periods, namely, near future (2025–2050) (NF) and far future (2075–2099) (FF). The classification was done based on the availability of the future climate data.

Figure 2 represents the overall methodological framework adopted for the study. The methodology consists of climate change projection, hydrological model development and climate change impact on hydrology.

The simulation from GCMs is usually found to deviate from the observed climatological data due to systematic and random model errors (Shrestha et al. 2017). Since RCMs are downscaled from GCMs, RCMs usually inherit a part of error from the GCMs. In order to match the original data and minimize the difference in means for each weather variable the extracted RCM data need to be bias corrected. Shrestha et al. (2017), conducted a study on bias correction techniques and found that simple methods like linear scaling can also be sufficient for hydrological analysis. Accordingly, future RCM climate data were bias corrected using two approaches, i.e. linear scaling and quantile mapping methods. Parameters such as correlation coefficient (R²), standard deviation (SD) and root mean square error (RMSE) were used for the evaluation of suitability of method. Based on the values of those parameters, linear scaling was adopted for temperature data while the quantile mapping method was used for rainfall bias correction, similar to Budhathoki et al. (2020).

Overall, comparing the values of standard deviation (SD) and root mean square error (RMSE) amongst two methods, the quantile mapping method was found to be suitable for rainfall and for temperature, the linear scaling method was found to be better. Therefore, this study opted for two different methods for the bias correction.

Due to limited observed data for the study area, the validated gridded global sets were used in place of observed data. A few models were built and run without any calibration and validation with different global datasets. The model results were than equated with the available observed data and displayed a good relationship with the available observed flow; they were then used for further analysis. Apart from model simulation with different datasets, the performance of gridded global datasets was also evaluated using statistical indicators and RClimdex tool.

Ten temperature and rainfall extreme indices were also analyzed in the study. The extreme indices for temperature are Frost days (FD0); Number of summer days greater than 25 °C (SU25); Coldest daily minimum temperature (TNn); Warmest daily minimum temperature (TNx); Coldest daily maximum temperature (TXn); and Warmest daily maximum temperature (TXx). The extreme indices for precipitation are: Consecutive dry days (CDD); Consecutive wet days (CWD); Wettest Consecutive five days (R5); and Precipitation from very wet days (R95). They were analyzed for the selected two future periods under RCP 4.5 and RCP 8.5.

The Soil and Water Assessment Tool (SWAT) was used to simulate the hydrology under climate change scenarios. SWAT is a physically based semi-distributed hydrological model which is widely used to quantify the impact of change in land use, climate and vegetation on the streamflow and water quality in watersheds. Several studies have reported the good performance of the SWAT model in the watersheds of the Himalayan regions (Rostamian et al. 2008; Bharati et al. 2014; Palazzoli et al. 2015; Rajib et al. 2016; Budhathoki et al. 2020).

The SWAT model has been used to study the impact of climate change on the hydrology of transboundary rivers. Kaini et al. (2020) applied the SWAT model to analyze climate change impacts on the hydrological regime of a river basin and its implications for future irrigation water availability in the Koshi River basin. They found that the average flow in the Koshi River is projected to increase in future. Čerkasova et al. (2018) applied the SWAT model to assess the impact of climate change hydrology and water quality in a transboundary river, Nemunas River in Europe. Hajihosseini et al. (2020) also used the SWAT model to analyze the impacts of land use changes and climate variability on the transboundary Hirmand River of Afghanistan. Similarly, the model has been extensively applied in the Lower Mekong Region to investigate the impact of climate change, land use change and demographic change on the hydrology and water quality (Oeurng et al. 2016; Trang et al. 2017).

The model has been successfully applied in the ungauged watersheds where there are few hydrological monitoring stations. The data needed for modelling SWAT are usually available from the local governments. The model makes use of DEM for delineation of watersheds and divides the delineated watershed further into sub basins and HRUs (hydrologic response unit) based on land and soil information. The water balance equation used by the model and steps involved during simulation are based on equations of Neitsch et al. (2005) for the SWAT model. To quantify the reliability of the models' output, model evaluation is necessary. The output of the model is considered to be reliable if the performance of the model falls within a certain limit. The criteria set by Moriasi et al. (2015) were used for evaluation of the model performance and accordingly NSE, R² and PBIAS.

The impact of climate change on streamflow was analyzed using the range of variability approach (RVA). The RVA is the most widely used multivariable approach for the assessment of flow regime alteration (Richter et al. 1998). In an RVA study, the full range of pre-impact data is divided into three different categories for each parameter:

• High Range: Values in a set of annual values above 67th percentile

• Middle Range: Values between 33rd and 67th percentile in a set of annual flow values

• Low Range: Values below 33rd percentile in a set of annual flow values

HAF allows the assessment of the range of flow regimes significantly affected by climate change. A positive hydrologic alteration value indicates an increase in frequency from pre-impact to post-impact period while negative indicates decrease in frequency. The degree of alteration for each index (D_i) is equal to the absolute value of HAF. Richter et al. (1998) further suggested that if the value of Di is between 0 and 33 percent, then it represents that there is little or no alteration (i.e. low alteration), 33–67 percent represents moderate alteration and 67–100 percent represents high alteration. The degree of hydrological change ‘D’ is a measure that quantifies the deviation of the post-climate change flow regime from the pre-climate change flow.

Figure 3(a) and 3(b) shows the NF and far FF average temperature change compared to the baseline (BL) (1980–2005) period. Under both scenarios the temperature has an increasing trend. The average temperature under RCP 4.5 and RCP 8.5 is projected to increase by 0.02 °C and 0.05 °C per year in the future, respectively. Overall, the average annual temperature in the study is expected to rise by 1.5 °C by the end of the century under RCP 4.5 and 3.6° under RCP 8.5. The increase in minimum and maximum temperature is consistent with the study conducted by Immerzeel (2007) in the Brahmaputra Basin from 2000 to 2100 based on the model of six downscaled GCMs. Likewise, the average annual rainfall under RCP 4.5 is projected to decrease by 1.9% and under RCP 8.5, it is projected to increase by 2.86% in the near future (2025–2050) and decrease by 1.3% in the far future (2074–2099).

It can be observed that the number of days with temperature less than 0 °C has decreased compared to baseline and the annual average maximum of daily maximum and daily minimum temperature are projected to increase compared to baseline with the highest increase in the far future under both the scenarios in Table 4. It can also be observed that the duration of extreme climate events is decreasing whereas the intensity of the extremes is increasing. The maximum number of consecutive days with rainfall less than 1 mm in a year was found to have reduced compared to baseline under the scenario for the near future whereas a slight increase was shown for the far future. The number of days with rainfall greater than 1 mm in a year have reduced compared to baseline irrespective of scenario and time period. However, it can be observed that the intensity of rainfall is high compared to baseline under both scenarios and time periods.

The SWAT model was calibrated and validated at Chimakoti station, and the result of its performance is presented in Table 5. The model was calibrated using discharge data of only one station due to the lack of data as the stations are in the other side of the border. The topographic data, land use and soil data along with climate and hydrological data were fed into the model with a warmup period of 5 years. Further elaboration of the data resolution has been stated in section 2.2. The calibration period was from 1990 to 2000 and validation period from 2001 to 2005. The calibration process of the model is the crucial part in the hydrological modeling. Identification of key parameters based on the study objectives is an essential step for model calibration (Ma et al. 2000). The model was calibrated using the automatic calibration software SWAT-CUP under SUFI-2 optimization algorithm (Abbaspour 2008; Arnold et al. 2012). To perform the uncertainty analysis, calibration, and validation of the hydrological simulation outputs, SWAT-CUP was used as an interface between the SWAT and calibration algorithms. The SWAT-CUP is a useful tool for automatic model calibration and sensitivity analysis. This program consists of five different procedures, e.g. GLUE, ParaSol, SUFI-2, IS, and MCMC for model calibration and sensitivity analysis (Khoi & Thom 2015). The SUFI-2 is the most popular and efficient procedure compared to others (Yang et al. 2008).

The determination of which modeling parameter values to calibrate is followed by parameter estimation. The most important factors in effective automated model calibration are parameter recognition, data quantity and accuracy, appropriate choice of a performance measure, and implementation of a proper search method.

Sensitivity analysis was performed for 36 parameters prior to calibration. The sensitivity analysis is a very important process to reduce the number of parameters prior to calibration. A total of 19 out of 36 parameters were found to be sensitive. The sensitive parameters and their rank, along with the range and fitted value are listed in Table 6.

The sensitivity of parameters can reflect the hydrological processes and their governing factors in the watershed. The top 10 ranked parameters are related to channel characteristics, groundwater properties and soil characteristics of the watershed. The top ranked parameter in the sensitivity analysis is CH_K2 (channel effective hydraulic conductivity) which regulates the transmission losses caused by surface runoff as it flows through the main channel and the fitted value lies within the range. This shows that the surface hydrology in this watershed is primarily governed by channel property. Followed by the parameters are: OV_N and CH_N2. The OV_N is an effective roughness coefficient that accounts for the effects of raindrop impact, channelization of flow into rills, obstacles such a ridges and rocks, friction over the land surface, and erosion and transport of sediment. The lower value of ALPHA_BF shows that the basin is slow in response to recharge.

One of the significant parameters is curve number, which is primarily influenced by the types of land cover in the basin. The response of the main land cover groups is the source of the high sensitivity to CN2. The dominating soil categories in the watershed belong to the soil hydrological groups A and B, which are characterized by high to medium infiltration rates, high to medium water transport rates, and high to incredibly high drainage capability. This explains the sensitivity to SOL_AWC, ESCO, and GWQMN. In summary, it can be generalized that the basin's response to surface runoff processes is faster than that of groundwater recharge and hence the basin stores a lower amount of groundwater.

Futher, the parameters focusing on discharge quantity were adopted from Muleta & Nicklow (2005); Neto et al. (2014); Rajib et al. (2016); and Rostamian et al. (2008).

According to the model performance evaluation classification suggested by by Moriasi et al. (2015), there was a good agreement between the observed and simulated flow at Chimakoti station. The plots between observed and simulated flow along with the calibration and validation results are shown in Figure 4(a) and 4(b). The base flow is well captured by the model, although there is minimal variability in peak flows.

The calibrated and validated model was then used for the flow simulation. The statistics on flow are shown in Table 7. The data shows that the average flow is slightly lower than the baseline flow except for RCP 8.5 in near future, where the average annual flow is 2.70% higher compared to the baseline flow. From the water availability point of view, the climate change will not have a significant impact on the water resources of the study area. The flows under the scenarios are projected to be higher in the near future compared to flows in the far future. The high variability in maximum and minimum flow, low value of Q90, i.e. the probability of exceeding the given flow by 90% of the time and high value of flood flows (Q5), suggests building storage facilities to meet downstream water requirements for irrigation, domestic and industrial uses, hydro-power production and controlling floods.

Storage systems such as dams or reservoirs are one of the many strategies for adapting to climate change. In the lack of such storage, the removal of water from the river or the operation of hydroelectric power plants relies on the monthly availability and flow variability. The variation in the projected average monthly flow with respect to the baseline flow is provided in Table 8. The change in monthly flow varies from +114.63% to –38.31%. Most of the flows are seen to decrease during the dry season and increase during the wet season. The variation shows a need for a storage system in the study region and changes in the operating regulations to address changes in monthly flows.

The change in flow peak with change in flow magnitude under different scenarios is shown in Figure 5. During the baseline period, the peak flow of 189 m³/sec is in September, whereas in the near future the peak flow of 185 m³/sec under RCP 4.5 and 200 m³/sec under RCP 8.5 and the peak flow of 159 m³/sec under RCP 4.5 and 150 m³/sec under RCP 8.5 in the far future is in August. The change in the peak flow is due to climate change and indicates the melting of snow.

The high and low flow at the study area outlet is demonstrated in Figure 6. Compared to baseline, under RCP 4.5, the high flow is expected to decline by 4.89% and 17.40% in the near and far future respectively, and under RCP 8.5, the high flow is projected to increase by 3.78% in the near future and decline by 16.72% in the far future. Low Q90 value and fluctuations in flow values justifies the need to construct storage within the study region.

For three RVAs ranges of monthly SWAT flows corresponding to two RCPs and two future periods the degree of hydrological alteration was calculated, and the result is shown in the Table 9.

It was found that the degree of alteration in the near future is always lower than that in the far future except for mid-range flow under RCP 4.5, where the degree of alteration is 2% higher in the near future compared to that in the far future. The degree of alteration for low RVA flow in the long term under both scenarios is relatively large, suggesting that the low flow range will mostly be changed with climate change.

Water balance components play a prominent part in the hydrological cycle and contribute to the river discharge. Comparison of the water balance components of the future period with the water balance components of the baseline period has been studied. Figure 7 shows the impact of climate change on the annual average water balance components at the outlet of the study area during baseline, near and far periods.

The model anticipated a decline in snowmelt relative to the baseline, which is consistent with research by Prasch et al. (2011), and that the ice melt will reduce as the size of glacier snow decreases. The rate of potential evapotranspiration is seen to decline in the near future, but then increase in the far future under both scenarios. Evaporation is seen to increase relative to baseline, but the rise is very minimal. The surface runoff appears to be more in the near future, but less in the far future. Furthermore, the water yield in the study area is found to decrease by 14 and 48 mm annually under RCP 4.5 in the near and far future respectively compared to baseline, and under RCP 8.5, it is expected to increase by 16 mm in the near future and decrease by 65 mm in the far future compared to baseline.

The key objective of the study was to assess the impact of changing climate on the hydrology of the transboundary river that flows from Bhutan to India. The river is known as Wangchu in Bhutan and Raidak in India. The future bias corrected climatic data (rainfall and temperature) under RCP 4.5 and RCP 8.5 scenarios for three different time periods i.e. baseline 1980–2005, near future from 2025–2050 and far future from 2074–2099 were forced into a SWAT model. The average annual temperature in the study area was projected to increase by 1.5 °C and 3.6 °C respectively under RCP 4.5 and RCP 8.5. Likewise, the average annual rainfall was projected to decrease by 1.9% under RCP 4.5, increase by 2.86% in the near future and decrease by 1.38% in the far future under RCP 8.5. The flow statistics for historical and expected flow indicated that, with one exception in the near future under RCP 8.5, the future flow will be lower than the baseline flows. The flow under both scenarios was projected to be greater in the near future (2025–2050) compared to the far future (2074–2099). Monthly flows were projected to vary from (+) 114.63% to (−) 38.31%. Building storage facilities to meet downstream water requirements for farming, domestic and industrial applications, generating the required hydropower and controlling floods looks like a viable solution to varying flows. The peak flow was also projected to shift from September to August in future, which indicates the impact of climate change on snow melt.

In addition to statistical flow analysis, the river's alteration in hydrological regime was also performed using the IHA tool by a range of variability methods. It was found that the low flows in the future will be highly altered compared to high and middle flows. Overall, it was found that the flow of the river would decrease in line with the projected decrease in rainfall. The outcome of this study would be useful in understanding the potential impact of climate change on Wangchu River in Bhutan and Raidak River in India, and in helping governments and individuals plan adaptation strategies accordingly.

Most of the rivers in Bhutan flow to India but apart from mutual understanding between the two countries, no treaty or agreement has been drawn that concerns the management of this transboundary river. Detailed research regarding the impact of climate change on both upstream and downstream can be conducted. Therefore, a study with mutual sharing of data between India and Bhutan will also help in research, planning and development of both the nations. Furthermore, in order to carry out any research related to climate, a large amount of observed meteorological and hydrological data is very important. The observed data in Bhutan for all the stations are really limited. Since the observed data are limited, a detailed study of global data can be used for all the rivers in Bhutan. Decision makers should advocate the need to invest in mutual collaboration and urge legislators to recognize that each territory presents its own series of diverse issues that enable transborder agencies to be flexible and take particular care in crisis management.

The authors would like to acknowledge the National Center for Hydrology and Meteorology, Bhutan for providing the required data, and the Asian Institute of Technology (AIT) and Thai Pipe Industry Co. Ltd for providing the necessary support to conduct the study.

Data cannot be made publicly available; readers should contact the corresponding author for details.