Daily Flow Simulation Using Wetspa Model with Emphasize on Soil Erosion (Study Area: The Neka Catchment in Mazandaran Province, Northern Iran)

1. Introduction

Soil erosion by water is one of the most important land degradation processes in Mediterranean environments. This process is strongly linked to problems of flooding and channel management. The relationship between land use and erosion in mountainous forested watersheds has been known in a qualitative sense for some time. Vegetation management, forest road construction, and forest fires, impact basin sediment yield by increasing the amount of sediment available for transport and the amount of surface water available to transport it.For early flood warnings as well to get time for planning and operation of civil protection measures it has become very important that forecasts are made and simulation of floods is carried out. With the ever increasing demand for water resources, it has become very important that the natural processes of floods be predicted, so that current and future environmental issues can be addressed well in time. A simplified representation of the natural hydrological system is the hydrological model. In this model, different physical processes are represented at different time scales and at a wide range of times. This has basically been associated with a lack of appropriate observational data to constrain model states, increase in the number of model outputs [18] and lastly, the complexity of the model. Basically, the distributed hydrological models give us the opportunity to deal with forcing implement these models.The models will help provide the means by which important information regarding existing and future stream flow conditions, very important information regarding hydrological state variables and the state of knowledge on basins of interest can easily be captured for use. Every entity in Iran has suffered great losses due o the floods together with socio-economic development. Among all the different kinds of natural disasters, flood is ranked first in terms of frequency, affected area, losses caused and the severity. On August 10th, 2001, a big flood took place in Golestan and Gorgan River with return period of 200 years and caused a lot of damage. The damages caused by the flood include 15,000 hectares damages to agricultural lands, 10,000 people rendered homeless, 10,000 hectares of damage to forests, and the greatest loss was the human death toll. 247 human beings had been killed in the disaster. The total damage to the entire province was a staggering 491 billion Rials. Once again on July 29th, 1999 Neka city in Mazandaran province was hit by a flood similar to the one which had hit the Golestan and Gorgan Rivers. About a billion dollar worth of damage was caused, including more than 4000 shops and homes damaged to about 50% to 100 %, 400 km railway damaged. 33km of road and about 100 people injured. The Neka river basin is located in northern Iran. It is frequently affected by storms and heavy rain which causes inundation. Flood forecasting modelling is the most important component of the real-time flood forecasting system. This system can mitigate such natural disasters. Flood warning and forecasting systems mostly use hydrologic/hydraulic models. These models, when optimally validated and calibrated can be very effective in minimizing flood damage through non-structural means. In the early years, these flood models were very simple with sophistication in technology comes effectiveness. With advances in geographic information systems and remote sensing theses models have now become more effective. The advantage of these models is that spatially distributed basin characteristics on stream flow can be reflected by these models. There are various studies in which this particular model has been applied, including the Alzette river basin in Luxembourg [9], Barebeek catchment in Belgium [3], the Hornad watershed in Slovakia [2], the Suoimuoi catchment in northwest Vietnam [8], the Simiyu river (Lake Victoria) in Tanzania [16] and the Suriname river basin, in central Suriname [14].

2. Methodology

2.1. WETSPA model

The WETSPA model which was proposed by Wang et al, (1996) [19] and predicts regional or basin Water and Energy Transfer between Soil, Plants and Atmosphere. The unique thing about it is that it has a physical basis and is a distributed hydrological model, laying down the concept of the composition of a basin hydrological system from atmosphere, canopy, root zone, transmission zone, and saturation zone. In order to briefly describe the model, it is vital to mention that heterogeneity in the basin is dealt with by its division into various grid cells which have further divisions for maintaining water and energy balance into a bare soil and vegetated portions. Another feature is that a vertical flow in a single dimension simplifies the movement of water in the soil and is inclusive of surface infiltration percolation and increase in capillary in the area which not saturated and underground water’s recharge. All this can be done by defining the proper and appropriate parameter classes of land use, topography and soil type, and also by using the available data base.

