Water Quality Modeling of Mahabad Dam Watershed–Reservoir System under Climate Change Conditions, Using SWAT and System Dynamics

The total phosphorus (TP) concentration, as the primary limiting eutrophication factor in the Mahabad Dam reservoir in Iran, was studied, considering the combined impacts of climate change, as well as the scenarios on changes in upstream TP loadings and downstream dam water allocations. Downscaled daily projected climate data were obtained from the Beijing Normal University Earth System Model (BNU-ESM) under moderate (RCP4.5) and extreme (RCP8.5) scenarios. These data were used as inputs of a calibrated Soil and Water Assessment Tool (SWAT) model of the watershed in order to determine the effects of climate change on runoff yields in the watershed from 2020 to 2050. The SWAT model was calibrated/validated using the SUFI-2 algorithm in the SWAT Calibration Uncertainties Program (SWAT-CUP). Moreover, to model TP concentration in the reservoir and to investigate the effects of upstream/downstream scenarios, along with forecasted climate-induced changes in streamflow and evaporation rates, the System Dynamics (SD) model was implemented. The scenarios covered a combination of changes in population, agricultural and livestock farming activities, industrialization, water conservation, and pollution control. Relative to the year 2011 in which the water quality data were available, the SD results showed the highest TP concentrations in the reservoir under scenarios in which the inflow to the reservoir had decreased, while the upstream TP loadings and downstream dam water allocations had increased (+29.9%). On the other hand, the lowest TP concentration was observed under scenarios in which upstream TP loadings and dam water allocations had decreased (−18.5%).


Introduction
Human activities have influenced the water quality in aquatic ecosystems by altering nutrient fluxes into receiving water bodies. Complex interactions of inflow, weather, and soil and land use practices impact external nutrient loadings, which are the main driving forces of eutrophication in reservoirs. Moreover, the recent global warming is expected to affect nutrient loss dynamics in watersheds by changing atmospheric and meteorological properties, such as precipitation patterns, atmospheric water vapor, and evaporation. This situation could make lakes and reservoirs more vulnerable to eutrophication [1][2][3]. In recent years, climate change and its impacts on the quantitative and qualitative aspects of water have been the focus of several studies [4][5][6]. Various climate models and scenarios were used to investigate possible climate change impacts on hydrological variables [7,8].
reservoirs. Previous studies provided a variety of scenarios for management practices, while the effects of global warming and how it might interfere with other possible future anthropogenic alterations have not received sufficient attention. Considering the crucial importance of the issue, there is a need for an accurate and more credible assessment of the cumulative effects of the climate change phenomena on watershed-reservoir systems.
In this regard, to study climate change impacts on the Mahabad Dam watershed-reservoir system, downscaled daily climate predictions from the Beijing Normal University Earth System Model (BNU-ESM), under moderate (RCP4.5) and extreme (RCP8.5) scenarios, were used as input to a calibrated SWAT Ver. 2012.10.21 model of the Mahabad Dam watershed. The SWAT model was used to simulate streamflow in the watershed. Moreover, the Sequential Uncertainty FItting Ver. 2 (SUFI-2) program of the Calibration and Uncertainty Procedures (SWAT-CUP) Ver. 5.1.6.2 program was implemented for the sensitivity analysis, calibration, and validation of the developed SWAT model. Finally, the effects of upstream/downstream scenarios (increase or decrease in pollution and water allocations), along with the climate change impacts on TP concentration in the Mahabad Dam reservoir were evaluated using the Stella Ver. 9.0.2 SD model. These simulations were carried out for the period 2020-2050, and the average TP concentration values in the reservoir were compared with the values during 2011 in which the water quality data were available. Figure 1 shows the framework of the present study.
Water 2019, 11, x FOR PEER REVIEW 3 of 17 alterations have not received sufficient attention. Considering the crucial importance of the issue, there is a need for an accurate and more credible assessment of the cumulative effects of the climate change phenomena on watershed-reservoir systems.
In this regard, to study climate change impacts on the Mahabad Dam watershed-reservoir system, downscaled daily climate predictions from the Beijing Normal University Earth System Model (BNU-ESM), under moderate (RCP4.5) and extreme (RCP8.5) scenarios, were used as input to a calibrated SWAT Ver. 2012.10.21 model of the Mahabad Dam watershed. The SWAT model was used to simulate streamflow in the watershed. Moreover, the Sequential Uncertainty FItting Ver. 2 (SUFI-2) program of the Calibration and Uncertainty Procedures (SWAT-CUP) Ver. 5.1.6.2 program was implemented for the sensitivity analysis, calibration, and validation of the developed SWAT model. Finally, the effects of upstream/downstream scenarios (increase or decrease in pollution and water allocations), along with the climate change impacts on TP concentration in the Mahabad Dam reservoir were evaluated using the Stella Ver. 9.0.2 SD model. These simulations were carried out for the period 2020-2050, and the average TP concentration values in the reservoir were compared with the values during 2011 in which the water quality data were available. Figure 1 shows the framework of the present study.

