Hydrological response to climate change of the Brahmaputra basin using CMIP5 general circulation model ensemble

The Ganges–Brahmaputra–Meghna (GBM) river system plays an important role in China, Bhutan, India, Nepal and Bangladesh. The GBM basin is the third-largest freshwater outlet to the world's oceans (Chowdhury & Ward 2004). The Brahmaputra River contributes 67% of the total annual water flow of Bangladesh (Immerzeel 2008). This river basin is the main source of water in Bangladesh. The assessment of streamflows through this river can play a vital role for the water management of the country. However, estimation of water scarcity or water availability depends on the understanding of the hydrological system that is the main governing backbone of all kinds of water movement and water pollution (Jha 2011).

Global warming due to increasing concentration of greenhouse gases is likely to have a significant impact on precipitation, runoff processes and water resources (Arnell & Reynard 1996; Zhang et al. 2010; Haddeland et al. 2012; Cuo et al. 2015; Pervez & Henebry 2015). This raises the question of whether climate change is a threat to human water security or not. Previous studies have shown that simulated climate change impacts vary substantially depending on the climate model and emission scenarios used (Arnell 1999; Bronstert et al. 2002; Bronstert 2004; Wentz et al. 2007; Hurkmans et al. 2008; Wang et al. 2008; Montenegro & Ragab 2010). However, most of them focused on shifts in the timing of hydrological regimes where runoff is dominated by snow melt (e.g. Adam et al. 2007; Hurkmans et al. 2010; Arnell & Gosling 2013) or long-term mean annual streamflow and water availability (e.g. Christensen et al. 2004; Kling et al. 2014). Few studies have reported that if climate change can alter the risk of hydrological extremes at regional scales (Lehner et al. 2006), then developed basins with a dense population can become highly vulnerable to hydrological extremes.

Hydrological models can be classified based on model input or the types of principle it follows. A hydrological model can be classified as a lumped or distributed model based on the model parameters as a function of space and time. A lumped model considers the entire basin as a single entity whereas a distributed model can provide outputs that are distributed in space (Moradkhani & Sorooshian 2008). On the other hand, models can be classified as stochastic or deterministic. For a given input, a deterministic model always produces the same output whereas different outputs can be produced by a stochastic model. A model can also be classified as an event based or continuous model (Moradkhani & Sorooshian 2008). Event based models produce output for a specific time while the continuous model produces a time series of output. Most importantly, hydrological models can be classified as empirical, conceptual or physically based. Empirical models are based on observation with consideration of the hydrological data and thus these models are also known as data driven models. Models can be based on artificial neural networks (Islam 2009; Wu et al. 2009; Gholami et al. 2015), machine learning techniques (Taormina & Chau 2015), auto-regressive moving average (Wang et al. 2015), or a hybrid of two statistical models (Chau & Wu 2010; Chen et al. 2015) or machine learning. Conceptual model consists of a number of interconnected reservoirs to represent the physical process of the catchment. These models require a large number of datasets for calibration and validation. The Stanford Watershed Model is the first major conceptual model.

