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

Assessing the impacts of climate change on a transboundary river plays an important role in sustaining water security within as well as beyond the national boundaries. At times, the unilateral decision taken by one country can increase the risk of negative effect on the riparian countries and if the impact is felt strongly by the other country, it can lead to international tension between them. This study examines the impact of climate change on hydrology between a shared river which is Wangchu river in Bhutan and Raidak river in India. The river is mainly used to produce hydropower in the two largest hydropower plants on which the majority of Bhutan’s economic development depends and is mainly used for agriculture in India. The Soil and Water Assessment Tool (SWAT) was used for future flow simulation. Future climate was projected for near future (NF) from 2025–2050 and far future (FF) from 2074–2099 using an ensemble of three regional climate models (ACCESS, CNRM-CM5 and MPI-ESM-LR) for two RCPs (Representative Concentration Pathways), RCP 4.5 and RCP 8.5 scenario. The ensemble results indicated that, in future, the study area would become warmer with temperature increase of 1.5 °C under RCP 4.5 and 3.6 °C under RCP 8.5. However, as per RCP 4.5 and RCP 8.5, rainfall over the study area is projected to decrease by 1.90% and 1.38%, respectively. As a consequence of the projected decrease in rainfall, the flow in the river is projected to decrease by 5.77% under RCP 4.5 and 4.73% under RCP 8.5. Overall, the results indicated that the degree of hydrological change is expected to be higher, particularly for low flows in both Wangchu and Raidak River. Since transboundary water is shared for economic growth, climate change adaptation and opportunities should also be considered by both the nations for better water management.


INTRODUCTION
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 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 2 and the length of the river is 160 km. Table 1 shows the main physiographical characteristics of the study area.

Historical data
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 2 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.

Future climate projection
The simulation from GCMs is usually found to deviate from the observed climatological data due to systematic and random model errors . 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 2 ), 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.

Hydrological modeling
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 2 and PBIAS.

Flow alteration analysis using range of variability approach
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: The software then determines the predicted frequency and the frequency at which post-impact values of IHA (Stefanidis et al. 2016;Brouziyne et al. 2021) parameters will fall respectively within each of the three categories. This expected frequency is equal to the number of values in the category during the pre-impact period multiplied by the ratio of post-impact years to pre-impact years. Finally, for each of the categories, the hydrological alteration factor (HAF) will be calculated as shown in the equation below: HAF ¼ (observed frequency-expected frequency) expected frequency 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.

RESULTS AND DISCUSSION
3.1. Future climate scenarios Corrected Proof 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 .
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.

Hydrological modeling performance
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 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).

Corrected Proof
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.

Corrected Proof
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.
3.3. Impact of climate change on hydrology

Changes in future streamflow
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.

Changes in future peak flow
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 3 /sec is in September, whereas in the near future the peak flow of 185 m 3 /sec under RCP 4.5 and 200 m 3 /sec under RCP 8.5 and the peak flow of 159 m 3 /sec under RCP 4.5 and 150 m 3 /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.

Changes in future high and low flows
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.

Changes in future water balance
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. 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.