The model for root zone water balance for each grid cell is obtained by Inputs and outputs equations:

$$D\frac{d\theta }{dt}=P-I-S-E-R-R$$ E1

where D (m) is root depth, θ (L³L⁻³] is soil moisture, t (m) is time, I (LT⁻¹) is initial abstraction including interception and depression losses, S (LT⁻¹) is surface runoff or rainfall excess, E (LT⁻¹) is evapotranspiration, R (LT⁻¹) is percolation out of the root zone, and F (ms⁻¹) is interflow. The assessment for excess rainfall is done by means of modification in a moisture-related rational process with a latent runoff coefficient with due consideration to factors such as land cover, slope, soil type, magnitude of rainfall, and pre soil moisture.

$$S=c\left(P-I\right){\left(\frac{\theta }{{\theta }_s}\right)}^{\alpha }$$ E2

where θ_s(L³L⁻³) is soil porosity or saturated water content, c (-) is potential runoff coefficient.A lookup table has been used for deriving the values of C, with associated values to slope, soil type and land use classes [10]. Literature was searched to obtain default values but they may be changed by the used if necessary as suitable according to a region’s specific situation. Liu.et al (2005) [8] did this by creating a look up table for catchments, associating latent rainfall excess coefficient to various arrangements of slope, soil type and land use. This study utilizes these values. Rainfall excess coefficient, given by :

$$\alpha =\max \left[1,K_{run}+\left(\frac{1-K_{run}}{P_{\max }}\right)\right]$$ E3

Where K _run (-) is surface runoff exponent and P _max (LT^-1) is a rainfall intensity scaling factor. Is impacted by rainfall intensity and the effect is reflected by α (−). Α is greater when the rainfall intensity is not high and as a result, surface runoff is lower, and the approach is shifted towards high rainfall intensity leading to runoff and soil moisture’s linear association. When the surface runoff exponent is 1, it implies P_max parameter, which is threshold rainfall intensity, leading to linearity between actual runoff coefficient and comparative soil moisture content. In order to measure the scale for this parameter, we can compare observed and computed peak discharges when floods are high. The calculation of evapotranspiration from soil and vegetation is done on the basis of a relationship explored by Thornthwaite and Mather (1955) [17] as defined by probable growth level, evapotranspiration, vegetation type, and soil moisture content:

 

