Agricultural Water Productivity-Based Hydro-Economic Modeling for Optimal Crop Pattern and Water Resources Planning in the Zarrine River Basin , Iran , in the Wake of Climate Change

For water-stressed regions/countries, like Iran, improving the management of agricultural water-use in the wake of climate change and increasingly unsustainable demands is of utmost importance. One step further is then the maximization of the agricultural economic benefits, by properly adjusting the irrigated crop area pattern to optimally use the limited amount of water available. To that avail, a sequential hydro-economic model has been developed and applied to the agriculturally intensively used Zarrine River Basin (ZRB), Iran. In the first step, the surface and groundwater resources, especially, the inflow to the Boukan Dam, as well as the potential crop yields are simulated using the Soil Water Assessment Tool (SWAT) hydrological model, driven by GCM/QM-downscaled climate predictions for three future 21th-century periods under three climate RCPs. While in all nine combinations consistently higher temperatures are predicted, the precipitation pattern are much more versatile, leading to corresponding changes in the future water yields. Using the basin-wide water management tool MODSIM, the SWAT-simulated water available is then optimally distributed across the different irrigation plots in the ZRB, while adhering to various environmental/demand priority constraints. MODSIM is subsequently coupled with CSPSO to optimize (maximize) the agro-economic water productivity (AEWP) of the various crops and, subsequently, the net economic benefit (NEB), using crop areas as decision variables, while respecting various crop cultivation constraints. Adhering to political food security recommendations for the country, three variants of cereal cultivation area constraints are investigated. The results indicate considerably-augmented AEWPs, resulting in a future increase of the annual NEB of ~16% to 37.4 Million USD for the 65%-cereal acreage variant, while, at the same time, the irrigation water required is reduced by ~38%. This NEB-rise is achieved by augmenting the total future crop area in the ZRB by about 47%—indicating some deficit irrigation—wherefore most of this extension will be cultivated by the high AEWP-yielding crops wheat and barley, at the expense of a tremendous reduction of alfalfa acreage. Though presently making up only small base acreages, depending on the future period/RCP, tomatoand, less so, potatoand sugar beet-cultivation areas will also be increased significantly.


Introduction
Food and water security will pose a great challenge in the near future due to rapid growth of population and often unsustainable water usage.The renewable water resources per capita in the Middle East and North Africa (MENA), as the most water-scarce regions of the world, are expected to decline from 750 to 500 m 3 by 2025, while the water withdrawals will increase by up to 50% [1].
Natural and anthropogenically induced climate change will act as an additional external driver threatening the future food security by exacerbating the water shortage, and, concomitantly, the decrease of crop production, as temperatures and irrigation water requirements increase [2].All of this holds particularly for the Middle East, including Iran, where groundwater reserves diminish at an alarming rate [3].
Improving the management of agricultural water use is of utmost importance, as irrigation water uses account for 70% of the global freshwater withdrawals, particularly due to the fact that the irrigated areas have dramatically increased in the 20th century, providing now about 40% of the world's food [4,5].In addition, agriculture has also an important role in the economy, in terms of the Global Gross Domestic Product (GDP), especially in developing countries, although its share has been decreasing over the last twenty years [6].Thus for Iran the GDP contribution of agriculture decreased from 23 to 9% [7], although the irrigated lands increased by 17% between 2003 and 2008 [8], confirming the low economic productivity of agriculture.Moreover, about 90% of the food demands are derived from agriculture supplies in Iran, but with a cost of exploitation of 92% of the available freshwater resources [5], which indicates that the economic agricultural return on water use is outstandingly low in Iran.
The agricultural production will also need to be increased globally by 70% up to year 2050, due to a 40%-projected population increase [9].The situation is even worse for developing countries where the food production should be doubled by that time [10].FAO predicts that only 10% of the global production growth (21% in developing countries) can be achieved by agricultural land expansions with the remainder coming from crop yield enhancements [11].
For a long-term sustainable water resources management for agriculture it is important to quantify and evaluate the possible impacts of climate change scenarios on the future water availability and crop production potentials.Previous publications evaluated the impacts of climate change on the water resources using hydrologic simulation models that are based on GCM-predictions [12,13].These and numerous other studies indicate that climate change will have undeniable impacts on the hydrology, namely, streamflow changes in a basin, which ultimately affects the water availability there.
For example, the impacts of climate change on the crop yield, food security, and crop water demands in sub-Saharan Africa and the North China Plain are investigated by Chijioke et al. [14] and Mo et al. [2], respectively, using different crop prediction and simulation models.These studies show that climate change may have either beneficial or harmful effects on the crop extent and productivity in irrigated or rain-fed agricultural lands.Other recent studies focus on the analyses of the impacts of a changing climate and agricultural demands on the water management and crop production [15,16] and indicate that climate change will lead to hydrologic changes and thus alter crop yields and crop water productivity (CWP), so that some adaptation strategies are required.
Over recent years, many publications on optimizing crop pattern and water allocation to maximize crop productions and economic benefits and to enhance the agricultural water management have appeared.A multi-crop planning (MCP) optimization model that is based on a nonlinear programming (NLP) algorithm was utilized for cropping pattern and water allocation by Bou-Fakhreddine et al. [17] to maximize the net financial return.Firstly, two linear formulations and a relaxed version were established from the NLP and then the MCP problem is solved by implementing two meta-heuristic algorithms, Simulated Annealing (SA) and Particle Swarm Optimization (PSO).Fazlali and Shourian [18] optimized water allocation by considering optimum cropping pattern for the Arayez plain in Iran, using the Shuffled Frog Leaping Algorithm coupled with MODSIM [19] and employing irrigations depths and cultivation areas as decision variables.However, the authors did not consider the CWP index, impacts of climate change, and management of the conjunctive water uses.In fact, few of these issues have been addressed by Fereidoon and Koch [20] who employed a MODSIM-LINGO-PSO algorithm to maximize the economic benefits of Karkheh Dam in Iran, in terms of water allocation for agriculture, under the impact of future climate change.The authors separated the optimization into a three-stage procedure, wherefore first the MODSIM allocates the available water of Karkheh Dam, with its inflow being simulated by the SWAT-hydrological model.Then, a linear optimization to maximize the crop yields in response to different assumed levels of available water is carried out, and, finally, a PSO algorithm is used to maximize the economic return.
The purpose of the current research is to jointly optimize the crop pattern and irrigation planning under climate change-and cropping pattern scenarios and so to maximize the net economic agricultural benefits and increase the crop production, altogether considering the constraints of limited available water resources.This objective is achieved by using an integrated hydro-economic model that consists of a combination of the CSPSO (Constrained Stretched Particle Swarm Optimization) method and MODSIM water management and planning model as a simulation-optimization approach.
The research area is the Zarrine River Basin (ZRB) belonging to the basin of Lake Urmia (LU), which has been shrinking tremendously over the recent decades.The impacts of climate change scenarios on the water resources and the crop production will be evaluated while considering the crop pattern scenarios for the irrigated croplands using the available water supply sources, namely, the Boukan Dam as the most important water management infrastructure of the region, as well as inter-basin discharges from river reaches and groundwater shallow aquifers.
To simulate the basin's water resources, i.e., the inflow of the Boukan Reservoir, interbasin flows, groundwater recharges, and other hydrologic variables of the ZRB in response to the changing climate, future (up to year 2098) downscaled climate predictors (min.and max.temperatures, precipitation) are taken from the recent study of Emami and Koch [13], and entered into the river basin hydrologic model, SWAT, which is firstly calibrated and validated for the discharge and then for the crop yields by adjusting the crop parameters and crop water requirements.Next, a water planning and management simulation model is prepared, using MODSIM for managing the conjunctive agricultural water uses (Q) in the river basin, by respecting the water-distribution priority constraints imposed for the basin as a whole.The model is then adjusted to allocate the available water to the major crop arable areas of the ZRB based on the crop water productivity (CWP) and the net economic benefit (NEB) of the crop production, both of which are combined in an index of agro-economic water productivity (AEWP).Finally, to enhance the crop production and the economic efficiency of the MODSIM-recommended water management policy, a Constrained Stretched Particle Swarm Optimization (CSPSO) algorithm is developed and fully coupled with the MODSIM in order to maximize the total NEB of the ZRB as the objective function and, consequently, the crop production, based on the AEWP-and Q-combination.The optimal arable crop areas and corresponding irrigation schedules are determined while using this CSPSO-MODSIM model under the constraints of arable areas for the irrigation plots and three different cereal crop pattern scenarios when considering three impact scenarios of climate change (RCP45, 60, and 85) for three future periods (near, middle, and far).

