Evaluation of Statistical-Downscaling/Bias-Correction Methods to Predict Hydrologic Responses to Climate Change in the Zarrine River Basin, Iran

Abstract

Modeling the hydrologic responses to future changes of climate is important for improving adaptive water management. In the present application to the Zarrine River Basin (ZRB), with the major reach being the main inflow source of Lake Urmia (LU), firstly future daily temperatures and precipitation are predicted using two statistical downscaling methods: the classical statistical downscaling model (SDSM), augmented by a trend-preserving bias correction, and a two-step updated quantile mapping (QM) method. The general circulation models (GCM) input to SDSM are climate predictors of the Canadian Earth System Model (CanESM2) GCM under the representative concentration pathway (RCP) emission scenarios, RCP45 and RCP85, whereas that to the QM is provided by the most suitable of several Climate Model Intercomparison Project Phase 5 (CMIP5) GCMs under RCP60, in addition. The performances of the two downscaling methods are compared to each other for a past “future” period (2006–2016) and the QM is found to be better and so is selected in the subsequent ZR streamflow simulations by means of the Soil and Water Assessment Tool (SWAT) hydrological model, calibrated and validated for the reference period (1991–2012). The impacts of climate change on the hydrologic response of the river basin, specifically the inflow to the Boukan Reservoir, the reservoir-dependable water release (DWR), are then compared for the three RCPs in the near- (2020–2038), middle- (2050–2068) and far- (2080–2098) future periods assuming (1) the “current” consumptive demand to be continued in the future, and (2) a more conservative “recommended” demand. A systematic future shortage of the available water is obtained for case (1) which can be mitigated somewhat for (2). Finally, the SWAT-predicted ZRB outflow is compared with the Montana-based estimated environmental flow of the ZR. The latter can successfully be sustained at good and fair levels for the near- and middle-future periods, but not so for the summer months of the far-future period, particularly, for RCP85.

1. Introduction

Climate predictions of general circulation models (GCMs) indicate a profound variability of future precipitation and temperatures in terms of trends and extremes which, eventually, will increase the likelihood of severe and irreversible impacts on the ecosystem, including rivers and lakes [1]. Historical observations over the last century have also confirmed the impacts of a changing climate, as the average surface temperature increased globally by about 0.74 °C between 1905 and 2005 and, moreover, the globally averaged warming rate has almost doubled over the last 50 years, 1965 to 2005 [2]. This global warming trend will most likely continue in the near future owing to the ongoing current emissions of greenhouse gases [3].

As climate change has direct and indirect impacts on the hydrologic cycle, studies in that regard are of utmost importance for securing the future sustainability and appropriate management of water resources in a region (e.g., [4,5]).

The pressures of climate change on hydrology and water availability, in particular, are especially high in countries across the Middle East and North Africa (MENA) region, including Iran, as the most water-scarce regions in the world, not least due to high population growth [6] which leads to increasing water demand, while often facing, at the same time, poor water resources management. Because of this situation, countries like Iran are highly vulnerable to the impacts of climate change [7], including a higher frequency of natural disasters and food insecurity.

Regional hydrological climate change impact studies need a fine-scale spatial resolution of the hydro-climatic parameters, which can be provided through the downscaling of GCM outputs, and so bridge the gap between the large-scale climate parameters of the GCM models, the predictors, and the catchment-scale of local hydro-climatic variables, the predictands.

Several statistical downscaling (SD) methods are available for that purpose and have been used in regional climate impact studies (e.g., [8,9,10]). In a SD method it is assumed that the derived empirical relationships between the raw GCM outputs (predictors) and the local past observational climate data (predictands) are also valid for the future [11]. The success of a SD method is dependent on the quality of the data used for calibration, selection of the effective predictors, and the choice of the empirical transfer function [8]. This holds particularly true for the prediction of rainfall, which is often known to be unsatisfactory, when the classical statistical downscaling models (SDSM) are used (e.g., [12,13]).

Furthermore, to achieve better water management it is vital to have a steady and dependable water supply for different uses. For a reservoir operation, a dependable water supply is defined as “the reservoir-dependable water release (DWR)” which is the maintainable yield of water that can be released in a guaranteed manner over a certain time period without shortages to the reservoir system and regional environmental requirements [14].

With this background, in the present study the impacts of future climate change signals on the hydro-climatic variables are assessed for the Zarrine River Basin (ZRB) in north-western Iran. To that avail, two SD methods, the classical SDSM [15] and the non-stationary quantile mapping (QM) method, with an empirical cumulative distribution function (ECDF) transfer function [10], have been used to predict downscaled climate variables of minimum and maximum temperatures and precipitation for the near future (2020–2038), middle future (2050–2068) and far future (2080–2098) for different representative concentration pathways (RCP) emission scenarios. These RCPs are a set of four new pathways adopted by the Intergovernmental Panel on Climate Change (IPCC) and describe four different 21st-century trajectories of GHG emissions and atmospheric concentrations, pollutant emissions and land use [3]. The ensuing climate parameters are then used as drivers in the basin-wide hydrologic simulation model, the Soil and Water Assessment Tool (SWAT), to assess the future hydrology and the available water resources of the ZRB, with emphasis on the prediction of the DWR of the Boukan Dam and the outflow of the basin to the Lake Urmia (LU), which has been desiccating enormously in recent years and is on the brink of a significant environmental disaster [16], as well as the inflow of the Boukan Dam, which is supplying most of the agricultural irrigation- and potable water demands of the region.

2. Study Area

The ZRB is the most important catchment and the main inflow source of the LU, the largest salt-water lake in the Middle East and Iran, and has been shrinking in recent decades. The ZRB, with a total area of 12,025 km², is located in the southern part of the lake basin, between 45°46′ E to 47°23′ W longitude and 35°41′ S to 37°44′ N latitude and includes parts of Kurdestan, and the West and East Azarbaijan provinces, with the four big cities Miandoab, Shahindej, Takab and Saghez (see Figure 1a).

The main channel of the Zarrine River has a total length of about 300 km and most of its course is within a mountainous area. The region’s climate varies from semi-wet cold or wet-cold in the mountain areas to semi-dry in the vicinity of LU, and the annual rainfall in the basin varies between 200 mm/yr in the lower catchment area and 800 mm/yr in the mountains. The snow height varies from 5–63 mm/yr on average, and the maximum snowfall is observed mostly in the south and west of the basin. The Boukan Reservoir is the most important and largest operating dam of the basin, with a gross storage capacity of 760 million m ³ (MCM) and a live storage capacity of 654 MCM. Its water is used for agricultural irrigation and the supply of drinking water (110 MCM/yr). The agricultural areas within the basin cover a total area of 74,318 ha, all irrigated by groundwater and surface water resources, including from the Boukan Reservoir, as the crop-growing season of the basin is mostly during the dry season lasting from mid-spring to mid-autumn.

3. Materials and Methods

3.1. Climate Change Scenarios and Predictions

Even though global climate models (GCM) may run at increasingly high resolutions in the future, i.e., could serve by themselves as dynamical down scaling models, present-day climate impact studies still require some classical downscaling. For the simulation of the future water availability of the ZRB, as well as of the Boukan Reservoir inflow and of the environmental flow to LU, climate predictors for minimum and maximum temperatures and precipitation of several GCM models are downscaled, and the corresponding predictands inputted subsequently to the hydrological model, SWAT. More specifically, two downscaling models, the well-known SDSM model and an updated two-step QM method of bias correction, are employed to that end and their results will be compared for selecting the best predictions of future climate from a wider range of climate change scenarios.

