Abstract

This study presents climate change impacts on streamflow for the Subarnarekha basin at two gauging locations, Jamshedpur and Ghatshila, using the Soil and Water Assessment Tool (SWAT) model driven by an ensemble of four regional climate models (RCMs). The basin's hydrological responses to climate forcing in the projected period are analysed under two representative concentration pathways (RCPs). Trends in the projected period relative to the reference period are determined for medium, high and low flows. Flood characteristics are estimated using the threshold level approach. The analysis of variance technique (ANOVA) is used to segregate the contribution from RCMs, RCPs, and internal variability (IV) to the total uncertainty in streamflow projections. Results show a robust positive trend for streamflows. Flood volumes may increase by 11.7% in RCP4.5 (2006–2030), 76.4% in RCP4.5 (2025–2049), 20.3% in RCP8.5 (2006–2030), and 342.4% in RCP8.5 (2025–2049), respectively, for Jamshedpur. For Ghatshila, increment in flow volume is estimated as 15.7% in RCP4.5 (2006–2025), 24.2% in RCP4.5 (2025–2049), 35.9% in RCP8.5 (2006–2030), and 224.6% in RCP8.5 (2025–2049), respectively. Segregation results suggests that the uncertainty in climate prediction is dominated by RCMs followed by IV. These findings will serve as an early warning for the alarming extreme weather events India is currently facing.

INTRODUCTION

A meticulous analysis of climate change impacts on the water resources for different radiative forcing is essential for establishing appropriate adaptive strategies and planning at the regional level (Tan et al. 2014; Krysanova et al. 2017). Climate models currently do not simulate the hydrological cycle at a sufficiently high resolution for describing the basin-scale hydrological impacts on anthropogenic climate change (Kundzewicz et al. 2018). Hence, integration of climate models with hydrological models is essential to analyse the hydrological responses at the catchment scale. Hydrological models are typically driven by the bias-corrected outcomes from the number of climate models (general circulation models (GCMs) or RCMs), fed with scenarios of future radiative forcing (RCPs). Distinct climate models may produce an entirely different projection of future climate due to model biases and internal variability (Kling et al. 2012). Therefore, an ensemble of climate models is preferred instead of a single model to facilitate the characterisation of structural uncertainty and to improve climate predictions. Furthermore, climate change impact assessment studies assimilate a cascade of uncertainties.

Two major types of uncertainties arise from the integration of hydrological models with the climate models: traditional uncertainties (Refsgaard et al. 2013) and climate prediction uncertainties (Hawkins & Sutton 2009). Traditional uncertainties arise from hydrological models due to uncertainties in the model inputs, model structure, and model parameters and the climate prediction uncertainties arise due to the model uncertainty, scenario uncertainty, and IV (Hawkins & Sutton 2009). Model uncertainty arises due to the response of different climate models to produce dissimilar changes in climate in the presence of the same radiative forcing. Scenario uncertainty arises due to imperfect knowledge of the external factors affecting the climate system (e.g. future emission of greenhouse gasses). Internal variability (IV) represents the natural variability of the climate system occurring in the absence of external forcing and contains processes inherent to the atmosphere, the ocean, and the combined ocean-atmosphere system (Deser et al. 2012). These uncertainties due to climate prediction do not equally contribute to the uncertainties in streamflow projections (Chawla & Mujumdar 2018). Therefore, segregation of uncertainty due to climate prediction is essential to identify the relative contribution of each factor. Several methodologies have been reported for the disintegration of climate associated uncertainties in hydrological sciences, including ‘fourth order polynomial fitting’ (Hawkins & Sutton 2009), ‘ANOVA’ (Yip et al. 2011), ‘maximum entropy theory’ (Lee et al. 2017), ‘polynomial fitting combined with initial condition method’ (Zhuan et al. 2018), and ‘cumulative uncertainty approach’ (Kim et al. 2019). ANOVA is a relatively simple and powerful approach for partitioning of uncertainty. As per Yip et al. (2011), ANOVA requires fewer assumptions than other methodologies. ANOVA does not assume IV as constant, like the methodology discussed by Hawkins & Sutton (2009). Unlike the methodologies proposed by Lee et al. (2017) and Kim et al. (2019), ANOVA can also quantify the uncertainty due to interaction of model and scenario. Overlooking the substantial interaction term, such an analysis could lead to a considerable impact on the interpretation of the data (Kim et al. 2019). Therefore, we have used ANOVA for the disintegration of climate associated uncertainties. Examination of such uncertainties is crucial for the mitigation of the hazardous impacts of climate change.

Climate change can notably alter the regional hydrological cycle via changes in various climatic elements (Gelete et al. 2019). The frequent heavy precipitation leads to a considerable increase in the frequency, duration and intensity of extreme events like floods, and the economic losses caused by extreme events may have dramatic impacts on society. Since the probability of occurrence of extreme events cannot be minimised, their impacts can be lessened through mitigation and adaptation measures. Hence, it is obligatory to identify how flood characteristics will alter in the coming decades. These characteristics (duration and volume) do not only have scientific applications but are also helpful in water management (Quesada-Montano et al. 2018). For short-term flood forecasting, machine learning methods are very effective (Fotovatikhah et al. 2018; Mosavi et al. 2018; Yaseen et al. 2019); however, methodologies are different for the projection of flood characteristics for future climate. A popular method for the depiction of the characteristics of hydrological extremes is the threshold level approach (Yevjevich 1967). The threshold level approach is an effective, simple, and well-established methodology and it has been widely used for drought analysis (Fleig et al. 2006; Van Loon 2015; Quesada-Montano et al. 2018).

To the best of our knowledge, the assessment of flood characterisation so far has focussed on the past and present-day climate (including under naturalized and human influenced conditions) and the future prediction of flood characteristics is still a rarity. Developing countries are more vulnerable to the detrimental impacts of floods due to low flood protection strategies and vulnerability of populations, therefore, such an analysis could be beneficial to countries such as India. To fill this gap, the current study provides a generalised framework for the comprehensive assessment of potential impacts of climate changes on the streamflow at regional scale. The Subarnarekha basin of India is selected as the case study to demonstrate the methodology.