The Zarrine River Basin
The Zarrine River is the main inflow source of the LU, the largest inland wetland of Iran which used to be the largest lake in the Middle East before it dwindled significantly in recent decades, with detrimental effects on the surrounding ecosystems of the lake.The Zarrine River Basin (ZRB) is located in northwestern Iran, south of LU between 45 • 46 E to 47 • 23 W longitude and 35 • 41 S to 37 • 44 N latitude (see Figure 1).The total length of the main channel is about 300 km and most of its course stretches through a mountainous area.The basin covers an area of about 12,000 km 2 , including parts of Kurdestan and the West and East Azarbaijan provinces, wherefore its larger portion is mountainous with an elevation of up to 3297 m and the smaller one is rather plain with an elevation going down to 1264 m.The big cities of the basin are namely Miandoab, Shahindej, Tekab, and Saghez.
The climate of the region varies from semi-wet cold or wet-cold in the mountain areas to semi-dry in the vicinity of LU.The average annual temperature varies between 8 and 12 • C while the annual precipitation (rain and snow) in the basin varies between 200 mm/year in the lower catchment area and 800 mm/year in the mountains.The maximum snowfall is recorded mostly in the south and west of the basin with snow heights varying from 5 to 63 mm/year.
The Boukan Dam/Reservoir is the largest operating dam of the ZRB, with a gross storage capacity of 760 Million m 3 (MCM) and a live storage capacity of 654 MCM.51% of its water is used for agricultural irrigation and the remainder for the supply of drinking and industrial water and environmental rights (826 MCM/year totally).
The agricultural areas within the basin cover a total area of 74,318 ha, all irrigated by both groundwater and surface water resources, including water from the Boukan Reservoir, as the crop-growing season there is mostly during the dry months between spring and autumn.
The current applied irrigation efficiency is about 38% for the areas that are irrigated by the surface water from the dam and the river, and 50% for the areas using groundwater resources; all numbers that are lower than the averages of most developing countries (45%) and developed countries (60%) [21], indicating a non-efficient use of surface water due to outdated irrigation methods and systems with a large loss of water through evaporation and seepage.However, advanced irrigation technologies that enhance the irrigation efficiency, can have negative effect on river flow, as a higher efficiency rarely decreases water consumption, which are usually associated with a reduction in recoverable return flows and an increase in crop yields and crop transpiration [22,23].
It should also be noted that the area of irrigated land has been increased by 36% from 1976 to 2013 [24], and this in spite of a catastrophic 88% decrease of the LU surface and ensuing environmental and ecological crises.
The irrigated croplands of the ZR Basin are demonstrated in Figure 2, including the current irrigated croplands of the ZRB (green) and the future agricultural development plan of ZRB, namely Rahimkhan (RK) Plain (light green) located in the downstream of the Boukan Dam.portion is mountainous with an elevation of up to 3297 m and the smaller one is rather plain with an elevation going down to 1264 m.The big cities of the basin are namely Miandoab, Shahindej, Tekab, and Saghez.The climate of the region varies from semi-wet cold or wet-cold in the mountain areas to semi-dry in the vicinity of LU.The average annual temperature varies between 8 and 12 °C while the annual precipitation (rain and snow) in the basin varies between 200 mm/year in the lower catchment area and 800 mm/year in the mountains.The maximum snowfall is recorded mostly in the south and west of the basin with snow heights varying from 5 to 63 mm/year.
The Boukan Dam/Reservoir is the largest operating dam of the ZRB, with a gross storage capacity of 760 Million m 3 (MCM) and a live storage capacity of 654 MCM.51% of its water is used for agricultural irrigation and the remainder for the supply of drinking and industrial water and environmental rights (826 MCM/year totally).
The agricultural areas within the basin cover a total area of 74,318 ha, all irrigated by both groundwater and surface water resources, including water from the Boukan Reservoir, as the crop-growing season there is mostly during the dry months between spring and autumn.
The current applied irrigation efficiency is about 38% for the areas that are irrigated by the surface water from the dam and the river, and 50% for the areas using groundwater resources; all numbers that are lower than the averages of most developing countries (45%) and developed countries (60%) [21], indicating a non-efficient use of surface water due to outdated irrigation methods and systems with a large loss of water through evaporation and seepage.However, advanced irrigation technologies that enhance the irrigation efficiency, can have negative effect on river flow, as a higher efficiency rarely decreases water consumption, which are usually associated with a reduction in recoverable return flows and an increase in crop yields and crop transpiration [22,23].
It should also be noted that the area of irrigated land has been increased by 36% from 1976 to 2013 [24], and this in spite of a catastrophic 88% decrease of the LU surface and ensuing environmental and ecological crises.
The irrigated croplands of the ZR Basin are demonstrated in Figure 2, including the current irrigated croplands of the ZRB (green) and the future agricultural development plan of ZRB, namely Rahimkhan (RK) Plain (light green) located in the downstream of the Boukan Dam.The main agricultural crops of the ZRB as well as the RK include alfalfa (ALFA), apple (APPL), barley (BARL), potatoes (POTA), sugar beet (SGBT), tomatoes (TOMA), and wheat (WWHT), and these are the ones that are considered in this study.

Data
The data needed for this study in the ZRB is gathered from different available sources.Most of the data is required for the set-up of the SWAT-hydrologic model, namely, various geospatial maps and hydro-climate time series.
The Digital Elevation Map (DEM) with a spatial resolution of 85 m was produced by the Iranian surveying organization.The land-use classification map of the basin, demonstrating the situation in year 2007, was obtained from the Agricultural Statistics and the Information Center of the Ministry of Agriculture [25] and it has a resolution of 1000 m and distinguishes 10 land use classes.The soil map of the watershed was extracted from the Food and Agriculture Organization (FAO) digital soil global map, with a spatial resolution of 10 km of and eight types of soils within two layers.
The daily climate input data includes maximum and minimum temperatures and precipitation over the period of 1987 to 2015 and was obtained from the Iranian Meteorological Organization (IRIMO) for six synoptic stations that were located in or close to the ZRB (see Figure 1).Missing data in the records were filled in using the inverse distance weighting (IDW) interpolation method.As the daily data for other climate variables, including solar radiation, wind, and relative humidity, were unavailable, they were generated using the weather generator (WGEN) module of SWAT model based on monthly averages of the synoptic stations of Iran.
Daily streamflow data for six gauging stations of the Zarrine River (Figure 1) were obtained from the Iran Ministry of Energy for the period 1987 to 2012.
The crop and irrigation data, namely, crop irrigation sources, planting, irrigation and harvesting dates, or water demands were taken from Ahmadzadeh et al. and MOE [22,26].The observed crop yields are gathered from MOA [25] and the additional economic crop data were gathered from SCI and MOA [27,28].The main agricultural crops of the ZRB as well as the RK include alfalfa (ALFA), apple (APPL), barley (BARL), potatoes (POTA), sugar beet (SGBT), tomatoes (TOMA), and wheat (WWHT), and these are the ones that are considered in this study.

Data
The data needed for this study in the ZRB is gathered from different available sources.Most of the data is required for the set-up of the SWAT-hydrologic model, namely, various geospatial maps and hydro-climate time series.
The Digital Elevation Map (DEM) with a spatial resolution of 85 m was produced by the Iranian surveying organization.The land-use classification map of the basin, demonstrating the situation in year 2007, was obtained from the Agricultural Statistics and the Information Center of the Ministry of Agriculture [25] and it has a resolution of 1000 m and distinguishes 10 land use classes.The soil map of the watershed was extracted from the Food and Agriculture Organization (FAO) digital soil global map, with a spatial resolution of 10 km of and eight types of soils within two layers.
The daily climate input data includes maximum and minimum temperatures and precipitation over the period of 1987 to 2015 and was obtained from the Iranian Meteorological Organization (IRIMO) for six synoptic stations that were located in or close to the ZRB (see Figure 1).Missing data in the records were filled in using the inverse distance weighting (IDW) interpolation method.As the daily data for other climate variables, including solar radiation, wind, and relative humidity, were unavailable, they were generated using the weather generator (WGEN) module of SWAT model based on monthly averages of the synoptic stations of Iran.
Daily streamflow data for six gauging stations of the Zarrine River (Figure 1) were obtained from the Iran Ministry of Energy for the period 1987 to 2012.
The crop and irrigation data, namely, crop irrigation sources, planting, irrigation and harvesting dates, or water demands were taken from Ahmadzadeh et al. and MOE [22,26].The observed crop yields are gathered from MOA [25] and the additional economic crop data were gathered from SCI and MOA [27,28].

Basic Concepts of an Optimal Hydro-Economic Model
Generally, the greater the water-use productivity and economic efficiency, the lower are the conflicts over scarce water resources and the additional financial and environmental burdens in an agriculturally exploited basin, like the ZRB.Enters the fundamental concept of hydro-economics, which stipulates that water demands can be represented as value-sensitive water demand functions, so that water-uses at different locations and times have varying economic benefits [29].The next step is then to set up a hydro-economic model, which is a solution-oriented model for investigating the water management tradeoffs and improving the economic efficiency of water allocation by incorporating the economic value of agricultural water in the heart of the water management model.This model represents a spatially distributed water resources system, infrastructure, management options, and economic values comprehensively [29].In the third and final step, such a hydro-economic model may be applied to simulate agricultural crop pattern strategies, with the goal to find that optimal multi-crop pattern that somehow maximizes the economic crop profits, while adhering to the various constraints of the limited water resources, hydrology, and various environmental regulations.Mathematically, this amounts to the set-up of a classical constrained (nonlinear) optimization problem, wherefore (1) the forward problem, i.e., the objective or cost function, is computed by the hydro-economic model, simulating the hydrological constraints by classical hydrological models, and (2) the minimization of that objective function is done by some kind of an optimization routine (e.g., [20]).
In this section, an innovative hydro-economic model for the ZRB is developed using such a simulation-based optimization approach to coordinate multiple factors, including water allocation, crop production pattern, and economic gains.More specifically, the optimization problem is defined as a constrained optimization (CO) problem which searches for the optimal allocation of irrigated crop pattern under the constraints of the limited water resources and other demands that should be satisfied.Eventually, the objective of the optimization search algorithm is to maximize the agro-economic productivity, i.e., the economic net benefit of a crop per unit water use, given that the latter is limited in the study region and it may be even more so in the future under the impacts of imminent climate change there [13].
The decision variables of the optimization are the cultivated areas of major crops for a particular, politically given combination of required crop distribution and agricultural demand regions.The constraints of the optimization algorithm in this study are defined based on the allowable range of arable areas and cereal crop pattern limits (for more details see Section 3.8).

