ABSTRACT

The spatially distributed hydrologic model WetSpa that works on daily, hourly, and minutely timescales is used to predict the flood hydrographs and spatial distribution of the hydrologic characteristics in a river basin by combining elevation, soil and land-use data within Geographical Information System. This model was applied in Ziarat river basin (95.15 km²) located in Golestan Province of Iran. Hourly hydro-meteorological data from 2008 to 2010 consist of precipitation data of two stations, temperature data of one station and evaporation data measured at one station, which were used as input data of the model. Three base maps namely DEM, land-use and soil types were produced in GIS form using 30 × 30 m cell size. Results of the simulations revealed a good agreement between calculated and measured hydrographs at the outlet of the river basin. The model predicted the hourly hydrographs with a good accuracy between 62% - 74% according to the Nash-Sutcliff criteria. To evaluate the model performance during the calibration and validation periods an Aggregated Measure (AM) was introduced that measures different aspects of the simulated hydrograph such as shape, size, and volume. The statistics of Ziarat river basin showed that the results produced by the model were very good in the calibration and validation periods.

Keywords:

Flood Prediction, Hydrological Modeling, WetSpa Model, Ziarat River Basin

1. Introduction

Rainfall-runoff models are used to perform various water resource assessments based on current and future changes in the watersheds and also in flood forecasting. Particularly, the Distributed Model Intercomparison Project (DMIP) provides a forum to examine the suitability of distributed rainfall-runoff models in flood forecasting using operational quality data in order to improve the flow modeling and prediction along the whole of the river system [1] . A distributed rainfall-runoff model which has the capability of simulating the heterogeneity of the spatial distribution of rainfall and watershed characteristics may be a better approach for flood hydrograph simulation [2] . In flood prediction and rainfall-runoff computation, distributed hydrological models, fitted to a type of data which can receive benefits of geographical information system techniques, have currently became a more feasible approach [3] . Digital elevation models, digital data of soil type and land use and also powerful GIS tools are new possibilities in hydrological researches that lead to a more data driven of modeling and understanding of the fundamental physical processes underlying the hydrological cycle [4] . In recent years, a number of methods have been developed to estimate the parameters of hydrologic model. One of the frequently used and relatively simple algorithms is the Parameter ESTimation, PEST method [5] . There are examples of the application of the PEST algorithm for the calibration of hydrologic models in the literature [6] - [14] . An application of the WetSpa model in a small watershed located in Iran on hourly time scale is presented in this study. The paper is organized as follow: Following this Section, ‘‘WetSpa model” briefly describes the WetSpa model. Section ‘‘Application” focuses on the site, data used in the study area and the model simulation. Details of the model calibrations and validation and comparisons of the results are given in Section ‘‘Results and Discussion”. Finally, conclusions of this work are presented in Section ‘‘Conclusions”.

WetSpa Model

The WetSpa model simulates the runoff and river flow in a watershed basin on an hourly time step [15] - [22] . WetSpa is able to perform spatially distributed calculations through the availability of spatially distributed data sets (digital elevation model, land use, soil and radar-based precipitation data) and GIS technology. Precipitation, interception, depression storage, surface runoff, infiltration, evapotranspiration, percolation, interflow, ground water flow, and water balance in each layer are the hydrological processes that are considered in the model. The total water balance which is considered for each raster cell is composed of a separate water balance for the vegetated soil, bare-soil, open water, and impervious part of each cell. This allows to consider the non-uniformity of the land use in per cell which depends on the resolution of the grid. A mixture of physical and empirical relationships is used to depict the hydrological processes in the model. The model can predict the peak discharges and hydrographs in each place of the channel network and the spatial distribution of the hydrological characteristics of each cell. Hydrological processes are represented in a cascading way. After the precipitation, incident rainfall first encounters the plant canopy which intercepts all or part of the rainfall until reaching the interception storage capacity and then excess water reaches the soil surface and may infiltrate in the soil zone, enter depression storage, or may be diverted as the surface runoff. Some of the infiltrated water percolates to the groundwater storage and the remained is diverted as interflow. Total runoff from a grid cell is computed as the summation of surface runoff, interflow and groundwater discharge. The root zone water balance for each grid cell is modeled continuously through equating inputs and outputs as follow:

$\frac{\text{d}\theta }{\text{d}t}=P-I-S-E-R-F$ (1)