Subarnarekha is an inter-state basin (Jharkhand, Orissa, and West Bengal) providing a life-line to about 15 million people. In recent years, the frequency of floods has increased to a substantial degree in Subarnarekha. Every year during the rainy season, the river causes large scale devastation in Orissa by submerging the downstream areas (Singh & Giri 2018). Hence, it is essential to determine how projected climate change scenarios will affect the hydrological processes in the Subarnarekha basin in the near future. For this purpose, the SWAT hydrological model is used. SWAT has been verified as an efficacious tool for evaluating climate change impacts on streamflow globally (Tan et al. 2014; Swain & Patra 2019). To evaluate the impact of climate change on the streamflow of the basin, we used high resolution dynamically downscaled projections obtained from the Coordinated Regional Climate Downscaling Experiment (CORDEX-East Asia (EA)) for under a medium stabilization scenario, i.e. RCP4.5, and a very high emission scenario, i.e. RCP8.5. RCP4.5 is typically considered as an optimistic case, while RCP8.5 as a business-as-usual case (Rajsekhar & Gorelick 2017). Two other emission scenarios are RCP2.6 (mitigation scenario) and RCP6 (medium stabilization scenario). RCP2.6 is a low forcing scenario that provides outcomes that are far from reality (Mora et al. 2013). Further, though both medium-stabilization scenarios, i.e. RCP4.5 and RCP6, are equally acceptable, RCP4.5 shows a better match of trends of the average annual CO2 emission growth rate than RCP6 (Peter et al. 2013; Sanford et al. 2014). Therefore, to cover the entire range of scenarios, we selected RCP4.5 and RCP8.5 (Van Vuuren et al. 2011). Furthermore, on CORDEX-EA, the datasets are available for RCP4.5 and RCP8.5 only. We examined the flood characteristics using the threshold level method. The traditional uncertainties (parametric, input data, and model structural uncertainties) are handled through the sequential and uncertainty fitting algorithm (SUFI-2) tool of SWAT-CUP (Soil and Water Assessment Tool–Calibration and Uncertainty Procedure) software. The quantification of uncertainties due to climate predictions in projected high, medium, and low flows is carried out through ANOVA.

The specific objectives of this study are: (1) evaluation of SWAT hydrological model at two gauging stations (Jamshedpur and Ghatshila) in the Subarnarekha basin; (2) the quantitative assessment of climate change impacts on streamflow by looking for robust trends in high, low and medium flows; (3) quantification of flood characteristics for projected flows; and (4) the segregation of uncertainty due to climate change predictions.

METHODS

Study area

The Subarnarekha basin (21°15′ to 23°34′ N and 85°80′ to 87°32′ E) drains an area of 19,296 km2. Subarnarekha is a rainfed river that flows through the three states of Eastern India, before its confluence to the Bay of Bengal. The climate in the basin is tropical with hot summers and mild winters (average temperature varying from 9 to 41 °C), and the mean annual rainfall is 1800 mm. A significant amount of annual flow spreads over four wet months (June–September) owing to the south-west monsoon. The dominant land-use types of the basin are agriculture and forest, and the major soil type is loamy. Figure 1 presents the index map of the study area along with the area drained by each gauge.

Figure 1

Index map of the study area.

Figure 1

Index map of the study area.

Data

Geospatial data

The topographic information was obtained from Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) at 30 m spatial resolution. DEM data was further used for watershed delineation, extraction of flow direction, and calculation of sub-basin attributes. A landuse map of the basin was prepared from the unsupervised classification of Landsat-TM imagery (30 m resolution) of 2011. Soil maps of 90 m resolution were accessed from the Food and Agriculture Organization. In order to match with the scale of all geo-spatial datasets, the soil map was interpolated to 30 m using bilinear interpolation.

Hydro-meteorological data

The rainfall (1980–2014) at 0.25 × 0.25° resolution and maximum and minimum temperatures (1980–2014) at 1 × 1° resolution were acquired from the India Meteorological Department (IMD), Pune. The datasets for minimum and maximum temperatures were interpolated to 0.25 × 0.25° resolution using bilinear interpolation. Daily streamflow data was obtained from the Central Water Commission (CWC), Bhubaneswar for Jamshedpur and Ghatshila stations.

ERA-Interim reanalysis data

Reanalysis data in CORDEX-EA, which is available as ERA-Interim (ERAINT), was taken from the European Centre of Medium-Range Weather Forecast for testing of regional models with observational data. The initial and boundary conditions for ERAINT data were obtained from the National Centres for Environmental Prediction reanalysis. Reanalysis data was produced by combining observations with models and was often used as a proxy for observations where data were not available (Wehner 2013). ERAINT data is available on CORDEX-EA for 1989–2007. The details of ERAINT data are provided in Table S1 in the supplementary material.

RCP scenario data

For climate change impact analysis, future climate variables (maximum and minimum temperatures and precipitation) were obtained from all four RCMs, namely, HadGEM3-RA, RegCM4, WRF, and RSM from CORDEX-EA experiments. CORDEX–EA datasets are publically available at http://cordex-ea.climate.go.kr/cordex/.

The initial and boundary conditions for historical RCM data were obtained from the Atmosphere-Ocean coupled Hadley Centre Global Environmental Model (version-2) (HadGEM2-AO). For each model, data for different variables belonging to two RCPs, RCP4.5 and RCP8.5, emission scenarios were extracted for 20th-century climate (historical runs; 1981–2005) and future climate (2006–2049). These RCMs were selected as a large number of studies (Park et al. 2008; Yhang & Hong 2008; Cha et al. 2011; Baek et al. 2013; Gu et al. 2018) have shown that each model can reasonably reproduce the regional climate over the EA region. Table S1 presents the details of the reanalysis and RCP scenario data used in this study. All RCM outputs were brought down to consistent resolution (0.25 × 0.25°) using bilinear interpolation. In the current study, the ensemble median of RCMs projection was used for the simulation of future streamflow, as the multi-model ensemble (MME) median performs better than the MME mean (Knutti et al. 2013). As per Kundzewicz et al. (2018), the ensemble median was a more outlier-robust measure than the ensemble mean. Relatively few ensemble works have been performed for RCMs because of the lack of long-term simulations with multiple RCMs.

Bias correction of rainfall and temperatures

RCM simulations of temperature and precipitation often include significant biases. The reasons for such biases include systematic model errors due to faulty conceptualisation, discretisation, and spatial averaging (Teutschbein & Seibert 2010). Hence, the bias correction of RCM data is necessary to bring it closer to the observation data.

