Impact assessment of climate change on hydro-climatic conditions of arid and semi-arid watersheds (case study: Zoshk-Abardeh watershed, Iran)

Climatic paleontology evidences that climate change has always been present throughout the history of the planet, but the climatic changes of the last century have two distinct features, as compared with past climatic changes. First, human activities play a greater role in the nature of the current climate change. Second, the speed of recent climatic change is greater so that a considerable number of changes will be occurring in the Earth's atmosphere over the short term (IPCC 2007). The average surface temperature of the Earth is increasing, which is mainly the result of greenhouse gas emissions (Sartori 2019). The recent scenarios of the Intergovernmental Panel on Climate Change (IPCC) show a global increase in temperature, i.e. 1.1 °C in the current century. In addition to the change in the average values of climate variables, the change in rainfall in relation to its past values will also be other outcomes, which means increasing the intensity of dry or wet conditions (IPCC 2007). Changes in the cloud and precipitation formation are the secondary aspects of climate change. These changes are not simply made due to an increase in greenhouse gas concentrations, and they occur in response to global warming (changes in the concentration of aerosols directly affect clouds and precipitation; however, this effect is most severe on temperature). Thus, in many areas, the greatest impact of climate change on hydrology is likely due not only to the relatively small changes in the turbulent behavior of clouds and precipitation but also to the direct effect of temperature rise on the hydrology cycle. In this case, the amount of evapotranspiration from the surface of the soil and plants will necessarily increase significantly, while the size of the snow grains will become smaller and thaw earlier in the warm season (Fung et al. 2011). The increase in evaporation will undoubtedly drastically reduce the availability of water and provide the basis for the occurrence and intensification of drought effects. Runoff is a component of the hydrology cycle that is influenced by many parameters. Certainly, increasing the temperature and decreasing snowfall in the mountainous basins along with the evaporation of the water bodies will also lead to a decrease in the volume of runoff and water reserves and an intensification of the hydrological drought phenomenon.

Nowadays, global warming has significant effects on precipitation and runoff yield and water resources due to the increased concentration of greenhouse gases (Zhang et al. 2017). Rainfall as a key factor in changing the frequency and range of hydrological cycle has serious consequences in social, economic, and agricultural developments (Zhang et al. 2017). The average of climatic variables, especially the components of temperature, precipitation, and runoff in the annual or seasonal scale, plays a predominant role in the hydrological cycle. The climate change is usually investigated through the assessment of the average of these parameters, as well (Afshar et al. 2017). Abnormalities of temperature and precipitation have substantial impacts on the agricultural industry and the rapid development of the country. Population growth, urbanization, and rapid economic development impose massive challenges on the available water resources of Iran, now and especially in the future (Madani et al. 2016; Afshar et al. 2017).

At present, general circulation models (GCMs) are the most frequently used models for the projection of different climatic change scenarios (IPCC 2007; Afshar et al. 2017). The IPCC collects and reviews global climate models as the assessment reports (ARs) of climate change. Up to now, the IPCC has released five different versions of GCMs, including First Assessment Report (FAR) models, Second Assessment Report (SAR) models, Third Assessment Report (TAR) models, Fourth Assessment Report (AR4) models, and Fifth Assessment Report (AR5) models. The IPCC AR5 has introduced a new way in the development of the scenarios. These scenarios span the range of plausible radiative forcing scenarios and are called representative concentration pathways (RCPs), whereas AR4 uses scenarios from the IPCC Special Report on Emission Scenarios (SERS). The RCPs span a wider range of possibilities than the SRES marker scenarios used in the modeling for the IPCC AR3 and AR4. In contrast to SRES, some of the RCPs also include mitigation and adaptation policies (IPCC 2013).

