WetSpa is a spatial hydrological model used for flood analysis and watershed management. It is used for forecasting the effects of the spatial variable inputs on outputs by linking with GIS. Such an input for the case of WetSpa is the Precipitation; physical control like Topography, Soil and Land use and the output is the discharge. It is very useful for setting up watershed management plans.

This model is applied to the Margecany-Hornad river basin (1131km²). Daily hydro-meteorological data from 1991-2000 including precipitation data, temperature data and evaporation were used as input to the model together with three base maps DEM, Land use and Soil type with spatial resolution of 100 x 100 m cell size. Model uses geo-spatially reference data as input for deriving model parameters which includes most data types supported by Arcview, such as coverage, shape file, grid and ASCII file.

WetSpa model uses multiple layer to represent the water and energy balance for each grid cell taking into account the processes of precipitation, interception, snowmelt, depression, infiltration, evaporation, percolation, surface runoff, interflow and groundwater flow. The main objective of this model is to simulate the stream flow using DEM, Soil type, Land use based on hydrometrological data.

2. Methodology

2.1 Master data for WetSpa model

2.1.1 Digital elevation model (DEM)

The spatial and elevation resolution should be fine enough to capture the essential information allowing to taking care of the effect of spatial variability of the watershed characteristics on its hydrological response. In practice, the chosen resolution must allow adequate representation of the actual topography and accurate determination of watershed area and its river network, and its sub watersheds. The slope, aspect, flow direction, flow length and flow accumulation of each grid cell was derived from DEM. The DEM may contain artificial pits from which no water can flow out. These specific problems have to be solved by modifying elevations artificially to lead to flow directions as accurate as possible on any of the cells. Next, the identification of river network was performed by assuming that all cells draining more than a specified upstream area are part of that network. More or less detailed river networks can be identified depending on the selected upstream threshold area. Finally, the stream links, stream orders and the sub watersheds corresponding to these river reaches were identified.

2.1.2 Land use

Land use information is an important input to the WetSpa Extension, which is normally obtained from remotely sensed data for the same area as the DEM, and with the same grid cell size. For hydrological simulation purpose, all land use classes initially determined are grouped together into 14 WetSpa classes significantly distinguished from each other on the basis of their effect on hydrological processes, namely crop, short grass, evergreen needle leaf tree, deciduous needle leaf tree, deciduous broad leaf tree, evergreen broad leaf tree, tall grass, irrigated crop, bog marsh, evergreen shrub, deciduous shrub, bare soil, impervious area and open water. Each of these classes is characterized by quantitative attributes. For simulation purpose, the percentage of bare soil and impervious areas were estimated for each grid cell based on the high-resolution land use map.

2.1.3. Soil type

Soil types of the catchment were obtained from the soil information furnished by soil maps. The soil code system used in WetSpa Extension is based on the soil texture triangle developed by the United States Department of Agriculture (USDA), which is characterized by its percentage of clay, silt and sand, ranging from the fine textures (clay), through the intermediate textures (loam); and the coarser textures (sand). Therefore, the original soil coverage map has to be converted to a raster map with WetSpa soil codes in the phase of data preparation. The grid must be adjusted to the same grid structure as the DEM. The soil properties and hydraulic characteristics of those soil types are considered constant in the present version of the model.

2.2 Parameters deriving using base maps

The parameter derivation using GIS begins with an input of the catchment’s characteristics DEM (digital elevation model), Land use map and Soil type map where DEM, Land use and Soil type are 3 base maps used in the model. The DEM was used for deriving stream characteristics like the slope, stream order, flow direction, stream link, hydraulic radius, flow accumulation and watershed among others. From the Land use map, parameters like infiltration, root depth and Manning’s coefficient were derived. From the Soil type map the characteristics of the soils like initial moisture, porosity, residual moisture, conductivity; wilting point and field capacity were obtained. Some other parameters like runoff coefficients were obtained from the combination of slope, Land use and Soil type. The parameters obtained were then used for the WetSpa model.