Study Area and Data
The Mahabad Dam watershed (808 km 2 ) is located in the West-Azerbaijan province in the northwest of Iran (3644 N, 4539 E) and is one of the Urmia Lake subbasins. According to the meteorological records, during 1988 to 2012, the average annual temperature and precipitation in the area were 12 °C and 350 mm, respectively [43]. This watershed is mostly covered by agricultural fields and grasslands. The Kauter and Beytas Rivers originate from the southern heights of the plain and run to the north in parallel. They join and create the Mahabad Dam reservoir and continue running as the Mahabad River. The soil map and land use map of the Mahabad Dam watershed are

Study Area and Data
The Mahabad Dam watershed (808 km 2 ) is located in the West-Azerbaijan province in the northwest of Iran (36 • 44 N, 45 • 39 E) and is one of the Urmia Lake subbasins. According to the meteorological records, during 1988 to 2012, the average annual temperature and precipitation in the area were 12 • C and 350 mm, respectively [43]. This watershed is mostly covered by agricultural fields and grasslands. The Kauter and Beytas Rivers originate from the southern heights of the plain and run to the north in parallel. They join and create the Mahabad Dam reservoir and continue running as the Mahabad River. The soil map and land use map of the Mahabad Dam watershed are presented in Figure 2. Based on the records from hydrometric stations, the average flow rate of the Kauter and Beytas Rivers during 1988-2012 were 6.18 m 3 /s and 1.82 m 3 /s, respectively [44].
Water 2019, 11, x FOR PEER REVIEW 4 of 17 presented in Figure 2. Based on the records from hydrometric stations, the average flow rate of the Kauter and Beytas Rivers during 1988-2012 were 6.18 m 3 /s and 1.82 m 3 /s, respectively [44].
(a) (b)  The Mahabad Dam has a storage capacity of 200 million cubic meters and provides water for agriculture (71%), industry (11%), drinking (7%), and other miscellaneous purposes (11%). In recent years, the Mahabad Dam has experienced various environmental issues due to excess nutrient loading from its watershed, causing year-round eutrophication in the reservoir. Based on field measurements, phosphorus is the rate-limiting nutrient in the reservoir, which generates from upstream land use practices (72%), residential areas (16%), and livestock farming activities (12%) [44].
In order to perform this research study, downscaled and bias-corrected future daily precipitation and temperature data forecasted by the BNU-ESM were downloaded from the National Aeronautics and Space Administration (NASA) Earth Exchange (NEX) archive under RCP4.5 and RCP8.5 scenarios [45]. Moreover, in order to delineate the watershed, a Digital Elevation Model (DEM) with 30 m spatial resolution was acquired from the Advanced Spaceborne Thermal Emission and  The Mahabad Dam has a storage capacity of 200 million cubic meters and provides water for agriculture (71%), industry (11%), drinking (7%), and other miscellaneous purposes (11%). In recent years, the Mahabad Dam has experienced various environmental issues due to excess nutrient loading from its watershed, causing year-round eutrophication in the reservoir. Based on field measurements, phosphorus is the rate-limiting nutrient in the reservoir, which generates from upstream land use practices (72%), residential areas (16%), and livestock farming activities (12%) [44].
In order to perform this research study, downscaled and bias-corrected future daily precipitation and temperature data forecasted by the BNU-ESM were downloaded from the National Aeronautics and Space Administration (NASA) Earth Exchange (NEX) archive under RCP4.5 and RCP8.5 scenarios [45]. Moreover, in order to delineate the watershed, a Digital Elevation Model (DEM) with 30 m spatial resolution was acquired from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global DEM data source [46]. The soil map was retrieved from the Food and Agriculture Organization of the United Nations (FAO) soils portal using the Harmonized World Soil Database (HWSD) V. 1.2 [47]. The land use data were provided by Mahab Ghodss Consulting Engineering Company [44], and the meteorological data were provided by the Islamic Republic of Iran Meteorological Organization (IRIMO) for the weather station of Mahabad [43].

