Curve number modifications and parameterization sensitivity analysis for reducing model uncertainty in simulated and projected streamflows in a Himalayan catchment

Ecological Engineering 108 (2017) 17–29

Contents lists available at ScienceDirect

Ecological Engineering journal homepage: www.elsevier.com/locate/ecoleng

Research paper

Curve number modifications and parameterization sensitivity analysis for reducing model uncertainty in simulated and projected streamflows in a Himalayan catchment Vishal Singh, Manish Kumar Goyal

MARK

⁎

Dept of Civil Engineering, Indian Institute of Technology, Guwahati 781039, India

A R T I C L E I N F O

A B S T R A C T

Keywords: Curve number Streamflow SWAT model Model optimization CMIP5 GCMs

Climate essentially controls the supply of ecosystems, species ranges, and process rates on Earth. Modeling hydrological processes for hilly catchments dominated by snow cover and glaciers is complex and relies on improved calibration and uncertainty analysis methods to standardize the watershed based ecosystem management practices. In this study, a modified curve number (CN) approach has been employed to simulate streamflow and water yield at sub-catchment scale over hundred years utilizing Soil and Water Assessment Tool (SWAT), which relies on the main physical factors of the ecosystem such as landuse/landcover and soil. The model calibration and validation strength was evaluated using coefficients of determination (R2) and NashSutcliffe Equation (NSE) objective functions. The uncertainties (percentage) occurred in the modeled streamflow were estimated using two optimization algorithms such as the Sequential Uncertainty Parameter Fitting Approach (SUFI2) and the Parametric Solution (ParaSol). A parameterization based sensitivity analysis was carried-out to recognize the most influencing model calibration parameters. Statistical downscaling of daily temperature and precipitation datasets was performed utilizing Coupled Model Intercomparison Phase Five (CMIP5) Global Circulation Models (GCMs) with their Representative Concentration Pathway (RCP) experiments. The downscaled temperature and precipitation were utilized to assess the climate change impact on streamflows at sub-catchment scale. The historical and projected scenarios of streamflow (at the outlet) and water yield (at sub-catchment scale) showed substantial variabilities in their amount in both temporal (1991–2100) and spatial scales (sub-catchment 1 to sub-catchment 7). The magnitude of change analysis confirmed a substantial increase in the water yield across all the sub-catchments over Himalaya. The percent of change analysis ensured that the magnitude of change of water yield is highly vulnerable in the Himalayan catchments. The GCMs based projected scenarios demonstrated a consistent increase in the streamflow at both outlets (e.g. Lachung and Chungthang). The overall results show a consistent increase in precipitation and water yield amount over the Himalayan catchments. The variable streamflow in terms of amount and intensity may disrupt the ecosystem services in the Himalayan catchments.

1. Introduction Due to increments in the air temperature, the vapor carrying capacity is significantly affected across the world (Collins et al., 2013). The exponential rate of change of air temperature has been identified as one of the main causes of climate change (Taylor et al., 2012). The variations in the ratio of precipitation to snowfall and rainfall are increasing non-uniformly around the world (Chamaille-Jammes et al., 2007). A non-uniformity in the precipitation events and amount observed (Khadka et al., 2014). Extreme events will be more frequent in the Himalayan region in coming years (Shivam et al., 2016; Singh and

⁎

Corresponding author. E-mail addresses: [email protected], [email protected] (M.K. Goyal).

http://dx.doi.org/10.1016/j.ecoleng.2017.08.002 Received 23 March 2016; Received in revised form 24 July 2017; Accepted 2 August 2017 0925-8574/ © 2017 Elsevier B.V. All rights reserved.

Goyal, 2016; Singh et al., 2015). The Himalayan glacier ranges are the largest and most dynamic in the world (Pachauri et al., 2014; Collins et al., 2013). Several studies reported that the glaciers have been retreating since the end of the Little Ice Age (LIA) and they have been accelerated since 1970s (e.g. Collins et al., 2013). The earlier and accelerated melting rate of permanent snow covers and glaciers are contributing more water to the rivers (Collins et al., 2013). The earlier melting of snow-glaciers due to climate change have largely affected the seasonality and quantity of streamflow over Himalayan watersheds (Collins et al., 2013; Kulkarni et al., 2010). Many studies reveal that the Himalayan glaciers have been losing their mass and are shrinking

Ecological Engineering 108 (2017) 17–29

V. Singh, M.K. Goyal

sub-catchment spatial scales in the monthly time intervals. Schuol et al. (2008) found a significant variability in the discharge amount during historical and future time domains. In this study, the rainfall-runoff relationship based Soil Conservation Services (SCS) curve number (CN) method utilized for the computation of accurate streamflow discharge (Mishra et al., 2014; Neitsch et al., 2011). As model uncertainties reduce the confidence level of projection scenarios, the prime objective of this study was to implement the modified CNs in SWAT at a HRU scale by the computation of fractional slopes. The fractional slopes are computed at each hydrological response unit (HRU) scale and then the CNs are adjusted for each HRU to avoid the slope steepness, especially over hilly catchment (Jain et al., 2010). The deterministic models do not compute randomness (Singh et al., 2013a,b). This study incorporates stochastic advance optimization algorithms such as Sequential Uncertainty Parameter Fitting Approach (SUFI2) and Parametric Solution (ParaSol) (Abbaspour, 2011; Abbaspour et al., 2007). The results from both methods are compared to find out the suitability of the optimization methods for hilly catchments. The parameter based sensitivity analysis was done to find out the most influencing parameter over snow glacier induced Teesta river catchment (Singh et al., 2013a,b; Abbaspour, 2011). The optimization algorithms are found helpful to reduce various kinds of model uncertainties (e.g. parameter uncertainty, input data related uncertainties and model structure uncertainty) in calibrated and projected modeling outcomes (Singh et al., 2015, 2013a,b; Abbaspour et al., 2011). Another objective of this study was to evaluate the effects of climate change in the projected water yield and streamflows over the snowglacier induced Teesta river Himalayan catchment. The projection scenarios were generated utilizing the Coupled Model Intercomparison Project Phase Five (CMIP5) Global Circulation Models (GCMs) with their Representative Concentration Pathway (RCP) experiments (Taylor et al., 2012). The projected streamflow and water yield are characterized at each sub-catchment scale to highlight the spatial variations due to climate change. The annual surface runoff and water yield were analyzed at each sub-catchment scale in the future time domain to assess increase and decrease in their amount. The percentage of change analysis was conducted to determine the magnitude of change of water yield at each sub-catchment scale, which may be useful to highlight the stress over watershed ecosystem and their components in a Himalayan region even at smaller scales. The Teesta river catchment has very diverse ecology and thus the any change in streamflow and water yield with respect to glaciered and non-glaciered portions may affect the biodiversity of the region.