2.3 Flow computation

WetSpa uses water balance equation to derive the different types of flow. Since the objective of this project is to simulate the stream flow, only those components which contribute on the stream flow were listed. They are interflow, percolation, groundwater flow and excess rainfall.

2.3.1 Determination of interflow

Interflow is defined as the water which infiltrates the soil surface and moves laterally through the upper soil layers until it enters a channel, which includes litter flow, return flow, unsaturated through flow, saturated through flow and so on, but excludes the saturated groundwater flow. Due to the delayed flow time, interflow usually contributes to the falling limb of a flood hydrograph, but it may also be a part of peak discharge at the basin outlet, particularly for the areas with steep slope and forest cover in humid or template regions. The following equation was used for computation of interflow in WetsSpa model.

Where,

RI(t)= Amount of interflow out of the cell over time ∆t.

D= Root depth

S= Cell slope

K[]= Cell hydraulic conductivity at moisture content

W= Cell width

C=Scaling factor depending upon land use.

2.3.2 Determination of percolation

Percolation or groundwater recharge refers to the natural process by which water is added from soil water zone to the saturation zone of the aquifer. Groundwater recharge is an important component in the root zone water balance, which connects the soil water and the saturated groundwater. The main influencing factors to the groundwater recharge are the hydraulic conductivity, root depth, and water content of the soil. In WetSpa, percolation out of root zone is assumed to pass directly to the groundwater reservoir, and estimated based on the Darcy’s law, being the product of hydraulic conductivity and the gradient of hydraulic potential. The hydraulic conductivity at effective saturation was used to define percolation (Brook et al. 1966)

RG(t)=

Where,

RG(t)= Percolation

K= Saturated hydraulic conductivity

(t) = Moisture content of soil at time t

= Saturated moisture content

= Residual moisture content

2.3.4 Determination of Ground Water flow

For the purpose of stream flow prediction, an estimate must be made of flow from the groundwater storage into the stream for each time step. Since little is know about the bedrock, the simple concept of a linear reservoir is used to estimate groundwater discharge on a small subcatchment scale. The groundwater outflow is added to any runoff generated to produce the total stream flow at the subcatchment outlet. The general groundwater flow equation can be expressed as.

QG(t)= Cg[SG(t)/1000]^m

Where

QG(t)= The average groundwater flow at the subcatchment outlet (m³/s)

SG(t)= The groundwater storage of the subcatchment at time t (mm),

m= Exponent; m = 1 for linear reservoir.

Cg= Groundwater recession coefficient taking the subcatchment area into account

2.3.5 Determination of Excess Runoff

Excess rainfall is that part of rainfall in a given storm, which falls at intensities exceeding the infiltration capacity of the land surface. It may stay temporarily on the soil surface as depression, or become direct runoff or surface runoff at the watershed outlet after flowing across the watershed surface under the assumption of Hortonian overland flow. Direct runoff forms the rapidly varying portions of watershed hydrographs and is a key component for estimating the watershed response. The excess rainfall is calculated by above equations.

F(t)= P(t)-I(t)- PE(t)

Where,

PE(t)= Infiltration excess on cell i over time interval

I(t)= Interception loss on i^th cell at time interval.

F(t)= Infiltration on i^th cell at time interval.

P(t)= Precipitation over i^th cell at time interval.

a= exponent related with rainfall intensity

The parameter Ci has rather stable regularity under ideal condition. Default rainfall excess coefficients for different slope, soil type and land use are taken from literature (Kirkby 1978).

2.3.6 Total runoff

The basic principle that WetSpa model uses is the equation of mass balance. The total runoff in the stream is equals to the sum of excess runoff, interflow and base flow.

2.4 Hydrological data preparation

