Assessment of the Spatiotemporal Effects of Land Use Changes on Runoff and Nitrate Loads in the Talar River

: This research surveyed the effects of land use changes on ﬂow nitrate pollution in the Talar River (northern Iran), using Landsat images of 1991 and 2013 and the Soil and Water Assessment Tool (SWAT). The results indicated that forest areas decreased by 14.9% and irrigated crops, dry land farming areas, range lands and residential areas increased by 46.8%, 31.1%, 4.7% and 17.5%, respectively. To calibrate and validate the studied period, the Nash Sutcliffe model efﬁciency (NSE) and coefﬁcient of determination ( R 2 ) were applied, ranging from 0.57 to 0.75 and from 0.62 to 0.76 for ﬂow simulation and 0.84 and 0.63 and 0.75 and 0.83 for nitrate simulation, respectively. The results of land use scenarios indicated that respective water ﬂow and nitrate loads increased by 34.4% and 42.2% in 1991–2013 and may even increase by 42.3% and 55.9% in the simulated period of 2013–2050 in all sub-basins. It is likely that the main reason for these results was due to the increase in agricultural activities and the decrease in forestry areas. Our ﬁndings showed the useful combination of modelling techniques (land cover changes and SWAT) to develop valuable maps able to design correct land management plans and nature-based solutions for water quality of runoff water harvesting systems in the future.


Introduction
Humans condition and modify natural ecosystems, which may lead to negative consequences in fluvial ecosystems [1][2][3]. Water is the most important source and its quality directly affects human health, soil fertility, biogeochemical cycles, welfare, food security and industrial development [4,5].
Recent studies claim that land use changes are one of the most important factors that affect surface water characteristics and soil fertility [6][7][8]. During the last decades, land use changes and water quality of rainwater harvesting systems have been associated with population growth, urban development, intensive agricultural activities, degradation of pastures and deforestation [9,10]. Although these factors are causing an increase in land degradation processes, little attention is paid to water resources and ecosystem management plan elaborations [11,12]. Drastic land use changes, which are not well-planed, are able to change water quality and quantity and subsequently, modify water balance in arid and semi-arid watersheds, increasing soil erosion and nutrient losses [13][14][15].
Several methods can be used to quantify the hydrological status after land use changes, but there is no consensus to standardize one method for regional or global contexts [16,17]. For example, the paired-watershed method is commonly used for small watersheds to detect hydrological differences among similar basins with different land uses [18,19]. Another widely used method is time series analysis, which is a statistical method that requires a long-term data base, which is often not available in non-developed countries [20]. Therefore, it is necessary that in order to get the maximum information from limited data, more comprehensive tools are developed [21]. In this way, hydrological models are tools that may be able to simulate spatiotemporal variations of hydrological processes and enhance the understanding of land use changes [22].
Several models have been developed to simulate hydrological processes for watershed areas such as ANSWERS (Areal Nonpoint Source Watershed Environment Response Simulation) [23], AGNPS (Agricultural Non-Point Source Pollution Model) [24], HSPF (Hydrological Simulation Program-FORTRAN) [25], WEEP (Water Erosion Prediction Project) [26] and SWRRB (Simulator for Water Resources in Rural Basins) [27]. One of the most commonly used models applied is the Soil and Water Assessment Tool (SWAT), developed by Arnold et al. [28]. Many previous studies have demonstrated the ability of SWAT to detect the negative impacts of land use changes on hydrology and water quality in different landscapes under distinct environmental conditions [29,30]. For example, Chaplot et al. [31] assessed the impact of land use changes on water discharge, sediment transport and NO3-N load at the outlet of the central Iowa watershed. Also, Kigira et al. [32] evaluated the impact of land use changes on water and sediment yield in the Thika River catchment in combination with satellite images. Lam et al. [33], using SWAT, investigated the effect of the best management practices (BMP) such as buffer strips in extensive land use managements such as grazing in the Kielstau catchment (north Germany). In China, Ouyang et al. [34], from 1976 to 2006, evaluated the effect of agricultural development on non-point source (NPS) pollution. In these studies, the models were calibrated and validated using daily data of flow, sediment, NO 3 -N and total phosphorus. Li et al. [35] classified three categories of factors for runoff using SWAT: climate, land cover change and direct human activities. Their result showed that climatic variations were the main cause of runoff declines over the time in the entire basin. Giri et al. [36] implemented ten best management practices (BMP) in agricultural areas in the Saginaw River Watershed using SWAT based on four targeting methods (Load per Subbasin Area Index -LPSAI-, Load per Unit Area Index -LPUAI-, Concentration Impact Index -CII-, and Load Impact Index -LII-). They reported strong relationships between the LPSAI and LPUAI methods both for the sediment and total nitrogen targeting scenarios. However, the SWAT model, even in watersheds with limited hydrological data, is able to show a satisfactory performance, if it is well calibrated and validated [37]. It is important to remark that this rapid development in geostatistical techniques has been possible due to the intensive improvement of remote sensing (RS) techniques and geographic information system (GIS) databases [38][39][40]. Thus, there is no point to obviate the discussion about the impact of rainfall scenarios and land-use change on hydrologic responses in degraded catchments in order to design correct management policies and implications of the population [41][42][43].
In Iran, a decrease in forest areas has become one of the critical problems in recent years, and the water quality of rainwater harvesting systems is being affected. One of the most important examples is the Hyrcanian forest, which contains one of the most important and significant natural habitats for in-situ conservation of biological diversity, including species from the Jurassic era such as unique endangered forest species, mountain rangeland with rare plant species, and endangered animal species. Unfortunately, in two decades this unique forest area has decreased, with the Talar watershed being one of the most affected areas due to fuel, timber smuggling, livestock grazing, road development, exploitation of mines and construction of factories. However, little attention has been paid by the scientific community to quantify this dynamic to help rural areas.
Because accurate and long-term water quality samples are not easily available in Iran and the establishment of new water quality monitoring stations and data collection are very costly and time-consuming, our main goal is to apply a hydrological model to assess the relationship between land use changes and their effects on the water quality of rainwater harvesting systems (water flow and nitrate transport) in the degraded areas of Hyrcanian forest. We intend to provide scientific evidence to promote sustainable landscape and water management plans for rural and natural areas. Therefore, to achieve this goal, we calibrate, validate and conduct a sensitivity analysis of the SWAT model in combination with remote sensing data. A gauging station in the Talar River and Landsat images were used to evaluate land use changes between 1991 and 2013. Finally, a hypothetical scenario for 2050 based on annual changes between 1991 and 2013 was also predicted.