A number of studies have been done to investigate the impact of climate change on the hydrological components of watersheds in Iran based on the IPCC's Fourth Assessment Report models (AR4) (Vaghefi et al. 2015; Eslamian et al. 2017). These models together with older emission scenarios (SRES) have less resolution, in comparison with the Fifth Assessment Report (AR5) models. Thus, climate change studies with higher resolution climate models under the new emission scenarios (RCPs) of the AR5 seem necessary in the watersheds of Iran. In a recent study by Afshar et al. (2017), the annual climatic components of the Kashafrood basin in Iran were assessed in the historical and future periods by using the AR5 models (Afshar et al. 2017). According to the results, the precipitation component showed a significant decreasing change trend, and the average temperature also indicated a remarkable increase with a confidence level of 90, 99, and 99.9%. On the basis of the obtained results, the average temperature of the basin will increase from 0.56 to 3.3 °C, and precipitation will decrease about 10.7% by the end of the 21st century.

The study watershed in this research (Zoshk-Abardeh) is one of the important sources of income for the regional villagers and has a high ecotourism potential in Khorasan Razavi Province, Iran. It is considered as an urban watershed and according to the historical evidence has a high flood potential, as well. Thus, this study aims to assess the impact of climate change on the most important hydroclimatological factors (rainfall, temperature, and runoff) affecting on the ecotourism potential of the watershed.

The Zoshk-Abardeh Watershed as a sub-watershed of the Kashafrood basin, with an area of 9,225.9 ha, is located in the west of Mashhad in Khorasan Razavi Province, Iran (Figure 1). The average slope of the basin is 52.46%, its average altitude is 2,235 m above sea level, and the length of the basin is 21.53 km. The average annual precipitation of the Zoshk-Abardeh watershed is 400 mm. The maximum monthly precipitation is 71 and 65 mm in April and May, and each of these months has about 18 and 16% of the total annual rainfall, respectively. June, August, and September are the driest months of the year. This basin has the average, minimum, and maximum annual temperatures of 8.7, 2.8, and 14.2 °C, respectively. In this study, daily observational data for the Zoshk Station were obtained from the Khorasan Razavi Meteorological Organization and the Khorasan Razavi Regional Water Authority (KRRWA). Land use, topography, and soil classification maps were obtained from the Watershed Management and Natural Resources Organization (WMNRO) and the National Cartographic Center (NCC). The Zoshk hydrometric and rain gauge station was established by the KRRWA in 1978, and it is already under the control of the KRRWA. This station is a flood alarming station, which is recently equipped with several instruments such as limnograph, data logger, teleferic bridge, and real-time data transmitters. A list of data sets used in this work and their sources are provided in Table 1.

In recent years, various models have been used to simulate water-related processes, as computer systems evolve and their simulations grow (Arnold et al. 1998; Besharat et al. 2015; Afshar & Hassanzadeh 2017; Dadfar et al. 2019). Among them, the Soil and Water Assessment Tool (SWAT), a semi-distributed conceptual model for watersheds, was developed by the United States Department of Agriculture–Agricultural Research Services (USDA–ARS) for assessing and forecasting the impact of different management scenarios on water quality, groundwater resources, soil erosion, and pollution loading (Arnold et al. 1998). Each sub-basin in the model is sub-divided into the Hydrological Response Units (HRUs), which are homogeneous units in terms of soil parameters, land use, and slope (Memarian et al. 2014; Afshar et al. 2017). In this study, runoff volume was estimated using the Soil Conservation Service (SCS) method. The Manning equation was utilized to calculate flow velocity. The Muskingum approach was used in the routing phase, as well. The growth parameters in the SWAT crop database were updated, and the required information records were also added to the SWAT soil database.

It should be noted that one of the limitations in simulating the impacts of climate change using the SWAT model is the fact that this model cannot account for all the changes caused by climate change in the new water cycle. Thus, we have to ignore the effects of climate change on some components of the water cycle in the simulation process.

The NS criterion is a common and valid coefficient in the calibration process of hydrologic models. Meanwhile, the NS is more reliable than other criteria and it can be applied to a different type of modeling context. This advantage indicates its flexibility as a goodness-of-fit statistic. On the other hand, NS efficiency can be used to quantitatively describe the accuracy of model outputs other than discharge. This indicator can be used to describe the predictive accuracy of other models as long as there is observed data to compare the model results to. For example, NS efficiency has been reported in the scientific literature for model simulations of discharge; water quality constituents such as sediment, nitrogen, and phosphorus loading. The evaluation of model prediction is very fine when the NS values are between 0.75 and 1, and it may be satisfactory if NS is greater than 0.36 (Moriasi et al. 2007). Moreover, the maximum value of R² is 1; the higher values of this coefficient indicate a better performance of the model.

