Upper Bear River Analysis using TOPMODEL

Final Report CEE 6440 GIS in Water Resources

Dan Bolke

December 12, 2005

Rainfall-runoff models come in many shapes and sizes, each of them are trying to accomplish the same end result: to model the behavior of water as it meets the earth from heaven and the path it decides to go to get down the hill. Each watershed is different but if broken down into small enough pieces, all the pieces can be categorized and labeled where an understanding is achieved of what water will choose to do given certain choices.

This project covers an explanation of TOPMODEL and its use in the framework of ArcGIS in determining runoff in an upper Bear River watershed from a particular event. Goals in selection of this project are to use functions and techniques learned in the Water Resources GIS course in a practical sense and incorporating these skills into a real issue.

TOPMODEL Introduction

TOPMODEL is a set of conceptual tools used for reproducing the behavior of water in a watershed rather than a hydrological modeling package (Beven et al, 1995). Streamflow in TOPMODEL comes from two sources: Precipitation on saturated areas which directly contribute to surface runoff and subsurface downhill flow consisting of baseflow and return flow (Tarboton, 2003). The foundation of TOPMODEL is simple approximate hydrological theory which introduces a minimal number of parameters for calibration while still considering the subsurface hydrological behavior of the catchment with minimal measurements of internal state variables. Parameters are meant to be physically interpretable. The two main objectives of TOPMODEL are to model pragmatic and practical forecasting and to develop a theoretical framework or perceived hydrological processes where realism and model procedures, as well as dealing with scale issues, may be researched (Bevan et al, 1995).

Simplification of catchment properties is required for the modeling process. The hydrological processes within a catchment are distributed dynamically and heterogeneously. Soil and bedrock heterogeneity along with vegetation distribution complicate runoff modeling. Water generally flows downhill. Considering saturated zone flow where the soils are shallow relative to the hillslope scale, the water table should be nearly parallel to the surface topography due to gravitational control when the conditions are wet enough for lateral downslope saturated flow. Water will then flow from areas of steep slopes to areas of slope convergence. Streamflow is a result of different types of paths to arrive at the stream such as infiltration excess overland flow, saturation excess overland flow, return flow, or a rise in the watertable. It is impossible to take into consideration all of the perceived heterogeneity of unsaturated flows in a set of relatively simple, macroscale concepts (Beven et al, 1995).

To simplify complex nature of hydrological processes into functional mathematical structures a series of assumptions are made in TOPMODEL (Beven et al, 1995):

A1. The dynamics of the saturated zone are approximated by successive steady state representations.

A2. The effective hydraulic gradient of the saturated zone is approximated by the local topographic surface gradient, tanβ.

From these assumptions a relationship is derived between catchment storage (or storage deficit) and local levels of the water table (or storage deficit due to drainage) where the topographic index (a/tanβ) is the main component. The topographic index describes if a certain point in a catchment will reach saturated conditions. High values are a result of long slopes or upslope contour convergence and low slope angles.

A3. The distribution of downslope transmissivity with depth is an exponential function of storage deficit of depth to the water table:

(eq. 1)

where To is the latteral transmissivity when the soil just reaches saturation (m²/h), S is local saturation deficit (m) and m is the model parameter (m). This can also be stated in term so the depth to the water table as:

(eq. 2)

where z is the depth to the water table (m) and f is a scaling parameter (m^-1). The relationship between the parameters f and m is thus:

(eq. 3)

where Δθ₁ is the effective water content change per unit depth in the unsaturated zone due to rapid gravity drainage. The decay parameter m controls the effective depth in the soil profile simultaneously with the parameter T₀.

With the effective watertable gradient and saturated flow assumed to be parallel to the local surface slope tanβ as assumed in A2 then the downslope saturated subsurface flow rate qi per unit contour length (m²/h) is shown as:

(eq. 4)

and assuming quasi-steady-state flow in the soil, and

A4. The recharge rate r [m/hr] entering the water table is spatially homogeneous shown in figure 1. Flow is described as:

(eq. 5)

where a is hillslope area per unit contour length (m) that flows through point i.

Figure 1: Flow through any point i (Beven, 2000)

Combining these equations relates depth to the water table at any point to the topographic index.

(eq. 6)

The average of these depths is given by:

(eq. 7)

Assuming that r is spatially constant due to surface flow moving relatively slow because of vegetation cover, ln r can be eliminated and the average depth to the water table becomes:

(eq. 8)

where ln(a/To tanB) is the soil-topographic index, and

(eq. 9)

A separate areal mean value for transmissivity is defined as

(eq. 10)

Now local depth to the water table is finally defined as:

(eq. 11)

where z_i is the depth to the watertable at any point i and:

(eq. 12)

This is then used to predict the water table depths over the catchment based on the spatial distribution of the soil topographic index a/To tanβ, where areas of the local water table is z<0 then these areas will produce immediate runoff along with areas of a shallow depth that will become saturated (Beven, 1995).

Development and use of TOPMODEL