The Climate Change Model and Scenarios
General Circulation Models (GCMs) have been widely used to make climate predictions on seasonal to decadal timescales. GCMs represent the physical atmospheric and oceanic processes. On the other hand, Earth System Models (ESMs) include physical, chemical, and biological processes. Therefore, they reach far beyond GCMs and simulate all relevant aspects of the earth's system. The Beijing Normal University Earth System Model (BNU-ESM) was developed at Beijing Normal University and can simulate several observed features of the earth's climate system with high accuracy. The model can be used to study climate variability at different timescales, interactions between the ocean and atmosphere, and carbon-climate feedback [48].
Representative Concentration Pathways (RCPs) are four greenhouse gas concentration trajectories adopted by the Intergovernmental Panel on Climate Change (IPCC) for its Fifth Assessment Report (AR5) in 2014 [49]. RCP2.6, RCP4.5, RCP6.5, and RCP8.5 are groups of several scenarios that consider future atmospheric conditions and land use/land cover changes. RCPs are used to investigate the response of the climate system to various anthropogenic greenhouse gas emissions. The numbers associated with them indicate the predicted amount of increase of radiative forcing to be reached by the year 2100 [50].
In this study, the BNU-ESM under RCP4.5 and RCP8.5 was implemented in order to obtain temperature, precipitation, and evaporation values for the period of 2020 to 2050.

The Soil and Water Assessment Tool
Among the most commonly used continuous-time, semi-distributed, and physically-based models is the SWAT model, which can address different pollution problems for watershed scales. To model processes within a watershed, the SWAT integrates weather, surface and groundwater hydrology, soil properties, plant growth, nutrient cycles, and land management practices [51]. Based on interior outlet points along the stream network, the SWAT divides the watershed area into several subbasins. Then in each subbasin, areas with similar soil types, land uses, slopes, and management practices are subdivided into Hydrologic Response Units (HRUs). HRUs represent percentages of the total subbasin area. In each HRU, yields are calculated and then summed to determine the total subbasin output. This means that at the subbasin scale, the SWAT uses spatially distributed parameterization, while at the HRU scale, it uses lumped parameterization [52]. To delineate the watershed, the SWAT incorporates four primary data files: a DEM, a land use map, a soil map, and meteorological data (precipitation, minimum and maximum temperatures, wind, relative humidity, and solar radiation). In cases where any of the meteorological data are not introduced into the model, the SWAT utilizes a built-in weather generator model (WXGEN) to stochastically generate daily weather values based on historical monthly averages of parameters such as temperature, precipitation, relative humidity, wind, and solar radiation.
In this study, the outputs of the BNU-ESM under RCP4.5 and RCP8.5 were used as the input of the SWAT model to simulate streamflow in the watershed during 2020 to 2050.