In the current study, RCM outputs (rainfall, and maximum and minimum temperatures) were bias-corrected with IMD gridded data (i.e. considered as observed data) through a distribution mapping method (Wood et al. 2004; Li et al. 2010; Piani et al. 2010) as described in Dhage et al. (2016) and Salvi et al. (2013). The distribution mapping method was chosen because Teutschbein & Seibert (2012) compared the performance of six bias correction methods and found that ‘distribution mapping’ outperformed all other methods. Precipitation correction was obtained by fitting Gamma distribution while temperature correction was obtained by fitting Gaussian distribution. Since RCMs simulate a large number of days with low precipitation instead of dry conditions, the precipitation distribution was usually distorted. Therefore, a precipitation threshold was introduced for frequency correction of rainy days before the distribution mapping.

To carry out bias correction using the distribution mapping method, the length of RCM data was divided into training (1981–2000) and testing (2001–2005) periods. For reanalysis data (1989–2005), 1989–2000 was used for training and 2001–2005 was used for testing. For the projected climate scenarios under RCP4.5 and RCP 8.5, the data ranging from 2006 to 2049 was divided into 2006–2030 and 2025–2049, considering six years, i.e. 2025–2030, as the overlapping interval. Fidelity of reanalysis driven RCM data (post-bias correction) and observed climate variables are presented through the Taylor diagrams (Taylor 2001). These diagrams show how closely a pattern matches the observation through the 2-dimensional plots in terms of three statistics. These three statistics are the correlation; root mean square difference and the ratio of variance between two patterns.

Hydrological model

The SWAT model is a physically based, semi-distributed, process-based river basin model (Arnold et al. 2012). It divides the watershed into subbasins, and these subbasins are further divided into hydrological response units (HRU) by a particular combination of land use, soil, topography, and slope. SWAT simulates the hydrology of the basin in land phase, using water balance that consists of rainfall, infiltration, interception, percolation, evapotranspiration, and return-flow components at subbasin (or HRU) levels. Further, the flow generated at the subbasin (or HRU) is routed through the network of streams at the outlet. SWAT simulates the surface runoff using the Soil Conservation Service Curve Number method and peak runoff based on the modified rational formula (Poertner 1974). Evapotranspiration is estimated using the Hargreaves method (Hargreaves & Samani 1985). For routing, the Muskingum method (Gill 1978) is used. The detailed description of these methods is provided in Arnold et al. (1998). The SWAT model simulates the water balance through the following equation:
formula
(1)
where SWf is the final soil water content at time f, SW0 is the initial soil water content in day i, and Pi, Qi, ETi, Ri, and Gi are the total precipitation, surface runoff, evapotranspiration, percolation to vadose zone, and baseflow on day i, respectively.
For model calibration, SUFI-2 of SWAT-CUP is used. In addition to automatic calibration, SUFI-2 also performs sensitivity analysis as well as uncertainty analysis. SUFI-2 outperforms other tools in SWAT-CUP, i.e. generalized likelihood uncertainty estimation, parameter solution, and Markov chain Monte Carlo (Yang et al. 2008; Xue et al. 2014). Global sensitivity analysis uses a multiple regression approach to determine the sensitivity of each parameter as:
formula
(2)
where g represents the objective function, α is the regression coefficient, β is the coefficient of parameters, and b is the parameter. It further uses the t-test to identify the relative significance of each parameter using hypothesis testing, where the null hypothesis states that the parameter being investigated has no effect at all. Therefore, assuming a 5% significance level, a p-value for a particular parameter lower than 0.05 indicates that the null hypothesis can be rejected. A p-value lower than 0.05 shows that there is a 95% chance of a parameter being sensitive. The t-stat is the test statistic obtained from the t-test. For a given t-stat and corresponding degree of freedom (i.e. (n–1); where n is the number of observations), the p-value can be obtained from the t-table.

This algorithm maps all uncertainties (parameter, model, input, etc.) on the parameters (expressed as uniform distribution and ranges) (Abbaspour et al. 2015). SUFI-2 quantifies the robustness of calibration through the estimation of 95% prediction uncertainty (95PPU), p-factor, and r-factor. The p-factor is the fraction of measured data enveloped by the 95PPU band and ranges from 0 to 1, where 1 indicates total enveloping of the measured data within model prediction uncertainty (i.e. a perfect model simulation considering the uncertainty) (Abbaspour et al. 2015). For discharge, a value of p-factor >0.75 is acceptable, subject to the adequacy of input data. The r-factor, on the other hand, is the ratio of the average width of the 95PPU band and the standard deviation of the measured variable. A value of <1.5 is desirable for this index (Abbaspour et al. 2015).

In this study, model setup is done through ArcSWAT 2012. The model is calibrated and validated at two gauging sites, i.e. Jamshedpur and Ghatshila, using SUFI-2. The model is calibrated for 2002–2009, using 2000–2001 as the warm-up period, and validated for 2010–2013.

For evaluating the performance of the model, Nash–Sutcliffe efficiency (NSE) (Equation (3)) is primarily used, i.e. NSE is maximised through the objective function. In addition to NSE, some other statistics, such as coefficient of determination (R2) (Equation (4)), and percent bias (PBIAS) (Equation (5)) are used to provide more insights to match/mismatch between the observed and simulated datasets.

NSE estimates the relative magnitude of residual variance compared to the variance of observed data. It lies between –∞ to 1. The coefficient of determination (R2) describes the proportion of the variance in measured data explained by the model. Its value lies between 0 and 1, with higher values indicating less error variance, and values greater than 0.5 are considered acceptable. PBIAS estimates the average tendency of the simulated data to be greater or smaller than their observed counterparts. PBIAS ranges between –∞ to ∞, with 0 as an optimum value.
formula
(3)
formula
formula
(4)
formula
(5)
where Qo and Qs denote observed and simulated flows, respectively and Qo(avg) and Qs(avg) denote average observed and simulated flows, respectively.

Analysis methods

Projected changes in mean, high, and low flows

Anticipated changes in RCM simulated rainfall, minimum and maximum temperatures are analysed with respect to reference periods under two climate change scenarios: RCP4.5 (optimistic case) and RCP8.5 (business-as-usual). Further, the annual discharge variation coefficient, i.e. standard deviations of all daily flows divided by mean annual values, as compared to the reference period are presented through the heat map. A heat map is a 2- or 3-dimentional representation of data in which values are presented by colours.