The SUFI-2 algorithm for each of the input parameters of the SWAT model presents a t-stat value and then determines the critical parameters (Abbaspour et al. 2015). On the other hand, the significance of differences in multiple runs for the parameter b_i is examined using the t-test. A t-test determines the relative significance of each parameter b_i, and the t-statistic is a measure of sensitivity (larger absolute values indicate more sensitive parameters) (Abbaspour 2011). The t-stat absolute values of each parameter shows that the parameters with more t-stat have more relative sensitivity. The algorithm also provides a p-value for each parameter in order to determine the relative significant sensitivity of them, which is closer to the zero number indicating a greater sensitivity of the parameter (Abbaspour et al. 2015). For validation, the calibrated parameter ranges without any further changes were used, and an iteration was run with the same number of simulations as utilized for calibration. The efficiency measures were computed during the validation process, as well.

In contrast to emission scenarios used in the AR4 models of the IPCC, the RCP is utilized in the AR5. New emission scenarios are based on the emission forcing level until 2100 (IPCC 2013). Radiative forcing (Watts per square meter; W/m²) is the extra heat at the lower atmosphere that will retain as a result of additional greenhouse gases. In order to investigate the future climate change, the Model for Interdisciplinary Research on Climate-Earth System Models (MIROC-ESM), which has been cooperatively developed by the University of Tokyo, National Institute for Environmental Studies (NIES), and the Japan Agency for Marine-Earth Science and Technology (JAMSTEC) (K-1 Model Developers 2004; Nozawa et al. 2007), was selected among the newest extracted models presented in the AR5. According to Afshar et al. (2017), this model showed the highest agreement with observational data (Figures 2 and 3) based on the evaluation criteria, i.e. NS efficiency coefficient, percent of bias (PBIAS), coefficient of determination (R²), and ratio of the root-mean-square-error to the standard deviation (RSR) (Afshar et al. 2017).

This model consists of four emission forcing scenarios, namely RCP2.6, RCP4.5, RCP6.0, and RCP8.5 (Table 2). These RCPs were included one mitigation scenario leading to a very low forcing level (RCP2.6), two medium stabilization scenarios (RCP4.5 and RCP6.0), and one very high baseline emission scenarios (RCP8.5) (Van Vuuren et al. 2011). For example, in RCP2.6, which has the lowest emission, the total radiation forcing level will reach by a maximum of 3 W/m² in 2050, and then it follows a downward trend (Van Vuuren et al. 2007). The RCP4.5 and RCP6.0 scenarios would be a stable scenario with an increased total emission forcing level in 2070 and a constant concentration after 2070 (Clarke et al. 2007; Hijioka et al. 2008). Furthermore, the RCP8.5 scenario steadily increases emissions forcing by the end of the 21st century and is approximately equal to 8.5 W/m² (Riahi et al. 2007).

The general atmospheric circulation models provide a good opportunity to predict future climate under different radiative forcing scenarios. These models have a relatively low resolution and are not suitable for the studies at the regional and local level. Therefore, the output of these models should be downscaled to be used on a regional scale (Sigdel & Ma 2016). The downscaling approaches are classified into dynamic and statistical techniques. In this study, a statistical downscaling model, i.e. Bias Correction and Spatial Disaggregation (BCSD), was used as the spatial disaggregation and bias correction method (Ahmed et al. 2013). The spatial disaggregation method is performed in three steps. In the first step, the correction of statistical deviation of rainfall and temperature components (resulted from climate change simulation) in the scale of the grid of GCMs is calculated in a monthly time step. In the second step, the spatial downscaling is carried out from the network scale of GCMs to the target network (local or regional). Then, in the third stage, the outputs are disaggregated from a monthly time step to a daily time step (Afshar et al. 2017). The historical data of four stations in the vicinity of each of the reference stations (observational stations) were extracted using ArcGIS with the intervals of 0.5° (the format of the climatic data is Network Common Data Form (NetCDF)). The historical data on the reference station were obtained having the information of four stations in the vicinity and using interpolation methods such as bilinear interpolation and inverse distance weight.