$$\{\begin{array}{l}E=0for\theta \prec {\theta }_{\text{w}}\\ E=\left[c_v{K}_{ep}{E}_p-I\right]\left[\frac{\theta -{\theta }_w}{{\theta }_f-{\theta }_w}\right]for{\theta }_w\le \theta \le {\theta }_f\\ E=c_v{K}_{ep}{E}_P-Ifor\theta \ge {\theta }_f\end{array}$$ E4

Where c _v (-) is vegetation coefficient which varies throughout the year depending on growing stage and vegetation type, K _ep (-) is a correction factor for adjusting potential evaporation E _p (LT^-1),θ _w (L³L^-3) is moisture content at permanent wilting point, and θ _f (L³L^-3 )is moisture content at field capacity. Only increase in groundwater capillary can bring about evapotranspiration when wilting point (θ <θ w ) is higher than the water content, and groundwater storage G(m) and a scaling parameter Go(m) controls groundwater capillary rise:

$$E_g=(c{K}_{ep}{E}_p-I)G/G_o forq{q}_w$$ E5

Where E _g (LT^-1) is the evaporation from groundwater, K _ep (-) is a correction factor for adjusting potential evaporation E_p (LT^-1). Pan measurement or Pemman-Monteith or other equations utilizing accessible weather data, referring to water surface or a grass cover in large fields provide the model’s latent evaporation E_p but the actual latent evaporation may be defined by native aspects that these methods do not attend to. There is need of correction factor K _ep (-) for calculating these effects and on average, the value is somewhere around 1, and the model can measure this through a water balance simulation over the long term. Various flows are obtained using water balance equation. List includes merely interflow, percolation, groundwater flow and excess rainfall as these components were adding to stream flow as the aim was its simulation. The model evaluates overland flow and channel flow routes by a linear diffusive wave estimate of the St. Venant momentum equation where the equation models the cell’s flow process as [12, 15]:

$$\frac{\partial Q}{\partial t}+c_i\frac{\partial Q}{\partial x}-d_i\frac{{\partial }^{2}Q}{\partial x^2}=0$$ E6

Where Q (m³/s) is the flow discharge at time t (s) and location x (m), c_i is the kinematic wave celerity at cell i (m/s), d_i is the dispersion coefficient at cell i (m²/s). Manning relation can estimate the parameters c_i and d_i, which are required for defining the function of cell response, as [5]:

$$c_i=\frac{5}{3}{v}_i$$ E7

$$d_i=\frac{v_i{R}_i}{2S_i}$$ E8

Where R_i is the average hydraulic radius or average flow depth of cell i (m), S_i is the cell slope (m/m), and v_i is the flow velocity of the cell i (m/s) calculated by the manning equation. De Smedt F. et al. (2000) [3] and Liu et al. (2002, 2003) proposed an amateur passage time distribution, an estimated numerical key to the equation of diffusive wave associating the discharge when flow path concludes to the initial accessible runoff.

$$U_i(t)=\frac{1}{\sqrt{2\pi {\sigma }_i^2{t}^3/t_i^3}}\exp \left[-\frac{{(t-t_i)}^2}{2{\sigma }_i^{2}t/t_i}\right]$$ E9

Where t_i is the mean flow time from the input cell to the flow path end (s), and σ_i ² is the variation of the flow time (s²).There is spatial distribution for t_i and σ_i ², flow celerity and dispersion coefficient functions define it through convolution integral and flow paths given by topography:

$$t_i={\displaystyle \sum _{j=1}^N\left(\frac{1}{c_j}\right)}{l}_j$$ E10

$${\sigma }_i^2={\displaystyle \sum _{j=1}^N\left(\frac{2d_j}{c_j^3}\right)l_j}$$ E11

Surface water moves faster than groundwater and therefore the latter is made easy in the form of a lumped linear reservoir on sub-catchment scale, derived from GIS. Sub-catchment outlet’s groundwater flow is connected to overland flow and interflow that is routed to the primary channel from every cell, keeping in view the impact of river damping on every part of the flow. This is followed by routing the total hydrograph to the basin outlet by Eq. 6 derived channel response functions. The sum of the discharge is attained by convoluting every cell’s flow response. Simulation of every hydrological process within GIS is a benefit of this methodology.

Calibration process includes simulated discharge’s comparison with respect to observed discharge after the model has been executed. Model calibration is done after adjustments of input parameters and evaluation of the output. The input parameters being used were: groundwater recession coefficient (Kg), scaling factor for interflow computation (Ki), temperature degree-day coefficient (K_snow), initial groundwater storage (G_0), rainfall degree-day coefficient (K_rain), rainfall intensity corresponding to a surface runoff exponent of 1 (P_max), surface runoff exponent for a near zero rainfall intensity (K_run), and correction factor for potential evapotransipiration (K_ep). For minimization of the differences between simulated and observed discharge, graphical technique was used. Certain parameters were highly influential as far as simulated flow is concerned. These included groundwater recession coefficient, scaling factor for interflow computation, maximum groundwater storage, and initial groundwater storage. When observed flow was in excess of the simulated flow, there was an increase in the initial groundwater storage. When simulated discharge was in excess of the observed discharge, there was a reduction in the interflow-scaling factor. Manual adjustments in the parameters of the model were done to achieve a good match between outlet’s observed and simulated flow.

A statistics series is used for evaluation of observed hydrograph reproduced by WetSpa. Various factors are taken into account, including model confidence, evaluation on the basis of visual comparison, model efficiency, the bias, time to the peak, and evaluation of peak flow rate. Statistical data provides quantitative estimations of goodness of fit between predicted and observed values. These also indicate the extent to which observations reflect predictions. Assessment of the model’s predictive capabilities is done on the basis of test’s outcomes. Evaluation of the goodness of fit in the time to peak or the peak discharge may be done with the help of its absolute or relative errors, and the other criteria for evaluation are explained below:

1) Model Bias