The SWAT Calibration and Uncertainty Procedures
The SWAT-CUP program has been developed for the calibration, validation, and sensitivity analysis of the SWAT model parameters. Prior to the model calibration, more sensitive parameters have to be identified. In the SWAT-CUP, t-stat values and p-values are used to measure the sensitivities of parameters. Parameters that show higher t-stat values (lower p-values) are more sensitive and a have more significant impact on the target variable [53].
Calibration means adjusting the model input parameters with the goal of achieving the best fit between the observed and simulated values. In the SWAT-CUP, the goodness of calibration is measured using the p-factor (the fraction of data in the range of a 95% prediction uncertainty (95 ppu)) and the r-factor (the average thickness of the 95 ppu band, divided by the standard deviation of the observed data). The p-factor is a value between 0 and 1, and the r-factor has a range of 0 to ∞. When the p-factor = 1 and the r-factor = 0, the simulated model is precisely in accordance with the observed data. In general, p-factors greater than 0.7 and r-factors smaller than 1.5 show satisfactory calibration and validation results [54].
The SWAT-CUP uses five different calibration procedures: SUFI-2, Particle Swarm Optimization (PSO), Generalized Likelihood Uncertainty Estimation (GLUE), Parameter Solution (ParaSol), and Markov Chain Monte Carlo (MCMC). For large-scale models in which the calibration process can be very time-consuming, the semiautomated SUFI-2 is quite efficient [55]. The SWAT-CUP adjusts model parameters within a user-defined range to achieve the desired value of an objective function, such as the coefficient of determination (R 2 ) (Equation (1)), the Nash-Sutcliffe (NS) model efficiency coefficient (Equation (2)), and percent bias (PBIAS) (Equation (3)): where "Q" is a variable such as streamflow; "m" and "s" stand for measured and simulated; the bar indicates the average; and "i" is the ith measured or simulated value. Higher R 2 values show that the model fits the observed values better. The NS function has a range of −∞ to 1. NS = 1 corresponds to a perfect match of simulated values to the observed data. The values between 0 and 1 indicate that the simulated and observed values are close to each other, whereas values less than 0 show that the model has no predictive power [56]. Moreover, percent bias measures the average tendency of the simulated data to be larger or smaller than the observations. The optimum value is zero, where low magnitude values indicate better simulations. Positive values indicate model underestimation, and negative values indicate model over estimation [57].

System Dynamics Modeling and Scenarios
The computer-based, object-oriented SD modeling is based on feedback control and nonlinear dynamics to model complex systems. It assumes that to obtain the behavior of a whole system, parameters such as time delays, nonlinearities, system feedbacks, amplifications, and structural relationships between a system's elements should be studied instead of the individual components themselves. SD provides a more in-depth understanding, as well as a dynamic view of the behavior of complex systems and how they evolve. The first step in developing an SD model is to create a conceptual model that is generally referred to as a causal loop diagram. This dynamic hypothesis is then quantified and simulated using stock and flow diagrams [58]. Stocks, flows, converters, and connectors are the basic building blocks for SD models. Stocks represent accumulations, while flows transport quantities into or out of a stock. In water resource science, a poor understanding of the interconnections among different subsystems inhibits researchers from developing sustainable solutions [59]. Hence, decision-making in water resource engineering should be based on a holistic view that considers the system's sophisticated dynamics and feedback processes, as well as the interdependencies between the different processes. A consideration of the combined effects of system dynamics can improve management decisions and reduce the possibilities of adverse side effects or unintended consequences of policy decisions [60].
In order to investigate how upstream/downstream scenarios, along with climate change impact, affect the TP concentration in the reservoir, six scenarios were proposed, as presented in Table 3. Although the change in precipitation patterns changes TP loadings in the watershed, only the streamflow was simulated using the SWAT model, and it was assumed that under each scenario, the TP loadings change in a range of 10-30% higher or lower than the observed values in 2011. This gives the researchers more freedom in studying the impacts of proposed anthropogenic changes on upstream TP loadings.  Although the change in precipitation patterns changes TP loadings in the watershed, only the streamflow was simulated using the SWAT model, and it was assumed that under each scenario, the TP loadings change in a range of 10-30% higher or lower than the observed values in 2011. This gives the researchers more freedom in studying the impacts of proposed anthropogenic changes on upstream TP loadings.   The results show that under climate change conditions and comparing the average values during 2020-2050 with 2011, temperature and evaporation will increase, while precipitation will decrease. Compared with the temperature in 2011, the RCP4.5 and RCP8.5 predicted an 18.11% and a 22.50% increase in the average temperature during 2020-2050, respectively (Figure 3a). Comparing the same periods, the RCP4.5 and RCP8.5 predicted a 2.93% and 2.25% decrease in precipitation, respectively (Figure 3b). Moreover, the evaporation will increase 19% under the RCP4.5, and 23% under the RCP8.5 (Figure 3c). The ultimate result of these changes will be less rainfall and higher evaporation rates, which can impact both the quantitative and qualitative aspects of the water resources in the study area.