Study Area
Talar River is a mountainous watershed with a total area of 210,088.7 ha, which is located between the eastern longitudes 52 •  , and this region is characterized by a Mediterranean rainfall regime. The region has a mean annual rainfall of 552.7 mm and annual minimum and maximum temperature of 7.7 and 21.1 • C, respectively [44]. Figure 1 shows the location of Talar, which is also characterized by gentle steep hillslopes (less than 12%) and is mainly distributed between the southeast highland plains and land near the outlet of the watershed. The sloping areas (more than 60%) occupy 14.2% of the whole territory. The remaining 38% have inclinations greater than 30%. The lithology and stratigraphy of Talar watershed are diverse. Most of the area is covered by sedimentary, igneous and continental rocks [44]. Most of the rock units (about 61%) are related to the Mesozoic Era. Moreover, 58 types of soils with different physical and chemical properties have been registered. The average discharge for Shirgah station located at 52 • 53 10 eastern longitude and 36 • 57 10 northern latitude, registered 7.95 m 3 /s and the highest discharge was 93.46 m 3 /s . The most important land uses of the studied area are forests, dry land farming, irrigated lands, rangelands and residential areas [44].

SWAT Model
The SWAT model is a distributed, physics-based and continuous model that is used to perform simulations at the catchment scale. This model was developed to predict the impact of land use and land management practices on water, sediment and pollutant yields in large watersheds using a daily time step [28]. In the SWAT model, the main watershed is divided into multiple sub-watersheds, based on soil types, land uses and slope classes, resulting in hydrological response units (HRU).
Soil water content, surface runoff, nutrient cycles, sediment yield, crop growth and management practices can be also simulated in each HRU [45].
The simulation of the hydrological cycle allows separation into a land phase and a water phase. The simulation of the land phase is based on water balance equations. The generated runoff in each HRU is summed to calculate the water reaching the main channel in each sub-basin. The water phase of the hydrological cycle is defined as the process of routing runoff in a river channel using the variable storage coefficient method. A detailed description of the model can be also found in Neitsch et al. [46].
SWAT is also able to simulate evapotranspiration, infiltration, percolation, runoff generation, nutrient cycling and transport for each HRU. Water and sediment routing as well as in-stream nutrient processes can be simulated along the channel length for each sub-basin [46]. The algorithms for nitrogen cycling and transport are based on the EPIC (Erosion-Productivity Impact Calculator) model [45]. Net mineralization is simulated with one active and one stable organic nitrogen pool. Plant uptake of nitrogen is estimated using a supply and demand approach. Finally, nitrates in soils can be removed from the soil via denitrification, mass flow of water as well as plant uptake [46].