Model bias is defined to be the relative mean difference between observed stream flows and predicted stream flows for a simulation sample that is sufficiently large. It should also be reflective of the ability to reproduce water balance.

$$MB=\frac{{\displaystyle \sum _{i=1}^N\left(Q{s}_i-Q{o}_i\right)}}{{\displaystyle \sum _{i=1}^{N}Q{o}_i}}$$ E12

where CR1 is the model bias, Q_si and Q_oi are the simulated and observed stream flows at time step i (m³/s), and N is the number of time steps over the simulation period.The equation 2.12 shows the criterion. It basically gives systematic over-prediction or under-prediction for prediction sets. The fit is considered better when the MB value is low. Observed flow volume’s perfect simulation is represented with the value 0.0.

2) Modified Version of Nash-Sutcliffe Efficiency for High Flow Assessment

Equation 13 describes a slightly different version of the Nash-Sutcliffe criterion. Actually, it is a fusion of Hoffmann, El Idrissi et al et al calibration condition (2004) [7] :

$$NSH=1-\frac{{\displaystyle \sum _{i=1}^N\left(Q{o}_i+\overline{Qo}\right){\left(Q{s}_i-Q{o}_i\right)}^2}}{{\displaystyle \sum _{i=1}^N\left(Q{o}_i+\overline{Qo}\right){\left(Q{o}_i-\overline{Qo}\right)}^2}}$$ E13

Where NSH is a modified version of the Nash-Sutcliffe criterion for gauging the aptitude with which the time evolution of high flows has been simulated. The formula shows how high discharges bear greater load as compared to low discharges. The value of 1 is the optimal value in NSH.

3) The Modified Correlation Coefficient (r _mod)

$$r_{\mathrm{mod}}=\left[\frac{\min \left\{{\sigma }_o,{\sigma }_s\right\}}{\max \left\{{\sigma }_o,{\sigma }_s\right\}}*r\right]$$ E14

Where $${\sigma }_o$$ and $${\sigma }_s$$ are the standard deviations of observed and simulated discharges respectively, r is the correlation coefficient between observed and simulated hydrographs.The Modified Correlation Coefficient (r _mod) is for the purpose of judging the accuracy of the replication of time progression of high flows. 1 is the perfect value for r. Model calibration allows one to discover the most suitable values for global modern parameters. The numerical information incorporated in the study is derived from Model Bias (MB), indicating the accuracy of water balance simulation, the Modified Correlation Coefficient (r_mod), which shows divergences in hydrograph shape and size [11], and the modified Nash-Sutcliffe efficiency for high flow (NSH), which weighs up the simulation of the stream flow hydrograph (13) provided by the functions below. The Aggregated Measure (AM) below is brought in to assess the efficacy of model workings during the calibration and confirmation phases (Eq 15). It will compute various features of the reproduced hydrograph such as shape, size and volume:

$$AM=\frac{r_{\mathrm{mod}}+Ns+(1-MB)}{3}$$ E15

AM requires a value of 1 in order to produce complete correlation. The time periods listed in Table 2 have been used in order to better sort out the proficiency of the model’s workings. (1; 6).