Originally developed to model small upland catchments in the UK, TOPMODEL has shown reasonable results with minimal calibration of parameter values. It has been successful also in simulating catchments in humid temperate regimes in the eastern US, New Zealand, and Scotland after calibration of the parameters. The parameters may be difficult to interpret physically, especially the transmissivity parameter T₀ which often yields very high values. These high values may be attributed to preferential flow pathway or zones of fractured regolith which would contribute to effective lateral downslope transmissivity values higher than expected for vertical hydraulic conductivity. Another reason is that the drainage upslope area, a, is assumed to extend to the divide. This is especially valid in drier catchments. Research has been done on the use of TOPMODEL in drier catchments. Successful simulations of runoffs from catchments at Mont-Lozere in the Cevenes, southern France and Gardon D'Anduze and Real Collobrier catchments also in southern France after calibrating the model to a small number of storms. On the contrary models of the Booro-Borotou catchmetn in the Cote d'Ivoire and cathcments in the Prades Mountains of Cataluna, Spain have resulted in satisfactory results only after the watershed has properly wetted up. In low precipitation catchments the soil may never reach a proper wetted state and have the tendency for short high intensity storms which mainly produce saturation excess overland flow which is not usually accounted for in TOPMODEL (Beven et al, 1995).

Watershed Selection

The Bear River watershed was selected because it is relatively undeveloped in the upper sections of the river. The outlet of the watershed selected is based on the location of a USGS stream gage used to evaluate the model also located in the basin are two SNOTEL sites for collecting precipitation data. Data for this basin is readily available.

Figure 2: Bear River basin with showing watershed of interest (source: E.M.R.G)

Figure 3: Head waters of Bear River showing SNOTEL sites and USGS stream gages with 12 unit HUCs

The precipitation event selected for the model occurred from August 14-18, 2003. Two days prior to the event 0.1 inches fell on the area to wet the soil after quite a long dry spell. The Lily Lake gage showed a total of 1.6 inches for the events and the Hayden Fork gage measured 1.3 inches. An average value of 1.5 inches was used for calculations covering the watershed.

Figure 4: Precipitation data from SNOTEL site at Lily Lake shown in blue box

The baseflow Qb used in calculations was 59 ft3/s or 6017 m3/hr from the figure

Figure 5: USGS stream gage showing runoff from precipitation events on august 13-18, 2003

Using TOPMODEL under the framework of ArcGIS

After the watershed has been selected and dems were obtained, the dems are clipped to an appropriate size to save time in the terrain and watershed processing. The terrain processing is done and an outlet point is added. The location of the outlet point for this watershed was selected as it is also the point of a USGS stream gage used to evaluate the model. The watershed processing is done based on this outlet point. The parameters of interest from the processing are:

demw – delineated watersheds draining through the outlet point

demsca – dinfinity contributing area

demslp – dinfinity slope

where the prefix dem indicates the input dem.

Figure 6: Watershed in upper Bear River after terrain and watershed processing from TauDEM with USGS stream gages and SNOTEL precipitation gages

A mask of the watershed is then created where each value of the raster which is in the shape of the watershed has a value of one. This is obtained by using the raster calculator in calculating:

(eq. 13)

where the expression demw>0 or where the watershed is covered by the demw will produce a true answer or a value of 1 and areas outside the demw will produce a false or 0 which makes the equation undefined, therefore producing a mask raster of the watershed where the value is one. When another raster is multiplied by the mask raster this produces a new raster which is only in the shape of the watershed. This is useful when making calculations were certain statistics are needed on the output raster such as mean or max or min values.

Figure 7: Shapefile of Individual watersheds created by TauDEM called demw with values ranging from 1 to 105

Figure 8: Showing the mask raster, each cell with value of 1 covering the watershed

This mask is then used to obtain only the areas of the outputs from the terrain processing rasters which are inside the watershed. The calculation of lambda uses the dinfinity contributing area and the dinfinity slope. The mean of the lambda calculation is the lambda value used in calculating the effective depth to the water table as explained prior.

(eq. 14)