where D [L] implicate the root depth, h [L³L⁻³] show the soil moisture, I, [LT⁻¹] is the initial loss consist of interception and depression storage, S [LT⁻¹] is the surface runoff, E [LT⁻¹] is the evapotranspiration from the soil, R [LT⁻¹] is the percolation out of the root zone, F [LT⁻¹] implicate the interflow, and t is the time [T]. The surface runoff is computed by a moisture-related modified rational method with a runoff coefficient dependent on the land cover, soil type and slope:

$S=C\left(P-I\right)\left(\frac{\theta }{{\theta }_s}\right)a$ (2)

where: θ_s = saturated soil moisture content [L³L⁻¹], Cr = potential runoff coefficient [−] depending on slope, land use and soil type, and α = empirical parameter [−]. Exponent α [−] in the formula is a variable reflecting the effect of rainfall intensity on runoff generation.

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

where Q [L³T⁻¹] implicates the discharge, t [T] shows the time, x [L] shows the distance along the flow direction, c [LT⁻¹] is the location dependent on the kinematic wave celerity, is interpreted as the velocity by which a disturbance travels along the flow path, and d [L²T⁻¹] is the location dependent on the dispersion coefficient, which measures the tendency of the disturbance to disperse longitudinally as it travels to the downstream. Assuming that the water level gradient equals the bottom slope and the hydraulic radius approaches the average flow depth for overland flow, c and d can be approximated by c = (5/3)v, and = (vH)/(2S0) [23] , where v [LT⁻¹] is the flow velocity computed by the Manning equation, and H [L] shows the hydraulic radius or the average flow depth. An approximate solution to the diffusive wave equation in the form of a first passage time distribution is applied [24] . That relates the discharge at the end of a flow path to the available runoff:

$u_i\left(t\right)=\frac{l_i}{2\sqrt{\text{π}{d}_i{t}^3}}\exp \left[-\frac{{\left(c_{i}t-l_i\right)}^2}{4d_{i}t}\right]$ (4)

where U(t) [T⁻¹] implicates the flow path unit response function, serving as an instantaneous unit hydrograph (IUH) of the flow path, which makes it possible to direct the excess water from any grid cell to the outlet of the basin or to the any downstream convergent point, t₀ [T] shows the flow time, and σ [T] is the standard deviation of the average flowtime. Two parameters t₀ and σ are spatially distributed and can be obtained through integration along the topographic determined flow paths as a function of flow celerity and dispersion coefficient.

$t_0={\displaystyle \int c^{-1}\text{d}x}$ (5)

${{\displaystyle \sigma ={\displaystyle \int \left(2d/c^3\right)\text{d}x}}}^{\text{}}$ (6)

As the groundwater movement is much slower than the surface water and near surface water system movements and the understandings about the bedrock is little, groundwater flow is simplified as a lumped linear reservoir in small GIS derived subwatershed scale. With considering the river damping effect for all flow components, overland flow and interflow are directed firstly from each grid cell to the main channel, and are joined with groundwater flow at the outlet of the subwatershed. Then the total hydrograph is routed to the outlet of the basin by the channel response function derived from Equation (4). The amount of total discharge is sum of the overland flow, interflow, and groundwater flow, and is obtained by convolution of the flow responses of all grid cells. One advantage of this approach is allowing to the spatially distributed runoff and hydrological parameters of the basin for using as inputs for the model. Inputs of the model consist of digital elevation data, soil type, land use data, and measured climatological data. Stream discharge data are optional for model calibration. All hydrological processes are simulated within a GIS framework. Because a large part of the annual precipitation is in the form of snow, snow melt simulating is done by a model based on hourly temperature data. The conceptual temperature index or degree-day method is used in this study because of its simplicity but it has not a strong physical foundation. The method replaces the full energy balance with a term linked to air temperature. It is physically sound in the absence of shortwave radiation when much of the energy supplied to the snowpack is atmospheric long wave radiation [25] [26] . The equation is as follow:

$M=\max \left[0,\left(Ksnow+KrainP\right)\left(Ta-To\right)\right]$ (7)

where M implicates the daily snowmelt [mm], Ta [˚C] shows the mean air temperature, To [˚C] shows a threshold melt temperature, Ksnow is a melt-rate factor [md⁻¹˚・C⁻¹], and Krain is a degree-day coefficient that shows the heat contribution from rainfall [d⁻¹˚・C⁻¹]. The critical melt temperature To is often intuitively set to 0˚C. The melt-rate factor Ks now is an effective parameter and may vary with location and characteristics of the snow. However, Ks now, To and K rain can be calibrated.

2. Application

2.1. Study Area