This includes the trial-and-error approach, the computerized method of numerical parameter optimization, or a mixture of the two. The PEST program was implemented for the process of design auto-calibration (4). There is a possibility of happening upon a local optimum instead of its global counterpart because the optimization subroutine is a local search technique. Clearly, in order to get suitable initial parameter values, a primary manual calibration is needed. The PEST program is used to regulate the WetSpa model with the obtained initial parameters values. One of the primary functions of the program is the forecasting of floods at the basin opening, and so we assume that the high flow reports are more significant than the low flow numbers in the model calibration routine.

Category	Aggregated measure(AM)
Excellent	>0.85
very good	0.70-0.85
good	0.55--0.70
poor	0.4-0.55
very poor	<0.40

Table 1.

Model performance categories to indicate the goodness fit.

2.2. Study area

The Neka basin is situated in Northern Iran and river is the central branch of Neka River. The Ablu station is situated next to the convergence of the Neka River and watershed covers the area of 1864.74km² up to the Ablu station. The Neka Basin is a large catchment with heights ranging from 36m to 3814m. Mean elevation is 1531m, mean slope is 24.81%. The Max flow length is 163km and stream order outlet is 133. Other characteristics are given in Table 1 and Figure 1.

Figure 1.

Location of the neka watershed in Iran.

A digital elevation model (DEM) was extracted from the topography map supplied by Iranian National Geographical Organization (IRNGO). The map was transformed into a 50m grid DEM.

Figure 2.

Elevation map of Neka basin.

Iranian Soil Conservation and Watershed Management Research Institute of Iran (SCWMRI) provided the data for ground cover (Figure 3) depicts the land use map for the study, presenting land cover for year 2000. 5 types of land cover can be observed, 45% of the basin is covered by forest land, 40% by rangeland, agriculture and villages account for 3% and the last 12% by short grass (figure3).The soil data was not available for the study, thus, soil map was obtained using a resource evaluation and land capability map (1973) from the Iranian Agriculture and Natural Resources Ministry, Soil sciences and Fertility Institute. 3 different types of soil are present in the catchment. Clay loam accounts for 89.4% of the basin, 5.6% is silt loam and the remaining 4.8% is sandy loam (Figure 3).

Figure 3.

Land use types in Neka basin.

Figure 4.

Soil types in Neka basin.

2.3. Inputs

For this study, precipitation, the potential evapotranspiration (PET) and discharge data were obtained from Water resources management of Iran(WRM.IR). The data included daily discharge data at 1 gauging station, precipitation for 7 stations and PET for 4 stations. All the data sets covered a period of 14 years i.e.1986-1990 and 2002-2004 Though daily discharge data is available for some locations in the catchment, only the Ablu station and the Neka station are being used for the study and model calibration.

2.4. Stream flow simulation by WetSpa model in Neka basin