Climate Change Impacts on Streamflow
The watershed was delineated within the ArcGIS interface using the ArcSWAT automatic watershed delineation tool. The streams were laid out based on the model recommended minimum drainage area of 1589.78 ha, and a total number of 45 subbasins and 165 HRUs were formed. HRUs were generated based on the land use and soil types that made up at least 20% of a given subbasin's area. Figure 4 shows the delineated Mahabad Dam watershed in the ArcSWAT. The watershed was delineated within the ArcGIS interface using the ArcSWAT automatic watershed delineation tool. The streams were laid out based on the model recommended minimum drainage area of 1589.78 ha, and a total number of 45 subbasins and 165 HRUs were formed. HRUs were generated based on the land use and soil types that made up at least 20% of a given subbasin's area. Figure 4 shows the delineated Mahabad Dam watershed in the ArcSWAT.   Using monthly streamflow data, the SWAT model was calibrated from 1988 to 2004. These data were obtained from the hydrometric stations located at the Kauter and Beytas Rivers' entry points to the reservoir. A three-year model warmup from 1988 to 1991 was set in order to stabilize base-flow conditions in the model. Following calibration, the streamflow was validated for the period of 2005-2012. Table 4 shows the calibrated values for the most sensitive parameters affecting streamflow in this watershed, and Figure 5a,b show the observed and simulated streamflow values for the Kauter and Beytas Rivers, respectively. According to criteria set by Moriasi et al. (2007) [56] for evaluating the model performance in calibration and validation, an NS value between 0.50 and 0.65 is considered "satisfactory," a value between 0.65 and 0.75 is rated as "good", and "very good" is attributed to values between 0.75 and 1.00. Therefore, except for the streamflow validation in the Beytas River, which fits into the satisfactory category, the other calibration and validation values were very good. Moreover, the PBIAS value for the streamflow calibration in the Kauter River is positive, which indicates model underestimation. On the other hand, the PBIAS value for the streamflow validation in the Kauter River and the streamflow calibration and validation in the Beytas River are negative, indicating model overestimation.
After performing the calibration and validation, the outputs of the BNU-ESM, under the RCP4.5 and RCP8.5 scenarios were used to run the SWAT model from 2020 to 2050 in order to simulate the future streamflow in the watershed. Figure 6 shows the average yearly streamflow in the watershed from 2020 to 2050 under climate change conditions. According to criteria set by Moriasi et al. (2007) [56] for evaluating the model performance in calibration and validation, an NS value between 0.50 and 0.65 is considered "satisfactory," a value between 0.65 and 0.75 is rated as "good", and "very good" is attributed to values between 0.75 and 1.00. Therefore, except for the streamflow validation in the Beytas River, which fits into the satisfactory category, the other calibration and validation values were very good. Moreover, the PBIAS value for the streamflow calibration in the Kauter River is positive, which indicates model underestimation. On the other hand, the PBIAS value for the streamflow validation in the Kauter River and the streamflow calibration and validation in the Beytas River are negative, indicating model overestimation.
After performing the calibration and validation, the outputs of the BNU-ESM, under the RCP4.5 and RCP8.5 scenarios were used to run the SWAT model from 2020 to 2050 in order to simulate the future streamflow in the watershed. Figure 6 shows the average yearly streamflow in the watershed from 2020 to 2050 under climate change conditions. As precipitation is forecasted to decrease and temperature and evaporation to rise, the SWAT model predicted reduced streamflow under both RCPs. The overall average streamflow of thirty years, simulated under the RCP4.5 and RCP8.5 scenarios, predicts a 5.1% and 3% decrease in As precipitation is forecasted to decrease and temperature and evaporation to rise, the SWAT model predicted reduced streamflow under both RCPs. The overall average streamflow of thirty years, simulated under the RCP4.5 and RCP8.5 scenarios, predicts a 5.1% and 3% decrease in streamflow, respectively, compared to the average streamflow in 2011. Figure 7 shows the developed stock and flow diagram. Three stocks were used in the model to represent the reservoir volume, TP load in the reservoir, and TP load in sediments. As illustrated, the reservoir volume is a function of inflows, outflows, and evaporation. The TP load in the reservoir fluctuates based on inflows and outflows, along with settling and resuspension processes. Furthermore, the TP load in sediments changes with settling, resuspension, and burial processes. In the next step, the TP concentration in the reservoir was calibrated by adjusting the initial TP load in the reservoir, initial TP load in sediments, TP settling rate, TP resuspension rate, and TP burial rate. Table 5 presents the parameters and their calibrated values used to calibrate the TP concentration in the reservoir.  In the next step, the TP concentration in the reservoir was calibrated by adjusting the initial TP load in the reservoir, initial TP load in sediments, TP settling rate, TP resuspension rate, and TP burial rate. Table 5 presents the parameters and their calibrated values used to calibrate the TP concentration in the reservoir.  Table 6 summarizes the simulation results of the average annual TP concentration in the reservoir. In each scenario, the inflow and evaporation rates represent the average change during 2020-2050, due to climate change conditions. Other parameters in each of the six scenarios were varied between 10 and 30 percent higher or lower than the observed value during 2011. Since the water quality records in the reservoir were only available for one year, the average projected TP concentration values were compared with the TP concentration in the reservoir during 2011 (84.13 µg/L). The results show that Scenario 2 (increase in agriculture and livestock farming activities) created the maximum TP concentrations in the reservoir under both RCPs. Relative to the year 2011, the TP concentration in the reservoir increased from 10.8% to 29.9% for RCP4.5 and 9.1% to 26.2% for RCP8.5. In this scenario, the inflow to the reservoir decreased, while the TP loading from upstream, water use downstream, and evaporation rate increased. Since 72% of the watershed phosphorus load generates from upstream land use practices, such as agriculture and livestock farming activities, the TP concentration in the reservoir showed a high sensitivity to changes in loadings from these pollution sources. Moreover, as agriculture takes up a significant portion of the dam's water allocation, the decreased volume of the reservoir, as a result of the increased agricultural water use, resulted in higher TP concentrations in the reservoir.