In WetSpa model, four types of data are needed. They are precipitation, potential evaporation, temperature and measured discharge. All these data should be in same time scale and same period. These data were stored in text format in subdirectory input under model. The unit of precipitation and potential evaporation is in mm. The unit of temperature is taken in degree centigrade and unit of discharge is used in m³/sec format.

2.5 Global Model Parameter

To make the calibration process simple and user friendly, WetSpa model has summarised calibration in 12 global parameters. They are correction factor of PET, interflow scaling factor, groundwater recession coefficient, initial soil moisture, initial groundwater storage, base temperature for snowmelt, temperature degree-day coefficient, rainfall degree-day coefficient, surface runoff exponent, and the rainfall intensity corresponding to surface runoff exponent of 1. These parameters have physical interpretations and are important in controlling runoff production and hydrographs at basin outlet, but difficult to assign properly on a grid scale. Therefore, calibration of these global parameters against observed flow data is preferable.

i. Correction factor for potential evapotranspiration (K_ep):

The potential evapotranspiration data used in the model were obtained from pan measurement or calculated by Pemman-Monteith or other equations using available weather data. These reference evapotranspiration rates refer to water surface or a grass cover in large fields. Actual reference or PET rates, however, may depend on local factors that are not addressed by these methods like effect of land use, elevation, as well as the micro-meteorological conditions for the grid to be simulated may be different from those prevailing at the site of the meteorological station whose data are being used. To account for these effects, a correction factor is required in the computed PET.

ii. Scaling factor for interflow computation (Ki):

Interflow is an important factor for runoff contribution in humid temperature region with sloping landscapes and well-vegetated cover. In WetSpa model, the interflow is assumed to occur when soil moisture exceeds the field capacity and there is sufficient hydraulic gradient to move the water. Due to change in conductivity, porosity with depth even in homogeneous soil, quick flow of water due to stream through root canals, animal tunnel, and subsurface erosion may become a critical component of peak flow. To account these phenomena, scaling factor is used.

iii. Groundwater recession coefficient (Kg):

The groundwater recession coefficient reflects the storage characteristics of the sub watershed and is same for all hydrographs at a given location. Calibration of this parameter is necessary to adjust the computed hydrograph equal to observed low flow hydrograph.

iv. Relative initial moisture content (K_ss):

Initial moisture is key element in the model for controlling the hydrological process for realistic starting point but is not important for the long-term flow simulation. The initial moisture content was expressed as a ratio of field capacity of the soil. It affects the hydrological processes only in the initial part of the simulation. It may also be estimated from the topographical wetness index (TWI) which is simply the logarthemic ratio of upslope drainage area to the local slope (Moore et al. 1993).

v. Initial groundwater storage (G_0):

In WetSpa model, groundwater balance is maintained on subcatchment scale. The groundwater contributes to the surface stream flow, which eventually reappears as base flow. During calibration, it was adjusted by comparing the computed and observed low flow for the initial phase.

vii. Temperature degree-day coefficient (K_snow):

It is a measure of effect of temperature on snow melt. It is measured in mm/degree/day. This in general varies with both in time and space and land use (Singh et al. 2000).

viii. Rainfall degree-day coefficient (K_rain):

The rainfall degree-day coefficient determines the rate of snow melting caused by condensation of humid air on the snow surface and advective heat transfer to snow pack by precipitation and is used for additional snowmelt due to rainfall. If zero value is given, the effect of rainfall on snowmelt will not be considered.

ix. Surface runoff exponent for a near zero rainfall intensity (K_run):

In WetSpa Extension, this exponent is assumed to be a variable starting from a higher value for a near zero rainfall intensity, and changing linearly up to 1 along with the rainfall intensity, when the predetermined maximum rainfall intensity is reached. If an exponent value 1 is given, the actual runoff coefficient is then a linear function of the relative soil moisture content, and the effect of rainfall intensity on the runoff coefficient is not taken into account.