Identification of spatial model parameters starts once the data is collected and processed for the use in the WetSpa model. Territorial features are taken out from the DEM including flow direction and accumulation, stream network, link and order, slope, elevation and hydraulic radius. To define the stream network the threshold is set to 10 i.e. when the upstream drained area exceeds 0.1km² then the cell is considered drained. To establish sub-catchments the threshold value is set at 1000, through which 265 sub-catchments were found with an average area of 6.971km².DEM extracted the slope map. The two are comparable with level slope (0.005 to 17.603%) in the middle area where as steeper slopes (17.603 to 158.389 %) were observed along the borders. The small area along the border has a very high gradient showing a ridge. Flow within the catchment is defined by the slope. In the above case the flow will be from the borders to the middle area of the catchment. The threshold value of 0.005% is considered for minimum slope while creating the grid of the surface slope. Should the calculated value be lower than the threshold value then to evade the inert water and low speeds, the value is taken as 0.005 %.If the two similar streams order come together then the stream order rises by one. It keeps on increasing as the river runs its course since more and more streams join it. Stream order ranges from 1 to 133 in this catchment.100 cells were chosen to classify first order stream. Meaning the first order stream comes into existence if the runoff from 100 cells provide to one cell. If one first order stream meets another one then they form a second stream order.Hydraulic radius was extracted from DEM (flow accumulation). Power law relationship with greater probability was employed which explains and uses the relationship between controlling area and hydraulic radius. It ranges between 0.009m for upland to the outlet of main channel at 6.098m in this catchment. With the progressing of river at the downstream the hydraulic radius is escalating, which is very obvious. Hydraulic radius is minute as the flow of water in the upstream area of the catchment is not fast. Beyond the rate of 0.5 (2-year return period). The hydraulic radius framework is produced. Subsequently, porosity, residual moisture; plant wilting point, field capacity, the grids of soil hydraulic conductivity and pore size distribution index are categorized on the basis of the soil texture grid with the help of attribute lookup Table 2.4.The first moisture map was developed with the help of soil map. Largely, in the region of catchments the moisture map differs as of 0.1-0.6 but, it usually has the high value of about 1 in the flow of the river. To measure the overspill from catchments the moisture map is used. The areas where the humidity or moisture is more, the overspill water is high as well. With the help of soil map the soil hydraulic conductivity was developed. In the flow of river, the difference of hydraulic conductivity is 0.6-1.51, in the majority regions of catchments. The map is utilized to find out the overspill from catchments.Another map developed with the help of soil map is the soil porosity. The difference of hydraulic porosity is 0.43-4.75 in the course of stream.Other things developed or derived from soil map include pore size distribution index, field capacity, Neka catchment, residual moisture and plant wilting point. The runoff coefficient is developed from the slope whereas land and soil type utilize maps. The runoff coefficient varies from 0.071 to 1 in the above catchment; conversely the runoff coefficient in the majority regions of the catchment is from 0.071 and 0.587. In the middle region of the catchment, the runoff coefficient is optimal because the land is used for cultivation and agriculture and due to the presence of clay soil. Because of the occurrence of steepness in the lofty regions, the overspill/runoff coefficient is highest. In the same way Manning’s roughness n coefficient, the grids of root depth and interception storage capacity are categorized again from land, making use of the grid. Manning’s n for channels is interposed on the basis of stream order grid/framework in which 0.050 m^-1/3s is set for lowest order, whereas 0.030 m^-1/3s is set for highest order (Figure 5).

Figure 5.

Manning coefficients map of Neka basin.

On the basis of hydraulic radius and Manning’s coefficient, the velocity map was formed. The Manning’s coefficient was acquired with the help of land use while hydraulic radius and the slope were developed from DEM. Manning’s equation helps in calculating velocity (Equation 16).

$$v_i=\frac{1}{n_i}{R}_i{}^{2/3}{S}_i{}^{1/2}$$ E16

Where, n_i is the Manning’s roughness coefficient (m³/s), R_i is the average hydraulic radius of cell i (m), S_i is the cell slope (m/m), and v_i is the flow velocity of the cell i (m/s).Velocity map shows the high velocity/speed regions (around m³/s) from the river/stream course. Besides the vicinity of high stream order, the low speed/velocity is found in catchments.The travel time map notifies the flow time of water in the catchments that differs from 0-54.91 hour. In the catchment inlet the travel time is optimal while in catchment outlet its zero. With the help of time deviation factor, the velocity of water is measured. The standard deviations are produced as a result of the movement of the flow velocity from time to a basin exit/outlet. It allows computing the IUH from every grid cell near the basin outlet. The projected standard flow time is shown in the Figure 14 between the grid cells and the basin outlet.Grids of depression storage capacity and potential runoff coefficient are attained with the help of attribute tables, which are formed by joining the grids of soil, elevation and the land use. 3 per cent is set as a proportion of impermeable region in villages. For the whole catchment the computed standard potential runoff coefficient is of 0.85. The gridirons for temperature, PET and precipitation are formed on the basis of environmental or geographical coordinates of every measuring location besides this on the catchment frontier utilizing the Thiessen polygon extension of the ArcView Spatial Analyst (Figures 6 and 7).

Figure 6.

Thiessen polygons for the precipitation stations.