The Ziarat watershed places in the north of Iran. The watershed has an area of 95.15 km² which is a small watershed, with elevation ranging from 488 to 3027 m. The mean elevation and slope of the watershed are 1714 m and about 41.4%, respectively, from which the drainage system and area were determined as shown in Figure 1. Part of the physiographical and hydrological characteristics of the Ziarat watershed is presented in Table 1. The digital maps of topography, land use and soil type are 3 base maps used in the model.

The grid size of DEM of the river basin was 30 m from which the drainage system and area were determined. The topographic data were provided from the numerical elevation data sets of National cartographic center. All GIS data are raster based with a 30 m grid size. Figure 2 shows the Ziarat watershed; topography, flow stations, the spatial distribution of the different land uses, and the soil texture. The land-use of the area, as shown in Figure 2, consist of crop (3.6%), grassland (10.7%), forest (77.6%), agriculture (6%), Urban (0.53%). clay loam (42.9%), sandy loam (26%) and clay (31.1%) are main soil types. Land use map was prepared using Landsat 7 ETM data of the study area in Golestan province.

For this study, precipitation, temperature, discharge and potential evapotranspiration (PET) data were provided from Water Research Institute of Golestan Province, The sets include hourly precipitation for two stations, temperature for one station, PET at one station, and hourly discharge data at one gauging station. All of the data are available for a 3-year period from 2008 to 2010.

2.2. Model Simulation

Identification of the spatial model parameters is undertaken after collecting and processing of required data for use in the WetSpa model. DEM is first used to extract the terrain features at each grid cell including elevation, flow direction, flow accumulation, stream network, stream link, stream order, slope, and hydraulic radius. The threshold of stream network delineating is set to 10 namely when the upstream drained area is greater than 0.1 km² a cell is considered be drained by a stream. The value of subwatersheds determining threshold is set to 3000, by which 17 subwatersheds are identified with an average subwatershed area of 5.5 km². A threshold value of minimum slope of 0.01% is considered in creating the grid of surface slope. If the calculated slope is less than the threshold value, the slope will set to 0.01% in order to avoid stagnant water or extreme low velocities.

The grid of hydraulic radius is generated with an exceeding frequency of 0.5 (2-years return period) which is resulted in an average hydraulic radius of 0.007 m for the upland cells and 0.5 m at the outlet of the main river channel (The related formula and calculation method are discussed in the documentation and user manual of the model by [27] ). In the next stage, the grids of soil hydraulic conductivity, porosity, field capacity, residual moisture, pore size distribution index, and plant wilting point are reclassified based on the soil texture grid through an attribute lookup table. Similarly, the grids of root depth, interception storage capacity, and Manning’s roughness coefficient are reclassified from the land use grid, in which the Manning’s coefficient for channels. The grids of potential runoff coefficient and depression storage capacity are obtained through attribute tables by combining the grids of elevation, soil and land use, for which the percentage of impervious area within an urban cell is set to 30%. The results are shown in Figure 3. From the figure it follows that non-forested and steeper areas share a very high potential runoff coefficient, but the forested and gentle slopes generate less surface runoff. The calculated average potential runoff coefficient is 0.58 for the entire watershed. The Thiessen polygon extension of the ArcView Spatial Analyst is used to create the grids for precipitation, temperature and PET based on the geographical coordinates of each measuring station and the watershed boundary.

Finally, the flow routing parameters include flow velocity, average travel time and its standard deviation from cells to the watershed outlet and to the subwatershed outlet. Figure 4 shows the calculated mean travel time from cells to the basin outlet for Ziarat watershed. Flow time for the most remote areas is around 10 h and the mean travel time for the watershed is 3 h.

3. Results and Discussion

Calibration and Validation

Model calibration and validation is done by the 3 years (2008-2010) measured hourly precipitation, temperature, PET, and discharge data. The 3-years period is divided into a 2-years and a 1-year period that the former is used for the model calibration and the second is used for the model validation. The calibration process is only performed manually for the global model parameters, whereas the spatial model parameters are kept as they are. Initial global model parameters are specifically chosen according to the basin characteristics as discussed in the documentation and user manual of the model [24] . Then, graphically and statistically comparison is done between the simulation results and the observed hydrograph in Ziarat.

Figure 5 presents a graphical comparison between the observed and calculated hourly flow in Ziarat for the year 2008 of the calibration period. Simulation

of the hourly stream flow for the second part of the data namely from 2010 is done using the calibrated global parameters in order to the model validation.

Figure 6 presents a graphical comparison between the observed and calculated hourly flow in Ziarat for the year 2010 of the validation period. Figure 5 and Figure 6 show that the model can well reproduce both the autumn-and winter flood hydrographs. The model performance provides satisfaction in both calibration and validation periods. Evaluation criteria for the calibration and validation period are presented in Table 3. Therefore, a primary manual calibration is needed to obtain proper initial parameter values. With the resultant initial