The SWAT (Arnold et al. 1998) was employed for hydrological simulation of the watershed based on the downscaled outputs (using the BCSD method) of climate change models for historical and future periods. The future projection was conducted in the near- (2006–2037), mid- (2037–2070), and far- (2070–2100) future periods related to historical records in the period of 1992–2005. The SWAT model was executed using the three years as the warm up period with the definition of 294 HRUs and the delineation of 12 sub-basins. Figure 4 depicts the modeling framework of the current study.

The SWAT-CUP (Abbaspour et al. 2007) was employed to run the SUFI-2 algorithm for SWAT calibration and sensitivity analysis. Seventeen parameters were chosen for SWAT sensitivity analysis using the 500 simulations of SUFI-2. Table 3 shows the t-stat and p-value for various parameters affecting the runoff volume. Comparing the obtained t-stat for each parameter shows that the parameters CN2, SOL_BD, and SOL_K have the highest relative sensitivity and CH_K2, SMFMX, and SOL_AWC parameters have the lowest relative sensitivity. The high sensitivity of SWAT to the alterations in parameters CN2, SOL_BD, and SOL_K was also reported by many research works such as Arnold et al. (1998) and Memarian et al. (2014).

The calibration and validation steps were carried out in accordance with the data of the Zoshk hydrometric station on a monthly time step. The calibration was performed using water flow data in the period of 2000–2006, and the model was validated based on the water flow records during the period of 2007–2010. The calibration process was accomplished using the 3,000 simulations of the SUFI-2 algorithm on the sensitive parameters detected through the sensitivity analysis. Finally, the optimal values of these parameters were determined for the simulation of monthly runoff (Table 4).

The model performance was evaluated using R² and NS coefficients between observed and simulated records. Based on the results, the coefficients R² and NS were estimated to be 0.52 and 0.47 for runoff simulation in the calibration stage and 0.46 and 0.42 in the validation stage, respectively (Table 5). Given that the amounts of NS are between 0.36 and 0.75 (Moriasi et al. 2007), the results of the model are acceptable, but the model performance is not significantly high, in general. A comprehensive literature review shows that the lack of water diversion data for agricultural and human consumption (as well as wastewater discharge), lack of irrigation use, agricultural management and crop planting data within the watershed, lack of complete understanding of surface and groundwater interactions, and lack of the knowledge on construction projects have a large negative impact on the efficiency of hydrological models (Kouchi et al. 2017).

The hydrograph of monthly observed and simulated flow rates during calibration and validation periods, which is shown in Figure 5, has been used in order to evaluate the efficiency of the model in the base and peak discharge simulation and to investigate their time harmonization to actual data, as well. The analysis of this hydrograph shows that the model has estimated the amount of peak discharge less than the actual amounts, which is confirmed by the average monthly simulated discharge during calibration and validation periods. The average monthly simulated discharge during calibration and validation periods is 0.33 and 0.38 m³/s, respectively, while the observed data in these two periods are 0.45 and 0.50 m³/s, respectively. The degree of uncertainty was calculated by two factors called r-factor and P-factor. As shown in Figure 6, the results showed that more than 50% of the observational data in both calibration and validation phases are bracketed by the 95PPU uncertainty estimation band and the r-factor fluctuates around 0.5, which indicates a rather acceptable degree of certainty in simulation (Abbaspour 2011; Ozdemir et al. 2017). It should be noted that if the results are of high quality, then 80–100% of the data should be bracketed by the 95PPU (P-factor), while low quality results may contain many outliers, and it may be sufficient to account for only 50% of the data in the 95PPU (Abbaspour et al. 2007). Some reasons that cause low uncertainty can be rooted in the development and human activities (human activities including water withdrawal for aquaculture and gardening activities, wastewater entry into the river, and partial narrowing of the river bed due to the occupation of stream bank) in the catchment area.

