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

Abstract

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%).

1. 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]. Watershed-scale models were implemented to study the effects of climate change and land management practices on watershed yields. Moreover, System Dynamics (SD) modeling, based on the notion of systems thinking, has emerged to facilitate a holistic analysis of complex human–environmental systems, such as water resource systems [9].

According to the literature, in a study by Woznicki and Nejadhashemi (2012) [10], the Soil and Water Assessment Tool (SWAT) outputs, under future climate change scenarios, showed increased sediment, Total Nitrogen (TN), and Total Phosphorus (TP) yields in agricultural and rangeland watersheds in Kansas and Nebraska. Michalak et al. (2013) [11] studied the Maumee River basin using SWAT and concluded that the precipitation intensity, timing of fertilizer application, and tillage practices affect the Dissolved Reactive Phosphorus (DRP) yields. Moreover, Bosch et al. (2014) [12] studied the impact of climate change on Lake Erie until 2100. The SWAT model of the Maumee River basin showed a 5–11% increase in streamflow, a 2–32% increase in sediment yield, a 1–5% increase in DRP, and a 0–7% increase in TP loading. The authors illustrated that Best Management Practices (BMPs) could significantly reduce the amount of climate change-induced increases in sediment and nutrient loads. In another research study by Malagó et al. (2017) [13], the SWAT model was used to model water and nutrient fluxes in the Danube River basin in Europe, which is under great nutrient pressure. The authors concluded that the main sinks of TN and TP diffuse emissions were plant uptake, soil retention, riparian filter strips, and river retention. More recently, Du et al. (2019) [14] implemented the SWAT model to study the impacts of different land use and climate change scenarios on the runoff yield of the Dagu River basin in China. The authors found climate change impacts more influential on runoff than land use change. In another study, Abbasi et al. (2019) [15] used the SWAT model to simulate the discharge of sediment and pesticides into the Malewa River Basin in Kenya. Based on the model results, higher concentrations of pesticides were found between May and mid-July. Finally, the authors concluded that considering the similarity between measured and simulated pesticides, the SWAT model could be used as an initial evaluation modeling tool for upstream to downstream suspended sediment and pesticide transport in catchments.

Furthermore, SD has been used for various water management [16,17,18,19] and environmental management studies [20,21,22,23]. More specifically, SD has been used for studies involving flood management [24,25,26,27], water allocation [28,29,30], climate change impact on water [31,32,33], carbon footprint of water [34,35,36], and energy planning [37]. Models have also been developed [38] for Lake Mead and the Las Vegas water supply system to educate the public about water conservation. A review of SD applications in water resources is provided by Mirchi et al. (2012) [39].

Dawadi et al. (2013) [32] investigated the effects of climate change and population increase on the water resources of the Las Vegas Valley (LVV) in southern Nevada. Using an SD model, the authors studied the impact of climate change on water demand and the water supply from the Colorado River. Mirchi and Watkins Jr. (2013) [40] developed an integrated SD model to simulate the natural processes driving eutrophication in Lake Allegan in Michigan and the interactions between socioeconomic subsystems. The model was used to characterize the lake’s recovery from its hypereutrophic state and to assess the effectiveness of a number of proposed Total Maximum Daily Load (TMDL) reduction policies. In another study, Liu et al. (2015) [41] developed an SD framework for managing the water quality of the Dianchi Lake in Yunnan Province in China under four scenarios: the business-as-usual, spatial adjustment of industries and population (S1), wastewater treatment capacity construction (S2), and structural adjustment of agriculture (S3). The authors concluded that S2 was the most effective for improving the lake’s water quality. Moreover, Duran-Encalada (2017) [42] studied the US–Mexico border region of the Rio Grande/Rio Bravo Water Basin, estimated the variation in the quality/quantity of water due to climate change, and assessed its impact on community development. The authors proposed an SD model to understand the complex interaction of factors affecting water quality and quantity and their impacts on social and economic conditions.

Based on the literature review, numerous studies were performed through isolated approaches to investigate the hydrological climate change impact solely or its impact on either watersheds or 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.

2. Materials and Methods

2.1. Study Area and Data

The Mahabad Dam watershed (808 km²) 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³/s and 1.82 m³/s, respectively [44].