parameter values, the PEST program is applied to calibrate the WetSpa model. The statistical criteria used in the Model analysis are ModelBias (MB), reflecting the ability of reproducing the water balance, the Modified Correlation Coefficient (r_mod), which reflects differences both in hydrograph size and in hydrograph shape [27] , and the Nash-Sutcliffe efficiency (NS), which evaluates the ability of reproducing the stream flow hydrograph [28] , given by following expressions:

$r_{\mathrm{mod}}=r\times \left[\frac{\min \left\{{\sigma }_o,{\sigma }_s\right\}}{\max \left\{{\sigma }_o,{\sigma }_s\right\}}\right]$ (8)

where ${\sigma }_o$ and ${\sigma }_s$ show the standard deviations of the observed and simulated discharges respectively, r indicates the correlation coefficient between observed and simulated hydrographs. The perfect value for this criterion is 1.

$MB=\frac{{{\displaystyle \sum }}_{i-1}^n\left(Qsi-Qoi\right)}{{{\displaystyle \sum }}_{i-1}^{n}Qoi}$ (9)

where, MB is the model bias, Qsi and Qoi are the simulated and observed streamflows at time step i (m³/s), and N is the number of time steps over the simulation period. Model bias measures the systematic under or over prediction for a set of predictions. A lower MB value indicates a better fit, and the value 0.0 represents the perfect simulation of observed flow volume.

$NS=1-\frac{{{\displaystyle \sum }}_{i-1}^N{\left(Qsi-Qoi\right)}^2}{{{\displaystyle \sum }}_{i-1}^N{\left(Qoi-Qo\right)}^2}$ (10)

where, NS is the Nash-Sutcliffe efficiency used for evaluating the ability of reproducing the time evolution of streamflows. The NS value can range from a negative value to 1, which 1 indicates a perfect fit between the simulated and observed hydrographs. Aggregated Measure (AM) is used to measure the different aspects of the simulated hydrograph such as shape, size and volume

$\text{AM}=\frac{r_{\mathrm{mod}}+NS+\left(1-\left|MB\right|\right)}{3}$ (11)

A value of 1 for AM shows a perfect fit. To categorize the goodness of the model performance, the intervals listed in Table 2 have been adopted [29] .

4. Conclusion

It has been proven that the distributed models are useful in scenario analyses because of their ability to predict the effect of spatially changing variables. To avoid the inherent complexity in estimating the surface runoff a simple but effective approach is presented wherein the whole basin is divided into grid cells, giving the possibility to simulate the hydrologic processes at reasonably small scale [30] . Coupling to GIS makes WetSpa a powerful tool to capture local

complexities of a watershed and temporal variation of river flows, especially peak discharges [3] . The mentioned model was tested in Ziarat Watershed of Golestan Province in Iran with 3 years of observed hourly rainfall, temperature and evaporation data. The results indicated a good agreement with the measured hydrograph at the outlet of the basin, also the model is able to consider the precipitation, antecedent moisture and runoff-generating processes in a spatially realistic manner based on topography, land use and soil type; the model is especially useful to analyze the effects of topography, soil type, and land use on the hydrologic behavior of a river basin resulting in a fairly high accuracy for both high and low flows, and the general hydrological trends are well captured by the model. Finding optimal values for the global model parameters can be accomplished by model calibration, e.g. manually through trial and error, automatically through numerical parameter optimization, combination of both techniques. For the model auto-calibration, the PEST program was used [5] . The objective function used by PEST is the sum of the squared differences between observations and predictions. Since the optimization algorithm is a local search method, there is a risk of locating a local optimum rather than the global optimum. The model results are very good based on the performance categories defined in Table 3. To evaluate the model performance during the calibration and validation periods, an Aggregated Measure (AM) is introduced that measures different aspects of the simulated hydrograph such as shape, size and volume. The statistics for Ziarat basin show that the model produces very good results for the calibration and validation period that it is in agreement with studies of [31] [32] [33] [34] [35] . Generally speaking, this study shows that the model can be useful for investigating the simulation and prediction of streamflows at the outlet of the basins with different rain condition and antecedent soil moisture.

Cite this paper

Azizi, M., Mohajerani, A. and Akhavan, M. (2018) Simulating and Prediction of Flow Using by WetSpa Model in Ziyarat River Basin, Iran. Open Journal of Geology, 8, 298-312. https://doi.org/10.4236/ojg.2018.83019

References