Projected streamflow are compared with the observed ones through the flow duration curves (FDCs) for high, medium, and low flows. FDC analysis is useful to discern what percentage of time the river flow exceeds a certain value known to cause flood damage.

Quantification of mean, high, and low flows is carried out through the estimation of mean annual flow (Q50) and two runoff percentiles representing low flows (Q90) and high flows (Q10). Further, the trends in medium, high, and low flows are calculated between the median of annual values for future periods (2006–2030 and 2025–2049) and the median of annual values for the reference period (1981–2005). The statistical significance of trends is estimated at the 0.05 significance level by nonparametric Wilcoxon signed rank test.

Estimation of flood volumes for RCP scenarios

For estimation of surplus volumes for both RCP scenarios, the threshold level method (Yevjevich 1967; Quesada-Montano et al. 2018) is used. As per this method, flood characteristics can be derived from time series of observed or simulated hydro-meteorological variables using a user-defined threshold level. We have used the symmetrical threshold approach (Quesada-Montano et al. 2018), by considering the low and high flow percentiles (Q90 for low flows and Q10 for high flows) as a threshold obtained from the reference period (1981–2005). A flood occurs when the flow is above the high-flow threshold, and drought occurs when the flow is below the low-flow threshold for any particular month or series of months. The surplus (flood)/deficit (drought) volumes are obtained as the deviation of flow above/below these thresholds for any month or sum of deviation of flow if it prolongs for more than one month. The duration of flood/drought events is obtained as the consecutive number of months having flow above/below these thresholds.

Uncertainty analysis

The contribution of three sources of uncertainties (model uncertainty, scenario uncertainty, and internal variability) in the simulated future streamflow is obtained through ANOVA (Yip et al. 2011; Bosshard et al. 2013; Chawla & Mujumdar 2018). The difference between future simulated values and reference values (Equation (6)) is considered as the input in the ANOVA method for partitioning of uncertainty due to individual sources in Q50, Q90, and Q10 as:
formula
(6)
where and represent the projected and observed values for a given variable Y.
An ANOVA model is constructed with different contributing factors, i.e. RCMs (represented by α; with four RCMs represented as i = 4 levels), RCPs (represented by β; with two RCPs represented as j = 2 levels), and the error term (ε) representing the internal variability (Equation (7)). A non-linear term (αβ) shows the interaction between these two terms and is also considered a factor for uncertainty. Hence, ΔY is written as the summation of the overall mean (μ) due to contributing factors and interaction between factors.
formula
(7)
ANOVA partitions the total variance into components due to different sources of variation. Therefore, the total variance in simulated future streamflow is given by the sum of squares from contributions and interactions as follows:
formula
(8)
formula
(9)
formula

Here, (Equation (8)) presents the total sum of squares, and present the sum of squares from contributors and interactions, respectively.

Uncertainty is estimated by the ratio of variance contribution due to individual factors to total variance ( or ). Therefore, the variance contribution due to each term (Equation (9)) is estimated. Relatively few studies have considered the interaction of model and scenario while segregating the uncertainty in climate predictions.

RESULTS AND DISCUSSION

Bias correction

In the current study, bias correction was carried out on a daily scale for both reanalysis data (1989–2005) and RCM data (1981–2005 and 2006–2049) with respect to the observed data. Frequency correction of rainfall days for all RCMs was done at a monthly scale prior to bias correction. For the frequency correction of rainy days, observed data (1981–2005) was used as a benchmark and grid wise threshold values were obtained for respective RCMs. These threshold values were obtained by arranging the observed and RCM rainfall data in descending order, reading the corresponding value of RCM rainfall for the day below which observed time series show zero rainfall. The obtained threshold values for the respective grids were then used for the future time periods (2006–2030 and 2025–2049). To assess the suitability of the bias correction method, the regional rainfall cycle produced by bias-corrected time series was compared with the observed time-series obtained from IMD data. Figure S1 in the supplementary material shows the spatial plot for observed data, bias corrected data, and bias uncorrected data for 1981–2005 for RegCM4 for the month of August. It is clear from Figure S1 that bias corrected data is able to capture the pattern of the observed data. The results show that distribution mapping performs well in the bias correction of rainfall and temperature, which is consistent with the findings of Piani et al. (2010), Teutschbein & Seibert (2010), Teutschbein & Seibert (2012), Smitha et al. (2017), and Heo et al. (2019).

Spatial patterns of the bias-corrected reanalysis data (precipitation, maximum and minimum temperatures) for all models, including MME, and observations at the regional level are presented through Taylor diagrams (Figure 2). From Figure 2, it can be concluded that RCM simulations perform well for Tmax and Tmin having high correlation (0.80–0.99) with the observed data. In the case of rainfall, the simulated data are slightly far-off from observed data; nevertheless, a reasonably fair correlation of 0.40–0.70 with the observed data is obtained. A relatively low correlation in precipitation responses may be because climate models often fail to represent the statistical properties of observed climate variables such as precipitation (Refsgaard et al. 2013). Furthermore, the precipitation projections obtained from numerous climate models and numerous assumptions about the future, e.g. socio-economic developments and RCP, can be largely different (Kundzewicz et al. 2018)

Figure 2

Taylor diagrams for inter-comparison of ERAINT RCM and observed data: (a) precipitation, (b) maximum temperature, and (c) minimum temperature. The reference point at a unit distance and 0° from the origin indicates the perfect performance. The distance to the reference point indicates the normalised root mean square error between the modelled and observations. The radial distance from the origin specifies normalised standard deviations, and cosine of azimuthal positions shows the spatial correlation between the models and observations.

Figure 2

Taylor diagrams for inter-comparison of ERAINT RCM and observed data: (a) precipitation, (b) maximum temperature, and (c) minimum temperature. The reference point at a unit distance and 0° from the origin indicates the perfect performance. The distance to the reference point indicates the normalised root mean square error between the modelled and observations. The radial distance from the origin specifies normalised standard deviations, and cosine of azimuthal positions shows the spatial correlation between the models and observations.