x. Rainfall intensity corresponding to a surface runoff exponent of 1 (P_max):

This parameter corresponds to threshold rainfall intensity over which the surface runoff exponent equals 1, and the actual runoff coefficient becomes a linear function of the relative soil moisture content. Calibration of this parameter was performed by comparison of the observed and computed surface runoff volume and the peak discharge for high floods.

2.6 Model Calibration

The calibration process involved the comparison of the simulated discharge with respect to observed discharge after running the model. The model was calibrated by adjusting the input parameters and evaluating the out put. The input parameters adjusted were scaling factor for interflow computation (Ki), groundwater recession coefficient (Kg), initial groundwater storage (G_0), temperature degree-day coefficient (K_snow), rainfall degree-day coefficient (K_rain), surface runoff exponent for a near zero rainfall intensity (K_run), rainfall intensity corresponding to a surface runoff exponent of 1 (P_max), correction factor for potential evapotranspiration (K_ep). The graphical technique was used to minimize the differences between the observed discharge and the simulated discharge. The main parameters that had a high influence on simulated flow were scaling factor for interflow computation, groundwater recession coefficient, initial groundwater storage and maximum ground water storage. For cases where the observed flow was above the simulated flow, initial ground water storage was increased. Interflow scaling factor was decreased when the simulated discharge was above the observed discharge. The K_snow factor was increased when the calibrated flow was below the observed flow but its effect was not much influential.

The model parameters were manually adjusted till a good match was obtained between the simulated flow and observed flow in the outlet.

2.7 Model evaluation

In addition to graphical comparison of result, the input parameters were also evaluated to check how close they were to the best values after running the model. The model evaluation coefficient gives the statistical assessment of the model performance. The results of this statistical assessment were shown by a set of output parameters. These are C1, C2, C3, C4 and C5.

C1 is the model bias and it is defined as relative mean difference between the predicted and observed flow for sufficiently large simulation. Lower value of C1 means better result. The ideal value is 0.

C1=

Where,

= simulated discharge

= observed discharge

n= number of time steps

C2 is model confidence which was calculated as the portion of the sum of square of the deviations of the simulated and observed discharge from the average observed discharge (Nash, et al. 1970). It varies in between 0 to 1 and its best value is 1.

C2=

Where,

= mean observed discharge.

C3 is the Nash-Sutcliffe efficiency. It shows how well the flow is simulated by the model. It involves the standardization of the residual variance. Its value can range from negative to 1 with 1 indicating a perfect fit between the simulated and observed hydrograph.

C3=1-

C4 is the logarithmic version of the Nash-Sutcliffe efficiency for low flow. It gives emphasize for evaluating the quality of low flow simulation and its best value is 1.

C4=1-

Where,

= arbitrarily chosen sufficiently low value to avoid the problems with nil observation or simulated discharge.

C5 is the adapted version of the Nash-Sutcliffe efficiency for high flow evaluation. It is the criteria for evaluating the ability of reproducing the time evolution of high flow as it can be seen that more weight is given to high discharge than low flow. Its perfect value is 1. After running the model the values of the parameters were obtained. In the beginning, these Ci values were very far from its ideal values. The global parameters were then adjusted and the model was run till a good match was obtained.

C5=1-

3.0 Results

3.1 Description of data and result of simulation

The DEM, Soil type and Land use maps are shown in Figure 1, Figure 2 and Figure 3 respectively. The basin has a hilly topography around the catchment boundary with elevation ranging from 332 to 1556 m. It is flat at the central and eastern part of the catchment while it has tall mountains in catchments boundary. The soil was classified as texture categories of loam, sandy loam, silt loam, silt clay loam, clay loam, sandy clay loam and loamy sand, sandy clay and clay loam. The soil mostly contains some fraction of loam and can support vegetation and that explains why it has vegetation in almost the entire catchments. The Land use map shows different types of land use on the catchments. The land use types in the catchments were more or less equally distributed. It consist of mix farming & grass, evergreen needle, deciduous broad tree, evergreen tree, tall grass, Bog marsh, Evergreen and deciduous shrub, Bare soil an stream & impervious areas. The evergreen needles were concentrated close to the catchments boundary and the mix farming and grass were concentrated around the flat central part of the catchment.