Before the actual calculation of Z or depth to the watertable, soil data must be obtained for the effective porosity (θ_e) and the hydraulic conductivity (K). From the STATSGO website (http://www.ncgc.nrcs.usda.gov/products/datasets/statsgo/) soil data is obtained for Utah which has over 700 divisions in the state which are classified with identifications called MUID's.

Figure 9: STATSGO data showing four MUID soil classifications covering the watershed

Tables accompany these STATSGO data which contain information about the soil classification, soil composition, wildlife habitat, plant and vegetation information and hydrologic information as well as many other attributes to the soil and land use. To find the effective porosity and hydraulic conductivity (both parameters were not readily given in the STATSGO data), the classification of texture shown in table 1 was used found in the comp table which categorized the soil areas of interest into loam and sandy loam soils. Parameters such as porosity, effective porosity, wetting front suction head and hydraulic conductivity for these textures are estimated by Rawls et al. (1983). These parameters can then be added to the attribute table (table 3) where a raster can be made from the shapefile using the feature to raster tool.

Table 1: Comp table from STATSGO data with the surface textures of the four MUIDs

Table 2: Soil parameters given from soil texture, effective porosity and hydraulic conductivity are used in TOPMODEL. Numbers in parentheses are one standard deviation from given parameter value.

Table 3: Attributes of soil classifications

Soil moisture content is the water volume / total volume. The Δθ₁ used in equation 3 is the difference between the saturated moisture content and the field capacity for the soil. Figure 10 shows the estimated values of water content for three classifications of soil. Values for the loam soil were given from the figure and values for the sandy loam were averaged between the sandy soil and the loam.

Figure 10: Soil moisture content for different soils

Source: Klocke and Hergert, (1990)

The f model parameter used in calibration of the model is related to m by equation 3. m is found by plotting the inverse of Q vs time, where the slope of the recession is 1/m. This is then calibrated to the individual model from a few watersheds. Calibration of models is beyond the scope of this project but the m modeling parameter is estimated in Figure 11.

Figure 11: time vs the inverse of Q. m=1000000

Once all the parameters have been found, calculations for the average depth to the water table (Z') can be made using the raster calculator as well as individual depths (Z) using equations 8 and 11.

Areas where the Z value is less than 0 have saturated conditions and runoff will be the same as the precipitation.

Figure 12: Z values or depth to the watertable and pink values are flat which is considered saturated

Runoff occurs where the soil is saturated or will become saturated during an event. Calculations for runoff are based on the following:

Z<0 100% runoff occurs in saturated soil

Z>P/θ_e Z>P/θ_e 0% runoff occurs because not enough precipitation for saturated conditions

0<Z<P/θ_e runoff=P-Z*θ_e soil will become saturated during the event

where P is the precipitation in meters and θ_e is effective porosity.

Results and Conclusions from TOPMODEL on the Bear River watershed

After working out the model of the event in the upper Bear River Watershed some of the values seemed to be far from expected values. The Z values or depth to the watertable values ranged from -11,829,997 meters to 146,347,200 meters with an average of 10,851,589 meters. These extreme values are a result from very low values of the f modeling parameter which was estimated from the slope of the recession curve of the hydrograph, which had values on the order of 10^-7 m^-1. The lowest f parameter value observed from research papers was 0.01 m^-1 from Franchini (1995) and the high was 66.67 m^-1. Values for the baseflow Q_b were calculated to be 322,081 m³/hr which is about 3140 ft³/s. This seems over predict the baseflow by quite a bit according to figure 5 of the USGS stream gage. The value of 59 ft³/s or 6015 m³/hr which was the flow in the stream before the event occurred.

The area under the curve of figure 5 showed the total runoff from the event to be about 184,000 m³ after subtracting off the baseflow. The event runoff as calculated using TOPMODEL theories where precipitation fell on areas where Z < 0 or saturated soil or the soil would become saturated produced an event total of 621 m³. This low value raises questions on the validity of the model and the use of TOPMODEL in an arid climate such as Utah. The flat areas or where tanβ or the slope is 0 produced a runoff of 122,963 m³, producing a total of 123,584 m³.

This project was an excellent opportunity to get introduced to TOPMODEL and advancing my capabilities using ArcGIS.

References

Beven K., R. Lamb, P. Quinn, R. Romanowicz, and J. Freer, (1995), "Topmodel," Chapter 18 in Computer Models of Watershed Hydrology, Edited by V.P. Singh, Water Resources Publications, Highlands Ranch, Colorado, p.627-668.

Beven, K.J., (2000), Rainfall Runoff Modelling: The Primer, John Wiley, Chichesster

E.M.R.G. (Environmental Management Resource Group), (2005), bearriverinfo.org .

Franchini, M., J. Wendling, C. Obled, E. Todini, (1995) "Physical interpretation and sensitivity analysis of the TOPMODEL," Journal of Hydrology, 175 (1996) 293-338.

Klocke, N. L., G. W. Hergert, (1990), "How Soil Holds Water," NebGuide, G90-964. http://ianrpubs.unl.edu/fieldcrops/g964.htm

Rawls, W. J., D. L. Brakensiek and N. Miller, (1983), "Green-Ampt Infiltration Parameters from Soils Data," J Hydraul. Div. Am. Soc. Civ. Eng., 109(1): 62-70.

Tarboton, D. G., (2003), "CEE6400 Homework 8: Topmodel and Unit Hydrographs" Utah State University, http://www.engineering.usu.edu/classes/cee/6400/hw8/ .

Tarboton, D. G., (2003), "Rainfall-Runoff Processes," Utah Water Research Laboratory, Utah State University

USDA, "SNOTEL Data" (accessed Nov. 2005), http://www.wcc.nrcs.usda.gov/snotel/

USDA, "State Soil Geographic (STATSGO) Database" (accessed Dec. 2005), http://www.ncgc.nrcs.usda.gov/products/datasets/statsgo/

USGS, "Daily Streamflow for Utah" (accessed Nov. 2005), http://nwis.waterdata.usgs.gov/ut/nwis/