Figure 7.

Thiessen polygons for the evaporation stations.

3. Results and evaluation

Forming validation and calibration, the 14 years (1986-1999) and (2002-2004) calculated PET, discharge data and the daily precipitation. From the 14 years time, the initial 12 years are selected for model calibration whereas the last two years are selected for model validation. Two models are used; the first one is the global model that is for the calibration processes while the second model being the spatial model whose factors remain unchanged. Global model factors are exclusively selected on the basis of basin traits, as talked about in the user manual of the model and the concerned documents. The imitation outcomes are evaluated to the experimental hydrograph on the Ablu station in the catchment basin statistically as well as graphically. By examining the base flow that is estranged from observed hydrograph, the first groundwater flow recession coefficient is anticipated. Amendments are mandatory of the parameters with respect to the total flow capacity and the fitting of base flow. The interflow scaling aspect is attuned for the recession and peak area of the flood hydrograph, which is receptive on behalf of high and low flows. Additionally, there are two parameters that are managing the quantity of surface runoff, which include the surface runoff exponent of about zero rainfall power and the rainfall intensity equivalent to the surface runoff exponent of about 1. These are attuned mostly for little storms as the real runoff coefficient is less because of the low rainfall intensity. With the evaluation of the water balance and hydrographs for the initial period, the active groundwater storage and soil water are attuned. The optimal active groundwater storage manages the quantity of vapour emerged from the groundwater and as a result it can be attuned by evaluating the flow amount in dry stages. As of the last two years of 14 year period the calibrated global parameters are utilized to reproduce the daily stream flow (Table 2).

No	Symbol	Parameter	Feasible range	Estimate	PEST(95%confidence interval)
1	Ki	Interflow scaling factor	0-10	1.263	1.2-1.3
2	kg	Groundwater recession coefficient (d-1)	0-0.05	0.0001	0.0079-0.0001
3	K_ss	Initial soil moisture (mm)	0-2	0.8424	0.98-1
4	K_ep	Correction factor for PET (-)	0-2	0.5373	0.53-0.54
5	G0	Initial active groundwater storage(mm)	0-500	1	1-2
6	G_max	Maximum active groundwater storage(mm)	0-2000	50	42.5-59.4
7	K_run	Moisture or surface runoff exponent (-)	0-8	7.8573	7.5-8
8	P_max	Maximum rainfall intensity (m)	0-500	133.9	130-250
9	T0	Threshold melt temperature (^oC)	-1-1	0	0-0.53
10	K_snow	Melt-rate factor (mm ^oC^-1d^-1)	-2 - 2	-2	-2-0.5
11	K_rain	Rainfall melt-rate factor (^oC^-1d^-1)	0-0.05	0	0-0.0211

Table 2.

Calibrated model global parameters.

The model is calibrated using daily stream flows observations for Neka basin. Calibration is first done manually ( trial-Error method) and then automatically using PEST to minimize the sum of square differences between observed and predicted stream flow. The validation as well as calibration periods found the performance of the model satisfactory. Tables 2 and 3 give the criteria for evaluation for the periods of calibration and validation. Hoffman et al. (2004)(7) gave four evaluation criteria, which are used here. The observations of high flow values are going to be more important in the model calibration than the low flow values because prediction of floods at the basin outlet is one of the primary purposes of the model. The evaluation in the process is carried out according to the criteria of Nash-Sutcliffe efficiency for high flow evaluation and the modified Correlation Coefficient r( mod) and Aggregated Measure(AM).The respective values were found for the validation period: the Nash-Sutcliffe efficiency was equal to 0.87, the flow volume was found to be 1.3% underestimated, whereas the modified Correlation Coefficient was 0.747 and Aggregated Measure was found to be 0.85 that performance is excellent for validation phase,Also The respective values were found for the calibration period: the Nash-Sutcliffe efficiency was equal to 0.78, the flow volume was found to be 2.2% underestimated, whereas the modified Correlation Coefficient was 0.734 and Aggregated Measure was found to be 0.83 that performance is very good for calibration phase. These results show that a lot of factors including the precipitation, antecedent moisture, and runoff generating processes were considered by the model in a manner that was spatially realistic and the basis for it was topography, land use and soil type. These factors resulted in producing highly accurate rates for the capture of both high flows as well as the general hydrological trends (Table 3).