The results show that the multi-model ensemble of RCMs performs better than individual models for all climate variables. The results are consistent with the underlying principle that ensembles perform better in representing the current and future conditions (Kay et al. 2006; Park et al. 2008; Teutschbein & Seibert 2010; Smitha et al. 2017). Therefore, incorporation of ensemble simulation is essential in hydrological responses to reduce uncertainties.

Assessment of hydrological processes through SWAT model

The sensitivity of the simulated discharge to the model parameters was checked through the global sensitivity analysis performed using SUFI-2 of SWAT-CUP. Based on the results obtained from the global sensitivity analysis, eight parameters were found to be sensitive (p-value < 0.05). Table S2 presents the list of sensitive parameters having p-values < 0.05 along with t-stats. The SUFI-2 algorithm is used for the model calibration, validation, and uncertainty analysis. Optimum parameter values were obtained, after three iterations, with 1000 model runs in each iteration. As mentioned earlier, the model calibration and validation were performed at two gauging stations, i.e. Jamshedpur and Ghatshila. Figure 3(a)–3(d) presents the time series plot for corresponding observed and calibrated/validated SWAT simulated streamflow, and the scatter plot for observed and calibrated/validated SWAT simulated streamflow, at Jamshedpur and Ghatshila gauging stations, respectively. The values of performance measures for the calibration at Jamshedpur stations were NSE = 07, R2 = 0.72, and PBIAS = –8.9%, respectively. The negative value of PBIAS shows the underestimation of simulated flow at Jamshedpur station during calibration (Moriasi et al. 2007). The values of uncertainty measures for the calibration at Jamshedpur station were p-factor = 0.98 and r-factor = 1.17. For the validation at Jamshedpur station the values of performance measures were NSE = 60, R2 = 0.61, and PBIAS = 3% and uncertainty measures were p-factor = 0.87 and r-factor = 1.15, respectively. The positive values of PBIAS shows the overestimation of simulated streamflow at Jamshedpur station during validation. At Ghatshila station, the values of performance measures during calibration were NSE = 71, R2 = 0.72, and PBIAS = –14% (underestimation) and uncertainty measures were p-factor = 0.91 and r-factor = 1.1, respectively. However, during validation the values of performance measures NSE = 0.68, R2 = 0.69, and PBIAS = –5% (underestimation) and uncertainty measures were p-factor = 0.93 and r-factor = 1.2, respectively.

Figure 3

Calibration and validation of SWAT model: (a) time series plot for calibration/validation at Jamshedpur station, (b) scatter plot for calibration/validation at Jamshedpur station, (c) time series plot for calibration/validation at Ghatshila station, and (d) scatter plot for calibration/validation at Ghatshila station.

Figure 3

Calibration and validation of SWAT model: (a) time series plot for calibration/validation at Jamshedpur station, (b) scatter plot for calibration/validation at Jamshedpur station, (c) time series plot for calibration/validation at Ghatshila station, and (d) scatter plot for calibration/validation at Ghatshila station.

This shows that the performance of SWAT is satisfactory (evaluation measures: R2 > 0.6, NSE > 0.50, PBIAS < ±15 (as per Moriasi et al. 2015) and uncertainty measures: p-factor > 0.75, and r-factor < 1.50) (as per Abbaspour et al. 2015) for both calibration and validation periods for Jamshedpur as well as Ghatshila station.

Impact of climate change on streamflow

Analysis of climate variables under RCP scenarios

Figure 4(a) presents the MME projections for mean monthly rainfall under two RCP scenarios and the reference period. Projected rainfall showed a general increase during August under both RCP scenarios, though in RCP8.5 (2025–2049), the increasing trend was found during July and August. During the rest of the months, projected rainfall had a decreasing trend with respect to the reference period, except for RCP8.5 during 2006–2030. MME projections for rainfall show that peak rainfall will increase by 5.2% in RCP4.5 (2006–2030), 7% in RCP4.5 (2025–2049), 7.5% in RCP8.5 (2006–2030), and 25% in RCP8.5 (2025–2049), respectively. Further, MME projections for rainfall show that mean annual rainfall will decrease by 27% in RCP4.5 (2006–2025), 22% in RCP4.5 (2025–2049), 10% in RCP8.5 (2025–2049) and increase by 8.5% in RCP8.5 (2006–2030), respectively.

Figure 4

Comparison of mean monthly observed and projected climate: (a) rainfall, (b) maximum temperature, and (c) minimum temperature.

Figure 4

Comparison of mean monthly observed and projected climate: (a) rainfall, (b) maximum temperature, and (c) minimum temperature.

Figure 4(b) and 4(c) shows the MME projections for mean monthly minimum and maximum temperatures under two RCP scenarios and the reference period, respectively. Minimum temperature showed consistent trends with the reference period (e.g. increasing during January–June and decreasing during July–December). For the maximum temperature also, the projected data follow the reference period trend. MME projections for minimum temperature show an average increase of 2.4% (0.42 °C) in RCP4.5 (2006–2025), 5.3% (1.03 °C) in RCP4.5 (2025–2049), 5.0% (1.07 °C) in RCP8.5 (2006–2030), and 16% (3.38 °C) in RCP8.5 (2025–2030), respectively. MME projections for maximum temperature show an average increase of 1.5% (0.50 °C) in RCP4.5 (2006–2025), 2.6% (0.84 °C) in RCP4.5 (2025–2049), 3.5% (1.08 °C) in RCP8.5 (2006–2030), and 4.5% (1.40 °C) in RCP8.5 (2025–2030), respectively.

Figure S2 presents the coefficient of variation of annual streamflow as compared to the reference period in the form of a heat map. It ranges from –48 to –19%, with negative values implying that the flow variability in the projected period is relatively stable as compared to the reference period (Zhang et al. 2016). Under RCP4.5 and RCP8.5, the coefficient of variation with respect to the reference period is relatively low for 2006–2030, indicating less deviation (standard deviation) in flow probability due to less variability in precipitation and temperatures (Zhang et al. 2016) (Figure 4).

Projected changes in mean, high, and low flows