The reason for the marginal performance of the model can be related to the following. (1) There is an error in the observational statistics of the hydrometric station and rain gauge, which is confirmed by the experts in the regional organization of water authority. (2) The amount of water extracted from the river by the orchards and upstream wells of the basin is not available. (3) There are also a number of valleys and contact springs in the ridge of the basin, with no accurate statistical data on the discharge of those springs in different seasons, and there is only one estimated discharge of them. The performance of the SWAT model and SUFI-2 algorithm have been evaluated in several watersheds in semi-arid regions of Iran. For example, Shafiei et al. (2013) evaluated their performance in hydrological simulation at one of the closest basins (Neyshabour Basin) to the study area (Shafiei et al. 2013). In this study, Hussein Abad and Andarab hydrometric stations were used in the calibration process. The results showed that the R² and NS coefficients at the calibration stage for the Hussein Abad station were 0.74 and 0.65, and for the Andarab station were 0.31 and 0.11, respectively. During the validation process, the R² and NS coefficients showed a significant change. Those were 0.42 and −0.82 for the Hussein Abad station, and 0.06 and −1.7 for the Andarab station, respectively. The P-factor and r-factor at the calibration phase were 67% and 1.41 for the Hussein Abad, and 23% and 0.25 for the Andarab station, respectively. During the validation process, the P-factor and r-factor for the Hussein Abad station were 53% and 2.8, and for the Andarab station were 3% and 0.38, respectively. Results confirmed that the model was not able to correctly simulate low flows. This problem was explained by the simplifications of the conceptual model and the complex interaction between runoff and subsurface flow in low rainfall events. The same approach was used to simulate the daily stream flow for the Roodan watershed (as a semi-arid basin) located in the southern part of Iran (Jajarmizadeh et al. 2012). For calibration, the values of NS and R² were reported at 0.66 and 0.68, respectively, and for validation, the values were 0.51 and 0.55. They reported that the performance of modeling was acceptable, but the SWAT model tended to underestimate the simulation of high stream flow for both calibration and validation periods. This issue was attributed to the application of the plant evapotranspiration method that imposes an overestimation of runoff for high flows and leads to better simulation of perennial flows. Moreover, they confirmed that the precipitation duration and intensity were not being considered by the SCS method for the simulation of stream flow in the SWAT model. This limitation is more profound for the watersheds that are located in arid and semi-arid area. The application of the SWAT model for estimating runoff in two mountainous semi-arid basins in central Iran was evaluated by Rostamian et al. (2008). Model calibration and uncertainty analysis were performed with the SUFI-2 algorithm. The measures P-factor and r-factor indicated a relatively acceptable model calibration and accounting of the uncertainties. The weakness of the model to simulate runoff for some months was probably due to the poor characterization of snowmelt processes in these mountainous watersheds, lack of sufficient discharge data, and lack of input data for the simulation of groundwater recharge and groundwater–river interaction. Based on the above literature review, it can be revealed that the SWAT model in many watersheds of arid and semi-arid regions of Iran due to different causes (most of them are reliable in this study and more than ever the lack of highly accurate data) cannot establish an enormous performance in a hydrological simulation of the basin. However, its acceptable robustness could be a suitable ground for its applications in other studies such as climate change.

The parameters of temperature and precipitation were downscaled using the BCSD method for RCP2.6, RCP4.5, RCP6.0, and RCP8.5 emission scenarios, and for the three time periods, i.e. near- (2014–2042), mid- (2043–2071), and far- (2072–2100) futures. The MIROC-ESM model was chosen among the latest extracted models from the AR5. Regarding the precipitation data as an important parameter in investigating the climate change, this model showed the highest compliance with observational rainfall data, based on the evaluation metrics such as NS (NS = 0.95), percent of bias (PBIAS = −2.88), coefficient of determination (R² = 0.97), and the ratio of root-mean-square-error to the standard deviation of the measured records (RSR = 0.33) (Afshar et al. 2017). The outputs of the downscaled MIROC-ESM model (rainfall and minimum and maximum temperature components in future periods) were initially fed into the validated SWAT model, and the predicted runoff for future intervals were extracted.