Modules of the CSPSO-MODSIM Integrated Hydro-Economic Model
The individual modules (hydrological, water management, agro-economic) of the hydro-economic model are developed in this research while using different simulation models (GCM, QM, SWAT, MODSIM) that were bound together with an optimization method (CSPSO).The flow chart of the connection of the models and processes in the integrated hydro-economic model is presented in Figure 3, from which the main steps are retrieved as follows: -

Predicting Future Climatic Scenarios by Updated Quantile Mapping
Regarding the projections of future climate change scenarios, they are mostly taken from Emami and Koch [13] who used several climate models (GCMs) of the CMIP5 archive [30], as mentioned in the fifth IPCC report [31] to select the most suitable one based on a climate model's skill to simulate the past climate in terms of minimum and maximum temperatures and precipitation.For assessing the impacts of climate change on the regional scale, the climate predictors of the selected GCM models were then downscaled using a recently-coming-to-the-fore statistical downscaling method, namely, QM (Quantile Mapping) [32,33], which proved to have better prediction performances than the more commonly used classical statistical downscaling model (SDSM).We forego a detailed description of the QM-downscaling method and refer the reader to [13].
In the first step of the QM method, the monthly biases of the future GCM-simulated climate variables between year 2020 to 2098 are removed while using a trend-preserving bias correction, namely, the ISI-MIP approach [32].In this approach, the GCM-simulated temperatures (min.and max.) are corrected applying an additive correction factor CTj and for the precipitation and a multiplicative correction factor CPj to each month of a year (j = 1, …, 12) of the GCM-simulated precipitation.In the second step a new, updated quantile mapping method [33] is used to correct the daily biases of the temperatures and precipitation in each month (e.g., all Januaries) using EDCDFm and CNCDFm CDF-(cumulative distribution function) matching methods, respectively.Emami and Koch [13] proved that these QM-variants perform better than other statistical downscaling methods in removing biases in the GCM-climate predictors for the study region by delivering much better correlations with the observed predictands (temperatures and precipitation) at the high-resolution, local scale, and of all this without almost no extra computational costs.

Predicting Future Climatic Scenarios by Updated Quantile Mapping
Regarding the projections of future climate change scenarios, they are mostly taken from Emami and Koch [13] who used several climate models (GCMs) of the CMIP5 archive [30], as mentioned in the fifth IPCC report [31] to select the most suitable one based on a climate model's skill to simulate the past climate in terms of minimum and maximum temperatures and precipitation.For assessing the impacts of climate change on the regional scale, the climate predictors of the selected GCM models were then downscaled using a recently-coming-to-the-fore statistical downscaling method, namely, QM (Quantile Mapping) [32,33], which proved to have better prediction performances than the more commonly used classical statistical downscaling model (SDSM).We forego a detailed description of the QM-downscaling method and refer the reader to [13].
In the first step of the QM method, the monthly biases of the future GCM-simulated climate variables between year 2020 to 2098 are removed while using a trend-preserving bias correction, namely, the ISI-MIP approach [32].In this approach, the GCM-simulated temperatures (min.and max.) are corrected applying an additive correction factor CT j and for the precipitation and a multiplicative correction factor CP j to each month of a year (j = 1, . . ., 12) of the GCM-simulated precipitation.In the second step a new, updated quantile mapping method [33] is used to correct the daily biases of the temperatures and precipitation in each month (e.g., all Januaries) using EDCDFm and CNCDFm CDF-(cumulative distribution function) matching methods, respectively.Emami and Koch [13] proved that these QM-variants perform better than other statistical downscaling methods in removing biases in the GCM-climate predictors for the study region by delivering much better correlations with the observed predictands (temperatures and precipitation) at the high-resolution, local scale, and of all this without almost no extra computational costs.

Model Setup and Calibration
The Soil Water Assessment Tool (SWAT) model is a physically-based, river basin-scale, time continuous simulation model that operates on a daily time step.Although this model has originally been developed to mainly simulate the impacts of land management practices in large and complex watersheds [34], it is widely used as a long-term rainfall-runoff model and efficient hydrologic simulator of water quantity and quality, so that it has increasingly being used to investigate climate change impacts on agro-hydrological systems [12,35].
The SWAT model requires quite a wide range of input data, as described in Section 2.2.To represent the large-scale spatial heterogeneity of the study basin more precisely, the SWAT modeled domain, i.e., the major basin, is divided into several sub-basins, which are usually delineated with the help of the topographic DEM using the ARCSWAT extension of ARCMap program.Then, each sub-basin is parameterized using a set of HRUs (Hydrologic Response Units), which are based on a unique combination of soil and land cover and management.For the SWAT-model of the ZRB, 11 sub-basins with a total of 908 HRUs have been defined.
After the parameterization of the SWAT-model's input data entries is done by using the stochastic sequential uncertainty fitting version 2 (SUFI-2) optimization algorithm, embedded in the SWAT-CUP decision-making framework [36].Various kinds of objective functions as measures for the goodness of the fit of the modeled to the observed streamflow are available in the SUFI-2 algorithm, wherefore Krause et al. [37] indicates that for a reliable calibration and validation of the model, a combination of different efficiency criteria, such as the coefficient of determination R 2 , the Nash-Sutcliff efficiency coefficient NSE, and bR 2 (b is the slope of the regression line between observed and simulated streamflow), should be considered.
In the present application, the calibrated input parameters of the model have been taken from Emami and Koch [13], where further details of the calibration/validation as well as of the model setup are presented.Basically, the optimal range of model input parameters was determined hierarchically sub-basin-wise, from the utmost upstream sub-basin (11) outlet down to the main outlet of the basin.
Based on the SUFI-sensitivity analysis, 24 model parameters were shown to be sensitive parameters for affecting the stream discharge, out of which the SCS curve number (CN2), the groundwater delay time (GW_DELAY), and the moist bulk density of the soil (SOL_BD) turned out to be the three most sensitive variables.With these optimized parameters, good fits of the modeled to the observed discharges were obtained at the six streamflow stations (sub-basin outlets) for the calibration (1998-2012) and the validation (1991-1997) periods, with average R 2 > 0.7, NSE > 0.6 and bR 2 > 0.5, which, according to the classification that was proposed by Moriasi et al. [38], is considered to be satisfactory.As the SUFI-computed uncertainty of the calibrated model that is quantified by the Pand the R-factor (see [13,39]) has average values of R > 0.75 and P close to 1, there is enough confidence in the calibrated SWAT-model for the ZRB.

Predicting the Impacts of Future Climate Change on the Hydrologic Cycle
It is important to evaluate the hydrologic responses to future changes of climate for improving adaptive water management, as the variability of precipitation and temperature in terms of trends and extremes will eventually increase the likelihood of severe and irreversible negative impacts on the ecosystem, including lakes and rivers.To that avail, the future downscaled climatic scenarios i.e., the QM-bias corrected predictions of the minimum and maximum temperatures and precipitation are employed as weather input drivers of the calibrated basin-wide hydrologic simulation model, SWAT.As a result, the future hydrologic cycle and the available water resources of the ZRB under the climate change can be assessed, especially the input of the dam, the Inter-basin discharge of the river reaches, the groundwater recharges and the withdrawals of the shallow aquifers.

Crop Yield Simulation and Calibration of the Potential Crop Yield
The SWAT model is also capable of simulating crop productions and yields efficiently, as has been shown in many publications (e.g., [40,41]).To do this in the ZRB, the current management operations of the various crops there are specified initially in the SWAT model, together with the corresponding planting and harvesting dates and the irrigation sources, based on information that was given by Ahmadzadeh et al. [22].The crop yields for seven major crops in the ZRB are then calibrated by adjusting a set of effective parameters in the model, until the averages of the simulated crop yields of the basin match those of the observed ones (gathered from MOA [25]) in a reasonable manner.These simulated crop yields are then extracted from the SWAT file output.hruto represent the potential crop yields of the ZRB.The reservoir characteristics, irrigation losses and essential demands are considered through the water use management (.wus) and reservoir (.res) SWAT-data files.