The flow duration curves (FDC) for Jamshedpur and Ghatshila gauging stations are drawn to compare the exceedance probability of daily streamflow, between the observed (1981–2005) and projected (2006–2049) climate to assess the flow patterns during 2006–2049 (Figure 5(a)). It is evident from the figure that under both emission scenarios the magnitude of streamflow increases, but the curve becomes relatively flatter during Q40–Q80 (medium flows). A flat curve generally indicates that groundwater contributions to the stream reach are significant which helps sustain the flow throughout the year (Chambers et al. 2017). FDCs further show a considerable increase in both high flows (10% exceedance) and low flows (90% exceedance) during 2006–2049 with respect to observed streamflows (1981–2005). Further, the Wilcoxon signed rank test is performed to analyse the trends in the mean (Q50), high (Q10), and low flows (Q90). Positive trends (statistically significant at 0.05 significance level) are found in Q10, Q50, and Q90 for 2025–2049 for RCP4.5 and RCP8.5 at both Jamshedpur and Ghatshila (Figure 5(b)). For 2006–2030, though positive trends are found for Q50, Q10, and Q90 at Ghatshila for both RCP4.5 and RCP8.5, these are not statistically significant. The same condition (positive trend but statistically insignificant) is also found for Q10 at Jamshedpur for 2006–2030 (RCP4.5). The positive trends obtained here in Q50, Q10, and Q90 under both RCP4.5 and RCP8.5 are consistent with the findings obtained from FDCs (Figure 5(a)). The results therefore show a high probability of floods in the Subarnarekha basin during 2025–2049.

Figure 5

(a) Comparison between exceedance probability of daily streamflow for observed and projected climate through flow duration curves (FDCs) at (i) Jamshedpur station, and (ii) Ghatshila station. (b) Quantification of difference in flow medians between projected and observed flows: (i) Q10 at Jamshedpur; (ii) Q50 at Jamshedpur; (iii) Q90 at Jamshedpur; (iv) Q10 at Ghatshila; (v) Q50 at Ghatshila; and (vi) Q90 at Ghatshila (I and II denote time periods 2006–2030 and 2025–2049, respectively).

Figure 5

(a) Comparison between exceedance probability of daily streamflow for observed and projected climate through flow duration curves (FDCs) at (i) Jamshedpur station, and (ii) Ghatshila station. (b) Quantification of difference in flow medians between projected and observed flows: (i) Q10 at Jamshedpur; (ii) Q50 at Jamshedpur; (iii) Q90 at Jamshedpur; (iv) Q10 at Ghatshila; (v) Q50 at Ghatshila; and (vi) Q90 at Ghatshila (I and II denote time periods 2006–2030 and 2025–2049, respectively).

Estimation of flood characteristics

The threshold level method is used in conjunction with the symmetrical threshold approach to estimate flood characteristics in the future period. Upon investigation of flood characteristics during the reference period 1981–2005, the method is able to capture the major flood in Subarnarekha during 1984 (Rakhecha 2002). This justifies the suitability of the threshold level method for flood characterisation in the Subarnarekha river basin. This finding is consistent with studies such as Kottegoda & Natale (1994), Zanchettin et al. (2008), Montanari (2012), Van Loon (2015), Baldassarre et al. (2017) and Quesada-Montano et al. (2018).

Consequently, the flow thresholds obtained from the reference periods are applied to future scenarios, and the volumes above/below these thresholds are estimated. In projected flows, no flow is found below the low flow percentile Q90, illustrating no risk of drought in the projected climate for the Subarnarekha basin. Figure 6 presents the magnitude of flood volumes and durations during the reference and the projected periods. The details of these attributes are given in Table S3. The results show that for a significant number of events flow exceeds the high flow percentile Q10 in both RCP4.5 and RCP8.5 scenarios, which justifies the high probability of occurrence of floods during the mid-century as obtained from the flow duration curves. Flood volumes may increase by 11.77% in RCP4.5 (2006–2030), 76.49% in RCP4.5 (2025–2049), 20.33% in RCP8.5 (2006–2030), and 342.49% in RCP8.5 (2025–2049) for Jamshedpur. For Ghatshila, the increment in flow volume is estimated as 15.75% in RCP4.5 (2006–2030), 24.20% in RCP4.5 (2025–2049), 35.92% in RCP8.5 (2006–2030), and 224.69% in RCP8.5 (2025–2030). These findings suggest that the increasing flood trends during the 2006–2049 period are going to follow the retrospective flood trends in SRB (Singh & Giri 2018).

Figure 6

