NAME r.hydro.CASC2D - GRASS raster command to execute fully integrated distributed hydrologic modeling. SYNOPSIS r.hydro.CASC2D, r.hydro.CASC2D help, r.hydro.CASC2D [-toepidbuq] [watershed_mask=name] elevation=name [initial_depth=name] [storage_capacity=name] [interception_coefficient=name] [roughness_map=name] [conductivity=name] [capillary=name] [porosity=name] [moisture=name] [pore_index=name] [residual_sat=name] [lake_map=name] [lake_elev=name] [radar_intensity_map=name] [links_map=name] [nodes_map=name] [channel_input=name] [table_input=name] [dis_profile=name] [wat_surf_profile=name] [hyd_location=name] [r_gage_file=name] outlet_east&north&slope=east,north,bedslope [Manning_n=value] [unif_rain_int=value] [num_of_raingages=value] time_step=value [gage_time_step=value] [radar_time_step=value] rain_duration=value tot_time=value [write_time_step=value] [unit_el_conv=value] [unit_lake=value] [unit_space=value] [d_thresh=value] discharge=name [dis_hyd_location=name] [depth_map=name] [inf_depth_map=name] [surf_moist_map=name] [rate_of_infil_map=name] [dis_rain_map=name] DISCLAIMER: The r.hydro.CASC2D model is Copyleft 1994 by B. Saghafian and F. L. Ogden. There are absolutely no warranties, expressed or implied, made regarding the accuracy of the results produced by this program. Use of this program is allowed as long as the original developers are properly acknowledged. No part of this document describing the use of the model may be used or duplicated without a complete citation. 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 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 notifying potential users. HISTORY: CASC2D originally started 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 while 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 and The University of Iowa, now Asst. Professor, Deptartment of Civil and Environmental Engineering, University of Connecticut, Storrs, Connecticut (email: ogden@eng2.uconn.edu). This version became known as r.hydro.CASC2D, part of the GRASS GIS for hydrologic simulations (Saghafian, 1993). DESCRIPTION: r.hydro.CASC2D simulates 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 being retained by vegetation and in this model is determined by specifying two parameters. Infiltration process causes surface water to soak into the permeable ground. 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, an explicit and an 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. Of particular importance are units of the input and output parameters. PARAMETERS/OPTIONS: The following input/output parameters/options control complexity of the simulation. Carefully read the NOTES section. watershed_mask=name map of watershed boundary (or mask). elevation=name map of elevation (DEM). initial_depth=name map of initial overland (not lakes) depth in mm. storage_capacity=name map of vegetation storage capacity in tenths of mm. interception_coefficient=name map of interception coefficient (values in 1000*actual coefficient). roughness_map=name map of surface roughness coefficient (values in 1000*Manning's n). conductivity=name map of soil saturated hydraulic conductivity in tenths of mm/hr. capillary=name map of soil capillary pressure head at the wetting front in tenths of mm. porosity=name map of soil effective porosity (values in 1000*porosity). moisture=name map of initial soil moisture (values in 1000*moisture). pore_index=name map of soil pore-size distribution index (Brooks & Corey lambda) in 1000*index. residual_sat=name map of soil residual saturation (values in 1000*residual saturation). lake_map=name map of lakes categories. lake_elev=name map of lakes initial water surface elevation (also see unit_lake). radar_intensity_map=name map of radar- (or otherwise-) generated time series of rainfall intensity in mm/hr. links_map=name map of channel network link numbers. nodes_map=name map of channel network node numbers. channel_input=name channel input data file name (ASCII). table_input=name look-up table file for links with breakpoint cross section (ASCII) (must use with -i option) dis_profile=name channel initial discharge profile file name (ASCII). wat_surf_profile=name channel initial water surface profile file name (ASCII). hyd_location=name file name containing link and node addresses of internal locations where discharge hydrographs are to be saved (ASCII). r_gage_file=name raingage rainfall input file name (ASCII). outlet_east&north&slope= easting, northing, and bed slope at the outlet value,value,value cell. Manning_n=value spatially uniform overland Manning's n roughness value. unif_rain_int=value spatially uniform rainfall intensity in mm/hr. num_of_raingages=value number of recording raingages. time_step=value computational time step duration in sec. gage_time_step=value time step of recorded raingage data in sec. radar_time_step=value time increment between radar- (or otherwise) generated rainfall maps in sec. rain_duration=value total rainfall duration in sec. tot_time=value total simulation time in sec. write_time_step=value time increment for writing output raster maps in sec. 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 the region easting and northing values must be DIVIDED to convert them into meters. d_thresh=value threshold overland depth, in meters, below which overland routing will not be performed (i.e. average depression storage). discharge=name outflow discharge output file name (ASCII). dis_hyd_location=name output file name for discharge hydrograph at internal locations (ASCII) depth_map=name output map of surface depth in mm. inf_depth_map=name output map of cumulative infiltration depth in tenth of mm. surf_moist_map=name output map of surface soil moisture in number of fractions of a thousand. rate_of_infil_map=name output map of infiltration rate in mm/hr. dis_rain_map=name output map of distributed rainfall intensity in mm. FLAGS: There are several flags whose utility is driven by data availability and/or user's choice. -t interpolates raingage rainfall intensities using Thiessen polygon technique. The default technique uses inverse-squared-distance 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 semi-arid regions with no to very low base flow discharge. Use of this option often limits the computational time step to small values. Channel bed smoothing is recommended. No hydraulic structures, except reservoir spillways, can be simulated. -p assumes uniform channel geometry in each link (needs -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. Supercritical slopes cannot be handled; a warning message will be printed if supercritical flow is encountered. -d performs only drainage of the basin for implicit routing (needs -i option). It is essential for implicit channel routing technique that a minmum 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. This option floods the watershed completely with water, and drains the water using a downstream depth vs. time relation. When the depth at the outlet reaches normal depth, the values of depth and discharge at each nodes are 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 drain the watershed. -b determines initial flow depths for implicit routing by performing standard step backwater computations. 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. -u prints 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 skips printing iteration, time, and discharge values to the screen. No status report is printed. IMPORTANT NOTES 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 thus it piles up and may create undesirable backwater effects within the watershed and eventually dictate use of shorter time step in order to accommodate such effects. It is recommended, in such instances, to try 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: 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. The DEMs almost always contain errors. Depending on the source of the DEM these errors may be attributed to various causes. Large flat area in the DEM may be due to the limited vertical resolution of elevation data. Routing over such flat areas sometimes create 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. As a rule of thumb, the user must cross check the DEM with topographic maps of the area. One way to find out where the potential sources of error in the DEM may be is to iteratively run r.hydro.CASC2D 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_con=10 or unit_el_con=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 values see Bras (1990) p. 233. 5) roughness_map: This map represents the spatial variation 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 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 in accordance with antecedent soil moisture condition. Based on soil texture classifications, tables of the first three parameters can be found in Rawls et al. (1983). Wherever reliable field measurements of such parameters are available they may substitute for the tables. Note that the conductivity map must be adjusted for the percent of rock fraction within the soil, if known. 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 paper 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). 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 has an outlet stream then the lake must also be numbered as a link in the links_map. See links_map 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 they have no outlet. The routing of the lakes is performed by linear reservoir technique. Rainfall is added directly to the lake inflow. No interception however is abstracted for the lake cells. 9) radar_intensity_map: This is the common name 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.tps 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 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 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 will follow, the links and nodes map would only contain man-made or fluvial (link types 1 and 8) reaches and thus no internal boundary condition link and node will show up in the links and nodes map (however they will be known by the channel_input file, see the attachment). The nodes map must contain the node numbers corresponding to the links map. Thus no node for links representing internal boundary conditions is 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 general rules for stream numbering are: A) Upstream links must have smaller link numbers than downstream links. B) Link numbers may not be skipped. C) Link numbers must start with 1. D) Node numbers increase in the downstream direction. E) Looped reaches cannot be simulated. F) Streams may run only in x-y directions and not diagonally. G) Link types cannot be mixed within a link. For more dicussion of links and nodes maps, see channel_input description and the last section of this document titled "ATTACHMENT for channel_input and table_file File Options". 11) channel_input, table_file: For a detailed description of channel data file set-up and requirements: see the attachment titled "ATTACHMENT for channel_input and table_file File Options", which is a modified version of the document presented at the CASC2D Workshop by F.L. Ogden, 9-10 June, 1994, 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. 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. When -d or -b flag is specified in conjunction with -i flag, these files will be output files at the end of the drainage or backwater run. The drainage or backwater run must precede the actual implicit routing (-i) simulation, for which these two files will be input files and describe the stream network initial conditions. 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. 13) hyd_location: If selected, this input ascii file shall contain the link and node addresses of the locations at which a discharge hydrograph is to be printed. The first column holds link numbers and the second column must be filled with the corresponding node numbers. For the output, see 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 coordinates. If a gage falls outside in a different zone to the left of the active region's UTM 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. 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: 150.0 212.0 45.0 104.0 320.0 173.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 (150.0,212.0), (45.0,104.0), and (320.0,173.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 coordinate system the eastings and northing will have to be in feet rather than meters. 15) outlet_east&north&slope: 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 will be 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. Also 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 chest-board pattern, i.e. oscillations are observed in the water surface level. This may eventually result in crashing the model. As such, the time step needs to be halved 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 thus 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. 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 time span of the 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 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.28. 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 one. 27) unit_space: This conversion factor is used to convert the horizontal grid spacing resolution to meters. The default is one. For state plane coordinate system in ft, set unit_space equal to 3.28. 28) d_thresh: This parameter represents an average step depression 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 often constitutes most of this time. 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 and results of mass balance calculations, are also printed in this file. 30) dis_hyd_location: Whenever hyd_location option is selected, the output will be saved under this option. 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 name 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 dischrge 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. In abnormal exits caused by oscillating depths, negative values may also be observed. 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 specifying the four, for generic event- based infiltration option, or the six, for continuous infiltration option, soil parameter maps. surf_moist_map output contains the soil moisture values only 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. REFERENCES: 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. 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, "Raster-Based Hydrologic Modeling of Spatially-Varied Surface Runoff", Water Resources Bulletin, AWRA, in press for 1995. 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. 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. ATTACHMENT for channel_input and table_file File Options GENERAL NOTES The naming convention associated with the Preissmann 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. Table 1. Link Types Link Description of Link Type Number of Number of Type Parameters Nodes # 1 Fluvial Link, Trapezoidal 5 >=2 Cross-Section 2 Overflow Weir 8 2 3 Culvert/Weir 8 2 4 Reservoirs 5 2 7 Bridge Crossings 8 2 8 Look-up Table (Breakpoint) Cross-Section 4 >=2 At 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 section 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. 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 conditionslink types, only weirs can be simulated at present. A major data entry in the channel_input file is the channel bed elevation, 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 (Ogden et al. 1994) to produce physically realistic longitudinal profiles because of errors inherent in any DEM. EXAMPLE OF INPUT FILE FORMAT AND CONSTRUCTION WARNING: The channel_input file and table_file both follow SI unit convention. In this example, a channel input data file is constructed for a fictitious watershed. The stream network will consist of trapezoidal and look-up table fluvial links, and an internal reservoir. There are no internal boundary conditions in this example other than a reservoir. 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 (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 7. Also links 4 and 5 drain into link 6, which also drains in to the reservoir (link 7). Link 10 originates at the outlet of the reservoir. Links 8 and 9 flow into link 11. Finally, links 10 and 11 merge to form link 12, the most downstream link. Note that the LINK map contains continuous sequence of link numbers, except for the internal boundary conditions which don't show up in the LINK map. Now, refer to the NODE map (nodes_map option). In 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 the 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 4 nodes, for example, it actually has 5. 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. Table 2. Number of Nodes in each Link for Example Watershed Stream Network Link Number Number of Nodes Number of Nodes as as in nodes_map in channel_input file 1 4 5 2 7 8 3 6 7 4 4 5 5 5 6 6 2 3 7 (reservoir) 0 2 8 5 6 9 12 13 10 5 6 11 2 3 12 (most downstream 9 9 link) The first portion of the 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". 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. In our example problem, the number of links is 12, and the maximum number of nodes is 13 (in link 9). Remember that all links which are not at the outlet or not immedietly upstream of a reservoir must have an extra node for connectivity purposes. The largest number of nodes is necessary for memory allocation. 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" 13 This data, which constitutes the first portion of channel_input file, 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 13 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 upstream dependencies (maximum 2, minimum 0) (INT). 4 Upstream dependency #1 (INT). 5 Upstream dependency #2 (INT). 6 Downstream dependency (zero for outlet) (INT). As far as link types in the current example, assume that links 1 through 9, except reservoir link 7, are well-described as trapezoidal (type 1), and we are using look-up table data to describe the cross-sections of links 10, 11, and 12 (type 8). 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 7 4 1 0 0 0 6 5 1 0 0 0 6 6 1 2 4 5 7 7 4 2 3 6 10 8 1 0 0 0 11 9 1 0 0 0 11 10 8 1 7 0 12 11 8 2 8 9 12 12 8 2 10 11 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. As an interpretation, examine link number 10 above. It is of link type 8 (look-up table cross-section data), has 1 upstream dependency (link 7), and flows into link 12. 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 coefficient 0.035 Bottom width (m) 2.50 Channel depth (m) 1.75 Side slope H:V varies Bed elevation dependent upon location Now, lets build the input file section for link #1 which has 4 nodes on the NODE map, plus an extra for connectivity at the upstream end of link 3 (total of 5 nodes). This input file section will look like: 1 5 0.035 2.50 1.75 2.00 265.0 0.035 2.50 1.75 2.00 264.0 0.035 2.50 1.75 2.00 263.0 0.035 2.50 1.75 2.00 262.0 0.035 2.50 1.75 2.00 261.0 The 1 and 5 on the first line indicate that it is link 1, with 5 nodes. The columns represent Mannning's n, bottom width, channel depth, trapezoidal side slope, and talweg elevation, respectively. Note that 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. 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 has a surface area of 0.1 square km, a spillway width of 12m, a spillway discharge coefficient of 0.97, an initial water surface elevation of 249.50m, and a spillway crest elevation of 250m. The resulting input file is shown below. 7 2 0.100 12.0 0.97 249.50 250.00 0.000 0.0 0.00 0.00 0.00 The number of nodes (line entries) per reservoir links is set at two in channel_input file, but the numbers entered in the second row entry are for possible future expansions and thus are currently irrelevant. Also note that the bed elevation at the downstream limit of link 3 must be lower than the assumed initial water surface elevation in the reservoir (link 7). This is mandatory at all downstream boundaries between channels flowing into the reservoirs. Assume that the channels described by look-up table hydraulic properties have essentially 4 different cross-sections. If look-up table data are used, the look-up tables are stored in a different 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 link 10 is approximated by one cross-section, link 11 is approximated by one cross-section, and link 12 is approximated by two cross-sections. So, lets say that link 10 uses table entry 1, link 11 uses table entry 2, and link 12 uses table entry 3 for the its first four nodes and entry 4 for the other five remaining nodes. Link 10 has 6 nodes including the downstream junction node so the channel input file construction for this link is shown below. 10 6 1 244.0 1 243.5 1 243.0 1 242.5 1 242.0 1 241.5 Similarly, link 11 uses look-up table entry #2, and has 3 nodes including the downstream junction node. Again with contrived talweg elevations, the data file looks like: 11 3 2 242.25 2 242.0 2 241.5 Link 12 is the most downstream link and it is described by data from two different look up tables. Link 12 also has 9 nodes. This portion of the file may look like: 12 9 3 241.51 3 241.37 3 241.03 3 240.75 4 240.54 4 240.44 4 240.18 4 239.89 4 239.50 Hence, 239.5 is the talweg elevation at the outlet of the catchment. A complete description of table entries in table_file file will be presented in the next section. If a weir existed in the network, it would not have any entry in links_map or nodes_map whereas connectivity data in the second portion of channel_input file must prove its extistance. The number of nodes in channel_input file for weirs must be set at two and thus two lines of entry should describe the weir parameters as in Table 6: Table 6. Parameter Description for Weirs Line Entry Column Parameter Description 1 1 forward direction weir dischagre coeficient 1 2 reverse direction weir dischagre coeficient 1 3 maximum in-bank height (channel depth), meters 1 4 crest length, meters 1 5 elevation of weir crest, meters 2 1-4 0. (reserved for future expansions) 2 5 bed elevation, meters No weir may be simulated at the very downstream of the network (i.e. outlet); in such cases it is suggested to add a short fluvial 2-node link after the weir and treat this extra link as the most downstream link. FINAL CHANNEL FILE FORMAT The example channel data file (typically called chn.dat) developed for this example looks like: ---------------------------CUT HERE---------------------------- 9.81 1.0 0.5 0.55 100.0 30.0 3600.0 0.07 12 13 1 1 0 0 0 3 2 1 0 0 0 3 3 1 2 1 2 7 4 1 0 0 0 6 5 1 0 0 0 6 6 1 2 4 5 7 7 4 2 3 6 10 8 1 0 0 0 11 9 1 0 0 0 11 10 8 1 7 0 12 11 8 2 8 9 12 12 8 2 10 11 0 1 5 0.035 2.50 1.75 2.00 265.0 0.035 2.50 1.75 2.00 264.0 0.035 2.50 1.75 2.00 263.0 0.035 2.50 1.75 2.00 262.0 0.035 2.50 1.75 2.00 261.0 2 8 0.035 2.50 1.75 2.00 267.0 0.035 2.50 1.75 2.00 266.5 0.035 2.50 1.75 2.00 266.0 0.035 2.50 1.75 2.00 265.0 0.035 2.50 1.75 2.00 264.0 0.035 2.50 1.75 2.00 263.0 0.035 2.50 1.75 2.00 262.0 0.035 2.50 1.75 2.00 261.0 3 7 0.035 2.50 1.75 2.00 261.0 0.035 2.50 1.75 2.00 258.5 0.035 2.50 1.75 2.00 256.0 0.035 2.50 1.75 2.00 254.0 0.035 2.50 1.75 2.00 252.0 0.035 2.50 1.75 2.00 250.0 0.035 2.50 1.75 2.00 248.0 4 5 0.035 2.50 1.65 1.70 255.2 0.035 2.50 1.65 1.70 254.3 0.035 2.50 1.65 1.80 253.2 0.035 2.50 1.65 1.80 252.1 0.035 2.50 1.65 1.90 251.0 5 6 0.035 2.50 1.75 2.10 256.2 0.035 2.50 1.75 2.10 255.0 0.035 2.50 1.75 2.15 254.0 0.035 2.50 1.75 2.20 253.0 0.035 2.50 1.75 2.20 252.0 0.035 2.50 1.75 2.25 251.0 6 3 0.035 2.50 1.75 2.00 251.0 0.035 2.50 1.75 2.00 249.0 0.035 2.50 1.75 2.00 247.4 7 2 0.100 12.0 0.97 249.50 250.00 0.000 0.0 0.00 0.00 0.00 8 6 0.035 2.50 1.75 2.00 245.71 0.035 2.50 1.75 2.00 244.24 0.035 2.50 1.75 2.00 243.55 0.035 2.50 1.75 2.00 243.1 0.035 2.50 1.75 2.00 242.6 0.035 2.50 1.75 2.00 242.25 9 13 0.035 2.50 1.75 2.00 250.1 0.035 2.50 1.75 2.00 249.1 0.035 2.50 1.75 2.00 248.0 0.035 2.50 1.75 2.00 247.0 0.035 2.50 1.75 2.00 246.1 0.035 2.50 1.75 2.00 245.0 0.035 2.50 1.75 2.00 244.0 0.035 2.50 1.75 2.00 243.8 0.035 2.50 1.75 2.00 243.5 0.035 2.50 1.75 2.00 243.2 0.035 2.50 1.75 2.00 242.7 0.035 2.50 1.75 2.00 242.5 0.035 2.50 1.75 2.00 242.25 10 6 1 244.0 1 243.5 1 243.0 1 242.5 1 242.0 1 241.51 11 3 2 242.25 2 242.0 2 241.51 12 9 3 241.51 3 241.37 3 241.03 3 240.75 4 240.54 4 240.44 4 240.18 4 239.89 4 239.50 ---------------------CUT HERE--------------------------- TABLE_FILE INPUT FORMAT Now, 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 4. Then we require 4 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 in each particular table, designated by "numhts". The second number is the vertical distance between hydraulic property points. 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, lets assume we know geometric variables for table entries 1, 2, 3 and 4. The file "table.dat" (table_file option) may look like: 4 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 4 6 1.0 0.0 0.0 0.0 0.0 1.0 4.2 9.4 11.9 2.0 10.1 42.3 112.4 3.0 17.2 101.2 326.5 4.0 23.4 147.7 889.6 5.0 24.5 226.6 1170.3 As a review, consider table 4 above, corresponding to cross sectional values of the five most downstream nodes of link 12. This table has 6 entries, each separated by 1.0 meter of depth. At a depth of 5.0 m, for example, the channel has a top-width of 24.5m, a cross-sectional area of 226.6 m2, and a conveyance of 1170.3 m3/s. REFERENCES 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., August 1-5, 1994, Buffalo, N.Y., 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.