Some preprocessing of three master maps were done to get the map of slope, stream order, hydraulic radius, initial moisture, runoff coefficient, velocity distribution and travel time to outlet. They were shown in Figure 4 to Figure 10 respectively. Similarly, thiessen polygons were created for the precipitation stations (Figure 11), evapotranspiration (Figure 12) and temperature (Figure 13). The result of comparison of simulated discharge verses observed discharge is shown in Figure 14.

Figure1: Elevation map. In the DEM, it can be see that flatter area is on the central part of the catchment while the elevated areas are close the boundary.

Figure 2: Soil map. In soil map, one can see the distribution of soil type in the catchment. The loam soil has dominancy than other type of soil. Sandy clay and clay loam has little fraction on the central part of the catchment.

Figure 3: Land use map. There is grass and mix farming on the central part of the catchment while deciduous tree and evergreen tree on the boundary of the catchment. Also, it can be seen the impervious area and bare soil near to the central part but its area is limited.

Figure 4: Slope map. It is derived from DEM. The slope map looks similar to DEM with the steep slopes (11.749 to 93.919 %) in the catchment boundaries and flat slopes (0.01 to 11.749%) in the central area of the catchment. The very high slope is limited to small area near to the boundary showing the ridge. The slope defines the flow within the catchment. In the given catchment above, the flow occurs from the boundary zone to the central part of the catchment.

Figure 5: Stream order map. Stream order increases by one level if two stream of same order meets each other. It goes on increasing as the river moves downstream cause of addition of more and more stream on it. In this catchment, the stream order is ranging from 1 to 295. In this catchment, 100 cells to define first order stream was selected. It means if the runoff from 100cell contributes to one cell, it starts forming first order stream. If first order stream meets another first order stream, it forms 2^nd order stream, and goes on like this.

Figure 6: Hydraulic radius map.

It was derived from the DEM (flow accumulation) using power law relationship with an exceeding probability which relates hydraulic radius to controlling area. In this catchment it is ranging from 0.005m to 1.461m. It can be clearly seen that hydraulic radius is increasing as the river moves toward downstream. In upstream part of the catchment area, the flow is small hence hydraulic radius is also small and vice versa.

Figure 7: Initial moisture map: The initial moisture map was derived from the soil map. In most of the area of the catchments, the initial moisture varies from 0 to 0.444 but has higher value (close to 1) in the stream path. This map was used for determining the runoff from the catchments. Runoff will be high in the zone having high initial moisture content.

Figure 8: Runoff coefficient map. This map was derived from the slope, Soil type and Land use maps. In catchment above, the runoff coefficient is ranging from 0.071 to 1, however most of the catchments area consist of the runoff coefficient between 0.071 and 0.587. The runoff coefficient is the highest in the central part of the catchment possibly due to the land use of crop and mixing farming in the area and existence of the clay soil. The runoff coefficient is also high in the high-elevated area around boundary because of steepness.

Figure 9: Velocity map.

The velocity map was based on the slope, hydraulic radius and Manning’s coefficient. The slope and hydraulic radius were derived from the DEM and the Manning’s coefficient was obtained from the Land use.

The velocity was computed using Manning’s equation.

Where, nis the Manning’s roughness coefficient (m-1/3s), R is the average hydraulic radius of cell i (m), S is the cell slope (m/m), and v is the flow velocity of the cell i (m/s).

From the velocity map, high velocity areas (around 3m/s) were found in the stream path. The low velocities were found in the catchment other than the location of stream of high order.