3.1.1. Selecting the Climate Model Intercomparison Project Phase 5 (CMIP5) General Circulation Model (GCM)

An assessment of the future impacts of climate change is heavily dependent on the projections of the climate model used. The first step of the selection of one (or more) appropriate GCM models is based on the availability of the data, the climate scenarios and the grid resolution of the model for a particular region. In a second step, the performances of suitable GCM models are compared to each other using different statistical approaches (e.g., [17,18]).

For use in the SDSM downscaling method, the GCM’s output data archive contains two sets of 26 large scale atmospheric variables: observed/assimilated climate data of the NCEP/NCAR reanalysis project (a joint project between National Centers for Environmental Prediction and National Center for Atmospheric Research) and future atmospheric predictors. In this study, the Canadian Earth System Model (CanESM2) with a latitude x longitude grid resolution of 2.81° × 2.81° has been selected, as it is the only GCM model from the Climate Model Intercomparison Project Phase 5 (CMIP5) archive whose daily large-scale atmospheric predictors can be directly fed into SDSM, as well as the 1961–2005 assimilated prediction data of the reanalysis project. In addition, the future predictor datasets (two RCP emission scenarios, RCP45 and RCP85, for the 2006–2100 period) taken from the Canadian Center for Climate Modeling and Analysis (http://climate-scenarios.canada.ca) are employed. The SDSM methodology is then applied three times in a similar manner to the maximum and minimum temperature and the rainfall at the six synoptic climate gauges located in the watershed (see Figure 1b).

For the application of the QM downscaling method, the performances of six different GCM- models from the CMIP5 archive (some of the available GCMs have different grid resolutions and sources, as listed in

Table 1. Details of the GCM models of the Climate Model Intercomparison Project Phase 5 (CMIP5) project used here for QM downscaling.

GCM Model	Atmospheric Grid		RCP45	RCP60	RCP85	Source
	Latitude	Longitude
CESM-CAM5	0.942	1.25	✓	✓	✓	National Science Foundation, Department of Energy, National Center for Atmospheric Research, USA
CCSM4	0.942	1.25	✓	✓	✓	National Center for Atmospheric Research, USA
MRI-CGCM3	1.121	1.125	✓	✓	✓	Meteorological Research Institute, Japan
MIROC5	1.401	1.406	✓	✓	✓	Atmosphere and Ocean Research Institute, University of Tokyo, Japan
CNRM-CM5	1.401	1.406	✓	✘	✓	Centre National de Recherches Meteorologiques, France
CMCC-CM	0.75	0.75	✓	✘	✓	Centro Euro-Mediterraneo per I Cambiamenti Climatici, Italy

) over the ZRB are analyzed for different RCP scenarios and, using a multi-criteria scoring approach (e.g., [19]), the most suitable GCMs are selected. More specifically, the past-performance of a GCM, i.e., its ability to simulate the climate of the reference period in the study region, is assessed based on six performance indicators, namely, system bias (B) (ideally close to 1 [20]), coefficient of determination (R²), (ideally close to 1), normalized root mean square error (NRMSE) (ideally close to 0), and quantile biases (QB_t) (ideally close to 1 [20]) at 25th, 50th and 75th percentile thresholds of the observations. These terms are defined in Equations (1)–(4) (e.g., [19,20]):

$$B=\frac{{\sum }_{i=1}^n\left(p_i\right)}{{\sum }_{i=1}^n\left(o_i\right)}$$

(1)

$$R^2={\left(\frac{\left({\sum }_{i=1}^n\left(o_i-\overline{o}\right)\times \left(p_i-\overline{p}\right)\right)}{\sqrt{{\sum }_{i=1}^n{\left(o_i-\overline{o}\right)}^2}\times \sqrt{{\sum }_{i=1}^n{\left(p_i-\overline{p}\right)}^2}}\right)}^2$$

(2)

$$NRMSE=\frac{RMSE}{\overline{o}}=\frac{\sqrt{\frac{{\sum }_{i=1}^n{\left(p_i-o_i\right)}^2}{n}}}{\overline{o}}$$

(3)

$$Q{B}_t=\frac{{\sum }_{i=1}^n\left(p_i|p_i\ge t\right)}{{\sum }_{i=1}^n\left(o_i|o_i\ge t\right)}$$

(4)

where o_i and p_i are the observed and predicted GCM weather data and $\overline{o}$ and $\overline{p}$ are their averages, respectively, t is the threshold (e.g., 25th percentile), and n is the number of years in the reference period. The raw climate data for the analysis were taken from the German Climate Computing Center (DKRZ) (https://esgf-data.dkrz.de). It should also be noted that other GCM models of CMIP5 were excluded from the analysis, due to either a lack of available data or to similarity of the grid resolution or the source. The schematic of the grid cells for the selected GCMs for use in the QM method is shown in Figure 1c.

Based on this approach, the climate predictors of the first four models (CESM-CAM5 to MIROC5) of Table 1 are employed for further QM downscaling, using RCPs 45, 60 and 85. The most benevolent RCP26 scenario is excluded from the scenarios, as there is general agreement in the scientific community that it is practically impossible to be reached in the near future (e.g., [21,22]). In fact, the current CO₂ emissions have globally been tracking marginally even above the RCP85 scenario over the past three decades (1980–2010), but could still be pulled down by sustained mitigation from the major emitters, to make the three emission scenarios feasible [23].

3.1.2. Statistical Downscaling Model (SDSM)

Downscaling techniques are divided into two main groups: Statistical and dynamical. In situations where low-cost, rapid assessments of localized climate change impacts are required, statistical downscaling is currently still the most effective method.

For assessing regional impacts of climate change, a hybrid of multi-linear regression (MLR) and stochastic weather generator (SWG) methods, as implemented in the well-known SDSM model [15,24], is used for downscaling the daily minimum and maximum temperatures and the precipitation, using the coarse-resolution CMIP5 CanESM2 GCM predictors (see Section 3.1.1).

In the SDSM downscaling procedure, a regression equation is developed between a few selected large-scale atmospheric predictors and the local scale predictands at the meteorological synoptic stations, using the efficient dual simplex algorithm. These large-scale relevant predictors are selected based on correlation analyses and scatter plots and also considering the physical sensitivity between selected predictors and relevant predictands for the region.

Ideally, predictor-variable candidates should be physically and conceptually sensible with respect to the predictand and accurately modeled by the GCM. For downscaling rainfall, in particular, it is also recommended that the selected predictors should describe atmospheric circulation properties, such as thickness, stability and moisture content. However, in practice, the choice of predictor variables is constrained by data availability from the GCM archives.

Specifically, the structure and operation of SDSM includes five distinct tasks: (1) preliminary screening of potential downscaling predictor variables; (2) assembly and calibration of SDSM; (3) synthesis of ensembles of current weather data using observed predictor variables; (4) generation of ensembles of future weather data using GCM-derived predictor variables; and (5) diagnostic testing/analysis of observed data and climate change scenarios.

The setup and calibration–validation of the SDSM model is done here using the SDSM decision support tool (version 4.2) (see Wilby and Dawson (2007) [24], for details). However, because of systematic differences between downscaled and observed temperatures, as an additional step a trend-preserving bias correction method [25], has been applied.

3.1.3. Updated Quantile Mapping (QM) Bias Correction Technique

The raw GCM climate predictors are generally subjected to biases in both the amount and the frequency of a climate variable when compared to local climate variable for a reference period. As mentioned, this holds particularly for precipitation. In view of these deficiencies, a new class of statistical downscaling methods, so-called QM bias-correction methods have been proposed in recent years (e.g., [26]) which appear to do a better job than classical SDSM in predicting local climate.

The fundamental basis of all bias-correction methods is to quantify possible biases by comparing the distribution of the GCM- simulated climate predictors with that of the local observed variables in the historical reference period. Moreover, the main assumption in the application of the classic QM is the stationarity of the biases and other distributional properties, i.e., variance and skewness, which means that these parameters should be the same for the reference and future periods, whereas only the mean of the distribution may change [27].

As there is now enough evidence that the statistics of the climate variables will not be necessarily stationary, non-stationary QM bias-correction techniques have recently been proposed [28,29] which have been demonstrated to be more effective tools for decreasing the bias of GCM predictors than the earlier stationary QM methods.

In this study a two-step updated version of the non-stationary QM method is employed. In the first step, the classical trend-preserving bias correction method [25] is applied to correct the GCM monthly biases of temperatures and precipitation. For the temperature, an additive correction constant $C{T}_j$ defined for each month of a year (j = 1 to 12), defined as:

$$C{T}_j=\frac{{\sum }_{i=1}^n{T}_i^{OBS}-{\sum }_{i=1}^n{T}_i^{GCM}}{n}$$

(5)

is applied resulting in corrected values ${\tilde{T}}_{i,j}^{GCM}$ of the temperature:

$${\tilde{T}}_{i,j}^{GCM}=T_{i,j}^{GCM}+C{T}_j,$$

(6)

whereas for the precipitation a multiplicative factor $C{P}_j$ defined as:

$$C{P}_j=\frac{{\sum }_{i=1}^n{P}_i^{OBS}}{{\sum }_{i=1}^n{P}_i^{GCM}}$$

(7)

is applied, i.e., the corrected precipitation is:

$${\tilde{P}}_{i,j}^{GCM}=P_{i,j}^{GCM}\times C{P}_j.$$

(8)

In the second step a new, updated quantile mapping method is used to correct the daily biases of the climate variable in each month of the year. The modified non-stationary cumulative distribution function (CDF)-matching (CNCDFm) technique and equidistant CDF matching (EDCDFm) [26,29] are used, respectively, for the bias correction of the daily precipitation and temperatures (minimum and maximum). For further details of the QM method, the reader is referred to Emami and Koch (2017) [10].

3.2. Integrated Hydrological Simulation Model, Soil and Water Assessment Tool (SWAT)

The SWAT model is a physically-based semi-distributed interdisciplinary river basin model that has proven to a be a very efficient tool for assessing quantity and quality of surface water resources for a different range of scales and conditions all over the world [30,31]. In recent years the SWAT hydrological model is increasingly being used to investigate the impacts and responses of climate change on the hydrologic cycles and water resources in various regions [32,33,34], including in Iran [10,13,35]. In such applications, the meteorological driver inputs to SWAT are usually downscaled climate predictors from GCMs.

In the following two sub-sections the set-up and the methodology of the calibration and validation of the SWAT-model for the ZRB study region is summarized.

3.2.1. Soil and Water Assessment Tool (SWAT) Model Setup and Process Simulation

In the first step of the SWAT model setup for the ZRB, the basic required inputs include the digital elevation model (DEM), the river stream network, land-use and soil maps, all of which were prepared with the ArcGIS interface of SWAT (ArcSWAT). To represent the large-scale spatial variability of the watershed, the model area is divided into sub basins which are further partitioned into a set of HRUs (hydrologic response units) representing unique combinations of land-use, and soil and management practices. In this study the spatial variability of the watershed is modeled in the SWAT model using 11 sub-basins (as shown in Figure 1b) and 908 HRUs.

The SWAT model requires quite a large number of input parameters and data which were gathered from different sources. The DEM (Figure 2a) was produced by the Iranian surveying organization and has a spatial resolution of 85 m. The land-use map of the watershed, reflecting the situation in the year 2007, was obtained from the agricultural statistics and the Information Center of the Ministry of Agriculture. It has a resolution of 1000 m and delineates 10 land use classes, indicating that the major land cover of the basin is 45% dry land, 37% pasture and 14% forest (see Figure 2b). The soil map of watershed was extracted from the Food and Agriculture Organization (FAO) digital soil map of the world (see Figure 2c) at a spatial resolution of 10 km and includes 8 types of soil with two layers for the ZRB study area.

The main crops in the study region contain wheat, barley, alfalfa, potato, tomato, sugar beet, and apple. The current applied irrigation efficiency is about 38% for the areas irrigated by the surface water resources from the dam and the river, and 50% for the areas using groundwater resources [36].

The climate input data required to drive the SWAT model include daily maximum and minimum temperatures and daily precipitation of six synoptic stations located in or close to the ZRB (see Figure 1b) and were obtained from the Iranian Meteorological Organization (IRIMO) over the period of 1987 to 2015. Missing data in the records were filled using the inverse distance weighting (IDW) interpolation method which has been found in [37] as an efficient method to that end.

The other climate drivers of the hydrologic cycle including solar radiation, wind and relative humidity are generated by the weather generator (WGEN) of the SWAT from the monthly averages, because of the unavailability of daily data. Also, daily streamflow data of six gauging stations of the Zarrine River (see Figure 1b), required for the calibration and validation of the hydrologic model, were obtained for the period 1987 to 2012 from the Iran Ministry of Energy.

Using the above mentioned data sets, SWAT solves the water balance equation for the land phase of the hydrological cycle on a daily time step by routing the accumulated water through the different sections of the surface and sub-surface compartments of the water cycle to the major stream tributary, based on the spatial characteristic for each HRU and aggregated for each sub-basin [31,38]. The various simulated hydrological outputs of the SWAT model, in physical process order, include then snowmelt (prevalent in the northern mountain areas of the ZRB), evapotranspiration (computed using the Hargreaves method), surface runoff (computed using the modified Soil Conservation Service (SCS) curve number (CN) procedure), lateral subsurface flow (the fast component of surface water infiltration), return flow from shallow aquifers (from percolated water) to the streams, deep groundwater recharge from both percolation and the shallow aquifer (=loss for the streamflows) and, eventually, the accumulated streamflow (=water yield) (routed along the stream-section using the Muskingum method), after subtraction of stream transmission losses.

3.2.2. Calibration and Validation

Although numerous publications deal with the calibration and validation of the SWAT model, these two steps remain still challenging tasks, given the large number of input and hydrological “tuning” parameters required to run the model [31]. The procedures to do this vary from the classic manual calibration by trial–error to the automated complex and global optimization algorithms, such as shuffled complex evolution (SCE), particle swarm optimization (PSO) or genetic algorithm (e.g., [39]).

In this study, a semi-automated approach, the stochastic sequential uncertainty fitting version 2 (SUFI-2) algorithm, embedded in the SWAT-CUP decision-making framework [40] is employed for calibration and uncertainty analysis of the SWAT model. This procedure has been widely used in recent SWAT-related publications [10,38].

In SUFI, the input parameters of the SWAT model are expressed as ranges to be treated as random variables, accounting for all possible sources of uncertainties. As initial parameters, default values and/or determined parameters from the manual calibration or deduced parameters from the literature are taken. Then the most effective parameters are defined based on the sensitivity analysis of the SWAT-CUP. In this study, a sub-basin-wise sensitivity analysis is done hierarchically from the last streamflow gauging station in the upstream sections (sub-basin 11, Sonnateh station) down to the outlet of the total basin (sub-basin 2, Nezamabad station).

The SUFI-2 procedure assigns an uncertain range to each parameter in each HRU and each sub-basin which will then be improved iteration by iteration. The best parameter ranges of the model are then found using the 95% prediction uncertainty (95PPU) of the objective function which is a measure of the goodness of the fit of the model to the observed streamflow. Different kinds of objective functions can be used in SUFI-2/SWAT-CUP, wherefore it has been indicated in [41], that, for a reliable calibration and validation of the model, a combination of different efficiency criteria, such as the R² (coefficient of determination), the Nash–Sutcliff efficiency (NSE) and bR², should be considered.

3.3. Reservoir Dependable Water and Environmental Water Rights

The main objective of a sustainable water-management planning is to provide a minimal dependable water supply for all stakeholders. To do this properly in a water-deficient region like the ZRB, the reservoir releases must meet the potential demand, even during extreme drought periods, i.e., the capacity of the former must be sufficiently high. In that regard, the criterion of “the reservoir dependable water release (DWR)” comes into play, which is the maintainable yield of water that can be released in a guaranteed manner over a certain time period, without shortages to the reservoir system and the regional environmental requirements [14]. These reservoir releases are simulated by the reservoir module of the SWAT model and an average annual release rate for the different future periods/scenarios is then computed and compared with the currently projected water demands in the ZRB—which, in decreasing order of priority, are potable (Pot), industrial (Ind), environmental (Env) and agricultural (Agr) demands [36]—as well as with a recommended alternative of reduced future demands [42]. The statistical analysis of the flow changes for the different future periods/scenarios is done for three percentiles (25th, 75th and 90th) and the long-term average.

Furthermore, the environmental water demands of the Zarrine River, as the main source of LU, are calculated using the Tennant (Montana) method [43] at two levels, “fair or degrading” and “good”. In this method, a relationship between habitat suitability and the proportion of the river streamflow allocated to satisfy the environmental demand is assumed, with the former defined from 20% (fair condition) to 40% (good condition) in spring and summer (April to September) and from 10% to 20%, respectively, for autumn and winter (October to March), as recommended by Tennant (1976) [43] and also used in Ministry of Energy (MOE) (2014) report [36].

4. Results and Discussion

4.1. SDSM Downscaling of Precipitation and Temperatures

Calibration and Validation of the SDSM Model

The SDSM-downscaling model is calibrated on observed data for the 13-year period 1987–1999 and validated for the 5-year period 2000–2005. Summaries of the statistical results for the calibration and validation of the model for the two (min and max) temperatures and the precipitation, using the standard error (SE) and the coefficient of determination (R²) as goodness-of-fit measures [24], are listed in

Table 2. Goodness-of-fit measures (based on standard error (SE) and R²) of the statistical downscaling model (SDSM), for min. and max. temperatures in the calibration and validation periods for the different climate stations in the ZRB.

Station	Tmp Max.				Tmp Min.				Effective Predictors
	CAL		VAL		CAL		VAL
	R2	SE	R2	SE	R2	SE	R2	SE
St1	0.77	2.1	0.94	2.8	0.63	2.1	0.90	2.6	500 hPa Geopotential 850 hPa Specific humidity 1000 hPa Specific humidity Air temperature at 2 m
St2	0.80	1.7	0.96	2.4	0.75	1.7	0.94	2.2	850 hPa Specific humidity 1000 hPa Specific humidity Air temperature at 2 m
St3	0.79	1.9	0.95	2.7	0.60	2.9	0.81	3.5	500 hPa Geopotential 1000 hPa Specific humidity Air temperature at 2 m
St4	0.81	1.7	0.96	2.3	0.63	2.4	0.87	2.9	500 hPa Geopotential 850 hPa Specific humidity 1000 hPa Specific humidity Air temperature at 2 m
St5	0.78	1.8	0.96	2.4	0.65	2.1	0.91	2.5	500 hPa Geopotential 850 hPa Specific humidity 1000 hPa Specific humidity Air temperature at 2 m
St6	0.82	1.7	0.95	2.8	0.71	1.7	0.90	2.6	500 hPa Geopotential 850 hPa Specific humidity 1000 hPa Specific humidity Air temperature at 2 m

and

Table 3. Similar to Table 2, but for the precipitation.

Station	CAL		VAL		Effective Predictors
	R2	SE	R2	SE
St1	0.33	4.8	0.35	3.2	500 hPa Meridional wind component 500 hPa Geopotential height 850 hPa Divergence of true wind 500 hPa Specific humidity Total precipitation
St2	0.28	4.7	0.26	2.6	500 hPa Meridional wind component 850 hPa Divergence of true wind 500 hPa Specific humidity Total precipitation
St3	0.3	6.3	0.30	3.6	500 hPa Meridional wind component 850 hPa Divergence of true wind Total precipitation
St4	0.26	2.7	0.20	2.0	1000 hPa Divergence of true wind 500 hPa Specific humidity 850 hPa Specific humidity Total precipitation
St5	0.25	5.3	0.28	3.1	1000 hPa Divergence of true wind 500 hPa Specific humidity 850 hPa Specific humidity Total precipitation
St6	0.49	2.7	0.15	2.8	500 hPa Specific humidity 850 hPa Specific humidity Total precipitation

, respectively.

Based on the results in the two tables, the derived predictor–predictand relationships are considered quite satisfactory, although the performance of the SDSM model is better for the min. and max. temperatures than for the precipitation, for which the biases are still considerable (a well-known common deficiency of the SDSM method [12]). For that reason, the latter is corrected in an additional step using the trend-preserving bias correction method [25], somewhat similar to the QM, but with different CP_j correction factors in Equations (3) and (4).

4.2. QM-Downscaling of Temperatures and Precipitation

4.2.1. Selecting the GCM Model for QM

Using the set of the first four GCMs of Table 1 as eligible candidates for GCM selection, their output of climate predictions are compared with the corresponding observations using the bias and goodness-of-fit indicators at the six weather stations of the basin (see Figure 1b) in the time period 1987–2005 and, based on their skill performances, the GCM models are ranked accordingly. The results of the scoring analysis for the average of the two temperatures (min. and max. temperature scores were very similar) and for the precipitation are summarized in

Table 4. Performance skill scores of the selected GCM models of CMIP5 for the max and min temperatures with the average score highlighted in bold.

Factor		Model	St1	St2	St3	St4	St5	St6	Skill Score	Final Score
B		CESM-CAM5	0.83	1.02	1.72	2.30	1.97	1.03	2	2.5
		CCSM4	0.88	0.91	1.70	2.17	1.93	0.94	4	3.0
		MRI-CGCM3	0.90	0.67	1.89	2.05	2.47	0.87	3	3.3
		MIROC5	1.53	1.14	3.05	2.80	3.38	1.16	1	1.2
NRMSE		CESM-CAM5	0.68	0.61	1.74	2.05	2.09	0.66	2
		CCSM4	0.64	0.59	1.57	1.90	1.96	0.63	3
		MRI-CGCM3	0.30	0.43	1.21	1.27	1.68	0.27	4
		MIROC5	1.85	1.14	5.03	4.69	5.86	1.30	1
R2		CESM-CAM5	0.58	0.58	0.58	0.56	0.53	0.55	4
		CCSM4	0.61	0.61	0.61	0.59	0.58	0.58	3
		MRI-CGCM3	0.92	0.90	0.92	0.92	0.88	0.88	1
		MIROC5	0.64	0.64	0.66	0.64	0.59	0.61	2
QBt	25th	CESM-CAM5	0.92	1.04	1.39	1.63	1.38	1.07	2
	50th		0.98	1.05	1.34	1.50	1.29	1.07	2
	75th		1.03	1.08	1.32	1.43	1.26	1.09	3
	25th	CCSM4	0.96	0.93	1.37	1.50	1.40	0.96	3
	50th		1.02	0.80	1.34	1.43	1.33	1.00	3
	75th		1.09	1.03	1.35	1.40	1.33	1.05	2
	25th	MRI-CGCM3	0.95	0.79	1.41	1.47	1.54	0.91	4
	50th		1.00	0.85	1.36	1.38	1.43	0.95	4
	75th		1.07	0.93	1.35	1.36	1.40	1.00	4
	25th	MIROC5	1.39	1.09	1.99	1.80	1.90	1.10	1
	50th		1.34	1.07	1.76	1.62	1.66	1.08	1
	75th		1.34	1.07	1.64	1.52	1.57	1.08	1

and

Table 5. Similar to Table 4, but for the precipitation.

Factor		Model	St1	St2	St3	St4	St5	St6	Skill Score	Final Score
B		CESM-CAM5	0.94	1.04	0.81	1.47	0.93	1.02	6	5.7
		CCSM4	1.23	1.27	0.95	1.76	1.09	1.24	4	3.7
		MRI-CGCM3	1.51	2.11	1.23	2.45	1.42	1.61	2	2.7
		MIROC5	0.77	1.58	0.68	1.91	1.10	1.21	5	4.7
		CNRM-CM5	0.94	2.14	0.83	3.13	1.81	1.63	1	1.8
		CMCC-CC	1.48	1.13	1.34	2.01	1.42	0.86	3	2.5
NRMSE		CESM-CAM5	1.25	1.19	1.26	1.62	1.30	1.11	6
		CCSM4	1.43	1.31	1.32	1.87	1.40	1.31	4
		MRI-CGCM3	1.57	1.95	1.34	2.56	1.40	1.50	3
		MIROC5	1.05	1.55	1.15	1.90	1.20	1.27	5
		CNRM- CM5	1.11	2.14	1.13	3.35	1.70	1.62	1
		CMCC-CC	1.65	1.52	1.77	2.42	1.83	1.40	2
R2		CESM-CAM5	0.08	0.08	0.08	0.05	0.06	0.17	4
		CCSM4	0.06	0.08	0.05	0.06	0.04	0.10	2
		MRI-CGCM3	0.12	0.16	0.08	0.10	0.11	0.12	5
		MIROC5	0.13	0.11	0.10	0.10	0.09	0.10	3
		CNRM	0.18	0.09	0.18	0.13	0.18	0.09	6
		CMCC-CC	0.16	0.00	0.00	0.00	0.00	0.00	1
QBt	25th	CESM-CAM5	0.94	1.00	0.80	1.46	0.92	1.00	6
	50th		0.95	0.95	0.81	1.47	0.95	0.95	6
	75th		0.98	0.92	0.86	1.44	1.02	0.92	6
	25th	CCSM4	1.23	1.24	0.95	1.76	1.09	1.23	4
	50th		1.21	1.16	0.94	1.73	1.10	1.16	4
	75th		1.22	1.10	0.94	1.67	1.11	1.09	4
	25th	MRI-CGCM3	1.51	2.08	1.22	2.46	1.41	1.59	2
	50th		1.45	1.88	1.15	2.37	1.35	1.44	2
	75th		1.43	1.74	1.08	2.23	1.29	1.31	2
	25th	MIROC5	0.77	1.55	0.68	1.90	1.09	1.18	5
	50th		0.75	1.48	0.65	1.84	1.05	1.12	5
	75th		0.73	1.40	0.63	1.74	1.01	1.06	5
	25th	CNRM-CM5	0.93	2.07	0.82	3.10	1.78	1.58	1
	50th		0.95	1.95	0.83	3.04	1.73	1.48	1
	75th		0.97	1.83	0.84	2.87	1.66	1.38	1
	25th	CMCC-CC	1.49	1.13	1.34	2.03	1.43	0.87	3
	50th		1.56	1.11	1.35	2.05	1.47	0.84	3
	75th		1.56	1.06	1.28	1.94	1.42	0.80	3

, respectively.

For the temperature, the assigned skill (rank) score (=4 for the model with the best performance) varies from 1 to 4 and for the precipitation from 1 to 6. The final skill score of each model is calculated as the average of the six scores (statistical factor, as listed in the first column of the table) and enumerated in the last column. Based on these average skill scores, the climate predictors of the CESM CAM5 and MRI-CGCM3 models are selected as the best choice for the QM downscaling of the future precipitation and temperatures, respectively.

4.2.2. Calibration and Validation of the QM Model

After applying the CT_j and CP_j correction factors (j = 1 to 12, number of months in a year; n = 19, total number of years) in Equations (5)–(7), the monthly biases of the GCM simulated climate variables are removed and then the QM-model is calibrated and validated for each month of the year on the daily time step during the period of 1987–2005, with 1987–1998 as the calibration and the remainder as the validation period. This is done for each of the six climate stations. The CDFs of the min. temperatures resulting from the various bias corrections steps, as epitomized in Equations (4) and (6), are shown representatively for St2 station for the months March and July in the two panels of Figure 3. Corresponding plots for the max. temperature at St4 station for January and October are presented in Figure 4.

One may notice from the two figures how the CDFs of the bias-corrected min. or max. temperature (Bias Correctd_Val) are very similar to the CDFs of the observed data during the validation period (Obs_Val), which clearly indicates that the performance of the QM model in correcting the biases for the min. and max. temperature is also quite satisfactory.

It should also be noted that the CDF of the bias-corrected variable in the calibration period is not shown in the figures, as it approximately overlaps with the CDF of “obs_Cal”, as it should.

Following the descriptions of Equations (3) and (4), the precipitation is downscaled/bias-corrected by the QM method in a similar manner. In analogy to the previous figures, the two panels of Figure 5 show the results of the calibration and validation of the precipitation at station St3 for the months of May and December. One can notice that the QM-downscaling model performs well in the validation period, as the CDF of the bias-corrected precipitation (Bias Correctd_Val) fits the CDF of the observed precipitation (Obs_Val) quite well.

4.3. 2006–2015 Projections of Temperatures and Precipitation

As the historic observed climate data at the synoptic stations (Figure 1b) were only available from 1987 to 2015, the reliability of the GCM/SDSM/QM predicted min./max. temperatures and precipitation for the period 2006–2015 is evaluated by comparing them with the observed data. To that avail, the goodness-of-fit measures, coefficient of determination (R²), standard error (SE) and index of agreement (IA) (defined as (1—standardized RMSE (Equation (3))) and so varying between 0 and 1 [44]) are computed for each climate station and the range of values is presented for the various combinations of downscaling methods/RCP-scenarios for the average of the min. and max. daily temperatures and the daily precipitation in the three plots of Figure 6. One may notice from these plots that the efficiency of the QM and SDSM methods is overall satisfactory, but the QM performs clearly much better than SDSM in the 2006–2016 evaluation period.

4.4. Future Projections of Temperature and Precipitation

Employing the various combinations of downscaling methods/scenarios as presented in Figure 6, the climate is predicted for three future periods: near (2020–2038), mid (2050–2068) and far (2080–2098).

In

Table 6. Relative changes compared with the observed period (1987–2005) of mean air temperature (°C) and precipitation (%) for the three future periods (near, mid, far) using the various downscaling model/scenarios.

Model	Future Period	Scenario	∆TMP (°C)	∆PCP (%)
SDSM	Near	RCP45	1.39	8%
		RCP85	1.58	15%
	Mid	RCP45	2.58	12%
		RCP85	3.83	22%
	Far	RCP45	3.25	13%
		RCP85	6.21	29%
QM	Near	RCP45	0.61	−4%
		RCP60	0.75	4%
		RCP85	1.18	3%
	Mid	RCP45	1.4	15%
		RCP60	1.51	3%
		RCP85	2.85	17%
	Far	RCP45	1.63	−3%
		RCP60	2.2	13%
		RCP85	4.7	−19%

the relative changes (∆) of the future predicted (average min/max) temperatures and precipitation, relative to the historical period (1987–2005), are listed for all downscaling methods/scenarios/periods combinations. As for the temperature changes in the table, in agreement with the general trend of the RCP scenarios, they are all positive and increase from RCP45 to RCP85, going beyond the 4 °C postulated globally by the year 2100 for the latter [45]. This holds particularly for the SDSM-predicted average temperature change of 6.21 °C which appears to be unrealistically pessimistic.

The situation is even more versatile for the SDSM-predicted precipitation changes. Although a global precipitation increase is likely in the future, it will not be uniform regionally [3]. Thus, it is expected that dry regions get drier and wet regions wetter [46]. The global mean temperature predicted for RCP85 will be nearly 3 °C higher than of RCP45 by year 2100, but all climate scenarios predict an increase of the precipitation by the doubling of radiation forcing in RCP85, relative to RCP45. So, according to the fluctuating trend of the historic precipitation in the ZRB, with even some dry years in the time-period 1998–2002 it may be expected that in some future periods/scenarios the precipitation will also decrease. That contradicts obviously the SDSM-predicted precipitation, but is more consistent the QM precipitation in a sensible manner.

As mentioned earlier, the classical SDSM method is often deficient in the downscaling of the precipitation [12], and for this reason a trend-preserving bias correction of the SDSM downscaled precipitation has been included here. However, Figure 7 shows that, although this add-on improves the predicted precipitation slightly, the observed precipitation is still significantly overestimated. It should also be noted that the better resolution of the GCM models employed in the QM method (CESM-CAM5 and MRI-CGCM3) than that of the SDSM-method (CanESM2) may also be the reason for the better performance of the former.

Therefore, because of the aforementioned larger discrepancies in the SDSM predicted and observed temperature and precipitation and owing to the better performance of the QM model in the recent 10-year period (2006–2015) (see Section 4.3), for the subsequent simulation of the future water resources and the hydrologic changes in the ZRB just the QM downscaled precipitation and temperature under the various scenarios are considered.

4.5. Calibration and Validation of the SWAT Model

Firstly, the initial parameters of the calibration are defined using the SWAT default parameters and information from a previous study of the authors [10]. Then the sensitivity analysis is done by means of the SWAT-CUP automatic procedure and the 24 most sensitive parameters of the SWAT model determined. The results are listed in

Table 7. The most sensitive ranked SWAT input parameters used for the different sub-basins.

Rank	Parameter Name	Dimension	Final Value/Range	Rank	Parameter Name	Dimension	Final Value/Range
1	CN2.mgt	dimensionless	59–90	13	ALPHA_BNK.rte	days	0.10–0.77
2	GW_DELAY.gw	Day	5–28	14	SNOCOVMX.bsn 1	mm	116
3	SOL_BD(1).sol	g/cm3	0.95–1.54	15	SOL_Z(1).sol	mm	281–392
4	GW_REVAP.gw 1	dimensionless	0.037	16	SNO50COV.bsn 1	dimensionless	0.57
5	RCHRG_DP.gw 1	dimensionless	0.014	17	SOL_K(1).sol	mm/h	7.22–14
6	REVAPMN.gw	mm H2O	1–450	18	TIMP.bsn 1	dimensionless	0.57
7	ALPHA_BF.gw	1/day	0.03–0.75	19	GWQMN.gw	mm H2O	750–3514
8	ESCO.hru	dimensionless	0.51–0.99	20	CH_N2.rte 1	dimensionless	0.016
9	SFTMP.bsn 1	°C	1.93	21	SMTMP.bsn 1	°C	1.08
10	CH_K2.rte 1	mm/h	0.5	22	GW_SPYLD.gw 1	m3/m3	0.05
11	SHALLST.gw	mm	943–5000	23	SMFMX.bsn 1	mm H2O/°C-day	7.95
12	SOL_AWC(1).sol	mm H2O/mm	0.14–0.24	24	SMFMN.bsn 1	mm H2O/°C-day	0.73

¹ Values of these parameters are the same or assumed to be the same for the whole basin.

.

To determine the optimal range of the named 24 model parameters, an initial uncertain range is given globally for each parameter, which is then fine-tuned in the calibration employing the SUFI-2 optimization algorithm hierarchically sub-basin-wise from the utmost upstream sub-basin (11) outlet down to the main outlet of the basin (sub-basin 2 outlet). We forgo a discussion of the details of the calibration results and refer to Emami and Koch (2017) [10]. Suffice to say that the SCS curve number (CN2), the groundwater delay time, and the moist bulk density of the soil turned out to be the three most sensitive parameters.

The SWAT-simulated monthly inflow to the Boukan Dam (the discharge at the outlet of sub-basin 8) and the main outflow of the basin (the discharge at the outlet of sub-basin 2, station Nezambad) are compared with the observed records in Figure 8 for the calibration (1998–2012) and validation period (1991–1997). It should be noted that the calibration and validation are done on the monthly scale because of the unavailability of daily data. Also, the first 4 years of the validation period (1987–1990) are used for the warm-up of the SWAT model. Figure 8 shows that the performance of the SWAT simulation for both the calibration and validation periods is quite satisfactory, although some peaks are underestimated which is mainly because of the occurrence of the some dry years with low precipitation, especially, during 1998–2002.

Quantitative statistical measures of the SWAT calibration and validation of the streamflow at the outlet stations of the six sub-basins (see Figure 1b) are presented in

Table 8. Statistical measures for the SUFI-2—optimized objective functions for calibration and validation periods for various sub-basin outlet stations.

Streamflow Station	Calibration			Validation
	R2	NS	bR2	R2	NS	bR2
Nezamabad	0.7	0.69	0.47	0.65	0.62	0.48
Chooblooche	0.50	0.48	0.40	0.56	0.43	0.44
Sarighamish	0.65	0.65	0.50	0.60	0.55	0.47
Boukan Dam	0.85	0.72	0.51	0.77	0.65	0.48
Safakhaneh	0.57	0.55	0.45	0.59	0.49	0.40
Sonateh	0.72	0.67	0.56	0.77	0.75	0.70
Average	0.67	0.63	0.48	0.67	0.58	0.50

in terms of the coefficient of determination (R²), Nash–Sutcliffe (NS) and the weighted bR² (coefficient of determination multiplied by the slope of the regression line between simulated and observed streamflow values). Again, all efficiency criteria are within the acceptable range, with the NS, as the most prominent measure, ranging between 0.5 and 0.7.

4.6. Climate Change Impacts on the Zarrine River Basin (ZRB) Hydrology and Water Resources

4.6.1. SWAT-Simulated Future Stream Inflow to the Boukan Dam

The QM-downscaled temperatures and precipitation predicted under the three RCPs for the three future periods 2020–2038 (near), 2050–2068 (mid) and 2080–2098 (far) are employed as input drivers in the calibrated and validated SWAT model to evaluate the impacts of possible climate change on the hydrologic cycle and water balance of the ZRB. The emphasis here is on the estimation of the inflow to the Boukan Dam as the main water infrastructure of the region and at the origin of the environmental outflow of the basin entering the LU.

In Figure 9, the statistical measures of the SWAT-simulated streamflow predictions at the inlet of the Boukan Dam obtained for the various QM/RCP/periods are compared with those of the simulated inflow of the dam during the historical period. One can notice that the low (25th percentile) quantile of the dam inflow will decrease for all future scenarios, whereas the high (75th percentile) and the extreme high (95th percentile) inflow quantile are more likely to increase in the future scenarios, except for the RCP85 scenario in the far future (2080–2098) period. Up to the end of the mid-future period (year 2068), the annual average of the inflow to the dam will be slightly increasing for all, but for RCP45 (near future), for which it is more likely to decrease. However, one can expect that this increase of available water will not cover that of the also increasing demands (Section 4.6.3), due to expanded agricultural and other development plans. Moreover, later in the far future, there will be a gradual decreasing trend for RCP45 and RCP85 again, although less so for RCP60, so that the water shortage in the study area will be further exacerbated.

In summary, based on Figure 9, one can conclude that the inflow to the Boukan Dam is generally decreasing in the low and average quantile range of the streamflow by the year 2100, whereas high and extremely high flow ranges exhibit slight increases by that time.

4.6.2. SWAT-Predicted Future Water Balance Parameters for the ZRB

The main SWAT-predicted ZRB water balance parameters for the historical and the three future periods are listed in

Table 9. Future (three periods) SWAT-predicted water balance parameters for the ZRB in comparison with the historical parameters for the three different RCP scenarios.

Water Balance Parameter	Historical	Near Future			Middle Future			Far Future			Average
		RCP45	RCP60	RCP85	RCP45	RCP60	RCP85	RCP45	RCP60	RCP85	RCP45	RCP60	RCP85
PCP	392	−5%	4%	3%	15%	3%	17%	−2%	13%	−19%	3%	7%	1%
SNM	58	5%	36%	−36%	−5%	−16%	−38%	−41%	−57%	−83%	−14%	−12%	−52%
ET	290	5%	3%	6%	10%	7%	19%	25%	36%	14%	13%	16%	13%
∆SW 1	17	−90	−28	−100	−61	−41	−33	−84	−46	−106	−78	−38	−80
SURQ	46	−13%	13%	−24%	11%	2%	2%	−37%	−24%	−63%	−13%	−2%	−28%
GWQ	84	−18%	5%	−2%	15%	−2%	4%	−46%	−32%	−67%	−17%	−10%	−21%
LATQ	38	−16%	−8%	21%	21%	3%	24%	−11%	13%	−26%	−3%	3%	5%
WYLD	168	−16%	4%	−2%	15%	1%	7%	−36%	−20%	−57%	−12%	−5%	−17%

¹ Difference of the averages of the soil water (mm) stored in the profile during the simulation period.

. One can notice from the table that the input precipitation (PCP) exhibits a wide range of change, varying from a −19% decrease (RCP85, far future) to a +17% increase (RCP85, middle future), together with a very slight growth of the average. Notwithstanding this, the future predicted snowmelt (SNM) and the soil water change (∆SW) will decrease remarkably on average, while the evapotranspiration (ET) is rising, i.e., altogether, there will be an unfavorable trend for the available renewable water resources in the ZRB. In fact, this can also be seen directly from the estimated streamflow in the main channel, called water yield (WYLD)—which is equal to the sum of the surface runoff (SURQ) + groundwater recharge to the reach (GWQ) + lateral subsurface flow (LATQ)—which will decrease by −5% (RCP60) and −17% (RCP85) on average over the three future horizons.

Generally, Table 9 indicates that the future impacts of climate change on the hydrological cycle in the ZRB will be non-uniform, depending on the different future periods/RCP scenarios. The RCP45, RCP60 and RCP85 can be called the “worst case” scenarios of the near, middle and far future periods, respectively, whereas the middle future and the other two periods (near and far future) RCP45 and RCP60, respectively, will be the “best case” scenarios. From the table one may also deduce that if one proceeds “as usual” (without taking climate mitigation policies), the available surface and groundwater resources of the ZRB will decline more and more (−28% for SURQ and −21% for GWQ for RCP85) on average over the three future periods.

4.6.3. Future Dependable Reservoir Water Release (DWR)

The annual dependable water releases (DWR) or safe yields of the Boukan Dam for the different future periods/RCP scenarios are determined by summing up the corresponding SWAT-simulated average monthly dam releases over the corresponding prediction period/scenarios.

The results are shown in the three panels of Figure 10. One may notice from the figure that the seasonal climate cycle in the region is clearly reflected in the variation of the dam’s water release. Furthermore, the smallest and largest differences of the future dam releases to the historical ones will be in the medium- and far-future periods, respectively. Regarding the effects of the different RCPs, the average monthly dam releases will have a slow positive (RCP45 and RCP85) or negative (RCP60) change, when going from the near to the middle future period; however, the monthly DWR decrease strongly for the far-future period. For the latter, there will be also much less seasonal variation of the dam water release than has been the case for the historical period, due to the higher future summer precipitation.

It should also be noted that in the above calculation of the future monthly DWR—as the SWAT-estimated evaporation from the reservoir seems to be underestimated and is quite low compared with the estimated values extracted from MOE (2014) [36]—the SWAT-estimated evaporation of the dam has been corrected with a historical correction factor, based on the MOE-reported average evaporation value of 59.4 MCM/yr in the historical period.

Summing up the monthly DWR of Figure 10, the annual DWR of the reservoir are obtained. The results are listed in

Table 10. Comparison of the future dependable water release (DWR) of the reservoir with the present and future water demands assuming two development plans.

Water Category		Historic	Near Future			Middle Future			Far Future
			RCP45	RCP60	RCP85	RCP45	RCP60	RCP85	RCP45	RCP60	RCP85
DWR (MCM/yr)		983	896	1169	1092	1312	1106	1212	644	857	355
Demands (MCM/yr)	Current	1030	1288			1373			1501
	Recommended		943			1007			1136
Supplied Demand (%)	Current	95%	70%	91%	85%	96%	81%	88%	43%	57%	24%
	Recommended		95%	>100%	>100%	>100%	>100%	>100%	57%	75%	31%

for the different future periods/RCP scenarios and are to be compared with the corresponding total demands [36,42]. Two kind of future demand scenarios are assumed in Table 10: (1) Current, based on predictions of up to year 2042 [36], considering a continuation of the current trends for the potable and industrial demands, together with the agricultural development plans, population growth and other characteristics of the region up to that time and linearly extrapolated for the middle and future periods; (2) Recommended, based on the recommended scenario of the potable and industrial predicted demands [36], together with a decreased agricultural demand, based on new executive strategies of demand management recommended by the Urmia Lake Restoration Committee [42]. These strategies include an increase of the irrigation efficiency, better river bed and bank management, deficit irrigation, improvement of the Zarrine irrigation network and completion of irrigation secondary networks with surface and modern techniques.

Table 11. Division into the water-use categories of the projected demands for the three future periods.

Water Category	Near Future				Middle Future				Far Future
Demands (MCM/yr)	Pot	Env	Ind	Agr	Pot	Env	Ind	Agr	Pot	Env	Ind	Agr
Current	229	383	20	656	300	383	34	656	442	383	20	656
Recommended	212	383	18	348	277	383	18	348	407	383	16	348

lists the division into the named categories of the projected demands for the three future periods. Obviously, as indicated in Table 10, the consumptive water demands (Pot, Ind and Agr) for scenario (2) will be less than for (1).

It can be seen from Table 10—and in line with the results of Table 9—that compared with the historic period, the DWR increases for the near and middle future periods, but decreases after that time for the RCP45 and RCP85 scenarios, which raises alarms about a severe water crisis in the far future.

This is mostly a consequence of the high evaporation from the Boukan Reservoir (see Table 9) and the ensuing reduction of the dam’s outflow into the ZRB. Although the future situation appears to be more positive for RCP60, for none of the RCPs investigated can the projected future water demands be met, particularly for the far-future period when the percentage ratio of DWR to demand is only 43% for RCP45 and even 24% for RCP85 for the “current” demand option.

On the other hand, for the “recommended” demanding option, these numbers go up to 57% and 31%, respectively, i.e., the water shortages will be decreased by 14%, respectively 7%, by year 2100. It is interesting to see that for the near and middle future the recommended water demands can be fully supplied (100%), with exception of the near future/RCP45 scenario, where the water demand is only satisfied by 95%.

Therefore, it is clear that not only the further agricultural development plans of the Boukan Dam should be stopped, but that the current agricultural demands should be reduced by 30% to 40% by implementing the executive strategies of the MOE (2016) report [42] demand management. In addition, the primary consumptive demands (potable and industrial) should be altered from a continuation of the current trends of MOE (2014) [36] to the recommended scenario of MOE (2016) [42], otherwise the ZRB will face a more water-scarce future.

4.6.4. Future Monthly Environmental Water Demands of the ZR

As the environmental flow is a non-consumptive demand and comes back to the river, for securing the river water rights it is necessary to control the latter at the most downstream point of the basin, i.e., the outlet of the basin. The monthly environmental water demands of the ZR, which eventually will also affect LU, are calculated using the Tennant (Montana) method (see Section 3.3) for “fair” and “good” levels and compared with the SWAT-simulated environmental water flow for the three future periods/RPCs in the three panels of Figure 11.

As can be seen from the top and middle panels of that figure, for the near- and middle-future periods the environmental rights will be satisfied at both “fair” and “good” levels for all RCPs, whereas in the far-future period (bottom panel), the river discharge will not meet the environmental demands at “fair” and “good” suitability levels during months March to June for the extreme RCP85 scenario and fails also to supply the demands in April and May for scenario RCP45. However, the demand is at least satisfied at the “fair” level for the middle future/RCP45 scenario. In contrast, for the more benevolent far future/RCP60 scenario, the regulated environmental flow will be satisfied at both “good” and “fair” levels.

In conclusion, the supply of the future environmental demands of the ZR entering the LU is overall satisfactory, but it will be of some concern in the far future.

5. Summary and Conclusions

In this study the hydrological and environmental impacts of climate change are assessed using two different statistical downscaling methods, the well-known SDSM model, but combined here with a trend-preserving bias correction methodology, and the two-step updated quantile mapping method (QM). Temperature and precipitation predictors of the CANESM2-GCM are employed for SDSM and of a set of CMIP5-GCMs for use in QM is evaluated employing a skill score multi-criteria method. Eventually, CGCM3 and CESM-CAM5 are found to be the most suitable ones for QM downscaling of temperatures and precipitation, respectively.

Both downscaling models are calibrated and validated for min./max. temperatures and the precipitation for the 1987–2005 historic base and the “future” 2006–2015 period. As the QM-predictands do better here than the SDSM ones, the QM future climate scenarios are selected for further evaluation of the hydrological impacts of climate change on the water resources of ZRB and the operation of the Boukan Dam.

To do this, the integrated hydrologic model SWAT is calibrated and validated quite satisfactory on streamflow records at various gauging stations for the 1991–2012 observation period, using the SUFI-2 algorithm. Then the model is run for the QM-predicted climate variables under RCP45, RCP60 and RCP85 scenarios for three future periods (near, middle and far future).

The results show that, compared with the historical period, the inflow of the Boukan Reservoir and the main water balance parameters of the ZRB will have mostly negative trends, especially, in the far future (2080–2098) period. Furthermore, low and mean quantiles of the flow entering the reservoir will decrease in the future periods, whereas the high and extreme high flows will slightly increase. The future surface runoff, groundwater contribution to the main reach and the total streamflow will also decline on average. Although water yield and reservoir inflow will increase slightly in the middle (2050–2068) period, they will fall gradually again in the future (2080–2098), especially for RCP45 and RCP85 scenarios.

Next, the annual dependable water release (DWR) of the Boukan Reservoir is calculated and compared with the expected future demands. Two scenarios for the latter are investigated, whereby the first is the continuation of the “current” policies and the second is a new “recommended” demand management scenario. The SWAT–DWR predictions indicate that for the “recommended” scenario all demand types will be successfully fulfilled in the near- and middle-future periods (up to the year 2068), whereas for the “current” scenario there will be a significant water shortage in this period, especially for RCP45 and RCP85. In fact, for the far future, the water supplied from the reservoir is alarmingly low for RCP45 and RCP85, as it will not be able to meet fully even the most essential demands (potable, industrial and environmental) for both demand scenarios.

Finally, the SWAT-predicted outflow of the basin and the Montana-based environmental rights of the Zarrine River are compared. The results demonstrate that the environmental flow can be successfully supplied at both “good” and “fair” levels for the near- (2020–2038) and middle- (2050–2068) future periods. Nevertheless, the river discharge cannot meet the rights for RCP45 at a good level and for RCP85 at “good” and “fair” levels during the summer months in far-future years.

In conclusion, the results indicate that the future impacts of a changing climate in the ZRB will be diverse for different periods and scenarios, with a high likelihood of shrinkage of the available freshwater resources in the basin, with adverse effects on LU, particularly, in the far-future period (2080–2098). Despite the fact that for the earlier future periods (near and middle) an increase of the water resources, i.e., precipitation, streamflow and the DWR in some scenarios/periods is predicted, the surface- and soil-water losses from evapotranspiration as well as evaporation from the Boukan Reservoir, in conjunction with a projected surge of demand, will more likely counteract this positive effect. This holds particularly when the “current” demand operational policy is perpetuated in the future. In contrast, by implementing instead the more conservative recent MOE “recommended” demand strategy, water shortages can be avoided, at least up to the middle-future period (2050–2068). Finally, it may be noted that enactment of further greenhouse reduction policies (mitigation) could be very effective for the conservation of water resources in the ZRB as well, as the results clearly show that the intermediate emission scenarios, RCP45 or RCP60, will have less future-adverse hydrological effects than the high-emission scenario, RCP85.