A number of physically based hydrological models are available and among them the most widely used models are the MIKE SHE (Système Hydrologique Européen) and SWAT (Soil and Water Assessment Tool) models. The MIKE SHE model has been developed by the Danish Hydraulic Institute which can simulate surface and groundwater movement and their interactions in a large basin. However, this model requires extensive physical parameters and is only commercially available. On the other hand, SWAT, a modelling tool developed by the US Department of Agriculture (USDA), is a semi-distributed complex physical model which has been developed to test and forecast the flow and nutrient circulation in ungauged basins. SWAT has been extensively used in many countries worldwide for flow prediction or to study the hydrological response to climate change (Spruill et al. 2000; Santhi et al. 2006; Zhang et al. 2010; Patel & Srivastava 2013; Fiseha et al. 2014). It has been extensively applied for hydrological water budgeting of the Ganges and 12 other major basins in India (Gosain et al. 2006). SWAT has several advantages including its computational efficiency, auto-calibration and uncertainty analysis techniques and is readily available as an open source tool. Moreover, the SWAT model is an effective tool in managing water resources (Tang et al. 2012) and achieving a high-quality calibration and uncertainty results (Abbaspour et al. 2004; Yang et al. 2008). SWAT has proven to be very successful in the watershed assessment of hydrology and water quality (Neitsch et al. 2002). It is a basin-scale, continuous-time model that operates on a daily time step and is designed to predict the impact of management on water, sediment, and agricultural chemical yields in ungauged watersheds. Major model components include weather, hydrology, soil temperature, plant growth, nutrients, pesticides, and land management. In SWAT, a watershed is divided into multiple sub-watersheds, which are then further subdivided into hydrologic response units (HRUs) that consist of homogeneous land use, management, and soil characteristics. The HRUs represent percentages of the sub-watershed area and are not identified spatially within a SWAT simulation. Alternatively, a watershed can be subdivided into only sub-watersheds that are characterized by dominant land use, soil type, and management. The water balance of each HRU in the watershed is represented by four storage volumes: snow, soil profile (0 to 2 metres), shallow aquifer (typically 2 to 20 metres), and deep aquifer (more than 20 metres). Flow, sediment, nutrient, and pesticide loadings from each HRU in a sub-watershed are summed, and the resulting loads are routed through channels, ponds, and/or reservoirs to the watershed outlet (Neitsch et al. 2002). The SWAT model was found to be very useful to study changes of flows of the semi-ungauged Ganges due to climate change (Narsimlu et al. 2013). Hence, it is expected that SWAT will be able to generate flows for the Brahmaputra River basin which is unique in hydro-morphological nature.

The Brahmaputra is a major transboundary river which drains an area of around 530,000 km² and crosses four different countries: China (50.5% of total catchment area), India (33.6%), Bangladesh (8.1%) and Bhutan (7.8%) (Gain et al. 2013). Figure 1 shows a map of the study area. Considering the climatology and elevation of the Brahmaputra basin, it can be classified into three physiographic units which will respond differently to the anticipated climate change (Immerzeel 2008). Elevations 3,500 m or above in the north basin, especially including China, are called the Tibetan plateau. The Himalaya belt, covers 28.6% of the basin, with elevations ranging from 100 to 3,500 m above mean sea level. The agricultural floodplains, of India and Bangladesh, cover 27.0% of the basin, with elevations up to 100 m above mean sea level, of which 40% are flood prone.

In the lower Brahmaputra, average temperature in winter is about 17°C and summer temperatures are on average as high as 27°C (Gain et al. 2011). Total annual precipitation is about 2,354 mm concentrated in the monsoon months of June, July, August and September (JJAS). The heaviest rainfall locations are Cherrapunji and Mawsynram, which are located just to the south of the Brahmaputra basin. Most of the runoff of this river is derived from heavy rainfall of 510–640 cm in the Abor and Mishmi hills in Arunachal Pradesh and 250–510 cm in the Brahmaputra plains. The annual precipitation over the Brahmaputra basin showed a decreasing trend in the Assam part and an increasing trend in the Bangladesh part (Mirza et al. 1998). A number of recent studies (Dhar & Nandargi 2002; Ghosh & Dutta 2012) suggested that the general circulation over the basin area undergoes abrupt seasonal changes during late spring and early summer, due to tropospheric warming over the Asian landmass, and causes summer precipitation. The mean annual value of such pre-monsoonal heavy rains shows a rainfall above 100 mm/day for 7.7 days and above 300 mm/day for 1.6 days for the observation period of 1993–2001. Overall, 66–85% of the annual rainfall occurs during the monsoon and 20–30% occurs during the pre-monsoon season, while a very small percentage of the annual rainfall occurs in winter. Rise of water in the Brahmaputra, with the increase in upstream rainfall, starts as early as April (Chowdhury & Ward 2004). The high stages last from June to August, and in some cases even until September. Several peaks occur during the monsoon with the first flooding event occurring in July and August and a second one of lesser volume occurring in September.

