Simulation of surface runoff in Karaj dam basin, Iran

A major part of the Iranian capital drinking water is supplied from Karaj reservoir 100 km northwest of Tehran. This reservoir collects water from 849 km2-catchment which is undergoing accelerated changes due to deforestation and urbanization. The main objective of this study is to develop a catchment modeling platform which translates ongoing land-use changes, soil data, precipitation and evaporation into surface runoff of the river discharging into the reservoir: Soil and Water Assessment Tool, SWAT, model along with hydro-meteorological records of 1997–2011. A variety of statistical indices were used to evaluate the simulation results for both calibration and validation periods; among them, the robust Nash–Sutcliffe coefficients were found to be 0.58 and 0.62 in the calibration and validation periods, respectively. This project has developed a reliable modeling platform with the benchmark land physical conditions of the Karaj dam basin.


Introduction
More detail information on the status of rainfall runoff increases decision-making power on future programs for watershed managers of natural resources for sustainable development. Recently, rainfall-runoff models are widely used by hydrologists to simulate watersheds runoff and play a key role in water resources management (Bilondi et al. 2013). On the other hand, several programs and techniques have been developed to reduce uncertainty of parameters to get best fit of parameters in the hydrological modeling (Singh and Dutta 2016). These hydrological models are utilized in assessing the land-use change, effects of climate change, water quantity and quality of river flow, forecasting of flood, development and protection of water resources, integrated management of water components and most of the water resource management measures Woolhiser 2002, 2016;Wurbs 1998). In most of the basins that require water resource planning, there is a lack of hydrometric stations for surface runoff measurement, or the data of measurement stations are incomplete; it is also unlikely to have measurement stations in all regions in the near future. Hydrological models, therefore, can evaluate the obtained runoff from precipitation in basins with no or incomplete surface runoff data with minimum cost and time by simulating the rainfall-runoff process.
Soil and Water Assessment Tools (SWAT) is one of the most popular hydrological models for runoff simulation in the world. The SWAT model is a continuous and distributed hydrological model to perform at the watershed scale (Hosseini and Ashraf 2015; Krysanova and Srinivasan 2015). This model introduced by Arnold et al. (1998) and further developed by Arnold and Fohrer (2005) can simulate hydrological processes in basins of small, medium and large scales. This model can be run on monthly, daily and hourly time scales and can simulate all components of water balance (Arnold et al. 1998).
SWAT model was selected among fifteen hydrological models to separate water components, which successfully fulfilled the annual and monthly water budget estimation as well as suspended sediment yield in the Taleghan catchment during 1987-2007, with a reliable accuracy (Hosseini 2010). The SWAT-CUP (Abbaspour et al. 2007b) is a computer program that links the sequential uncertainty fitting (SUFI-2) algorithm to SWAT model. Up to now, researchers used SUFI-2 algorithm for model calibration and uncertainty analysis of parameters of SWAT model as a popular method in the world. Fukunaga et al. (2015) investigated the application of the SWAT hydrologic model in a tropical watershed at Brazil. Narsimlu et al. (2015) in the Kunwari river basin applied SUFI-2 algorithm for model calibration, sensitivity and uncertainty analysis for 19 years from 1987 to 2005. Nyeko (2015) assessed the capabilities and limits of SWAT model in modeling watershed that has limited field and hydrologic data for possible use in water resources management. Krysanova and Srinivasan (2015) assessed five projects of different applications of SWAT covering the impacts of land cover change, climate change and combined impacts of climate change and human activities in water management. Vilaysanea et al. (2015) applied SWAT model to test the capability of the model for predicting stream flow and used SUFI-2 algorithm for calibration and uncertainty analysis in Xedone river basin. Lin et al. (2015) in their study found out the effects of land-use and land cover changes on runoff response using SWAT model. They used two different land-use scenarios (1985 and 2006, with reduced forest and increased cropland and urbanized area) in Jinjiang catchment. Hosseini et al. (2012) studied the effects of land-use changes on the water balance components in Taleghan basin using SWAT model. Both SUFI-2 and Parasol programs were used to calibrate and verify the model using the data from the period 1987-2007. A comparison of the water budget components performed by this model in Taleghan catchment from 1995 to 2004 using the three land-use maps of 1987, 2001by Hosseini and Ashraf (2015 showed reliable results for both calibration and validation periods.
The main objective of this study is to perform possibility of the SWAT model and SUFI-2 algorithm for runoff simulation of the study area, which will contribute to the reservation of natural resources in the Karaj dam basin.

Study area
Karaj dam basin is located in the southern slopes of the Mountain Alborz between the geographical coordinate's longitude 51°E-51°21′E and latitudes 35°30′N-36°6′N (Fig. 1). The area of this basin is 84,894 hectares and average annual water flow for Karaj Dam is 472 million cubic meters. The average daily discharge is estimated 154.54 m 3 s −1 . There are ten hydrometry stations to measure discharge in the study area. The Sira hydrometry station, located above the Karaj dam, selected to runoff simulation.

Description the SWAT model
Soil and Water Assessment Tools (SWAT) is a physically based model for assessing the impact of management and climate on water supplies, sediment and agricultural chemical yields in catchments (Narsimlu et al. 2015). This model is a basin-scale, spatially distributed watershed delivery model developed by the Agricultural Research Service (ARS) at the US Department of Agriculture (USDA). Its purpose is to simulate water, sediment and chemical yields on large river basins and possible impacts of land use, climate changes and watershed management. Output from SWAT can be daily, monthly or yearly, but in all cases is based on a daily model time step. SWAT can be applied in watersheds up to several 1000s of km 2 , using a two-level disaggregation scheme. Preliminary sub-basin identifications are carried out based on topographic criteria, followed by further discretization using land use and soil type. The physical properties inside each sub-basin are then aggregated with no spatial significance.
In SWAT, a catchment is divided into multiple sub-catchments with hydrologic response units (HRUs) that consist of homogeneous land use, management, topographic, and soil characteristics (Abbaspour et al. 2007a). SWAT-CUP is a computer program that provided sensitivity analysis, calibration, validation, and uncertainty analysis of SWAT models (Abbaspour et al. 2007b).
The hydrology model is based in a water balance equation, comprising surface runoff, precipitation, evapotranspiration, infiltration and sub-surface runoff. Potential evapotranspiration is estimated using either the Priestley or Taylor (1972) or the Penman-Monteith (Monteith 1965) methods.
The soil profile is represented by up to 10 soil layers, a shallow aquifer and a deep aquifer. When the field capacity in one layer is exceeded, the water is routed to the lower soil layer. If this layer is already saturated, lateral flow occurs. From the bottom soil layer, percolation goes into the shallow and deep aquifers. Water reaching the deep aquifer is lost, but in turn flow from the shallow aquifer due to deep aquifer saturation that is added directly to the sub-basin channel.
Surface runoff is computed by the SCS curve number method (Arnold et al. 1990) and is therefore a non-linear function of precipitation and a retention coefficient. Once all hydrological processes are calculated for a homogeneous part of the sub-basin, the resulting flows are considered to contribute directly to the main channel. SWAT includes a routing module based on the Routing Outputs to Outlet model (ROTO; Arnold et al. 1990), which takes into account the connectivity of the river network.

Description the SUFI-2 algorithm
A common factor in almost all the published papers is the calibration/validation and uncertainty analysis of the models (Abbaspour et al. 2018). In SUFI-2, parameters uncertainty accounts for all sources of uncertainties. These sources include variables (e.g., rainfall), the conceptual model, model parameters and measured data. To evaluate such uncertainties, SUFI-2 offers two criterion factors, the P-factor and the R-factor. The P-factor indicates the percentage of measured data bracketed by the 95% prediction uncertainty (95PPU), whereas the R-factor calculates the average thickness of the 95PPUb and divided by the standard deviation of the measured data. Theoretically, the value of the P-factor ranges from 0 to 100%, while that of the R-factor ranges from 0 to infinity. A P-factor of 1 and R-factor of zero indicate a simulation that exactly complies with measured data.

Input data
Meteorological data with topography, soil and land-use maps are prepared to simulate model. Digital elevation model with resolution 90 meters is used to create topographic data (Fig. 2). Twenty-one hydrologic sub basins were identified in study area. SWAT model needs soil-type and land-use properties to simulate hydrological parameters within each sub-basin (Winchell et al. 2009).
The volume of surface runoff is simulated using the curve number (CN) method. This method was created by Soil Conservation Service (USDA 1972) which calculates surface runoff as a function of soil type, slope, initial soil humidity and management practices (Neitsch et al. 2011). The evapotranspiration potential is determined using the Penman-Monteith equation and is corrected for the soil cover and simulated plant growth to get the real evapotranspiration rate (Neitsch et al. 2011;Monteith 1965).
The soil type classified into four classes by the FAO (1997) is presented in Fig. 3. The land-use map provided by SCWMRI (2002) is illustrated in Fig. 4.
The climatic data including daily precipitation, maximum and minimum temperature, relative humidity, solar radiation and wind speed were obtained from the Karaj meteorological station. The data used were for the period January 1, 1997-December 30, 2011.

Performance evaluation of the model
The performance of the model was evaluated by visual and statistical comparison of the measured and simulated data. The graphical technique provided an initial general overview (American Society of Civil Engineers 1993). Interpretation of the hydrograph first focused on the peak flows and then on the base flow. There are no existing standards describing the range of values of the statistical parameters that would show acceptable performance of the model (Loague and Green 1991). In this study, seven criteria including R 2 , NSE, PBIAS, RMSE, RSR, P-factor and R-factor are applied to analyze the model (Table 1). American Society of Civil Engineers (1993) has highlighted the Nash-Sutcliffe efficiency coefficient (NSE) (Nash and Sutcliffe 1970) among the various statistical parameters to evaluate hydrological models.
The P-factor (percentage of measured data bracketed by the 95% prediction boundary) is often named 95PPU (Percentage Prediction Uncertainty). The 95PPU is calculated at the 2.5% and 97.5% levels of the cumulative distribution of an output variable obtained through Latin hypercube sampling (Abbaspour 2011). The range of the P-factor varies from 0 to 1, with values close to 1 indicating good fitness between simulated and observed values (Yang et al. 2008a, b).
Another measure quantifying the strength of a calibration and uncertainty analysis is the R-factor, which is the average thickness of the 95PPU band divided by the standard deviation of the measured data. The calibrated parameter ranges can be generated with an acceptable value of the R-factor and P-factor.

Sensitivity analysis of the model parameters
The sensitivity analysis procedure employed measurement data for the period January 1, 1997-December 30, 2011, to evaluate the fit between the measured and modeled time series data. This enabled the identification of the parameters that were influenced by the characteristics of the hydrographic basin, and those to which the model was most sensitive. Evaluation was then made of the way in which adjusting the value of a parameter affected the model output, in order to identify parameters that might improve the characteristics of the model (Veith and Ghebremichael 2009). Sensitive parameters were calculated using the SUFI-2 algorithm and SSE (sum of square errors) equation between the measured and simulated daily runoff data.
The parameters used for the flow were selected based on the experience in the study area. The initial simulation to determine the sensitivity of the model to different parameters was performed using default parameter values. The values were then varied within upper and lower limits established according to the characteristics of each parameter, using three methods. In the first procedure, the initial value of the parameter is modified by adding an increment. The second method consists of multiplying the initial value by a set amount. In the third method, the initial value is substituted by a different value (Van Griensven et al. 2007).
The sensitivity of the model to a parameter is determined using the percentage difference between the output values of the objective function for simulations performed immediately before and after changing the value of a parameter (Veith and Ghebremichael 2009). A higher t-stat value indicates greater sensitivity for a given parameter.
Thirty-five parameters selected for sensitivity analysis that eight parameters were more sensitive for simulation of the flow (Table 2). These parameters were CN2, ALPHA_ BF, CH_K2, ESCO, CH_N2, REVAPMN, GW_REVAP and SOL_BD.

Calibration and validation of the parameters
This research passed through three stages which are as follows: (1) the setup (also known as warm-up) period from 1997 till the end of the year 1999 (3 years), (2) the  (5 years). Usually in modeling, few initial years are selected for setup or warm-up of the model. In general, one-third and two-thirds of the original data sets are selected for model calibration and validation, respectively. Visual analysis of hydrographs shows overestimating of peak runoff with slightly wide uncertainty band (95PPU). Simulated runoff closed to observe one by selecting proper range of sensitive parameters.
The results of calibration period were obtained with four iterations and 500 runs during 2000-2006 (Fig. 5). The calibrated ranges of parameters were performed by model during validation period from 2007 to 2011 (Fig. 6). The results demonstrated that the performance of the calibrated and validated model was satisfactory, with NSE > 0.75, R 2 > 0.5, PBIAS < ± 10, RMSE ≤ (SD/2), and 0.06 ≤ RSR ≤ 0.70. Main statistical criteria in calibration periods show R 2 , P-factor, R-factor and NSE are 0.58, 0.36, 0.36 and 0.57, respectively. This analysis shows 0.62, 0.63, 0.49 and 0.59 in validation period, respectively (Table 3). Coefficient of determination (R 2 ) and Nash-Sutcliffe efficiency coefficient (NSE) show a reliable correlation between observed and simulated values in both calibration and validation periods.
Statistical analyses are used to evaluate the simulation accuracy of the model performance. The observed and estimated mean monthly discharge at Sira hydrometry station in calibration and validation periods is plotted in Fig. 7a, b.
The results obtained for the statistical criteria demonstrated that the model was able to describe the hydrological processes in the watershed of the Karaj River.

Conclusion
Many urban and rural societies within Karaj River catchment are dependent on surface water resources. The longterm fluctuation in surface flow may have serious impacts on the daily lives of upstream residents as well as the reservoir and the capital residents in both terms of flooding and yield. SWAT model reproduced surface runoff in Karaj River reliably for both calibration and validation periods using the entire records of years 1997-2011. Seven main statistical criteria including R 2 , NSE, PBIAS, RMSE, RSR, P-factor and R-factor are used to evaluate the modeling results with the observed records of stream flow with minimum uncertainty. The use of sensitivity analysis enabled the identification of the most important parameters required to model the hydrological processes in the Karaj   Steel and Torrie (1960) 0.5 ≤ R 2 ≤ 1 Satisfactory Nash-Sutcliffe efficiency      River catchment, hence reducing the quantity of model parameters that needed to be calibrated. In this study, eight sensitive parameters including CN2, ALPHA_BF, CH_K2, ESCO, CH_N2, REVAPMN, GW-REVAP and SOL-BD were used in a successful effort to decrease uncertainty. Using the developed modeling and database platform, planners and decision makers will have access to a set of reliable tools to test any proposed scenarios of development within the different gauged or ungauged sub-catchments of the main basin and evaluate the impacts on surface water.