Period	MB(CR1)	NSH(CR5)	r(mod)	AM	Performance
23/9/1986-31/12/1998; (calibration period)	0.094	0.76	0.733	0.82	very good
0.041	0.74	0.729	0.81	very good
0.118	0.75	0.725	0.78	very good
0.018	0.78	0.734	0.83	very good
22/9/2002-21/9/2004; (validation period)	0.012	0.87	0.746	0.85	Excellent

Table 3.

Model performance of WetSpa model for calibration and validation period.

Simulated flows for the whole amount of time (1986-1999) and (2002-2004). Figures 8 shows a Graphical comparison between observed and simulated river discharge in Neka basin at 2002-2004. In general, one can notice a reasonable agreement between model results and observations. Peaks in the hydrograph are rather well predicted. As well for size as for time of occurrence.

Figure 8.

Observed and simulated daily flow at Neka basin for validation phase.

Figure 9.

observed and simulated flow versu ferqunecy.

Figure 10.

scatter plot of observed versus simulated flow for validation phase.

 

4. Conclusion and recommendation

This research allowed for the model to be examined for the Neka catchment in Iran for a period of 14 years. The daily rainfall and evaporation input data was tested. A fine comparison with the calculated hydrograph was obtained. The precision of this reproduction is around 0.87, 0.974 and 0.747 according to the Nash-Sutcliffe criteria, square correlations, and modified correlation coefficient respectively. The aggregated measure (AM) was found to be 0.85 that the model performance is excellent for the validation phase. This shows that the model is suitable for using rainfall data, antecedent moisture and evapotranspiration etc and generating runoff in a spatially rational based on the topography, land use and soil type and creating a nominally higher precision simulation catchment for the high flows. The graphical comparison amongst observed and estimated daily flows for the 14 years of model simulation proves that the season flood hydrographs have been properly reproduced by the distributed hydrological WETSPA model and this is achieved via reaching conclusions from using parameters based on the topography and other features of basin such as soil and land use. The consistency of the model estimation depends on how well the model‟s construction is created along with how well it is parameterized. Model calibration is also essential to fix the model‟s workings (8). Manual and mechanical calibrations are the two kinds of parameter estimations which are used. Mechanical calibration uses a search algorithm to check the perfect parameters and allows for numerous benefits using a physical approach. WETSPA model parameterization and the spatial format of the model parameters are estimated by employing the accessible field information to define the most essential differences. This approach ensures that the model uses that data that has been represented in the catchment. In this research, the WETSPA model was first manually calibrated through a trial and error process more using more than 500 tests for finding best parameters. This approach also allows for the use of automatic calibration methods to enhance the working of the model. The research recommendations are:

• An improved investigation can also presented the impact of coordinated alterations to the weather and land use on the hydrological procedures within the region, such as impact on flooding and soil moisture distribution.

• Using practical procedures to allow for the growth which accounts for a joint effect of elevation, slope, general movement of the atmosphere etc on the spatial division of rainfall, temperature and PET. This can also lead to an increased consistency of the model inputs and lower the indecision of the model outputs. This is important when designing for a large mountainous catchment. The radar data can also be added to the WetSpa model to simulate the spatial division of rainfall at each step at a time.

• Some of the model errors are caused due to lack of accurate and efficient input data such as rainfall, temperature and evapotranspiration. Hence, for increasing the efficiency of the model, increased number of rainfall, temperature and evapotranspiration stations, as well suitable distribution of the measuring stations over the watershed is required.