The focus of this study is to assess the changes of river flow in the lower Brahmaputra River basin where the hydrological impact of climate change is expected to be particularly strong due to glacier melt, extreme monsoon rainfall and sea level rise (Immerzeel et al. 2010; Gain et al. 2011). Average discharge of the Brahmaputra is approximately 20,000 m³/s. The climate of the basin is monsoon driven with a distinct wet season from June to September, which accounts for 60–70% of the annual rainfall (Immerzeel 2008).

Several types of data are required as input for SWAT to develop the model using the ArcSWAT interface. Topographic data have been obtained from the Shuttle Radar Topographic Mission (SRTM) with a spatial resolution of 90 m. Figure 1 shows the GBM basin. The sub-basin parameters, such as slope gradient, slope length of the terrain and the stream network characteristics (channel slope, length and width) have been derived from the analysis of the digital elevation model (DEM). The DEM has been masked for the Brahmaputra basin area as shown in Figure 1. A soil map of the study area was collected and extracted from the Food and Agriculture Organization (FAO) Digital Soil Map of the World (Figure 2(a)). A ‘usersoil’ database table was created for the study area from the available interpretations and lookup tables. Landuse maps are required for the delineation of the HRU of the model. A 300 m resolution landuse map called GlobCover based on data from 2009 to 2010 has been collected from the European Space Agency (Figure 2(b)). This data has been reclassified to match the SWAT land classes. One key assumption about the land use/land cover data is that in the future it would be the same as the present period. This assumption was considered as future land use/land cover projections were unavailable. The projection of the future land use and land cover pattern over the study area needs detailed modelling work which is beyond the scope of this study.

Watersheds have been delineated through automatic watershed delineation techniques. A threshold for minimum sub-basin area has been selected for the stream network and outlet calculation. As Bahadurabad station on the downstream of the Brahmaputra basin has been selected as the watershed outlet, this delineation resulted in a watershed area of 510452.91 km² and a total of 39 sub-basins.

In developing countries, there is a lack of full and realistic long-period climatic data. Therefore, the weather generator solves this problem by generating data from satellite data. The model requires the daily values of weather variables consisting of precipitation, maximum temperature and minimum temperature, solar radiation, wind speed and relative humidity over the study area. Due to unavailability of the ground-based observed data, satellite-based data products are the only available source of the data. The daily observed precipitation data can be obtained from three different sources: APRHODITE (Asian Precipitation – Highly-Resolved Observational Data Integration Towards Evaluation) (Yatagai et al. 2012), GPCP (Global Precipitation Climatology Project) (Adler et al. 2003), TRMM (Tropical Rainfall Measuring Mission) (Kummerow et al. 1998). Resolution of these precipitation data products are 25 km and they are available from 2000 onward. Paul et al. (2015) showed that out of the three data products, TRMM has shown better performance while comparing simulated flow with the observed flow. Temperature data were extracted from the ERA-Interim reanalysis data product. Due to the unavailability of any other high resolution (25 km) gridded data product over this region, ERA-Interim data has been chosen to derive the model. ERA-Interim is a global atmospheric reanalysis data product which has been continuously updated in real time since 1979. Basically, it uses a four-dimensional variational analysis (4D-Var) to produce gridded data products including a large variety of parameters related to describing weather, ocean and land surface conditions.

A weather generator model called WXGEN is integrated with SWAT. WXGEN generates climatic data or fills missing data using monthly statistics which are calculated from existing daily data. Meteorological data are analysed to determine the various statistical parameters, like, mean monthly maximum and minimum temperature, average monthly wind speed, number of rainy days, average daily solar radiation, standard deviation for air temperature, precipitation, skewness for daily precipitation, and probabilities of wet day following a dry day or a wet day, required as input by the weather generator file in SWAT.

The landuse and soil maps, as derived from the above-mentioned methods, have been imported and linked with the respective database table to create appropriate lookup tables. The ‘Multiple HRU’ option has been selected which generates 304 HRUs for the whole watershed. For TRMM datasets, a period of 5 years (2000–2004) has been selected for calibration and 5 years (2005–2009) for validation. In addition, 1 year has been kept as a warm-up period for both calibration and validation. A warm-up period allows the model to get a fully operational hydrological cycle and thus helps to stabilize the model. The main methods used in modelling the hydrologic processes in SWAT were the curve number method for runoff estimating, Penman-Monteith method for Potential Evapotranspiration (PET) and Muskingum method for channel routing.

