r.hydro.CASC2D - GRASS raster command to execute fully
integrated distributed hydrologic modeling.
(GRASS Raster Program)
Note: This is beta software: version 1.01.
r.hydro.CASC2D, r.hydro.CASC2D help, r.hydro.CASC2D [-smtoepibduq] elevation=mapname time_step=value tot_time=value discharge=mapname outlet_eastNnorthNslope=east,north,bedslope rain_duration=value [watershed_mask=mapname] [initial_depth=mapname] [storage_capacity=mapname] [interception_coefficient=mapname] [roughness_map=mapname] [Manning_n=value] [conductivity=mapname] [capillary=mapname] [porosity=mapname] [moisture=mapname] [pore_index=mapname] [residual_sat=mapname] [lake_map=mapname] [lake_elev=mapname] [radar_intensity_map=mapname] [links_map=mapname] [nodes_map=mapname] [channel_input=mapname] [table_input=mapname] [dis_profile=mapname] [wat_surf_profile=mapname] [hyd_location=mapname] [r_gage_file=mapname] [unif_rain_int=value] [num_of_raingages=value] [gage_time_step=value] [radar_time_step=value] [write_time_step=value] [unit_el_conv=value] [unit_lake=value] [unit_space=value] [d_thresh=value] [dis_hyd_location=mapname] [depth_map=mapname] [inf_depth_map=mapname] [surf_moist_map=mapname] [rate_of_infil_map=mapname] [dis_rain_map=mapname]
In above, the command line arguments have been rearranged so that the required parameters are placed first. When doing "r.hydro.CASC2D help", the user will see a different sequence of parameters.
Users with no background in or understanding of distributed hydrology are strongly advised against using this code in any mode, particularly in operational mode. Besides knowledge of basic hydrology, experience with typical numerical techniques used in physically-based hydrodynamic models is recommended as it will help the user grasp capabilities and limitations of this model. This manual is significantly condensed for electronic distribution and is in no way comprehensive. Users are encouraged to experiment with the model and venture in hydrology textbooks and journal papers to learn more about the topics touched upon in this manual. The r.hydro.CASC2D code is continuously being improved. Changes in the source code of r.hydro.CASC2D may be made at any time, without notification. No claims are made regarding the suitability of r.hydro.CASC2D for any particular purpose. The model r.hydro.CASC2D was written for research and educational purposes.
CASC2D originally began with the two-dimensional overland flow routing algorithm developed in APL by Prof. P.Y. Julien at Colorado State University. The overland flow routing module was converted from APL to FORTRAN by Bahram Saghafian, then at Colorado State University, with the addition of Green & Ampt infiltration and explicit channel routing (Julien and Saghafian, 1991, Saghafian, 1992, and Julien et al., 1995). The Fortran version was reformulated, significantly enhanced, and re-written in the C programming language by Bahram Saghafian at the U.S. Army Construction Engineering Research Laboratories. Implicit channel routing code was developed and added to r.hydro.CASC2D by Fred L. Ogden (Ogden, 1994), formerly at Colorado State University, now Asst. Professor, Department of Civil and Environmental Engineering, University of Connecticut, Storrs, Connecticut. This version became known as r.hydro.CASC2D, part of the GRASS GIS for hydrologic simulations (Saghafian, 1993).
r.hydro.CASC2D is a physically-based, distributed, raster hydrologic model which simulates the hydrologic response of a watershed subject to a given rainfall field. Input rainfall is allowed to vary in space and time. Major components of the model include interception, infiltration, and surface runoff routing. Interception is a process whereby rainfall is retained by vegetation. Interception is estimated using an empirical three parameter model. Infiltration is the process whereby rainfall or surface water is pulled into the soil by capillary and gravity forces. The Green and Ampt equation with four parameters is applied to model the event-based infiltration. For continuous soil moisture accounting, redistribution of soil moisture can also be simulated whenever the non-intercepted rainfall intensity falls below the saturated hydraulic conductivity of the soil. The redistribution option requires two more soil hydraulic parameters. Excess rainfall becomes surface runoff and is routed as overland flow and subsequently as channel flow. The overland flow routing formulation is based on a two-dimensional explicit finite difference (FD) technique, while two different FD techniques, one explicit and one implicit, provide options for routing one-dimensional channel flow. Through a step function, a depression depth may be specified, below which no overland flow will be routed.
The following sections describe various aspects of the model. In executing the model, the likelihood of making a serious mistake by an inexperienced user is unfortunately very high due to complexity of input and output options. Thus the user must thoroughly read all the details of this manual and then select appropriate options for his or her needs. Since at this time GRASS is an integer GIS, all input maps are stored in integer form. This requires conversion from integer to floating point representation via a scaling factor. The scaling factor for all floating point values is fixed. Therefore, the integer maps of parameter values are created by multiplying floating point values (e.g. volumetric water content) by a multiplier.
The following input/output parameters/options control complexity of the simulation. Map and file names in square brackets [] are optional. Some maps are mutually exclusive ( logical -or- ), while some maps require other maps to enable proper function ( logical -and- ). Carefully read the NOTES section.
TOPOGRAPHY
elevation=mapname map of elevation (DEM).
outlet_eastNnorthNslope=value,easting, northing, and bed slope at the outlet
value,value (comma delimited)
TIME
time_step=value computational time step duration in seconds. (typ. 1 to 30 seconds)
tot_time=value total simulation time in sec.
OVERLAND FLOW
[Manning_n=value] spatially uniform Manning's n roughness value for overland flow.
-or-
[roughness_map=mapname] spatially varied map of Manning's n roughness coefficient (values in 1000*Manning's n).
[watershed_mask=mapname] map of watershed boundary (or mask). This option is recommended, as it speeds execution greatly.
[d_thresh=value] threshold overland depth, in meters, below which overland routing will not be performed (i.e. average depression storage).
[initial_depth=mapname] map of initial overland (not lakes) depth in mm.
RAINFALL
rain_duration=value total rainfall duration in sec.
[unif_rain_int=value] spatially uniform rainfall intensity in mm/hr.
-or-
[r_gage_file=filename] raingage rainfall input file name (ASCII).
and-
[num_of_raingages=value] number of recording raingages.
-and-
[gage_time_step=value] time step (temporal resolution) of recorded raingage data in sec.
or-
[radar_intensity_map=mapname] prefix of time series of maps of radar- (or otherwise-) generated rainfall intensities in mm/hr.
and-
[radar_time_step=value] time increment between radar- (or otherwise-) generated rainfall maps in sec.
INTERCEPTION
[storage_capacity=mapname] map of vegetation storage capacity in tenths of mm.
-and-
[interception_coefficient=mapname] map of interception coefficient (values in 1000*actual coefficient).
INFILTRATION
[conductivity=mapname] map of soil saturated hydraulic conductivity in tenths of mm/hr (Req'd for G&A and Redist).
[capillary=mapname] map of soil capillary pressure head at the wetting front in tenths of mm (Req'd for G&A and Redist).
[porosity=mapname] map of soil effective porosity (values in 1000*porosity) (Req'd for G&A and Redist).
[moisture=mapname] map of initial soil moisture (values in 1000*moisture) (Req'd for G&A and Redist).
[pore_index=mapname] map of soil pore-size distribution index (Brooks & Corey lambda) in 1000*index (Req'd for Redist).
[residual_sat=mapname] map of soil residual saturation (values in 1000*residual saturation) (Req'd for Redist).
LAKES
[lake_map=mapname] map of lakes categories.
[lake_elev=mapname] map of lakes initial water surface elevation (also see unit_lake).
CHANNEL ROUTING
[channel_input=filename] channel input data file name (ASCII), required for explicit (EX) and implicit (IM) channel routing methods.
[links_map=mapname] map of channel network link numbers. (EX & IM)
[nodes_map=mapname] map of channel network node numbers. (EX & IM)
[table_input=filename] look-up table file for links with breakpoint cross section, link type 8, (ASCII) (IM)
[dis_profile=filename] channel initial discharge profile file name (ASCII). (IM)
[wat_surf_profile=filename] channel initial water surface profile file name (ASCII). (IM)
[hyd_location=filename] file name containing link and node numbers of internal locations where discharge hydrographs are to be saved (ASCII).
UNITS
[unit_el_conv=value] unit conversion factor by which the values in elevation must be DIVIDED to convert them into meters.
[unit_lake=value] unit conversion factor by which the values in lake surface elevation map must be DIVIDED to convert them into meters.
[unit_space=value] unit conversion factor by which all region easting and northing values must be DIVIDED to convert them into meters.
OUTPUT
discharge=filename outlet hydrograph file name (ASCII).
[dis_hyd_location=filename] output file name for discharge hydrograph at internal locations (ASCII)
[write_time_step=value] time increment for writing output raster maps in seconds.
-and-
[depth_map=mapname] output maps of surface depth in mm.
[inf_depth_map=mapname] output maps of cumulative infiltration depth in tenths of mm.
[surf_moist_map=mapname] output maps of surface soil moisture in number of fractions of a thousand.
[rate_of_infil_map=mapname] output maps of infiltration rate in mm/hr.
[dis_rain_map=mapname] output maps of distributed rainfall intensity in mm.
FLAGS
There are several flags whose utility is driven by data availability and/or user's choice.
-s do not check square tolerance of a grid cell. With this option non-square grid cell can be simulated using width of east-west resolution.
-m r_gage_file units in mm/hr. Without this option r_gage_file units in in/hr.
-t spatially interpolates raingage rainfall intensities using Thiessen polygon technique. The default technique uses inverse-distance-squared proportionality for interpolation of rainfall intensity over space.
-o routes edge-accumulated overland flow out of active region (ONLY when no mask is specified). Often the surface water accumulated at the edges of the current region creates severe backwater effects and may limit the use of longer computational time steps.
-e performs one-dimensional explicit finite difference channel routing. May be suitable for low- to medium-intensity rainstorms over small arid and semi-arid watersheds with no base flow discharge. This option often limits the computational time step to small values (<10 seconds). Channel bed smoothing is recommended to eliminate adverse slopes. No hydraulic structures, except reservoir spillways, can be simulated.
-p assumes uniform channel geometry in each link (requires -e option). If this option is chosen, the channel input file (channel_input) must have only one line per fluvial link rather than the default of one line per node.
-i performs Preissmann double sweep implicit channel routing. Particularly suitable for watersheds with some base flow to avoid dry-bed condition with channel slopes less than 1%. Supercritical slopes cannot be handled; a warning message will be printed if supercritical flow is encountered.
-b initializes the channel depth and discharge files (similar to -d, requires -i option) using the standard step backwater method. This option must be used with -i flag and replaces -d option to write "dis_profile" and "wat_surf_profile" files. At present, only link types 1, 2, and 8 may exist in the channel network.
-d performs channel initialization for implicit routing (similar to -b, requires -i option) by flooding the entire channel network with a horizontal water surface and draining down to normal depth using a y(t) outlet boundary condition (similar to -b). It is essential for implicit channel routing technique that a minimum initial base flow discharge exist in the channels and also the corresponding initial water surface profile at each node in the channel network have a realistic value. When the depth at the outlet reaches normal depth, the values of depth and discharge at each node is written to "dis_profile" and "wat_surf_profile" files for use in start up of actual simulations. With this option no other component of the model is executed; only the implicit channel routing is performed to create initial depth and discharge files necessary for start up of actual simulations.
-u writes discharges in cubic feet per second (cfs) and volumes in cubic feet in "discharge" file. The internal calculations are all in SI units regardless of this flag. The default SI unit prints the discharges in cubic meters per second (cms) and volumes in cubic meters (m3).
-q quietly skips printing iteration, time, and discharge values to the screen. No status report is printed.
The user must pay close attention to the following notes prior to any simulation:
1) watershed_mask: Although optional, preparation of this map is highly recommended as it cuts down on the memory requirements by the amount directly proportional to the ratio of mask area over the region area (also see -o flag). If the basin (mask) delineation has not been performed correctly, surface water may accumulate near the edges of the watershed. This excess water has no way out of the watershed and will accumulate, creating undesirable backwater effects within the watershed which eventually dictate use of shorter time step in order to accommodate such effects. It is recommended, in such instances, to re-delineate the watershed along the edges. An alternative would be to lower the elevation of the cells being flooded by excessive surface runoff near the edges such that no artificial backwater is created.
2) elevation: The elevation map in the form of Digital Elevation Model (DEM) is undoubtedly one of the most important inputs for distributed modeling. Thus the quality of the DEM plays a major role in success of distributed hydrologic simulations. DEMs almost always contain errors. Large flat areas in the DEM may be due to the limited vertical resolution of elevation data. Routing over such flat areas usually creates problems for the numerical techniques used in distributed physically-based models. Unreal pits in the DEM may be artifacts of interpolation scheme used to rasterize digitized contours, or due to coarse resolution in areas of concave topography. As a rule of thumb, the user must cross check the DEM with topographic maps of the area. One way to discover potential errors in the DEM is to run r.hydro.CASC2D iteratively with no other option except uniform rainfall and a relatively short time step. Writing surface depth maps should also be selected (see depth_map and write_time_step options). The simulation may terminate normally, at which case the surface depth maps must be examined in order to determine where most water accumulation occurred and whether such accumulation areas is justified by the topographic map of the area. Also stream network may be checked for proper connectivity via depth maps. Alternatively, if the model crashed at a certain location (whose address is printed on the screen as well as at the bottom of discharge file) the user must zoom in and double check the DEM with the topographic map. Often times some manual editing of the DEM is necessary in order to impose the actual drainage trend of the topo map. Prior to editing it is recommended that the DEM be multiplied by 10 or 100, particularly if the original vertical resolution of the DEM is also a concern. To account for this multiplication factor use unit_el_conv=10 or unit_el_conv=100, whichever appropriate. Note that only the unreal pits and depressions must be removed since they most likely trap surface runoff which would have otherwise contributed to the outlet. The real lakes and reservoirs may be simulated if delineated (see lake_map). In any event, non-smoothed DEMs require short computational time steps, while properly smoothed DEMs, particularly those with coarser grid space resolution, allow use of longer time steps. Another source of concern, which was briefly mentioned in above paragraph, about the quality of the DEMs is that a nicely connected stream network cannot be derived from non-smoothed DEM. If no delineated network (in the form of a vector file, for example) is available, an approach similar to what was described before may be taken in order to edit the DEM. However where a network has been delineated independent of the DEM then the elevation of stream cells should be checked so that they are not higher than those of the surrounding cells. Otherwise the stream cells will not properly collect the surface runoff.
3) initial_depth: This is a map which contains initial overland depth values, if any. Rarely used since prior to the storm overland planes are often dry. For channel initial depth see dis_profile and wat_surf_profile options.
4) storage_capacity & interception_coefficient: In r.hydro.CASC2D, the interception rate (i) is expressed as:
i(t)=r(t) while I < a
i(t)=b*r(t) while I > a
where r(t) denotes rainfall intensity at time t; a is the storage capacity; b is the interception coefficient; and I is the cumulative interception depth. Storage capacity map, as well as interception coefficient map, are usually a reclass of vegetation map. For table of storage capacity and interception coefficient values see Gray (1970, section 4.6) or Bras (1990) p. 233.
5) roughness_map: This map represents the spatial distribution of overland Manning roughness coefficient n. This map could be a reclass of vegetation cover map and/or the land use map. Tables of Manning's n are available in most hydrology textbooks. By using Manning resistance equation it is assumed that the overland flow over watersheds satisfies the conditions of turbulent flow over rough surfaces.
6) conductivity, capillary, porosity, & moisture: Four soil property maps are needed for modeling infiltration process using the Green-Ampt technique. These maps respectively contain saturated hydraulic conductivity, capillary suction head, effective porosity, and initial soil moisture content. The first three maps may be reclasses of soil texture map and the last one must be prepared considering antecedent soil moisture conditions. Based on soil textural classifications, tables of estimates of the first three parameters can be found in Rawls et al. (1983). Wherever reliable field measurements of such parameters are available they may be substituted for table values. Note that the saturated hydraulic conductivity map must be adjusted for the percent of rock fraction within the soil, if known. If these four maps are specified, Green & Ampt infiltration calculations are performed. If one or more of the maps is not specified, no infiltration calculations are made.
7) pore_index & residual_sat: These two maps are required when continuous soil moisture accounting is of primary interest. The model is capable of redistributing the soil moisture during periods of no- or low- intensity rainfall, over which infiltration capacity may recover for the next burst of storm intensity. The technique used herein for hiatus and post-hiatus stages is primarily based on the method by Smith el al. (1993). In this model Green-Ampt equation is used for post-hiatus stage. Tables of pore-size distribution index (Brooks & Corey lambda) and residual saturation are given by Rawls et al. (1982). As in (6), measured values should be substituted where available. When the four maps specified in (6) as well as the two maps specified here are given as command line arguments, the redistribution infiltration routine is used. All six maps must be specified to enable redistribution infiltration calculations.
8) lake_map & lake_elev: The first map is the delineated lake map with different categories corresponding to individual lakes and the second map holds the initial water surface elevation in the lakes. It is best to number the lakes categories sequentially. If a lake or reservoir connects to the implicit channel network at its outlet, then the lake must also be numbered as a link as reflected in the channel_input file. See channel_input and links_map options for more details. However isolated pond areas may be simulated as well; in fact such ponds are recommended to be delineated as part of the lake_map. Such isolated lakes should not be present in the links map since their outlet is not connected to the channel network. The routing of the lakes is performed by linear reservoir technique. Rainfall is added directly to the lake inflow. No interception or infiltration is abstracted from the lake cells.
9) radar_intensity_map: This is the common prefix of radar- or otherwise- generated rainfall maps. One possible source of this type of data could be the NEXRAD system (Crum and Albery, 1993). Alternatively, such time series of raster maps may come from the output of an interpolation software, which takes in the rainfall site data corresponding to different times and generates the time series of raster maps. For instance, where raingage rainfall data is available, one could run s.surf.rst of GRASS, one time for every recording time, and save the interpolated rainfall maps to be used as input for r.hydro.CASC2D. In any case one must pay attention to the unit of intensity which is in integer millimeters per hour for map values. Also see radar_time_step, and rain_duration options. Example: If one sets radar_intensity_map=rain.map and there are a total of eight maps in regular intervals of 10 minutes (radar_time_step=600 and rain_duration= 4800, both in seconds), then the actual name of the maps in the current mapset would have to be: "rain.map.1", "rain.map.2" through "rain.map.8", respectively corresponding to time periods of 0-10, 10-20 through 70-80 minutes.
10) links_map & nodes_map: A link is a channel segment, of finite length, which is comprised of two or more computational nodes placed at the center of each grid cell. Any internal boundary condition (weirs and lakes/reservoirs, for example) is also considered a link with only two nodes, one node upstream of the internal link and one downstream. Although all links and nodes must be numbered as discussed below, the link and node maps only contain trapezoidal or look-up table cross sections (link types 1 and 8). Internal boundary condition link types (e.g. weirs) do not appear in the link or node maps. The internal boundary condition link types do appear in the channel_input file (see attachment). The purpose of the link and node maps is to provide connectivity information between overland flow grid cells and channel (link,node) pairs. This information is used to pass overland flow to the 1-D channel model as lateral inflow. The nodes map must contain the node numbers corresponding to the links map. Thus nodes corresponding to internal boundary condition links are not present in the nodes map. At a cell where a link terminates and another link begins (junctions for example) the node number must comply with the downstream link; i.e. this cell is assigned a value of 1 in the nodes map. However one must note that the last node of the upstream link also lies at such junction cells. This last node number must be accounted for in the main ascii channel data file (channel_input), in which the total number of nodes per link must be given. An exception to this rule is the most downstream link leading to the outlet. The links map essentially describes topology of the stream network and its strict conformation with the stream numbering conventions is vital. Output from the GRASS command r.watershed cannot be directly used in conjunction with r.hydro.CASC2D. The link numbers assigned by r.watershed must be re-numbered in accordance with the following rules. Additionally, the node map may be constructed by re-numbering the link map. The general rules for link and node numbering are:
A small example network with three links is given below:
LINK MAP NODE MAP
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
0 0 0 1 1 0 0 0 0 0 2 0 0 0 0 0 3 2 0 0 0 0 0 1 0 0
0 0 0 1 0 0 0 0 0 2 2 0 0 0 0 0 4 0 0 0 0 0 3 2 0 0
0 0 0 1 1 0 0 0 2 2 0 0 0 0 0 0 5 6 0 0 0 5 4 0 0 0
0 0 0 0 1 0 0 2 2 0 0 0 0 0 0 0 0 7 0 0 7 6 0 0 0 0
0 0 0 0 1 1 3 2 0 0 0 0 0 0 0 0 0 8 9 1 8 0 0 0 0 0
0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0
0 0 0 0 0 0 3 3 0 0 0 0 0 0 0 0 0 0 0 3 4 0 0 0 0 0
0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0
This network has three fluvial links, link 1 has 9 nodes in the node map, while link 2 has 8 and link 3 has 5. Note that the junction of links 1 and 2 is labeled link 3. However, links 1 and 2 have hidden nodes at this junction, which must have the same bed elevation as the first node in link 3, but different cross-sections. The junction is really not the same point in space for all three links, but really represents a short distance from the confluence. In the channel_input file, link 1 really has 10 nodes and link 2 has 9. Link 3 would only have 5 nodes, because its downstream end is not a junction, rather it is the watershed outlet. The link and node maps tell r.hydro.CASC2D that the overland flow in the grid cell at row 5, column 4 is passed to link 1, node 5. For more discussion of links and nodes maps, see channel_input description and the last section of this document for channel input and table file options".
11) channel_input, table_file: For a detailed description of channel data file set-up and requirements: see the last section of this document which is a modified version of the document presented at the CASC2D Workshop by F.L. Ogden, June 9-10, 1994, at the University of Memphis, Tennessee (Revision 2, 10 January 1995). Also it is recommended for the user to refer to Ogden (1994) for an overview of the implicit channel routing formulation and numerical scheme. Note that for explicit channel routing, the same channel input file may be used provided that only link types 1 and 4 are present. If explicit channel routing is used, the channel talweg elevations are taken from the channel_input file, not from the DEM.
12) dis_profile & wat_surf_profile: These two files hold, respectively, the discharge profile and the water surface profile of all nodes within the stream network, including internal boundary condition nodes. These files are created by running r.hydro.CASC2D with the -i option AND either -d or -b. After creation, these files are used to initialize the implicit channel routing scheme for actual watershed simulations. They must be present for actual simulations. The -d flag creates initial discharge and depth files by flooding the entire channel network with a horizontal water surface, and draining the network using a depth vs. time relation at the watershed outlet. This draining proceeds until normal depth is reached at the outlet. At this time, the dis_profile and wat_surf_profile files are written for use in channel initialization for later simulations. The -b flag creates the initialization files using the standard step backwater method. Note that in dis_profile the unit of discharge is in cubic meters per second and the unit of depth in wat_surf_profile is in meters. Except for the reservoirs, the water surface profile file actually holds the flow depth values and thus is measured relative to the bed. CAUTION: The implicit channel routing routine must have good first approximation of depth and discharge at each node in the network. Inaccurate values (guesses) will surely lead to a crash. Additionally, the implicit channel routine is for SUBCRITICAL (Froude number < 1) flow only. A newly created channel input file will almost always crash when first run. It is then the task of the user to identify the location of the problem. Usually, there are regions of hydraulically steep slopes which cause supercritical flow. The channel code will warn the user of supercritical flow, including node and link number when it occurs, and then exits. The user should look at that node/link combination, and alter the data file to eliminate the reach of steep slopes. One workable solution to prevent stability problems at low flows is to use the so-called "Preissmann Slot" (see Cunge et al., 1980, for example). Be patient, as getting the channel network running is a time consuming process. There are a host of possible errors including: abrupt changes in cross-section, DEM-error, etc., which can cause problems. The implicit routing method is not applicable to steep (mountainous or upland) streams. If supercritical flow is encountered in upland (1st order) streams, they can be eliminated from the network. The user should verify the suitability of this approach on each watershed.
13) hyd_location: If a filename is specified, this input ascii file contains the link and node numbers of the locations at which discharge hydrographs are to be written to the file named in the dis_hyd_location option. The first column in this file should contain link numbers and the second column must be filled with the corresponding node numbers. For a description of the output, see the dis_hyd_location option.
14) r_gage_file: This is an ascii file which must be provided when raingage rainfall data is being simulated. At the top of the file, two-column lines hold the easting and the northing for each raingage. The number of such lines is determined by the total number of recording raingages (num_of_raingages). The location of any of the gages does not have to be within the current region nor within the current watershed mask as long as the easting and the northing are not specified relative to the current region, but are based on absolute values in the UTM or SP coordinates. If a gage falls outside in a different zone to the left of the active region's zone, then negative values are also acceptable. Note however that raingages well outside the watershed under analysis generally provide poor rainfall estimates. The subsequent lines at the bottom of the raingage file must reflect temporal variation of rainfall intensity. The number of columns per line is equal to the number of raingages (specified via num_of_raingages). The columns are separated by space. The number of lines in this lower portion is equal to number of instances, separated by a constant time interval, the raingages have made a recording. As usually is the case, the unit of rainfall intensity for raingage data must be in inches per hour. Example: For three raingages, each recording rainfall intensity every 2 minutes for a total duration of, say, 10 minutes, a file called "rain.inp" may look like this:
205150.0 750212.0
20545.0 750104.0
205320.0 750173.0
0.0 0.0 0.55
1.75 2.25 0.80
1.00 1.80 1.50
0.65 0.90 0.70
0.0 0.50 0.30
In above example, the eastings and northing of the first, second, and third raingages are (205150.0,750212.0), (20545.0,750104.0), and (205320.0,750173.0) respectively. The intensities recorded by the first gage, for example, are 0.0, 1.75, 1.00, 0.65, and 0.0 inches per hour, respectively, over 0-2, 2-4, 4-6, 6-8, and 8-10 minutes, etc. For this example r_gage_file=rain.inp, num_of_raingages=3, gage_time_step=120, and rain_duration=600. For state plane (SP) coordinate system the eastings and northing will have to be in feet rather than meters.
15) outlet_eastNnorthNslope: These three values determine the location of the outlet, in terms of its easting and northing, and the outlet bed slope. One needs to make sure that the outlet described by its easting and northing is not only within the active region but also inside watershed mask. Often times the region is not set to its original settings after zoom-in operations (d.zoom) are performed and this may put the outlet outside the active region thus causing the model to eventually crash. The bed slope is equal to tangent of the angle which is made between the bed profile at the outlet and horizontal plane. This slope is primarily used to calculate the outflow overland discharge at the outlet, if any, based on normal depth boundary condition, when no channel routing is performed (all surface flow treated as overland flow and channels essentially assumed wide). Normal depth is also assumed to prevail at the outlet when explicit channel routing (-e) has been selected.
16) Manning_n: The alternative to simulating spatially varied roughness coefficient is to provide a single value of Manning roughness coefficient n via Manning_n parameter. The user is warned if both roughness_map and Manning_n have been specified.
17) unif_rain_int: This option represents the intensity of the spatially- uniform temporally-constant (up to the rainfall duration) rainfall. Therefore this option may replace r_gage_file and radar_intensity_map options. The unit is in millimeters per hour.
18) num_of_raingages: The total number of recording raingage. If this variable is set to one, no spatial interpolation is performed and the rainfall is treated as uniform is space but could vary in time. The temporal variation will be then provided by r_gage_file. Note that when num_of_raingages is set to one, the easting and northing of the gage is irrelevant although two (arbitrary) values must still be provided at the top of raingage file (r_gage_file).
19) time_step: This represents the duration of computational time step in seconds and is a critical variable determining the total execution time for a particular simulation. These is no firm guide for the selection of the time step; it comes with experience and strongly depends on watershed and rainfall characteristics. The general rule for overland routing and explicit channel routing is that shorter time steps must be used for higher intensity storms, finer horizontal grid resolution (grid spacing), steeper watershed slopes, larger watershed areas, and smoother surfaces. Stability of explicit routing depends upon Courant number. Unfortunately, the critical condition for Courant number limits the length of the computational time step which must be used for the entire simulation unless variable time step algorithm is implemented in the future. Shorter time steps must be used when backwater effects are generated, mainly in flat areas which are not part of the lake map. If the time step is too long for any particular simulation the surface water depth in completely flat areas may show a checker-board pattern, i.e. oscillations are observed in the water surface level. This eventually results in a crash. As such, the time step should be decreased and the simulation repeated; or the flat areas be delineated within the lake map.
20) gage_time_step: This variable represents the time interval, in seconds, between recording instances of rainfall intensities by raingage(s). It is implied that the rainfall data has been recorded in regular intervals. See the notes under r_gage_file for an example.
21) radar_time_step: A value expressing the time interval between consecutive rainfall maps. A uniform time step is implied.
22) rain_duration: This is the total duration of rainfall in seconds. If multistorm events are simulated, the time from the beginning of the first storm to the end of the last storm constitutes the total rainfall duration. As such, selection of soil moisture redistribution capability is recommended via specifying pore index (pore_index) and residual saturation (residual_sat) maps.
23) tot_time: The total simulation time in seconds. If the falling limb of the discharge hydrograph is of particular interest the total simulation time must be set to a value greater than total rainfall duration plus the expected recession time.
24) write_time_step: This time parameter determines the frequency of writing output raster maps and is equal to the time interval, in seconds, at which output raster maps are saved. Also see depth_map, inf_depth_map, surf_moist_map, rate_of_infil_map, and dis_rain_map output options.
25) unit_el_conv: This conversion factor is used to convert elevation values in DEM (elevation) map to meters. This parameter doesn't need to be specified if the DEM is already in meters since the default is one. For DEM units in cm or ft, unit_el_conv must be respectively set to 100 or 3.281.
26) unit_lake: This conversion factor is used to convert the water surface elevation values in the lake_elev map to meters. The default is 1.0.
27) unit_space: This conversion factor is used to convert the horizontal grid spacing resolution to meters. The default is 1.0. For state plane coordinate system in ft, set unit_space equal to 3.281.
28) d_thresh: This parameter represents an average step retention storage below which the overland depth will not be routed. Another words, all depressions less than or equal to d_thresh are filled before any overland flow could begin. The unit must be in meters so for 2 mm of depression storage set d_thresh=0.002. Higher values of depression storage would reduce the total execution time of the model as the overland routing consumes most of the CPU effort.
29) discharge: The discharge hydrograph computed at the outlet will be saved in this ascii file under the current directory. Some other information, such as peak discharge, is also printed in this file.
30) dis_hyd_location: Whenever the hyd_location option is selected, the discharge at individual node/link pairs will be saved in this file. The discharge hydrographs at the (link,node) locations specified in hyd_location file are grouped in columns. The number of lines in this file is determined by the total number of iterations equal to tot_time divided by time_step.
31) depth_map, inf_depth_map, surf_moist_map, rate_of_infil_map, & dis_rain_map: depth_map is the common filename prefix given to the time series of output raster maps containing surface depth values in millimeters. At the channel cells the channel flow depth will be recorded in the maps. The surface depth maps, and all other output maps, are saved at regular intervals determined by write_time_step option. If write_time_step is not set, no output raster map will be saved. The first map always corresponds to the initial condition and naturally shows the water surface profile corresponding to the base flow discharge within the channel network when implicit channel routing is performed. Similarly the last depth map corresponds to the end-of-simulation time, or to the time at which the program finished abnormally, for example due to selection of a long step which generated oscillations leading to a crash. Abnormal program termination caused by oscillating depths may show negative depths in the overland plane. There is a hard-coded limit of 2000 output raster maps for each simulation. inf_depth_map and rate_of_infil_map are two options to save the output raster maps of cumulative infiltration depth, in tenth of mm, and rate of infiltration in mm/hr, respectively. These two options can only be selected when infiltration is being computed via either the Green & Ampt or redistribution methods. surf_moist_map output contains the soil moisture values at the soil surface and it may solely be selected with continuous infiltration option (pore_index and residual_sat maps are specified). dis_rain_map option saves the instantaneous rainfall intensities, in mm/hr, at write_time_step intervals and may be set for the simulations involving spatially distributed rainstorms; i.e. raingage rainfall data or radar (or other sources of rainfall raster maps) data. Example: If one sets depth_map=depth.out, tot_time=1000, and write_time_step= 100, then a total of 11 raster depth maps will be saved in a normally-terminated simulation: depth.out.00, depth.out.01, depth.out.02 through depth.out.10, respectively corresponding to times equal to 0, 100, 200, through 1000 seconds. These maps may be viewed sequentially in animation form using the GRASS animator program xganim.
The naming convention associated with the Preissmann double-sweep 1-D implicit channel routing method is based on the concept to links and nodes. A link is a channel segment, or an internal boundary condition, which is comprised of two or more computational nodes. All internal boundary conditions contain two nodes, while fluvial reaches may be of any size greater than or equal to two nodes. The following discussion of input file format must first distinguish between different link types. A few of the possible link types are presented below in Table 1. As of August 1995, only link types 1, 2, 4, and 8 are supported. Development is continuing on link types 3 and 7.
Table 1. Link Types
Link Description Number of Number of Type of Link Type Parameters Nodes 1 Fluvial Link, Trapezoidal Cross-Section 5 >=2 2 Overflow Weir 8 2 3 Culvert/Weir (not yet supported) 8 2 4 Reservoir 5 2 7 Bridge Crossing (not yet supported) 8 2 8 Fluvial Link, Look-up Table (Breakpoint) 4 >=2 Cross-SectionAt present the channel model formulation accepts cross-sectional input for only two different channel geometries, namely trapezoidal, and breakpoint via look-up table. Trapezoidal channel parameters which must be provided as input at each computation node include: Manning's n, bottom width, channel depth, side slope, and bed elevation. Look-up tables of cross sectional properties must include cross-sectional area, top width, and conveyance at equal depth intervals. Smooth transition in channel cross sectional properties within links of type 8 and between all connecting fluvial links often plays a vital role in the success of simulations. Abrupt changes in cross-sections can lead to numerical errors in mass conservation.
As far as handling the reservoirs (link type 4), flow is not routed through reservoirs in this version of the CASC2D channel routing code. Instead the linear reservoir approximation is used. Among other internal boundary conditions link types, only weirs can be simulated at present.
The channel_input file contains the channel bed elevation at each node, which constructs the longitudinal profile of channels in the drainage network. Ideally, the modeler should use surveyed cross sections and talweg profiles of the channel network. However, extensive surveys are often impossible for the entire drainage networks on large watersheds. In lieu of an extensive survey, there are existing tools for extracting channel topology from digital elevation models (DEM). However, if the channel network is extracted from a DEM, a smoothing algorithm must be applied (e.g. Ogden et al., 1994) to produce physically realistic longitudinal profiles because of errors inherent in any DEM.
EXAMPLE OF INPUT FILE FORMAT AND CONSTRUCTION
NOTE: The channel_input file and table_file both follow SI unit convention. All units of length are in meters, including elevations.
In this example, a channel input data file is constructed for a fictitious watershed. Note that this particular example is merely for demonstration and has not been tried in any simulation run. The stream network will consist of trapezoidal and look-up table fluvial links, an internal reservoir, and an overflow weir. Regarding the MASK map (watershed_mask option), note that all watershed cells must be marked with 1, while all raster cells outside the watershed boundary are marked with 0. Also all the cells in the lake should carry the value 1 on the LAKE map (lake_map option). If there was another lake in this example, it would be denoted with the number 2, etc.
The channel link numbering scheme typically employed in double-sweep routing was explained in the main text. With reference to the LINK map shown below, (links_map option), assume that links 1 and 2 drain into link 3, which in-turn drains in to the lake, which is assigned link number 4. Also links 7 and 8 drain into link 9. Link 5 originates at the outlet of the reservoir. Link 5 is split into link 6 because of a change in some cross-section property. Links 6 and 9 flow into link 10, which is immediately upstream from a weir which must be considered link 11. Finally, link 12, the most downstream link, lies below the weir (link 11). Note that the LINK map contains continuous sequence of link numbers, except for the internal boundary conditions such as the weir (link 11) and lakes (link 4) which does not appear in the LINK map.
EXAMPLE LINK MAP
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0
0 0 0 1 0 0 0 0 0 0 0 0 0 2 2 0 0 0 0 0
0 0 0 1 1 0 0 0 0 0 0 0 2 2 0 0 0 0 0 0
0 0 0 0 1 1 1 0 0 0 0 0 2 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 0 0 0 2 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 0 2 2 2 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 3 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 7 7 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 7 7 0 0 0 0 5 0 0 0 0 0 0 0 0 0
0 0 0 0 0 7 7 0 0 0 5 5 0 0 0 0 0 0 0 0
0 8 8 8 0 0 7 0 0 0 0 5 6 0 0 0 0 0 0 0
0 0 0 8 8 8 9 9 0 0 0 0 6 6 0 0 0 0 0 0
0 0 0 0 0 0 0 9 9 9 0 0 0 6 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 9 9 9 9 10 10 10 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 12 12 12 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 0
Now, refer to the NODE map (nodes_map option) shown below. In the NODE map, each link is numbered from 1 to the number of grid cells spanned by that link, with the exception of the internal boundary conditions. Conceptually, internal boundary conditions (including reservoirs) have two nodes but they are not given any node numbers in the NODE map. This is how the program recognizes internal boundary conditions. Furthermore, all links except the one leading to the watershed outlet must have an extra node to provide connectivity to the downstream link. Therefore, even though the NODE map may show link 1 to have 13 nodes, it actually has 14. This implied extra node for link 1, shared between the most downstream node of link 1 and the most upstream node of link 3, provides the connection between link 1 and link 3. The number of nodes in each link in this example are shown below in Table 2. Note that the node map entries for the lake (link 4) are 0.
EXAMPLE NODE MAP
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 2 0 0 0 0 0 0 0 0 0 3 2 0 0 0 0 0
0 0 0 3 4 0 0 0 0 0 0 0 5 4 0 0 0 0 0 0
0 0 0 0 5 6 7 0 0 0 0 0 6 0 0 0 0 0 0 0
0 0 0 0 0 0 8 9 10 0 0 0 7 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 11 0 10 9 8 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 12 13 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 2 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 4 5 0 0 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 6 7 0 0 0 2 3 0 0 0 0 0 0 0 0
0 1 2 3 0 0 8 0 0 0 0 4 1 0 0 0 0 0 0 0
0 0 0 4 5 6 1 2 0 0 0 0 2 3 0 0 0 0 0 0
0 0 0 0 0 0 0 3 4 5 0 0 0 4 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 6 7 8 9 1 2 3 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 3 4 5 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0
Table 2. Number of Nodes in each Link for Example Watershed Stream Network
Link Number Number of Nodes Number of Nodes
as in nodes_map in channel-input file
1 13 14
2 10 11
3 3 4
4 0 2
(reservoir)
5 4 5
6 4 5
7 8 9
8 6 7
9 9 10
10 3 4
11 0 2
(weir)
12 6 6
(outlet link)
The first portion of the channel input file is used to pass physical constants and simulation parameters to the model. These include the gravitational acceleration "g", kinetic energy correction factor "alpha", the friction slope weighting factor "beta", the spatial derivative weighting coefficient "theta", the length of each node "dx", the computational time step "dt" (seconds), the total simulation time "tt" (seconds), and the discharge in all first order streams "qmin" (m3/s). At present, the time step "dt" must be identical to the computational time step used in the overland flow routing portion of r.hydro.CASC2D. The program must be told the number of links, and the largest number of nodes in any link in the network for dynamic memory allocation. In our example problem, the number of links is 12, and the maximum number of nodes is 14 (in link 1). Remember that all links which are not at the outlet or not immediately upstream of a reservoir must have an extra node for connectivity purposes. The total number of links is called "nlinks", and the largest number of nodes in any link in the network is called "maxnodes". In this example, we will use the constants and parameters given in Table 3:
Table 3. physical constants and simulation parameters
"g" 9.81 m/s2
"alpha" 1.0
"beta" 0.5
"theta" 0.55
"dx" 100.0 m
"dt" 30.0 s
"tt" 3600.0 s
"qmin" 0.07 cms
"nlinks" 12
"maxnodes" 14
This data constitutes the first portion of channel_input file and is arranged into a header which must have the form (note floating point and integers):
9.81
1.0 0.5 0.55
100.0
30.0 3600.0
0.07
12
14
The second portion of the input file describes link types as well as the network topology and connectivity. This is accomplished using a line-input format, one line for each link in the network. Each line contains 6 values arranged in columns, as shown in Table 4.
Table 4. Connectivity information format
Column Data Description
1 Link number (INT)
2 Link type (INT)
3 Number of links upstream from this link (max. 2, min. 0) (INT)
4 Upstream link #1 (INT) (0 if no upstream links)
5 Upstream link #2 (INT) (0 if no or only one upstream link)
6 Downstream link (0 for outlet ) (INT)
As far as link types in the current example, assume that links 1 through 9, except reservoir link 4 (link type 4), are well-described as trapezoidal (type 1), and we are using look-up table data to describe the cross-sections of links 10 and 12 (type 8). The weir (link 11) is of link type 2. Thus the second portion of input file for describing link types and connectivity looks like:
1 1 0 0 0 3
2 1 0 0 0 3
3 1 2 1 2 4
4 4 1 3 0 5
5 1 1 4 0 6
6 1 1 5 0 10
7 1 0 0 0 9
8 1 0 0 0 9
9 1 2 7 8 10
10 8 2 6 9 11
11 2 1 10 0 12
12 8 1 11 0 0
Note that only the outlet link has 0 as a downstream dependency. It is important that this be the only link with 0 as a downstream dependency. Also, the location and topology of internal boundary condition links (4 and 11) with respect to other links are known via above connectivity information in the channel_input file. As an interpretation, examine link number 10 above. It is of link type 8 (look-up table cross-section data), has 2 upstream dependencies (links 6 and 9), and flows into link 11.
The third portion of the input file contains the individual hydraulic property information for each node in the network. Assume for now that all the trapezoidal channels in the network have the properties shown in Table 5.
Table 5. Cross-Sectional Properties of Example Trapezoidal Channels
Manning roughness coefficient varies
Bottom width (m) varies
Channel depth (m) 1.75
Side slope H:V varies
Talweg elevation dependent upon location
Now, build the input file section for link #1 which has 13 nodes on the NODE map, plus an extra for connectivity at the upstream end of link 3 (total of 14 nodes). This input file section will look like:
1 14
0.035 2.10 1.75 2.00 264.40
0.035 2.10 1.75 2.00 263.85
0.035 2.10 1.75 2.00 263.41
0.035 2.10 1.75 2.00 262.98
0.035 2.10 1.75 2.00 262.75
0.035 2.10 1.75 2.00 262.54
0.035 2.10 1.75 2.00 262.37
0.035 2.10 1.75 2.05 262.02
0.035 2.20 1.75 2.05 261.87
0.035 2.20 1.75 2.10 261.75
0.035 2.20 1.75 2.10 261.67
0.035 2.20 1.75 2.15 261.52
0.035 2.30 1.75 2.15 261.25
0.035 2.40 1.75 2.15 261.00
The 1 and 14 on the first line indicate that it is link 1, with 14 nodes. The columns represent Manning's n, bottom width, channel depth, trapezoidal side slope H/V, and talweg elevation, respectively.
Link 2 might have an entry which looks like:
2 11
0.035 1.80 1.75 1.80 264.21
0.035 1.80 1.75 1.80 264.10
0.035 1.90 1.75 1.90 263.81
0.035 1.90 1.75 1.90 263.56
0.035 2.00 1.75 1.90 263.34
0.035 2.10 1.75 1.90 262.86
0.035 2.10 1.75 2.00 262.24
0.035 2.10 1.75 2.00 261.76
0.035 2.20 1.75 2.10 261.48
0.035 2.30 1.75 2.10 261.26
0.035 2.30 1.75 2.10 261.00
The cross-section entry for link 3 would then look like:
3 4
0.035 2.60 1.75 2.20 261.00
0.035 2.60 1.75 2.40 260.20
0.035 2.80 1.75 2.90 259.67
0.035 2.80 1.75 2.90 259.17
The elevation of the downstream end of links 1 and 2 must be equal to the bed elevation of the upstream end of link 3. It is required that the bed elevation of all channel inverts at each junction be equal. Even though they are not exactly the same points in space, they are assumed close enough to have the same elevation.
Flow Flow
\ /
\ \ \ / / /
\ v \ / v /
\ \ / /
\ link 1 \ /link 2 /
\ node 14\ / node 11/
\ /
\ /
| |
| link 3 |
| node 1 |
The input for other links with trapezoidal cross section is similar to the input for link 1 (see the next section for the complete input file).
Additionally, assume that the reservoir (link 4) has a surface area of 0.498 square km, a rectangular spillway width of 12m, a spillway discharge coefficient of 0.97, an initial water surface elevation of 260.10m, and a spillway crest elevation of 260.50m. The resulting input file entry is shown below.
4 2
0.498 12.0 0.97 260.10 260.50
0.000 0.0 0.00 0.00 0.00
The number of nodes (line entries) in the channel_input file per reservoir link is two, but the numbers entered in the second row entry are for further improvements and for now their value is irrelevant. Also note that the bed elevation at the downstream limit of link 3 MUST be lower than the initial water surface elevation in the reservoir (link 4). This is mandatory at all downstream boundaries between channels flowing into the reservoirs. Reservoirs serve as downstream boundary conditions at upstream links. The elevation of reservoirs should therefore never be allowed to fall below critical depth in any channels which flow into the reservoir. If this happens, the code will crash.
The input for links 5, 6, 7, 8, and 9 will be similar to links 1, 2, and 3. See the constructed input file at the end of this text. Trapezoidal channel links are simple in terms of input. The channel depth field is intended to represent the average bank-full depth of the channel. This number is not used in the present version of the code, but will be in the future when the overland flood plane and channel flows are coupled iteratively. At present the flows are not fully coupled. Water can flow only from the overland flow plane into the channels, not from the channels back to the flood plane.
If look-up table data are used (link type 8), the look-up tables are stored in a separate file. The GRASS command line option for this look-up table file is table_file. If there are any link type 8 in the channel_input data file, the program will attempt to open the file table_file (typically called "table.dat"). If this file is not present, the program will terminate. In channel_input file, all that is needed is the table number for a particular cross section, and the talweg elevation. Assume that the 4 nodes in link 10 are approximated by one cross-section, which is given as table entry 1, then the channel_input file entry for link 10 would look like:
10 4
1 244.0 1 243.5 1 243.0 1 242.5
In this format, the first column represents the look-up table number for the given node, and the second column represents the bed elevation for that node. Table numbers can be mixed within a link. However, abrupt changes in cross-section should be avoided because they cause significant continuity errors in the formulation. The format of the table_file is discussed later in this document.
Link 11 is an overflow weir link. While this link has no nodes which appear in the node map for overland flow connectivity, it does have two nodes in terms of channel topology. The meaning of the columns in the weir node input section is presented in Table 6.
Table 6. Parameter Description for Weirs
Line Entry Column Parameter Description
1 1 forward direction weir discharge coefficient
1 2 reverse direction weir discharge coefficient
1 3 0.0 (reserved for future use)
1 4 crest length, meters
1 5 elevation of weir crest, meters
2 1-4 0.0 (reserved for future use)
2 5 downstream bed elevation, meters
The first (upstream) node represents the crest of the rectangular weir. The second (downstream) node represents the bed of the channel just downstream from the weir. Assume that this weir has a crest elevation of 244.4 m, a crest width of 8 m, a discharge coefficient in the forward direction of 0.92, a discharge coefficient for flow in the reverse direction equal to 0.85, and the channel bed elevation immediately downstream from the weir is 239.8m. The entry for this weir in the channel_input file would thus appear as:
11 2
0.92 0.85 0.00 8.00 244.4
0.00 0.00 0.00 0.00 239.8
Link 12 uses table entry 2 for the its first four nodes and entry 4 for the remaining two nodes. Note that link 12 has 6 nodes in the node map, and 6 nodes in the channel_input file. Because link 12 is the outlet link, we do not add an extra node at the downstream end for connectivity purposes. The portion of the channel_input file for link 12 looks like:
12 6
2 239.80
2 239.47
2 239.33
2 239.05
3 238.54
3 238.44
Hence, 238.44 is the talweg elevation at the outlet of the catchment. Also notice that the upstream end of link 12 and the downstream node of link 11 have the same bed elevation.
Weirs cannot be used as the outlet boundary condition for this channel model. If your watershed has a weir at the outlet, just place a short link downstream from the weir with appropriate characteristics.
FINAL CHANNEL FILE FORMAT
The example channel data file (typically called chn.dat) developed for this example looks like:
----------------------BEGIN channel_input FILE HERE----------------------------
9.81
1.0 0.5 0.55
100.0
30.0 3600.0
0.07
12
14
1 1 0 0 0 3
2 1 0 0 0 3
3 1 2 1 2 4
4 4 1 3 0 5
5 1 1 4 0 6
6 1 1 5 0 10
7 1 0 0 0 9
8 1 0 0 0 9
9 1 2 7 8 10
10 8 2 6 9 11
11 2 1 10 0 12
12 8 1 11 0 0
1 14
0.035 2.10 1.75 2.00 264.40
0.035 2.10 1.75 2.00 263.85
0.035 2.10 1.75 2.00 263.41
0.035 2.10 1.75 2.00 262.98
0.035 2.10 1.75 2.00 262.75
0.035 2.10 1.75 2.00 262.54
0.035 2.10 1.75 2.00 262.37
0.035 2.10 1.75 2.05 262.02
0.035 2.20 1.75 2.05 261.87
0.035 2.20 1.75 2.10 261.75
0.035 2.20 1.75 2.10 261.67
0.035 2.20 1.75 2.15 261.52
0.035 2.30 1.75 2.15 261.25
0.035 2.40 1.75 2.15 261.00
2 11
0.035 1.80 1.75 1.80 264.21
0.035 1.80 1.75 1.80 264.10
0.035 1.90 1.75 1.90 263.81
0.035 1.90 1.75 1.90 263.56
0.035 2.00 1.75 1.90 263.34
0.035 2.10 1.75 1.90 262.86
0.035 2.10 1.75 2.00 262.24
0.035 2.10 1.75 2.00 261.76
0.035 2.20 1.75 2.10 261.48
0.035 2.30 1.75 2.10 261.26
0.035 2.30 1.75 2.10 261.00
3 4
0.035 2.60 1.75 2.20 261.00
0.035 2.60 1.75 2.40 260.20
0.035 2.80 1.75 2.90 259.67
0.035 2.80 1.75 2.90 259.17
4 2
0.498 12.0 0.97 260.10 260.50
0.000 0.0 0.00 0.00 0.00
5 5
0.035 2.50 1.75 2.10 249.20
0.035 2.50 1.75 2.10 248.65
0.035 2.50 1.75 2.15 247.94
0.035 2.50 1.75 2.20 247.31
0.035 2.50 1.75 2.20 246.97
0.035 2.50 1.75 2.25 246.42
6 5
0.035 2.50 1.75 2.00 246.42
0.035 2.50 1.75 2.00 245.97
0.035 2.50 1.75 2.00 245.24
0.035 2.50 1.75 2.00 244.73
0.035 2.50 1.75 2.00 244.00
7 9
0.035 1.75 1.75 2.00 254.20
0.035 1.75 1.75 2.00 253.55
0.035 1.75 1.75 2.00 252.92
0.035 1.90 1.75 2.00 252.29
0.035 1.90 1.75 2.00 251.71
0.035 1.90 1.75 2.00 251.24
0.035 2.00 1.75 2.00 250.82
0.035 2.00 1.75 2.00 250.34
0.035 2.00 1.75 2.00 250.08
8 7
0.035 1.60 1.75 2.00 253.71
0.035 1.90 1.75 2.00 253.21
0.035 1.90 1.75 2.00 252.94
0.035 2.00 1.75 2.00 252.55
0.035 2.10 1.75 2.00 251.80
0.035 2.10 1.75 2.00 251.11
0.035 2.10 1.75 2.00 250.08
9 10
0.035 3.10 1.75 2.90 250.08
0.035 3.20 1.75 2.90 249.16
0.035 3.30 1.75 3.00 248.63
0.035 3.40 1.75 3.00 248.02
0.035 3.40 1.75 3.10 247.33
0.035 3.20 1.75 3.10 246.65
0.035 3.10 1.75 3.10 246.05
0.035 3.20 1.75 3.10 245.88
0.035 3.20 1.75 3.10 245.52
0.035 3.50 1.75 3.10 245.24
0.035 3.70 1.75 3.40 244.74
0.035 3.50 1.75 3.40 244.51
0.035 3.10 1.75 3.40 244.00
10 4
1 244.00
1 243.50
1 243.00
1 242.50
11 2
0.92 0.85 0.00 8.00 244.40
0.00 0.00 0.00 0.00 239.80
12 6
2 239.80
2 239.47
2 239.33
2 239.05
3 238.54
3 238.44
--------------------- END OF channel_input FILE ---------------------------
TABLE_FILE INPUT FORMAT
If link type 8 (look-up table x-sections) are present in the channel_input file, we require a file typically called "table.dat", to store all look-up table values. This file must begin with an integer equal to the number of tables contained within the file. In our example this number is 3. Then we require 3 tables. The first line of each table is an integer equal to the table number. The second line in the table entry contains two numbers, an integer, and a floating point (real) value. The first number is equal to the number of entries (rows) in each particular table, designated by "numhts". The second number is the vertical distance between hydraulic property points, which must be a constant for each table. For instance, if you describe the variation of cross-sectional area, topwidth, and conveyance at 0.5 m vertical intervals, this number will be 0.5.
Each table entry must then contain "numhts" lines, each with four entries. The first line must be:
0.0 0.0 0.0 0.0
subsequent lines must contain the following:
height topwidth x-section_area x-section_conveyance.
For instance, assume we know geometric variables for table entries 1, 2, and 3. The file "table.dat" (table_file option) may look like:
3
1 5 1.0
0.0 0.0 0.0 0.0
1.0 2.1 5.9 13.5
2.0 6.9 22.3 122.2
3.0 14.7 72.1 312.8
4.0 19.6 145.6 789.4
2 8 0.5
0.0 0.0 0.0 0.0
0.5 1.1 1.5 7.5
1.0 1.7 3.4 17.1
1.5 3.1 11.2 43.9
2.0 5.2 32.4 98.5
2.5 7.2 49.2 187.4
3.0 9.2 86.4 312.5
3.5 12.1 143.1 624.9
3 5 1.0
0.0 0.0 0.0 0.0
1.0 2.2 5.4 15.5
2.0 6.6 15.3 132.2
3.0 15.7 81.8 352.8
4.0 21.4 155.2 839.6
For further clarification, consider table 2 above, corresponding to cross sectional properties of the four most upstream nodes of link 12. Table 2 has 8 entries, each separated by 0.5 meter of depth. At a depth of 3.0 m, for example, the channel has a top-width of 9.2m, a cross-sectional area of 86.4 m2, and a conveyance of 312.5 m3/s. The conveyance K is used to calculate the discharge using:
Q=K*sqrt(Slope)
Where using Manning's equation, K is defined (in SI units) as:
K = 1/n * (Area) * (Hydraulic_radius)^(2/3)
These tables can be produced easily using a spreadsheet using measured channel cross-section data.
Bras, R. L., 1990, Hydrology: An introduction to hydrologic science, Addison-Wesley, Reading, Mass., 643 p.
Crum, T. D., and R. L. Alberty, 1993, The WSR-88D and the WSR-88D operational support facility, Bulletin of the American Meteorological Soc., 74(9), pp. 1669-1687.
Cunge, J.A., F.M. Holly, and A. Verwey, 1980, Practical Aspects of Computational River Hydraulics, Iowa Institute of Hydraulic Research, 404HL The University of Iowa, Iowa City, IA 52242. 420 p.
Gray, D.M., 1970, Handbook on the Principles of Hydrology, National Research Council of Canada, Water Information Center Inc., Water Research Building, Manhasset Isle, Port Washington, N.Y., 11050.
Julien, P. Y., and B. Saghafian, 1991, CASC2D users manual - A two dimensional watershed rainfall-runoff model, Civil Engr. Report, CER90-91PYJ-BS-12, Colorado State University, Fort Collins, CO.
Julien, P. Y., Saghafian, B., and F. L. Ogden, 1995, "Raster-Based Hydrologic Modeling of Spatially-Varied Surface Runoff", Water Resources Bulletin, AWRA, 31(3), pp. 523-536.
Ogden, F.L., 1994, de-St Venant channel routing in distributed hydrologic modeling., Proc. Hydraulic Engineering `94, ASCE Hydraulics Specialty Conference, G.V. Cotroneo and R.R. Rumer, eds., Vol. 1, pp. 492-496.
Ogden, F.L., Saghafian, B., and W.F. Krajewski, 1994, GIS-based channel extraction and smoothing algorithm for distributed hydrologic modeling, Proc. Hydraulic Engineering `94, ASCE Hydraulics Specialty Conference, G.V. Cotroneo and R.R. Rumer eds., August 1-5, 1994, Buffalo, N.Y., pp. 237-241.
Rawls, W. J., Brakensiek, D. L., and N. Miller, 1983, Green-Ampt infiltration parameters from soils data, J. of Hydraulic Engineering, ASCE, 109(1), pp. 62-70.
Rawls, W. J., Brakensiek, D. L., and K. E. Saxton, 1982, Estimation of soil water properties, Trans. of ASAE, pp. 1316-1320.
Saghafian, B., 1992, Hydrologic analysis of watershed response to spatially varied infiltration, Ph.D. Dissertation, Civil Engr. Dept., Colorado State University, Fort Collins, CO.
Saghafian, B., 1993, Implementation of a distributed hydrologic model within Geographic Resources Analysis Support System (GRASS), Proceedings of the Second International Conference on Integrating Environmental Models and GIS, Breckenridge, CO.
Smith, R. E., Corradini, C., and F. Melone, 1993, Modeling infiltration for multistorm runoff events, Water Resources Research, 29(1), pp. 133-144.
B. Saghafian and F. L. Ogden, Colorado State University