(Pachauri et al., 2014; Khadka et al., 2014; Krishna, 2005) and that the flow of Himalayan rivers is increasing (Singh et al., 2015). Thus the downstream portions of the rivers are flooding (Radic et al., 2014; Kulkarni et al., 2010). The eastern Himalayas is recognized as one of the most fragile ecosystems in the world. The poor and marginalized people living in this region face a pressing challenge of adapting to changing climate (Chettri et al., 2009). The shifts in rainfall and snowfall pattern over Himalayas due to climate change will affect the habitat of many species (Chakraborty et al., 2016). In particular, Chettri et al. (2009) and Huss (2011) conducted studies on the Himalaya and found that glaciers in this region have shrunk and the permanent snow cover is significantly decreased. The relative contribution of discharge from glaciers has decreased and the yearly peak discharge has shifted towards spring. Collins et al. (2013) worked on the climatic variations in Himalayan region and deduced that discharge of Ganges river, especially in Nepal, significantly fluctuated. However, the level of flow was usually preserved from the 1960s to 2000s. Bhutiyani et al. (2008) found that the discharge of the Satluj river significantly increased from the year 1991–2004 due to enhanced temperature warming. Bhambri et al. (2011) concluded that the snowmelt induced runoff will be higher in the 21st century over the Himalayan region and that could affect the aquatic ecosystems in Himalayas. Distributed and deterministic eco-hydrological models have been effectively utilized for the management of water resources. The physical and real time meteorological parameters based hydrological models were found more reliable for the computation of water balance than other modeling approaches (Singh et al., 2013a,b; Abbaspour, 2011; Arnold et al., 2012). The hydrological response unit (HRU) scale based semi-distributed hydrological model Soil and Water Assessment Tool (SWAT) and stochastic optimization tool SWAT Calibration and Uncertainty Program (SWATCUP) have already proven their utilities in various water balance studies around the world (Singh et al., 2013b; Abbaspour, 2011; Schuol et al., 2008; Gosain et al., 2006). To incorporate the present climate projections (e.g. Global Circulation Models: GCMs), we emphasized that the present modeling framework should be defined and explained (Li et al., 2015; Vermote, 2015; Slater et al., 2009; Molotch and Margulis, 2008). We highlighted the need to parameter optimization algorithms and uncertainty analysis to generate more accurate projection scenarios of main hydrological components such as streamflow, water yield and precipitation. The Teesta river in the eastern Himalayas is characterized by extremely high elevation ranges (e.g. 1400 m to 7400 m, approximately) (Singh and Goyal, 2016). The topography is dominated by snow cover, glaciers and steep slopes (Singh and Goyal, 2016). In previous SWAT based studies, several authors reported unusual high peaks in modeled streamflows, especially when modelling hilly terrain. The high peak in the compute discharge may be occurring due to steep slopes. It is found that this problem is very common in the curve number (CN) based runoff computation (Rahman et al., 2013; Schilling et al., 2008). The unreliable high flow peaks may enhance the uncertainty level, which further makes calculating of water balance unreliable (Singh et al., 2015; Mishra et al., 2014). This will tend to increase an uncertainty level in the projection of water balance components such as water yield and snowfall/snowmelt (Schulz et al., 2008). This could affect the downstream movement of streamflow and may affect the aquatic life of the living species. The SWAT model has proven its suitability for the calculation of in the projection of water balance components in large river basins across the world (Singh et al., 2015; Schuol et al., 2008; Gosain et al., 2006; Arnold et al., 2012). Gosain et al. (2006) used SWAT for simulating twelve large river basins in India to quantify the climate change impact on hydrology and concluded that the future discharge volume has been enhanced for those rivers originating from the Himalayan glaciers due to the higher melting rate of snow and glaciers. Schuol et al. (2008) used SWAT to analyze the hydrology of the entire Africa continent at