It is important to downscale global circulation model (GCM) output from global scale (>150 km) to a regional scale (∼50 km) using regional climate models (RCMs). Our study area is located in the CORDEX South Asia domain in which RCMs are simulated. Outputs from only 11 RCMs are available over this region and all of them are considered for this study. Multiple regional models were chosen to capture the uncertainties of future flows due to the uncertainties inherent in the available climate models. RCMs have been simulated for the historical period from 1980 to 2009 which is considered as the baseline for historical simulations. Future river discharge for the Brahmaputra River has been estimated for three time windows, i.e., 2020s (2010–2039), 2050s (2040–2069) and 2080s (2070–2099). Table 1 shows the South Asia CORDEX climate scenarios. The calibrated model has been simulated to predict future discharge of the river under different climate models. Out of the four emission pathways adopted by the Fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC), the upper extreme scenario of RCP8.5 indicates the extreme climate change scenarios in the future. To understand the most extreme hydrological changes in the study area that may happen due to climate change, RCP8.5 scenario has been chosen. This will provide the opportunity to study the global warming level beyond 2°C which was motivated by the Paris Agreement of the United Nations Framework Convention on Climate Change (UNFCCC). The RCP2.6, RCP4.5, RCP6.0 scenarios have projected less global warming than the RCP8.5 scenario. Therefore, adopting the high emission RCP8.5 is advantageous for the study of the extreme climate change impact over the study area.

Many hydrological models contain parameters that cannot be determined directly from field measurements. Therefore, model calibration is used to adjust such parameters to optimize the agreement between observed and simulated values (Tolson & Shoemaker 2007; Zhang et al. 2009; Wu et al. 2012). There are numerous parameters in hydrological models which can be classified as physical parameters (i.e., parameters that are physically measurable from the properties of watershed) and process parameters (i.e., parameters representing properties which are not directly measurable) (Sorooshian & Gupta 1995). The options for calibration and uncertainty analysis techniques included in the SWAT are the MCMC (Markov Chain Monte Carlo) method (Vrugt et al. 2008), GLUE (Generalized Likelihood Uncertainty Estimation) (Beven & Binley 1992), ParaSol (Parameter Solution) (Yang et al. 2008), and SUFI-2 (Sequential Uncertainty Fitting-II) (Abbaspour et al. 2004). However, the SUFI-2 algorithm has been extensively applied (Abbaspour et al. 2004; Wagener & Gupta 2005; Yang et al. 2008) as it needs a minimum number of model simulations to attain a high-quality calibration and uncertainty results. Therefore, the SUFI-2 algorithm of SWAT-CUP has been applied in this study for model calibration. Under this algorithm, sensitivity analysis of parameters has been performed by regressing Latin hypercube generated parameters against objective function values. SWAT parameters were selected for model performance (calibration and validation) and included the following: ALPHA_BF, ALPHA_BNL, CH_N2, CN2, EPCO, ESCO, GWQMN, GW_DELAY, GW_REVAP, HRU_SLP, OV_N, REVAPMN, SFTMP, SLUSUBBSN, SOL_AWC, SOL_BD, SOL_K.

The ability of the SWAT model to simulate streamflow has been evaluated based on standard regression analysis which determines the strength of the linear relationship between simulated and measured data. For this, dimensionless indices were used which provide a relative model evaluation assessment, and error indices having the same units as those of the data of interest, which quantifies the deviation. Four widely used indicators were used to measure performance of the model at the calibration and validation stage as mentioned by Moriasi et al. (2007). These are coefficient of determination (R²), Nash-Sutcliffe efficiency (NSE), percent bias (PBIAS) and ratio of Root Mean Square Error (RMSE) and Standard Deviation termed as RSR.