Figure 10: Travel time map to catchments outlet. In this map, it can be seen that the travel time of flow varies from 0 to 79.106 hr. The travel time is highest for the catchment inlet and zero at the catchment outlet. The travel time map was used in the flow routing process of the model. This ensures that the flow velocity was estimated as a time variant variable.

Figure 11: Thiessen polygon for precipitation stations. There were 44 precipitations stations including both inside and outside of the catchment. Thiessen polygons were used to represent the area corresponding to particular rain gauge stations.

Figure 12: Thiessen polygon for Evapotranspiration. There were four-evapotranspiration stations in total but only two of them have been used for the water balance analysis of the catchmenst.

Figure: 13 Thiessen polygons for temperature . Five temperature stations have been used in this catchment.

After large number of trials for adjusting the input parameters and running the model, the best result of output parameters obtained are shown below. It can be seen a good match between the simulated stream flow and observed stream flow.

Figure14: Comparison of the simulated and the measured discharge. In this graph the simulated discharge was found to be close to the observed discharge.

The corresponding values of model output evaluation coefficients were:

Model bias (C1) =0.025

Model confidence (C2) =0.949

Nash Sutcliffe coefficient (C3) =0.738

Nash Sutcliffe coefficient for low flow (C4) =0.724

Nash Sutcliffe coefficient for high flow (C5) =0.762

And corresponding global input parameters used were:

Interflow scaling factor (Ki) = 4

Groundwater recession coefficient (Kg) =0.005

Initial relative moisture content (K_ss) =1.0

Correction factor for potential evapotranspiration (K_ep) =0.803

Initial ground water storage (G_0) =5.15

Maximum groundwater storage (G_max) =35.33

Temperature degree day coefficient (K_snow) =0.8

Rainfall degree day coefficient (K_rain) = 0.15

Surface runoff exponent for a near zero rainfall intensity (K_run) =6

Rainfall intensity corresponding to surface runoff exponent 1 (P_max) =0.762

4.0 Discussion

The spatial hydrological model Wetspa was used to simulate the stream flow in the Margecany-Hornad river basin. Three base maps DEM, Soil type and Land use map were used together with hydro-metrological data mainly precipitation, evapotranspiration and temperature. The base maps were used to derive other maps like runoff coefficient, slope, stream order etc. With all data prepared, the model was run. The 12 global parameters were adjusted till a good match was obtained between the simulated flow and observed flow. The result of the simulated result obtained from the model was evaluated by the statistical assessment. The results of this statistical assessment were shown by a set of output parameters C1, C2, C3, C4 and C5. The result showed that a good match was obtained between the simulated and observed flow.

5.0 Conclusion

Wetspa model was used to simulate the stream flow in the Margecany-Hornad catchment. This model used three base maps together with some hydro-metrological from 1991 to 2000 as input data. They were precipitation, evapotranspiration, temperature and measured discharge. The calibration process was simplified by use of 12 global input parameters. The model was run and its performance was evaluated in terms of C1, C2, C3, C4 and C5 coefficient whose values were 0.025, 0.949, 0.738, 0.724 and 0.762 respectively. After running the model several times with adjusted global parameters, a good match was obtained between the simulated and observed flow. Thus this model can be used to predict the flow in the stream thus it is very useful for the flood prediction and watershed manangement.

6.0 References

Brooks, R.H., and Corey, A.T., 1966. Properties of porous media affecting fluid flow.

Kirkby, M.J. (ed.), 1978. Hill Slope Hydrology

Maidment, D. R. 2002. Arc Hydro, Gis for Water Resources.

Moore, I.D., Gessler, P., Nielsen, G. and Peterson, G., 1993. Soil attribute prediction using terrain analysis.

Nash, J.E. & Sutcliffe, J.V., 1970. River flow forecasting through conceptual model

Singh, P., Kumar, N., and Arora, M., 2000. Degree-day factors for snow and ice for the Dokriani Glacier, Garhwal Himalayas