The annual trend analysis of the AR5 model outputs under four RCPs is presented in the Table 6 within the three future periods for hydro-climatic variables. Utilizing the MK test, it can be established that the alterations in rainfall records are decreasing in all RCPs, but their trend is not statistically significant at the confidence level of 95%. According to Table 6, mean temperature variations follow a significant trend in some intervals under different scenarios. In the near- and mid-future periods, there are increasing changes under the RCP2.6 scenario, but the trends are not statistically significant at the 5% level. In the far-future, a significant increasing trend is observed under the RCP2.6 scenario, meanwhile in far-future under the RCP4.5 scenario, there are increasing changes, but the trends are not statistically significant at the level of 5%. In the mid- and far-future periods under the RCP6.0 scenario, a significant increasing trend has been observed. Finally, in the mid-future, under the RCP8.5 scenario, there is a significant increasing trend. However, the increasing changes in the near- and far-future periods are not statistically significant at the confidence level of 95%. The analysis of the annual trend of variables indicates that the amount of precipitation will decrease in this watershed during the future periods by the end of the 21st century. The most decreasing alterations in the rainfall and the highest increase in the temperature are achieved under the highest concentration of greenhouse gases, i.e. RCP8.5. Moreover, in the near-, mid-, and far-futures, the runoff changes are decreasing under the RCP2.6 scenario, but the trend is not statistically significant at the level of 5%. In the mid- and far-future periods under the RCP4.5 scenario, there is a statistical significant decreasing trend in runoff volume; however, the decreasing variation in the near-future is not significant at the confidence level of 95%. In the near-, mid-, and far-futures under the RCP6.0, runoff variations are declining, but the trend is not statistically significant. In the far-future period, under the RCP8.5, there is a significant decreasing trend; however, in the near- and mid-futures, runoff declining changes is not statistically significant at the 5% level. Reduced rainfall and also increased temperature in the watershed will reduce the rate of runoff in the future periods in such a way that the security of the inhabitants of the region will be severely affected. The results of the studies by Afshar et al. (2017) and Gebre & Ludwig (2015) also support and confirm these findings (Gebre & Ludwig 2015; Afshar et al. 2017).

The performance evaluation of the SWAT model for runoff simulation was done using the R² and NS metrics. The R² and NS coefficients during the calibration step were estimated to be 0.52 and 0.47, respectively, while these measures during the validation process were 0.46 and 0.42, respectively. Results showed that the SWAT performance for the simulation of monthly runoff in the Zoshk-Abardeh watershed was not satisfactory, but it was in an acceptable range. Then, the SWAT model was employed to simulate runoff volume based on the downscaled rainfall and temperature values of the validated MIROC-ESM model. The outputs of this AR5 model demonstrated the highest agreement (with the R², NS, PBIAS, and RSR of 0.97, 0.95, −2.88, and 0.33 for rainfall, respectively) with observed historical records. The BCSD downscaling approach was used to extract rainfall and temperature values. The climate change simulation indicated a decreasing trend for precipitation variations in all future periods, but this trend was not statistically significant at the 5% level. The temperature variable in all RCPs had an increasing trend. The minimum alterations of temperature were obtained under the RCP2.6 scenario during the near- and far-future periods. However, temperature trend analysis under the RCP4.5 scenario during the near- and mid-futures and under the RCP6.0 scenario during the near-, mid-, and far-futures showed a significant upward trend at the 95% of confidence level. Runoff variations (simulated by the SWAT model) under the RCP4.5 scenario during the mid- to far-future and under the RCP8.5 scenario during the far-future period followed a significant downward trend. Runoff volume during the near-future period under the RCP4.5 scenario and throughout the near- to mid-future under the RCP8.5 scenario had declining variations, but its trend was not statistically significant at the 5% level. In general, these results indicated that the amount of temperature would follow an increasing tendency, while precipitation and runoff volume would follow a decreasing movement in the Zoshk-Abardeh watershed by the end of the 21st century. The increasing drift of temperature, and in particular the minimum temperature, can affect the evapotranspiration rates as well as future snowfall in the region. Thus, the amount of runoff will be impacted and decreased. Consequently, these hydro-climatic changes will impact the watershed biodiversity. Finally, the investigation of climate change impacts on groundwater, land use, and land cover condition can be suggested for further researches in the Zoshk-Abardeh watershed, as well.

The authors acknowledge the Regional Water Authority of Khorasan Razavi province for supplying necessary hydro-climatic data.

All authors declare no conflict of interest.