Pearson's correlation coefficient (r) and coefficient of determination (R²) describe the degree of collinearity between simulated and measured data. The correlation coefficient, which ranges from −1 to 1, is an index of the degree of linear relationship between observed and simulated data. If r = 0, no linear relationship exists. If r = 1 or −1, a perfect positive or negative linear relationship exists. Similarly, R² describes the proportion of the variance in measured data explained by the model. R² ranges from 0 to 1, with higher values indicating less error variance, and typically values greater than 0.5 are considered acceptable (Van Liew et al. 2003; Santhi et al. 2006; Moriasi et al. 2007). Although r and R² have been widely used for model evaluation, these statistics are oversensitive to high extreme values (outliers) and insensitive to additive and proportional differences between model predictions and measured data (Legates & McCabe 1999).

The deviation term (DV) is used to evaluate the accumulation of differences in streamflow volume between simulated and measured data for a particular period of analysis. PBIAS has been selected for recommendation for several reasons: (1) DV has been recommended by ASCE (1993), (2) DV is commonly used to quantify water balance errors and its use can easily be extended to load errors, and (3) PBIAS has the ability to clearly indicate poor model performance (Gupta et al. 1999). PBIAS values for streamflow tend to vary more, among different auto-calibration methods, during dry years than during wet years (Gupta et al. 1999). This fact should be considered when attempting to do a split-sample evaluation, one for calibration, and one for validation.

RSR incorporates the benefits of error index statistics and includes a scaling/normalization factor, so that the resulting statistic and reported values can apply to various constituents. RSR varies from the optimal value of 0, which indicates zero RMSE or residual variation and therefore perfect model simulation, to a large positive value. The lower the RSR, the lower the RMSE, and the better the model simulation performance.

According to Moriasi et al. (2007), if model simulation has a NSE > 0.50, RSR < 0.70, and PBIAS ± 25% this is satisfactory. Ramanarayanan et al. (1997) suggested that if R² > 0.5 or NSE > 0.4, the model simulates natural phenomenon well. Therefore, these criteria can be used to evaluate the results of calibration, validation, and simulation obtained in this study.

The SWAT model has been simulated using TRMM gridded rainfall data sets for the Brahmaputra River basin. Each simulation was independently calibrated and validated against one discharge station. Calibration was carried out by adjusting the parameters until a good match was obtained between calculated and observed flows for each rainfall estimator. Then the same adjusting parameters were used for other rainfall estimators. Critical parameters identified during calibration were the curve number (CN2), base flow recession constant (ALPHA_BF), groundwater delay (GW_DELAY), minimum threshold depth of water required in shallow aquifer for ground-water flow to occur (GWQMN), and baseflow alpha factor for bank storage (ALPHA_BNK). The sensitivity ranks for SWAT parameters for the Brahmaputra River watershed are presented in Table 2.

A comparison of hydrographs for the calibration and validation periods between observed and model-generated flow using TRMM precipitation data is shown in Figure 3. Statistical indicators for evaluating model performance are shown in Table 3. R² values are found to be 0.88 and 0.69 for the calibration and the validation periods respectively. NSE values are found as 0.84 and 0.66 for calibration and validation periods, respectively. RSR for simulated discharge using TRMM data was quite good (values are less than 0.5) for both calibration and validation periods. An RSR value less than 0.5 is reasonably accepted (Moriasi et al. 2007). Here, TRMM data have been found more accurate for the both calibration and validation as PBIAS values are less than ±25%.

Box and whisker diagrams of the changes of the mean seasonal discharge under RCP 8.5 scenarios forced by 11 different RCM model generated climate data for four periods (historical, 2020s, 2050s, 2080s) are shown in Figure 4. The percentage of changes of the monthly mean discharge in the future scenarios of the Brahmaputra River at the Bahadurabad station are shown in Figure 5. Due to the uncertainty of the climate predictions, as represented by the different climate scenarios, the predicted discharge of Brahmaputra River at Bahadurabad station also varies and depends on the climate scenarios.