TP Concentration in the Reservoir
Among all scenarios, Scenario 6 (water conservation and pollution control) resulted in the minimum TP concentrations in the reservoir under both RCPs. Relative to the year 2011, the TP concentration in the reservoir decreased from 3.4% to 17.5% for RCP4.5 and 4.6% to 18.5% for RCP8.5. In this scenario, the phosphorus loading from all pollution sources decreased, as well as the outflow from the dam.
Based on the findings of this study, even under the most optimistic scenarios, serious conservation tactics will be necessary to lower the reservoir TP concentrations to oligotrophic (0-12 µg/L) or mesotrophic (12-24 µg/L) levels. To counteract this situation, adaptations to reduce upstream nutrient loading and restrictions on dam water allocations are required. Adjustments can include nutrient and soil management practices, as well as procedures to minimize nutrient loss to surface waters, such as the establishment of wetlands and riparian buffer zones, and restrictions on dam water allocations, such as less intensive agriculture and reduced domestic water use.

Conclusions
Using the BNU-ESM under the RCP4.5 (moderate) and RCP8.5 (extreme) scenarios, as well as the SWAT and an SD approach, the combined impacts of climate change scenarios on upstream TP loadings and dam water allocations were investigated on the TP concentration in the Mahabad Dam reservoir in Iran. This study indicated that climate change would have a significant impact on the meteorological conditions and the water yield in the Mahabad Dam watershed during the period of 2020-2050. Comparing the average value during this period with the average value in 2011, under the RCP4.5 and RCP8.5, the BNU-ESM predicted an 18.11% and a 22.50% increase in temperature, and a 2.93% and 2.25% decrease in precipitation, respectively. Consequently, as a result of altered climatic conditions, the SWAT predicted a 5.1% and 3.0% decrease in streamflow under RCP4.5 and RCP8.5, respectively. Scenarios on population increase, changes in upstream pollution rates, and dam water allocations, along with streamflow and evaporation rate alterations due to climate change, showed that the trophic state of the Mahabad Dam reservoir would deteriorate under scenarios in which the TP loadings from upstream and water use downstream increased. Therefore, the scenario of increased agriculture and livestock activities yielded the highest TP concentration in the reservoir (109.3 µg/L). On the other hand, the TP concentration in the reservoir showed the lowest values under scenarios in which the TP loadings from upstream and water use downstream decreased (68.6 µg/L). This situation was most evident in the scenario of water conservation and pollution control. As the results of this study indicated, even the most optimistic scenarios of TP loadings and dam water allocations still created eutrophic conditions in the reservoir. This situation demands serious bans or limiting strategies on activities that generate pollution upstream and the precise management of dam water allocations to improve the trophic state of the reservoir.