2. Study area description The Teesta river catchment (up to Chungthang gauge) part of Indus river basin has been selected for the analysis, which is a part of eastern Sikkim Himalaya India. The flow in the Teesta river originates from many small and large tributaries. Lachung river is one of the major tributaries of Teesta river joins the Teesta river at Chungthang gauge and also contributed a significant amount of flow in the downstream of Teesta river. The Teesta river originates from the Chhombo Chhu and Khangchun Chho glacial lakes. These lakes are situated at an elevation of 5280 m (Krishna, 2005). Many small tributaries meet the Teesta river on the eastern flank and several larger tributaries meets Teesta on the western flank. These tributaries contributed a large amount of drainage density. The Chungthang gauge point has been selected as the final outlet of the catchment (Fig. 1). The two most important meteorological factors i.e. air temperature and precipitation show significant variation from downstream to upstream portions. The temperature decreases with increasing height referred to temperature lapse rate (TLR) and precipitation gradient (changes in precipitation from one location to another) significantly vary over the Teesta river catchment (Singh and Goyal, 2016). The downstream portion corresponds to moderate elevation ranges (e.g. 18

Ecological Engineering 108 (2017) 17–29

V. Singh, M.K. Goyal

Fig. 1. Study area map.

precipitation datasets were adjusted at all the climate stations such as Chopta valley (SB1), Thangu (SB2), Muguthang (SB3) Lachen (SB4), Yumthang (SB5), Lachung (SB6) and Chungthang (SB7) by considering the topographical variations (Fig. 1). Each climate station represents local climates (e.g. temperature and precipitation lapse rates) and topographical variations (e.g. elevation variations) across the whole selected study area (Singh and Goyal, 2016). These datasets have been adjusted spatially at each climate station utilizing the methodology presented by Singh and Goyal (2016). The other meteorological parameters such as daily humidity, daily wind speed and daily net solar radiation are necessarily required for SWAT modeling. These variables have been downloaded from the SWAT weather data portal (http:// swat.tamu.edu/software/links/).

around 1400 m) and the upstream portion corresponds to high elevation ranges (e.g. around 7400 m) (Singh and Goyal, 2016). In winters, the upstream catchment part of Teesta river is mostly covered with permanent and temporary snow covers, and glaciers. The downstream portion of Teesta catchment receives a significant amount of rainfall during summer seasons and the river catchment receives around 75% precipitation during the monsoon season (June to October). The Teesta river catchment receives the highest precipitation in July and August, while the maximum snowfall is during January to March (Krishna, 2005).

3. Input data 3.1. Climate data

3.2. GIS data layers

The daily observed gridded (0.5° × 0.5°) precipitation and temperature datasets for the year 1980–2005 (around 26 years) have been used for the analysis. These datasets were collected from the Indian Meteorological Department and Indian Institute of Tropical Meteorology. These gridded datasets are prepared using more than 1800 observed precipitation and temperature gauges distributed across the India. These gridded precipitation and temperature datasets significantly utilized in various research and scientific studies (Shivam et al., 2016; Subash and Sikka, 2014; Singh et al., 2013a,b; Bhutiyani et al., 2008). These gridded datasets were utilized for the downscaling of daily temperature and precipitation at each sub-catchment scale. Additionally, the measured daily precipitation datasets (1980–2005) are also available at Lachung and Chungthang stations and thus utilized for the analysis. The Teesta river catchment was divided into seven sub-catchments (SB) mainly to spread the grid points across the whole catchment. The observed gridded daily minimum and maximum temperature and

The physical thematic data-sets such as digital elevation model (DEM), landuse/landcover (LULC) and soil map have been used in SWAT modeling to define the drainage and physiographical characteristics of the selected study area. The Cartographic Satellite mission (CARTOSAT) derived DEM (30 m spatial resolution) has been used in the delineation of streams, catchments, slopes and stream linking path channels. The fractional slope at each HRU scale has been computed using DEM. The HRUs are defined as per the unique combination of LULC, soil and slope. The LULC map (at 1:50,000 scale) is obtained from the BHUWAN portal of Indian Space Research Organization (http://bhuvan.nrsc.gov.in/bhuvan_links.php). The soil map (at 1:2, 50,000 scale) of the selected study region was obtained from the SWAT Waterbase portal (http://www.waterbase.org/).

19

Ecological Engineering 108 (2017) 17–29

V. Singh, M.K. Goyal

4. Methodology used

et al., 2007). In this study, the Global Sensitivity Analysis (GSA) approach has been employed to rank the model parameters using p-value and t-stat tests (Singh et al., 2013a,b; Abbaspour, 2011). If the p-value corresponds to 0, then the parameter can be identified as the most sensitive parameter (Abbaspour, 2011; Chamaille-Jammes et al., 2007). In GSA approach, a significance test (t-stat) is evaluated based on the significance level “alpha” (α). The value of alpha is considered as 0.05 at 95% confidence interval (Abbaspour, 2011). If the parameter values are recorded greater than +1.96 or lower than −1.96, then a significant (p < 0.05) positive or negative trend can be expected, respectively (Singh et al., 2013b; Abbaspour, 2011). In this study, total 11 calibration parameters such as R__CN2.mgt (curve number), V__ALPHA_BF.gw (baseflow), V__GW_DELAY.gw (groundwater delayed flow), V__GWQMN.gw (minimum depth of groundwater aquifer), V__ESCO.hru (Soil evaporation compensation factor), V__SFTMP.bsn (snowfall temperature), R__SOL_AWC.sol (soil available water content), R__SOL_K.sol (soil hydraulic conductivity), V__CH_N2.rte (manning roughness coefficient of the main channel), V__CH_K2.rte (effective conductivity) and V__ALPHA_BNK.rte (bank erosion), relevant to surface water, groundwater and snowmelt are considered. A SUFI2 algorithm works based on the iterative process, which updates the parameter values in each iteration and thus a new updated value is resulted after each iteration. The R prior to parameter represents relative change and V prior to the parameter represents the absolute change in the parameter value (for details see Abbaspour (2011)). SUFI2 algorithm computes the parameter uncertainty to ensure the observed data fall into the 95% Prediction Uncertainty (95PPU) bands. In SUFI2, the uncertainty bands have physically meaningful range that can be occurred due to use data inputs such as LULC, soil, precipitation, temperature and other data. From first iteration to the next iteration, SUFI2 decreases the uncertainty in steps while monitoring the p-factor and r-factor for next several iterations. The pfactor and r-factor are very crucial in SWAT modeling to determine the model uncertainty level in the modeling outcome. The p-factor actually gives the information of the percentage of modeled data is bracketed by observed data. The r-factor defines the thickness of uncertainty bands (Abbaspour, 2011). The 95PPU bands defined by upper and lower uncertainty bands, which can be computed between 2.5% and 97.5% levels of the cumulative distribution of an output variable. The 95PPU bands compute the uncertainty using Latin Hypercube Sampling (LHS) method (Singh et al., 2013b; Abbaspour, 2011). The 95PPU prohibits 5% of the very bad simulations or outliers in both upper and lower bounds (Abbaspour, 2011). Theoretically, the value of p-factor ranges between 0–100% and r-factor ranges between 0 to infinity. The value of pfactor = 1 and r-factor = 0 corresponds to measured data equivalent to observed data. For the goodness-of-fit computation between the calibrated/validated and observed datasets, the two objective functions such as the coefficient of determination (R2) (Zhang et al., 2014) and Nash-Suttcliffe Equation (NSE) (Moriasi et al., 2007; Nash and Sutcliffe, 1970) objective functions have been applied. The equations are given below:

4.1. Soil and water assessment tool (SWAT) SWAT is defined as semi-distributed, deterministic, and continuous time scale based hydrological model, which is fully capable of simulating and forecasting various water balance components at Hydrological Response Unit (HRU) and sub-catchment scale. Around 2400 peer reviewed publications are available related to SWAT model applications (Neitsch et al., 2011). The modeling outcomes can be simulated on the daily, monthly and yearly time steps (Neitsch et al., 2011; Arnold et al., 2012). In SWAT modeling, each sub-catchment (SB) corresponded with main channel path and various HRUs. In this study, the modified curve number (CN) method has been used to enhance the accuracy of streamflow and other water balance components. In SWAT modeling, the main hydrological component such as streamflow simulated at each HRU scale and then the aggregation of streamflow was done at each SB scale. Then water is routed to the main outlet of the catchment through Muskingum routing method (Arnold et al., 2012). A modified CN equation, initially evaluated by Neitsch et al. (2011) and Mishra et al. (2014), has been implemented in SWAT modeling. The SWAT derived CNs are adjusted as per the respective HRU slope and LULC class. The CNs modifications have been done as per the following steps (Neitsch et al., 2011):

Qs =

(R − Ia)2 R − Ia + S

(1)

Here, Qs is defined as surface runoff (mm), R is defined as daily rainfall (mm), Ia is defined as the initial abstraction (e.g. infiltration and interception prior to runoff) and S is defined as the retention parameter for the LULC and soil.

1000 S = 25.4 ⎛ − 10⎞ ⎠ ⎝ CN

(2)

Here, CN stands for the curve number. Ia is the initial abstractions and it is generally similar to as 0.2S and the equation will be as:

Qsurf =

(Rday − 0.2S )2 (Rday + 0.8S )

(3)

Rday is the Rainfall of the day. Runoff will be occurred if Rday > Ia. The slopes are generated for each HRU and then a new modified curve number generated for each HRU as per the given equation (Mishra et al., 2014).

CN 2α =

CN 2*322.79 + 15.63(α ) α + 323.52

(4)

Here, CN2α is the moisture condition II CN adjusted for slope α is the fractional slope for each HRU. Here, the value of α varies from 0.14 to 1.4 (Mishra et al., 2014). Likewise, CN2 stands for moisture condition II CN and CN3 stands for moisture condition III CN. In SWAT model, the LULC classes are defined as Forest (FRST), Shrubs (SHRB), Snow cover (SNOW), Glaciers (ICE) and pasture or grazing land (PAST).

2

⎡ (Qm, i − Qm)(Qs, i − Qs ) ⎤ ⎢∑i ⎥ ⎣ ⎦ Maximize R2 = 2 ∑ (Qm,i − Qm) ∑ (Qs,i − Qs )2

4.2. Model calibration and validation The stochastic optimization analysis tool SWAT Calibration and Uncertainty Program (SWATCUP) has been used for the model calibration and parameterization sensitivity analysis (Abbaspour, 2011). The SUFI2 and ParaSol algorithms inbuilt in SWATCUP used for the optimization of parameter values by adopting the concept of aggregate parameter selection (Yang et al., 2007). The sensitivity analysis has been performed to recognize the importance of a specific parameter to the given process (Fig. 2). In sensitivity analysis, one can determine the behavior of the parameters as per the given objective function. The sensitivity of the parameter can be observed when the parameter values influence the defined objective function (e.g. R2 and NS) (Abbaspour

i

i

(5)

Q is a variable (e.g. discharge), m and s stand for measured and simulated, respectively. i is the ith measured or simulated data.

Maximize NSE = 1 −

∑i ∑i

(Qm − Qs ) i 2 (Qm, i − Qm)2

(6)

Here, Q is defined as a variable (e.g. discharge), m and s stand for measured and simulated, respectively and bar stands for average. The ParaSol method utilizes a global optimization criteria, which 20

Ecological Engineering 108 (2017) 17–29

V. Singh, M.K. Goyal

Fig. 2. Methodology chart.

The Representative Concentration Pathway (RCP) experiments explore an estimate of the radiative forcing (RF) concentration level into the atmosphere that could be expected in the 21st century (Taylor et al., 2012). Different RCP experiments (e.g. RCP2.6, RCP4.5 and RCP8.5) have been generated to highlight the climate extremity from the baseline scenario to an extremely high level scenario. The RCP8.5 categorized as the extreme highest emission scenario and thus it represents the maximum RF concentration level. Likewise, the RCP2.6 and RCP4.5 are the mitigation (baseline scenario) and intermediate RF scenarios (Singh et al., 2016; Taylor et al., 2012). A statistical downscaling model (SDSM) (Wilby et al., 2014) has been utilized for the downscaling of daily temperature and precipitation. To minimize the errors in downscaled GCM variables, the bias correction method which is already evaluated by Singh and Goyal (2016) has been applied. After the spatial interpolation of GCM predictors (1980–2005) at each climate station (or SB scale), the SDSM model has been setup for the downscaling of temperature and precipitation data (2006–2100) using multimodel CMIP5 GCMs such as ESM–2 M, CM2P1 and CM3 with multiple Representative Concentration Pathway (RCP) experiments (e.g. RCP2.6, RCP4.5 and RCP8.5) (Singh and Goyal, 2016; Shivam et al., 2016). SDSM adds or removes precipitation days randomly based on prevailing likelihood of events in each month. If a month is mostly wet, it has a greater chance of a precipitation day being added. On the other hand, if a month is comparatively dry it has a bigger chance of a precipitation day being removed. SDSM model is able to perform linear, exponential and logistic regression analysis for making a regression equation between the predicted (e.g. observed precipitation and temperature) and predictors (e.g. GCM variables) (Shivam et al., 2016; Wilby et al., 2014; Taylor et al., 2012). In this study, a linear regression method has been followed for the model calibration and future projection (2006–2100) of minimum-maximum temperature and precipitation. In SDSM, the Reanalysis climate datasets based predictors during the historical time (Taylor et al., 2012) used for the model calibration (see Sivam et al. (2016)). The bias correction analysis has been performed by adopting the guidelines presented in earlier research articles (Shivam et al., 2016; Singh et al., 2016). The bias corrected daily temperature and precipitation data sets, then utilized for the projection of water yield and streamflows in SWAT model.

can be optimized using objective functions (Abbaspour, 2011). In this method, an uncertainty analysis can be performed under the choice of two statistical methods (Abbaspour, 2011). ParaSol adopts Shuffle Complex (SCE-UA) algorithm, which minimizes a function for up to 16 variables (Duan et al., 1992). The SCE-UA algorithm used widely for watershed model calibration. This algorithm has been also used for various hydrological applications such as subsurface hydrology, soil erosion, remote sensing and land surface modeling (Abbaspour, 2011). The two major steps of ParaSol method are given below: 1. In this method, the modeled outcomes have been categorized into “good” simulations and “poor simulations” by defining a threshold value. The threshold value can be determined by Bayesian statistics using probability density region of the model results. 2. The prediction uncertainty built correspondingly from the good simulations 4.3. Downscaling of CMIP5 GCMs and projected scenarios The six GCM grids (2.5° × 2.5°) neighboring the selected Teesta river catchment (up to Chungthang gauge) have been selected for the downscaling of daily temperature and precipitation data. The downscaling has been done to incorporate the global effect of climate change at local scale, especially in two main climate influencing parameters such as air temperature and precipitation. For the downscaling, total 16 most relevant predictors (Singh and Goyal, 2016; Taylor et al., 2012) have been chosen to cover the circulation domains of the predictors (Fig. 1). The six GCM predictors are spatially interpolated at each climate station (point specific) and also corresponded to each SB. Therefore, for the point scale downscaling of selected meteorological variables (e.g. minimum-maximum temperature and precipitation) under the CMIP5 GCMs, an Inverse Distance Weightage Approach (IDWA) (Lu and Wong, 2008; Snell, 1998; Snell et al., 2000) has been performed. The IDWA method has already evaluated by various authors in similar studies across the world (Singh and Goyal, 2016; Kharin et al., 2013; Goyal and Ojha, 2010). The CMIP5 GCM variables (e.g. daily minimum and maximum temperature, precipitation, etc.) are characterized by the concentration level of radiative forcing (RF) in the atmosphere (Taylor et al., 2012). 21

Ecological Engineering 108 (2017) 17–29

V. Singh, M.K. Goyal

Fig. 3. (a) Comparison of SWAT simulated streamflow using single slope, without CN adjustment and adjusted CNs, (b) SUFI2 based calibrated and validated outcomes of streamflow (scalar) at Lachung and Chungthang gauge stations, and (c) ParaSol based calibrated and validated outcomes of streamflow (scalar) at Lachung and Chungthang gauge station.

22

Ecological Engineering 108 (2017) 17–29

V. Singh, M.K. Goyal

5. Results

per SUFI calibration results, the Chungthang gauge has performed better correlation with observed datasets than Lachung. This is because the model has been calibrated from upstream (e.g. Lachung) to downstream (e.g. Chungthang). The calibration results (Fig. 3c) for ParaSol method on a monthly time step show a good correlation between computed and observed datasets. The Chungthang gauge station results (Fig. 3c) showed a better prediction than Lachung station. Table 1 provided a comprehensive assessment of the optimization results obtained from both methods. As per ParaSol results (Table 1), the R2 is recorded fairly well during calibration and validation. However, the NSE coefficient is computed little lower in case of ParaSol during calibration and validation. The lower NSE illustrated that the result obtained from a PraSol method shows higher uncertainty in their extreme values, because NSE is more sensitive for the extreme peak values (Abbaspour, 2011). In the ParaSol method, a Chungthang gauge station secured better R2 and NSE than Lachung gauge on a monthly calibration step. The overall results from the ParaSol method found satisfactory but less predictive than SUFI2 method. The modeled uncertainties have been determined through the objective functions such as p-factor and r-factors resulted in Table 1. The p-factor and r-factor computed fairly well for both the methods, which demonstrates a moderate uncertainty level in the simulated flows (Table 1). Most of the calibration steps procured p-factor > = 0.55 (Table 2a and b). The modeling outcomes show that the simulated data are mostly falling within the uncertainty band (upper 95PPU), as already explained in the methodology. As per the daily and monthly computational outcomes, around 55% simulated data (in some cases it is more than 60%) sets are matched with the observed datasets, while as per the monthly basis observations, around 60% simulated discharge values are matched with the observed data values. The r-factor statistics showed an optimal uncertainty level in the simulated discharge. The uncertainty evaluation results, utilizing p-factor and r-factor, suggested that the predicted outcomes are fairly optimized and comparable to previous studies (Zhang et al., 2014; Schilling et al., 2008; Abbaspour et al., 2007). In this study, total 11 parameters which are found most relevant to the streamflow and snowmelt process as per the literature survey have been selected to evaluate model uncertainty issues. A total of 11 parameters were taken for SUFI 2 calibration process and 8 parameters were selected for the ParaSol based optimization process. The sensitivity of the modeling calibration parameters has been tested using global sensitivity optimization. The selection of the most relevant and influencing parameters has improved the accuracy of the modeling outcomes. During sensitivity evaluation, the NSE objective functions has been considered in both the optimization methods. The best fitted value of all calibration parameters was obtained after running multiple iterative processes. The fitted best coefficient values computed based on a monthly calibration have been shown in Table 2 for both the optimization methods. Out of 11 calibration parameters (Table 2), eight parameters such as R_CN2, V_GW_DELAY, R_SOL_K, V__ESCO.hru, R__SOL_AWC, V__CH_N2, V__CH_K2 and V__ALPHA_BNK are recorded as significant sensitive parameters for SUFI calibration and V__GW_DELAY.gw is recorded as significantly sensitive for ParaSol method. The description of each parameter has been explained in Table 2. For significant sensitive parameters, the p-values are recorded almost zero (< 0.05) and t-stat values are computed > +1.96 or < −1.96. These parameters show a significant response in the optimization of best coefficient values. The significant sensitive parameters help to optimize the best model parameter values, which found helpful to reduce the model uncertainty during simulation and projection of main water balance components. In both optimization methods, the parameters such as R_CN2 and V_GW_DELAY have found most significant sensitive parameters. Results (Fig. 4a and b) show the sensitivity of R_CN2 and V_GW_DELAY as per the given objective function such as NSE and R2. Results (Fig. 4a) shows

Initially, the streamflow was computed as per the single slope (or constant slope) method and it figured highly uncertain at both catchment outlets shown in regression plots (Fig. 3a). Based on a single slope, the computed R2 is found less corresponding to observed and simulated stream flows and the fitted trend line look highly uncertain (Fig. 3a). Thus, multiple HRU slopes were computed and incorporated into SWAT modeling. The multiple HRU slope method improved the correlation between observed and modeled streamflows (Fig. 3a). The streamflows showed scalar values because as per the data dissemination policy of the Central Water Commission India, the Teesta river comes under classified river basin and we cannot publish the original observed values. The R2 value is not much improved between observed and modeled streamflow at both stations. Several extreme peaks can be observed in the modeled discharge (Fig. 3a). As per the basic approach of SWAT, the slope factor cannot be incorporated with CNs and thus the effect of slope may be misunderstood in the streamflow computation. Therefore, the SWAT model based computing CNs for each HRU class have been adjusted. The fractional slopes are computed at each HRU scale as per the fractional area of HRU class as per the methodology explained. Fig. 3(a) clearly showed a significant improvement in the computation of streamflow in modified CNs based regression plots. In the modified CN method, the R2 value is significantly improved as per the comparison made between the daily observed and simulated discharge at both the outlets such as Lachung and Chungthang (Fig. 3a). The modified CN method significantly reduced the unreliable extreme high peaks in the SWAT simulated discharge and thus the R2 value is significantly improved. The occurrence of high peaks due to slope steepness makes modeling outcomes unreliable and the increases the model uncertainties. The calibration and validation have been performed on both daily and monthly time scales. A linear regression plots (Fig. 3b) show calibration/validation results on a daily time step of SUFI2 method. Table 1 shows the calibration results for SUFI2 and ParaSol method. In SUFI2 method, the NSE coefficient shows better correlation at Chungthang gauge than Lachung gauge. As per SUFI2 results, for Chungthang station, the NSE values are computed as 0.76 and 0.78 on a monthly time step during calibration and validation, respectively. Whereas, for Lachung station, the NSE values are computed as 0.61 and 0.0.83 on a monthly time step during calibration and validation, respectively. The Lachung is situated at the upstream of Chungthang and thus the uncertainty level can be expected higher at Lachung gauge. As per SUFI2 results, the performance of SWAT model is fairly recorded for the monthly time step. Around 300 simulations were performed to get the highly optimized parameter values in each iteration process. Total 4 iterations were performed to get the best optimization results. It has been illustrated that in a hilly terrain environment, the SUFI2 algorithm performs better during monthly time step than a daily time step. However, it has not recorded much variation between daily and monthly calibration results. The daily calibration results were not discussed in this study. Overall, as Table 1 Calibration and validation results of SUFI2 and ParaSol methods as per monthly basis analysis. Time step

Calibration/ Validation

Outlet

p-factor

r-factor

R2

NSE

SUFI2

Calibration

Lachung Chungthang Lachung Chungthang

0.67 0.62 0.57 0.55

0.31 0.26 0.39 0.49

0.61 0.76 0.83 0.78

0.46 0.48 0.64 0.79

Lachung Chungthang Lachung Chungthang

0.58 0.52 0.57 0.55

0.46 0.44 0.46 0.42

0.62 0.72 0.71 0.72

0.46 0.44 0.52 0.43

Validation ParaSol

Calibration Validation

23

Ecological Engineering 108 (2017) 17–29

V. Singh, M.K. Goyal

Table 2 Calibration results of SUFI2 and ParaSol methods on the monthly basis analysis. Time Step

Parameter Name

1st Iteration

Final Iteration

Sensitivity Results

Min Value

Max Value

Fitted Value

Min Value

Max Value

Fitted Value

t-stat

p-value

SUFI2

R__CN2.mgt V__ALPHA_BF.gw V__GW_DELAY.gw V__GWQMN.gw V__ESCO.hru V__SFTMP.bsn R__SOL_AWC.sol R__SOL_K.sol V__CH_N2.rte V__CH_K2.rte V__ALPHA_BNK.rte

−0.20 0.05 10.00 0.00 0.80 −3.00 −0.20 −0.80 0.00 5.00 0.00

0.20 0.40 50.00 2.00 1.00 5.00 0.40 0.80 0.30 130.00 1.00

0.02 0.12 21.87 1.33 0.86 0.26 −0.08 −0.40 0.17 109.13 0.87

−0.63 0.16 8.65 1.58 0.75 −2.03 −0.23 −0.96 0.13 82.09 0.67

0.33 0.15 30.57 2.20 0.83 0.04 0.05 −0.50 0.18 113.82 1.10

0.04 0.13 26.00 1.81 0.83 −1.49 −0.17 −0.59 0.16 82.36 0.73

6.28 1.94 −7.37 0.88 2.18 0.80 −2.05 −4.72 7.71 −16.76 3.99

0.00 0.05 0.00 0.38 0.03 0.42 0.04 0.00 0.00 0.00 0.00

ParaSol

V__GW_DELAY.gw V__GWQMN.gw V__ESCO.hru V__ALPHA_BNK.rte R_SOL_AWC.sol R_SOL_K.sol R__SOL_BD.sol V__SFTMP.bsn

8.00 1.50 0.74 0.66 −0.20 −0.95 −0.70 −2.00

30.00 2.20 0.82 1.09 0.05 −0.50 −0.30 0.04

21.47 1.80 0.80 0.79 −0.20 −0.89 −0.38 −1.77

8.64 1.58 0.75 0.70 −0.22 −0.95 −0.80 −2.02

30.00 2.20 0.80 1.00 0.05 −0.51 −0.30 0.04

28.80 1.65 0.78 0.71 −0.13 −0.94 −0.47 −0.41

29.23 0.69 0.37 0.57 0.81 −0.65 −1.55 0.77

0.00 0.49 0.71 0.57 0.42 0.51 0.12 0.44

(RCP4.5 and RCP8.5) (SFig. 1a–c). The precipitation scenarios showed a significant increase in precipitation amount from the year 2008–2009 (SFig. 1–c). The black trend line (SFig. 1a–c) shows the average of all RCPs and it clearly shows an increasing trend from the year 2008–2100. All GCMs show a consistent and progressive increase in the precipitation amount in the 21st century. Each GCM model and their RCPs have distinctive, which are consisted with different radiative forcing levels. The RCP8.5 shows maximum increase in precipitation across all the SBs and corresponded with the extreme radiative forcing level. Among all the SBs, SB1, SB2 and SB3 are maximally fluctuated their peaks. Several trend lines (SFig. 1a–c) show a very extreme fluctuations in their beaks, especially after 2050s. The projected trends (SFig. 2a) for streamflow simulated at the Chungthang and Lachung outlets. It shows the monthly average projection of streamflow at Lachung gauge. All months show a significant increase in streamflow amount except January and August (SFig. 2a). Among all the months, The April, May, September and October show maximum increase. The higher runoff amount has observed during these months, mainly due to higher contribution of snowmelt runoff. The maximum increase was observed in RCP8.5 experimental scenarios. The monthly average projection of streamflow at Chungthang gauge is also computed ((SFig. 2b). In this, most of the month shows a significant increase in streamflow amount except July and August. In July, the trend shows decrease in streamflow runoff, illustrating a shift in the monsoon precipitation pattern. The RCP4.5 and RCP8.5 based projected scenarios showed a maximum increase in streamflow and a significant fluctuation in their peaks are also observed. Among all the GCMs, the CM3 model based projected scenarios found most vulnerable to climate change. The SB wise intra-decadal time series (TS) variations (TS1: 2011–2040, TS2: 2041–2070 and TS: 2071–2100) in the projected water yield is calculated (SFig. 3a–c). The results show that the calculated water yield has consistently increased from TS2011-2040 to TS2071-2100 (SFig. 3a–c). The largest increase has been recorded in RCP8.5 water yield scenarios as it represents the extreme level changes. The highest difference in water yield has been computed between TS2011-2040 and TS2071-2100 among all the GCMs scenarios. To quantify the magnitude of change in water yield at each SB, the percentage (%) of change has been calculated. The percentage of change in the water yield has also been evaluated in intra-decadal TS sets (Fig. 6 and Table 3). The Red colour in the bar chart (Fig. 6) shows the variations between TS1 and TS3 and the Yellow colour in the bar chart

the sensitivity results of R_CN2 and V_GW_DELAY from SUFI2 method. The optimized coefficient ranges of R_CN2 and V_GW_DELAY mostly vary from −0.03 to 0.03 and −2 to 28, respectively. The parameter values corresponded to 0.55–0.62 and 0.58–0.62, respectively, for R_CN2 and V_GW_DELAY (Fig. 4a). Similarly, results (Fig. 4b) shows the sensitivity results of R_CN2 and V_GW_DELAY from ParaSol method. In this, the optimized coefficient ranges of R_CN2 and V_GW_DELAY mostly vary from −0.05 to 0.10 and 10–20, respectively. The parameter values mostly corresponded to 0.48–0.53 and 0.47–0.52, respectively, for R_CN2 and V_GW_DELAY (Fig. 4a). Results (Fig. 5) of the main water balance components such as precipitation (PRECIPmm), snowmelt (SNOWmm), surface runoff (SURQmm) and water yield (WYLDmm) displayed at sub-catchment (SB) wise on a yearly basis (1991–2005). During model simulation, the initial 2 years (1989–1990) were taken as the warm up period. These variable showed significant variations in their amount from SB1 to SB7 (Fig. 5), which highlights the variable topographical responses of main water balance components over the Teesta river catchment. The water balance components during historical time significantly varied from the year 1991–2005 (Fig. 5). The spatial variations in the amount of selected variables illustrating that the topography of Teesta catchment is highly variable from SB1 (upstream) to SB2 (downstream). The SB7 corresponded to a maximum amount of precipitation than other SBs, while snowmelt was induced maximally at SB1 than other SBs. The SB1 and SB3 have categorized with extreme elevation ranges and thus these SBs corresponded to maximum snowfall during winter and snowmelt during summer. The water yield was computed maximally at SB7. The projected scenarios (2006–2100) of the main water balance components such as precipitation, water yield and streamflows have been used to highlight the significant impacts of climate change over the Teesta river Himalayan catchment. The details about projection scenarios were provided in an additional Supplementary material (see Supplementary Figs. 1–3). The severity of changing climate has been presented as per the CMIP5 GCMs corresponded with low (RCP2.6) to extreme emission scenarios (RCP8.5). The temperature and precipitation downscale at each sub-catchment scale. The projected scenarios (2008–2100) of main balance components have been generated in SWAT based on the downscaled temperature and precipitation. In model projection, the initial two years (2006–2007) were taken as the warm up period. The annual projected scenarios of precipitation calculated for each SB using multiple GCMs and RCPs such as 2M (RCP2.6, RCP4.5 and RCP8.5), CM3 (RCP2.6, RCP4.5 and RCP8.5) and CM2P1 24

Ecological Engineering 108 (2017) 17–29

V. Singh, M.K. Goyal

Fig. 4. Sensitivity analysis results as per the R2 objective function for the most significant sensitive parameter such as R_CN2, V_GW_DEL from SUFI2 and ParaSol methods; (a) results from SUFI2 method and (b) results from ParaSol method.

the percentage of change in the water yield (Fig. 6). This shows the long term climate change effect on the watershed components. The RCP4.5 and RCP8.5 based projected water yield scenarios showed a consistent increase in their amount as these experiments are defined as the extreme experiments. The outcomes (Table 3) show that the downstream SBs are more vulnerable to climate change and thus computed maximum percentage of change among all the GCMs and RCP experiments. The downstream portion of the catchment corresponds to moderate elevation ranges and revealed a significant variation in precipitation as observed from the GCM based different RCP experiments (SFig. 1–c). The observations which have projected as per the CMIP5 GCMs strongly

shows the variations between TS2 and TS1. Based on the results, the Red colour in the bar chart (e.g. TS1 and TS3) shows the maximum variation in water yield. The minimum percentage of change has been recorded between TS1 and TS2 among all the TS sets for all the GCMs. The results of percentage of change analysis at each sub-catchment scale validated the outcomes of the projection of water yield and precipitation at each sub-catchment scale (SFig. 1a–c). Results (Table 3) showed that most of the SBs show an enormous increase (e.g. positive values) in the water yield except few observations. The results (Table 3) clearly indicated that the higher percentage of change (e.g. 49.67% at SB7) has been computed between TS1 and TS3, similarly highlighted in 25

Ecological Engineering 108 (2017) 17–29

V. Singh, M.K. Goyal

Fig. 5. Sub-catchment wise variation in main water balance components such as (a) precipitation (mm), (b) snowmelt (mm), (c) surface runoff (mm), and (d) water Yield (mm) simulated by SWAT model during historical time (1991–2005).

two optimization algorithms gave a broad specification of the parameters at sub-catchment and HRU scale. The sub-catchment scale wise distribution of parameters enhanced the calibration strategy (Abbaspour, 2011). Similarly, our results showed a significant improvement in the overall calibration strategy and thus a good NSE observed at both the outlets during validation (Table 1). In both methods, the SUFI2 performs better over snow and glacier induced Himalayan catchment than the ParaSol algorithm (Table 1). Zhang et al. (2014) and Singh et al. (2013a,b) revealed the same and found that the SUFI2 is performing better in reducing parameter based modeling uncertainties. The reduction in model uncertainties significantly enhanced the future projection of main water balance components. The influence of climate change is found to have a more projecting effect on the water yield while impact of LULC change affect the surface runoff conditions (Kundu et al., 2017; Tong et al., 2012). The CN and groundwater delayed flow parameters relevant to water yield and surface runoff were identified as the most significant parameters for the optimization processes (Fig. 4). Many few studies resulted the CN as a significant sensitive parameter in hilly catchments (Singh et al., 2015; Rahman et al., 2013). The sensitivity of CN can be expected due the presence of steep slopes in hilly catchments. The CN is a function of LULC and soil (Mishra et al., 2014), and thus effective management practices in hilly catchments will protect the additional surface runoff. The accurate estimation of CNs at smaller scale (or HRU scale) through hydrological modeling tools will provide a significant feedbacks to the eco-hydrologist, which could be helpful to make an effective basin management plan over hilly catchments. The other soil related parameters such as groundwater delayed flow (V__GW_DELAY.gw) also found significant sensitive parameter. The snowfall temperature (V__SFTMP.bsn) is recognized as the significant parameter to control snow accumulation in temperature index model (Fontaine et al., 2002). However, in this study, the snowfall temperature did not show their sensitivity in the snowmelt induced streamflow calibration process. The reason is that the snowfall temperature is optimized closer to the real value. This study presented the correlation between the most influencing climate factors (e.g. temperature and precipitation) and watershed

supports that the future climatic conditions can be more vulnerable and the water yield amount will be enhanced in the 21st century over Himalayan catchment. 6. Discussion The CN method extensively used in various rainfall runoff and agricultural watershed management studies but CN method found uncertain in the streamflow simulation for hilly catchments. Many authors reported inconsistencies in the streamflow simulation over hilly catchments using SWAT inbuilt CN method (Singh et al., 2015; Zhang and Pan, 2014; Rahman et al., 2013; Huang et al., 2006). Previous Jha et al. (2014) and Huang et al. (2006) concluded that the CN cannot account the effect of slope. We observed several extreme high peaks in the modeled discharge (Fig. 3a) and a poor correlation between the observed and modeled discharge values. Unlike to previous studies, the present study based results showed a significant improvement in streamflow simulation using modified curve numbers (Mishra et al., 2014) adjusted at HRU scale. The HRU based fractional slopes helped to increase the accuracy of streamflow simulation at smaller scale (Fig. 3a). The improved streamflow results over a very hilly terrain catchment clearly provide an advantage to the eco-hydrologists in similar topographical geographies across the world. The small scale accurate computation of streamflow and other water balance components would be advantageous to eco-planners to maintain the biodiversity of the region, especially over Himalayas. Singh et al. (2015) and Abbaspour (2011) concluded that the model uncertainties can enhance in extreme high peak flow and low flow conditions due to under and overestimation of simulated discharge. This will reduce the accuracy of long term prediction of hydrological components. We observed several high peaks and overestimation of streamflow in the simulated outcome in case of without slope condition. The uncertainty in the water yield computation can increase due to wrong parameterization in the calibration process (Singh et al., 2013a,b; Abbaspour, 2011). This study applied two advance optimization algorithms such as SUFI2 and ParaSol already evaluated by few researchers (Singh et al., 2015, 2013a,b; Abbaspour et al., 2007). These 26

Ecological Engineering 108 (2017) 17–29

V. Singh, M.K. Goyal

Fig. 6. Percentage of change (%) in water yield (an average of all sub-catchments) in intra-decadal time series (TS) durations; TS1: 2011–2040, TS2:2041-2070 and TS3:2071-2100.

27

Ecological Engineering 108 (2017) 17–29

V. Singh, M.K. Goyal

has depleted over the Himalayas and thus the additional melt water could be the reason of enhancing streamflow amount in Himalayan rivers. This study results revealed the same (Table 3). The observations (Fig. 6) clearly highlights that the water yield and streamflow significantly increased across the Teesta river catchment. The SWAT based simulated snowmelt, surface water runoff and water yield showed substantial changes across the whole catchment during historical time (Fig. 5). The CMIP5 GCMs based climatic projections significantly showed an increase in the runoff at each subcatchment scale (Fig. 6). The above results highlighted the severity of climate change in the watershed components. Shivam et al. (2016) conducted a study on Subansiri river Himalayan catchment utilizing CMIP5 GCMs and revealed that the precipitation will be enhanced over Himalayas during monsoon season. This could be a major cause for the enhancement of the streamflow in the Himalaya. The projected changes in streamflow and precipitation may alter the natural balance of watershed components and their ecosystem services in the Himalayan catchments (Shivam et al., 2016; Kadka et al., 2014). This will cause a serious threat to the Himalayan biodiversity. The high intensity precipitation events and increasing amount of runoff will be dangerous for the ecosystem planning and management in Himalayan catchments. The combined effects of altered streamflow and high precipitation will likely reduce the reproductive success of the catchment. These modeling assets will definitely provide a great scope to the eco-hydrological planners to build a watershed management plan in the data scarce region.

Table 3 Percentage of change analysis results in the intra-decadal time series (TS) (e.g. TS1 2011–2040, TS2 2041–2070 and TS3 2071–2100) Water Yield. General Circulation Models (GCMs)

Coupled Climate Model version 3 (CM3)

Earth System Model (ESM2 M)

Coupled Model version 2 Phase 1 (CM2P1)

Subcatchments (SB)

Water Yield (mm)

TS1-TS3

TS1-TS2

TS2-TS3

RCP2.6

SB1 SB2 SB3 SB4 SB5 SB6 SB7

2.89 4.00 1.70 5.96 5.57 −0.48 6.20

5.68 6.04 2.75 9.37 4.60 −0.51 3.51

−2.63 −1.93 −1.02 −3.12 0.93 0.04 2.60

RCP4.5

SB1 SB2 SB3 SB4 SB5 SB6 SB7

9.59 22.48 17.22 7.54 14.27 7.52 19.30

6.44 22.54 9.74 1.61 9.28 4.43 14.12

2.96 −0.05 6.82 5.83 4.57 2.96 4.54

RCP8.5

SB1 SB2 SB3 SB4 SB5 SB6 SB7

21.32 27.73 20.04 10.99 25.65 16.18 49.67

7.47 21.65 8.90 7.04 14.67 29.82 20.72

12.89 5.00 10.23 3.69 9.57 −10.51 23.98

RCP2.6

SB1 SB2 SB3 SB4 SB5 SB6 SB7

−2.15 2.88 0.02 −0.95 −1.74 −0.98 0.09

−0.37 1.33 0.94 −3.10 1.41 0.58 −1.62

−1.78 1.53 −0.92 2.22 −3.11 −1.56 1.73

RCP4.5

SB1 SB2 SB3 SB4 SB5 SB6 SB7

3.36 8.58 4.41 6.72 7.50 7.19 11.40

0.14 1.84 5.61 −0.85 2.17 0.48 4.16

3.21 6.62 −1.14 7.63 5.22 6.67 6.95

RCP8.5

SB1 SB2 SB3 SB4 SB5 SB6 SB7

13.80 14.14 14.75 16.35 12.79 5.89 20.90

3.79 5.60 5.72 5.67 2.00 1.93 11.48

9.64 8.08 8.54 10.11 10.57 3.88 8.45

RCP4.5

SB1 SB2 SB3 SB4 SB5 SB6 SB7

0.31 6.85 1.06 4.99 1.08 2.87 2.19

1.66 5.61 −0.27 8.02 3.24 −0.73 3.20

−1.33 1.18 1.33 −2.81 −2.09 3.63 −0.98

RCP8.5

SB1 SB2 SB3 SB4 SB5 SB6 SB7

10.12 8.50 6.93 16.82 7.28 7.41 6.78

0.09 5.13 3.26 5.17 5.74 4.33 2.86

10.02 3.21 3.56 11.07 1.46 2.95 3.81

7. Conclusion A hydrological model SWAT has been applied for the projection and evaluation of streamflows and water yield in a snow/glacier induced Teesta river Himalayan catchment. An advanced stochastic optimization tool SWATCUP inbuilt SUFI2 and ParaSol algorithms have been utilized for the model calibration, validation and parameter uncertainty analysis. The SWAT model performed well in a very complex hilly catchment, as all the objective functions (e.g. p-factor, r-factor, R2 and NSE) reached satisfactory values. The sensitivity analysis helped to identify the most influencing model calibration parameters, such as CNs and groundwater delayed flow. It was observed that in the extreme hilly catchment, sometimes the models fail to predict the actual streamflow. Our results clearly showed that the modified and adjusted CNs based on HRU slopes improved the overall simulation process and projected more accurate water yield and other water balance components. The uncertainties occurred due to extreme high flow peaks successfully minimized by the modifications of CNs. The HRU based computation of fractional slopes helped to obtain the accurate estimation of streamflow in this extreme elevation topography catchment. In this research, we evaluated the impact of climate change on streamflows in the historical and 21st century time utilizing most relevant CMIP5 GCMs with their multiple RCP experiments. Historical water balance components were simulated at each sub-catchment scale and the annual results of streamflow showed a significant increase in their amount at both outlets. The sub-catchment wise assessment of water yield also illustrated that their extremity will be enhanced in future time. The multiple RPCs experiments based projected scenarios showed the extremity of precipitation and water yield at each subcatchment. The monthly projected scenarios (2008–2100) of streamflows clearly show a significant increase in their amount, especially in March, April, October and November months, which will add more water into the valley stream and thus the spring discharge quantity will be increased in future. The RCP4.5 and RCP8.5 experimental scenarios were shown a maximum increase in stream flows. It revealed that the precipitation variability will be enhanced in coming years. The precipitation amount will be increased over the Himalayan terrains as projected by the GCMs. The above demanding assessment of eco-hydrological components such as water yield, precipitation and

ecosystems through the projection of water balance components. Our results showed that the precipitation amount and total water yield are significantly increasing over the Himalayan river catchment due to climate change (see SFig. 1). A very few precious studies are conducted in the same context and revealed the same (Singh et al., 2016; Wambura et al., 2015). Rahman et al. (2013) concluded that due to climate changes the discharge over snow-glacier induced catchments is increasing. Kulkarni et al. (2010) found that the permanent snow cover 28

Ecological Engineering 108 (2017) 17–29

V. Singh, M.K. Goyal

climate and land use changes on the water balance. Ecol. Eng. 105, 42–57. Li, L., Vrieling, A., Skidmore, A., Wang, T., Muñoz, A.R., Turak, E., 2015. Evaluation of MODIS spectral indices for monitoring hydrological dynamics of a small, seasonallyflooded wetland in Southern Spain. Wetlands 35 (5), 851–864. Lu, G.Y., Wong, D.W., 2008. An adaptive inverse-distance weighting spatial interpolation technique. Comput. Geosci. 34, 1044–1055. Mishra, S.K., Chaudhary, A., Shrestha, R.K., Pandey, A., Lal, M., 2014.. Experimental verification of the effect of slope and land use on SCS runoff curve number. Water Resour. Manag. 28 (11), 3407–3416. Molotch, N.P., Margulis, S.A., 2008.. Estimating the distribution of snow water equivalent using remotely sensed snow cover data and a spatially distributed snowmelt model: a multi-resolution, multi-sensor comparison. Adv. Water Resour. 31 (11), 1503–1514. Moriasi, D.N., Arnold, J.G., Van Liew, M.W., 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. Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models part I—a discussion of principles. J. Hydrol. 10 (3), 282–290. Neitsch, S.L., Arnold, J.G., Kiniri, J.R., Williams, J.R., 2011. Soil Water Assessment Tool Theoretical Documentation Version 9. Texas Water Resource Institute of Technical Report No. 406. Texas A and M University. Pachauri, R.K., Allen, M.R., Barros, V.R., Broome, J., Cramer, W., Christ, R., Dubash, N.K., 2014. Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Radic, V., Bliss, A., Beedlow, A.C., Hock, R., Miles, E., Cogley, J.G., 2014. Regional and global projections of twenty-first century glacier mass changes in response to climate scenarios from global climate models. Clim. Dyn. 42 (1–2), 37–58. Rahman, K., Maringanti, C., Beniston, M., Widmer, F., Abbaspour, K., Lehmann, A., 2013. Streamflow modeling in a highly managed mountainous glacier watershed using SWAT: the Upper Rhone River watershed case in Switzerland. Water Resour. Manag. 27 (2), 323–339. Schilling, K.E., Jha, M.K., Zhang, Y.K., Gassman, P.W., Wolter, C.F., 2008. Impact of land use and land cover change on the water balance of a large agricultural watershed: historical effects and future directions. Water Resour. Res. 44 (7). Schuol, J., Abbaspour, K.C., Srinivasan, R., Yang, H., 2008. Estimation of freshwater availability in the West African sub-continent using the SWAT hydrologic model. J. Hydrol. 352 (1), 30–49. Shivam, Goyal, M.K., Sarma, A.K., 2016. Analysis of the change in temperature trends in Subansiri River basin for RCP scenarios using CMIP5 datasets. Theor. Appl. Climatol. 1–13. Singh, V., Goyal, M.K., 2016. Analysis and trends of precipitation lapse rate and extreme indices over north Sikkim eastern Himalayas under CMIP5ESM-2M RCPs experiments. Atmos. Res. 167, 34–60. Singh, S.K., McMillan, H., Bardossy, A., 2013a.. Use of the data depth function to differentiate between case of interpolation and extrapolation in hydrological model prediction. J. Hydrol. 477, 213–228. Singh, V., Bankar, N., Salunkhe, S.S., Bera, A.K., Sharma, J.R., 2013b. Hydrological stream flow modeling on Tungabhadra catchment: parameterization and uncertainty analysis using SWAT CUP. Curr. Sci. 104 (9), 1187–1199. Singh, V., Goyal, M.K., Chu, X., 2015. Multicriteria evaluation approach for assessing parametric uncertainty during extreme peak and low flow conditions over snow glaciated and inland catchments. J. Hydrol. Eng. 21 (1). http://dx.doi.org/10.1061/ (ASCE)HE.1943-5584.0001217. Slater, A.G., Clark, M.P., Barrett, A.P., 2009. Comment on “Estimating the distribution of snow water equivalent using remotely sensed snow cover data and a spatially distributed snowmelt model: a multi-resolution, multi-sensor comparison” by Noah P. Molotch and Steven A. Margulis (Adv. Water Resour. 31 (2008) 1503–1514). Adv. Water Resour. 32 (11), 1680–1684. Snell, S.E., Gopal, S., Kaufmann, R.K., 2000. Spatial interpolation of surface air temperatures using artificial neural networks: evaluating their use for downscaling GCMs. J. Clim. 13, 886–895. Snell, S.E., 1998. Spatial interpolation of surface air temperatures using artificial neural networks: evaluating their use for downscaling GCMs. J. Clim. 13, 886–895. Subash, N., Sikka, A.K., 2014. Trend analysis of rainfall and temperature and its relationship over India. Theor. Appl. Climatol. 117 (3–4), 449–462. Tong, S.T., Sun, Y., Ranatunga, T., He, J., Yang, Y.J., 2012.. Predicting plausible impacts of sets of climate and land use change scenarios on water resources. Appl. Geogr. 32 (2), 477–489. Vermote, E., 2015. MOD09A1 MODIS/Terra Surface Reflectance 8-Day L3 Global 500m SIN Grid V006. NASA EOSDIS Land Processes DAAChttp://dx.doi.org/10.5067/ MODIS/MOD09A1.006. Wambura, F.J., Ndomba, P.M., Kongo, V., Tumbo, S.D., 2015. Uncertainty of runoff projections under changing climate in Wami River sub-basin. J. Hydrol.: Reg. Stud. 4, 333–348. Wilby, R.L., Dawson, C.W., Murphy, C., O’Connor, P., Hawkins, E., 2014. The statistical downscaling model- decision centric (SDSM-DC): conceptual basis and applications. Clim. Res. 61, 259–276. Zhang, S., Pan, B., 2014. An urban storm-inundation simulation method based on GIS. J. Hydrol. 517, 260–268. Zhang, Z., Lu, W., Chu, H., Cheng, W., Zhao, Y., 2014. Uncertainty analysis of hydrological model parameters based on the bootstrap method: a case study of the SWAT model applied to the Dongliao River Watershed, Jilin Province, Northeastern China. Sci. China Technol. Sci. 57 (1), 219–229.

streamflow at smaller scale (or sub-catchment scale) in both historical and future time domains would provide great assets to eco-hydrologists, especially over Himalayas. Based on the above finding, eco-hydrologists can develop efficient strategies to manage watershed based ecosystem services. Acknowledgments This present research work has been carried out under DST research project No. SB/DGH-66/2013 and financial support is gratefully acknowledged. We also thank to GFDL, New Jersey USA, Earth Explorer NASA, CWC New Delhi India, IMD Pune India and IITM Pune India for providing the required datasets for the completion of this research work. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ecoleng.2017.08.002. References Abbaspour, K.C., Yang, J., Maximov, I., Siber, R., Bogner, K., Mieleitner, J., Zobrist, J., Srinivasan, R., 2007. Modelling hydrology and water quality in the pre-alpine/alpine Thur watershed using SWAT. J. Hydrol. 333, 413–430. Abbaspour, K.C., 2011. SWAT-CUP4: SWAT Calibration and Uncertainty Programs—A User Manual. Swiss Federal Institute of Aquatic Science and Technology, Eawag. Arnold, J.G., Moriasi, D.N., Gassman, P.W., Abbaspour, K.C., White, M.J., Srinivasan, R., Jha, M.K., 2012. SWAT. Trans. ASABE 55 (4), 1491–1508. Bhambri, R., Bolch, T., Chaujar, R.K., 2011. Mapping of debris-covered glaciers in the Garhwal Himalayas using ASTER DEMs and thermal data. Int. J. Remote Sens. 32, 8095–8119. Bhutiyani, M.R., Kale, V.S., Pawar, N.J., 2008. Changing streamflow patterns in the rivers of northwestern Himalaya: implications of global warming in the 20th century. Curr. Sci. 95 (5), 618–626. Chakraborty, A., Joshi, P.K., Sachdeva, K., 2016. Predicting distribution of major forest tree species to potential impacts of climate change in the central Himalayan region. Ecol. Eng. 97, 593–609. Chamaille-Jammes, S., Fritz, H., Murindagomo, F., 2007. Climate-driven fluctuations in surface-water availability and the buffering role of artificial pumping in an African savanna: potential implication for herbivore dynamics. Austral Ecol. 32 (7), 740–748. Chettri, N., Eriksson, M., Jing, F., Mool, P.K., Sharma, E., Shrestha, A.B., Tse-ring, K., 2009. Climate Change Impacts and Vulnerability in the Eastern Himalayas. ICIMOD. Collins, D.N., Davenport, J.L., Stoffel, M., 2013. Climatic variation and runoff from partially-glacierized Himalayan tributary basins of the Ganges. Sci. Total Environ. 468, S48–S59. Duan, Q., Sorooshian, S., Gupta, V., 1992. Effective and efficient global optimization for conceptual rainfall-runoff models. Water Resour. Res. 28 (4), 1015–1031. Fontaine, T.A., Cruickshank, T.S., Arnold, J.G., Hotchkiss, R.H., 2002. Development of a snowfall-snowmelt routine for mountainous terrain for the soil water assessment tool (SWAT). J. Hydrol. 262 (1), 209–223. Gosain, A.K., Rao, S., Basuray, D., 2006. Climate change impact assessment on hydrology of Indian river basins. Curr. Sci. 90 (3), 346–353. Goyal, M.K., Ojha, C.S.P., 2010. Robust weighted regression as a down- scaling tool in temperature projections. Int. J. Glob. Warm. 2 (3), 234–251. http://dx.doi.org/10. 1504/IJGW.2010.036135. Huang, M., Gallichand, J., Wang, Z., Goulet, M., 2006.. A modification to the soil conservation service curve number method for steep slopes in the Loess Plateau of China. Hydrol. Process. 20 (3), 579–589. Huss, M., 2011. Present and future contribution of glacier storage change to runoff from macroscale drainage basins in Europe. Water Resour. Res. 47 (7). Jain, S.K., Tyagi, J., Singh, V., 2010. Simulation of runoff and sediment yield for a Himalayan watershed using SWAT model. J. Water Resour. Prot. 2 (3), 267. Jha, R.K., Mishra, S.K., Pandey, A., 2014. Experimental verification of the effect of slope, soil, and amc of a fallow land on runoff curve number. J. Indian Water Resour. Soc. 34 (2). Khadka, D., Babel, M.S., Shrestha, S., Tripathi, N.K., 2014. Climate change impact on glacier and snow melt and runoff in Tamakoshi basin in the Hindu Kush Himalayan (HKH) region. J. Hydrol. 511, 49–60. Kharin, V.V., Zwiers, F.W., Zhang, X., Wehner, M., 2013. Changes in temperature and precipitation extremes in the CMIP5 ensemble. Clim. Change 119 (2), 345–357. Krishna, A.P., 2005. Snow and glacier cover assessment in the high mountains of Sikkim Himalaya. Hydrol. Process. 19 (12), 2375–2383. Kulkarni, A., Rathore, V., Singh, B.P., Ajai, S.K., 2010. Distribution of seasonal snow covers in central and western Himalaya. Ann. Glaciol. 51, 123–128. Kundu, S., Khare, D., Mondal, A., 2017. Individual and combined impacts of future

29