From Figure 4, it is seen that the simulations show mixed results for the RCP8.5 scenarios of the 11 RCMs. The dry period (DJF) shows no significant change up to year 2099. The variations observed are mostly negligible. However, from Figure 5, it can be seen that in a few cases the monthly changes during the winter season exceed 40%. As the flow in this season during the base period is so low, a small change in the flow will show greater change in percentage. During the pre-monsoon period (MAM), some of the models show significant increases of the discharge peaks, while most of the models show that the peak during this season will remain relatively unchanged. It is seen from Figure 4 that monthly mean flow in the pre-monsoon period varies from −12% to 29% in the 2020s (2010–2039). In the 2050s (2040–2069) and 2080s (2070–2099), the monthly mean flow varies from −3% to 24% and −12% to 31%, respectively.

The flows in the monsoon period (JJAS) however, show mixed results. Six (model no. 3, 5, 8, 9, 10 and 11 of Table 1) from among the 11 simulations show that the peak discharge, as well as their median, will increase as time passes throughout the century. Two of the models (model no. 10 and 11 of Table 1) show discharges exceeding 100,000 cms during the 2080s. On the other hand, models 1, 2, 4, 6 and 7 of Table 1 show discharges which are relatively unchanged with a decreasing tendency over time. Three of these models show a relatively low rise up to the middle of the century then the flow starts falling towards the end century. From Figure 4, the change of monthly mean flow during monsoon (JJAS) is found to be varying randomly. It also indicates a gradual increase of monsoon flow towards the end of the century for most of the RCMs. It is seen that in the 2020s (2010–2039) the mean monthly flow in monsoon will vary from −12% to 14%. In the 2050s and 2080s, mean monthly flow in this season will vary from −10% to 22% and −19% to 36%, respectively.

The changing flow pattern during the post-monsoon period (ON) is similar to the monsoon season (JJAS). The models that show increasing tendency during JJAS also show the same during ON. The same is true for the other models as well. The monthly flow during the post-monsoon season is following a gradual increase in most of the RCMs, as indicated by Figure 4. Figure 5 indicates that, in the 2020s the mean monthly discharge in monsoon will vary from −11% to 17%. In the 2050s and 2080s, it will vary from −16% to 17% and −17% to 31%, respectively.

A SWAT model has been setup over the ungauged Brahmaputra basin using the satellite-based DEM product, land use and soil information. The model has been successfully calibrated and validated using the gridded precipitation product from the TRMM satellite and temperature data from the ERA-Interim reanalysis product. Although the model is able to capture the average flow of the basin, it underestimates the peak of the discharge during the flood years.

To assess the possible changes of the future flow of the Brahmaputra River, the model has been simulated, forced by the RCP8.5 scenarios of the 11 RCMs. The flow during the dry period remains almost the same throughout the 21st century, while the flow during the pre-monsoon period exhibits an increasing tendency towards the end of the century. Most of the ensemble agreed regarding the increment of change of monthly flow from June to December throughout the end of the century. However, there might be major uncertainty of the changes of future flow during March, April, and May. The inter-quartile ranges of the change of flow during March, April and May have been found to be higher than other months of the year. On the other hand, the inter-quartile ranges of the change of flow during JJAS are significantly low, indicating less uncertainty about increasing future flow. The amount of monsoon flow will be always more prominent than the amount of pre-monsoon flow. The probability of the range of uncertainty in pre-monsoon flow will be more dominant than the probability of the range of uncertainty in monsoon flow. Having an increasing trend of monsoon flow through the century is very threatening for the lower riparian country of Bangladesh considering the severity of recent floods.

Although this study is able to show the opportunity for using SWAT as a policy advisory tool for assessing impact of climate change on water resources of the Brahmaputra basin, the performance of the model can be further improved. Availability of the hydro-meteorological data from ground-based observatories is essential to improve the performance of the model during calibration and validation. High resolution DEM, soil maps, and land use data will further improve the quality of the assessment. High resolution global circulation models (GCMs) can resolve the sub-grid scale processes which can improve the regional climate prediction of rainfall over this region.

The research leading to these results has received funding from the European Union Seventh Framework Programme FP7/2007-2013 under grant agreement no. 603864 (HELIX: High-End cLimate Impacts and eXtremes; http://www.helixclimate.eu).