Abstract
The river situation in a dry or semi-dry area is extremely affected by climate change and precipitation patterns. In this study, the impact of climate alteration on runoff in Kashafrood River Basin (KRB) in Iran was investigated using the Soil and Water Assessment Tool (SWAT) in historical and three future period times. The runoff was studied by MIROC-ESM and GFDL-ESM2G models as the outputs of general circulation models (GCMs) in the Coupled Model Intercomparison Project Phase 5 (CMIP5) by two representative concentration pathway (RCP) scenarios (RCP2.6 and RCP8.5). The DiffeRential Evolution Adaptive Metropolis (DREAM-ZS) was used to calibrate the hydrological model parameters in different sub-basins. Using DREAM-ZS algorithm, realistic values were obtained for the parameters related to runoff simulation in the SWAT model. In this area, results show that runoff in GFDL-ESM2G in both RCPs (2.6 and 8.5) in comparing future periods with the historical period is increased about 232–383% and in MIROC-ESM tends to increase around 87–292%. Furthermore, GFDL-ESM2G compared to MIROC-ESM in RCP2.6 (RCP8.5) in near, intermediate, and far future periods shows that the value of runoff increases 59.6% (23.0%), 100.2% (35.1%), and 42.5% (65.3%), respectively.
INTRODUCTION
The impact of climate change on the availability of water resources has been attracting the attention of researchers throughout the world (Faramarzi et al. 2009; Rodrigues et al. 2014; Zuo et al. 2015). In the context of precipitation frequency changes and climate change, numerous studies indicate that global warming leads to ecosystem degradation and water crisis (Falkenmark & Rockström 2006; Zuo et al. 2015) and exacerbated water scarcity in dry and semi-dry regions. About 600 million people live in areas with less than 500 m3 water per capita, therefore enduring a severe shortage of water (Pereira et al. 2009). Thus, there is a great deal of concern among researchers about the increasing threat of water scarcity not only in the region but also at global scale (Vörösmarty et al. 2010). As a result, in the global warming context, comprehensive assessment of water resources is a vital part of understanding the availability of water and improving water management towards maintainable, effective, and instrumental usage of scarce freshwater resources. Furthermore, climate change has a significant impact on water availability in river basins and there is an interaction between water resources and climate change (Faramarzi et al. 2009). Spatial and temporal variabilities of water resources in basins are severely impacted by climate change (Vörösmarty et al. 2010). Due to the high variation and low rate of precipitation in dry and semi-dry areas, many sectors can suffer direct or indirect effects on the economy and can lead to loss of people (Anyamba et al. 2014) and, most likely, effects on runoff variation (Chiew et al. 2009). At present, general circulation models (GCMs) are the most frequently used models for the projection of different climatic change scenarios. The Intergovernmental Panel on Climate Change (IPCC) gathers and reviews global climate models as a part of the international climate change assessment reports for policymakers and public users (IPCC 1990). Up to now, IPCC has released five different versions of GCM. The IPCC Fifth Assessment Report of the Intergovernmental Panel on Climate Change (AR5) has introduced a new way of developing scenarios. Numerous studies have been reported about runoff assessment using representative concentration pathway (RCP) scenarios of AR5 on future climate conditions at global and regional scales.
Tan et al. (2017) used five GCMs under RCP2.6, RCP4.5, and RCP8.5 for two future periods (2015–2044 and 2045–2074) and concluded an increase of future annual runoff in Kelantan River Basin in Malaysia. Sun & Peng (2020) using ten GCMs of the CMIP5 under RCP4.5 estimated the future annual runoff in the north–south transect of eastern China for a future period (2011–2100). Yang et al. (2019) simulated precipitation and runoff in Luanhe River Basin in north China using GCMs under four RCPs (i.e., RCP2.6, RCP4.5, RCP6.0, and RCP8.5) for a future period (2020–2030). The result showed that runoff may increase between 39 and 58% compared to the observed runoff depth.
Iran, with a dry and semi-dry climate in the Middle East, is recognized to have low precipitation, high temperature, and high evapotranspiration demand. The water leveling is relatively fragile and endangered by water scarcity that can significantly exacerbate because of climate change (IPCC 2007). Arguably, climate change is one of the big challenges facing Iran; it has an especially harmful impact on water resources which creates problems for the environment and the economy, especially the agricultural sectors. There is a need to predict the potential climate change impacts on the interval and amount of precipitation that can be used for supporting and controlling water resources suitably and decreasing the water scarcity problem (Ternik et al. 2013). A semi-distributed hydrological model, Soil and Water Assessment Tool (SWAT), was used to supply more insights into the mechanisms of hydrological procedures and to assess hydrological processes and effects of climate change on water resources spatiotemporally (developed by United States Department of Agriculture (USDA; Arnold et al. 1998). In northeast Iran, the Kashafrood River Basin (KRB) is the longest river and supplies surface water for Mashhad city. Owing to its extreme dry weather, this basin has been suffering from water scarcity and pollution (Afshar et al. 2017). Up to now, water scarcity related to climate change in the KRB has not been well measured in the analysis of climate change and climate policy construction. Thus, the main purpose of this study was to assess the potential climatic changes on future runoff in KRB according to the CMIP5; Taylor et al. 2012).
The correct values of the parameters in the hydrological model (i.e., SWAT) cannot be directly measured (Li et al. 2010) and should be computed through a calibration process (Laloy et al. 2010; Laloy & Vrugt 2012; Leta et al. 2015; Schoups & Vrugt 2010). The calibration process in a modeling scheme plays an important role in extracting the optimal parameter values which accord reasonably with reality (Leta et al. 2016).
Formal and informal Bayesian frameworks have been developed to calibrate and analyze uncertainty in different river basin models. Formal Bayesian approaches apply likelihood with employing the real assumption to measure the total predictive uncertainty and the posterior probability densities of the parameters (Kuczera & Parent 1998; Schoups & Vrugt 2010; Muñoz et al. 2014; Rivera et al. 2015). The Markov Chain Monte Carlo (MCMC) algorithms can be mentioned as a good example of formal Bayesian approaches. DiffeRential Evolution Adaptive Metropolis (DREAM-ZS), as the most commonly used MCMC algorithm, originated from the DREAM (Vrugt et al. 2009a), which is used as a robust and effective sampler (Vrugt et al. 2009b; Laloy & Vrugt 2012; Laloy et al. 2012). The DREAM-ZS relies on the samplings from the past state and produces candidate points for individual chains. DREAM-ZS provides some advantages: (i) posterior sampling requires only three parallel chains, leading to reduced time for burn-in; (ii) requirement of less function assessments than DREAM, leading to suitable limiting dispensation; (iii) comprising a snooker updater which induces jumps beyond the parallel direction updates and boosts the candidate points' variety (Ter Braak & Vrugt 2008). Although many techniques exist to analyze the uncertainty (Makowski et al. 2002; Mantovan & Todini 2006; Yang et al. 2008; Zhang et al. 2015; Samadi et al. 2017), only a small number of studies have used the DREAM-ZS algorithm to analyze and calibrate the uncertainty in different hydrological models (Laloy & Vrugt 2012; Joseph & Guillaume 2013; Surfleet & Tullos 2013; Laloy et al. 2015; Leta et al. 2015; Zahmatkesh et al. 2015; Nourali et al. 2016; Engeland et al. 2016; Liu et al. 2017; Athira et al. 2018; Han & Zheng 2018).
According to the literature review, the assessment of uncertainty analysis on the SWAT model with DREAM-ZS algorithm for predicting runoff with CMIP5 models in multisite calibration has not been reported yet. Therefore, to fill this gap, the current study was conducted to survey the uncertainty prediction capabilities of the DREAM-ZS algorithm to extract posterior ranges of parameters to predict runoff in future time periods.
In this research, the DREAM-ZS algorithm is applied as a formal Bayesian method for uncertainty analysis, calibration, and measurement of SWAT model parameters in the KRB. First, the main goal of this research is to adjust and verify the SWAT model for the KRB as a large mountainous watershed in Iran with a new uncertainty analysis algorithm (DREAM-ZS). Second, the quantity of monthly runoff is calculated at the sub-basin level considering the effect of water resources management activities in the KRB in future times until the end of the 21st century and with new climatic models (GCMs models) of CMIP5 and RCP scenarios.
MATERIALS AND METHODS
The study area
The study area, Kashafrood Basin (KB), is located in Razavi Khorasan Region (RKR), in northeastern Iran. KB is a massive watershed in RKR with an area of 16,870 km2 and is known as a large-scale watershed. The longest river in this area is Kashafrood River (KR) and its total length is about 285 km. Geographically, KB is located approximately between longitudes 58° 15′ E and 61° 13′ E and latitudes 35° 35′ N and 37° 07′ N (Figure 1). The topographical elevation of the basin varies from 390 to 3,302 m above sea level. The average slope is 4.7% and the average elevation is 1,846 m above sea level. Figure 1 shows the topography of KB, and the location of gauging and meteorological stations. Due to mountainous topographical conditions, KB has a cold climate, semi-dry with high evapotranspiration in summer and low annual precipitation. The annual precipitation in this area is about 340 mm with relatively high spatial fluctuations. In addition, the maximum and minimum temperatures in a year are 20.6 and 7.1 °C, respectively (Afshar et al. 2018). The data set used in the study consists of a digital elevation model (DEM), land-use map (see Figure 2) and soil map from IRS satellite imagery, daily temperature and precipitation data and five monthly runoff data of the basin for 1997–2011 (see Figure 1). Torogh and Kardeh dams are located respectively south and north of Mashhad and provide water for agricultural and drinking sectors (Figure 1). The main soil type in the center and north of the basin is a heterogeneous blend of silt, sand, and clay, and in the southern and eastern parts of the basin it is clay and silt (Afshar et al. 2018).
DEM (digital elevation model) map and the location of gauging and meteorological stations in KB (adopted from Afshar et al. 2018).
DEM (digital elevation model) map and the location of gauging and meteorological stations in KB (adopted from Afshar et al. 2018).
SWAT model and data sets
SWAT, a semi-distributed conceptual model, was improved in United States Agriculture Department–Agricultural Research Services to set a rate and predict the effect of management design scenarios on groundwater resources, water quality, pollution loading, and soil erosion (Arnold et al. 1998; Narsimlu et al. 2015; Afshar & Hassanzadeh 2017; Afshar et al. 2018). The SWAT model was used at a daily time step for simulating the runoff in this study. To delineate the basin, the DEM of the watershed was used in the SWAT model. Each sub-basin is also divided into unique land-use/slope/soil hydrologic response units (HRUs). A description of each HRU is prepared by superimposing the land use, slope, and soil maps (Neitsch et al. 2011). In this research, KRB was divided into 217 sub-basins and 635 HRUs with 20, 10, and 20% threshold values for land-use, slope, and soil classes, respectively. The hydrometric stations are located at the outlet of the sub-basins and are selected in terms of area and topographic status to obtain a homogeneity between the sub-basins. By examining the overlapping of the maps, we concluded that by defining these thresholds, we can integrate small polygons with an area of less than 25 hectares, depending on the scale of the maps, into their larger adjacent polygons. Afterwards, a multiple slope option with four classes of 0–5, 5–10, 10–15 and >15% was selected for slope discretization (Afshar et al. 2018). Moreover, the modified Soil Conservation Service Curve Number (SCS-CN) was used to evaluate the surface runoff. The potential evapotranspiration (PET) was also calculated by the Hargreaves method. The Muskingum method was employed to route runoff. Flow was simulated from 1997 to 2011, and the first three years were selected as the model warm-up period. Finally, the monthly data set was used for calibration purposes, which related to an 11-year time span (2001–2011).
DREAM-ZS algorithm
In this study, a prior set of DREAM-ZS algorithm has been considered via the calibrated parameter set of calibration and uncertainty procedures (SWAT-CUP) and sequential uncertainty fitting program phase 2 (SUFI-2) process (Afshar et al. 2018). The flowchart of SWAT modeling with the DREAM-ZS sampler used for predictive uncertainty and calculation of parameter values is illustrated in Figure 3.
Sensitivity and uncertainty assessment
To reach a better understanding about the effect of sensitive parameters on runoff, the sensitivity analysis was done by SWAT-CUP tool. The dominant variables of SWAT (see Table 1) could be determined via Latin hypercube sampling (LH-OAT) sensitivity analyzing method (Van Griensven et al. 2006). LH-OAT was applied to explain the initial population of parameters. Posterior distribution of parameters (which were applied to generate model outputs) with uniform prior distributions was obtained based on the preselected ranges (Table 1). Therefore, 95% estimation uncertainty (95PPU) limits (including parameters, model framework, etc.) were achieved via computing 97.5 and 2.5% of cumulative distribution with regard to the outputs.
Sensitive parameters (adopted from Afshar et al. 2018)
Parameters . | Range of parameter . | Parameters . | Range of parameter . | ||
---|---|---|---|---|---|
Minimum . | Maximum . | Minimum . | Maximum . | ||
CN2 | −0.4 | 0.4 | ESCO | 0.01 | 1 |
GW_DELAY | 0 | 400 | SLSUBBSN | 10 | 150 |
ALPHA_BF | 0 | 1 | CH_N2 | 0 | 0.3 |
SOL_AWC | −0.3 | 0.3 | CH_K2 | 0 | 150 |
SOL_K | −0.8 | 0.8 | SFTMP | −5 | 5 |
SOL_BD | −0.3 | 0.3 | SMTMP | −5 | 5 |
GW_REVAP | 0.02 | 0.2 | SMFMN | 0 | 10 |
SHALLST | 0 | 1,000 | TIMP | 0.01 | 1 |
RCHRG_DP | 0 | 1 | SURLAG | 1 | 24 |
EPCO | 0.01 | 1 | PCPMM | −0.5 | 0.5 |
Parameters . | Range of parameter . | Parameters . | Range of parameter . | ||
---|---|---|---|---|---|
Minimum . | Maximum . | Minimum . | Maximum . | ||
CN2 | −0.4 | 0.4 | ESCO | 0.01 | 1 |
GW_DELAY | 0 | 400 | SLSUBBSN | 10 | 150 |
ALPHA_BF | 0 | 1 | CH_N2 | 0 | 0.3 |
SOL_AWC | −0.3 | 0.3 | CH_K2 | 0 | 150 |
SOL_K | −0.8 | 0.8 | SFTMP | −5 | 5 |
SOL_BD | −0.3 | 0.3 | SMTMP | −5 | 5 |
GW_REVAP | 0.02 | 0.2 | SMFMN | 0 | 10 |
SHALLST | 0 | 1,000 | TIMP | 0.01 | 1 |
RCHRG_DP | 0 | 1 | SURLAG | 1 | 24 |
EPCO | 0.01 | 1 | PCPMM | −0.5 | 0.5 |
DREAM-ZS algorithm is evaluated by three indices: P-factor shows the percentage of the observed data bracketed by 95PPU (Vrugt et al. 2009b; Abbaspour 2011; Alazzy et al. 2015), d-factor is the average thickness of 95PPU band divided by standard deviation of the related calculated variable (Abbaspour 2011), and Nash–Sutcliffe efficiency (NSE) represents the goodness-of-fit between the observed and simulated data (Nash & Sutcliffe 1970). P-factor, NSE, and d-factor are in the range of 0–100%, −∞–1, and 0–∞, respectively. In the calibration process, P-factor close to 100% (observations bracket via prediction uncertainty width) and d-factor value of zero (to reach a smaller uncertainty limit) show good agreement between the observed and simulated values (Yang et al. 2008). Also, the model prediction evaluation is facilitated when the NSE values are 0.75–1 and it may be satisfactory if NSE is greater than 0.36 (Moriasi et al. 2007). The parameter ranges can be determined at acceptable values of NSE efficiency, d-factor, and P-factor.
CMIP5 models
By CMIP5 outputs of GCMs models, 61 climatic time periods can be simulated all over the world by climate change models. In this research, about 14 models were selected to be used for selecting the appropriate CMIP5 models which are most adapted to the observational data (Chaturvedi et al. 2012; Xu et al. 2012; Jiang & Tian 2013; Ou et al. 2013; Iqbal & Zahid 2014; Salzmann et al. 2014; Woldemeskel et al. 2016; Afshar et al. 2017). Table 2 summarizes the total information and names of associated institutions of the models. CMIP5, unlike CMIP3, uses the new RCPs based on radiative forcing level by 2,100, shown in Table 3. The four RCP scenarios are RCP2.6, RCP8.5, RCP4.5, and RCP6.0 (Van Vuuren et al. 2011a). RCP2.6 is representative of the low end of the scenario literature in terms of emissions and radiative forcing (Van Vuuren et al. 2011b), while RCP8.5 includes high emission scenarios which cover the entire range of radiative forcing (Rao & Riahi 2006; Riahi et al. 2007, 2011). However, RCP2.6 and RCP8.5 were utilized in this research to investigate the future climate impacts in runoff. The bias corrected spatial downscaling method was chosen to downscale the precipitation data (for more information, see Afshar et al. 2017).
GCMs (global climate models) used in this research
Model name . | Organization . |
---|---|
HadGEM2-ES | MOHC, Met Office Hadley Centre, UK |
GFDL-CM3 | NOAA, Geophysical Fluid Dynamics Laboratory, USA |
GFDL-ESM2M | |
GFDL-ESM2G | |
MIROC5 | MIROC, Centre for Climate System Research, Japan |
MIROC-ESM | |
MIROC-ESM-CHEM | |
IPSL-CM5A-LR | IPSL, Institute Pierre Simon Laplace, France |
IPSL-CM5A-MR | |
NorESM1-M | NCC, Norwegian Climate Centre, Norway |
BCC-CSM1.1 | BCC, Beijing Climate Centre, China |
CSIRO-Mk3.6.0 | CSIRO-QCCCE, Commonwealth Scientific and Industrial Research Organization, Australia |
CCSM4 | NCAR, National Center for Atmospheric Research, USA |
CESM1(CAM5) | NSF-DOE-NCAR, Community Earth System Model Contributors |
Model name . | Organization . |
---|---|
HadGEM2-ES | MOHC, Met Office Hadley Centre, UK |
GFDL-CM3 | NOAA, Geophysical Fluid Dynamics Laboratory, USA |
GFDL-ESM2M | |
GFDL-ESM2G | |
MIROC5 | MIROC, Centre for Climate System Research, Japan |
MIROC-ESM | |
MIROC-ESM-CHEM | |
IPSL-CM5A-LR | IPSL, Institute Pierre Simon Laplace, France |
IPSL-CM5A-MR | |
NorESM1-M | NCC, Norwegian Climate Centre, Norway |
BCC-CSM1.1 | BCC, Beijing Climate Centre, China |
CSIRO-Mk3.6.0 | CSIRO-QCCCE, Commonwealth Scientific and Industrial Research Organization, Australia |
CCSM4 | NCAR, National Center for Atmospheric Research, USA |
CESM1(CAM5) | NSF-DOE-NCAR, Community Earth System Model Contributors |
Information about representative concentration pathways (RCPs) (Van Vuuren et al. 2011a)
RCPs . | Descriptions . |
---|---|
2.6 | Peak in radiative forcing at ∼3:1 Wm2 (∼490 ppm CO2 equivalent) before 2,100 and then decline (the selected pathway declines to 2.6 Wm2 by 2,100) |
4.5 | Stabilization without overshoot pathway to 4.5 Wm2 (∼650 ppm CO2 equivalent) at stabilization after 2,100 |
6.0 | Stabilization without overshoot pathway to 6 Wm2 at (∼850 ppm CO2 equivalent) stabilization after 2,100 |
8.5 | Rising radiative forcing pathway leading to 8.5 Wm2 (∼1,370 ppm CO2 equivalent) by 2,100 |
RCPs . | Descriptions . |
---|---|
2.6 | Peak in radiative forcing at ∼3:1 Wm2 (∼490 ppm CO2 equivalent) before 2,100 and then decline (the selected pathway declines to 2.6 Wm2 by 2,100) |
4.5 | Stabilization without overshoot pathway to 4.5 Wm2 (∼650 ppm CO2 equivalent) at stabilization after 2,100 |
6.0 | Stabilization without overshoot pathway to 6 Wm2 at (∼850 ppm CO2 equivalent) stabilization after 2,100 |
8.5 | Rising radiative forcing pathway leading to 8.5 Wm2 (∼1,370 ppm CO2 equivalent) by 2,100 |
The best CMIP5 models with the highest fit to the observed data were selected on the following steps, after downscaling the climatic parameters of the precipitation and temperature: The historical data of four stations which were located on the network nodes (with intervals of 0.5 degrees (longitude and latitude) to each other) around each station were extracted by ArcGIS. Then, the historical data of the station were calculated via inverse distance weight (IDW) method. Thereafter, according to the precipitation data and Thiessen's polygon averaging method, the CMIP5 model was assessed by Nash–Sutcliffe efficiency (NSE), percent of bias (PBIAS), and ratio of the root mean square error to the standard deviation of measured (RSR) criteria and coefficient of determination (R2) (Moriasi et al. 2007). The results showed that GFDL-ESM2G and MIROC-ESM might be chosen as the best climate models for comparing and predicting the effects of future climate change on the runoff component in the KRB. These models showed the highest compliance with observational precipitation data, based on evaluation metrics such as NSE (MIROC-ESM) = 0.95 and NSE (GFDL-ESM2G) = 0.92), PBIAS (MIROC-ESM) = −2.88 and PBIAS (GFDL-ESM2G) = −2.93), R2 (MIROC-ESM) = 0.97 and R2 (GFDL-ESM2G) = 0.94), and RSR (MIROC-ESM) = 0.33 and RSR (GFDL-ESM2G) = 0.37).
RESULTS AND DISCUSSION
Since the observed data are very restricted, only the results related to the calibration data set are illustrated and the results of the validation data set are excluded.
Results of DREAM-ZS based optimization
DREAM-ZS was employed for separating the effects of data, model structure, and parameter set uncertainties from the total predictive uncertainty (Vrugt et al. 2008). To choose candidate points for the individual chains, this classic method uses samples of previous steps. Three parallel chain runs were assigned to the DREAM-ZS algorithm (to compute the posterior distribution) initializing prior parameter set distribution samples for 50,000 runs (Table 4). DREAM-ZS convergence was constantly checked by R statistic for each parameter (see Figure 4) and could be achieved when all the sampling parameters fell under the threshold of 1.2 (see Figure 4, dashed line). In the current research, the algorithm reached the favorite convergence rate after performing a few iterations (about 5,250 runs) for all sensitive parameters.
Fitted value of posterior parameters with ranges of parameters
Parameter . | Range of parameter . | Fitted value . | Parameter . | Range of parameter . | Fitted value . | ||
---|---|---|---|---|---|---|---|
Minimum . | Maximum . | Minimum . | Maximum . | ||||
CN2 | 0.197 | 0.281 | 0.23 | ESCO | 0.41 | 0.501 | 0.50 |
GW_DELAY | 132.7 | 159 | 133.55 | SLSUBBSN | 96.5 | 109.5 | 99.54 |
ALPHA_BF | 0.0 | 0.08 | 0.042 | CH_N2 | 0.102 | 0.13 | 0.12 |
SOL_AWC | 0.21 | 0.28 | 0.28 | CH_K2 | 101.5 | 114 | 101.72 |
SOL_K | 0.15 | 0.32 | −0.16 | SFTMP | −1.0 | −1.81 | −1.45 |
SOL_BD | 0.095 | 0.167 | −0.11 | SMTMP | −0.95 | −1.92 | −1.70 |
GW_REVAP | 0.085 | 0.115 | 0.11 | SMFMN | 7.5 | 8.95 | 8.16 |
SHALLST | 460.8 | 510.5 | 485.30 | TIMP | 0.50 | 0.70 | 0.51 |
RCHRG_DP | 0.229 | 0.32 | 0.24 | SURLAG | 9.60 | 12.2 | 10.22 |
EPCO | 0.29 | 0.41 | 0.29 | PCPMM | −0.177 | −0.325 | −0.25 |
Parameter . | Range of parameter . | Fitted value . | Parameter . | Range of parameter . | Fitted value . | ||
---|---|---|---|---|---|---|---|
Minimum . | Maximum . | Minimum . | Maximum . | ||||
CN2 | 0.197 | 0.281 | 0.23 | ESCO | 0.41 | 0.501 | 0.50 |
GW_DELAY | 132.7 | 159 | 133.55 | SLSUBBSN | 96.5 | 109.5 | 99.54 |
ALPHA_BF | 0.0 | 0.08 | 0.042 | CH_N2 | 0.102 | 0.13 | 0.12 |
SOL_AWC | 0.21 | 0.28 | 0.28 | CH_K2 | 101.5 | 114 | 101.72 |
SOL_K | 0.15 | 0.32 | −0.16 | SFTMP | −1.0 | −1.81 | −1.45 |
SOL_BD | 0.095 | 0.167 | −0.11 | SMTMP | −0.95 | −1.92 | −1.70 |
GW_REVAP | 0.085 | 0.115 | 0.11 | SMFMN | 7.5 | 8.95 | 8.16 |
SHALLST | 460.8 | 510.5 | 485.30 | TIMP | 0.50 | 0.70 | 0.51 |
RCHRG_DP | 0.229 | 0.32 | 0.24 | SURLAG | 9.60 | 12.2 | 10.22 |
EPCO | 0.29 | 0.41 | 0.29 | PCPMM | −0.177 | −0.325 | −0.25 |
Convergence rate for Gelman–Rubin scale-reduction metric of 20 parameters of SWAT according to the DREAM-ZS method.
Convergence rate for Gelman–Rubin scale-reduction metric of 20 parameters of SWAT according to the DREAM-ZS method.
After these chains converged, the posterior parameter distributions of the SWAT model were produced with the last 20% of the parameters' set (i.e., 10,000 parameters' set). These plots are illustrated by Figure 5. Also, the best values (maximum likelihood) for all parameters are represented by (×) (see Table 5) where the x-axis represents the posterior ranges for all parameters. Results show that some parameters are efficiently optimized by the DREAM-ZS process for a complex model and resemble less uncertainty as a simply identifiable parameter. Furthermore, EPCO, TIMP, and ALPHA_BF values maintain a uniform distribution of shape and show the low uncertainty in the SWAT model. In other parameters, distributions of CH_N2, GW_DELAY, and ESCO are markedly varied to be concentrated around their lower or upper limits. The ranges of other parameters based on the results must be relaxed as well. Some parameters reached similar distributions or very close distributions. For instance, SOL_AWC, SOL_BD, ALPHA_BF, SURLAG, SLSUBBSN, and CN2 had almost similar values of 0.26, −0.14, 0.04, 11.2, 99, and 0.21, respectively.
Values of criteria of total uncertainty (U1) and parameter uncertainty (U2)
Hydrometry gauges . | U1 . | U2 . | ||||
---|---|---|---|---|---|---|
NSE . | P . | d . | NSE . | P . | d . | |
SARASSHA | 0.72 | 0.91 | 1.96 | 0.72 | 0.13 | 0.19 |
ZIRBAGOL | 0.70 | 0.89 | 1.48 | 0.70 | 0.14 | 0.22 |
GOLJAGHR | 0.57 | 0.98 | 5.00 | 0.57 | 0.14 | 0.15 |
HESDEHB | 0.54 | 0.98 | 3.43 | 0.54 | 0.15 | 0.16 |
KARTIAN | 0.56 | 0.93 | 2.79 | 0.56 | 0.11 | 0.17 |
Hydrometry gauges . | U1 . | U2 . | ||||
---|---|---|---|---|---|---|
NSE . | P . | d . | NSE . | P . | d . | |
SARASSHA | 0.72 | 0.91 | 1.96 | 0.72 | 0.13 | 0.19 |
ZIRBAGOL | 0.70 | 0.89 | 1.48 | 0.70 | 0.14 | 0.22 |
GOLJAGHR | 0.57 | 0.98 | 5.00 | 0.57 | 0.14 | 0.15 |
HESDEHB | 0.54 | 0.98 | 3.43 | 0.54 | 0.15 | 0.16 |
KARTIAN | 0.56 | 0.93 | 2.79 | 0.56 | 0.11 | 0.17 |
To evaluate the uncertainty of overland flow parameters for five considered hydrometry gauges following the convergence, the last 20% of the SWAT samples which fitted the model well were extracted to produce the output values and 95% confidence interval was defined by computing 97.5 and 2.5% levels of the cumulative distribution of the output variables. The DREAM-ZS algorithm plotted the parameter uncertainty (dark shaded area), the time series data related to total uncertainty (light shaded area), the best simulation obtained by SLS likelihood function (black line), and the monthly flow discharge (filled circles) observed in five gauges from 2001 to 2011, as shown in Figure 6. The model for KRB was calibrated according to the data from one downstream (KARTIAN), one upstream (SARASSHA), and three mid-stream (GOLJAGHR, HESDEHB, and ZIRBAGOL) gauging stations. Also, in Figure 5, the x-axis of plots is demonstrated time from 0 to 132 months. The obtained values for three evaluation measures (NSE, d-factor, and P-factor), which show the performance of uncertainty analysis, are presented in Table 5. The uncertainty origins, such as U1 and U2, might be considered as the total (and parameter) uncertainty.
As can be seen in Figure 6, the most observation points (over 90%) of DREAM-ZS take place within 95% predictive uncertainty limits and P-factor demonstrates high performance of the DREAM-ZS algorithm. Also, the results illustrate that the total uncertainty (and parameters' set uncertainty) band covers 89–98% (and 11–34%) of the measurements' data by identical d-factor of 1.36–5 (and from 0.15 to 0.29), respectively. The minimum and maximum P and d-factors are for Gand ZIRBAGOL (1.48 and 89%) and GOLJAGHR (5 and 98%) stations, respectively. Considering only the SWAT parameter uncertainty, 11–15% of the observations are bracketed. On the contrary, limits of the lower and upper parameters' uncertainty involve most of the observation of the band. Also, the gap between observed data and the parameter uncertainty may result from the model structural, input and output insufficiency (Laloy et al. 2010). Although the total prediction uncertainty bounds are quite wide, the parameter uncertainty bounds are quite narrow. On average for all stations, the value of NSE obtained was around 0.62 (Table 5). Also, the NSE values show relatively good agreement between simulated and observed runoff values. Using the DREAM-ZS algorithm, the peak runoff values are mostly underestimated, while in only a few months are they overestimated (especially in the end part of the calibration period) in all runoff stations. The accuracy of the model is good in peak discharge estimations, although not excellent. In wet periods, the rainfall gauge stations are not a very good representative of precipitation over the basin, because the map shows there are no rain gauges in some high and most low altitude parts of the study basin. Also, the KRB is a mountainous watershed and is exposed to high spatial and temporal variability in rainfall distribution that cannot be captured exactly using rain gauges in the basin. Therefore, such limitation could have contributed to some degree of simulation uncertainty (Cho et al. 2009; Memarian et al. 2014). Overall, the results show a good agreement between the simulated and measured data providing a realistic modeling.
Results of runoff predictions for the future via CMIP5 models
The spatial variations of runoff in the KRB during the historical (1997–2011) and future time periods (near: 2014–2042, intermediate: 2043–2071, and far: 2072–2100) (Afshar et al. 2018) in GFDL-ESM2G and MIROC-ESM models under both scenarios (RCP8.5 and RCP2.6) are respectively illustrated in Figures (7)–(9).
Spatial variation of runoff in the KRB for both models and under RCP2.6 in future periods (left panels belong to MIROC-ESM, right panels belong to GFDL-ESM2G).
Spatial variation of runoff in the KRB for both models and under RCP2.6 in future periods (left panels belong to MIROC-ESM, right panels belong to GFDL-ESM2G).
Spatial variation of runoff in the KRB for both models and under RCP8.5 in future periods (left panels belong to MIROC-ESM, right panels belong to GFDL-ESM2G).
Spatial variation of runoff in the KRB for both models and under RCP8.5 in future periods (left panels belong to MIROC-ESM, right panels belong to GFDL-ESM2G).
The simulation results show that in the northern part of the watershed, the runoff values would be higher than the southern part with almost the same pattern for different time periods. However, there are some differences in the sub-basin scale in the amount of runoff (Figures (7)–(9)). The mean annual runoff in the near future under RCP2.6 (RCP8.5) in GFDL-ESM2G and MIROC-ESM models, respectively, is about 173 (182) and 108 (147) mm/year. Also, in the intermediate future period, the mean value of annual runoff under RCP2.6 and RCP8.5 in GFDL-ESM2G (MIROC-ESM) is 141 (70) and 160 (119) mm/year, respectively. Simulated annual runoff in the far future period via the GFDL-ESM2G (MIROC-ESM) model under RCP2.6 and RCP8.5 is about 176 (124) and 125 (75) mm/year. The results show that in the near and intermediate periods, by both models under RCP8.5, simulated runoff will be higher than RCP2.6, while in the far future under RCP2.6, the mean annual of the simulated runoff will be higher than that in the other scenario. Furthermore, the average annual runoff obtained by the GFDL-ESM2G model in all future periods and under both scenarios is higher than the other selected model (MIROC-ESM). By comparing the runoff values in historical and future periods it is clear that the hydrology condition of the watershed will be improved due to the increase of runoff in all future periods. More details are presented in Table 6 with regards to comparing runoff in future time periods and via both CMIP5 models with historical time.
Comparison of runoff in historical and future time periods in both CMIP5 models (mm/year)
Period times . | GFDL-ESM2G . | MIROC-ESM . | ||
---|---|---|---|---|
RCP2.6 . | RCP8.5 . | RCP2.6 . | RCP8.5 . | |
Near future | 172.72 | 181.53 | 108.19 | 147.50 |
Mid future | 140.65 | 160.42 | 70.24 | 118.71 |
Late future | 176.51 | 124.74 | 123.86 | 75.44 |
Historical | 37.54 |
Period times . | GFDL-ESM2G . | MIROC-ESM . | ||
---|---|---|---|---|
RCP2.6 . | RCP8.5 . | RCP2.6 . | RCP8.5 . | |
Near future | 172.72 | 181.53 | 108.19 | 147.50 |
Mid future | 140.65 | 160.42 | 70.24 | 118.71 |
Late future | 176.51 | 124.74 | 123.86 | 75.44 |
Historical | 37.54 |
The minimum, maximum, and mean annual runoff values are shown in Tables 7 and 8, respectively, for upstream (in the higher part of a stream) and downstream (the direction that a river flows) sub-basins in the KRB. In general, sub-basins located in the upstream show higher runoff value, while downstream sub-basins show smaller annual runoff values under both scenarios and models. The average annual runoff values for upstream in the MIROC-ESM model under RCP2.6 for near, intermediate, and far future are about 138.8, 94.8, and 163.6 mm/year, while in the downstream it decreases to around 43.5, 51.5, and 47.9%, respectively. In the mentioned model under RCP8.5 scenario almost the same pattern is seen in the future periods and in downstream runoff the three future periods show 41.4, 38.2, and 45.3% less values than those for the upstream. In the MIROC-ESM model under both scenarios and future periods, the range of mean annual runoff variation in upstream and downstream is between 38.2 and 51.1%, while in the GFDL-ESM2G model the range of this variation is between 14.0 and 44.0%. The lowest difference between upstream and downstream will happen in the far future under RCP8.5 in the GFDL-ESM2G model and under RCP2.6 (RCP8.5), the variation in the near, intermediate, and far future is about 40.9 (44.0), 41.0 (30.2), and 38.7 (14.0)%, respectively. On the other hand, runoff declined about 60.7% from the upstream to downstream in the historical time period.
Variations of runoff for upstream and downstream of the KRB in the MIROC-ESM model
Range . | Historical . | RCP2.6 . | RCP8.5 . | ||||
---|---|---|---|---|---|---|---|
Near future . | Mid future . | Late future . | Near future . | Mid future . | Late future . | ||
Upstream sub-basins of KRB | |||||||
Min | 1.57 | 46.39 | 23.87 | 48.84 | 73.13 | 38.08 | 23.21 |
Max | 125.95 | 249.20 | 179.78 | 278.57 | 353.00 | 467.46 | 206.22 |
Mean | 54.22 | 138.84 | 94.84 | 163.58 | 186.68 | 147.26 | 97.91 |
Downstream sub-basins of KRB | |||||||
Min | 0.70 | 6.63 | 4.78 | 9.54 | 15.48 | 10.11 | 7.46 |
Max | 96.23 | 215.18 | 155.08 | 254.13 | 277.05 | 465.54 | 203.75 |
Mean | 21.32 | 78.38 | 46.32 | 85.23 | 109.40 | 90.94 | 53.59 |
Range . | Historical . | RCP2.6 . | RCP8.5 . | ||||
---|---|---|---|---|---|---|---|
Near future . | Mid future . | Late future . | Near future . | Mid future . | Late future . | ||
Upstream sub-basins of KRB | |||||||
Min | 1.57 | 46.39 | 23.87 | 48.84 | 73.13 | 38.08 | 23.21 |
Max | 125.95 | 249.20 | 179.78 | 278.57 | 353.00 | 467.46 | 206.22 |
Mean | 54.22 | 138.84 | 94.84 | 163.58 | 186.68 | 147.26 | 97.91 |
Downstream sub-basins of KRB | |||||||
Min | 0.70 | 6.63 | 4.78 | 9.54 | 15.48 | 10.11 | 7.46 |
Max | 96.23 | 215.18 | 155.08 | 254.13 | 277.05 | 465.54 | 203.75 |
Mean | 21.32 | 78.38 | 46.32 | 85.23 | 109.40 | 90.94 | 53.59 |
Variations of runoff for upstream and downstream of the KRB in the GFDL-ESM2G model
Range . | Historical . | RCP2.6 . | RCP8.5 . | ||||
---|---|---|---|---|---|---|---|
Near future . | Mid future . | Late future . | Near future . | Mid future . | Late future . | ||
Upstream sub-basins of KRB | |||||||
Min | 1.57 | 74.61 | 60.56 | 83.19 | 74.67 | 71.52 | 49.78 |
Max | 125.95 | 397.13 | 330.62 | 403.39 | 434.39 | 335.27 | 213.59 |
Mean | 54.22 | 217.99 | 177.61 | 219.55 | 233.74 | 189.47 | 134.29 |
Downstream sub-basins of KRB | |||||||
Min | 0.70 | 47.94 | 19.94 | 53.87 | 36.12 | 41.31 | 27.96 |
Max | 96.23 | 281.34 | 236.77 | 290.47 | 276.64 | 268.97 | 239.37 |
Mean | 21.32 | 128.68 | 104.70 | 134.64 | 130.75 | 132.16 | 115.46 |
Range . | Historical . | RCP2.6 . | RCP8.5 . | ||||
---|---|---|---|---|---|---|---|
Near future . | Mid future . | Late future . | Near future . | Mid future . | Late future . | ||
Upstream sub-basins of KRB | |||||||
Min | 1.57 | 74.61 | 60.56 | 83.19 | 74.67 | 71.52 | 49.78 |
Max | 125.95 | 397.13 | 330.62 | 403.39 | 434.39 | 335.27 | 213.59 |
Mean | 54.22 | 217.99 | 177.61 | 219.55 | 233.74 | 189.47 | 134.29 |
Downstream sub-basins of KRB | |||||||
Min | 0.70 | 47.94 | 19.94 | 53.87 | 36.12 | 41.31 | 27.96 |
Max | 96.23 | 281.34 | 236.77 | 290.47 | 276.64 | 268.97 | 239.37 |
Mean | 21.32 | 128.68 | 104.70 | 134.64 | 130.75 | 132.16 | 115.46 |
On a sub-basin scale (Figure 10), the highest rate of runoff in both models, RCPs, and in all future periods occurred in sub-basin number 18 and also under RCP8.5 respectively happened in sub-basin number 100 (MIROC-ESM in the intermediate and far future) and sub-basin 170 (GFDL-ESM2G model in far future). On the other hand, in the GFDL-ESM2G model under both scenarios and all future periods, the lowest value of runoff was seen in sub-basin number 190, while in the MIROC-ESM model in all future periods and RCPs sub-basin number 203 (except in far future and RCP2.6 in sub-basin 216), experienced the minimum value of this component.
Sub-basins' number map of KRB and selected sub-basin for runoff component.
CONCLUSIONS
In this study, formal Bayesian frameworks (DREAM-ZS) were investigated to simulate runoff discharge over a semi-dry basin system (with an area of 16,870 km2) in Iran by the SWAT model. Twenty parameters out of initial parameters were determined as sensitive parameters to assess uncertainty and calibration processes at five gauging stations. The DREAM-ZS algorithm was used to optimize the parameters by taking SLS as a likelihood function. After 50,000 runs, the best posterior parameters were obtained. Investigation of marginal density implied that both parameter ranges and shape had almost a uniform posterior distribution for some parameters. The distribution of ALPHA_BF, CN2, and CH_N2 could be well-defined, whereas the remaining parameters were either positively or negatively skewed. Results of calibration and uncertainty analysis indicated that at all stations the peak runoff in most months predicts underestimation via the calibrated model, while the falling limb had overestimated values in some months. The base flow had a high correspondence between simulated and observation data for all calibration stations. For the total uncertainty, NSE was obtained in the range of 0.54–0.72, most data were bracketed within 95PPU bounds (89–98%), and various d-factors (1.48–5.00) were noticed. In DREAM-ZS, the peak runoff points were slightly overestimated and underestimated as compared to the observed data in some months. Therefore, uncertainty analysis approved the high efficiency of DREAM-ZS in predicting parameters uncertainty via the KRB as a future research plan. The mean annual runoff in the MIROC-ESM model under RCP2.6 (RCP8.5) from near to far future periods, respectively, was 108 (147), 70 (119), and 124 (75) mm/year while in the other selected model the results of the simulation were 173 (181), 141 (160), and 176 (125) mm/year. Comparison of GFDL-ESM2G and MIROC-ESM models for near, intermediate, and far future periods under RCP2.6 (RCP8.5) showed the mean annual runoff in GFDL-ESM2G was about 59.6% (23.0%), 100.2% (35.1%), and 42.5% (65.3%) more than the other model. The results of this research can provide reliable information about runoff at sub-basin scale for policymakers.
A follow-up work to this study would be an investigation on water resource security (demand, availability, scarcity, reliability, and vulnerability) based on the blue and green water concepts by applying the mentioned models in this region.