Agro-Economic Water Productivity and Net Economic Benefit
For a better management of the future water resources in a water-scarce region, such as the ZRB, it is necessary to make the water supply-and/or the irrigation system as efficient as possible in economic and water resources perspectives.Because of the competition of the different stakeholders for the scarce freshwater resources in the region, not only for agriculture, a paradigmatic policy shift is required from (a) maximizing productivity per unit of area to (b) maximizing productivity or economic value per unit of consumed water [42], as both the irrigated agriculture's land base and the water supplies are continuously being depleted and reallocated, in order to produce even more agricultural crops.To achieve this policy shift, the net benefits of the water used, should be enhanced by prioritizing more economic efficient crops i.e., higher productivity per unit of water, while the water use should be constrained to the allocated water resources for agricultural water-use.
The crop-water productivity CWP c (kg/m 3 ) is defined as the ratio of amount of crop yield that is produced Y c (kg/ha) to the amount of water applied per unit crop area Q c (m 3 /ha) during the crop's production [43].The next step is then to define the agro-economic water productivity AEWP c (USD/m 3 ) as the ratio of net total economic value of crop NEB c (USD) to the total amount of water Q t c (m 3 ) supplied under the priority constraints that were provided by MODSIM (see following subsection) [44]: where Price c (USD/kg) is the selling price of the crop, Cost c (USD/ha) is the total production costs of the crop, A c (ha) is the crop cultivation area, Irr c = Q t c /A c (m) is the irrigation water height, and the other variables are as defined above.
Based on the definitions above, optimizing the water resources allocation can be an important measure for increasing agro-economic water productivity AEWP c in areas where the water is scarce like ZRB.In this study, a water management model for deficit irrigation is developed to determine the allowable agricultural water, which at the first stage the potable, industrial, and environmental water demands (essential demands) are supplied according to their higher priorities as a constraint and then the remaining available water will be allocated to the agricultural water-uses, i.e., different crop demands based on their AEWPs.Finally, for boosting the net economic benefits NEB c as the ultimate objective of this study, the combination of agro-economic water productivity AEWP c and total allocated water Q t c is maximized through an optimization algorithm.As in Equation (1), the Price c and Cost c are constant for each crop, the total crop economic benefit will also be increased.

MODSIM Water Resources Management and Planning Module
The MODSIM simulation model is a generalized river basin network model for developing basin-wide strategies of short-term water management, long-term operational planning, drought contingency planning, water rights analysis, and conflict resolution between different water users [19,45].This model has enjoyed widespread application across the world to simulate operations, as mentioned [20,[46][47][48].
The core idea behind MODSIM is to represent a complex river basin system by a flow network consisting of coupled sequences of nodes and links, with the former symbolizing storage components, such as reservoirs and aquifers, points of inflow, demands, diversions, and river confluences, and the latter representing river reaches, pipelines, canals, and stream-aquifer interconnections defining stream depletions from pumping and return flows from seepage and other water applications.Further network elements of the model consist of unregulated inflows, reservoir operating targets, consumptive and instream flow demands, evaporation and channel losses, reservoir storage rights and exchanges, and stream-aquifer modeling components.In addition, various surface and ground water resources with their inter-relationships can be represented by highly nonlinear, non-convex, or discontinuous equations [49].Details on how each of the components is modeled and calculated in the model can be found in Fredericks et al. [19].
The model sequentially solves a linear optimization problem within the confines of mass balance throughout the network over the planning period by means of a highly efficient network flow optimization (NFO) algorithm solved with the Lagrangian relaxation algorithm RELAX-IV [50].More specifically, the following constrained flow optimization problem is solved for each time interval (t = 1, . . ., T) over the planning horizon: where A is the set of all links in the network; N is the set of all nodes; q k is the integer valued flow rate in link k; c k are cost weighting factors, i.e., the water right priorities per unit flow rate in link k; b it is the (positive) gain or (negative) loss at node i at time t; I i the set of all links terminating at node i (inflow links); O i the set of all links originating at node i (outflow links); and, l k and u k are lower and upper bounds, respectively.The cost factor c k for accounting for active storage and demand links priorities are generally calculated while using the following formula: where PR k is the integer priority ranking that ranges between 1 and 5000 for the reservoirs or the demands, where the negative sign states that high-rank nodes (1, 2, etc.) are given more weights in the minimization of the cost function (Equation ( 2)).The cost factor may also include economic factors that are defined in this study for the crop demand links as the cost of the crop water supply (cc c ) based on the agro-economic water productivity index AEWP c (Equation ( 1)), i.e., cc c = −AEWP c , which means that more water should be allocated to a crop that provides more economic benefits than others for the same amount of water supplied.The schematic network of the MODSIM ZRB model with the various demand nodes are represented in Figure 4.These demand nodes can be categorized into four groups: (1) the network of the dam croplands (ZRB), suffixed _ZRB, (2) the future development croplands of the RK, suffixed _RK, (3) the crop demands supplied from the reach, suffixed _rch, and (4) the demands supplied from the groundwater aquifers, wherefore two aquifer storages are defined cumulatively for the upstream-, suffixed _gwup, and downstream, suffixed gwdown areas of the Boukan Dam.The conjunctive water uses, i.e., surface-and groundwater irrigation are linked by connecting the river network nodes and the shallow aquifer storages with the corresponding demand nodes.Other non-agricultural water demands include the potable demand of Boukan and Tabriz (about 60% of its potable demand) cities, the industrial demands of the dam, the Legzi water transfer, the other agricultural water rights, the environmental rights of the Boukan Dam, and the LU environmental water demand.The hydrologic inputs for the MODSIM-model are captured from the flow results of the SWAT model, including the inflow of the dam, the river discharges and the inter-basin flows (the inter-basins 7-3 and 2-1 describe the generated water in sub-basins 7 to 3 and sub-basins 2 to 1, respectively), the storage of the shallow aquifers, the recharges (GW_Rchg), the contributing flows to the surface water and the drainage (GWQ + Drainage).The different types of the water demands, water transfers (such as the transfer to the RK Plain from the Simineh River), the return flows from the agricultural water use (about 10% to the same source, groundwater aquifers or the Zarrine River) are all initially defined based on values that are given by MOE [25] and the irrigation and water requirements of the various crops and irrigation losses are adjusted based from provisional values of Ahmadzadeh et al. [22].The MODSIM ZRB simulation model is customized on the MODSIM 8.5 platform using the custom coding module in VB.NET routine, with further details being provided in Section 3.8.

Constrained Stretched Particle Swarm Optimization (CSPSO) Search Algorithm
The Particle Swarm Optimization (PSO) method that was proposed first by Kennedy and Eberhart [51] is a stochastic evolutionary social behavior-based optimization algorithm for solving nonlinear global optimization problems in an efficient way.The main idea behind the development of the PSO is social sharing of information among individuals of a population in nature (the flock or swarm), in order to provide an evolutionary advantage for all individuals to move towards some optimum [52][53][54].
A PSO model consists of a number of N particles (the swarm) moving around in the D-dimensional search space, with each particle representing a possible solution to a numerical problem.In this D-dimensional search space, the actual location of the i-th particle (i = 1, …, N) can be represented by a D-dimensional vector of position Xi = (xi1, xi2, …, xiD) and the position change (velocity) by another D-dimensional vector Vi = (vi1, vi2, …, viD).The best candidate solution that is the best previously visited position of the particles of the swarm is denoted as Pi = (pi1, pi2, …, piD).
Assuming that the g-th particle is the best and denoting the iteration by the superscript n, the swarm is manipulated according to the following two equations [55,56]: The hydrologic inputs for the MODSIM-model are captured from the flow results of the SWAT model, including the inflow of the dam, the river discharges and the inter-basin flows (the inter-basins 7-3 and 2-1 describe the generated water in sub-basins 7 to 3 and sub-basins 2 to 1, respectively), the storage of the shallow aquifers, the recharges (GW_Rchg), the contributing flows to the surface water and the drainage (GWQ + Drainage).The different types of the water demands, water transfers (such as the transfer to the RK Plain from the Simineh River), the return flows from the agricultural water use (about 10% to the same source, groundwater aquifers or the Zarrine River) are all initially defined based on values that are given by MOE [25] and the irrigation and water requirements of the various crops and irrigation losses are adjusted based from provisional values of Ahmadzadeh et al. [22].The MODSIM ZRB simulation model is customized on the MODSIM 8.5 platform using the custom coding module in VB.NET routine, with further details being provided in Section 3.8.

Constrained Stretched Particle Swarm Optimization (CSPSO) Search Algorithm
The Particle Swarm Optimization (PSO) method that was proposed first by Kennedy and Eberhart [51] is a stochastic evolutionary social behavior-based optimization algorithm for solving nonlinear global optimization problems in an efficient way.The main idea behind the development of the PSO is social sharing of information among individuals of a population in nature (the flock or swarm), in order to provide an evolutionary advantage for all individuals to move towards some optimum [52][53][54].
A PSO model consists of a number of N particles (the swarm) moving around in the D-dimensional search space, with each particle representing a possible solution to a numerical problem.In this D-dimensional search space, the actual location of the i-th particle (i = 1, . . ., N) can be represented by a D-dimensional vector of position X i = (x i1 , x i2 , . . ., x iD ) and the position change (velocity) by another D-dimensional vector V i = (v i1 , v i2 , . . ., v iD ).The best candidate solution that is the best previously visited position of the particles of the swarm is denoted as P i = (p i1 , p i2 , . . ., p iD ).
Assuming that the g-th particle is the best and denoting the iteration by the superscript n, the swarm is manipulated according to the following two equations [55,56]: ) where d = 1, . . ., D, with D the number of decision variables; i = 1, . . ., N, with N the size of the swarm; r 1 , r 2 are uniformly distributed random numbers in [0, 1]; n = 1, 2, . . ., the iteration number; c 1 , c 2 acceleration coefficients that are, respectively, the cognitive and social components of the particle velocity, representing the impact of self-knowledge and the collective effect of the population; χ a constriction parameter that is employed, alternatively to the inertia parameter ω to limit the velocity to the range [v min , v max ].By incorporating a recently proposed technique called Function Stretching into the classical PSO, Parsopoulos and Vrahatis [53] arrived at the SPSO method, which has the capability to alleviate the attractions of local minima of the objective function and so to rise the success rates for finding a truly global solution of the problem.
The basic idea of SPSO is to use a two-stage transformation of the original objective function f (x), which can be applied immediately after a local minimum x of the function f (x) has been detected, defined as follows: where γ 1 , γ 2 , and µ are arbitrary user-defined constants.The first transformation stage, G(x) elevates the function f (x), eliminating so most of the local minima, whereas the second stage, H(x) stretches the neighborhood of x upwards, as it assigns higher function values to those points.The location of the global minimum is left unchanged, as both stages do not alter the local minima located below x.
Further details of the implementation of the of CSPSO-methodology and its extension for solving the constrained optimization (CO) problem of the maximization of the net economic benefit of the total cultivated crops in the area under the constraints of limited water resources, water right priorities, and crop area limitations are present in the following sub-section.