Boxplots showing the distribution flood characteristic (surplus volumes and durations: (a) flood volume at Jamshedpur station, (b) flood duration at Jamshedpur station, (c) flood volume at Ghatshila station, and (d) flood duration at Ghatshila station (I and II denote time period 2006–2030 and 2025–2049 respectively).

Figure 6

Boxplots showing the distribution flood characteristic (surplus volumes and durations: (a) flood volume at Jamshedpur station, (b) flood duration at Jamshedpur station, (c) flood volume at Ghatshila station, and (d) flood duration at Ghatshila station (I and II denote time period 2006–2030 and 2025–2049 respectively).

Uncertainty analysis

ANOVA is used to decompose the uncertainty in the simulated future streamflows (Q10, Q50 and Q90) to identify the contribution of three sources of uncertainties, i.e. RCMs, RCPs, and IV. Figure 7 presents the results of the uncertainty analysis. As evident from the figure, RCMs are the primary source of uncertainty for all flows, followed by the internal variability. This suggests that uncertainty in the climate projections is primarily controlled by the selection of RCM. Our findings are consistent with a large number of studies where, in general, climate models are found to be the major contributor to uncertainty in future streamflow projections (Yip et al. 2011; Bosshard et al. 2013; Sellami et al. 2016; Krysanova et al. 2017; Chawla & Mujumdar et al. 2018; Kim et al. 2019). The magnitude of uncertainty due to IV implies its significant importance over the regional scales. Though IV is irreducible in nature, Smith et al. (2007) suggested that IV could be slightly reduced by appropriate initialization of climate predictions with observation datasets. Furthermore, interaction of the climate model and emission scenario (especially for high and medium flows) is found to be the next significant contributor to uncertainty. The contribution of the interaction of climate model and emission scenario to the uncertainty implies that the RCM and scenario should both be considered together (or in different combinations) for the future assessment of climate impacts on water resources (Yip et al. 2011; Bosshard et al. 2013; Chawla & Mujumdar et al. 2018). RCP is also found to be a contributor to uncertainty; however, its contribution is less in comparison to other sources.

Figure 7

Contribution of different sources to uncertainty in annual streamflow projections in: (a) high flows, (b) medium flows, and (c) low flows (I and II denote time period 2006–2030 and 2025–2049 respectively, while J and G denotes Jamshedpur and Ghatshila stations respectively).

Figure 7

Contribution of different sources to uncertainty in annual streamflow projections in: (a) high flows, (b) medium flows, and (c) low flows (I and II denote time period 2006–2030 and 2025–2049 respectively, while J and G denotes Jamshedpur and Ghatshila stations respectively).

SUMMARY AND CONCLUSIONS

This work seeks to address the two main objectives through the integration of RCM simulated outputs to the hydrological model: the first is to analyse the potential impacts of climate change on future streamflow simulations along with quantification of projected flood characteristics, and the second objective is related to quantification of uncertainty due to RCMs, emission scenario, and internal variability in streamflow projections. A semi-distributed SWAT hydrological model is used to carry out the analysis and to account for the spatial variability of streamflow. Model calibration/validation is carried out through the SUFI-2 algorithm. Relative to the reference period, increasing trends are found in future streamflow (low, high, and medium flows). The assessment of flood characteristics is carried out through the threshold level method by using symmetrical high/low flow thresholds. This method is found to be appropriate in the identification of flood characteristics for SRB. Upon the segregation of uncertainty, RCMs are found to be the major source of uncertainty in the streamflow projections followed by IV.

The major finding of the study is that the Subarnarekha basin may experience floods in the near future. The crux of the study, however, is that even though the results presented in this work may change quantitatively, we do not expect them to change qualitatively. Therefore, we suggest the incorporation of more RCMs (or GCMs) and emission scenarios for the analyses of climate change impact on SRB to improve confidence in the quantitative aspects as well.

These findings could provide solutions to the extreme weather events that India is currently facing, by providing an early warning. Further, the flood characteristics estimates would be helpful for the water managers and policy makers in forming adaptation measures for sustainable water resources management and disaster risk reduction.

SUPPLEMENTARY MATERIAL

The Supplementary Material for this paper is available online at https://dx.doi.org/10.2166/wcc.2020.254.

REFERENCES

Arnold
J. G.
Srinivasan
R.
Muttiah
R. S.
Williams
J. R.
1998
Large area hydrologic modeling and assessment part I : model development
.
J. Am. Water Resour. Assoc.
34
,
73
89
.
Arnold
J. G.
Moriasi
D. N.
Gassman
P. W.
Abbaspour
K. C.
White
M. J.
Srinivasan
R.
Santhi
C.
Harmel
R. D.
Griensven
A. v.
Liew
M. W. V.
Kannan
N.
Jha
M. K.
2012
SWAT:model use, calibration and validation
.
Trans. ASABE
55
,
1491
1508
.
Baek
H. J.
Lee
J.
Lee
H. S.
Hyun
Y. K.
Cho
C.
Kwon
W. T.
Marzin
C.
Gan
S. Y.
Kim
M. J.
Choi
D. H.
Lee
J.
Lee
J.
Boo
K. O.
Kang
H. S.
Byun
Y. H.
2013
Climate change in the 21st century simulated by HadGEM2-AO under representative concentration pathways
.
Asia Pac. J. Atmos. Sci.
49
(
5
),
603
618
.
Baldassarre
Di. G.
Martinez
F.
Kalantari
Z.
Viglione
A.
2017
Drought and flood in the Anthropocene: feedback mechanisms in reservoir operation
.
Earth Syst. Dyn.
8
(
1
),
225
233
.
Bosshard
T.
Carambia
M.
Goergen
K.
Kotlarski
S.
Krahe
P.
Zappa
M.
Schär
C.
2013
Quantifying uncertainty sources in an ensemble of hydrological climate-impact projections
.
Water Resour. Res.
49
(
3
),
1523
1536
.
Chambers
B. M.
Pradhanang
S. M.
Gold
A. J.
2017
Simulating climate change induced thermal stress in coldwater fish habitat using SWAT model
.
Water
9
(
10
),
26
28
.
Deser
C.
Phillips
A.
Bourdette
V.
Teng
H.
2012
Uncertainty in climate change projections: the role of internal variability
.
Clim. Dyn.
38
(
3–4
),
527
546
.
Dhage
P. M.
Raghuwanshi
N. S.
Singh
R.
Mishra
A.
2016
Development of daily temperature scenarios and their impact on paddy crop evapotranspiration in Kangsabati command area
.
Theor. Appl. Climatol.
128
(
3–4
),
983
997
.
Fleig
A. K.
Tallaksen
L. M.
Hisdal
H.
Demuth
S.
2006
A global evaluation of streamflow drought characteristics
.
Hydrol. Earth Syst. Sci.
10
(
4
),
535
552
.
Fotovatikhah
F.
Herrera
M.
Shamshirband
S.
Chau
K. W.
Ardabili
S. F.
Piran
M. J.
2018
Survey of computational intelligence as basis to big flood management: challenges, research directions and future work
.
Eng. Appl. Comput. Fluid Mech.
12
(
1
),
411
437
.
Gelete
G.
Gokcekus
H.
Gichamo
T.
2019
Impact of climate change on the hydrology of Blue Nile basin
.
Ethiopia: A Review. J. Water Clim. Change
1
12
.
doi: 10.2166/wcc.2019.014
.
Gill
M. A.
1978
Flood routing by Muskingum method
.
J. Hydrol.
36
,
353
363
.
Hargreaves
G. H.
Samani
Z. A.
1985
Reference crop evapotranspiration from temperature
.
Appl. Eng. Agric.
1
,
96
99
.
Hawkins
E.
Sutton
R.
2009
The potential to narrow uncertainty in regional climate predictions
.
Bull. Am. Meteorol. Soc.
90
(
8
),
1095
1107
.
Kay
A. L.
Jones
R. G.
Reynard
N. S.
2006
RCM rainfall for UK flood frequency estimation. II. Climate change results
.
J. Hydrol.
318
(
1–4
),
163
172
.
Kim
Y.
Ohn
I.
Lee
J. K.
Kim
Y. O.
2019
Generalizing uncertainty decomposition theory in climate change impact assessments
.
J. Hydrol.
3
.
doi: 10.1016/j.hydroa.2019.100024
.
Knutti
R.
Masson
D.
Gettelman
A.
2013
Climate model genealogy: generation CMIP5 and how we got there
.
Geophys. Res. Lett.
40
(
6
),
1194
1199
.
Kottegoda
N. T.
Natale
L.
1994
Two-component log-normal distribution of irrigation-affected low flows
.
J. Hydrol.
158
(
1–2
),
187
199
.
Krysanova
V.
Tobias
V.
Stephanie
E.
Shaochun
H.
Ilias
P.
Michael
S.
Alexander
G.
Rohini
K.
Valentin
A.
Berit
A.
Alejandro
C.
Ann van
G.
Dipangkar
K.
Anastasia
L.
Vimal
M.
Stefan
P.
Julia
R.
Ousmane
S.
Xiaoyan
W.
Michel
W.
Xiaofan
Z.
Fred
F. H.
2017
Intercomparison of regional-scale hydrological models and climate change impacts projected for 12 large river basins worldwide – a synthesis
.
Environ. Res. Lett.
12
(
10
),
105002
.
Kundzewicz
Z. W.
Krysanova
V.
Benestad
R. E.
Hov
Ø.
Piniewski
M.
Otto
I. M.
2018
Uncertainty in climate change impacts on water resources
.
Environ. Scie. Pol.
79
,
1
8
.
Lee
J. K.
Kim
Y. O.
Kim
Y.
2017
A new uncertainty analysis in the climate change impact assessment
.
Int. J. Climatol.
37
(
10
),
3837
3846
.
Mora
C.
Frazier
A. G.
Longman
R. J.
Dacks
R. S.
Walton
M. M.
Tong
E. J.
Sanchez
J. J.
Kaiser
L. R.
Stender
Y. O.
Anderson
J. M.
Ambrosino
C. M.
Fernandez-Silva
I.
Giuseffi
L. M.
Giambelluca
T. W.
2013
The projected timing of climate departure from recent variability
.
Nature
502
(
7470
),
183
187
.
Moriasi
D. N.
Arnold
J. G.
Liew
M. W. V.
Bingner
R. L.
Harmel
R. D.
Veith
T. L.
2007
Model evaluation guidelines for systematic quantification of accuracy in watershed simulations
.
Trans. ASABE
50
(
3
),
885
900
.
Moriasi
D. N.
Gitau
M. W.
Pai
N.
Daggupati
P.
2015
Hydrologic and water quality models: performance measures and evaluation criteria
.
Trans. ASABE
58
(
6
),
1763
1785
.
Park
E.
Hong
S.
Kang
H. Ã
.
2008
Characteristics of an East-Asian summer monsoon climatology simulated by the RegCM3
.
Meteorol. Atmos. Phys.
100
,
139
158
.
Peters
G. P.
Andrew
R. M.
Boden
T.
Canadell
J. G.
Ciais
P.
Le Quéré
C.
Marland
G.
Raupach
M. R.
Wilson
C.
2013
The challenge to keep global warming below 2°C
.
Nat. Clim. Chang.
3
(
1
),
4
6
.
Poertner
H
.
1974
Practices in Detention of Urban Storm Water Runoff
.
APWA Special Rep. No. 43
.
American Public Works Association
,
Washington, DC
, p.
206
.
Quesada-Montano
B.
Di Baldassarre
G.
Rangecroft
S.
Van Loon
A. F.
2018
Hydrological change: towards a consistent approach to assess changes on both floods and droughts
.
Adv. Water Resour.
111
,
31
35
.
Rakhecha
P. R.
2002
Highest floods in India
.
IAHS-AISH Publ.
271
,
167
172
.
Refsgaard
J. C.
Drews
M.
Jeppesen
E.
Madsen
H.
Markandya
A.
Olesen
J. E.
Porter
J. R.
Christensen
J. H.
2013
The role of uncertainty in climate change adaptation strategies – A Danish water management example
.
Mitig. Adapt. Strateg. Glob. Chang.
18
,
337
359
.
Sanford
T.
Frumhoff
P. C.
Luers
A.
Gulledge
J.
2014
The climate policy narrative for a dangerously warming world
.
Nat. Clim. Chang.
4
(
3
),
164
166
.
Singh
A. K.
Giri
S
.
2018
Subarnarekha River: The Gold Streak of India
.
The Indian Rivers
.
Springer
,
Singapore
, pp.
273
285
.
Smith
D. M.
Cusack
S.
Colman
A. W.
Folland
C. K.
Harris
G. R.
Murphy
J. M.
2007
From a global climate model
.
Paleoceanography
796
(
2007
),
796
799
.
Van Loon
A. F.
2015
Hydrological drought explained
.
Water
2
(
4
),
359
392
.
Van Vuuren
D. P.
Edmonds
J.
Kainuma
M.
Riahi
K.
Thomson
A.
Hibbard
K.
Hurtt
G. C.
Kram
T.
Krey
V.
Lamarque
J. F.
Masui
T.
Meinshausen
M.
Nakicenovic
N.
Smith
S. J.
Rose
S. K.
2011
The representative concentration pathways: an overview
.
Clim. Chang.
109
(
1
),
5
31
.
Wood
A. W.
Leung
L. R.
Sridhar
V.
Lettenmaier
D. P.
2004
Hydrologic implications of dynamical and statistical approaches to downscaling climate model outputs
.
Clim. Chang.
62
,
189
216
.
Yang
J.
Reichert
P.
Abbaspour
K. C.
Xia
J.
Yang
H.
2008
Comparing uncertainty analysis techniques for a SWAT application to the Chaohe Basin in China
.
J. Hydrol.
358
(
1–2
),
1
23
.
Yevjevich
V.
1967
An Objective Approach to Definitions and Investigations of Continental Hydrologic Droughts. Hydrology Papers
.
Colarado State University
,
Fort Collins, CO
,
USA
.
Yip
S.
Ferro
C. A. T.
Stephenson
D. B.
Hawkins
E.
2011
A simple, coherent framework for partitioning uncertainty in climate predictions
.
J. Clim.
24
(
17
),
4634
4643
.
Zanchettin
D.
Traverso
P.
Tomasino
M.
2008
Po River discharges: a preliminary analysis of a 200-year time series
.
Clim. Chang.
89
(
3–4
),
411
433
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY-NC-ND 4.0), which permits copying and redistribution for non-commercial purposes with no derivatives, provided the original work is properly cited (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Supplementary data