Table 1. The land use classification of the Mahabad Dam watershed.

Land Use Type	Watershed Area (%)
Dryland Farming	66.37
Dense Pasture	13.82
Non-Dense Pasture	11.09
Irrigated Farming	4.49
Forest	2.95
Water	1.03
Urban (Medium Density)	0.13

and

Table 2. The soil types in the Mahabad Dam watershed.

Soil Type	Sand (%)	Silt (%)	Clay (%)	Watershed Area (%)
Taconic	43	35	23	72
Benson	35	37	30	28

represent the land use and soil classification of the watershed, respectively.

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].

2.2. 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.

2.3. 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.

2.4. 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 $\infty$. 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²) (Equation (1)), the Nash–Sutcliffe (NS) model efficiency coefficient (Equation (2)), and percent bias (PBIAS) (Equation (3)):

$${\mathrm{R}}^2=\frac{[{\sum }_{\mathrm{i}}({\mathrm{Q}}_{\mathrm{m},\mathrm{i}}-{\overline{\mathrm{Q}}}_{\mathrm{m}})\left({\mathrm{Q}}_{\mathrm{s},\mathrm{i}}-{\overline{\mathrm{Q}}}_{\mathrm{s}}\right){]}^2}{{\sum }_{\mathrm{i}}{\left({\mathrm{Q}}_{\mathrm{m},\mathrm{i}}-{\overline{\mathrm{Q}}}_{\mathrm{m}}\right)}^2{\sum }_{\mathrm{i}}{\left({\mathrm{Q}}_{\mathrm{s},\mathrm{i}}-{\overline{\mathrm{Q}}}_{\mathrm{s}}\right)}^2}$$

(1)

$$\mathrm{NS}=1-\frac{{\sum }_{\mathrm{i}}{\left|{\mathrm{Q}}_{\mathrm{m}}-{\mathrm{Q}}_{\mathrm{s}}\right|}_{\mathrm{i}}^2}{{\sum }_{\mathrm{i}}{\left|{\mathrm{Q}}_{\mathrm{m},\mathrm{i}}-{\overline{\mathrm{Q}}}_{\mathrm{m}}\right|}_{\mathrm{i}}^2}$$

(2)

$$\mathrm{PBIAS}=100\times \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Q}}_{\mathrm{m}}-{\mathrm{Q}}_{\mathrm{s}}\right)}_{\mathrm{i}}}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{Q}}_{\mathrm{m},\mathrm{i}}}$$

(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² 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].

2.5. 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. System Dynamics (SD) model scenarios.

#	Scenario	Description
1	Population increase	Increase in phosphorus loading from residential sources and increase in domestic water use
2	Increase in agriculture and livestock farming activities	Increase in phosphorus loadings from agricultural fields and livestock farming activities and increase in agricultural water use
3	Industrialization	Increase in industrial water use
4	Water conservation	Decrease in water allocations
5	Pollution control	Decrease in phosphorus loadings from pollution sources
6	Water conservation and pollution control	Decrease in water allocations and phosphorus loadings from pollution sources (combination of Scenario 4 and 5)

. 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.

3. Results and Discussion

3.1. Climate Change Impacts on Temperature, Precipitation, and Evaporation

Figure 3 shows the average future temperature, precipitation, and evaporation values for the period of 2020–2050 and for the year 2011, under both scenarios.

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.

3.2. 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.

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. The calibrated values for parameters affecting streamflow in the Mahabad Dam watershed.

Parameter	Description	Subbasins	Calibrated Value
DEEPST.gw	Initial depth of water in the deep aquifer (mm)	All	7566.89
ALPHA_BF.gw	Baseflow alpha factor (1/days)	All	0.87
REVAPMN.gw	Threshold depth of water in the shallow aquifer for “revap” or percolation to the deep aquifer to occur (mm)	All	702.66
SNO_SUB.sub	Initial snow water content (mm)	All	28.308
CH_K1.sub	Effective hydraulic conductivity in tributary channel alluvium (mm/h)	11, 13, 14, 21, 24–26, 28, 36	215.298
		15–20, 22, 23, 27, 29–35, 37–45	5.06