Data Collection
In this study, SWAT version 2012 was used to simulate water quality and quantity affected by land use changes. The input data for the SWAT model for this study were weather data, a digital elevation model (DEM), soil and land-use, streamflow, and water quality data. A detailed description of these input data and their sources are presented in Table 1.
The land use data were generated from Landsat satellite images (Landsat Thematic Mapper (TM) for 1991 and Landsat Enhanced Thematic Mapper Plus (ETM+) for 2006, which were downloaded from the USGS website ( Figure 2). Both images correspond to June. Land-use maps were generated using a supervised classification based on the maximum likelihood algorithm. Then, the overall accuracy and kappa coefficient (κ) were used to assess the classification accuracy. Land use data were classified into five different land use types including forest, irrigated land, dry land farming, rangeland and urban areas. To achieve these goals, the software ENVI Version 5.1 (Boulder, CO, USA) was used. In this study, for determining the land use scenario for 2050, parameters such as elevation classes, slope, distance from the river and roads as well as residential areas were used as indicators of land use change. After the image classification was completed for 1991 and 2006, the differences among land use classes were calculated between the two paired-periods. Then, by dividing this difference over a period of 15 years (1991-2006), the annual rate of change was obtained for each year. Finally, the land use polygons in the GIS environment were reduced by considering parameters.  The watershed was delineated into 23 sub-watersheds using a digital elevation model. Then, the area based on special conditions of soil types, land uses and topography was divided into HRU units. Daily weather data (precipitation, maximum and minimum temperature, relative humidity, and wind speed and sunshine hours) were obtained from three local weather gauges from

Model Performance Evaluation
Due to limitations related to the input data such as slope, elevation, soil type, land cover and complexity of hydrologic processes, the calibration of SWAT output was considered as complex. There are many statistical methods to access and evaluate the accuracy of model results. In this study, the accuracy of SWAT simulation was evaluated with the observed values using the coefficient of correlation (R 2 ), Nash-Sutcliffe coefficient (ENS) and percent bias (PBIAS).

Coefficient of Determination (R 2 )
This method is considered a good indicator to determine the data uniformity observed and simulated. R 2 ranges from 0 to 1, where higher values indicate less error variance [47] and values greater than 0.50 are considered acceptable [48]. R 2 was calculated using Equation (1):

Nash-Sutcliffe Efficiency (NSE)
NSE is a normalized statistical method that assesses the relative magnitude of the residual variance compared to measured data variance [49]. NSE was calculated using Equation (2): NSE ranges between −∞ and 1.0; a value of 1 is the optimal value. Values between 0.0 and 1.0 are generally viewed as acceptable levels of performance. Generally, the model simulation was considered satisfactory if NES > 0.5 [47].

Percent Bias (PBIAS)
PBIAS allows calculation of the average tendency of the simulated values to be larger or smaller than their observed values [50]. A model with low-magnitude values from PBIAS indicates accurate simulation (optimal value is 0). Positive and negative values of PBIAS indicated that the model registered an underestimation or overestimation bias, respectively [50]. PBIAS was obtained as follows from Equation (3): where PBIAS is the deviation of data being evaluated, expressed as a percentage. Table 2 shows the classification of statistical results of NSE and PBIAS for flow and nutrients [47]. Table 2. General performance for recommended statistics at a monthly time step [47].

Sensitivity Analysis
In this study, Swat Cup software and the SUFI-2 program were used to conduct the sensitivity analysis, calibration and validation. Many parameters are able to affect the model results; therefore, it is necessary to specify the key parameters that will register more impact. In our case, 24 parameters were selected using monthly observational discharge (2001)(2002)(2003) and nitrate (2012-2013) data. Sensitivity analysis was automatically performed with 500 iterations and 24 sensitivity parameters.

Accuracy Assessment of Land Use Classification
It is well known that the accuracy assessment is the most important part of remote sensing classification and land use change classifications. Lu et al. [51] stated that the most common accuracy assessment elements include overall accuracy, producer's accuracy, user's accuracy and the Kappa coefficient. In this study, the Kappa coefficient and overall accuracy were also used to assess the verification of land use classification (Tables 3 and 4).  Forest, irrigated land, dry land farming and residential areas in 1991 were accurately classified using the error matrix (Table 3). Residential areas and irrigated land are sometimes mistaken for range land. This may be due to the proximity of these land uses to rangeland areas. The classification error matrix in 2006 shows that dry land farming was accurately classified, followed by forest, residential areas, irrigated land and range lands (Table 4). In 1991 and 2013, the classification efficiency reached 96.7% and 87.2%, respectively, and the corresponding kappa coefficient was 0.94 and 0.80.

Land Use Changes from 1991 to 2013
A total of five land use classes were identified from the satellite images, including forest, irrigated land, dry land farming, range land and residential areas.
In Table 5, land use changes are summarized. Land use maps obtained using Landsat TM and ETM images of Talar River for 1999 and 2013 are displayed in Figure 4.  It can be observed that the dominant land uses are forest and range lands with total areas of about 80% ( Table 5). The 1991 land use map shows that the forest, irrigated land, dry land farming, range land and residential areas had respective areas of 40.8%, 0.5%, 11.3%, 46.6% and 0.8% (Figure 3a), but the land use map shows that forest, irrigated land, dry land farming, range land and residential areas had areas of 34.8%, 0.7%, 14.8%, 48.8% and 0.9% in 2013 (Figure 3b). From 1991 to 2013, forest decreased from 40.8 to 34.8% but other land uses increased (Table 5 and Figure 3). This dynamic has been reported in several areas around the world such as in the Mediterranean landscapes or tropical territories, showing critical land degradation processes such as soil erosion or a decrease in soil fertility after these changes [52][53][54]. Therefore, land planners and policy makers should detect these drastic land use changes and apply restoration measures or strict plans to protect forestry areas and rural areas that use water for traditional agricultural activities [55,56].

Land-Use Scenario in Year 2050
During the last twenty years, a loss of forestry areas has been registered. Thus, we try to evaluate the impact of land use changes on water flow and nitrate transport for a future scenario. To achieve this goal, three land use scenarios are chosen: 1991 (S1), 2013 (S2) and 2050 (S3). The image classification is presented in Table 6. Figure 5 shows land-use scenarios for the year 2050. As other authors observed, parameters such as slope, distance from roads and residential areas affected land use polygon size [57]. We can observe that the outlook for future water quality of rainwater harvesting systems is not so hopeful. Our results show that a decrease in natural areas and an increase in human expansion would be higher than nowadays. It is important to remark that several authors claimed that the loss of forestry areas is affecting the levels of carbon, ozone and biodiversity [58,59]. Moreover, the risk of floods and slides may also increase, as demonstrated other authors in other Iranian areas [60][61][62].

Sensitivity Analysis
Because of different characteristics of watersheds such as soil physical and chemical properties, land use and physiographic characteristics, the effects of parameters can vary with respect to hydrology, sediment and nutrient discharges [63,64]. Hence, parameter sensitivity analysis must locally performed although none of the sensitivity analysis methods is clearly superior [65].
Since water quality models have comprehensive parameters and there are various data to compare with simulation results (e.g., flow, suspended sediment, nitrogen and phosphorus), sensitivity analysis methods are needed in order to adjust parameters considering the model outputs [66]. In this study, the selected parameters for sensitivity analysis based on relevant literature [67] and SWAT documentation [27,46] are presented in Table 7. The soil conservation service curve number (SCS-CN) was considered in this study as the most sensitive parameter for the three stations and the main parameter for flow simulation. The other most sensitive parameters are SMTMP, SOL_BD, ALPHA_BNK, SOL_K, SOL_AWC, SMFMN, SLSUBBSN and REVAPMN. The last-ranked flow parameters are LAT_TTIME, CANMX, TIMP, CH_N2, GWQMN, SURLAG, TLAPS, ESCO, GW_REVAP and ALPHA_BF, (meaning in Table 7). Also, for nitrogen discharge, the most sensitive parameters are CN2, and GWQMN and ALPHA_BF and also water quality parameters SDNCO, ERORGN, SOL_NO3 and SOL_ORGN (meaning in Table 7).
Lenhart et al. [68] and Zhai et al. [69] reported that CN2 (SCS runoff curve number for moisture condition II) was also the most sensitive parameters in the simulation of flow and water quality. They explained that this is because the land use and soil types often have a great influence on the CN2 value [70]. As a consequence, this parameter has a major influence on the composition of the water balance and the simulated flow [71].

Stream Flow Simulations
The model simulation was evaluated for both the calibration and validation periods using graphical comparisons and three statistical measures (Table 8): R 2 , NSE and PBIAS. Prior to model calibration, the mean of the simulated stream flow data for Shirgah, Kerikkola and Polsefid relative to the mean observed data was approximately 1.4, 2.6 and 0.02%, respectively. For calibration (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) and validation (2009-2013) periods, NSE and R 2 ranged from 0.57 to 0.75 and 0.62 to 0.76, respectively, for the three monitored stations. The monthly flow simulations show a good calibration and validation based on the performance ratings of Moriasi et al. [47] with rates of 0.65 < NSE < 0.75 and ±10 < PBIAS < ±15. The PBIAS coefficient for calibration and validation period ranges from 1.3% to 13%. These results show an underestimation of the simulated mean.
In Figure 6, hydrographs of observed and simulated monthly flow at Shirgah, Kerikkola and Polseffid stations during 2001-2008 (calibration period) and 2009-2013 (validation period) are presented. We can appreciate that time and the volume of flow peaks are simulated with high accuracy, but the maximum flow volume of most peaks is underestimated. In general, the model statistical evaluations indicate a good model performance for flow simulation in Talar watershed. Our findings agree with those presented by Cao et al. [72], who, using SWAT, also noted that the number and location of rainfall gauges strongly affected the accuracy of the flow simulation. Wu and Chen [73] also detected similar results, but they had better coefficients in their studies. They obtained an NSE and R 2 of 0.75 and 0.80 for the calibration period, and 0.74 and 0.87 for the validation period using SWAT-cup and SUFI-2, respectively. Moreover, they were able to report reasonable consistency among the simulated and observed runoff values as well as the responses to the precipitation dynamic. In another case, Xu et al. [74] compared climate data measured at Xingshan climate station and gridded climate data for the accuracy of the simulation result. They stated that the gridded dataset had lower accuracy of stream flow predictions, and for this reason, gridded climate data were not used for improving simulation results.

Nitrate Simulations
The concentration of nitrates in water is also considered an important indicator to evaluate the water quality of rainwater harvesting systems. Figure 7 shows the monthly measured and simulated nitrate loads for calibration and validation periods (2012 and 2013).
The model has calibration and validation values close to the measured values for most of the monthly nitrate loads (Figure 7). These results are also confirmed by our statistical measures in Table 8. The model efficiency NSE for calibration and validation periods is 0.74 and 0.62, respectively, and the R 2 coefficient for calibration and validation periods is 0.75 and 0.83, respectively. PBIAS is 5.4 (for calibration period) and −39.3 (for validation period). These results show that the model has a general underestimation (in first period) and a high overestimation (in second period). Similarly, the model performance in simulating nitrate loads for both calibration and validation periods can be considered suitable (with 0.65 > NSE > 0.75, ±10 < PBIAS < ±40) [46] and acceptable in stream flow and nitrate simulations. In this study, the result of nitrate simulation is close to that for other catchments like the Lam Takong River Basin in Thailand [75]. These authors obtained R 2 values higher than 0.6 and PBIAS values lower than 25% for nitrate simulation. Moreover, they established fair relationships between observation and simulation results.  Table 9 shows the percentile changes of discharge and nitrate at the sub-basin scale under three land use change scenarios. For both the first period (1991-2013) and the second period (2013-2050), the discharge increases in most sub-basins because of deforestation and the development of rainfed and irrigated agricultural land. The water discharges do not change in some sub-basins such as 3, 8 and 21 (in the first period) and sub-basins 3 and 13 (in the second period). Additionally, the highest changes in discharge are registered in sub-basins 22, 7 and 2, (in three periods). Nitrate output shows an increase, except in sub-basin 3 and 8 for the S1-S2 scenario. The highest nitrate output observed in the sub-basin is located in sub-basin 13, with a rate of 613.4% (Table 9). In addition, nitrate output increases for almost all sub-basins for the S2-S3 scenario. Here, the highest nitrate rate is calculated for sub-basins 23, 15 and 16, with rates of 1277.6, 1240.54 and 3820.7, respectively. The third scenario (S1-S3) shows a sharp increase in nitrate output in several sub-basins, except 8 and 16. Figures 8 and 9 show the changes in average annual discharge and nitrate yield in all sub-basins for the three land use scenarios. The results also show that forest destruction and the development of dry land farming are causing a higher concentration of nitrates in surface water. In general, the results state that the nitrate loss originates from agricultural lands and rangelands, which can be considered as driving factors of water quality [76]. We conclude that these areas are the main sources of nitrate pollution in Talar watershed, being located in the hotspot of dry land farming, which affects irrigated and forestry areas. The obtained results correspond with the conclusions of Ahearn et al. [77], Behera and Panda [78] and Mohammad and Ahmad [79], who claimed that dry land farming is one of the most important land uses that is able to affect nitrate mobilization. Castillo et al. [80] also reported that nitrate is significantly related to the water quality of rainwater harvesting systems and the elimination of forestry areas in the catchments. In fact, apart from other many factors, forest land and irrigated areas are identified as a sink for NPS pollution [81,82] because decayed organic matter from forest land can produce organic nutrients and cause nutrient loss during rainfall. Table 9. Percent changes in simulated average annual discharge and nitrate transport for three land use scenarios (S1: 1991, S2: 2013 and S3: 2050).

Broader Impact
In recent years, the Mazandaran province has experienced widespread land use changes due to its favorable climatic conditions and fertile soils, population growth and food demand needs. A large amount of forest land has been turned into garden and agriculture land uses. Due to a recent water shortage, a large amount of riverside land is devoted to rice cultivation due to easy access to water. These extensive land use changes in Mazandaran province endanger water and soil resources in watersheds and sustainable development. These results can be extrapolated to other areas with similar environmental problems applying the same methods and increasing the scales. Thus, we consider that this study can show a broader impact and be useful for land managers and decision makers in the future.

Conclusions
In this study, an evaluation of the impact of land use changes on water quality of runoff water harvesting systems of the Talar River watershed, using satellite data and the SWAT model, was conducted. We conclude that the combination of both techniques is able to provide useful information to analyze spatiotemporal land use dynamics and water quality of rainwater harvesting systems. The simulated flow and nitrate loads were in reasonable agreement with the measured values using statistical coefficients such as R 2 , NSE and PBIAS. The results of land use scenarios for all sub-basins indicated that the flow and nitrate increased, with rates of 34.4% and 42.2% in the first period (1991-2013) and 42.3%, and 55.9% in the second period (2013-2050), respectively. The spatial modeling and assessment of water and nitrate yield presented may allow better understanding of the land use change challenges in Talar River watershed. Moreover, the results will allow researchers, managers and policy makers to identify critical areas, appropriate land use management plans and pollution control measures.