Formulation of the Constrained Optimization Problem
The ultimate objective of the hydro-economic optimization model is to maximize the economic productivity of the sum of all the crops in an irrigation plot under the constraints of limited water resources and crop areas available.By summing up Equation (1) for all cmax crops, the objective function is: (with the notations as given for Equation (1), so that the final optimization problem can be stated, as follows: wherefore, for application of classical optimization (minimization) routines, like CSPSO, the objective function Z is replaced by its negative.
This, yet unconstrained optimization formulation will be extended to a constrained one by adding appropriate constraints, as discussed below: First of all, the actual crop yield Y c entering Equation ( 12) is dependent on the amount of irrigated water available, which, in turn, depends on the actual water allocation (simulated by MODSIM).Thus, if the irrigation water amount applied, Q c , is less than the crop water requirement, WR c , the actual crop yield, Y c , will not reach its potential (maximum) crop yield Ymax c .This is epitomized by the following FAO-equation for the water production function [57]: where Ymax c is the potential crop yield, estimated by SWAT as mentioned earlier, WR c is the crop water requirement, Q c is the MODSIM-simulated allocated irrigation water, Ky c is a FAO-based yield response factor describing the effect of a reduction of irrigation water on the crops yield losses, with values being gathered from Steduto et al. [58].The crop water requirements WR c were determined using NetWat software from CropWat application series of the Iranian Water Directive (IWD) that was developed by the Iran government for estimating the effects of future climate change [59].Equation (14) shows clearly that unless the crop water requirement WR c is fully met by irrigation water Q c , the actual crop yield Y c will be less than its potential crop yield Ymax c .The second kind of constraints for the optimization problem ( 14) model arises then firstly from the fact that the sum of all croplands cannot exceed the total arable area of the irrigation plot, i.e., ∑ where At ZRB = 684.2km 2 , At RK = 125 km 2 , and At ups = 250 km 2 are the future total arable areas supplied by the Boukan Dam irrigation network (ZRB), of the RK plain and of the agricultural areas upstream of the dam, respectively, according to the SWAT-land use map and based on information of MOE [25].
The sums in Equations ( 15)-( 17) run over seven different areas, as this is the number of the major crops cultivated in the region, which are, in alphabetic order, alfalfa, apple, barley, potatoes, sugar beets, tomatoes, and wheat, although the apple croplands are assumed to be constant to current cultivated areas because it is irreplaceable as orchard and have very low economic value.Further constraints that are varied later in the modelled scenarios pay attention to the fact that a high-benefit (low water costs) crops cannot be solely cultivated over the whole ZRB, but there are constraints on areas attributed to the various crops, not to the least to satisfy the population's food demands.More specifically, for each arable crop an allowable range of area size A c is defined, wherefore the maximum values are taken from MOE [25] and for the minimum areal sizes two alternatives are investigated.In both of them, denoted as Smin 1 and Smin 2 in Table 1, the minimum areas devoted to the two cereals (barley, wheat) of the Boukan Dam network plot (BARL_ZRB and WWHT_ZRB) are set to the current arable area, whereas for the other crop/area demand nodes, the minimum is assumed to have (1) no limitation (=0) for Smin 1, and (2) 60% of the maximum arable area for Smin 2 .Based on these two numbers, the corresponding areal sizes for the different crops have been computed for the different irrigation plots and listed in Table 1.
Table 1.Area constraints (km 2 ) for the two minimal-area scenarios (Smin 1 and Smin 2 ) devoted to the different crops for the three irrigation demand areas (ZRB, RK, and upstream areas of the Boukan Dam).Finally, as the production of cereal crops, i.e., barley (BARL) and wheat (WWHT) has strategic importance in the ZRB as well as for Iran overall, three cereal crop pattern scenarios are considered further in this study, namely, that the sum of the cropland areas devoted to barley (A BARL ) and wheat (A WWHT ) will be a portion X of the maximum-minimum range above the minimum areas, i.e.,

Crop/Area Name
where X is the limiting production factor of the cereal crop pattern scenario and is set to, respectively, X = 0.35, 0.5, and 0.65 for the three scenarios investigated, and A min and A max are the sum of the area limits of the different demand plots for the two crops (barley and wheat) (Table 1).

Integration of MODSIM and CSPSO Algorithms
To solve the optimization problem, Equation (13), which is subject to the four constraints, Equations ( 15)- (18), by means of the CSPSO method, a penalty function method is used (e.g., Parsopoulos and Vrahatis [54]).The latter allows for converting a constrained optimization problem to an unconstrained optimization one by adding the constraints as a weighted penalty to the objective function (Equation ( 13)) in the form: where nc is the number of constraints, PF i is the penalty factor of i-th constraint (Equations ( 15)-( 19)), which takes the binary values 0 or 1, depending on whether the constraint is satisfied or violated, respectively, and h is the static penalty weighting factor value, which is found to be 1,000,000 by trial and error based on a convergence analysis as well as PSO-literature recommended values.Other parameters to be specified in the various CSPSO-equations of the previous sub-section are taken in agreement with recommended literature values as (e.g., [54,55]): N = 40 (swarm size), n = 500 (maximum iteration number); c 1 = 1.2 and c 2 = 0.8 (acceleration coefficients), γ 1 = 5000, γ 2 = 0.5, µ = 10E-10 (three constraints of function stretching) and χ = 1 (constriction parameter).
Finally, to arrive at the CSPSO-MODSIM integrated hydro-economic model (see Figure 3), the MODSIM water management tool is embedded in the CSPSO-optimization method as an inner layer of the iteration process.Thus, in the first iteration, the search algorithm generates the decision variables of the arable crop areas that should meet the ranges specified in Table 1.In the next step, the water demands and their priorities are designated to the MODSIM model based on the initial irrigation depth (estimated with SWAT) and the AEWPs (Equation ( 1)).Then, the optimal irrigation depths of the crops are estimated in MODSIM using Equations ( 2)-( 4), with flow inputs from SWAT.The actual crop yields are predicted using Equation ( 14) and returned to the CSPSO model, together with the optimal irrigation depths, to calculate the fitness/penalty function, Equation (19).This procedure, coded in the MATLAB© environment, is repeated for each iteration of the CSPSO, using the swarm-intelligence, until the penalized objective function converges to the maximum net economic benefit of the total crop production in the three irrigation plots.Usually, after 400-500 iterations, no noticeable further improvement in the objective function was obtained.

Historical and Future Climate Projections
The min.and max.temperatures-as well as the precipitation-predictors of the CGCM3-and the CESM-CAM5 GCM from the CMIP5-GCMs archive were found, respectively, by virtue of the skill score multi-criteria method, to be the most suitable climate predictors and used then for further QM-downscaling.This task was carried out for three RCP-scenarios (RCP45, RCP60, and RCP85).
The QM-downscaling model was calibrated and validated for each month of the year (e.g., January) for the time periods 1987-1998 and 1999-2005, respectively.The validation of the model was considered satisfactory, as the CDFs of the bias-corrected min/max temperatures and precipitation fit the CDF of the corresponding observed variables.In addition, the reliability of the QM-downscaled predictors was also evaluated by comparing them with the observed data for the period 2006-2015, using the goodness of fit measures coefficient of determination (R 2 ), standard error (SE), and index of agreement (IA), which are summarized for the three RCPs that were investigated in Table 2 (for further details see [13]). Figure 5 shows the spatial distributions of the observed/QM-downscaled average annual temperatures and precipitation in the ZRB for the historical reference period .
GCM/QM-downscaled climate predictions for three RCP scenarios (RCP45, RCP60, and RCP85) are shown for three future periods, near (2020-2038), middle (2050-2068), and far future (2080-2098), in Figures 6-8, respectively.One may notice from these figures that the trends in all future RCP-scenarios are approximately the same, such that, as compared with the historical references period (see Figure 5), both the temperature and the precipitation are mostly increased.
In particular, for the near future period (Figure 6), RCP45-and RCP60-scenarios turn out to be wetter than RCP85, with RCP45 having at the same time a moderate temperature rise, RCP60 practically no temperature rise, and RCP85 exhibiting a rather high temperature increase.
in Figures 6-8, respectively.One may notice from these figures that the trends in all future RCP-scenarios are approximately the same, such that, as compared with the historical references period (see Figure 5), both the temperature and the precipitation are mostly increased.
In particular, for the near future period (Figure 6), RCP45-and RCP60-scenarios turn out to be wetter than RCP85, with RCP45 having at the same time a moderate temperature rise, RCP60 practically no temperature rise, and RCP85 exhibiting a rather high temperature increase.For the middle future period (Figure 7), and as compared with the near future period (Figure 6), the RCP45-scenario will be drier again, going hand in hand with a small temperature rise, whereas RCP60 will become wetter while temperature will rise moderately.For RCP85, the trend is straightforward, with both a temperature and precipitation increase.

Observed/bias-corrected predictors
Finally, the far-future period (Figure 8) must be highlighted as most critical, as, when compared with the middle-future period (Figure 7), on one hand, the temperature increases another 3% to 14%, and, on the other hand, the precipitation decreases by another 4% to 22%, depending on the RCPs, with RCP85, as the most extreme, at the upper ends of these ranges.For the middle future period (Figure 7), and as compared with the near future period (Figure 6), the RCP45-scenario will be drier again, going hand in hand with a small temperature rise, whereas RCP60 will become wetter while temperature will rise moderately.For RCP85, the trend is straightforward, with both a temperature and precipitation increase.
Finally, the far-future period (Figure 8) must be highlighted as most critical, as, when compared with the middle-future period (Figure 7), on one hand, the temperature increases another 3% to 14%, and, on the other hand, the precipitation decreases by another 4% to 22%, depending on the RCPs, with RCP85, as the most extreme, at the upper ends of these ranges.

Climate Change Impacts on the Hydrology of the ZRB
To evaluate the future climate impacts on the hydrology of the ZRB, the calibrated SWAT model was driven by the QM-downscaled GCM-temperatures and precipitation predictors of the previous section to simulate, among other hydrological components, the inflow to the Boukan as the major irrigation water reservoir in the ZRB.The annually aggregated SWAT-modelled inflow hydrograph of the Boukan Dam over the historical and future predicted periods under the three RCPs, together with their linear trend lines, is presented in Figure 9.
One may notice from the figure that there is a consistent decreasing inflow trend over the whole, i.e., historic and future time periods, though the trend is stronger for the former-where the inflow has decreased by a factor of two from a maximum of ~2000 MCM/year in the 1990s to less than 1000 MCM/year at the present time-than the latter.Also, the negative trend is more pronounced for the more extreme RCP85 climate scenario than for the other two, RCP45 and RCP60.

Climate Change Impacts on the Hydrology of the ZRB
To evaluate the future climate impacts on the hydrology of the ZRB, the calibrated SWAT model was driven by the QM-downscaled GCM-temperatures and precipitation predictors of the previous section to simulate, among other hydrological components, the inflow to the Boukan as the major irrigation water reservoir in the ZRB.The annually aggregated SWAT-modelled inflow hydrograph of the Boukan Dam over the historical and future predicted periods under the three RCPs, together with their linear trend lines, is presented in Figure 9.
One may notice from the figure that there is a consistent decreasing inflow trend over the whole, i.e., historic and future time periods, though the trend is stronger for the former-where the inflow has decreased by a factor of two from a maximum of ~2000 MCM/year in the 1990s to less than 1000 MCM/year at the present time-than the latter.Also, the negative trend is more pronounced for the more extreme RCP85 climate scenario than for the other two, RCP45 and RCP60.
To evaluate the future climate impacts on the hydrology of the ZRB, the calibrated SWAT model was driven by the QM-downscaled GCM-temperatures and precipitation predictors of the previous section to simulate, among other hydrological components, the inflow to the Boukan as the major irrigation water reservoir in the ZRB.The annually aggregated SWAT-modelled inflow hydrograph of the Boukan Dam over the historical and future predicted periods under the three RCPs, together with their linear trend lines, is presented in Figure 9.
One may notice from the figure that there is a consistent decreasing inflow trend over the whole, i.e., historic and future time periods, though the trend is stronger for the former-where the inflow has decreased by a factor of two from a maximum of ~2000 MCM/year in the 1990s to less than 1000 MCM/year at the present time-than the latter.Also, the negative trend is more pronounced for the more extreme RCP85 climate scenario than for the other two, RCP45 and RCP60.In addition, the results of the statistics of the average annually accumulated of the SWAT-simulated monthly inflow to the Boukan Dam, starting from the minimum over various percentiles to the maximum, together with the resulting water balance components, surface runoff (SWQ), lateral subsurface flow (LWQ), groundwater inflow (GWQ), and water yield (WYLD = SWQ + GWQ + LWQ), are presented for the three future periods under the three RCPs in Table 3.As can be seen from the table, the highest increase and decrease of the annual dam inflow are predicted for RCP60 an RCP85, respectively.When compared with the minimum and the mean historical dam inflow, those values will be increased in the near and middle future, except for RCP45 in the near future, whereas they will decrease by 2% to 23% for the far future period.The maximum dam inflow will augment by 22% to 66% in the time period 2020-2068, but decrease again by 6% to 31% after 2080.The low (25%)-quantile of the dam inflow will also experience a decrease, except for RCP45 in the middle future, whereas the high (75%) quantile will mostly increase, except of the far future period.The trends of the water balance components, e.g., of the water yield, are the same as the mean dam inflow.
The last column of Table 3 lists the Water Supply Stress Index WaSSI is defined as the ratio of the total water demand WD (extracted from Emami and Koch [13]) (second to last column) of all sectors to the total water supply from surface and groundwater sources, i.e., the water yield WYLD (third to last column).Obviously, low values of WaSSI < 1 mean low water stress and WaSSI-values reaching 1 and above indicate high stress.Based on the WaSSI-indices that were computed in this way, Table 3 shows that the water stress in the ZRB will be higher than that of the historical period for the near-and middle-future periods and rise to an even alarming level (WaSSI > 1) in the far-future period.

Crop Yield Simulation
As a major ingredient of the CSPSO-MODSIM crop pattern optimization model, the potential and actual crop yields of the major crops in the ZRB must be correctly known as these simulated and adjusted in the iterative optimization process, based on the prioritized water allocation and while using the FAO Equation ( 14) for the crop yield.
To simulate the crop production processes with the SWAT model, firstly the scheduled irrigation management operations are entered in its management module (.mgt).The most important management operations include planting, irrigation operation, and harvest operation, for which the needed data has been taken from Ahmadzadeh et al. and MOE [22,25], i.e., the irrigation operations are defined for the HRUs with the major crops and using their monthly crop water requirements (WR c ), in terms of their irrigation depths, listed for the seven crops in Table 4.These crop water depths are then later employed in the FAO equation, together with the available, prioritized water allocation Q c , to update the actual crop yield in the CSPSO-MODSIM model iteration process.The crop yields are computed in the SWAT-model, based on crop parameters that are specified in the crop.datinput file.A set of initial crop yield effective parameters, as described in Table 5, are adjusted based on literature values [22,60] and fine-tuned during the calibration/sensitivity analysis of the model to minimize the residuals of observed-simulated annual crop yields.The resulting final values of the crop parameters are also listed in Table 5.The SWAT-simulated crop yields are later applied as potential crop yields, Ymax c in the FAO Equation ( 14), as it was found that the amount of irrigation water dispersed to a crop was more than compensating its crop water need so that the crop yields remained the same when an unlimited source of irrigation water was applied.Using the calibrated crop parameters, the crop water requirements and the irrigation time scheme, as specified in Tables 4 and 5 Using the calibrated crop parameters, the crop water requirements and the irrigation time scheme, as specified in Tables 4 and 5, the average crop yields over the time period 1987-2012 are simulated and compared with the observed average crop yields in the ZRB.The regression-and the bar-plot of Figure 10 indicate that the SWAT-model simulates the observed actuals crop yields-which as mentioned are in fact potential crop yields, Ymaxc-in a satisfactory manner.For the economic analysis of the crop productions, various crop economic values including the guaranteed selling price by the government, production costs, crop yields, crop areas, irrigation depths, and potential AEWP of the seven major crops in the ZRB, gathered from [26,27] and listed in Table 6, are required.The table includes also potential future crops to be discussed later.
From Table 6 one may notice that the average total potential AEWPc for the current crop-areas situation is about 0.37 USD/m 3 , excluding APPL, grown in orchards and which are, for that reason, not considered in the multi-crop optimization and just applied as a demand node in the MODSIM model.In fact, APPL has a negative AEWPc, which means that its production cost is more than its selling price.In contrast, TOMA and WWHT have the highest-and BARL and POTA the lowest-yet positive agricultural economic water productivities AEWPc.If all the crop water requirements can be fulfilled completely and the potential crop yields be reached, the potential For the economic analysis of the crop productions, various crop economic values including the guaranteed selling price by the government, production costs, crop yields, crop areas, irrigation depths, and potential AEWP of the seven major crops in the ZRB, gathered from [26,27] and listed in Table 6, are required.The table includes also potential future crops to be discussed later.From Table 6 one may notice that the average total potential AEWP c for the current crop-areas situation is about 0.37 USD/m 3 , excluding APPL, grown in orchards and which are, for that reason, not considered in the multi-crop optimization and just applied as a demand node in the MODSIM model.In fact, APPL has a negative AEWP c , which means that its production cost is more than its selling price.In contrast, TOMA and WWHT have the highest-and BARL and POTA the lowest-yet positive agricultural economic water productivities AEWP c .If all the crop water requirements can be fulfilled completely and the potential crop yields be reached, the potential economic net benefit of the six crops (except APPL) will be ~32.2Million USD, but with an expense of 735 MCM total water use.In this section the CSPSO-MODSIM hydro-economic algorithm is executed under the two minimum arable cultivation-area constraints (Smin 1 and Smin 2 , see Table 1) and three different possible levels of cereal production rates (X = 50%, 35%, and 65%) to find the most suitable crop pattern, i.e., the one providing the average annual maximum net economic benefit, NEB t .This is done for the three RCP climate scenarios and the three (near, middle, far) future target periods.
The models starts by firstly ascribing initial values for 18 decision variables, i.e., the cultivated areas (A c ) of the 6 = 7 − 1 major crops (the apple demand node is excluded, because an orchard cannot be replaced by other crops) in each of the three demand areas (see Figure 4) (3 × 6 = 18 variables)-depending on the ranges of the individual crop areas for the two minimum area scenarios (Smin 1 and Smin 2 ), as specified in Table 1.
In the subsequent step, the MODSIM model is run, starting with the initial agricultural water productivities, AEWP c (Equation ( 12)), based on the current conditions (see Table 6).Then, the allocated water for each crop is captured from the MODSIM model under the different, and, while using the FAO Equation ( 14), the actual crop yield, the actual AEWP, and the main objective function Z, i.e., the total NEB t (Equation ( 12)), is estimated (updated) and the process repeated as part of the iteration scheme of the CSPSO-MODSIM process, wherefore the decision variables (crop areas) are adjusted based on the CSPSO-swarm information, until the maximum (=negative minimum) net annual economic benefit (=total NEB c ) (averaged over the simulation horizon and based on the allocated water supply) will be reached.

Optimal Future Net Annual Economic Benefits NEB t and AEWP
Results for the optimal future annual economic benefits NEB t obtained for the two minimum arable cultivation-area constraints (Smin 1 and Smin 2 ) under the various scenarios/time periods are shown in Table 7, wherefore, for simplification and strategic implications, the results in each group have been averaged over all three future time periods.
As can be seen from the table, for the Smin 1 -minimum crop area scenario, for both of the Smin 1 and Smin 2 -minimum crop area scenarios, the 65%-rate for cereal production is recommended over the total future period (up to year 2100) and this holds for all three RCPs and it is also in line with the increasing food demands of the country.However, as in the Smin 1 -scenario the minimum cultivated area for all non-cereal crops is set to zero (see Table 1), such a drastic extension of the cereal cultivation area might not generally be acceptable and will more likely increase the social dissatisfaction of the farmers and stakeholders in that region.
In contrast, for the Smin 2 -scenario, the optimal NEB t turns out to be quite higher than for the Smin 1 -scenario.Therefore, the results of Smin 2 are favored here, and for this reason they only will be discussed further in the following paragraphs.Thus, one may notice from the table that the total optimal NEB t s will be particularly high for the medium-emission scenarios, RCP45 and RCP60, with an annual net income of ~35.3 and 37.4 Million USD, respectively, which is up to 16% higher than the average potential net economic benefit of the historical period (=32.2Million USD, see Table 6), whereas for the high-emission scenario, RCP85, there will be no production gain.
Moreover, the highest annual economic benefits of the recommended crop patterns are expected to be in the middle future period, and this holds for all three RCPs, whereas the least crop production gain will be in the far future period, 2080-2098, as expected.
Table 7. Optimal Z (=total annual NEB t )-values (Million USD) for three kinds of cereal production-rates (X = 35%, 50% and 65%) under the two minimum arable cultivations areas constraints (Smin 1 and Smin 2 ) for three future periods and three RCPs.In Table 8, the results for the now-on considered Smin 2 -scenario are analyzed further in terms of the average optimal total AEWP for the six major crops of the ZRB and total agricultural supplied water Q per year.As can be seen from this table, for the selected scenario Smin 2 , the total AEWP is with 0.37 USD/m 3 the highest for the 65%-rate cereal-production scenario.Interestingly, for RCP85 the total AEWP will be the highest, i.e., it will have the lowest required water supply.This confirms that the algorithm will increase the AEWPs in the case of future water shortage.In fact, the total agricultural supplied water Q in this selected 65%-scenario will be decreased from the current 735 MCM/year (see explanations for Table 6) to 517 MCM/year for RCP45 and RCP60, and to 426 MCM/year for RCP85.
However, if the ultimate agricultural crop cultivation objective were to be shifted from maximum economic value to minimum agricultural water use, the 35% and 50% cereal production-rate scenario should be selected, respectively, for the near-middle and far future periods, which means that the "water-hungry" cereal production areas should not be extended that much.9, which lists the percentile optimal arable areas (Ac) of the six crops considered for the three irrigation plots in the ZRB, i.e., the 18 decision variables of the CSPSO-MODSIM model.Several conclusions can be drawn from the inspection of the multitude of numbers in the table.
Table 9. Optimal crop areas (km 2 and %, relative to those of the historical period) for the Smin2-scenario under the 65% cereal area constraint for the three irrigation plots (see Figure 4), the three future periods/RCPs.9, which lists the percentile optimal arable areas (Ac) of the six crops considered for the three irrigation plots in the ZRB, i.e., the 18 decision variables of the CSPSO-MODSIM model.Several conclusions can be drawn from the inspection of the multitude of numbers in the table.
Table 9. Optimal crop areas (km 2 and %, relative to those of the historical period) for the Smin2-scenario under the 65% cereal area constraint for the three irrigation plots (see Figure 4), the three future periods/RCPs.9, which lists the percentile optimal arable areas (A c ) of the six crops considered for the three irrigation plots in the ZRB, i.e., the 18 decision variables of the CSPSO-MODSIM model.Several conclusions can be drawn from the inspection of the multitude of numbers in the table.
Table 9. Optimal crop areas (km 2 and %, relative to those of the historical period) for the Smin 2 -scenario under the 65% cereal area constraint for the three irrigation plots (see Figure 4), the three future periods/RCPs.* ∆Ac: difference of the average crop acres (in km 2 or in %) over the future periods and of those of the historic period. 1 Irrigation plot downstream of the dam; 2 Irrigation plot upstream of the dam.

Crop
First of all, the total future arable area is, when compared with its historical value of 632 km 2 , increased to up to 900 km 2 , depending slightly on the future period and the RCP assumed, with the lowest extension of 36% being recommended for the far future/RCP60-and the near future/RCP85-time/RCP-, but of about 45% in the near future/RCP45-and the middle future/RCP85-scenarios.
Secondly, most of this arable area increase is made up by an extension of the cultivation areas of the two cereals, namely, WWHT, with ∆A c = 203 km 2 , i.e., nearly a doubling from its historical value, and, less so, of BARL, with ∆A c = 53 km 2 , in line with their importance as strategic crops in the ZRB.In fact, this is not surprising, as these cereals have, as compared with the other crops considered, rather low crop water demands, which, together with their relative high prices, lead to high AEWPs, particularly, for WWHT (see Table 6).
Thirdly, there is a tremendous reduction of the total ALFA crop area of more than −30%.This is the reason why, for fulfilling the demands of the alfalfa production in the region, Naraqi et al. [61] recommended to cultivate alfalfa in the riverside of Aras River, which is more suitable for this crop, and then transport it to the feeding stocks.
Fourthly, also significantly augmented by several 100%, though from a small base acreage, are the TOMA-cultivation areas, essentially, for the same reasons as for the cereals above.Somewhat smaller acreage increases, though still a doubling of the trifling base values, are experienced for the POTA and SGBT crop areas, namely, for the wetter RCP45-and RCP60-scenarios, with more water becoming available then.

Implications and Concluding Remarks
Interestingly, in connection with the afore-mentioned agro-econometric predictions for the ZRB, it has been recommended by the Lake Urmia Restoration Committee [62] to decrease the agricultural demands of the LU basin by about 40%, based on new executive strategies of demand management to avoid an imminent ecological and environmental disaster of the lake and the surrounding ecosystem.In fact, such a reduction of the future agricultural water irrigation is approximately found here, as, following Table 10, about 62% of the total agricultural demands of ZRB can be supplied for the RCP45 and RCP85 scenarios.
The recommended LU-disaster mitigation strategies [62] 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, and last, but not least, suggestions to replace some high water-consuming crops, like sugar beet and potato, with some less consuming ones, such as canola or sorghum.
Moreover, an increase of the crops' AEWP c can be generally also achieved by adopting proven agronomic and water management practices, such as deficit irrigation or modern irrigation technologies (e.g., pressured systems and drip irrigation).Furthermore, to improve the economic yield one may need to (a) switch from low-to high-value crops, for example, from apple to pistachio, (b) lower the costs of inputs (labor, water technologies), and (c) attempt to get multiple benefits of the irrigation water, e.g., using (cheaper) recycled wastewater [62].However, these adaptation strategies may have some further impacts that should be investigated, e.g., wastewater recycling will adversely increase the total water depletion.
It has also been proposed by Naraqi et al. [61] to replace parts of the low economic-value field crops and orchards with crops of higher economic benefits of less water use, such as canola, pistachio, or saffron.The initial AEWP c and other economic characteristics of these crops are also listed in Table 6 [26,27] and one can notice that these are indeed several times higher than the AEWP C of the ZRB major crops (see Table 6) analyzed heretofore.It may also be recommendable to develop further greenhouse vegetable cultivations.

Summary and Conclusions
In this paper, a complex, integrated hydro-economic model is developed as a crop pattern and irrigation planning tool consisting of a sequential combination of QM-climate downscaling, SWAT-hydrological modelling, MODSIM optimal water allocation, and a Constrained Stretched Particle Swarm Optimization (CSPSO) optimization search algorithm for maximizing the economic benefits of a multi-crop cultivation pattern.More specifically, the ultimate objective is then to optimize the future crop pattern of the major crops in the ZRB and their crop water irrigation depths in terms of net total economic benefits of the crop production, while considering the impact scenarios of the climate change and the cereal crop pattern scenarios.
To do this, climate change in the ZRB is first predicted using a statistical downscaling method, two-step updated QM (Quantile Mapping) bias correction technique that is based on the most suitable GCM outputs for the min./max.temperatures and precipitation namely CGCM3 and CESM-CAM5 of CMIP5 archive to assess the impact scenarios of the climate change (RCP 4.5, 6.0, and 8.5) in three 19-years future periods (near, middle, and far).In the next step, the downscaled weather variables are applied to the calibrated and validated basin-wide hydrologic model, SWAT, to simulate the future available water resources including the reservoir inflow and hydrologic changes.The crop yield potentials are also simulated using the SWAT model with applying the irrigation depths based on their monthly crop water requirements and adjusting the crop parameters.
According to the predicted impacts of climate change, RCP60 an RCP85 are expected to have the highest increase and decrease, respectively, of the inflow to the Boukan Dam.For all RCPs, the Boukan Dam inflow will be increased in the near and middle future, in comparison with the minimum and the mean historical dam inflow, except for RCP45 in the near future, whereas it is predicted that they will decrease by 2% to 23% for the far future period.The lowest available water resources are predicted for the far future regarding rather low precipitation and high temperature, especially for RCP85, which has the highest decrease of freshwater.The performance of the SWAT model for the streamflow and the crop yield simulation is quite satisfactory.The ratio of water demand across the water supply, i.e., the WaSSI (water stress) indices of the ZRB-simulated scenarios are predicted to be higher than the historical period for the near and medium periods, while the highest water stress is expected for the far future period.
To set up a water planning simulation module, the MODSIM model is then customized to allocate the available water based on priority constraints set for the ZRB.Finally, the CSPSO-MODSIM hydro-economic simulation-optimization model is developed by coupling this customized MODSIM model with a Constrained SPSO optimization algorithm, wherefore the objective function is here the total net economic benefit (NEB t )-which basically is the product of the total agro-economic water productivity (AEWP t ) and the total amount of water supplied Q t -under the constraints of different cereal crop pattern scenarios and the total arable areas existent for the three irrigation plots considered, i.e., the Boukan Dam downstream irrigation network, the plot upstream of the dam, and the RK plain.A penalty function method is employed to convert the constrained optimization problem to an unconstrained one, which can then be solved while using the classical SPSO algorithm.
The results of the application of the hydro-economic CSPSO-MODSIM model indicate that the optimal total net economic benefit (NEB t ) from all major crops in the ZRB basin (not including APPL, which is grown in orchards and cannot easily be altered over the short run) can be improved by appropriate adjustments of the different crop cultivation areas from that of the historical period (32.2 Million USD), to 35.3 and 37.4 Million USD for the two medium emission scenarios RCP45 and RCP60, respectively, but remain invariant for the extreme RCP85 (32.6 Million USD), though with an average reduction of the agricultural water use Q of 38%.
This overall increase of the future maximal NEB t is based on optimal crop pattern areas as decision variables that hint of a significant extension of the total arable area from its historical value of 632 km 2 to up to 900 km 2 , i.e., a relative increase that varies, depending on the future period/RCP, between 36% and 46% and so helps to guarantee the food security in the ZRB.Moreover, the strongest average increases of acreage of about 75%, amounting to an extension of approximately 250 km 2 , are reserved for the strategic crops BARL and WWHT, owing to the 65%-devoted areal model constraint for these cereals and also to their rather high estimated AEWPs and low water demands Q.
However, this large extension of the two cereal acreages is made primarily at the expense of the crop area for ALFA, which, because of its low AEWP and rather high water demand, will be decreased by more than 30% in the future, so that it becomes economically more efficient to import portions of that crop from a neighbour catchment, such as the Aras River Basin.
POTA-, SGBT-, and TOMA-crop areas, in particular, will, relative to their historic values, also be also tremendously increased, though less in absolute numbers than those of the cereals above.In this vegetable group, the TOMA share of crop pattern exhibits the largest extension of about 250%, owing to its high water productivity (AEWP) and this holds for all three RCPs investigated.For the wetter RCP45-and RCP60-scenarios, the POTA-and SGBT-crop areas will, because of more available water, be doubled as well.On the other hand, for the most extreme RCP85-scenario, BARL-and SGBT-areal shares will augment by an average 85%, whereas the POTA-acreage will then have the lowest percentile increase of only 50%.
Thus, in summary, to improve the total net income from agriculture in the ZRB, it is recommended to replace the high-consuming water and low agro-economic water productivity crops, such as ALFA, APPL, and SGBT with crops of less water demand and higher economic benefits that, in addition to the other major crops of the ZRB investigated here (BARL, WWHT, TOMA, POTA) could be some new region-specific crops, such as canola, saffron, and pistachio, which all have high absolute selling prices and relatively high AEWPs.
In conclusion, the sequential SWAT-MODSIM-CSPSO hydro-economic optimization model developed here turns out to be an efficient tool for quantitative studies of agricultural economics and water resources management, namely, in semi-arid basins with agricultural cultivation under deficit irrigation.Moreover, by driving the model by future climate predictors, climate change impacts on agricultural productivity can be evaluated, and the latter eventually be optimized under the vagaries of climate change.Using a model, as the one proposed here, will allow water resources authorities and other stakeholders, not only in agriculture, to find the most suitable regional water-and land management strategies and optimal land-use planning scenarios for future years.
As a caveat, though, it may be noted that the accuracy of the present hydro-economic optimization model could be improved by implementing more sophisticated crop yield forecasting methods, such as to allow the estimation of different crop-yield response factors Ky c for each stage of the crop growth; optimization of the crop production function; and, consideration of all the spatial and temporal changes of land use classes and their related water demands in the model.

Figure 1 .
Figure 1.Map of the Zarrine River Basin (ZRB) with the Boukan Dam and the climate and streamflow stations.

Figure 1 .
Figure 1.Map of the Zarrine River Basin (ZRB) with the Boukan Dam and the climate and streamflow stations.
Sustainability 2018, 10, x FOR PEER REVIEW 11 of 31 agricultural water rights, the environmental rights of the Boukan Dam, and the LU environmental water demand.

Figure 4 .
Figure 4. Schematics of the MODSIM model for the ZRB with related demand nodes, separated into three demand area groups, as indicated by the corresponding suffixes.

Figure 4 .
Figure 4. Schematics of the MODSIM model for the ZRB with related demand nodes, separated into three demand area groups, as indicated by the corresponding suffixes.

Figure 6 .
Figure 6.Spatial distributions of average precipitation (top) and temperature (bottom) in the ZRB during the near future period (2020-2038) for three RCPs.

Figure 9 .
Figure 9. Annual Soil Water Assessment Tool (SWAT)-simulated inflow to the Boukan Dam in the historical (1991-2012) and the three future periods, for the three RCPs, together with corresponding linear trend lines.
, the average crop yields over the time period 1987-2012 are simulated and compared with the observed average crop yields in the ZRB.The regression-and the bar-plot of Figure 10 indicate that the SWAT-model simulates the observed actuals crop yields-which as mentioned are in fact potential crop yields, Ymax c -in a satisfactory manner.Sustainability 2018, 10, x FOR PEER REVIEW 20 of 31

Figure 11 .
Figure 11.Current/historic crop pattern of the Boukan Dam (downstream) network-(left pie) and the dam upstream demand area (right pie) in km 2 .

Figure 12 .
Figure 12.Optimal crop pattern for RCP45 for the three future periods: the near (left panel), middle (middle panel) and far (right panel) future in km 2 .

Figure 12 .
Figure 12.Optimal crop pattern for RCP45 for the three future periods: the near (left panel), middle (middle panel) and far (right panel) future in km 2 .Sustainability 2018, 10, x FOR PEER REVIEW 24 of 31

Figure 14 .
Figure 14.Similar to Figure 11, but for RCP85.The pie charts of Figures 11-14 are quantified more clearly in Table9, which lists the percentile optimal arable areas (A c ) of the six crops considered for the three irrigation plots in the ZRB, i.e., the 18 decision variables of the CSPSO-MODSIM model.Several conclusions can be drawn from the inspection of the multitude of numbers in the table.
Details of each of the above modeling steps are described in the following sub-sections.

Table 2 .
Goodness of fitness measures of the QM-model for the different RCPs for the time period 2006-2015.

Table 3 .
Statistics of the inflow of the Boukan Dam, basin water balance components, total water demand (WD) (MCM/year), and WaSSI-index for historic and future periods under different RCPs, with % -values denoting relatives to the historic reference period.

Table 4 .
Monthly crop water requirements for the major crops in the ZRB and their irrigation intervals.

Table 5 .
List of effective crop yield parameters adjusted in the SWAT-calibration process.

Table 6 .
Current selling prices, production costs, crop yields, crop areas, crop water requirements and initial agro-economic water productivity (AEWP) of the presently cultivated major crops as well as of potential future crops in the ZRB.