CH_N1.sub	Manning’s “n” value for the tributary channels	11, 13, 14, 21, 24–26, 28, 36	26.46
		15–20, 22, 23, 27, 29–35, 37–45	17.279
CH_N2.rte	Manning’s “n” value for the main channel	11, 13, 14, 21, 24–26, 28, 36	0.235
		15–20, 22, 23, 27, 29–35, 37–45	0.116
CH_K2.rte	Effective hydraulic conductivity in main channel alluvium (mm/h)	11, 13, 14, 21, 24–26, 28, 36	160.026
		15–20, 22, 23, 27, 29–35, 37–45	371.048
EPCO.hru	Plant uptake compensation factor	All	0.825
CANMX.hru	Maximum canopy storage (mm)	All	3.417
OV_N.hru	Manning’s “n” value for the overland flow	11, 13, 14, 21, 24–26, 28, 36	10.78
		15–20, 22, 23, 27, 29–35, 37–45	22.0
SOL_AWC(1).sol____TACONIC	Available water capacity of the soil layer (mm H2O/mm soil)	All	0.12
SOL_K(1).sol____TACONIC	Saturated hydraulic conductivity (mm/h)	All	15.41
SOL_BD(1).sol____TACONIC	Moist bulk density (g/cm3)	All	0.9
SOL_ZMX.sol____TACONIC	Maximum rooting depth of soil profile (mm)	All	1124.15
SOL_AWC(1).sol____BENSON	Available water capacity of the soil layer (mm H2O/mm soil)	All	0.21
SOL_K(1).sol____BENSON	Saturated hydraulic conductivity (mm/h)	All	19.91
SOL_BD(1).sol____BENSON	Moist bulk density (g/cm3)	All	1.26
SOL_ZMX.sol____BENSON	Maximum rooting depth of soil profile (mm)	All	2036.35

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.

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.

3.3. TP Concentration in the Reservoir

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. The calibrated values for the parameters affecting the total phosphorus (TP) concentration in the reservoir.

Parameter	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
TP settling rate (×10−6 m/s)	0.35	0.59	1.12	3.54	2.8	2.06	2.31	0.63	3.39	2.39	0.94	1.52
TP resuspension rate (×10−7 m/s)	0	0	0.5	0	0	0	0	2.2	0.5	0.22	0.18	0.24
TP burial rate (×10−8 m/s)	1	1	1	1	1	1	1	0	0	0	1	1
Initial TP load in the reservoir (ton)	4.85
Initial TP load in sediments (ton)	100

presents the parameters and their calibrated values used to calibrate the TP concentration in the reservoir.

Table 6. The average annual TP concentration in the reservoir under climate change and TP loading scenarios.

Scenario	Upstream (%TP Change)			Downstream (%Water Allocations Change)				TP Concentration (μg/L)	% Change
	Agriculture	Livestock Farming	Residential	Agriculture	Industry	Domestic	Other
RCP4.5 (Inflow = −5.1%, Evaporation = +19%)
1			↑			↑		88.5–90.4	(+5.2)–(+7.5)
2	↑	↑		↑				93.2–109.3	(+10.8)–(+29.9)
3					↑			87.8–88.2	(+4.4)–(+4.8)
4				↓	↓	↓	↓	83.9–86.1	(−0.3)–(+2.3)
5	↓	↓	↓					73.0–82.8	(−13.2)–(−1.6)
6	↓	↓	↓	↓	↓	↓	↓	69.4–81.3	(−17.5)–(−3.4)
RCP8.5 (Inflow = −3.0%, Evaporation = +23%)
1			↑			↑		87.3–89.1	(+3.8)–(+5.9)
2	↑	↑		↑				91.8–106.2	(+9.1)–(+26.2)
3					↑			86.6–87.0	(+2.9)–(+3.4)
4				↓	↓	↓	↓	82.9–85.1	(−1.5)–(+1.2)
5	↓	↓	↓					72.1–81.7	(−14.3)–(−2.9)
6	↓	↓	↓	↓	↓	↓	↓	68.6–80.3	(−18.5)–(−4.6)
↑: +(10–30); ↓: −(10–30)

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.

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.

4. 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.