Revision as of 19:55, 23 July 2008

Preface

PREFACE

The work described in this report was authorized by Headquarters, U. S. Army Corps of Engineers (USACE). Funding for this report was provided by the Hydrologic Systems Branch, Coastal and Hydraulics Laboratory (CHL), Engineer Research and Development Center (ERDC) and the System Wide Water Resources Program (SWWRP). At the time of preparation, Mr. Earl Edris was the chief, Hydrologic Systems Branch, CHL, ERDC.

This report was prepared by Dr. Charles W. Downer, Coastal and Hydraulics Laboratory (CHL) Engineer Research and Development Center (ERDC), Dr. Fred L. Ogden, Department of Civil and Environmental Engineering, University of Connecticut, and Mr. Aaron Byrd USACE-ERDC-CHL.

This report was prepared under the general supervision of Mr. Earl Edris, Chief, Hydrologic Systems Branch, CHL, ERDC. Mr. Tom Richardson was Director of CHL. Dr. Steve Ashby (EL) was the SWWRP program manager. The report was reviewed by Dr Mark Jourdan, CH-HW and Dr Jeffery D. Jorgeson, CH-HW.

At the time of publication Dr. Jim Houston was ERDC Director.

This document and the software GSSHA are products of the Watershed Systems Group, Hydrologic Systems Branch, Coastal and Hydraulics Laboratory, Engineer Research Development Center. For more information about GSSHA, contact:

Barbara Parsons

Hydrologic Systems Branch

Coastal and Hydraulics Laboratory

Engineer Research Development Center

3909 Halls Ferry Road

Vicksburg, MS 39180

http://chl.wes.army.mil/software

This report should be cited as follows:

Downer, C. W., Ogden, F. L., and Byrd, A.R. 2008, GSSHAWIKI User’s Manual, Gridded Surface Subsurface Hydrologic Analysis Version 4.0 for WMS 8.1, ERDC Technical Report, Engineer Research and Development Center, Vicksburg, Mississippi.

The Watershed Systems Group (WSG) within the Coastal and Hydraulics Laboratory of the US Army Engineer Research and Development Center (ERDC) supports the US Army and the US Army Corps of Engineers (USACE) in both military and civil operations through the development, modification and application of surface and sub-surface hydrologic models. The Department of Defense (DoD) is also charged with managing approximately 200,000 km² of land within the United States on military installations and flood control and river improvement projects. The WSG provides the Army with predictions of stream flow and stage, inundated areas, saturated areas, soil moistures, groundwater levels, and contaminant fate and transport. Predictions are provided for anticipated changes in weather conditions, project alternatives and land-use changes. The WSG uses a variety of models that are supported by the DoD graphical user interfaces (GUI) Watershed Modeling System (WMS) (Nelson, 2001), Groundwater Modeling System (GMS) (Jones, 2001), and Surfacewater Modeling System (SMS) (Zundel, 2001). These GUIs are commonly referred to the XMS system. The XMS interfaces support a variety of model classes, from simple lumped-parameter runoff models, to 2-D overland, and 3-D unsaturated groundwater models.

For many problems the distributed modeling approach may offer substantial potential improvement in capability compared with traditional lumped-parameter hydrologic models such as the USACE surface hydrologic model HEC-1 (USACE, 1985). The US Army, with additional support from the US Environmental Protection Agency (EPA), funded the development of the physically-based, distributed parameter, Hortonian runoff model CASC2D (Ogden and Julien, 2002; Downer et al., 2002a). Past experience with CASC2D has been favorable when the model has been properly applied, i.e. when Hortonian flow is the dominant process (Doe and Saghafian, 1992; Doe et al. 1996; Ogden et al., 2000; Senarath et al., 2000; Downer et al., 2002a). CASC2D Version 1.18b is linked with WMS Version 5.1 (BYU, 1997a; 1997b), which greatly simplifies model setup, results analysis and visualization. The WSG and the US Army no longer support the development or application of the CASC2D model. CASC2D development continues at Colorado State Univerity.

While Army experience with CASC2D has generally been favorable, there are many instances where the assumptions inherent in the CASC2D model limit its applicability (Senarath et al., 2000; Downer et al., 2002a). Figure 1 illustrates hillslope hydrology with an emphasis on the different runoff and streamflow generating processes. When saturation excess runoff, groundwater discharge to stream, exfiltration, etc., contribute significantly to the stream flow, the application of Hortonian runoff models is ill advised and can lead to erroneous results (Loague and Freeze, 1985; Loague, 1990; Grayson et al., 1992; Smith et al., 1994; Loague and Kyriakidis, 1997; Downer et al., 2002a).

1.1 History

The GSSHA model is a significant reformulation and enhancement of the CASC2D model. The CASC2D runoff model began with a two-dimensional overland flow routing algorithm developed and written in APL (A Programming Language) by Professor P.Y. Julien at Colorado State University. The overland flow routing module was converted from APL to FORTRAN by Dr. Bahram Saghafian, then at Colorado State University, with the addition of Green & Ampt infiltration and explicit diffusive-wave channel routing (Julien and Saghafian, 1991; Julien et al., 1995). The FORTRAN code was reformulated, significantly enhanced, and re-written in the C programming language by Dr. Bahram Saghafian at the U.S. Army Construction Engineering Research Laboratory (CERL). Implicit channel routing was added to CASC2D by Fred L. Ogden (Ogden, 1994), formerly at Colorado State University, then Associate Professor, Department of Civil and Environmental Engineering, University of Connecticut, Storrs, Connecticut, now Cline Distinguished Chair of Engineering, Environment and Natural Resources, Department of Civil & Architectural Engineering and Haub School of Environment and Natural Resources, University of Wyoming. This version, named r.hydro.casc2d, was part of the GRASS GIS for hydrologic simulations (Saghafian, 1993). Work began in 1995 to re-formulate CASC2D with the addition of continuous simulation capabilities, including an interface with the Watershed Modeling System (WMS) interface developed by Brigham Young University (BYU). This version, known as CASC2D for WMS, is distinguished from its predecessors by the addition of a number of new capabilities, numerous improvements and bug fixes, and a more stringent copyright. Johnson et al. (2000) added overland and channel transport to the CASC2D model.

The GSSHA model is a direct result of the dissertation work of Charles W. Downer, USACE-ERDC-CHL (Downer, 2002), and was developed from a need to simulate watersheds with runoff producing processes other than Hortonian flow. While the capability of the CASC2D model was included in GSSHA (many of the processes were taken directly from CASC2D) the continuous nature of the GSSHA model resulting in a need to develop an entirely new model. The first release of the model respresents a fully coupled surface water/groundwater simulator with sediment transport capability (Downer and Ogden, 2006). Since the original development of GSSHA, a myriad of improvements and capabilities have been added to model: coupling of the Green and Ampt with redistribution (GAR) model to the saturated groundwater (Downer et al, 2002), improved channel routing including non-orthogonal stream networks, reservoirs, detention basins and hydraulic structures (Downer et al., 2008), improved soil mositure accounting for use with GAR (Downer, 2008), constituent transport (Downer and Byrd, 2007), and coupling of constituent transport with the Nutrient Simulation Model (NSM) (Johnson and Gerald, 2007). The new features have been tested in a variety of watersheds (Downer et al., 2002, Downer 2008a, Downer 2008b, and others).

1.2 Purpose

GSSHA is intended to be a complete physics based wateshed analysis model and includes important processes related to the generation of runoff, stream routing, overland and stream sediment processes, and constituent transport. The primary purpose for the original version of GSSHA model is to correctly identify and realistically simulate the important hydrologic processes in watersheds. The model is intended to simulate different types of runoff production and determine the controlling physical processes in watersheds, i.e. infiltration excess, saturated source areas, and groundwater discharge. In addition, the model is to physically simulate soil erosion, transport and deposition, as well as constituent transport. Development of the model has been directed by the following requirements:

• model must be capable of explicitly calculating flows, stream depths, and soil moistures in a variety of hydrologic regimes and conditions including non-Hortonian watersheds;
• formulation must account for sub-surface effects on stream flow;
• numerical algorithms must be robust;
• model must conserve mass;
• model simulates soil erosion, transport, and deposition;
• model simulates contaminant transport problems;
• simulation times must be short enough to allow real-time predictions for use at DoD training facilities;
• model must be supported by the standard DoD graphical user interface (GUI) WMS;
• source code must be available to the U.S. Army without restrictions or limitations on modification or publication of results.

1.3 Differences Between GSSHA and CASC2D

Introduction:Differences Between GSSHA and CASC2D

GSSHA is a physically-based, distributed-parameter, structured grid, hydrologic model that simulates the hydrologic response of a watershed subject to given hydrometeorological inputs. The watershed is divided into cells that comprise a uniform finite difference grid. Processes that occur before, during, and after a rainfall event are calculated for each grid cell and then the reponses from individual grid cells are integrated to produce the watershed response. Major components of the model include precipitation distribution, snowfall accumulation and melting, precipitation interception, infiltration, evapo-transpiration, surface water retention, surface runoff routing, channel flow routing, unsaturated zone modeling, saturated groundwater flow, overland sediment erosion, transport and deposition, channel routing of sediments, and constituent fate and transport on the overland and in channels.

During an event, rainfall is spatially and temporally distributed over the watershed. Rainfall may be intercepted by vegetation before reaching the land surface. Once an initial interception demand is reached, a fraction of the precipitation will reach the land surface. Upon reaching the land surface, precipitation may infiltrate due to gravity and capillary forces. Water remaining on the land surface may runoff as two dimensional (2-D) overland flow, after a specfied retention depth representing micro-topography has been reached. This water may eventually enter a stream and be routed to the watershed outlet as one dimensional (1-D) channelized flow. Between precipitation events, soil moisture accounting, evapo-transpiration (ET), and 2-D lateral groundwater flow may be occuring. When precipitation falls in the form of snowfall, the water equivalent volume remains on the land surface and is released as water according to an energy budget calculation.

On the overland flow plane sediment is detached due to rainfall impact and shear stresses due to overland flow. Sediments are routed overland along with the 2-D overland flow. Erosion and deposition continuously occurs on the overland plane as sediments are transported. Sediments may eventually be routed to the stream network where fines (silt and clay) are routed according the advection dispersion equation. Coarse materials are treated as bed load, which is computed according to Yang's method (ref).

Constituents may be assumed to within the soil column or on the land surface. In either case, constituent uptake occurs when water is ponded on the soil surface. Constituents move along in the 2-D overland flow, with reactions occuring as water moves across the watershed. Constituents may ultimately be deposisted into the stream network where they are transported according to the reactive advection dispersion equation.

2.1 Processes Simulated

GSSHA is a process-based model. Hydrologic processes that can be simulated and the methods used to approximate the processes with the GSSHA model are listed in Table 1. With the exception of channel routing, all processes and approximations in the original CASC2D model are also contained in the GSSHA model. The Preissmann channel routing routine (Cunge et al., 1980) was excluded because of known stability problems with the scheme when simulating trans-critical flows (Mesehle and Holly, 1997). Also, the upwind explicit channel routing method was replaced with a similar up-gradient explicit method.

Table 1 – Processes and approximation techniques in the GSSHA model

GA – Green and Ampt (1911), GAR – Green and Ampt with Redistribution (Ogden and Saghafian, 1997), RE – Richards’ equation (1931), ADE – Alternating Direction Explicit (Downer et al., 2000), ADEPC – Alternating Direction Explicit, Predictor-Corrector (Downer et al, 2000).

2.2 Time Steps and Process Updates

In GSSHA the user specifies the overall model time step, in seconds, that the model uses to loop through the processes, check update times, and update processes. To avoid missing updates of processes, such as rainfall, that may be specified at 1 minute intervals, the overall model time step should be integer divisible into 60 seconds or an integer multiple of 60 s (i.e. 5,10,15,20,30,60,120,180,300). Time steps such as 7,9,13,16,21,45,90,270 should not be used, as they may result in unexpected internal model behavior. The model time step also must not be greater than the finest resolution of inputs, such as rainfall. Typical time steps for GSSHA range from 10 to 300 seconds. Smaller time steps may be required for particularly difficult problems.

The computational time step is an important parameter affecting the performance of GSSHA. In addition to setting the pace of the model, the overall model time step is used to set or initialize the temporal discretizations of many model processes. While many processes, such as channel routing, saturated and unsaturated groundwater flow, have internal model stability checks, some methods of overland flow routing do not. If the time step is too large the program may crash or produce inaccurate results. Very small time steps result in inordinately long simulation times. The best way to determine the most efficient time step is through a temporal convergence study, where the time step is varied and the model behavior is observed. This allows the user to determine the maximum time step that can be used with acceptable accuracy. As the time step is increased, the outlet hydrograph will begin to shift in position in relation to simulations with smaller time steps. As the time step is further increased, the discharge at the outlet will oscillate; further increases in the computational time-step will result in program crashes. Results of a temporal convergence study, featuring hydrograph shifting, oscillations, and model crash, are shown in Figure 4.

The appropriate time step strongly depends on watershed and rainfall characteristics. In general, 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.

Shorter time steps must be used when backwater effects are generated in flat areas in the digital elevation model (DEM). If the time step is too long for any particular simulation the surface water depth in very flat areas may develop a checkerboard pattern due to oscillations in the water surface level. This eventually results in a crash. If this occurs the time step should be decreased and the simulation repeated.

At the time of this publication, only the time step for the saturated groundwater flow is specified in addition to the overall time step. The ET time step is fixed at 1 hr, the usual interval of hydrometeorological data available. To maintain stability the time step may be reduced internally for the explicit channel routing code, the unsaturated zone RE solver, the groundwater solver, and the explicit and ADE solutions for overland flow. Internal time step limitations in the model are described under the appropriate process sections. Rainfall updates are specified in the rainfall gage file and the interval between updates can vary as needed. Thus the overall time step is limited by:

1. stability issues in the overland flow scheme,
2. the smallest rainfall interval,
3. the groundwater time step, and
4. the need to be integer divisible into the groundwater time step and the smallest rainfall interval.

Timesteps for sediment and constituent fate and transport are based on the underlying hydrologic proecesses and do not have to be specified.

Guidance for time steps is shown in Table 2.

Table 2 – Recommended time-steps and stability criteria used in the GSSHA model

2.3 Inputs

GSSHA is a distributed-parameter, process-based model that requires the user to select the processes to be simulated and then provide the model with the data necessary to drive the selected options. Three types of input data are used. An ASCII text project file is used to provide the basic project information, select processes to be simulated, assign simulation parameters, and locate data files, tables and maps. Spatially distributed parameters can be assigned with maps of ASCII gridded data with a parameter value in each grid cell, with index maps and tables of parameter values that relate to the index maps, or with uniform values in every cell. Typically the data required to assign parameter values in every cell is not available. Standard practice in the application of GSSHA has been to develop index maps based on available data sources of land use, soil type, and vegetation. Typically these maps are combined to create a master land use/soil type/vegetation index map that can be used to assign all parameter values.

Parameters for each index map are then assigned using tables that reference the values in the index maps. A detailed description of using the index maps and Mapping Tables to build a GSSHA model is provided in Section 11. If available, the detailed maps containing parameter values in each cell, as described in Ogden (2000), may be input in lieu of the index maps and table.

Since distributed parameters may be assigned with a single uniform value in the project file, a table value linked to an index map, or with an ASCII map with a parameter value for every grid cell the model has been developed to prioritize parameter specification. While internally assigning parameter values the GSSHA model looks for the most detailed information first, GRASS ASCII maps, the second most detailed information second, table values linked to index maps, and finally a single uniform value from the project file. Once data from one of these sources is located, the search ends and the parameter values are assigned inside the GSSHA model. While this rule is generally applicable it is prudent not to specify multiple sources of the same parameter value in the project file. This avoids possible confusion and improper assignment of parameter values.

The Watershed Modeling System (WMS) interface, developed at the Environmental Modeling Research Laboratory (EMRL) at Brigham Young University, is recommended for developing input files and viewing output from the GSSHA model. The WMS produces GSSHA specific files from general Geographic Information System (GIS) data. WMS does not replace the functions of a GIS, though it can accept information in a variety of GIS formats. GSSHA relies on the GRASS ASCII data file format for storing spatially distributed data. The GRASS GIS is very helpful in the preparation of GSSHA data sets. Users of ARC/INFO and ARCVIEW can export data to GSSHA through the WMS interface. For more information about WMS, DoD and EPA personnel should contact:

XMS Model Support

Hydrologic Systems Branch

Coastal and Hydraulics Laboratory

Engineer Research Development Center

3909 Halls Ferry Road

Vicksburg, MS, 39180

(601) 634-4286

http://chl.wes.army.mil/software

Other users seeking information about WMS, should contact:

Aquaveo

3210 N Canyon Road

Provo, Utah 84604

801-691-5530

http://www.aquaveo.com/technical-support

email: support@aquaveo.com

GSSHA simulations require a project file that contains command line instructions or “cards” which pass options to GSSHA for a particular simulation. The name of the project file is given at run time as a command line argument for GSSHA. The following section presents all project file cards followed by a brief description of each. The project file consists of a single card on each line, followed by its argument, if any. While some cards require no argument, others require values, character strings, file names, table names, or map names. Tables are files that contain ASCII input data in a tabular format. A map name is simply the name of a floating point GRASS ASCII file that contains raster data. An index map refers to a similar file that contains integers indexed to tabular values. Throughout the manual, project card names will be presented in BOLD CAPITAL letters; arguments appear as CAPITAL ITALICS.

The project file cards may appear in the project file in any order, except the first card must be GSSHAPROJECT. Extraneous or misspelled cards are ignored. A card may be commented out by preceding the card with a pound sign (hash mark) “#”. When using the Mapping Table to assign parameter values for any process, the MAPPING_TABLE card is used. For project cards related to input parameters the units, if any, of the input argument are presented. For optional inputs, the default value, if any, is also presented.

3.1 Required Inputs

The following table lists all required project file card, which must be present in a project file. In addition to these cards, a method of assigning rainfall and overland roughness values must also be selected from the choices below.

GSSHAPROJECT

GRIDSIZE ##.##

ROWS  ##

COLS  ##

TOT_TIME  ##

TIMESTEP  ##

OUTSLOPE  ##.##

OUTROW  ##

OUTCOL  ##

ELEVATION  ##.##

WATERSHED_MASK "filename.msk"

HYD_FREQ  ##

SUMMARY  "filename.sum"

OUTLET_HYDRO  "filename.otl"

3.2 Mapping Table – Optional

When using the Mapping Table to assign any of the distributed parameters the Mapping Table card must be included.

MAPPING_TABLE  "filename.cmt"

ST_MAPPING_TABLE  "filename.smt"

3.3 Overland Flow – Required

3.3.1 Required Inputs

Overland roughness coefficients must be assigned to every cell in the active grid. Overland roughness coefficients may be assigned as either a single uniform value, from a GRASS ASCII grid map, or from the Mapping Table. Only one option should be selected. Do not specify either of these cards if roughness values are to be read from the Mapping Table file.

MANNING_N  ##.##

ROUGHNESS  "filename.ext"

3.3.2 Optional Inputs

The user can select the overland flow routing scheme. An overland retention depth may also be specified at the discretion of the user. Retention depth can be assigned as a single uniform value, from the Mapping Table, or from a GRASS ASCII grid map. Only one option should be specified.

OVERTYPE   [type]

RETEN_DEPTH  ["filename.ext"]

RETENTION ##.##

READ_OV_HOTSTART  "filename.ext"

WRITE_OV_HOTSTART  "filename.ext"

OV_BOUNDARY

INIT_ELEV_HEAD ##.##

OVERLAND_MOMENTUM 

OVERLAND_STRICT_DT 

To include wetlands, which are optional, simply set up a wetlands table. The wetlands are turned on if the table is present and has a valid index map associated with it.

3.4 Interception – Optional

Interceptions parameters may be assigned using the Mapping Table, or by assigning GRASS ASCII maps.

INTERCEPTION

STORAGE_CAPACITY  "filename.ext"

INTERCEPTION_COEFF  "filename.ext"

3.5 Rainfall Input and Options – Required

The following table lists all project file cards pertaining to rainfall input. Either, but not both, PRECIP_UNIF or PRECIP_FILE must be included in each GSSHA project file.

PRECIP_UNIF

RAIN_INTENSITY  ##.##

RAIN_DURATION  ##

START_DATE [yr mo day]

START_TIME [hr min]

PRECIP_FILE  "filename.gag"

RAIN_INV_DISTANCE

RAIN_THIESSEN

3.6 Infiltration – Optional

Infiltration may be calculated using four different infiltration options. Green and Ampt (GA) (Green and Ampt, 1911), a multi-layered Green and Ampt Model, Green and Ampt with Redistribution (GAR) (Ogden and Saghafian, 1995), and Richards’ equation (RE) (Richards, 1931). Only one of these four methods should be selected.

GREEN_AMPT

INF_REDIST

INF_LAYERED_SOIL

INF_RICHARDS

3.6.1 Green and Ampt (GA)

When GREEN_AMPT is selected, values of hydraulic conductivity, wetting front suction head, porosity, and initial moisture are required. Parameter values may be input using the Mapping Table file or with the series of GRASS ASCII maps using the cards below.

CONDUCTIVITY  "filename.ext"

CAPILLARY  "filename.ext"

POROSITY  "filename.ext"

MOISTURE  [value]

3.6.2 Green and Ampt with Redistribution (GAR)

When GAR is specified, the GA parameters plus two additional parameters must be provided: the pore-distribution index, and residual saturation. These may be input in the Mapping Table file or by providing two GRASS ASCII maps using the cards below.

PORE_INDEX  "filename.ext"

RESIDUAL_SAT  "filename.ext"

FIELD_CAPACITY  "filename.ext"

3.6.3 Multi-layered Green and Ampt

The multi-layered GA parameters are assigned with the table described in Section 11, and referenced to an index map. The required project cards are listed below.

SOIL_TYPE_MAP "filename.ext"

SOIL_LAYER_INPUT_FILE  "filename.ext"

3.6.4 Richards’ Equation

When INF_RICHARDS is selected additional cards are required. There are also a number of optional input and output cards, as described below. Richard' Equation parameters may be specified in the MAPPING_TABLE, as well as in separate ASCII files.

3.6.4.1 Required Inputs

INF_RICHARDS

RICHARDS_C_OPTION  [value]

3.6.4.2 Parameter Assignment - Required

Select either SOIL_TYPE_MAP and SOIL_LAYER_INPUT_FILE or use the Mapping Table.

SOIL_LAYER_INPUT_FILE "filename.ext"

SOIL_TYPE_MAP  "filename.ext"

3.6.4.3 Optional Inputs

WATER_TABLE  "filename.ext"

AQUIFER_DELTA_Z  ##.##

SEASONAL_RS

RICHARDS_UPPER_OPTION  [value]

GW_ASSIGN_THETA

RICHARDS_ITER_MAX  ##

RICHARDS_WEIGHT  ##.##

RICHARDS_K_OPTION  [value]

RICHARDS_DTHETA_MAX  ##.##

3.6.4.4 Optional Output

IN_THETA_LOCATION  "filename.ext"

3
1 1
10 9
100 90  

OUT_THETA_LOCATION  "filename.ext"

3.6.7 Special Infiltration Cards

READ_SM_HOTSTART  "filename.ext"

WRITE_SM_HOTSTART  "filename.ext"

3.7 Channel Routing – Optional

3.7.1 Required Inputs

The following cards are REQUIRED for the simulation of channel routing.

DIFFUSIVE_WAVE

CHAN_EXPLIC

CHANNEL_INPUT  "filename.cif"

STREAM_CELL  "filename.gst"

OUTLET_HYDRO  "filename.otl"

SECTION_TABLE  "filename.ext"

3.7.2 Initial Condition and Boundary Condition - Optional

EXPLIC_BACKWATER  "filename.ext"

EXPLIC_HOTSTART  "filename.ext"

WRITE_CHAN_HOTSTART  "filename"

READ_CHAN_HOTSTART  "filename"

OVERBANK_FLOW

OVERBANK_MAX_DV

OVERLAND_BACKWATER

CHAN_POINT_INPUT

MAX_COURANT_NUMBER

HEAD_BOUND

BOUND_DEPTH

BOUND_TS

CHANNEL_MOMENTUM

3.7.3 Stream Losses/Gains - Optional

3.7.3.1 General

STREAM_LOSS

3.7.3.2 Parameters

Used with WATER_TABLE or STREAM_LOSS card. These parameters may also be distributed along the channel by assigning in the CHAN_INPUT file.

K_RIVER  ##.##

M_RIVER  ##.##

3.7.4 Optional Output

IN_HYD_LOCATION  "filename.ihl"

OUT_HYD_LOCATION  "filename.ohl"

OUT_DEP_LOCATION  "filename.odl"

IN_SED_LOC  "filename.isl"

OUT_SED_LOC  "filename.osl"

IN_GWFLUX_LOCATION  "filename.igf"

OUT_HYD_LOCATION  "filename.ogf"

STRICT_JULIAN_DATE

CHAN_DEPTH  "filename.cdp"

CHAN_STAGE  "filename.cds"

CHAN_DISCHARGE  "filename.cdq"

CHAN_VELOCITY  "filename.cdv"

LAKE_OUTPUT    "filename.lel"

3.8 Continuous Simulations – Optional

Continuous simulations require general information about the watershed location, selection of a method to calculate evapo-transpiration (ET), hydrometeorological (HMET) data in one of three available formats, and the appropriate distributed data either from the Mapping Table file or from GRASS ASCII maps.

3.8.1 Required Inputs

LONG_TERM

LATITUDE  ##.##

LONGITUDE  ##.##

GMT  ##.##

SOIL_MOIST_DEPTH  ##.##

EVENT_MIN_Q  ##.##

ET_CALC_PENMAN

ET_CALC_DEARDORFF

3.8.2 Seasonal Canopy Resistance - Optional

SEASONAL_RS

SEASONAL_RS_SPRING

SEASONAL_RS_SUMMER_START

SEASONAL_RS_SUMMER_END

SEASONAL_RS_FALL

3.8.3 Format of Hydrometeorological (HMET) Data – Required, Select One Format

HMET_SURFAWAYS  "filename.hmt"

HMET_SAMSON  "filename.hmt"

HMET_WES  "filename.hmt"

3.8.4 ET Parameter Assignment – Required, Select Mapping Table or GRASS ASCII maps

Long-term simulation parameters must be assigned using either the Mapping Table or providing the GRASS ASCII maps as described below. Albedo, wilting point, transmission coefficient, vegetation height and canopy resistance are also required for ET_CALC_PENMAN.

ALBEDO  "filename.alb"

WILTING_POINT  "filename.wtp"

TCOEFF  "filename.tcf"

VHEIGHT  "filename.vht"

CANOPY  "filename.cpy"

3.8.5 Optional Inputs

TOP_LAYER_DEPTH ##.##

END_TIME [yr mo day hr min]

START_DATE [yr mo day ]

START_TIME [hr min]

3.8.6 Snow Card Inputs - Optional

Prior to version 6.1 there are no snow options. Additional snow capability has been added in v6.1 and beyond. Please note the GSSHA version numbers when using these cards.

Cards calling which snow melt algorithm to use

default (no card required)

NWSRFS_SNOW

EB_SNOW

Cards Associated with All Three Melt Methods

NWSRFS_SCF     ##.##

SNOW_TEMP_BASE ##.##

SNOW_NO_INFILTRATE

INIT_SWE_DEPTH #.# or File

SNOW_SWE_FILE  ***.swe

Example file when using INIT_SWE_DEPTH

Cards Associated with BOTH Hybrid Energy Balance and Temperature Index Methods

NWSRFS_FR_USE  ##.##

NWSRFS_TIPM       ##.##

NWSRFS_NMF        ##.##

NWSRFS_FUA        ##.##

NWSRFS_PLWHC      ##.##

NWSRFS_ELEV_SNOW  File

Example file when using NWSRFS_ELEV_SNOW

Elevations are in meters, all other values are in their standard formats.

Cards Associated with JUST Temperature Index Method

NWSRFS_MF_MAX  ##.##

NWSRFS_MF_MIN  ##.##

Cards Associated with Vertical Melt Water Transport (Vertical MWT)

SNAP_RETENTION

VERT_SNOW_RETENTION

Cards Associated with Lateral Melt Water Transport (Lateral MWT)

ROUTE_LAT_SNOW     ##.##

SNOW_DARCY     ##.##

Cards Associated with Orographic Effects Orographic effects are available in v6.1 and beyond. Note that the cards change between versions 6.1 and v6.2

HMET_OROG_GAGES  ***.txt

HMET_ELEV_GAGE   ##.##

OROGVAR_HMET

HMET_LAPSE_RATE  ##.##

YES_DALR_FLAG  ##.##

3.8.7 Distributed Hydrometeorology Data - Optional

This function allows GSSHA to read in raster-based hydrometeorology data. The input files must be in the same projection as the GSSHA model and must be larger than the model domain. Please see Distributed HMET Data for more details.

HMET_ASCII  ***.txt

3.8.8 Continuous Frozen Ground Index (CFGI) Index Model - Optional

These options are only valid with verions 6.2 and later. In previous versions the CFGI model is applied for any LONG_TERM simulations.

CFGI 

CFGI_INDEX  

CFGI_K  

GTFSM

3.9 Saturated Groundwater Flow – Optional

3.9.1 Required Inputs

These cards are REQUIRED to perform 2-D lateral groundwater simulations.

WATER_TABLE  "filename.wte"

GW_SIMULATION

GW_TIME_STEP  ##.##

GW_LSOR_CON  ##.##

AQUIFER_DELTA_Z  ##.##

GW_RELAX_COEFF  ##.##

AQUIFER_BOTTOM  "filename.bot"

GW_BOUNDFILE  "filename.bnd"

3.9.2 Optional Inputs

GW_ASSIGN_THETA

GW_UNIF_POROSITY

GW_POROSITY_MAP  "filename.por"

GW_UNIF_HYCOND  ##.##

GW_HYCOND_MAP  "filename.hyd"

K_RIVER  ##.##

M_RIVER  ##.##

GW_FLUXBOUNDTABLE  "filename.flx"

GW_LEAKAGE_RATE  ##.##

SINGLE_UNSAT_SAT  ##.##

• Values for porosity (θ_s) and saturated hydraulic conductivity (K_s) may be specified in one of four (4) ways:

1. Without specifying one of the two types of cards above, θ_s and/or K_s values for each cell will be assigned the value of bottom layer of the SOIL_LAYER_INPUT_FILE for the soil type as specified in the SOIL_TYPE_MAP.
2. Uniform values of θ_s or K_s may be specified using the appropriate card above.
3. GRASS ASCII maps of θ_s or K_s values for each cell may be assigned by using the appropriate card above.
4. Values may be assigned in the MAPPING_TABLE

3.9.3 Optional Output

GW_OUTPUT  "filename.gw"

OUT_WELL_LOCATION  "filename.iwl"

GW_WELL_LEVEL  "filename.owl"

IN_GWFLUX_LOCATION  "filename.igf"

OUT_GWFLUX_LOCATION  "filename.ogf"

3.10 Soil Erosion – Optional

3.10.1 Soil Erosion Simulation - Required

SOIL_EROSION  ##

OUTLET_SED_FLUX  "filename.osf"

3.10.2 Soil Erosion Parameters - Required

Soil erosion parameters must be specified in the Mapping Table. See description of inputs in Mapping Table Section - Chapter 12

3.10.3 Optional Inputs

SAND_SIZE  ##.##

WATER_TEMP  ##.##

SED_POROSITY  ##.##

ADJUST_ELEV  "filename.ele"

ADJUST_CHANNEL_BED  "none"

IN_SED_LOC  "filename.isl"

OUT_SED_LOC  "filename.osl"

OUTLET_SED_TSS "filename.oss"

OUT_TSS_LOC "filename.tss"

3.10.4 Optional Outputs

NET_SED_VOLUME "base_filename.ext"

VOL_SED_SUSP "base_filename.ext"

3.11 Constituent Transport – Optional

3.11.1 Contaminant Transport - Required

OV_CON_TRANS

3.11.2 Contaminant Transport Parameters – Required

Contaminant transport parameters must be specified in the Mapping Table. See description of inputs in Mapping Table Section - Chapter 12

3.11.3 Optional Inputs

SOIL_CONTAM

MIXING_LAYER_DEPTH  ##.## 

SOIL_STATIC_CONC 

SOIL_NOFLUX 

CONTAM_MAP  filename

CHAN_CON_TRANS

CHAN_DECAY_COEF  ##.##

CHAN_DISP_COEF  ##.##

OUT_CON_LOCATION  "filename.ocl"

OUT_MASS_LOCATION  "filename.oml"

ALL_CONTAM_OUTPUT

3.12 Subsurface Drainage Network – Optional

Subsurface drains can be simulated with either SUPERLINKS and/or a simple link/node routing package. The inputs are similar

3.12.1 Specifying the Methods

The following card is REQUIRED for the simulation of SUPERLINK subsurface drains.

STORM_SEWER

The following card is REQUIRED for the simulation of non-SUPERLINK subsurface drains.

STORM_DRAIN

3.12.2 Optional Inputs

The following cards are optional for the simulation of SUPERLINK subsurface drains.

GRID_PIPE

HIGH_HEAD_RELEASE

SUPERLINK_C_OPT

SUPERLINK_DRAINMOD

The following cards are optional for the simulation of non-SUPERLINK subsurface drains.

DRAIN_GRID_PIPE

SUPERLINK_DRAINMOD

3.12.3 Output Cards

The following cards are optional for the simulation of SURERLINK subsurface drains.

SUPERLINK_JUNC_LOCATION

SUPERLINK_NODE_LOCATION

SUPERLINK_JUNC_FLOW

SUPERLINK_NODE_FLOW

The following cards are optional for the simulation of non-SURERLINK subsurface drains.

DRAIN_OUTPUT

3.13 Output Files – Required

3.13.1 Required Output

SUMMARY  "filename.sum"

OUTLET_HYDRO  "filename.otl"

3.13.2 Optional Output

QOUT_CFS

QUIET

SUPER_QUIET

STRICT_JULIAN_DATE

OPTIMIZE  "filename.opt"

OPTIMIZE_SED  

PROJECT_PATH  

REPLACE_PARAMS  "filename.in"

REPLACE_VALS  "filename.in"

IN_HYD_LOCATION  "filename.ihl"

OUT_HYD_LOCATION  "filename.ohl"

OUT_DEP_LOCATION  "filename.odl"

OUTLET_SED_FLUX  "filename.osl"

 Requires SOIL_EROSION.

OUTLET_SED_TSS  "filename.oss"

 Requires SOIL_EROSION.  

IN_SED_LOC  "filename.isl"

OUT_SED_LOC  "filename.osl"

OUT_TSS_LOC  "filename.tss"

OVERLAND_DEPTH_LOCATION  "filename.odi"

OVERLAND_DEPTHS  "filename.odo"

OVERLAND_WSE_LOCATION  "filename.owi"

OVERLAND_WSE  "filename.owo"

OVERLAND_Q_CUM_LOCATION  "filename.cdl"

CUM_DISCHARGE  "filename.cds"

IN_THETA_LOCATION  "filname.itl"

OUT_THETA_LOCATION  "filename.otl"

OUT_WELL_LOCATION  "filename.igw"

GW_WELL_LEVEL  "filename.ogw"

IN_GWFLUX_LOCATION  "filename.igf"

OUT_GWFLUX_LOCATION  "filename.ogf"

OUT_CON_LOCATION  "filename.ocl"

OUT_MASS_LOCATION  "filename.oml"

3.13.3 Optional Output Maps

MAP_FREQ  ##

MAP_TYPE  ##

DISCHARGE  [value]

CUM_DISCHARGE_MAP  [value]

OVERLAND_CHAN_DISCHARGE  [value]

DEPTH  [value]

INF_DEPTH  [value]

SURF_MOIST  [value]

RATE_OF_INFIL  [value]

DIS_RAIN  [value]

CHAN_DEPTH  [value]

CHAN_DISCHARGE  [value]

MAX_SED_FLUX  [value]

NET_SED_VOLUME  [value]

GW_OUTPUT  [value]

GW_RECHARGE_CUM  [value]

GW_RECHARGE_INC  [value]

** if running on a linux/unix machine, then the path uses a forward slash instead of a backward slash.

4.1 Units

With few exceptions the units used for input and output in GSSHA are in SI units. Units for particular processes are consistent with those found in the literature. For instance, hydraulic conductivity is specified in units of cm/h, while rainfall rates are specified in mm/h. The particular units expected of each input are given in the project file card descriptions. The discharge hydrograph may be output in cfs by specifying the QOUT_CFS flag.

Furthermore, the model was developed based on the Universal Tranverse Mecator (UTM) coordinate system, as described in Section 4.4.

4.2 Grid Size

Distributed models are used over a very wide range of grid sizes, from 10 to 1000 m. The selection of an appropriate grid size for a GSSHA model requires consideration of both the available data and computational effort required. Typical grid sizes range from 10 to 250 m. The selection of the grid size for a given watershed determines the total number of grid cells used to describe the watershed, setting the computational effort and memory required. Note that if the grid size is halved, the memory required and computational time increase by a factor of 4. In general, smaller grids are less sensitive to sub-grid variability for Hortonian runoff (Ogden and Julien, 1993). Therefore, smaller grids are generally “better”, if the data exists to assign relevant watershed characteristics to each grid cell. However, smaller grid-sizes do not guarantee superior model performance. As with the time step, the most appropriate grid size can be determined with a convergence study, where the effects of increasing or decreasing the grid size can be observed on model output.

4.3 Total Event Simulation Time

This project file card is used to specify the total GSSHA event simulation time in minutes. If the volume of water remaining on the surface at the end of the simulation is greater than 5% of the rainfall volume, the total simulation time is too short to capture the entire runoff hydrograph, and a warning is printed in the run summary (SUMMARY card) file. The TOT_TIME card is ignored for continuous simulations. For continuous simulations the end of an event is the time at which the outflow discharge falls below the discharge specified by the EVENT_MIN_Q project file card. For continuous events the EVENT_MIN_Q is used only for accounting purposes, except for the last event when it is also used to stop the simulation. For some simulations the EVENT_MIN_Q value may never be reached. In that case the user should specify the END_TIME card. The END_TIME card can be used to stop any long term simulation at any desired point. The simulation may end before the specified END_TIME if conditions dictating the end of the simulation are encountered before the specified END_TIME.

4.4 Coordinate System

GSSHA performs calculations on raster grids and data can be input with GRASS ASCII raster data. While any logical, and consistent, coordinate system can be used, the preferred coordinate system for raster based data is the UTM map projection. The UTM System breaks the entire earth into zones 6 degrees of longitude wide. The zones that cover the continental United States are:

The standard specification of the UTM system includes (Davis et al., 1981):

• The reference ellipsoid is Clarke 1866 in North America.
• The origin of longitude is the central meridian.
• The origin of latitude is the equator.
• The unit of measure is the meter.
• A false easting of 500,000 m is used for the central meridian of each zone.
• The scale factor at the central meridian is 0.9996.
• The zones are numbered beginning with 1 for the zone between 180° W and 174° W meridians and increasing to 60 for the zone between meridians 174° E and 180° E.
• The latitude for the system varies from 80° N to 80° S.
• In the southern hemisphere, a false northing of 10,000,000 m is used.
• The scale error is 1/2500 on the central meridian.

The UTM coordinate system was chosen because of its global applicability and widespread acceptance. Note that watershed data that lie in two zones must be merged into one zone. The data should be merged into the zone that contains the majority of the watershed area. This is accomplished using a GIS.

4.5 Map Headers

Many spatially varied inputs for GSSHA can be input as ASCII GRASS GIS raster maps. Raster maps are defined as maps that assign attributes to areas, as opposed to vector maps, which assign attributes to lines or polygons. GSSHA requires that all raster grid cells are square.

Each raster map, including index maps, must have a header that conforms to the following example:

The entries for north, south, east, and west are the coordinates of the bounding rectangle that contains the entire watershed. The rows and cols entry in the header respectively contain the number of rows and columns in the watershed. Note that for this particular example header, that (north-south)/rows = 200 m, and (east-west)/cols=200 m. Therefore, the grid size of this particular model is 200 m, and the raster cells are square, as required by GSSHA. As the coordinate values in the map header are not used to determine geographic position on the globe, any values for north, south, east, and west may be used in the map headers provided that the grid size is accurate, the grids are square, and your maps are consistent. If the grid size calculated from the bounding rectangle, rows, and columns is not equal to the value specified in the project file using the GRID_SIZE project file card, GSSHA will not run. The header in each GRASS ASCII map must be IDENTICAL. This requirement forces the user to ensure that all maps are of the same geographic region and grid size.

The header is followed by rows and columns of space delimited data values. These values can be either integer or real, depending on the data type. For example, the mask, channel link, channel node, and index maps are by definition integer maps; maps containing the land surface and bedrock elevations contain real values.

4.6 Watershed Mask

The required WATERSHED_MASK project file card is used to input the name of the file containing the mask map. The mask map is used to define the watershed boundaries within the rectangular grid and reduces memory requirements by an amount directly proportional to the ratio of mask area over the bounding rectangle area. The watershed mask is a map containing only 0s and 1s, where 0 and 1 represent grid cells that are outside and within the watershed, respectively. For instance, a watershed that is approximately the shape of Texas would have the following a mask map file:

north: yyyyyyy.y1

south: yyyyyyy.y2

east:  xxxxxx.x1

west:  xxxxxx.x2

rows:  15

cols:  12

0 0 0 0 1 1 1 0 0 0 0 0

0 0 0 0 1 1 1 0 0 0 0 0

0 0 0 0 1 1 1 0 0 0 0 0

0 0 0 0 1 1 1 0 0 0 0 0

0 0 0 0 1 1 1 1 1 1 0 0

0 0 0 0 1 1 1 1 1 1 1 1

0 0 0 0 1 1 1 1 1 1 1 1

0 0 0 0 1 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1 1

0 1 1 1 1 1 1 1 1 1 1 1

0 0 0 1 1 1 1 1 1 1 1 0

0 0 0 0 1 0 1 1 1 1 0 0

0 0 0 0 0 0 1 1 1 0 0 0

0 0 0 0 0 0 1 1 1 0 0 0

0 0 0 0 0 0 0 1 0 0 0 0

If the watershed mask contains a 0 in a particular grid cell, GSSHA ignores data in all input maps for that particular grid cell. In the above map, there are a total of 15*12=180 grid cells. There are 85 grid cells with a 1 in the watershed mask. In this example, GSSHA would only allocate memory for 85 of the 180 grid cells for each map, representing a decrease of 53% of memory required to store the data for each map. Errors in watershed delineation (hence mask creation) will propagate through all data sets input to and output from GSSHA. The watershed should be delineated with care. WMS creates the watershed mask as part of the watershed delineation processes. GRASS users can use the r.watershed command to create the mask. During each simulation, the GSSHA model writes out the watershed mask with the internal number assigned to each grid cell in GRASS ASCII format in a file called “maskmap”. This file is useful for deciphering error messages from GSSHA that refer to the grid cell number where the problem occurred.

4.7 Elevation Map

Elevation data are input for every active grid cell in a GRASS ASCII map file specified with the ELEVATION project file card. The elevation data are perhaps the most important inputs for GSSHA modeling. The quality of the elevation data plays a major role in success of GSSHA simulations. Elevations in the grid are derived from digital elevation model data (DEM).

DEMs always contain errors. Large flat areas in the DEM may be due to the limited vertical resolution of elevation data from which the DEM was derived. Extensive flat areas usually cause problems for the 2-D explicit diffusive-wave overland flow routing used in GSSHA. Digital dams, pits, and depressions in the DEM may be artifacts of the interpolation scheme used to rasterize digitized contours, or due to coarse resolution in regions of concave topography.

In addition, the grid size used in a GSSHA model is normally coarser than the available DEM data. Elevation data in the grid must be somehow interpolated from the DEM data and lumped into larger areas. This process can introduce additional error in the elevation of grid cells. As a rule, the user must cross check the elevation values with in-field observations or topographic maps of the area. Digital topographic maps are often available and can be displayed as a background image in WMS.

One way to discover potential errors in the elevation data is to perform a simulation with the most basic GSSHA model: a single event with uniform rainfall, overland flow, a relatively short time step, and no other options. Surface depth output maps should be written frequently (see DEPTH and MAP_FREQ project file cards). If the simulation finishes without an error, the surface depth maps should be examined to determine where most water accumulates and whether such accumulations are justified by the topographic map of the watershed. Alternatively, the model may crash. The location of grid cells where problems occurred will be printed on the screen and also at the bottom of the run summary file.

Editing of the elevations in the grid is often necessary to impose the actual drainage trend observed in the topographic map. Digital dams, pits and depressions must be removed since they trap surface runoff that would otherwise contribute to the outlet discharge. Using grids with raw elevations requires shorter computational time steps, while properly prepared grids, particularly those with coarser resolution, allow use of longer time steps. If channel routing is performed, care should also be taken to ensure that overland flow runoff reaches grid cells that contain channel links. If the stream network is delineated independently from the DEM, e.g. from a digital line graph (DLG), then the elevations of the grid cells containing channel nodes should be checked to insure they are not higher than those of the surrounding cells. Otherwise the overland flow will not be correctly passed to the channels.

The cleandam.exe routine can be used to do much of the above. Cleandam.exe is described in a another section of the manual.

WMS offers the option of importing and displaying vector stream location files (DLG) to aid in stream channel delineation. WMS users can also automatically delineate streams from the DEM using the tools in the WMS software. GRASS users may use the r.watershed command to automatically delineate the streams. As either automatic delineation relies on the DEM to determine stream locations, there may be substantial differences from the DLGs. If an automated method is used to locate the stream locations from the DEM data and the grid resolution used in the GSSHA model is coarser than the available DEM data the stream may not fall in the lowest elevation grid cells. However the stream is located, the elevations of the cells containing the stream may need to be manually adjusted to ensure proper overland-channel interaction.

4.8 Optimizations

Senarath et al. (2000) demonstrated that the CASC2D model could be effectively parameratized by use of an automated calibration procedure, such as the shuffled complex evolution (SCE) method (Duan et al., 1992). Output that allows optimizations based on event peak discharge and event discharge volume can be produced by using the OPTIMIZE card. This card specifies a file that contains the peak discharge and discharge volume for each event in the rainfall file. Peak discharge (cms) and volume (m³) are written for the outlet and also at any locations in the stream network specified in the IN_HYD_LOCATION file. The last line contains the total outlet discharge volume for the entire simulation. The output format is:

GSSHA uses similar two-step explicit finite volume schemes to route water for both 1-D channels and 2-D overland flow, where flows are computed based on heads, and volumes are updated based on the computed flows. Compared with more sophisticated implicit finite difference and finite element schemes, the algorithms used in GSSHA are simple. The friction slope between one grid cell and its neighbors is calculated as the difference in water surface elevations divided by the grid size. Compared with the kinematic wave approach, this diffusive wave approach allows GSSHA to route water through pits or depressions, and regions of adverse slope. The Manning formula is used to relate flow depth to discharge. Use of the Manning formula implies that the flow is both turbulent and that the roughness is not dependent on flow depth. Neither of these assumptions may be valid on the overland flow plane. While being simple, the method is powerful because it allows calculations to proceed when only portions of the stream network or watershed are flowing. This is an important attribute as rainfall may occur on only a portion of the watershed. Though simple, the methods used in GSSHA are very robust allowing model time steps in the range of 1 to 3 minutes (Downer, 2002a; Downer et al., 2000) and minimal processing of the DEM and channel thalweg elevations.

5.1 Channel Routing

5.1.1 Explicit Channel Routing Formulation

The 1-D channel routing scheme is depicted in Figure 5. Inter-cell flows, Q_i-1/2 and Q_i+1/2 (m³/s) in the longitudinal, x, direction are computed from depths, d, at the n time level using the Manning equation for the head discharge relationship:

(1)

where: n is roughness coefficient, A is the area (m²), R is the hydraulic radius, and S_f is the friction slope, calculated in the x direction as:

(2)

where: S_ox is the land surface slope in the x direction. If negative flow occurs (flow in the upstream direction), the head in the downstream cell is used to calculate the flow as:

(3)

Since the flow direction may change at any point in the stream, especially in ephemeral streams near the beginning of rainfall events, the flow direction is determined around each node and the locally upstream cell properties are used to compute the flow. This simple local determination of the upstream cells prevents crashes in channels with adverse slopes when little or no water is present in the upstream cell. This method also allows better simulations of backwater effects.

Inter-node fluxes are used to calculate the volume, V, in each node as:

(4)

where: q_lat (m²/s) is the amount of lateral inflow from the overland flow cells adjacent to the node, and q_recharge (m²/s) is the exchange between the groundwater and channel. These new volumes are used to compute nodal values of A, d, and wetted perimeter at the n+1 time level. Calculations proceed from the upstream boundary to the downstream boundary.

Several modifications were made in the implementation of the channel routing scheme to accommodate groundwater/channel interactions. These modifications permit continuous interaction between channel nodes and the saturated groundwater cells. The channel routing scheme was modified to allow water to remain in the channel after channel routing ends, and for water to be present in the channel when channel routing begins. Because groundwater may discharge to the stream at anytime, channel routing is initiated anytime a minimum amount of water is in the channel network. If the channel routing scheme indicates there is no flow in the channel, channel routing is halted during periods outside precipitation events. Fluxes between the stream and the groundwater are still computed and adjustments to the stream volumes are made without routing. If groundwater discharges to the stream, channel routing will resume, but at the groundwater time step, which is typically larger than the channel routing time step.

Because GSSHA uses a finite volume representation of channel flow the standard stability criterion, Courant number < 1.0, does not strictly apply. Maintaining stability is dependent on volume changes during each time step. Experience with the scheme indicates that stability can be maintained with a time step limitation that keeps the maximum Courant number everywhere in the network less than 1/6. The user can specify this value using the MAX_COURANT_NUM card, which requires a real number value. The default is 0.04. Groundwater and over bank fluxes can induce instability and additional controls in the channel routing scheme are added to further reduce instability. If the channel routing scheme becomes unstable (negative depth occurs in one or more cells), despite the more restrictive control on the Courant number, the time step is reduced and the channel routing calculations are repeated. The channel routing time step may be automatically reduced to a value as low as 1/1000 of a second. Allowing the time step to become very small during periods of sharp transition allows a larger overall model time step to be used. For each call of the channel routing function the overall model time step is used as the beginning channel routing time step. The time step is only reduced when the stability controls are activated, and then only for that call of the channel routing routine. The current GSSHA version employs a predictor corrector scheme with Picard iterations on the flow area. This has made the scheme much more stable when reservoirs or overbank flooding are simulated.

Several new channel routing features have been added to GSSHA since the publication of the original manual. These features are described in the following SWWRP Tech Note.

Non-ortho Tech Note

Figure 5 – Explicit channel routing scheme

In GSSHA version 7.1.4 and higher, an inertial formulation of the flow equations, as described in Bates et al. (2010), is included as an option. To use this option include the CHANNEL_MOMENTUM card in your project file. This option may be helpful for simulations where channel hydrodynamics are important, such as flooding due to a storm surge. For these cases stability may be increased due the inclusion of inertial terms and a better accounting of the effects of friction. The method may allow longer time steps and potentially reduce simulation times. The user can control the time step coefficient for the momentum formulation by specifying a real value between 0.0 and 1.0 after the CHANNEL_MOMENTUM card. The default value is 0.2. Larger values increase speed, lower values increase stability. The user is referred to Bates et al. (2010) for details: Bates, P. D., Horritt, M. S., and T. J. Fewtrell (2010) A simple inertial formulation of the shallow water equations for efficient two-dimensional flood inundation modelling. J. Hydro. 387 (2010) 33-45.

5.1.2 Boundary Conditions

The upstream boundary condition in each first order link is a no flow condition. The default watershed outlet condition is normal flow, calculated using the channel slope at the watershed outlet. The downstream boundary condition can also be a specified head. The outlet boundary type is changed to a head by using the HEAD_BOUND project card. The boundary depth, m, is specified with the BOUND_DEPTH project card. When the head boundary is specified the depth at the outlet remains at the specified depth for the entire simulation period. Flows entering the outlet cell exit at the same rate. The head boundary condition is desirable when the condition at the outlet of the basin is a known head instead of normal depth. This might occur when the basin empties into a larger water body such as a river, pond, or lake, or when a hydraulic structure is near the watershed outlet. The up-gradient method of computing inter-node discharges allows the head boundary condition influence to propagate upstream.

5.1.3 Initial Conditions

The default for the explicit scheme is to start simulations from the dry bed condition. There is also the ability to save water surface profiles and flows from one simulation and use them as the initial condition of another simulation. The last set of depths and flows from the channel routing code can be saved in one file using the EXPLIC_BACKWATER project card to specify the file name to write depth (m) and discharge (m³ s^-1) for every node in the stream network. Alternately, the values of the water surface and the discharge can be saved in separate files by using the WAT_SURF_PROFILE_OUT and DIS_PROFILE_OUT project cards to provide the file names for output of the surface water profile (m), and the discharge profile, m³ s^-1), respectively. To begin a simulation from the saved values, the EXPLIC_HOTSTART project card is used to prescribe the name of the file containing the initial values of depth and discharge, or the WAT_SURF_PROFILE and DIS_PROFILE project cards are used to prescribe the files of initial water surface profile and discharge, respectively. The options can be used togather; both initial conditions can be read in and ending conditions can be written out. This feature is useful when breaking a simulation into multiple time periods. Initial conditions are also useful for starting groundwater simulations, as well as many other applications.

5.1.4 Describing the Stream Network

Channel routing is specified in the Project file by including the CHAN_EXPLIC or DIFFUSIVE_WAVE cards. Along with one of these two cards, the CHAN_INPUT and STREAM_CELL cards are used to specify the stream input files as described below.

The stream network in GSSHA is described with a series of links and nodes. A node is a single computational element in the stream network. A link is a channel segment comprised of two or more computational nodes. Two files are required to describe the channel network in GSSHA:

• Channel input file (CHAN_INPUT) – ASCII text file describing link connectivity and assigning attributes to every channel node.
• Stream cell input file (STREAM_CELL) - ASCII text file describing the channel and grid connectivity.

The steam cell file describes the topology of the stream network and contain the information needed to provide channel connectivity with the overland flow plane and the saturated groundwater grid. The CHAN_INPUT file contains the information necessary to connect the stream network and assign attributes to each of the stream nodes. WMS can create the STREAM_CELL and CHAN_INPUT file in the required format. Because the stream network is not confined by the grid, the STREAM_CELL input file is particulary complicated and difficult to build without the WMS interface. The rules that govern the creation of channel input file are listed below.

General Rules

1. Looped reaches cannot be simulated.
2. Link types may not be mixed within a link.
3. Internal boundary condition links contain two nodes, one upstream, and one downstream, of the internal boundary condition.

Required Numbering Scheme

1. First order stream links are numbered first.
2. The first link must be numbered 1.
3. Link numbers must always increase in the downstream direction. Upstream links must have smaller numbers than downstream links.
4. Link numbers must change at junctions.
5. Link numbers may not be skipped.
6. Node numbers increase in the downstream direction.
7. The most upstream node in a link is node 1.

Upstream and downstream links overlap one node. The first node in the downstream link is the same as the last node in the upstream link. In the CHAN_INPUT file this node must be explicitly contained in both the upstream and downstream links, such that the upstream link in the CHAN_INPUT will contain one more node computational node than actually exist. An exception to this rule is the most downstream link containing the outlet.

As an example, we'll consider the highly conceptualized stream network shown in the figure below.

Figure 6 – Conceptualized Stream Network

This stream network has 9 stream links, including a lake, link 7, with an outlet weir, link 8.

5.1.4.1 – Channel Input File

The channel input file name is specified with the CHAN_INPUT card in the project file. The channel input file consists of four parts:

• channel input file declaration
• channel routing constants
• channel network connectivity
• link (reach) information

These four sections are assembled in the listed order to produce the channel input file.

5.1.4.1.1 Channel input file declaration

The first line of the channel input file must contain the card "GSSHA_CHAN". This card tells GSSHA that it is reading a channel input file. In the CHAN_INPUT file the first line will look like:

GSSHA_CHAN

5.1.4.1.2 Channel routing constants

The second portion of the channel input file is used to pass physical constants and simulation parameters to the model. There are three channel routing parameters that must be prescribed. The program must be told the number of links, and the maximum number of nodes in any link in the network. In our example, the number of links is 9, and the maximum number of nodes is 9 (in link 6). Recall that all links except the outlet link must have an extra node for connectivity purposes. The total number of links is called "LINKS", and the largest number of nodes in any link in the network is called "MAXNODES". The three paramets, "ALPHA" "BETA" 'THETA" are not currently used. A value of 1.0 is suggested for these parameters.

Table 5 - Channel file input constants

This data constitutes the second portion of channel input file and is arranged into a header which must have the form (note floating point and integers):

ALPHA       1.000000
BETA        1.000000
THETA       1.000000
LINKS       9
MAXNODES    9

5.1.4.1.3 – Channel network connectivity

The third portion of the input file describes the network topology. This is accomplished using a line-input format, one line for each link in the network. Each line contains the identifier "CONNECT" followed by three or more integer values arranged in columns, as shown in Table 6. The total number of columns depends on the number of upstream dependencies, with one for each. The links should appear in ascending numerical order.

Table 6 - Connectivity information format

For our example, the stream connectivity information looks like.

CONNECT    1    4    0
CONNECT    2    4    0
CONNECT    3    4    0
CONNECT    4    5    3    1    2    3
CONNECT    5    6    1    4
CONNECT    6    7    1    5
CONNECT    7    8    1    6
CONNECT    8    9    1    7
CONNECT    9    0    1    8

Note that only the outlet link has 0 downstream dependencies. Links 1,2, and 3 have no upstream dependencies, while link 4 has 3.

5.1.4.1.4 - Link (Reach) information

Following the channel connectivity information is channel link and node information. Parameters are defined for each link in the stream network. Links are listed in ascending order.

5.1.4.1.4.1 - Link declaration

In the CHAN_INPUT file each link is declared with the "LINK" card, which is followed by the link number. For any CHAN_INPUT file the first line after the channel connect information will be:

LINK           1

5.1.4.1.4.2 - Link types

The input that follows the link declaration depends on the link type. Link types may be basic fluvial links, structures, or reservoirs.

5.1.4.1.4.2.1 - Basic fluvial channel links

Channel cross-sections may be either trapezoidal, or natural channels represented by breakpoint cross-sections. Different cross-section types may be mixed in the channel network. Typically, small streams in the upland areas are modeled with trapezoidal cross-sections; larger streams in the lower reaches of the catchment are modeled with break-point cross-sections. Data for the development of break-point cross-sections can be derived from field surveys of cross-section, detailed topographic maps, or high resolution digital elevation models. Coarse resolution DEMs should not be relied upon to develop channel cross-section geometry.

The hydrodynamic channel routing in GSSHA is very sensitive to channel cross-sections, more so than the longitudinal slope. More time should be spent developing accurate cross-sectional description than to getting a highly accurate longitudinal profile. Smooth transitions in channel cross-sectional properties between all connecting fluvial links often play a vital role in the success of simulations. Abrupt changes in cross-sections can lead to numerical mass conservation errors. It may be necessary to create transition links between measured break-point and trapezoidal cross-sections when adjoining links vary greatly in cross-section. The Manning’s roughness parameter used in the model is largely a calibration parameter. Physically realistic literature values may be used as initial guesses and as bound on the calibration (e.g. Chow, 1959; Barnes, 1967).

Basic fluvial links contain nodes. For a basic fluvial link the next line after the link declaration card is the channel node increment, or channel grid size. This is declared by using the "DX" card, followed by the node spacing in the current reach. The node spacing can be different in every link. For our example, link 1 has a node spacing of 107.840396 m. The input in the CHAN_INPUT file is:

DX             107.840396

After setting the node spacing for the link, the type of link is specified. Basic channel reaches may be defined as either "TRAPEZOID" or "BREAKPOINT" for trapezoid and natural cross sections, respectively. These basic reach types may be further described as "SUBSURFACE", for surface water/groundwater interaction, and as "ERODE", for erodible channel if performing sediment transport calculations. Both "SUBSURFACE" and "ERODE" may be combined with with either "TRAPEZOID" or "BREAKPOINT" with the use of underscores. They may used singly or combined. Such as

TRAPEZOID_ERODE
TRAPEZOID_SUBSURFACE

or,

TRAPEZOID_SUBSURFACE_ERODE

The basic channel type TRAPEZOID or BREAKPOINT should always be listed first, followed by the TRAPEZOID and/or SUBSURFACE descriptors.

For our example, Link 1 is a trapezoidal channel. So the input is:

TRAPEZOID

After the link type is defined, the number of nodes is specified with the "NODE" card. For our example, Link 1 has 7 nodes, so the input is:

NODES 7

After defining the number of nodes, information for the first node, node 1, in the link is specified. Required node information is the x and y location of the node and the elevation of the node. These are specified with "X_Y" and "ELEV" cards. For our example the information for link 1 is:

NODE 1
X_Y  1825.200000 933.400000
ELEV 112.694141

After describing the first node, the channel cross section must be defined. How the cross section is defined depends on the type of cross section. In either case, the "X_SEC" card is included next in the CHAN_INPUT file.

X_SEC

5.1.4.1.4.2.1.1 Trapezoidal cross-section

Trapezoidal cross-sections are symmetrical, assumed to be infinitely deep, and have a constant side slope. Trapezoidal cross-sections are defined by:

• Manning’s N - Manning’s roughness coefficient, n (dimensionless)
• bottom width (m)
• channel depth (m)
• side slope, z (change in X with a change in Y of 1).

This information is entered for each trapezoidal link with a series of cards. These cards are:

MANNINGS_N

BOTTOM_WIDTH

BANKFULL_DEPTH

SIDE_SLOPE

Each card is followed by the value. For our example Link 1 is a trapezoidal channel with Manning roughness of 0.03, a bottom width of 2 m, a bankfull depth of 2 m, and a side slope of 2.0. The X_SEC input for this reach is:

XSEC
MANNINGS_N     0.030000
BOTTOM_WIDTH   2.000000
BANKFULL_DEPTH 2.000000
SIDE_SLOPE     2.000000

5.1.4.1.4.2.1.2 Natural cross-section

"BREAKPOINT" cross sections are used to describe natural cross sections. Breakpoint cross sections are defined by the points along the channel cross section. These are entered for each link after the "X_SEC" declarations. Internal to GSSHA, these points and the specified Manning roughness will be used to develop a table of channel information for a specified number of increments between the channel bottom and the bank full depth.

The information in the CHAN_INPUT file is the "MANNINGS_N", the number of XY pairs "NPAIRS" that describe the channel cross section, the number of increments in the channel property table used in GSSHA "NUM_INTERP", followed by the series of XY pairs, specified with the "X1" card.

For our example link 9 is a "BREAKPOINT" cross section, with a roughness of 0.03, and 4 xy points that define the cross section. We are specifying that the cross section be divided into 2000 depth increments for computation purposes. The "X_SEC" input for link 9 appears as:

XSEC
MANNINGS_N     0.030000
NPAIRS         4
NUM_INTERP     2000
X1   -4.000000 2.000000
X1   0.000000 0.000000
X1   2.000000 0.000000
X1   6.000000 2.000000

5.1.4.1.4.2.1.3 Additional Channel Properties (X_SEC qualifiers)

Trapezoidal and breakpoint cross sections represent the two basic fluvial cross section types in GSSHA. For more complex models additional channel properties can be defined after the X_SEC card. The table below describes the additional X_SEC cards.

5.1.4.1.4.2.1.4 Node information

Each fluvial channel link must contain two or more nodes. The nodes are the computational elements used in GSSHA to perform the actual channel routing calculations. Nodes are defined by there XY location, as well as the thalweg elevation, the elevation of the lowest point in the channel cross section.

After describing the channel cross section, XYZ information must be entered for each computational node in the link. For each node the "NODE" card is used to define the node; the x and y position are entered with the "X_Y" card; followed by the node elevation is set with the "ELEV" card. This information is entered for each node in ascending order.

In our example, link 1 has 7 nodes. Information for node 1 was previously described, information for nodes 2 through 7 are listed as follows.

NODE 2
X_Y  1734.200000 991.266667
ELEV 105.977735
NODE 3
X_Y  1643.200000 1049.133333
ELEV 99.186421
NODE 4
X_Y  1552.200000 1107.000000
ELEV 92.323354
NODE 5
X_Y  1461.200000 1164.866667
ELEV 85.407261
NODE 6
X_Y  1370.200000 1222.733333
ELEV 78.464655
NODE 7
X_Y  1279.200000 1280.600000
ELEV 73.846634

It is required that the bed elevation of all channel inverts at each junction be equal. In our example the elevation of the downstream node of links 1, 2 and 3 are equal to the bed elevation of the upstream end of link 4 (Figure 5).

5.1.4.1.4.2 - Structure channel links

Structures are internal boundary conditions within the GSSHA stream routing network. GSSHA allows several types of structures to be defined. The 7 types of structures recognized in GSSHA are:

• Horizontal Weir
• Sag Vertical Curve (Parabolic) Weir
• Circular Culvert
• Box Culvert
• Rating Curve
• Scheduled Release
• Rule Curve

All culverts, for example, lie below a weir formed by a roadway or other embankment. In GSSHA, a hydraulic structure link can consist of any combination of the above. Other structures consist of a number of culverts, some with different dimensions, invert elevations, etc. This allows the user to simulate complex hydraulics from multiple drainage elements.

Hydraulic structures are point models that exist at a channel intersection, at the last node for the upstream link and at the first node of the downstream link, which occupy the same point. The upstream conditions are taken from the second to last node in the upstream link. The downstream conditions are taken from the second node in the downstream link. It is important to note that the bed elevation of the upstream condition is the bed elevation at the junction plus half the difference between the bed elevation at the junction and the bed elevation at the next to last node of the upstream link. The same applies to the downstream condition. To avoid numerical issues, it is advisable to the structure invert be at or above the bed elevation at the respective end.

Horizontal Weir

The "weir" card is used to specify a horizontal broad- or sharp-crested weir of specified length. The figure below shows the geometry of the weir looking downstream. The CREST_LENGTH card specifies the length L shown in the figure in meters.

Other required inputs are the elevation of the crest above the model elevation datum CREST_LOW_ELEV, and the forward- and reverse-flow discharge coefficients, DISCHARGE_COEFF_FORWARD and DISCHARGE_COEFF_REVERSE.

STRUCTTYPE     WEIR
CREST_LENGTH             6.200000
CREST_LOW_ELEV           61.500000
DISCHARGE_COEFF_FORWARD  0.600000
DISCHARGE_COEFF_REVERSE  0.600000

Sag Vertical-Curve Weir

The "sag_weir" card is used to specify a non-symmetric parabolic weir of the type that is commonly created by flow over a roadway crossing a stream, when the flood plain is small. Parabolic vertical curvature of roadways is very common. The geometrical parameters shown in the Figure below are based on the slopes of the two road segments that are to be joined by the parabola, the horizontal distance between the point of vertical curvature (PVC) and point of vertical tangency (PVT), and elevation of the low point on the curve.

<math>Q=C_d L H^{3/2}</math>

The required inputs for sag vertical curve weirs are: CREST_LENGTH L, which is the horizontal distance from PVC to PVT, CREST_LOW_ELEV or E_inv in the figure, the forward and reverse discharge coefficients DISCHARGE_COEFF_FORWARD and DISCHARGE_COEFF_REVERSE, the steepest approach slope entered as a negative number STEEP_SLOPE, and the shallowest approach slope entered as a positive number, SHALLOW_SLOPE.

STRUCTTYPE     SAG_WEIR
CREST_LENGTH             25.000000
CREST_LOW_ELEV           370.000000
CREST_LOW_LOC            0
DISCHARGE_COEFF_FORWARD  0.700000
DISCHARGE_COEFF_REVERSE  0.700000
STEEP_SLOPE              -0.10
SHALLOW_SLOPE            0.08

Circular Culvert

The "round_culvert" card specifies a culvert of circular cross-section. The following six flow types are considered in GSSHA for both circular and rectangular box culverts:

Required inputs for circular culverts are the upstream invert elevation UPINVERT (m), the downstream invert elevation DOWNINVERT (m), the inlet discharge coefficient INLET_DISCH_COEFF, the reverse flow discharge coefficient REV_FLOW_DISCH_COEFF, culvert length LENGTH (m),the diameter of the culvert DIAMETER (m), and the Manning roughness coefficient ROUGH_COEFF. The SLOPE card is not a require input, and is calculated by the code. However, if the length of the culvert is not known, the UPINVERT and SLOPE are used to calculate the DOWNINVERT elevation.

STRUCTTYPE     ROUND_CULVERT
UPINVERT                 409.600000
DOWNINVERT               409.500000
INLET_DISCH_COEFF        0.800000
REV_FLOW_DISCH_COEFF     0.900000
SLOPE                    0.000000
LENGTH                   10.000000
ROUGH_COEFF              0.013000
DIAMETER                 1.220000

As noted above, numerical issues can be avoiding by locating you inverts above the bed elevations at the upstream and downstream condition locations, which are at the mid-point between the junction (structure) node and the next upstream or downstream node, depending on whether it is the UPINVERT or the DOWNINVERT.

Box Culvert

The "rect_culvert" card is used to create a rectangular cross-section culvert. The same flow regimes shown in the figure above are considered.

Required inputs for rectangular culverts are the upstream invert elevation UPINVERT (m), the downstream invert elevation DOWNINVERT (m), the inlet discharge coefficient INLET_DISCH_COEFF, the reverse flow discharge coefficient REV_FLOW_DISCH_COEFF, culvert length LENGTH (m),the width of the rectangular culvert WIDTH (m), the height of the rectangular culvert HEIGHT (m), and the Manning roughness coefficient ROUGH_COEFF. The SLOPE card is not a require input, and is calculated by the code. However, if the length of the culvert is not known, the UPINVERT and SLOPE are used to calculate the DOWNINVERT elevation. If the UPINVERT and DOWNINVERT cards are specified, the slope of the culvert is calculated, and any value input through the SLOPE card is not used.

STRUCTTYPE     RECT_CULVERT
UPINVERT                 413.700000
DOWNINVERT               413.600000
INLET_DISCH_COEFF        0.800000
REV_FLOW_DISCH_COEFF     0.900000
SLOPE                    0.010000
LENGTH                   10.000000
ROUGH_COEFF              0.013000
WIDTH                    0.910000
HEIGHT                   0.670000

As noted above, numerical issues can be avoiding by locating you inverts above the bed elevations at the upstream and downstream condition locations, which are at the mid-point between the junction (structure) node and the next upstream or downstream node, depending on whether it is the UPINVERT or the DOWNINVERT.

Rating Curve

The "rating_curve" card is used to specify a series of discharge vs. water surface elevation points, and can be used to simulate the hydraulic performance of "unusual" hydraulic structures such as irregular weirs or roadway crossings.

Scheduled Release

The "scheduled_release" card is used to enter specified discharges as a function of time, such as those from a reservoir. This functionality is particularly valuable for simulations of past events, where reservoir releases are known.

Rule Curve

The "rule_curve" card is used to enter discrete discharges as a function of water level. These rule curves are often associated with the manual operation of flood control reservoirs. At different times of the year, flood control reservoirs are operated to try to maintain a certain "target" water level in the reservoir. If the water level exceeds the target by a certain amount, the reservoir operator is authorized by the rule curve to release a certain flow rate of water from the reservoir in an attempt to restore the reservoir level to the desired target level.

Composite Structure Link

Suppose that at a certain point a channel is crossed by a parabolic-shaped roadway, and three culverts, one square and two circular are used to pass the normal flow under the roadway. The situation is shown in the following figure.

This hydraulic structure consists of four elements: sag vertical-curve weir, box culvert, and two circular culverts. Notice that as shown in the figure, each element has a different invert elevation, which is allowed.

For the sake of discussion, assume that within the GSSHA channel network, this hydraulic structure is link number 38. The channel input file for this link would appear as:

LINK         38
STRUCTURE
NUMSTRUCTS     4
STRUCTTYPE     RECT_CULVERT
UPINVERT                 2369.790000
DOWNINVERT               2369.540000
INLET_DISCH_COEFF        0.800000
REV_FLOW_DISCH_COEFF     0.900000
SLOPE                    0.010000
LENGTH                   18.000000
ROUGH_COEFF              0.013000
WIDTH                    2.40000
HEIGHT                   1.270000
STRUCTTYPE     ROUND_CULVERT
UPINVERT                 2370.830000
DOWNINVERT               2370.080000
INLET_DISCH_COEFF        0.800000
REV_FLOW_DISCH_COEFF     0.900000
SLOPE                    0.000000
LENGTH                   19.000000
ROUGH_COEFF              0.013000
DIAMETER                 1.220000
STRUCTTYPE     ROUND_CULVERT
UPINVERT                 2370.410000
DOWNINVERT               2369.870000
INLET_DISCH_COEFF        0.800000
REV_FLOW_DISCH_COEFF     0.900000
SLOPE                    0.000000
LENGTH                   19.000000
ROUGH_COEFF              0.013000
DIAMETER                 1.220000
STRUCTTYPE     SAG_WEIR
CREST_LENGTH             45.000000
CREST_LOW_ELEV           2372.860000
CREST_LOW_LOC            0
DISCHARGE_COEFF_FORWARD  0.700000
DISCHARGE_COEFF_REVERSE  0.700000
STEEP_SLOPE              -0.05
SHALLOW_SLOPE            0.035

Note that the order in which the hydraulic elements appear in the STRUCTURE block when there is more than one structural element is not important.

5.1.4.1.4.3 - Reservoir channel links

Reservoirs are defined by specifying that a reservoir is present with the LAKE card in the .cif file. The minimum, initial, and maximum reservoir elevations are specified with the MINWSE, INITWSE, and MAXWSE, cards respectively. The number of cells within the maximum reservoir boundary is specified with the NUMPTS card, followed by a list of the cells in the reservoir domain, defined by row and column values. Reservoirs are dynamic features that can increase/decrease in size, take/release stream nodes, and overland flow cells as defined in Non-ortho Channels Technical Note Non-ortho Tech Note. The resevoir will not go into any cells not listed in the NUMPTS list, but water may flow out of the reservoir back on the overland flow plane. All reservoirs must have a hydraulic structure defined as the reservoir outlet. Other important restrictions are that the reservoir cannot rise above the MAXWSE, move into cells not defined as part of the reservoir, join with other reservoirs, or take in stream nodes not directly in the same branch as the reservoir. Any of these conditions will cause the program to stop, crash, or provide non-sensical results. Large mass balance errors are indicative of problems with reservoirs violating one or more of these conditions. Reservoirs can be kept within the desired confines using embankments. Water can be passed through an embankment holding a reservoir using lowspots along the embankment or using overland rating curves. Water can also pass from one lake to another using these features. Although lakes cannot exceed the maximum water level, they can go dry. If the water level falls below the MINWSE it will stop discharging but water can still seep out of the lake and/or evaporate.

Output from all the reservoirs can be specified with the LAKE_OUTPUT card specified in the project file. Lake elevations and volumes, respectively, will be written for each lake, in the order they are listed in the .cif file, for each time period specified in the HYD_FREQ card, in the file specified with the LAKE_OUTPUT card.

Within the .cif file, a standard lake input would look like the example below:

LINK         7
LAKE           
MINWSE       55.000000
INITWSE      62.100000
MAXWSE       95.000000
NUMPTS       210
18  3     18  2     19  3     17  4     17  3     18  4     16  5     17  5     17  2     19  4     
16  4     16  6     15  6     18  5     17  1     15  7     16  3     20  4     17  6     15  5     
14  8     14  7     19  5     16  2     16  7     15  4     14  6     18  6     16  1     15  8     
13  8     20  5     13  9     15  3     17  7     14  5     14  9     13  7     19  6     15  2     
16  8     12  9     14  4     18  7     13  6     12  10     15  9     15  1     20  6     12  8     
14  3     17  8     13  10     13  5     11  10     19  7     12  7     14  10     14  2     16  9     
11  9     13  4     11  11     18  8     12  6     15  10     10  11     14  1     20  7     12  11     
11  8     13  3     17  9     12  5     13  11     10  10     19  8     11  7     10  12     16  10     
13  2     9  12     14  11     12  4     10  9     18  9     11  12     11  6     20  8     15  11     
13  1     9  11     10  8     17  10     12  3     12  12     8  13     11  5     9  13     19  9     
9  10     13  12     16  11     10  7     12  2     10  13     8  12     18  10     11  4     9  9     
14  12     10  6     7  14     20  9     11  13     12  1     8  11     17  11     15  12     8  14     
11  3     9  8     12  13     10  5     7  13     19  10     8  10     9  14     16  12     9  7     
6  15     11  2     13  13     18  11     7  12     10  4     10  14     6  16     8  9     14  13     
20  10     9  6     7  15     6  14     11  1     17  12     7  11     11  14     5  16     10  3     
8  8     15  13     8  15     19  11     9  5     6  13     12  14     7  10     5  17     5  15     
16  13     8  7     9  15     10  2     18  12     4  17     13  14     6  12     9  4     7  9     
7  16     20  11     10  15     5  14     8  6     4  18     10  1     17  13     14  14     4  16     
6  11     9  3     11  15     7  8     3  18     19  12     8  16     8  5     5  13     15  14     
6  17     6  10     4  15     12  15     9  2     7  7     18  13     3  17     16  14     9  16     
3  19     5  12     8  4     7  17     2  19     13  15     20  12     6  9     4  14     5  18     

Additional Options

Seepage

Seepage, and/or groundwater exchange, is controlled by including the M_LAKE card in the project file, followed by the seepage rate (cm/hr). If multiple lakes are present and different seepage values are desired for each lake then the M_LAKE is included in the .cif file for each lake. Include the M_LAKE card, with the value, before the NUMPTS card. See example below.

Dynamic verses Static Lakes

The default is for lakes to be dynamic features in the channel network and in the overland flow plan. During simulations, channels and overland cells overtaken by the lake are taken out of the stream/overland domain and put in the lake. As the lake level recedes, these nodes/cells are put back. More details are provided in the document linked above. This feature can be turned off by including the STATIC card in the .cif for each lake. Include the STATIC card before the NUMPTS card (See example below). When a lake is STATIC it means that the lake is static in relation to the grid and the stream network, not that the stage, area, and volume are static. For a STATIC lake the footprint of the lake is assumed to be the maximum lake size, defined by the cells specified with NUMPTS. Using STATIC lakes may speed up simulations because complexity is reduced. However, STATIC lakes may increase infiltration and lead to underestimation of backwater effects on both the overland and in the stream network.

Elevation/Volume/Surface Area

The default is to get the stage/storage/area relationship from the grid elevations. However, this may lead to errors if the grid size is large. In any case, the user can specify the stage/storage/area relationship by including the NUMRATE card followed by the number of points in the rating curve. This is followed by the values of elevation (m), volume (m3), and area (m2), one line for each NUMRATE. These data are included AFTER the NUMPTS card and the cells that specify the lake imprint on the grid. See example below.

If all these options were selected for previous example, the input in the .cif file might look like this:

LINK         7
LAKE           
MINWSE       55.000000
INITWSE      62.100000
MAXWSE       95.000000
M_LAKE       10.000000
STATIC
NUMPTS       210
18  3     18  2     19  3     17  4     17  3     18  4     16  5     17  5     17  2     19  4     
16  4     16  6     15  6     18  5     17  1     15  7     16  3     20  4     17  6     15  5     
14  8     14  7     19  5     16  2     16  7     15  4     14  6     18  6     16  1     15  8     
13  8     20  5     13  9     15  3     17  7     14  5     14  9     13  7     19  6     15  2     
16  8     12  9     14  4     18  7     13  6     12  10     15  9     15  1     20  6     12  8     
14  3     17  8     13  10     13  5     11  10     19  7     12  7     14  10     14  2     16  9     
11  9     13  4     11  11     18  8     12  6     15  10     10  11     14  1     20  7     12  11     
11  8     13  3     17  9     12  5     13  11     10  10     19  8     11  7     10  12     16  10     
13  2     9  12     14  11     12  4     10  9     18  9     11  12     11  6     20  8     15  11     
13  1     9  11     10  8     17  10     12  3     12  12     8  13     11  5     9  13     19  9     
9  10     13  12     16  11     10  7     12  2     10  13     8  12     18  10     11  4     9  9     
14  12     10  6     7  14     20  9     11  13     12  1     8  11     17  11     15  12     8  14     
11  3     9  8     12  13     10  5     7  13     19  10     8  10     9  14     16  12     9  7     
6  15     11  2     13  13     18  11     7  12     10  4     10  14     6  16     8  9     14  13     
20  10     9  6     7  15     6  14     11  1     17  12     7  11     11  14     5  16     10  3     
8  8     15  13     8  15     19  11     9  5     6  13     12  14     7  10     5  17     5  15     
16  13     8  7     9  15     10  2     18  12     4  17     13  14     6  12     9  4     7  9     
7  16     20  11     10  15     5  14     8  6     4  18     10  1     17  13     14  14     4  16     
6  11     9  3     11  15     7  8     3  18     19  12     8  16     8  5     5  13     15  14     
6  17     6  10     4  15     12  15     9  2     7  7     18  13     3  17     16  14     9  16     
3  19     5  12     8  4     7  17     2  19     13  15     20  12     6  9     4  14     5  18     
NUMRATE 4
55.0 0.0 0.0
60.0 10.0 100.0
80.0 1000.0 10000.0
95.0 10000.0 100000.0  

5.1.4.1.5 Assembling the Channel Input File

In this section, the three portions of the channel input file are put together to form the channel input file for our example project.

GSSHA_CHAN
ALPHA       1.000000
BETA        1.000000
THETA       1.000000
LINKS       9
MAXNODES    9
CONNECT    1    4    0
CONNECT    2    4    0
CONNECT    3    4    0
CONNECT    4    5    3    1    2    3
CONNECT    5    6    1    4
CONNECT    6    7    1    5
CONNECT    7    8    1    6
CONNECT    8    9    1    7
CONNECT    9    0    1    8

LINK           1
DX             107.840396
TRAPEZOID
NODES          7
NODE 1
X_Y  1825.200000 933.400000
ELEV 112.694141
XSEC
MANNINGS_N     0.030000
BOTTOM_WIDTH   2.000000
BANKFULL_DEPTH 2.000000
SIDE_SLOPE     2.000000
NODE 2
X_Y  1734.200000 991.266667
ELEV 105.977735
NODE 3
X_Y  1643.200000 1049.133333
ELEV 99.186421
NODE 4
X_Y  1552.200000 1107.000000
ELEV 92.323354
NODE 5
X_Y  1461.200000 1164.866667
ELEV 85.407261
NODE 6
X_Y  1370.200000 1222.733333
ELEV 78.464655
NODE 7
X_Y  1279.200000 1280.600000
ELEV 73.846634

LINK           2
DX             97.619562
TRAPEZOID
NODES          7
NODE 1
X_Y  920.200000 1743.400000
ELEV 109.559260
XSEC
MANNINGS_N     0.030000
BOTTOM_WIDTH   2.000000
BANKFULL_DEPTH 2.000000
SIDE_SLOPE     2.000000
NODE 2
X_Y  980.033333 1666.266667
ELEV 102.920801
NODE 3
X_Y  1039.866667 1589.133333
ELEV 96.282341
NODE 4
X_Y  1099.700000 1512.000000
ELEV 90.355169
NODE 5
X_Y  1159.533333 1434.866667
ELEV 84.165096
NODE 6
X_Y  1219.366667 1357.733333
ELEV 77.843573
NODE 7
X_Y  1279.200000 1280.600000
ELEV 73.846634

LINK           3
DX             119.589339
TRAPEZOID
NODES          3
NODE 1
X_Y  1445.400000 1452.600000
ELEV 80.258959
XSEC
MANNINGS_N     0.030000
BOTTOM_WIDTH   2.000000
BANKFULL_DEPTH 2.000000
SIDE_SLOPE     2.000000
NODE 2
X_Y  1362.300000 1366.600000
ELEV 75.890504
NODE 3
X_Y  1279.200000 1280.600000
ELEV 73.846634

LINK           4
DX             97.528863
TRAPEZOID
NODES          6
NODE 1
X_Y  1279.200000 1280.600000
ELEV 73.846634
XSEC
MANNINGS_N     0.030000
BOTTOM_WIDTH   2.000000
BANKFULL_DEPTH 2.000000
SIDE_SLOPE     2.000000
NODE 2
X_Y  1214.524434 1207.600341
ELEV 68.100627
NODE 3
X_Y  1149.848867 1134.600682
ELEV 64.624511
NODE 4
X_Y  1085.173301 1061.601023
ELEV 61.082366
NODE 5
X_Y  1020.497734 988.601364
ELEV 57.443587
NODE 6
X_Y  955.822168 915.601705
ELEV 53.724277

LINK           5
DX             102.425766
TRAPEZOID
NODES          3
NODE 1
X_Y  955.822168 915.601705
ELEV 53.724277
XSEC
MANNINGS_N     0.030000
BOTTOM_WIDTH   2.000000
BANKFULL_DEPTH 2.000000
SIDE_SLOPE     2.000000
NODE 2
X_Y  881.349951 845.281535
ELEV 49.777956
NODE 3
X_Y  806.877734 774.961364
ELEV 45.831635

LINK           6
DX             102.425766
TRAPEZOID
NODES          9
NODE 1
X_Y  806.877734 774.961364
ELEV 45.831635
XSEC
MANNINGS_N     0.030000
BOTTOM_WIDTH   2.000000
BANKFULL_DEPTH 2.000000
SIDE_SLOPE     2.000000
NODE 2
X_Y  732.405518 704.641194
ELEV 42.106178
NODE 3
X_Y  657.933301 634.321023
ELEV 38.368474
NODE 4
X_Y  583.461084 564.000853
ELEV 34.624646
NODE 5
X_Y  508.988867 493.680682
ELEV 30.880817
NODE 6
X_Y  434.516650 423.360512
ELEV 30.234555
NODE 7
X_Y  360.044434 353.040341
ELEV 27.028553
NODE 8
X_Y  285.572217 282.720171
ELEV 21.302121
NODE 9
X_Y  211.100000 212.400000
ELEV 14.343808

LINK         7
RESERVOIR           
MINWSE       55.000000
INITWSE      62.100000
MAXWSE       95.000000
NUMPTS       210
18  3     18  2     19  3     17  4     17  3     18  4     16  5     17  5     17  2     19  4     
16  4     16  6     15  6     18  5     17  1     15  7     16  3     20  4     17  6     15  5     
14  8     14  7     19  5     16  2     16  7     15  4     14  6     18  6     16  1     15  8     
13  8     20  5     13  9     15  3     17  7     14  5     14  9     13  7     19  6     15  2     
16  8     12  9     14  4     18  7     13  6     12  10     15  9     15  1     20  6     12  8     
14  3     17  8     13  10     13  5     11  10     19  7     12  7     14  10     14  2     16  9     
11  9     13  4     11  11     18  8     12  6     15  10     10  11     14  1     20  7     12  11     
11  8     13  3     17  9     12  5     13  11     10  10     19  8     11  7     10  12     16  10     
13  2     9  12     14  11     12  4     10  9     18  9     11  12     11  6     20  8     15  11     
13  1     9  11     10  8     17  10     12  3     12  12     8  13     11  5     9  13     19  9     
9  10     13  12     16  11     10  7     12  2     10  13     8  12     18  10     11  4     9  9     
14  12     10  6     7  14     20  9     11  13     12  1     8  11     17  11     15  12     8  14     
11  3     9  8     12  13     10  5     7  13     19  10     8  10     9  14     16  12     9  7     
6  15     11  2     13  13     18  11     7  12     10  4     10  14     6  16     8  9     14  13     
20  10     9  6     7  15     6  14     11  1     17  12     7  11     11  14     5  16     10  3     
8  8     15  13     8  15     19  11     9  5     6  13     12  14     7  10     5  17     5  15     
16  13     8  7     9  15     10  2     18  12     4  17     13  14     6  12     9  4     7  9     
7  16     20  11     10  15     5  14     8  6     4  18     10  1     17  13     14  14     4  16     
6  11     9  3     11  15     7  8     3  18     19  12     8  16     8  5     5  13     15  14     
6  17     6  10     4  15     12  15     9  2     7  7     18  13     3  17     16  14     9  16     
3  19     5  12     8  4     7  17     2  19     13  15     20  12     6  9     4  14     5  18     

LINK         8
STRUCTURE
NUMSTRUCTS     1
STRUCTTYPE     WEIR
CREST_LENGTH             0.200000
CREST_LOW_ELEV           61.500000
CREST_LOW_LOC            0
DISCHARGE_COEFF_FORWARD  0.600000
DISCHARGE_COEFF_REVERSE  0.600000

LINK           9
DX             97.920932
BREAKPOINT
NODES          4
NODE 1
X_Y  211.100000 212.400000
ELEV 14.343808
XSEC
MANNINGS_N     0.030000
NPAIRS         4
NUM_INTERP     2000
X1   -4.000000 2.000000
X1   0.000000 0.000000
X1   2.000000 0.000000
X1   6.000000 2.000000
NODE 2
X_Y  140.866667 144.166667
ELEV 7.691532
NODE 3
X_Y  70.633333 75.933333
ELEV 4.305107
NODE 4
X_Y  0.400000 7.700000
ELEV 2.116070

5.1.4.2 - Grid connectivity STREAM_CELL file

The topology relationships between the channel and the overland grid are defined in the STREAM_CELL file. This file contains two types of information, the location of the channel nodes in the grid plane in terms of I (row) and J (column) grid location, and the percentage of each stream node in the listed grid cell. Our example stream network is shown over the grid in the following figure.

Figure 7 – Stream network over grid

As is seen in the figure, every node falls into a grid cell. In addition, every node length falls in one or more grid cells. A certain percentage of each computational node falls within any grid cell that contains a channel node. This information is contained in the STREAM_CELL file. The STREAM_CELL file contains the following information.

5.1.4.2.1 - STREAM_CELL file identifier

Every STREAM_CELL file begins with the card

GRIDSTREAMFILE

5.1.4.2.2 - Number of stream cells

The second line in the STREAM_CELL file is the total number of grid cells with streams "STREAMCELLS". In our example the total number of cells with at least part of stream node is 45. So the second line in our STREAM_CELL file is:

STREAMCELLS 45

5.1.4.2.3 - Link/node/grid matching information

Every cell with some part of a stream node in it is identified with the "CELLIJ" card, which gives the row and column of the cell to be described. This is followed by the number of nodes in the cell with the "NODES" card. Followed with the "LINKNODE" card for every node in the grid cell. The "LINKNODE" card specifies the link number, node number, and the percentage of the node within the grid cell. In our example, about 27% of link 1, node 1, falls in grid cell with row 11 column 19. The remaining portion falls in grid cell 11 18, along with about 15% of link 1 node 2. This information in the STREAM_CELL file looks like:

CELLIJ        11 19
NUMNODES      1
LINKNODE      1 1 0.276923
CELLIJ        11 18
NUMNODES      2
LINKNODE      1 1 0.723077
LINKNODE      1 2 0.150922

This information is repeated for every cell containing portions of a channel node. The information is listed is ascending link/node order. For our example, the entire STREAM_CELL file looks like:

GRIDSTREAMFILE
STREAMCELLS   45
CELLIJ        11 19
NUMNODES      1
LINKNODE      1 1 0.276923
CELLIJ        11 18
NUMNODES      2
LINKNODE      1 1 0.723077
LINKNODE      1 2 0.150922
CELLIJ        10 18
NUMNODES      1
LINKNODE      1 2 0.224903
CELLIJ        10 17
NUMNODES      2
LINKNODE      1 2 0.624176
LINKNODE      1 3 0.474725
CELLIJ        10 16
NUMNODES      1
LINKNODE      1 3 0.404307
CELLIJ        9 16
NUMNODES      2
LINKNODE      1 3 0.120968
LINKNODE      1 4 0.573626
CELLIJ        9 15
NUMNODES      2
LINKNODE      1 4 0.426374
LINKNODE      1 5 0.607143
CELLIJ        8 15
NUMNODES      1
LINKNODE      1 5 0.065385
CELLIJ        8 14
NUMNODES      2
LINKNODE      1 5 0.327473
LINKNODE      1 6 0.771429
CELLIJ        8 13
NUMNODES      5
LINKNODE      1 6 0.228571
LINKNODE      2 6 0.251513
LINKNODE      3 2 0.225581
LINKNODE      4 1 1.000000
LINKNODE      4 2 0.104115
CELLIJ        3 10
NUMNODES      1
LINKNODE      2 1 0.562662
CELLIJ        4 10
NUMNODES      2
LINKNODE      2 1 0.437338
LINKNODE      2 2 0.333705
CELLIJ        4 11
NUMNODES      1
LINKNODE      2 2 0.525414
CELLIJ        5 11
NUMNODES      3
LINKNODE      2 2 0.140882
LINKNODE      2 3 1.000000
LINKNODE      2 4 0.005014
CELLIJ        5 12
NUMNODES      1
LINKNODE      2 4 0.150561
CELLIJ        6 12
NUMNODES      2
LINKNODE      2 4 0.844425
LINKNODE      2 5 0.452031
CELLIJ        7 12
NUMNODES      1
LINKNODE      2 5 0.224292
CELLIJ        7 13
NUMNODES      3
LINKNODE      2 5 0.323677
LINKNODE      2 6 0.748487
LINKNODE      3 2 0.024719
CELLIJ        6 15
NUMNODES      1
LINKNODE      3 1 0.546330
CELLIJ        6 14
NUMNODES      1
LINKNODE      3 1 0.065298
CELLIJ        7 14
NUMNODES      2
LINKNODE      3 1 0.388372
LINKNODE      3 2 0.749699
CELLIJ        9 13
NUMNODES      1
LINKNODE      4 2 0.120459
CELLIJ        9 12
NUMNODES      2
LINKNODE      4 2 0.775426
LINKNODE      4 3 0.473984
CELLIJ        10 12
NUMNODES      1
LINKNODE      4 3 0.296769
CELLIJ        10 11
NUMNODES      2
LINKNODE      4 3 0.229247
LINKNODE      4 4 0.843854
CELLIJ        11 11
NUMNODES      2
LINKNODE      4 4 0.156146
LINKNODE      4 5 0.316932
CELLIJ        11 10
NUMNODES      2
LINKNODE      4 5 0.683068
LINKNODE      5 1 0.221867
CELLIJ        12 10
NUMNODES      1
LINKNODE      5 1 0.527704
CELLIJ        12 9
NUMNODES      2
LINKNODE      5 1 0.250430
LINKNODE      5 2 0.643934
CELLIJ        13 9
NUMNODES      2
LINKNODE      5 2 0.356066
LINKNODE      6 1 0.092353
CELLIJ        13 8
NUMNODES      2
LINKNODE      6 1 0.907647
LINKNODE      6 2 0.066001
CELLIJ        14 8
NUMNODES      1
LINKNODE      6 2 0.369135
CELLIJ        14 7
NUMNODES      2
LINKNODE      6 2 0.564864
LINKNODE      6 3 0.488068
CELLIJ        15 7
NUMNODES      1
LINKNODE      6 3 0.289850
CELLIJ        15 6
NUMNODES      2
LINKNODE      6 3 0.222082
LINKNODE      6 4 0.910135
CELLIJ        16 6
NUMNODES      2
LINKNODE      6 4 0.089865
LINKNODE      6 5 0.120701
CELLIJ        16 5
NUMNODES      2
LINKNODE      6 5 0.879299
LINKNODE      6 6 0.332202
CELLIJ        17 5
NUMNODES      1
LINKNODE      6 6 0.131281
CELLIJ        17 4
NUMNODES      2
LINKNODE      6 6 0.536516
LINKNODE      6 7 0.754269
CELLIJ        18 4
NUMNODES      1
LINKNODE      6 7 0.051997
CELLIJ        18 3
NUMNODES      3
LINKNODE      6 7 0.193734
LINKNODE      6 8 1.000000
LINKNODE      9 1 0.158045
CELLIJ        18 2
NUMNODES      1
LINKNODE      9 1 0.023685
CELLIJ        19 2
NUMNODES      2
LINKNODE      9 1 0.818271
LINKNODE      9 2 0.581870
CELLIJ        19 1
NUMNODES      1
LINKNODE      9 2 0.065419
CELLIJ        20 1
NUMNODES      2
LINKNODE      9 2 0.352711
LINKNODE      9 3 1.000000

WMS builds the STREAM_CELL file from GIS like inputs developed in the WMS interface. Building this file without the WMS interface is extremely tedious. Without the interface, it would be most practical to employ a stream network that conforms to the grid on a one to one basis. That is, the stream nodes and the grids cells are the same size, and the stream nodes do not cut across grid cells.

5.1.5 Longitudinal Channel Profile Smoothing

When available, channel cross-sectional geometry and longitudinal profile data collected from extensive field surveys can be input into GSSHA. In lieu of detailed stream bed profiles, the DEM may be used as an indicator of channel slope, assuming that over large reaches, the slope of the channel thalweg is on the order of the same slope as the land surface. This assumption is not true in the case of extremely sinuous channels. Channel slopes extracted from the grid elevations are skewed by errors in the DEM and from errors in interpolating from the DEM to the resolution of the grid. These errors may result in regions of adverse channel slope, undulating streambed elevations, or sharp transitions in channel slope. These conditions will require a reduction in computational time step for the explicit diffusive wave routing method. It may be impossible to simulate channel flow with in a channel with many uncorrected DEM derived errors.

Extensive regions of adverse slope are rare in natural channels. It is recommended that regions of adverse slopes be removed from the GSSHA channel network, unless field observations or obvious geologic controls indicate that they are justified. Smoothing of the channel thalweg profile should be done to remove stream sections of adverse slope, undulating streambed elevations, zero slope, and sharp slope transitions.

Ogden et al (1994) developed a method of smoothing stream channel profiles based on the realization that in regions of concave topography near streams, DEMs are likely to be positively biased. For such conditions filling depressions along stream paths is probably erroneous. A more accurate approach is to remove higher regions to connect depressions. WMS provides tools to smooth the channel thalweg profile both manually and with an automated version of the Ogden et al. (1994) algorithm.

Changes to the stream thalweg elevation alone will not result in proper simulations of channel-groundwater exchange or channel-overbank exchange. To properly simulate these processes, the correct elevations, or at least the correct difference in elevation, of the channel thalweg and overbank area are needed. Upon initialization GSSHA will check in channel thalweg profile for problems and report them, if any, in the "check_thalweg" file.

5.1.6 Losing and Gaining Streams

The calculation of stream losses from streams in arid regions with coarse textured substrate can be specified by including the STREAM_LOSS card in the project file. When simulating either a static or moving water table, streams may be either losing or gaining streams, and the STREAM_LOSS card need not be included in the project file. In either case, the hydraulic conductivity (K_rb) (cm hr^-1) and thickness (M_rb) (cm) of the substrate must be specified. Uniform values can be specified in the project file with the K_RIVER and M_RIVER cards. Including these values in the project file results in every node in every stream link, regardless of specified type, having these uniform values of stream bed parameters. Alternately, the values can be specified in the CHAN_INPUT file for each link specified as as "SUBSURFACE", as described above.

Only channel links identified in the CHAN_INPUT file as "SUBSURFACE" will gain and lose. Links can be identified as "SUBSURFACE" with the "SUBSURFACE" card. Links can also be identified as "SUBSURFACE" by appending "SUBSURFACE" to the link type declaration, i.e. "TRAPEZOID_SUBSURFACE".

Any "SUBSURFACE" link without a specified "K_RIVER" and "M_RIVER" card will be assigned the uniform values of K_RIVER and M_RIVER. The default values for M_RIVER and K_RIVER are 1.0 and 0.0, respectively, such that unless values are specified in either the project file or the CHAN_INPUT file, there will be no channel/subsurface interaction.

Calculation of the flux term is discussed in Section 8.5.1. The explicit channel routing code was modified to allow water to remain in the channel after channel routing ends, and for there to be a starting volume of water in each cell at the beginning of channel routing. When simulating groundwater, channel routing continues after a storm event ends until no water is moving in the channel or only a very small amount of water is left in the stream channels. If the stream discharges to the ground; style="ater after channel routing has ceased, the volume of the stream is updated each groundwater update call. No channel routing is performed under these conditions. However, if the groundwater is discharging to the stream, channel routing will resume, but at the groundwater time step.

5.1.7 Sediment Transport in Channels

Sediment transport in the channel network can be performed as described in Section 10. Sediment transport in the channel network requires links be identified as "ERODE" in the CHAN_INPUT file. For each link identified as "ERODE", the maximum erosion (m) should be specified with the "MAX_EROSION" card in the CHAN_INPUT file. The default value is 0.0. Typically, if sediment routing is desired, all stream links should be identified as "ERODE", with either the "ERODE" card, or appending "ERODE" to the channel type,i.e. "TRAPEZOID_ERODE".

5.2 Overland Flow Routing

5.2.1 Overland Flow Routing Formulation

Overland flow in GSSHA employs the same methods described for 1-D channel routing, except the calculations are made in two dimensions. Flow is routed in two orthogonal directions in each grid cell during each time step. The watershed boundary represents a no flow boundary for the overland flow routing and when a grid cell lies on the watershed boundary, flow is not routed across the boundary. There is an option in GSSHA to impose head boundaries along the watershed boundary and compute flows across the overland boundary. In GSSHA, Δx = Δy. Inter-cell fluxes in the x and y directions, p and q, respectively, are computed in cell ij from the depth, d_ij, at the n^th time level using the Manning equation for the head discharge relationship in the x and y directions, respectively, as

(6)

(7)

Depths in each cell are calculated at the n+1 time level based on the flows for each cell (Julien and Saghafian, 1991):

(8)

In addition to this original formulation in the CASC2D model, two additional methods of solving the equations have been added, an alternating direction explicit scheme (ADE) and an ADE scheme with an additional predictor-corrector step (ADE-PC) (Downer 2002a; Downer et al., 2000). Both the ADE and ADE-PC methods employ the up-gradient difference technique, Equation 3, for flows in the upstream direction (Downer, 2002a). The explicit scheme is being deprecated in v8.0 and will not be available in future verions. Fluxes other than inter-cell fluxes, direct evaporation (DET), infiltration, exfiltration, are accounted for before overland routing is computed.

In the ADE method, inter-cell flows are first calculated in the x direction according to Equation 6. Depths in each row are updated based on the flows in the x direction:

(9)

Inter-cell flows in the y direction are computed using the updated depths:

(10)

Depths in each column are updated based on the flows in the y direction:

(11)

With the ADE-PC method additional steps are added to improve accuracy and stability. As before, during each sweep, by rows or by columns, an estimate of heads is made based on the calculated flows, Equations 9 and 11. Next, using the updated depths, updated estimates of flow are computed at the n+1 time level

(12)

The original flows and the updated flows are then averaged to come up with an estimate of flows for the time step:

(13)

These flows are then used to update the original depths, Equations 9 and 11. This procedure is essentially the MacCormack method (MacCormack, 1969) except up-gradient differences are used in both the predictor and corrector steps. A similar method was successfully implemented by Wang and Hjelmfelt (1998).

In benchmark tests using the three methods: original explicit EXPLICIT, ADE, and ADE-PC, in an contrived watershed consisting of two converging planes (open book), the Goodwin Creek Experimental Watershed (GCEW) (Senarath et al., 2000), and Poplar Creek (Downer et al., 2002a), the ADE and ADE-PC methods ran with significantly larger time steps (Downer et al., 2000). Depending on the test case, time steps could be increased from 20% to 240% with commensurate decreases in simulation times.

The routing scheme is selected using the OVERTYPE card. The default value is ADE. The most efficient scheme to use depends on the particular watershed. The ADE-PC scheme can generally handle rougher terrain and typically requires less smoothing of the DEM, but the additional computational steps result in greater computation time, unless use of the ADE-PC scheme permits substantially greater time steps than with one of the other two methods. For smoothed DEMs or in watersheds with smoother terrain, the ADE and EXPLICIT methods usually can be employed, with a resulting savings in execution time.

To improve stability, the timestep in both the EXPLICIT and ADE schemes is variable. The time step is not variable in the ADE-PC method. The addition of a variable timestep has allows for increased overall model timestep when using either the EXPLICIT and ADE methods. These new methods have not been benchmarked against the ADE-PC method.

in GSSHA version 8.0 the EXPLICIT option is deprecated and will not be available in future versions.

In GSSHA version 7.1.1 and higher, an inertial formulation of the flow equations, as described in Bates et al. (2010), is included as an option. To use this option include the OVERLAND_MOMENTUM card in your project file along with ADE specified with the OVERTYPE card . This option may be helpful for simulations where overland hydrodynamics are important, such as flooding due to a storm surge. For these cases stability may be increased due the inclusion of inertial terms and a better accounting of the effects of friction. The method may allow longer time steps and potentially reduce simulation times. In addition, in GSSHA v7.1.1 and higher, a strict overland time step for the diffusive wave formulation, as described in Bates et al. (2010) has been included. This option is specified by including the OVERLAND_STRICT_DT card in the project file. This option also only works with the ADE method. While use of this option may increase stability for high overland flows, the timestep is very restrictive compared to the normal criteria used for the ADE method. If instability arises with use of the ADE method, it may be better to reduce the overall timestep, use the inertial formulation for flow OVERLAND_MOMENTUM, or use the ADE-PC method, in lieu of the OVERLAND_STRICT_DT option. For versions 7.14 and beyond, the user can control the time step coefficient for the momentum formulation by specifying a real value between 0.0 and 1.0 after the OVERLAND_MOMENTUM card. The default value is 0.2. Larger values increase speed, lower values increase stability. The user is referred to Bates et al. (2010) for details: Bates, P. D., Horritt, M. S., and T. J. Fewtrell (2010) A simple inertial formulation of the shallow water equations for efficient two-dimensional flood inundation modelling. J. Hydro. 387 (2010) 33-45.

Special overland flow cell types that may defined and result in a special type of overland routing are wetlands described in the GSSHA primer here: http://www.gsshawiki.com/Wetlands:Wetlands

And overland lakes, which are lakes not connected to the stream network. Overland lakes are specified with the OVERLAND_LAKE_MAP card that specifies the overland lake mask. The overland lake mask is a map with zeros and integer values, like the MASK_MAP except the integer values occur where overland lakes are desired. This mask map can be created in WMS by making a coverage in WMS for overland lakes and then assigning integer values to polygons that define the lakes, and then creating an index map from that WMS coverage. The significance of an overland lake is that no overland routing occurs within the overland lake cells. Other overland processes, such as infiltration, still occur.

5.2.2 Overland Flow Hydraulic Roughness

The GSSHA model requires that Manning roughness coefficients be assigned to every cell in the watershed mask. There are three ways to specify the hydraulic roughness of the overland flow planes in GSSHA. The first method is to apply a constant value over the entire watershed through the MANNING_N project file card. The second method is to use the MAPPING_TABLE to assign roughness coefficients using tabled values referenced to an index map. The third method is to produce a GRASS ASCII map of roughness coefficient, and provide the name of this map to GSSHA using the ROUGHNESS project file card. Table 8 provides typical values of the Manning roughness coefficient for overland flow over various surfaces:

Table 8 - Values of overland flow roughness coefficient Notes: ^aEngman (1986), ^bDowner (2002b), ^cSenarath et al (2000), ^dHEC (1985)

These values should be considered guidelines. Manning roughness coefficients are typically assigned from literature values and then adjusted through calibration. Additional sources of literature values are Liong et al. (1989), Engman (1986), and Ree et al. (1977). If calibrated values differ significantly from published values, there may be appropriate justification.

Typical Manning roughness coefficient values for open channel flow are considerably smaller than overland flow values because of deeper flow depths in the channel. In cases where the overland flow may become very deep, such as simulating a tidal surge, the user may want to use depth varying overland roughness. GSSHA is currently formulated to calculate the overland roughness at any depth (d) using the formula

n_d=n₀e^-Bd

where n_d is the Manning roughness at depth (d), n₀ is the specified roughness in the ROUGHNESS mapping table, and B is an exponent specified in the ROUGH_EXP mapping table. The default value for B is zero, which results in a static value roughness with depth. To modify this relationship, a value of B is specified for each roughness category in the ROUGH_EXP table. So that the tables will look like

ROUGHNESS "unif"
NUM_IDS 1
ID    DESCRIPTION1                            DESCRIPTION2                                   ROUGH   
1     Roughness ID                                                                           0.300000 

ROUGH_EXP "unif"
NUM_IDS 1
ID    DESCRIPTION1                            DESCRIPTION2                                   ROUGH_EXP
1     Roughness ID                                                                           0.5

Positive values of B result in the roughness decreasing with depth, common; negative values result in roughness increasing with depth, uncommon. A typical value is 0.5. The exponent may require calibration.

5.2.3 Runoff Retention

Natural land surfaces contain micro-topography, small depressions, that retain water prior to runoff. The water held in the grid cell, or retention storage, never becomes direct runoff and can only be removed from the land surface as infiltration or direct evaporation. In certain regions, the retention storage can be significant. Retention storage is input as a depth (mm) in each grid cell and may be optionally input to GSSHA as:

• a uniform value using the RETENTION card,
• as a table of values related to index maps by using the MAPPING_TABLE project card, and the RETEN_DEPTH card without specifying a file name or
• as an ASCII GRASS map through the use of the RETEN_DEPTH project file card with a specified file name.

5.2.4 Specifying Initial Depths on the Watershed

Initial depths on the overland flow plane at the beginning of the simulation may be specified by use of the OV_INIT_DEPTH project file card. This card is used to specify an ASCII GRASS map containing initial overland depth values, m. The initial depth can also be specified with the READ_OV_HOTSTART card. A map of the final overland flow values from a simulation can be written out using the WRITE_OV_HOTSTART card. This map can then be used to hotstart subsequent simulations. A single value of initial water surface elevation (m) can be specified using INIT_ELEV_HEAD card. The water surface elevation in every cell with an elevation less the specified elevation will be set to this value and the depths computed from the specified elevation and the cell elevation.

5.2.5 Simulations without Channel Routing

It might be desirable to perform simulations without channel routing in very small watersheds lacking a defined channel network, or in the beginning stages of developing a GSSHA model of a watershed. It is always prudent to build a GSSHA model one process at a time. Getting just the overland flow portion of the model to run is always the first step in building a model. For these reasons, GSSHA may be used without channel routing. Because the channel network normally provides the outlet point, the overland cell containing the outlet must be specified during simulations without channel routing. Normally this grid cell will have the lowest elevation in the watershed. The row and column containing the outlet grid cell are specified using the OUTROW and OUTCOL cards, while the slope of the outlet grid cell is specified using the OUTSLOPE card. If channel routing is enabled through the inclusion of CHAN_EXPLIC in the project file, the OUTROW, OUTCOL, and OUTSLOPE cards are not required, and ignored if present.

Rainfall is always a required input. Rainfall may be input as spatially and temporally uniform, at a specified rate for a specified duration, for a single event, or rainfall may be input as spatially and temporally varying for any number of rainfall events. The second method requires a separate rainfall data input file.

6.1 Spatially and Temporally Uniform Precipitation

Spatially and temporally uniform rainfall is specified with project file cards. Place the PRECIP_UNIF card in the project file, and specify the rainfall rate (mm/hr) and duration (minutes) with the RAIN_INTENSITY and RAIN_DURATION cards, respectively. Since GSSHA normally determines the starting date and time of a simulation from the rainfall input file, this information must be provided in the project file using the START_DATE and START_TIME cards. For example, the following five lines entered in a project file specify a rainfall rate of 27 mm/h, for a duration of 60 minutes, starting on 12 June 1996 at 13:30. Note that the date/time format in GSSHA input is always year, month, day, hour, and minute.

The use of spatially- and temporally-uniform rainfall is mutually exclusive with continuous simulations using the LONG_TERM project file card. Do not specify both PRECIP_UNIF and LONG_TERM in the same project file.

6.2 Spatially and Temporally Varied Precipitation

Assignment of spatially and temporally varied rainfall requires the creation of a rainfall input file. Like the project file, the rainfall input file is card based. The precipitation input file can contain multiple “events’. The following rainfall cards, Table 9, are recognized by GSSHA, and valid for use only in the precipitation input file (not in the project file):

Table 9 – Types of rainfall inputs

The following information on developing rainfall input should also be noted:

• In a given EVENT, the rainfall data source type (GAGES, RADAR, RATES, ACCUM) may NOT change.
• The data source type may change from one EVENT to the next.
• The number and location of rain gages may change from one event to the next.
• If only one gage is present, rainfall interpolation is impossible. The location of the gage is irrelevant and the gage coordinates are ignored. Rainfall is applied uniformly in space. This provides a means to apply a temporal distribution of rain in a spatially uniform fashion. Such temporally varying, spatially uniform rainfall distributions are commonly used in flood frequency analysis, i.e. TP40.
• A separate line with its own time of recording is used to input each instance of rainfall, allowing varying temporal resolution rainfall data to be input. This feature is particularly useful when using radar-rainfall estimates, since the temporal resolution can vary considerably.
• Avoid using precipitation data with temporal resolution coarser than 1 hour.
• The finest temporal resolution of GSSHA rainfall input is 1 minute. Rainfall rates change on integer minutes. Seconds are not allowed in the time field.
• The COORD card must be followed by the easting and northing of the rain gage. Easting and northing must be in coordinates in the same frame of reference as the header of all ASCII GRASS input maps. If the gage and map coordinates are not in the same system, the gages will not be placed at the correct location in relation to the watershed.
• If used, the optional strings must be enclosed in double quotes “like this”
• Tabs or spaces are used to delimit the file. DO NOT USE COMMAS.
• To improve readability, line feeds are allowed in the data file.
• If radar rainfall files are used they should be in the same projection as the model, which should be a UTM projection.

The following is an example of a multiple-event precipitation input file using radar-rainfall estimates for the first event, rain gage data for the second event, and radar ascii text files for the third event. (unlikely but illustrative):

EVENT “Event of 30 June 1995- rainfall stops on July 1st”
NRPDS 5 
NRGAG 3
COORD 205150.0 4750212.0 "center of radar pixel #1" 
COORD 205045.0 4750104.0 "center of radar pixel #2" 
COORD 205320.0  4751173.0 "center of radar pixel #3" 
RADAR 1995 06  30  22  56  0.00  0.00  0.00 
RADAR 1995 06  30  23  18  10.75  2.25  5.80 
RADAR 1995 06  30  23  39  21.16  1.80  41.50 
RADAR 1995 06  30  23  57  12.13  20.90  20.70 
RADAR 1995 07  01  00  09  11.71  16.50  2.30 

EVENT "Event of 4 July 1995- new raingage network data" 
NRPDS  4 
NRGAG  4 
COORD 204555.0  4751268.0 "location of raingage #1" 
COORD 205642.0  4750491.0 "location of raingage #2" 
COORD 205921.0  4750330.0 "location of raingage #3" 
COORD 206170.0  4749611.0 "location of raingage #4" 
GAGES 1995 07  04  09  47  0.0  0.0  0.0  0.0 
GAGES 1995 07  04  10  01  38.0  2.0  0.0  0.0 
GAGES 1995 07  04  10  16  16.0  14.0  3.0  0.0 
GAGES 1995 07  04  10  35  19.0  20.0  16.0  8.0 

EVENT "Event of 6 July 1995- new radar data" 
NRPDS  8 
NRGAG  -1 
RADAR 1995 07  06  10  54  kzmx_07_06_10_54.asc
RADAR 1995 07  06  11  00  kzmx_07_06_11_00.asc
RADAR 1995 07  06  11  06  kzmx_07_06_11_06.asc
RADAR 1995 07  06  11  12  kzmx_07_06_11_12.asc
RADAR 1995 07  06  11  18  kzmx_07_06_11_18.asc
RADAR 1995 07  06  11  24  kzmx_07_06_11_24.asc
RADAR 1995 07  06  11  30  kzmx_07_06_11_30.asc
RADAR 1995 07  06  11  36  kzmx_07_06_11_36.asc

6.3 Interpolation Between Gages

The rainfall interpolation technique for spatially varied rainfall is specified with either the RAIN_INV_DIST or RAIN_THIESSEN project cards. No interpolation method can create information without creating uncertainty. All interpolation methods “estimate” the spatially varied field from point measurements, introducing uncertainty. The Thiessen polygon method is simply a nearest-neighbor approach, while the inverse distance squared method produces smooth fields based on the assumption that the influence of a measured value decreases with the distance from the point of measurement squared. The spatial variability of instantaneous rainfall is correlated at a scale of only several km at most. Use caution. It is not appropriate to use rain gage data at distances greater than the correlation length of the rainfall rate field. For example, the use of 15-minute rain gage data from a gage that is 30 km from the catchment is completely unrealistic, making calibration impossible. If you have no rain gages in the catchment separated by less distance than the correlation length of the rainfall field, you do not have enough data to calibrate GSSHA.

For example, Ogden and Julien (1993) calculated the correlation length of radar-estimated rainfall rates from multi-parameter observations of convective rainfall. Putting aside the discussion of the appropriateness of the correlation length as indicator of spatial structure due to the anisotropy and non-stationary of rainfall, the calculated correlation length was on the order of 2.5 km. This means that a rain gage network with inter-gage distances greater than this distance will not capture the true spatial variability of rainfall. There is also a chance that significant rainfall will be completely missed by such a network.

The U.S. National Weather Service network of WSR-88D next generation (NEXRAD) weather radars offer the potential of providing rainfall estimates in locations where there are no rain gages. There are numerous error sources that affect the conversion of radar observations into rainfall rate estimates. Discussions of radar errors are beyond the scope of this manual. NEXRAD precipitation estimates can be used in GSSHA, by formatting the data into a GSSHA precipitation file using the RADAR precipitation type card. When using NEXRAD rainfall estimates, GSSHA assigns a rain gauge at the center of each radar data pixel. When combined with Thiessen polygon rainfall interpolation, this reproduces the original radar pixels. The use of inverse-distance squared interpolation should not be used with radar data.

6.4 Interception

The interception of rainfall by the vegetation is modeled in GSSHA using the two parameter method published by Gray (1970). An initial quantity of rainfall (mm), entirely intercepted by foliage, is specified with the STORAGE_CAPACITY project card. The fraction of rainfall retained as intercepted water after satisfying the storage capacity is specified with the INTERCEPTION_COEFF card. These two cards are used to specify GRASS ASCII maps of the required parameters. Alternatively, both storage capacity and the interception coefficient may be assigned use the Mapping Table and an index map based on vegetation.

In GSSHA, 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 (mm),
• b is the interception coefficient (unitless, between 0.0 and 1.0),
• and I is the cumulative interception depth.

Storage capacity and interception values are usually inferred from vegetation or land cover. Values of storage capacity and interception coefficient values can be found in Gray (1970, section 4.6) or Bras (1990) p. 233.

To reflect the seasonal nature of interception in temperate regions, the storage capacity variable can be seasonally adjusted with the SEASONAL_RS card. This card is specified in the project file. If this card is specified, the storage capacity is varied based on the month. This option is only valid if LONG_TERM simulations are being conducted. If using this option, peak growing season values should be specified for the storage capacity. These values will be modified inversely with the canopy resistance values, as described in Section 9 Continuous Simulations. While this method has not been independently validated, it is based on the reasonable assumption that the storage capacity for interception is dependent on the available leaf area. Canopy resistance is the inverse of the leaf area index, thus storage capacity is inversely related to the canopy resistance. Factors were developed for two regions in the northern hemisphere, basically north and south, north is specified as latitudes above 37.0 degrees. The storage capacity is divided by the following factors for each month. The methods were employed in the calibration and validation of flow and sediment data in the Eau Galle Model (Downer, 2008) Eau Galle TN.

Table - Seasonal storage capacity adjustment factors

Water ponded on overland flow plane cells will infiltrate into the soil as conditions permit. Infiltration is dependent upon soil hydraulic properties and antecedent moisture conditions, which may be affected by previous rainfall, run on, ET, and the location of the water table. In GSSHA, the unsaturated zone that controls infiltration may be simulated with a 1-D formulation of Richards’ equation (RE), which simulates infiltration, ET, and soil moisture movement in an integrated fashion. Infiltration may also be simulated using traditional Hortonian Green and Ampt (GA) (Green and Ampt, 1911) approaches which are simplifications of RE. There are three optional GA based methods to calculate infiltration for Hortonian basins: 1) traditional GA infiltration, 2) multi-layer GA, and 3) Green & Ampt infiltration with redistribution (GAR) (Ogden & Saghafian, 1997). The traditional GA and multi-layer GA approaches are used for single event rainfall when there are no significant periods of rainfall hiatus. The GAR approach is used when there are significant breaks in the rainfall, or for continuous simulations.

RE is a general equation and can be applied in any type of watershed or conditions. However, the simpler methods based on the GA equation are preferred when runoff is Hortonian, i.e. occurs due to infiltration excess, where the rainfall/run-on of water is greater than the possible infiltration rate. For fine textured soils the GAR method has been shown to closely mimic the RE solution (Ogden and Saghafian, 1997) and when applied in basins identified as Hortonian, the GAR method has been shown to produce results comparable with the RE (Downer and Ogden, 2003a).

However, when Hortonian flow is not the predominant stream flow producing mechanism, application of GA type models is ill advised and can result in erroneous results (Downer et al., 2002a). For cases where Hortonian flow is not the predominate process generating stream flow the RE should be used, and coupled with the saturated groundwater solution as appropriate. The saturated groundwater model can also be coupled with the GAR model. This represents a rather crude approximation of the actual processes, but has proven useful for numerous studies Eau Galle TN , JD31 TN. The advantage of this approach is a savings in computation time.

Representation of the soil column below each cell with Richards’ equation is presented. Formulation and application of the GA model is well described in other sources (i.e. Maidment, 1993) as well as the GAR method (Ogden and Saghafian, 1997). Formulation, solution, and application of the multi-layered GA model as applied in the GSSHA model is presented in Section 7.3.

7.1 Richards’ Equation

Detailed modeling of the soil water profile in the vadose, or unsaturated zone, is a key addition in the GSSHA model. It is this dynamic area that controls the flux of water between the surface and groundwater and partitions rainfall into infiltration, runoff, groundwater recharge and ET. The most rigorous way to model these complicated and integrated phenomena is to use Richards’ equation in the solution of the problem. While many have described RE as being infeasible for use in hydrologic predictions the equation has been successfully used in field scale and watershed scale models of soil moisture and runoff (Lappala et al., 1987; Hutson and Waggenet, 1989; Dawes and Hatton, 1993; Refsgaard and Storm, 1995; and others). Simpler methods, such as GA and GAR, which are approximations of RE do not provide detailed soil moisture profiles or simulate the movement of water from the groundwater to the unsaturated zone. Accurate representation of layered soils or soils with a water table difficult is also difficult with these approximations (Short et al., 1995). To simulate infiltration with Richards equation use the INF_RICHARDS option in the project file.

Processes in the unsaturated zone important to surface water hydrology, infiltration, ET, and groundwater recharge, are largely oriented in the vertical direction (Refsgaard and Storm, 1995). GSSHA solves the one-dimensional (vertical direction) head-based formulation of Richards’ equation

(14)

where:

C is the specific moisture capacity,

is the soil capillary head (cm),

z is the vertical coordinate (downward positive) (cm),

t is time (hr),

K () is the effective hydraulic conductivity (cm),

and W is a flux term added for sources and sinks (cm/hr).

The head-based formulation allows solution of Richards’ equation in both saturated and unsaturated conditions (Haverkamp et al. 1977). The head-based approach has been successfully implemented in the field scale unsaturated flow models VS2D (Lappala et al., 1987), LEACHM (Hutson and Wagenet, 1989), and SWAP (van Dam and Feddes, 2000) and in the surface water hydrology models TOPOG_IRM (Dawes and Hatton, 1993) and MIKE SHE (Refsgaard and Storm, 1995). With a head-based formulation, rainfall, evaporation and groundwater recharge can all be modeled without a change in variable. Known mass balance problems associated with solution of the head-based formula (Celia et al. 1987; Ross, 1990; Pan and Wierenga, 1995; Celia et al. 1990) have largely been eliminated by the development of new solution techniques (Kirkland et al., 1992; Rathfelder and Abriola, 1994; Pan and Wierenga, 1995).

In GSSHA, RE is solved using an implicit finite difference approach, which maps directly to the overland flow finite difference grid. The soil column below each overland flow cell is sub-divided into multiple unsaturated zone cells, Figure 9.

Figure 9 - GSSHA representation of the unsaturated zone

A 1-D approximation of the unsaturated zone imposes some limitations on the application of the model. Lateral flow in the unsaturated zone may be significant for a perched water table that does not extend to the soil surface. Lateral flow of unsaturated water may also occur in the absence of a perched water table (Zaslavsky and Sinai, 1981). Lateral flow in the unsaturated zone under either saturated or unsaturated conditions is most significant for steep slopes. Solution of the multi-dimensional RE is difficult and the additional computational burden is not commonly justified if the purpose of the model is to simulate surface water hydrology. Aspect ratio problems are also avoided because of the 1-D formulation employed. Other conditions in the unsaturated zone that are not explicitly simulated are macro-pores and hysteresis.

7.1.1 Discretization

The soil column must be sub-divided into cells to numerically solve Richards’ equation. The discretization used in the GSSHA model is shown in Figure 9. The soil column is defined in three layers, namely A, B, and C horizons. The soil properties and discretization can vary in each layer. Any size cell can be used in the discretization, however, if very large cells are used the solution may be poor. As discussed in Downer (2002a) and Downer and Ogden (2002b), cells on the order of 1 cm are generally needed in the top 10 cm of soil to accurately model the infiltration process. If larger cells are used in dry soils a large volume of water must be added to the cell before significant infiltration can occur because the hydraulic conductivity will be very low. This can result in significant underestimation of infiltration, as discussed in Downer (2002a) and Downer and Ogden (2002b). The grid size for each soil in the SOIL_TYPE_MAP in each soil layer is specified in the SOIL_LAYER_INPUT_FILE or in the Mapping Table.

The solution scheme employed is implicit, central difference in space and forward difference in time, and is thus second order accurate in space, first order accurate in time. For j being the cell numbering scheme, increasing downward, such that j is the cell of interest, j-1 is the cell above, and j+1 is the cell below, and n represents the time level, Figure 9. The following discretization is used for non-boundary cells

(15)

where:

C = water capacity,

= pressure head in the cell (cm),

K = hydraulic conductivity in the cell (cm/hr),

W = source term (cm/hr),

Δt = time step (hr), and

Δz = grid size (cm).

Values for variables represented between cells, or at the cell interface, are represented by j-1/2 and j+1/2 for the interfaces above and below, respectively. Inter-nodal distances, Δz_j-1\2, Δz_j+1\2, are defined as the distance between the center of cells j and j-1 and j and j+1, respectively, computed as:

(16)

Inter-cell hydraulic conductivities may be calculated two ways, an arithmetic weighting of the values or a geometric average. The inter-cell hydraulic conductivity weighting method is selected with the RICHARDS_K_OPTION project card with the argument being either ARITHMETIC or GEOMETRIC; the default is ARITHMETIC. In using the arithmetic weighting, a weighting factor, , is used to determine how much weight is placed on upper and lower cells. The value of is specified with the RICHARDS_WEIGHT project card; the default value is 1.0. If α = 1, then only the value from the j-1 cell is used; this is commonly referred to as backwards difference, or upwinding. If α = 0 only the j+1 value is used. This is commonly referred to as forward difference, or downwinding. If α is 0.5 then equal weight is applied to both the j-1 and j+1 cells. Lappala (1981) recommends that backward difference be used for modeling the infiltration process. The geometric average is best when there are large changes in the hydraulic conductivity between cells and is generally applicable in all situations.

7.1.2 Non-linear Coefficients

In Richards’ equation, hydraulic conductivity and water capacity values are dependent on the water content of the cell. The method to describe these relationships must be specified with the RICHARDS_C_OPTION project card. The possible arguments are BROOKS for the Brooks and Corey method (1964) and HAVERCAMP for the Havercamp method (Havercamp et al., 1977). The Brooks and Corey method (1964) as extended by Hutson and Cass (1987) may be used to estimate relative hydraulic conductivity from the soil pressure head as:

if < _b (17)

and if > _b K_r = 1.0 if

where: K_r is the relative hydraulic conductivity, _b is the air entry or bubbling pressure, and λ is the pore-size distribution index, which is the inverse of the ratio of the length of flow path through the soil matrix to the straight line length.

With the Havercamp method (Havercamp et al., 1977) as modified by Lappala et al. (1987) the relative hydraulic conductivity is calculated as:

(18)

where: A and B are parameters fitted to laboratory determinations of hydraulic conductivity at different soil pressure heads.

The water capacity C is also dependent on the pressure head . Water capacity is defined as the change in moisture with respect to head . For Brooks and Corey if < b the relationship between moisture content and head is expressed as

(19)

where: , θ is the moisture content, θ_r is the residual saturation, and θs is the saturated moisture content, or porosity. Rearranging and taking the derivative with respect to head yields

(20)

In the original formulation of Brooks and Corey, if > b, θ = s and C=0. This formulation permits changes in pressure head without a resulting change in water content for pressures greater than the bubbling pressure. Hutson and Cass (1987) extended the Brooks and Corey equation into the wet region where > b using the following formulation

(21)

where: , and θ_c is the critical moisture content, the moisture content where the curve changes from the standard Brooks and Corey curve to the modified curve of Hutson and Cass. The critical moisture content is defined as with the corresponding critical pressure head:

(22)

The moisture content can then be represented as

(23)

And, the water capacity is found by taking the derivative of the above equation.

(24)

For the Havercamp equations the relationship between the water content and the pressure head is

(25)

where: α and β are fitted parameters. According to Lappala (1981) the form of the equation can be expressed as:

(26)

Differentiating with respect to pressure yields the water capacity

(27)

The relationship between soil moisture and suction head as represented by the different methods is shown in Figure 10.

Figure 10 – Water retention curves (BC – Brooks and Corey).

7.1.3 Evapo-transpiration Source Term

During long-term simulations potential evapo-transpiration (PET) is calculated by either the Penman Monteith (Monteith, 1965; 1981) or Deardorff (1977; 1978) equations and is applied to each soil column below the overland flow plane. Any water ponded on the surface of the overland flow cell is reduced up to amount of the PET. Any remaining PET demand is applied to cells in the unsaturated zone down to the specified root depth (Figure 9). The actual evapo-transpiration (AET) is distributed over the cells in the specified root zone in proportion to the size of each cell. The AET is computed from the PET by adjusting the PET for the soil moisture in each cell. AET depends on the soil moisture, hydraulic properties of the soil and plant characteristics. At water contents at or above the field capacity (θ_fc) there is no stress on plants and AET is equal to PET (Shuttelworth, 1993; Dingman, 1994). The field capacity is the water content at which the suction pressure prevents gravity drainage of the soil. Field capacity is not an actually physically measurable quantity and there are many descriptions of the conditions that correspond to the concept of field capacity. At water contents below the wilting point, θ_wq, plant transpiration ceases and plants wilt (Dingman, 1994), and AET is zero. At intermediate points AET will depend on PET and the water content. Many relationships to relate AET to PET have been suggested (Dyck, 1983). In GSSHA AET is calculated from PET for water contents greater than the wilting point using the relationship

(28)

The appropriate equation depends on vegetation, climate, and soil type. If the exponent, P, is 1 the relationship is linear, above 1 the curve is convex, and below 1 concave. Currently P is set to 1.0. Future versions will allow the power to be specified by the user to reflect the local conditions. The AET for each cell is added to the source term, W, for that cell in the RE solution.

7.1.4 Upper Boundary Condition

The upper boundary condition varies depending on the condition of the top cell: specified flux for no surface ponding, and specified pressure (head) when infiltration excess results in surface ponding. The first cell in the column, j=0 (Figure 9), is located above the ground surface, and a pressure is always specified for this cell. For a flux boundary condition, the pressure in the top cell is zero and the flux is added to cells via the source term, W. For a head boundary the pressure in the top cell is equal to the depth of ponded water. The typical sequence of events is described below.

At the beginning of simulations and at times between rainfall events there is typically no ponded water on the overland flow plane. The upper boundary condition for the unsaturated zone is a negative flux, equal to the AET of the top cell, cell 1 in Figure 9. When rainfall, runoff or some other source of water is added to the overland flow cell a flux boundary condition is initially specified for the RE solution. A flux equal to the depth of ponded water divided by the current time step is added to the source term, W, of the first non-boundary cell, cell 1 in Figure 9. New heads are calculated with this assumed flux. These heads are used to compute the inter-cell fluxes. The computed flux in cell 1 is compared to the source term. If the calculated flux is less than the specified (assumed) flux, then all the water cannot infiltrate into the soil column during the current time step. In this case, water ponds on the soil surface and the upper boundary is changed to a specified head, the head being the depth of ponded water. Heads are re-computed using the head boundary. To save computation time, the upper boundary condition remains a specified head until the overland flow cell becomes dry. At this point the boundary condition changes back to a specified flux, and remains a specified flux until infiltration capacity of the soil column is again exceeded. Water that enters the top cell is infiltration. Any water that does not infiltrate can become runoff or direct ET (DET).

As discussed in Downer (2002a) and Downer and Ogden (2002b), the assignment of the hydraulic conductivity at the ground surface, K_1/2 between j=0 and j=1 (Figure 9), has important implications in determining infiltration, infiltration excess, and ET. Three methods, specified with the RICHARDS_UPPER_OPTION project card and the appropriate argument, can be used to compute the hydraulic conductivity at the ground surface. The default option, NORMAL, is to use the cell-centered value of hydraulic conductivity in cell 1 (Figure 9). The cell-centered value of hydraulic conductivity is always used for a flux boundary condition. If the soil surface boundary condition is a head, hydraulic conductivity at the soil surface may be assumed to be the saturation value, selected by using the GREEN_AMPT argument, or an average of the saturation value and the cell-centered value, selected with the AVERAGE argument. The assumption is that there is always a very thin layer of saturated material at the soil surface any time water ponds. Either method results in increased infiltration compared to using the cell-centered value of hydraulic conductivity. As discussed in Downer (2002a) and Downer and Ogden (2002b), testing at two watersheds indicated either alternative method allows the use of larger cell sizes in the unsaturated zone without seriously affecting calculated hydrologic fluxes.

7.1.5 Lower Boundary Condition

Three different lower boundary conditions can be specified. When the groundwater table is far from the ground surface the lower boundary condition is a zero head gradient. Water entering the N-1 cell in the soil column exits at the incoming rate (Figure 9). This boundary condition is valid when the water table is so deep that its effect on processes in the upper soil column is negligible, and is the default option when a WATER_TABLE is not specified by the user. The lower boundary can also be a fixed water table, specified with the WATER_TABLE project card that specifies the name of a map containing starting groundwater elevations. When a fixed water table is simulated the top of the last cell, N, is just at the surface of the saturated groundwater (Figure 9). The pressure at the top of the cell is zero, and the pressure at cell’s center is positive, calculated as 0.5ΔZ_N (Figure 9). This groundwater boundary fluctuates when saturated groundwater is simulated. For a moving water table, the size of the last non-boundary cell, N-1, and the number of cells, N, changes as the water table rises and falls (Figure 9). A moving water table is specified by using the GW_SIMULATION card in the project file, and then supplying the required saturated groundwater inputs.

The unsaturated zone does not include the saturated zone. The flux between the saturated and unsaturated zones is calculated as part of the overall unsaturated solution. This separation of the saturated and unsaturated zones helps in determining mass balance errors and maintaining mass balance in both the saturated and unsaturated zones because water is either in one zone or the other. When water crosses the boundary between the saturated and unsaturated zones it is removed from one zone and placed in the other. Mass balance errors can occur for a variety of reasons including: model formulation, spatial and temporal discretization used, solution technique, and abrupt changes in material properties. In GSSHA the mass balance for each compartment, overland, stream, unsaturated, saturated, are calculated independently and integrated to compute an overall mass balance. Mass balance errors in any compartment may be controlled by reducing time steps or increasing stability criteria, as needed. The methods used to link the saturated and unsaturated zones are further described in later sections.

7.1.6 Solution

Richards’ equation is highly non-linear because both the water capacity and the hydraulic conductivity depend on the pressure, or water content, of the soil. Numerical solution of RE requires some type of linearization. In GSSHA the RE is linearized by making the water capacity and inter-cell hydraulic conductivity constant during each time step. With flux updating of the heads, as described by Kirkland et al. (1992), this provides a stable, accurate, and mass conserving solution for most conditions.

For N cells including the upper and lower boundary cells, N-2 equations are needed. The well-known Thomas algorithm (Thomas, 1949) is efficient for solution of the resulting tri-diagonal matrix. After solving for heads in each cell, flux updating (Kirkland et al., 1992) synchronizes the heads and soil moistures, and improves the mass balance. Fluxes (cm hr-1) across the top face of each cell, f_j-1/2, are computed as:

(29)

These are used to compute the change in water content of each cell as

(30)

The new water content is

(31)

Once the water content in each cell is updated, pressure heads in each cell are calculated based on the head-moisture content relationship and used in the next time step. Kirkland (1991) found it necessary to restrict flux updating to cells that were both unsaturated and not immediately adjacent to saturated cells. Testing by the authors confirms this needed restriction. The reason for this requirement is that in saturated cells, head changes can occur without a change in water content. This deficiency in the formulation can cause errors in the solution and the mass balance. In this case, iterating on K_r and C can produce substantial improvements in accuracy and mass balance.

The maximum number of iterations is specified by the user using the RICHARDS_ITERMAX project card, the default being one. While iterating on the hydraulic conductivity and water capacity will almost always improve the solution and mass balance to some extent, iterating is particularly useful anytime there are saturated cells in the unsaturated zone. Cells can become saturated when large changes in material properties exist between adjacent layers, when the upper boundary condition is a head, and when the saturated groundwater is rapidly rising. Experience indicates that when the top boundary is the head condition, the accuracy and mass balance are always improved by iterating on the non-linear coefficients. When the top boundary condition is a head the maximum number of iterations changes from the default to five, unless the user has specified more iterations.

Picard iterations (Celia et al., 1990) are used anytime the user specifies iterations or the upper boundary condition is a head. Heads at the n+1 time level are first calculated using values of K_r and C from the n time level. Water contents are calculated based on the fluxes, and heads are updated based on the water contents. K_r and C are then calculated based on the updated values of head and water content. The water contents and heads are set back to the n time level values and K_rⁿ⁺¹ and Cⁿ⁺¹ are used to compute values of head at the n+1 time level. The procedure is repeated until the convergence criterion is met or the maximum number of iterations is reached. The convergence criterion is applied independently to each soil column. The convergence criterion is the difference in head between iterations. The error is adjusted for the head in each cell because the error is most important for wet cells, where a small absolute error in head can cause a large error in the solution. The error for each cell is expressed as

(32)

where k expresses the iteration number, which should not be confused with n, which represents the time level. Iterations for a particular soil column cease once the maximum change in head is less than 1 mm.

7.1.7 Time Step Limitation

Because the discretization is implicit there is no inherent stability criteria for the scheme. However, a time step limitation is desirable to keep the scheme accurate and mass conserving. The time step limitation in GSSHA is adapted from Belmans (1983). The time step is limited such that a maximum change in water content, Δθ_allow, is not exceeded. If the maximum change in water content in any cell, Δθ_max, exceeds Δθ_allow the time step is reduced so that an exceedance during the next time step is not likely. The following limitation is used

(33)

The maximum allowable change in water content is specified by the user, with a suggested range of 0.002 to 0.03 (Belmans, 1983). Smaller limitations will result in longer simulation times. In GSSHA each soil column has its own time step limitation computation. This can greatly increase the speed of the model when rapid changes in water content occur in only a few soil columns. In this case, Δt for these few soil columns may be very small while Δt for the bulk of the soil columns in the watershed is still very large in comparison. The time step is limited by specifying the maximum water content change allowable with the RICHARDS_DTHETA_MAX project card; the default value is 0.025.

7.2 Green and Ampt (GA)

The Green-Ampt (GA) equation, developed for infiltration into uniform soils, has proven very useful in the field of hydrologic engineering. In the GA model water is assumed to enter the soil as a sharp (completely saturated) wetting front. The soil is assumed to homgeneous and infinite, with no restrictions due to the water table. The water can be thought of as entering the soil as a piston, shown in the figure below.

GA Infiltration Model

The GA method is selected by placing the GREEN_AMPT card in the project file. Four soil hydraulic parameters are needed for modeling the infiltration process using the G&A technique:

• GA hydraulic conductivity, K_ga (cm/h), half the saturated hydraulic conductivity K_s
• GA wetting front suction head, Φ_f (cm),
• effective porosity, θ_e,
• and initial soil moisture content, θ_i.

These values may be assigned by using the Mapping Table and an index map or with four GRASS ASCII maps. If using the Mapping Table and an index map the MAPPING_TABLE card is placed in the project file. The MOISTURE card is used with no argument. If using the GRASS ASCII maps to assign values to every cell in the watershed the CONDUCTIVITY, CAPILLARY, POROSITY, and MOISTURE cards are used to specify the GRASS ASCII map names as arguments. Assignment of parameters for all infiltration methods is discussed in Section 7.5.

7.3 Multi-layer Green and Ampt

While the basic Green and Ampt model has been proven effective for modeling infiltration into soils, several important common natural phenomena cause the assumption of a vertically uniform soil to be overly restrictive. Conditions such as layered soils, non-uniform initial soil moistures, surface crust, lenses, and high water tables all violate this basic assumption, yet are routinely encountered in the field. Therefore it becomes necessary to deal with these situations in some manner. The method used can have significant, if not overwhelming, effects on the hydrograph produced. This section discusses a methodology to incorporate a three layered soil profile into the GSSHA hydrologic model using Green-Ampt infiltration. Layering in soil may also be simulated with the Richards’ equation, Section 7.1.

In GSSHA v5 and above, the multi-layer GA model can be used in continuous mode, as well as with groundwater and in simulated constituent transport, including soils.

7.3.1 Layered Soils

In nature, layered soils are the norm rather than the exception. Everyone is familiar with the three layer soil system consisting of a top, A, middle, B, and bottom, C, horizon (Figure 9). Typically, soils in the A horizon are loose and high in organic material. High biological activity increases the porosity and hydraulic conductivity of soils in this layer. Soils in the B horizon are typically less permeable, with lower organic content and reduced biological activity. Soils in the C and lower horizons tend to be even less permeable, with minimal biological activity. What is found in the field could vary greatly from this simple model. Erosion and sedimentation, both recent and in geologic time frames, can cause very different layering effects. This can result in layers that become more porous with depth, or alternating porous and impermeable layers. Also, thin impermeable surface crusts can exist where compaction, raindrop impact or tillage occurs.

Regardless of the layering, when a wetting front moves from one layer to the next, the infiltration rate will be reduced. This is true regardless of the orientation of the layering. If a wetting front moves from a more porous to a less porous layer, the infiltration is reduced due to the reduced hydraulic conductivity. If the wetting front moves from a less porous layer to a more porous layer, the reduced wetting front suction causes the reduction in the infiltration rate. If the soils making up the lower layers are much coarser than the top layers, the wetting front may become unstable, and infiltration will be greatly reduced.

In a single layer infiltration model, the effect of layering must be addressed in some manner. Two possible approaches are: 1) use the properties of the top soil to represent the entire column, or 2) use some type of average values based on the depths of the layers.

7.3.2 Varying Moisture Contents

Soils do not have to be layered to cause problems with the Green-Ampt assumption of uniform soil properties. Soil moisture can also vary with depth. The variation in soil moistures with depth could be caused by layering or due to other reasons, such as interaction with the ground water table. In the latter case, the soil will tend to be driest at the surface and saturated near the water table. During periods of evaporation, the soil moisture content may vary considerably with depth. This type of condition can be roughly approximated with the three soil layer model, where the only parameter that changes through the layers is the initial soil moisture.

7.3.3 Model Formulation

For GA infiltration in GSSHA, infiltration occurs in what is assumed to a uniform soil profile with constant hydraulic parameters and initial soil moisture throughout the profile. The Green-Ampt routine has been expanded to allow for three different soil layers to be modeled. These layers can be of any thickness and hydraulic parameters. The assumption of a sharp wetting front still applies. As the wetting front crosses the layers there is assumed to be an instantaneous change in the initial moisture, porosity, and wetting front suction head. However, the effective hydraulic conductivity is calculated at each time step based on the depth of the leading edge of the wetting front. The hydraulic conductivity, K (cm/hr), is calculated as the harmonic mean of the wetted layers. While the front is contained in the first layer, the hydraulic conductivity at any time n, K_n, is equal to K₁, the hydraulic conductivity of the A soil horizon.

After the wetting front passes into the second layer, the hydraulic conductivity is defined as:

(34)

where:

L₁=thickness of the A horizon (cm),

L₂=thickness of the B horizon (cm),

MD₁=moisture deficit of the A horizon (dimensionless),

MD₂=moisture content of the B horizon (dimensionless),

K₁=saturated hydraulic conductivity of the A horizon (cm/hr),

K₂=saturated hydraulic conductivity of the B horizon (cm/hr), and

Fⁿ= cumulative infiltration at the n time level (cm).

The moisture content MD is defined as the effective porosity θ_s minus the initial moisture content θ_i.

When the wetting front reaches the third layer, the effective hydraulic conductivity becomes:

(35)

where MD₃ (cm), and K₃ (cm/hr) are the moisture content and hydraulic conductivity, respectively, of the third soil layer.

Therefore the change in effective hydraulic conductivity is gradual, as would be expected in nature. The soil is assumed to be saturated at all points behind the wetting front.

The cumulative infiltration at the n+1 time level is (Adapted from Chow et al., 1988)

(36)

where: _ƒ, the wetting front suction head, and MD are for the layer containing the leading edge of the wetting front. This equation is solved iteratively by the Newton Raphson method. Once the cumulative infiltration is determined the infiltration rate, ƒ(cm/hr) is calculated by

(37)

Actually, for this method the calculation of the effective hydraulic conductivity lags the calculation of the infiltration by one time step. However, since the time steps for GSSHA are small, on the order of 1 minute, this should present minimal concern, because the hydraulic conductivity will change very little over such a short period.

In the case where rainfall rate is limiting, the cumulative infiltration is calculated as

(38)

where: dn+1 is the depth of water ponded on top of cell before infiltration is calculated.

The infiltration rate is then

(39)

7.3.4 Inputs

The multi-layered GA input represents an intermediate step in transforming from a map based input system to an index map and table based input system. The multi-layer GA approach is selected by placing the INF_LAYERED_SOIL card in the project file. Multi-layer GA requires the same inputs as the traditional GA approach except that the information must be specified in a table for all three layers. The table is specified with the SOIL_LAYER_INPUT_FILE card. This table is referenced to an index map that contains a unique integer number, starting with 1, for each soil type in the watershed. This GRASS ASCII map is specified with the SOIL_TYPE_MAP project card. The soil layer input file contains the number of soils on the first line, followed by the soil number and parameter values for each soil in the table. Soil numbers should start with soil 1 and increase monotonically. Every soil number in the table should correspond to a number in the SOIL_TYPE_MAP.

In GSSHA version 5.0 and above, the multi-layer GA model has been modified to work in long term simulation mode. The multi-layer GA model can also be coupled to saturated groundwater, and can be using in contaminant modeling, including simulation of contaminants in the soil column. The input table for the multi-layer GA model reflect these changes.

As shown below, values in the table are separated by spaces or tabs, not commas:

Where: Ks is the saturated soil hydraulic conductivity (cm hr^-1),
   λ is the pore distribution index,
   Sf is the moisture front suction head (cm),
   θ is the soil moisture (fraction of total): e - saturated, f - field capacity, wp - wilting point, r - residual, i - initial,
   and d is the layer depth (cm).
The numbered subscripts refer to the soil layer number, 1 – top layer, 2 – middle layer, 3 – bottom layer.

For a watershed with three soil types, soil one a uniform clay, soil two a uniform sand, and soil three a sand with an embedded clay layer, the table might look like:

The corresponding SOIL_TYPE_MAP would contain the integers 1, 2 and 3 corresponding to the locations of the soils 1, 2, and 3 in the SOIL_LAYER_INPUT_FILE.

Newer versions of GSSHA, 5.7 and above, support the multi layer inputs in the MAPPING_TABLE_FILE. For these versions the inputs can be included as stated above, or more effciently, included in the MAPPING_TABLE_FILE. In this case, the above inputs are put into the MULTI_LAYER_SOIL table, as shown below. The inputs and the format are the same as above. WMS v9.0 and above support this ability. The table below was developed with WMS v9.1. It should be noted that the term BUB_PRESS is actually the wetting front suction head Sf.

7.4 Green and Ampt with Redistribution (GAR)

The GAR method expands the capability of the GA method by redistributing soil moisture during periods of no- or low- intensity rainfall. This allows infiltration capacity to recover for the next burst of storm intensity, and makes the GAR method suitable for simulating multiple rainfall events in series. When using the GAR method, multiple wetting fronts occur as burst of rainfall occur. During rainfall haitus these fronts are redistributed to produce a new value of initial moisture for calculations during the next burst of rainfall, Fig 12.

This method is selected with the INF_REDIST card in the project file. The technique for hiatus and post-hiatus infiltration was developed by Ogden and Saghafian (1997), and is similar to the method developed by Smith et al. (1993).

The GAR technique requires three additional inputs to the four required for the GA method, residual saturation and pore-size distribution index (Brooks & Corey, 1964), and field capacity. These must be input with in the Mapping Table. Assignment of parameters is discussed in Section 7.5.

Ogden and Saghafian (1997) performed a comparison of the GAR approach against a numerical solution of RE on each of the 11 USDA soil textural classifications for two pulses of rainfall separated by a period of no rainfall. It is important to note that the traditional GA is invalid for this case, unless the period of no rainfall is long enough in duration for the soil moisture profile to equilibrate. These results and other analysis presented by Ogden and Saghafian (1997) indicate that the GAR approach is a reasonably good approximation to the numerical solution of RE for well-drained soils, subjected to multiple pulses of rainfall. The agreement between the two approaches is better for finer textured soils. Coarse textured soils dry faster at the surface because of moisture profile “bulging”. The GAR approach cannot predict this phenomena, and therefore, overpredicts near soil-surface water contents after rainfall.

Figure 12. GAR multiple wetting fronts in soil profile.

7.5 Parameter Estimates

It is best to use soil property and Green-Ampt infiltration parameters derived from field and laboratory measurements of infiltration on the study watershed. Even under controlled conditions hydraulic soil property measurements are very difficult. Hysteresis effects and the extremely non-linear behavior of soil water retention make it very difficult to uniquely identify soil infiltration parameters. Hydrologic studies seldom have budgets sufficient to determine the needed parameters in the field.

Considerable prior research has been performed to relate soil infiltration parameter values to textural classification. Some highly relevant references are Rawls and Brakensiek (1983) and (1985), and Rawls et al. (1982) and (1983). Table 9 summarizes Rawls and Brakensiek soil parameter estimates as a function of United States Department of Agriculture (USDA) textural classification. It is important to note that the values listed in Table 9 were derived from the geometric mean of tests on a large number of soil samples. Hydraulic conductivies for all GA based approaches are half of the saturated values listed in Table 9 (Rawls, et al., 1982). The variance of these values is large, indicating significant uncertainty or low correlation between textural classification and soil texture. However, these values are useful because they provide an initial estimate of infiltration parameters. The variances of the values in Table 9 are listed in the original papers, and are published in Maidment (1993).

USDA Textural Classification	Total Porosity/Saturation   θs   (cm3/cm3)	Effective Porosity/Saturation   θe   (cm3/cm3)	Field Capacity Saturation   θf   (cm3/cm3)	Wilting Point Saturation   θwp   (cm3/cm3)	Residual Saturation   θr   (cm3/cm3)	Bubbling Pressure Geometric Mean   b (cm)	Pore Size Distribution Arithmetic Mean   λ   (cm/cm)	Saturated Hydraulic Conductivity (multiply by 0.5 for GA methods)   Ks (cm/h)	Wetting Front Suction Head (Capillary Head)   ƒ (cm)
Sand	0.437	0.417	0.091	0.033	0.02	7.26	0.694	23.56	4.95
Loamy sand	0.437	0.401	0.125	0.055	0.035	8.69	0.553	5.98	6.13
Sandy loam	0.453	0.412	0.207	0.095	0.041	14.66	0.378	2.18	11.01
Loam	0.463	0.434	0.27	0.117	0.027	11.15	0.252	1.32	8.89
Silt loam	0.501	0.486	0.33	0.133	0.015	20.79	0.234	0.68	16.68
Sandy clay loam	0.398	0.330	0.255	0.148	0.068	28.08	0.319	0.30	21.85
Clay loam	0.464	0.390	0.318	0.197	0.075	25.89	0.242	0.20	20.88
Silty clay loam	0.471	0.432	0.366	0.208	0.040	32.56	0.177	0.20	27.30
Sandy clay	0.430	0.321	0.339	0.239	0.109	29.17	0.223	0.12	23.90
Silty clay	0.479	0.423	0.387	0.250	0.056	34.19	0.150	0.10	29.22
Clay	0.475	0.385	0.396	0.272	0.090	37.30	0.165	0.06	31.63

Table 9 - Rawls & Brakensiek soil parameter estimates.

In the table soil moistures θ, are listed for saturation (s), effective saturation (e), field capacity (f), wilting point (wp), and residual (r). These values are applicable for all approaches. These are followed by the bubbling pressure _b (used for RE), the pore distribution index (used for RE, GAR and multi-layer GA in continuous mode), saturated hydraulic conductivity K_s (used directly for RE, halved for all GA approaches), the wetting front suction head f (used for all GA approaches).

Standard practice in developing GSSHA models is to obtain digital soil textural classification data and use these data to develop an index map of soil types. Soil textural maps may be combined with land use or vegetation maps. Land use and vegetation can strongly influence soil hydraulic properties. The Mapping Table is used to assign initial parameters to the soil types in the index map. One or more of these parameters, typically K_s and f or _b, are used as calibration parameters. As discussed by Senarath et al (2000), calibration is best done using an automated calibration method, such as SCE (Duan et al, 1992), combined with long term simulations. The possible parameter values are bounded by the range found in literature values, unless other factors, such as land use or vegetation, dictate otherwise. The range of values may be narrowed by making field and laboratory measurements of parameters.

Impervious areas can be accounted for by modifying the infiltration parameters as described above, or the area available for infiltration into a cell can be adjusted by specifying the impervious fraction using the AREA_REDUCTION process and including an index map for impervious area. Include the AREA_REDUCTION card in the project file and an AREA_REDUCTION mapping table with the fraction impervious for each value in the associated index map. Also see the section on Mapping Tables.

8.1 General

Simulation of saturated 2-D groundwater flow is selected by placing the GW_SIMULATION card in the project file. When groundwater simulations are performed, the starting water surface elevations (m) must be provided in a file specified with the WATER_TABLE card; the file containing bedrock elevations (m) is specified with the AQUIFER_BOTTOM card. The groundwater time step must also be specified (s) using the GW_TIMESTEP card, typically 1 hour, 3600 s. A variety of boundary conditions can be specified with the GW_BOUNDFILE card. Saturated groundwater may be simulated using either the GAR or RE methods to calculate infiltration. In GSSHA v5 and above, multi-layer GA can also be used to simulate infiltration for groundwater models. The saturated groundwater solution was developed in conjunction with the RE unsaturated groundwater solution, and is intended to be used with RE. The GAR method was linked to the saturated groundwater model to provide a simple, quick, and approximate groundwater recharge estimate. This method is a useful screening tool, and can be used to help define parameter values before applying the more rigorous RE method. The multi-layer GA model was extended for use with continuous simulations and as a result, was also extended for use with the saturated groundwater model. Whether GAR or RE or multi-layer GA is used affects how parameter values are assigned, as described below.

8.2 Formulation

Trescott and Larson (1977) described the solution to the two-dimensional free surface groundwater problem, and the efficiency of various solvers. Their methods were largely followed in the development of this portion of the code; an exhaustive coverage need not be presented here. The overall approach, differences in approach, and integration into the GSSHA model are presented.

The controlling equation, as developed by Pinder and Bredehoeft (1968), is:

(40)

where:

T is the transmissivity (m²/s),

h is the hydraulic head (m),

S is the storage term (dimensionless),

and W is the flux term for sources and sinks (m/s).

It is assumed that off diagonal terms are not important and that transmissivity can be expressed as the product of the saturated hydraulic conductivity of the media (K) and the depth of the saturated media (b). For the free surface problem the head is the surface water elevation (E_ws).

(41)

This equation can be represented using a block-centered finite difference five-point implicit scheme as:

(42)

This representation varies from the original representation of Trescott and Larson (1977) in that the transmissivities and the storage terms are both time dependent and calculated implicitly using Picard iteration.

8.3 Solution

The equation is solved by successive overrelaxation by lines (LSOR) (for example Tannehill et al., 1997). LSOR was shown by Trescott and Larson (1977) to be capable of solving a variety of difficult groundwater problems, though not necessarily being the most efficient method. With LSOR, the two-dimensional problem is linearized by solving by rows or by columns. The user specifies solution by rows or by columns with the GW_LSOR_DIR project card, and the choice is made to align the direction of solution with the principal direction of flow, x (argument - 1) or y (argument - 2) (Trescott and Larsen, 1977). LSOR is an iterative method; the user selects the groundwater convergence criteria (m) with the GW_LSOR_CON project card. A typical value, the default, is 10^-5 m. The user also selects an over relaxation coefficient, ω, with the GW_RELAX_COEF card. During each iteration the head in each cell is adjusted by the over relaxation coefficient, ω, such that:

(43)

where k denotes the iteration number. For ω greater than 1.0, the next head is projected out on a line determined from previous iterations. This speeds the convergence process but may also reduce stability. Typically, a value of ω of about 1.2 results the fastest solution. For very difficult problems ω may need to smaller than 1.0, and the solution is said to be underrelaxed.

For each iteration the transmissivities in both the x and y directions are calculated based on the updated saturated depth, b, determined from the heads, so that:

(44)

where K_gw is the lateral hydrualic conductivity of the porous media. This is a variation from Trescott and Larson (1977) who calculate the transmissivity based on the value of b from the last, n^th, time step such that Tⁿ⁺¹ = K_gwbⁿ. The storage term, S, in the equation is used to represent the change in volume with respect to the change in head:

(45)

where V is the volume. In two-dimensional applications it is common practice to represent the storage as the porosity of the saturated groundwater media. When RE is used to define the unsaturated zone the storage term, S, is not a constant but also depends on the head. The storage term is updated each iteration. The storage term is set to the porosity of the unsaturated cell above or below the groundwater surface elevation, depending on whether the groundwater is rising or falling. For time steps when the groundwater moves more than one unsaturated zone cell, the storage term is calculated as a bulk storage term:

(46)

where ΣΔz is the distance that the groundwater surface is anticipated to move during the n+1 time step or k+1 iteration. The anticipated change in the groundwater surface elevation is determined from the source term, W, containing all fluxes in the cell and the storage capacity, S, of surrounding cells. An internal time step limitation attempts to limit the movement of the groundwater surface to a single unsaturated zone cell. This added step helps maintains overall mass conservation. The actual storage available is Θ_s--Θ, but because water content is a discrete value in each cell, calculating storage as Θ_s-Θ over ΣΔz can result in a lack of convergence of the groundwater solution, which is implicitly calculating h, S, b, and T simultaneously. As discussed below in the section on coupling of processes, the difference in water volume between Θ_s and Θ_s-Θ is added to the source term during the next groundwater update, such that mass is conserved. This may result in a small time lag in the groundwater response.

When GAR is used to provide recharge estimates to the groundwater model, the storage term is constant, equal to the value moisture deficit of the soil, effective porosity minus soil moisture. The moisture deficit is updated every rainfall event.

8.4 Assignment of Parameter Values

The value of K_gw (cm/hr) may be specified as a constant value with the GW_UNIF_HYCOND or a map of spatially distributed values specified with the GW_HYCOND_MAP card. The porosity of the groundwater media below the specified unsaturated zone can be a uniform value, specified with the GW_UNIF_POROSITY card, or can be distributed by specifying a map of spatially distributed values using the GW_POROSITY_MAP. If neither uniform or distributed values of these parameters are specified with the above cards, the values will be set to the those of the last soil layer of the soils in the unsaturated zone specified in the SOIL_TABLE_INPUT_FILE or Mapping Table files. When GAR is used to provide estimates of recharge, only hydrualic conductivity need be specified, as the storage is computed using the GAR infiltration parameters.

8.5 Boundary Conditions

Each cell in the active domain is assigned a boundary condition with a map of integer values specified with the GW_BOUNDFILE project card. Boundary conditions around the border of the watershed may be constant head (type 2), constant flux (type 0) or a combination of both. Flux boundaries are represented by setting the transmissivity of flow boundary cells to zero and entering the flux in the source/sink term, W. For all non-head boundary cells (type 1), the unsaturated zone provides a flux to the source/sink term, W. The lower boundary is a flux term (cm/hr) that represents a leaky aquifer, specified with the GW_LEAKAGE_RATE project card. The default is zero. This single user-specified value is used in all cells.

Possible internal boundary conditions include constant head (type 2), constant flux (type 6), variable flux (type 3), river head (type 5) and river flux (type 4), and reservoir cell (type 7). The constant flux and variable flux boundaries are used to represent pumping wells of constant and variable rate, respectively. Water added/subtracted by pumping is added to the source/sink term, W. The possible boundary types are:

0 – no flow

1 – regular infiltration cell (no special boundary condition)

2 – specified head, taken from WATER_TABLE file

3 – dynamic flux (well)

4 – stream cell with calculated flux between stream node and groundwater cell

5 – stream cell with a specified head, value from stream routing solution

6 – static flux (well) value taken from the GW_FLUX_BOUNDTABLE

7 - reservoir cell

8.5.1 Stream Boundary Cells – Type 4 & 5

When 1-D channel flow is simulated, grid cells containing stream channels may be represented as either head (type 5) or flux boundaries (type 4). For a head boundary, the elevation of the water surface in the stream node is used as a specified head in the groundwater solution. For a river flux boundary condition the flux between the groundwater cell and stream node is calculated during the channel routing update. Fluxes accumulated over multiple channel routing updates are accumulated and then added to the groundwater cell in the source/sink term, W. Additional inputs are required to simulate groundwater/stream interactions:

• thickness of the riverbed material, M_rb, (cm), and
• vertical hydraulic conductivity of riverbed material, K_rb (cm/hr).

These may be specified as uniform values using the M_RIVER and K_RIVER project cards, respectively. These may also be distributed along the stream network in the channel input file (cip file), specified with the CHAN_INPUT card, as described in Section 5.1 Channel Routing.

The flux throught the river bed is calculated based on the difference in elevation between the groundwater surface and the water surface elevation in the stream node. McDonald and Harbaugh (1988) developed a simple method to compute flows between the channel network and saturated groundwater based on the Darcy equation. If the groundwater surface elevation is above the riverbed elevation but below the river water surface elevation, the river discharges to the groundwater:

(47)

where:

f = per unit area discharge (m/s),

K_rb = hydraulic conductivity of the river bed material (cm/hr),

M_rb = depth of the river bed material (cm),

E_r = elevation of the river water surface (m), and

E_ws = elevation of the groundwater water surface (m).

The negative flux indicates that the discharge is into the groundwater. If the groundwater surface elevation is below the riverbed elevation the river leaks to the groundwater at the rate

When the STREAM_LOSS card is used in the absence of a known or calculated water table location, Equation 48 is used to compute the stream loss.

If the groundwater surface elevation is above the river water surface elevation, the groundwater discharges to the river according to Equation 46. In this case, the flux will be positive, indicating flow is from the groundwater to the river. The unit flux, ƒ, is multiplied by the top width and length of the channel segment. In calculating the top width an effective depth is defined. If the water surface elevation in the channel is higher than the water surface elevation of the groundwater, the channel depth is used in the calculation. When the groundwater level is higher than the water level in the stream, the effective depth used is the groundwater elevation minus the bed elevation.

The channel side slope properties are assumed to be the same as those of the channel bed. If properties vary significantly between the bed and the side slope an average or weighted average values should be input. This approach is highly empirical and the riverbed property values will require calibration to the system being modeled.

8.5.2 Static Pumping Wells – Type 6

Static pumping rate wells may be placed in any cell in the watershed domain. Pumping wells are located by placing a value of 6 in the GW_FLUX_BOUNDTABLE at each desired well location. The GW_FLUX_BOUNDTABLE specifies the name of a file with the pumping rates. The GW_FLUX_BOUNDTABLE has the following format.

I location well 1          J location well 1          Pumping rate (m³/d) well 1
I location well 2          J location well 2          Pumping rate (m³/d) well 2
I location well 3          J location well 3          Pumping rate (m³/d) well 3
…
…
…
I location well N-1    J location well N-1     Pumping rate (m³/d) well N-1 I location well N        J location well N         Pumping rate (m³/d) well N

Values are separated by spaces or tabs, not commas. Groundwater extractions have a positive pumping rate; injections have a negative pumping rate. The number and I, J locations of the wells in table must match the number and location of wells in the GW_BOUNFILE file.

8.5.3 Reservoir Cells - Type 7

Any cell that is within an active reservoir is assigned a boundary condition of 7. As the reservoirs change size within the GSSHA model, type 7 boundary conditions are assigned internally to GSSHA. If a cell has a type 7 boundary a flux between the reservoir and the groundwater will be calculated as the Darcy flux based on the head between the reservoir and the groundwater, the depth of sediments in the bottom of the reservoir, and the hydraulic conductivity of the sediments. The equations are the same as those described for groundwater channel interaction. Values of hydraulic conductivity for each cell are taken as the hydraulic conductivities prescribed for the infiltration model. The thickness of the bed sediments is specified with the MLAKE card in the CHAN_INPUT file. A uniform value, for all lakes, can be specified in the project file with the M_LAKE card. The value is specified in cm. There is no default value. A different value can be entered for any or all of the lakes. Groundwater interaction will not be calculated for reservoir links without the MLAKE card in the CHAN_INPUT file or the M_LAKE card in the project file.

8.5.4 Static and Dynamic Wells using the Mapping Table

In addition to specifying static wells above, both static and dynamic wells can be created by using the WELL_TABLE table of the mapping table file. An index map of wells with unique pumping characteristics is created (0's elsewhere.) This index map and the associated table is then used to modify the groundwater boundary map. This table has two parameters. The first is a flag (1 for yes and 0 for no) to indicate whether the well is dynamic (1) or static (0). If the well is static, the second parameter is the pumping rate (positive values) or injection rate (negative values). If the well is dynamic then the second value refers to a time series index; you will need to also set up the time series index table.

To create a dynamic well, there are four steps that need to be followed.
1) set up a time series of the pumping rates (positive values) or injection rates (negative values) (see Time series file format)
2) reference the time series in the project file (TIME_SERIES_FILE "filename.ts", as indicated in the above link to time series) preferably before the mapping table.
3) set up the time series index mapping table
4) set up the WELL_TABLE table, with the is_dynamic value set to 1 and the value set to the time series index for that particular time series.

8.6 Coupling of the Saturated Zone Model with the Richards’ Equation Model of the Unsaturated Zone

Groundwater:Coupling of the Saturated Zone Model with the Richards’ Equation Model of the Unsaturated Zone

8.7 Coupling of the Saturated Zone Model to the GAR Infiltration Model

The GAR routine described in Section 7.4 calculates infiltration and soil moisture based on the movement of multiple sharp wetting fronts that occur with rainfall burst followed by rainfall hiatus. The multi-layer GA model described in Section 7.5 calculates infiltration and soil moisture in a user specified distributed three layer soil model. The infiltration calculated by GAR or multi-layer GA can be used to provide a rough estimate of groundwater recharge to the 2-D saturated groundwater flow simulation. The groundwater recharge for each time step is equal to the infiltration computed by GAR or multi-layer GA. In the model the groundwater recharge is computed each time step as

(52)

where: R is equal to the recharge m³, F is the total infiltration (m), and A is the area of the overland flow cell (m²). The simple two layer soil moisture accounting routine is used to calculate soil moisture. When the groundwater elevations exceed the ground surface elevation, infiltration calculations for the cell cease, and the groundwater surface exchange is calculated as described above. Anytime exfiltration occurs, the infiltration and overland flow processes are initiated if they are not already active. These processes remain active as long as exfiltration occurs and until all water on the overland flow plane stops flowing and infiltrating. Infiltration is not calculated for cells in which exfiltration is occurring.

When using the GA models to couple the surface soils to the groundwater table the properties of the unsaturated groundwater media existing between the specified soil column, SOIL_DEPTH for GAR, and water table are assigned the groundwater porosity. The soil saturation (fraction) for this unsaturated soil may be specified by the user with the SINGLE_UNSAT_SAT card. The default value is 0.75, meaning the soil is 75% saturated. This corresponds, roughly, to field capacity, the value at which the soil can no longer gravity drain. For comparison, the equivalent value in MODFLOW, would be 0.0. To obtain a MODFLOW like simulation the user would specify:

SINGLE_UNSAT_SAT 0.0

When using the RE solution to couple to saturated groundwater, soil moistures are calculated along the entire soil column to the free groundwater surface, regardless of the depth of the groundwater. This provides a variable soil moisture along the entire soil column and represents an advantage of the RE method.

More detail on the coupling of GAR and the saturated groundwater model can be found in these documents:

JD31 TN
Soil Moisture Accounting

GSSHA can be used to perform continuous simulations by specifying the LONG_TERM card in the project file, and providing the necessary input data. Continuous simulations require that one of two optional methods for calculating evapo-transpiration (ET) be selected, Deardorff (1978) or Penman-Monteith (Monteith, 1965, 1981), as described in Section 9.1. Hourly hydrometeorological observations are required for either method. The file containing the hourly values of hydrometeorological values is specified with the HMET_TYPE card. The TYPE can be WES, SAMPSON, or SURFAWAYS, as described in Section 9.3. The hydrometeorological inputs are used to compute potential evapo-transpiration (PET) with one of two ET methods. This PET is then sent to the chosen soil moisture accounting method as a demand. The amount of the demand that can be satisfied is the actual ET (AET). The AET is met based on the soil hydraulic properties, soil mositure, balance with other sources (infiltration, flow from groundwater) and demands (drainage). Anytime long-term simulations are performed snowfall accumulation and melting are also calculated, Section 9.2. To perform long term simulations, infiltration must be calculated with either GAR or RE by use of the INF_REDIST or INF_RICHARDS cards, respectively. In GSSHA v5 and above the multi-layer GA model can also be used for continuous simulations. Multi-layer GA is specified with the INF_LAYERED_SOIL card. When the GAR or multi-layer GA method of infiltration is specified, soil moistures are calculated with a two layer soil moisture model, as described in section 9.2. When RE is chosen, soil moisture is calculated as part of the overall solution.

9.1 Computation of Evaporation and Evapo-transpiration

The evaporation and evapo-transpiration models incorporated in GSSHA allow calculation of the loss of soil water to the atmosphere, improving the determination of soil moistures. Two different evapo-transpiration options are included:

• bare-ground evaporation from the land-surface using the formulation suggested by Deardorff (1978), and
• evapo-transpiration from a vegetated land-surface utilizing the Penman-Monteith equation (Monteith, 1965, 1981).

Variants of these two representations are widely used in land-surface schemes of climate and distributed hydrologic models (e.g., Dickinson et al., 1986; Beven, 1979).

To accurately compute fluxes of soil water to the atmosphere, energy fluxes between the atmosphere and the ground must first be computed. These energy, or radiation, fluxes, discussed in Section 9.1.1, are the forcing terms in the evapo-transpiration calculations. Ground temperature is an important component in both the Deardorff and Penman-Monteith formulations, and the fluxes computed in Section 9.1.1 are used in the computation of the ground temperature as described in Section 9.1.2. A detailed presentation of the Deardorff method in Section 9.1.3 describes how the computed fluxes and ground temperatures can be related to bare-soil evaporation. Section 9.1.4 details how the computed fluxes are adjusted and applied for computation of evapo-transpiration with the Penman-Monteith equation. Finally, Section 9.1.5 describes the additional inputs needed for ET calculations.

9.1.1 Computation of Auxiliary Energy and Heat Fluxes

Realistic estimates of incoming and outgoing radiation fluxes are needed for accurate estimates of ET. The important components of the energy budget are long-wave radiation, discussed in Section 9.1.1.1, short wave radiation, discussed in Section 9.1.1.2, and heat conduction into the soil, Section 9.1.1.3. Influences of soil, water, vegetation and cloud albedos, atmospheric emissivities, and sloping-terrain effects on energy fluxes must be included in the calculations. The effects of albedos and emissivities are of great importance in determining radiation fluxes and warrant detailed representation (Pielke, 1984).

9.1.1.1 Net Incoming Short-wave Radiation

The net incoming short-wave radiation can be represented as:

(53)

where: A is the albedo, and R_{s, direct} and R_{s, diffuse} are direct and diffuse contributions of short-wave radiation, respectively.

It is also imperative that short-wave energy influxes be modified to account for sloping terrain effects. As described by Young (1972), even in the Great Plains region of North America (one of the flattest places on earth), only 7% of the land area can be classified as flat in relation to solar radiation calculations. Despite the significance of terrain, most hydrologic and Global Climate Model (GCM) land-surface schemes ignore the sloping terrain effects when estimating the soil energy budgets. Pielke (1984) found substantial differences between the solar radiation values from north and south facing slopes. Pielke and Mehring (1977) found that “the eastern slope of a 1-km mountain (with a slope of about 2°) to be about 1° to 2°C warmer in the morning and cooler by the same amount in the afternoon than the same location in the western slope.” The direct incoming solar irradiance can be represented as:

(54)

where: λ is angle of incoming radiation, ζ is the local land surface slope in the azimuthal direction of the sun, and R_s, horiz,direct,ground is the direct radiation on a horizontal ground surface, computed as:

(55)

where: the transmission coefficient, K_t (dimensionless) is a function of density, type and condition of vegetation, N is the fraction of sky covered by cloud cover (0-1), n is a turbidity factor of air that varies from 2.0 (clear air) to 5.0 (smoggy urban areas), currently fixed at a value of 2.0 in the code. K is the fraction of cloudless sky insolation received on a day with overcast skies, and is given as:

(56)

where z is the height to the cloud-base (km), fixed in the current formulation at 1.5 km. The molecular scattering coefficient, a, is defined as:

(57)

The optical air mass, m, is calculated as:

(58)

where: λ is the angle from the observer’s horizon to the center of the solar disk, and is often referred to as the “solar elevation angle”. Calculation of λ is discussed below.

If no measurements of R_{s,horz,direct} are available, then they may be estimated based on the time of day and year, and the location of the watershed from:

		, if λ < 90°, and	(59a)
			(59b)

where the solar constant, S_o, has a value of 1376 W m^-2 (Hickey et al., 1980). Following Paltridge and Platt (1976) the ratio of the actual earth-sun distance squared (a²) to the average earth-sun distance squared (r²) on any given Julian day m* can be estimated using the following relationship (60):

where:

(61)

where: m* is the number of the day [i.e., ranging from 0 (January 1) to 364 (December 31)], and d_o is the Julian day fraction of the year converted to radians. The angle (i) between the direct solar radiation and the normal to the slope is defined as (Kondrat’yev, 1969) :

(62)

The slope of the terrain (α) is given by:

(63)

where: ∂_zG/∂_x and ∂_zG/∂_y represent the incremental slopes in the x and y directions. The zenith angle (λ) is defined by:

(64)

where: angle φ is the latitude at the site of interest. The orientation of the sun’s azimuth β relative to the azimuth of the terrain slope γ is given by β–γ. The trigonometric relationship for γ is:

(65)

The trigonometric relationship for β is:

(66)

where: δ_sun is the declination angle of the sun (varying between +23.5? on June 21 to -23.5? on December 22), and can be obtained by using the following formulation (Paltridge and Platt, 1976):

(67)

where: δ_sun is in radians. The hour angle, h_r (0°≡noon) is represented mathematically when the sun is east of the observer’s longitude as follows (Curtis and Eagleson, 1982):

(68)

where: T_s is the standard time at the site of interest (counted from midnight; i.e. from 0.00 to 23.59), ΔT₂ is the difference between true solar time and mean solar time in hours (small; hence is neglected in this analysis), and ΔT₁ is the difference between the standard and local longitudes (in hours) given as:

(69)

where i^* = 1 for longitudes located to the east of Greenwich and i^* = -1 for those located to the west, and θ_S and θ_L are the longitudes of the standard- and the observer-meridian, respectively. The standard-meridian is defined as the meridian where the observer’s time zone is centered. When the sun is west of the observer’s longitude the following relationship is used:

(70)

The diffuse short-wave radiation is obtained by using the following relationship (Kondrat’yev, 1965):

(71)

Lee (1978) noted that the difference between R_s,diffuse and R_{s,horiz.,diffuse} is only 2% for slopes less than 30%. Lee (1978) also showed that for slopes less than 36%, the contribution to the total solar radiation by the reflection of total solar from surrounding sloped terrain is only 3% or less. The diffuse horizontal radiation, R_{s,horiz.,diffuse} data are extracted from NOAA-NREL CD-ROMs, and hence atmospheric scattering and absorption effects are neglected.

9.1.1.2 Net Incoming Long-wave Radiation

The net incoming long-wave radiation can be represented as follows:

(72)

where: A is the long-wave surface albedo, σ is the Stephan-Boltzmann constant with a value of 5.67 x 10-8 Jm^-2s^-1K^-4, T_g and T_a are the absolute ground and the air temperatures (Kelvins), respectively, K^* is the cloudiness-correction factor (dimensionless), and E_a and E_g are the emissivities of the air and the ground surfaces, respectively (dimensionless). Following Bras (1990), it is assumed that the ground surface acts as a blackbody, with an emissivity of 1.0. The long-wave radiation from clear skies is strongly related to the atmospheric water content. The following formulation is used to evaluate the atmospheric emissivity (Idso, 1981):

(73)

where e is the vapor pressure (mb), and can be obtained from the following relationship:

(74)

where: rh is the relative humidity (dimensionless, expressed as a decimal number between 0 and 1), and e_s is the saturated vapor pressure (mb). Using the Teten’s formula (Teten, 1930) the saturated vapor pressure can be approximately evaluated as:

(75)

The presence of clouds gives rise to an increase in long-wave radiation. The following relationship, suggested by the TVA (1972), is used to evaluate K^*:

(76)

where N is the fraction of sky covered by clouds.

9.1.1.3 Heat Conduction into the Soil

Many models neglect the soil heat flux (Manabe et al., 1974; Gates, 1974), however, Deardorff (1978) has found the assumption of an insulated soil surface to be especially poor under random atmospheric forcing-conditions. Following Kasahara and Washington (1971), in GSSHA the soil heat flux at the surface (positive when directed into the soil) is represented as a function of sensible heat flux. Deardorff (1978) found the following relationship to be of “intermediate but surprisingly acceptable accuracy” in determining the sum of energy fluxes into the soil (W m^-2) for short time-steps:

(77)

where H_s is the sensible heat flux (W/m²) (positive when directed upward), represented mathematically as:

(78)

where: c_p is the specific heat at constant pressure, equal to 1.013 kJ kg^-1 °C^-1, u_a is the wind speed (m/s) at the reference level, z, and c_H is the dimensionless heat or moisture transfer coefficient applicable to bare soil. Following Deardorff (1978) a value of 0.0025 is selected for c_H. The air density ρ_a is calculated as (kg m³):

(79)

where: P_a is the atmospheric pressure in kPa, and T_a is the air temperature in °C.

9.1.2 Estimation of Ground Temperature

Because ground temperature appears explicitly in the outgoing long-wave energy flux-term, and implicitly in the Deardorff (1978) and Penman-Monteith evaporation formulae, it is important for radiation and evaporation calculations. Ground temperature varies considerably during the diurnal cycle and is not generally measured or provided by meteorological or weather stations. In GSSHA the ground temperature is obtained numerically by solving the surface energy balance equation with the Newton-Raphson iterative method as described by Deardorff (1978).

Assuming that the energy lost due to temporary storage, advection and biochemical usage are negligible for non-vegetated surfaces, the surface energy balance formulation can be written as:

(80)

where: the terms H_s, R_l, R_s and G have been previously defined, and E is the evaporation latent heat flux (W/m²) from bare soil, as described in Section 9.1.3.

The saturation specific humidity (kg/kg), a function of ground temperature, q_sat (T_g) (mb) can be calculated from:

(81)

where: ε is the ratio of molecular weight of water vapor to that for dry air (0.622), and q_a is the specific humidity, estimated by substituting e = rhe_s for e_s, where r_h is the relative humidity.

The Clausius-Clapeyron equation is used to calculate the derivative of q_s with respect to T_g as:

(82)

where: R_water is the gas constant for water vapor (461 J kg^-1 °K^-1) and L is the latent heat of evaporation (2.50036 MJ kg^-1). Following Williamson et al. (1987), the latent heat of evaporation, known to vary with ground temperature, is a constant. The ground temperature can be obtained iteratively using the following relationship:

(83)

where: K indicates the iteration count, ƒ(T_g) is obtained from Equation (80) and ƒ’ (T_g) is obtained from Equation (82). The initial value for the iterative procedure is the surface temperature from the previous time step, T_g^n-1. The procedure is repeated until the following condition is satisfied:

(84)

A convergence criterion, δ_e of 0.001 °K is utilized in this analysis. When the convergence criterion is satisfied the Clausius-Clapeyron equation is solved, and T_gⁿ = T_g^k+1.

9.1.3 Evapo-transpiration

Evaporation from a bare surface can be computed using the Deardorff equation by specifying the ET_CALC_DEARDORFF card in the project file. When using the Deardorff method only the physical components of evaporation are considered. Plants are not considered. Compared to the Penman-Montieth method (1971), use of the Deardorff method can result in higher evaporation and lower soil moistures because residual soil moisture is the lower limit on soil moisture and resistance of the plant canopy to evaporation is not considered. The Deardorff method is not appropriate in vegetated areas and have limited applicability. Inputs to the mapping table for ET calculations, described below, are the same for the two methods.

The ET_CALC_PENMAN project card is used to select the Penman-Monteith method for evapo-transpiration. The Penman-Monteith equation is one of the most advanced resistance-based models available for the prediction of evapo-transpiration from a vegetated land-surface (Shuttleworth, 1993). Although Monteith (1965) points out several simplifying assumptions utilized to derive the Penman-Monteith evaporation formula, Lemeur and Zhang (1990) have found that for arid watershed the performance of the Penman-Monteith is better than both the CRAE (Morton, 1983) and the Advection-Aridity (Brutsaert and Stricker, 1979) models. In the Penman-Monteith model, the potential evapo-transpiration estimates are obtained by using the following relationship:

(85)

where:

• Δ is the slope of the specific humidity/temperature curve between the air temperature and the surface temperature of the vegetation (kPa °C^-1)
• λ is the latent heat of vaporization of water (2.50036 MJ/kg)
• c_p is the specific heat of air at constant pressure (1.013 kJ kg^-1 °C^-1)
• γ is the psychrometric constant kPa °C^-1)
• r_a is the aerodynamic resistance to the transport of water vapor from the surface to the reference level z (s/m)
• r_c is the (Monteith) canopy resistance (s/m) to the transport of water from some region within or below the evaporating surface to the surface itself. The canopy resistance is expected to be a function of the stomatal resistance of individual leaves. Under wet-canopy conditions r_c = 0.
• A_* is the available energy given by A_* = (R_l + R_s) - G (W/m²)
• R_s is the net incoming short-wave radiation at the reference level z, (W/m²) Section 9.1.1.1.
• R_l is the net long-wave radiation at the reference level z, (W/m²) Section 9.1.1.2.
• G denotes the sum of energy fluxes into the ground, to adsorption by photosynthesis and respiration and to storage between ground level and z (in W/m²), Section 9.1.1.3.

The standard reference level, z, is taken as 2 m. The current formulation does not explicitly calculate the leaf temperature of the vegetation canopy, and the ground temperature is substituted for the leaf temperature when calculating Δ. The psychrometric constant, γ, is defined as:

(86)

where: P is the atmospheric pressure (kPa). The gradient of the saturation vapor pressure curve with respect to temperature, Δ, is given by:

(87)

The rate of water diffusion from the ground surface due to turbulence is controlled by the aerodynamic resistance term, r_a. This term is a function of wind speed and the height of the vegetation cover. Mathematically, r_a is represented as:

(88)

where: z_u and z_e are the heights of the wind speed and humidity measurements (m), respectively, and U_z is the wind speed (m/s). Following Brutsaert (1975), z_om is assumed to be 0.123h_c and z_ov is assumed to equal 0.0123h_c, where h_c is the mean height of the crop (m). According to Monteith (1981) d = 0.67 h_c.

9.1.4 Parameter Values

Calculation of evapo-transpiration requires additional parameter values be assigned to every active grid cell. These parameters may be assigned with either the Mapping Table File, Section 11, or GRASS ASCII maps specified with project cards described in Section 3.8. ET parameters are typically assigned with a combination soil texture/land use (STLU) index map. The Deardorff method requires values of land surface albedo. For the Penman-Monteith method, values of land surface albedo, vegetation height, vegetation canopy resistance, and vegetation transmission coefficient are needed.

9.1.4.1 Land Surface Albedo

Ground surface temperature calculations, discussed in Section 9.1.2, require values of land surface albedo, which describe the fraction of long-wave radiation reflected back to the atmosphere. Values range from 0.0 to 1.0. Literature values for a variety of land covers compiled from a number of sources are listed in Table 11.

Table 11 – Land surface albedo values Notes: ^aThe smaller value is for high zenith angles; the larger value is for low zenith angles. From ^bSellers (1965); ^cMunn (1966); ^dRosenburg (1974); ^eOke (1973); ^fLee (1978); gde Jong (1973); ^hAtwater and Ball (1981); and ⁱEagleson (1970).

9.1.4.2 Canopy Stomatal Resistance

Plants lose water through their leaves through a process called transpiration. This water is lost through small opening in the leaves called stomata. Inside the leaf the air is near saturation, and the difference between air saturation and leaf saturation causes the plant to lose moisture due to diffusion. The plant can control the loss of water by opening and closing the stomata. The loss of water vapor from plants can be calculated from a resistance law, where the difference in saturation between the interior of the leaf and the atmosphere is the forcing term, and the canopy stomatal resistance is the resisting term.

The canopy stomatal resistance values entered in GSSHA must represent the resistance of the canopy for an entire grid cell. The canopy resistance is affected by coverage, time of day, and the type and condition of plants in the cell. Greater leaf coverage means lower resistance values. In GSSHA, noontime canopy resistance values are entered. Canopy resistance has a strong diurnal variation and the GSSHA model corrects the noontime canopy resistance for the time of day. Although stomatal resistances are not commonly available for mixed-vegetation types, data are available for a variety of single vegetation types (e.g. Szeicz and Long, (1969) for grass; Stewart and Thom (1973) for pine forest; Allen et al. (1989) for clipped grass and alfalfa). Noontime values of canopy resistance (s m^-1) for different types of plant under different conditions are listed in Table 12.

Table 12 - Canopy Stomatal Resistance Notes: ^aPielke (1984); ^bGash and Stewart (1975), ^cMonteith (1975); and ^dSceicz and Long (1969).

The Penman-Monteith equation has been found to be very sensitive to the value of canopy resistance. As pointed out by Lemeur and Zhang (1990), a 10% error in canopy resistance will result in a 10% error in the estimated evapo-transpiration. Senarath et al (2000) found that the discharges calculated during long-term simulations with CASC2D were also sensitive to the canopy resistance.

In the northern hemisphere temperate zone ET is subject to strong seasonal variations. ET is dependent on both climatic conditions and the vegetative cover. The seasonal variability of climatic conditions is reflected in the model with the hourly hydrometeorological inputs (Senarath et al., 2000). Vegetative cover is represented with simple land use/land cover indexes, such as forest, pasture, etc. Seasonal changes in vegetation cover are best simulated with plant growth models. Comprehensive plant growth models, such as SWAT (Arnold et al., 1999), and hydrologic models that include comprehensive plant growth models, TOPOG_IRM (Dawes et al., 1997) require extensive information on plant communities and growing conditions. Such detailed data are not routinely available.

Senarath et al. (2000) determined that of the ET parameters used in the Penman-Monteith equation, evapo-transpiration is most sensitive to the value of canopy resistance, and is quite insensitive to the other ET parameters. Leaf area and canopy resistance can vary by as much as several hundred percent during the year for crops, grasses, and deciduous forest in temperate regions (Monteith, 1975; Doorenbos and Pruitt, 1977; Federer and Lash, 1978). Based on this information, the canopy resistance was chosen as the vehicle to incorporate seasonal variability of ET in the GSSHA model.

Figure 13 – Seasonal Variation of Canopy Resistance in GSSHA

Mid-growing season values of canopy resistance are input as the starting point. For each month an amplification factor is used to represent the change in the canopy resistance related to plant growth. This mid-season value is then applied directly for the months of May through September. Thus, the amplification factor for May-September is 1.0. For the months November through February, the amplification factor is 4.0. This high factor relates to the death of crops, the browning of grasses and loss of leaves in deciduous trees. The months of March and October are considered transition months, and an amplification factor of 2.5 is used during these months. The timing of these events corresponds to values used for the transpirational leaf area of deciduous forest in North Carolina (Federer and Lash, 1978), though similar timing in conductance and ET is seen in other areas of the continental US (Nixon et al., 1972 for example). The implied assumptions of this simplistic approach, depicted in Figure 13, are then:

1. Beginning of spring - at the end of March grasses begin to green, crops to sprout and trees leaf;
2. Beginning of summer - by the first of May crops are near mature, and trees have full foliage;
3. End of summer - by the beginning of October crops begin to die, grasses to brown and trees lose leaves; and,
4. Fall/Winter - for the period November through February, plants are dormant and transpire little water to the atmosphere.

Use of this simple method at the Goodwin Creek Experimental Watershed (GCEW) in north Mississippi resulted in improved predictions of soil moistures and outlet discharge for periods outside the summer growing season (Downer, 2002a). The timing of seasonal changes should be adjusted to make the method applicable to climates different from that of southeast region of North America.

In GSSHA versions 5.x and higher, the above noted default values are applied for latitudes of 37 degrees and below. For latitudes greater than 37 degrees spring comes later and summer comes earlier; the four inflection points are April, July, September, and October. In versions 6.1 and above, these inflection points can be specified in the project file.

Seasonal variability of canopy resistance is selected by placing the SEASONAL_RS card in the project file along with the LONG_TERM card. In versions 6.1, and above, the inflection points on the curve can be specified by the user by using the cards SEASONAL_RS_SPRING SEASONAL_RS_SUMMER_START SEASONAL_RS_SUMMER_END and SEASONAL_RS_FALL. These cards are followed by the month, integer value, which the inflection point occurs. This capabilty allow users in the southern hemisphere or in regions of seasonal extremes to specify the timing on the inflection points. As shown in Figure 13, there is a linear transition between inflection points. ALL, or none, of the inflections points must be specified. For example, for the southern portion of the United States, below 37 latitude, the following cards describe the default conditions:

LONG_TERM

SEASONAL_RS

SEASONAL_RS_SPRING 3

SEASONAL_RS_SUMMER_START 5

SEASONAL_RS_SUMMER_END 9

SEASONAL_RS_FALL 10

9.1.4.3 Vegetation Height

The vegetation height is required to compute the aerodynamic resistance term in the calculation of turbulent diffusion. Vegetation heights are entered in meters. Sample values are listed in Table 13 (from Eagleson, 1970). The values may not be the representative, expected vegetation height-values of these vegetation/forest types.

Table 13 - Sample values of vegetation height

9.1.4.4 Vegetation Transmission Coefficient

The plant canopy can prevent radiation from reaching the ground surface, reducing the amount of radiation available to heat the ground surface and produce evaporation. The vegetation transmission coefficient describes the fraction of light that penetrates the vegetation canopy and reaches the ground. Values can range between 0.0, total canopy blocking of sunlight, to 1.0, total light penetration on bare soil. Table 14 list measured values of canopy resistance for grass.

Table 14 - Vegetation transmission coefficient values for grass Notes: Data from O.G. Sutton, “Micrometeorology,” McGraw Hill, New York, 1953.

9.2 Hydrometeorological Data

In order to perform continuous simulation in GSSHA hydrometeorological data are required for the entire period of the continuous simulation. The required data are hourly values of:

• barometric pressure,
• relative humidity,
• total sky cover,
• wind speed,
• dry bulb temperature,
• direct radiation, and
• global radiation.

These data are available from a wide variety of sources and three different input file formats can be used. The GSSHA model also contains provisions to synthetically produce solar radiation estimates and fill in short term data gaps.

9.3.1 Hydrometeorological Input File Formats

The three file formats that GSSHA supports are discussed below in order of preference. Simple file formats are preferred because the likelihood of a data format error is reduced. The WES format is the simplest because it uses no alphabetical characters. Users are encouraged to produce and use hydrometeorological data files in the WES format for continuous simulations. The other formats contain data quality flags, extraneous data, and occasionally errors. Such errors can be very difficult to locate. Reducing data down to the WES format through independent software (such as a spreadsheet) ensures that the data are properly formatted and contains no extraneous values.

9.3.1.1 HMET_WES (Recommended format)

This file contains only hourly values of the required hydrometeorological parameters: barometric pressure, relative humidity, total sky cover, wind speed, temperature, direct radiation, and global radiation. The HMET_WES format contains the following numbers in columns 1-11 (Table 15):

Table 15 – HMET_WES file format

Data format is important; the year is stored as a 4 digit number; the hour varies between 0 and 23. The data in columns 1-4, 6-9 are stored as integers, while the data in columns 5, 10, and 11 are stored as real numbers. The data in this file should be space delimited. The number of spaces between columns does not matter, provided that the total width of the line is less than 256 characters. There must be 1 hourly record for each hour of the intended duration of the continuous simulation, even if there are no data available for that hour. Gaps in the hourly record will cause the GSSHA simulation to terminate prematurely.

The following are three example lines from a HMET_WES formatted hydrometeorological data file (note the distinct presence of both real and integer numbers):

1995    4   14   23   29.625    86   100   000    53   9999.99   9999.99
1995    4   15    0   29.625    83   100   003    53   9999.99   9999.99
1995    4   15    1   29.615    90    80   003    52   9999.99    9999.99

This data input format is strongly recommended because it is compact, easy to read and write and contains no extraneous information. The no data flags are important. GSSHA has a set procedure for dealing with missing data.

9.3.1.2 HMET_SAMSON (Recommended over surface airways format)

In 1993 the National Climatic Data Center (NCDC) published a CD-ROM containing surface meteorological observations from 1961-1990. This CD-ROM is the result of collaboration between the National Oceanic and Atmospheric Administration (NOAA) and the U.S. Department of Energy, National Renewable Energy Laboratory. The SAMSON CD-ROM represents an excellent nationwide dataset for the 1961-1990 period. Unfortunately, this CD-ROM is not being updated as new data become available.

The following is an example of data from the SAMSON CD-ROM. A reader program provided with the SAMSON data produces this data format. The ASCII Samson data format includes a one line station header (in this case the data are from the Memphis, Tennessee station), and a one line descriptive header with column ID’s. These two header lines are followed by a user selectable number of hourly data points (in this example two hourly lines are shown).

~13893 MEMPHIS                TN  -6  N35 03  W089 59    87
~YR MO DA HR I    1    2       3       4        5       6  7       8     9   10     11   12    13     14       15     16   17         18   19        20  21
   82   4      1      1 0    0    0    0 ?0    0 ?0    0 ?0  0  0  12.8   6.1  64 1011 100  2.1   24.1  77777 999999999   16 99999.    0  26
   82   4      1      2 0    0    0    0 ?0    0 ?0    0 ?0  0  0  12.2   6.1  67 1011   80  2.1   24.1  77777 999999999   16 99999.    0  26

The first header line consist of the Weather Bureau Army Navy (WBAN) number for the station (13893), followed by 26 character descriptive string (MEMPHIS), the state name where the station is located (TN), the time difference between local standard and UTC in hours (-6), degrees and minutes of latitude, N|S, (N35 03), degrees and minutes of longitude ,E|W, (W089 59) and the station elevation in meters above sea level (87). The first header record containing the station information was written by NCDC with the following FORTRAN format statement:

(1X,A5,1X,A22,1X,A2,1X,I3,2X,A1,I2,1X,I2,2X,A1,I3,1X,I2,2X,I4)

The second line is for the identifier record, which shows the numerical identifier for the data elements in the National Solar Radiation Data Base (NSRDB) Synoptic Format (see: http://rredc.nrel.gov/solar/pubs/NSRDB/).

The observed data occupy the third record and all ensuing records (lines) in the file until a new year of data occurs. When a new year of data is encountered the header record and identifier record are repeated, followed by more data records. The data records include year, month, day, hour (1-24) in local standard time. The number below the “I” column is an “observation indicator” for which 0 means the weather observation was made, and a 9 indicates that the observation is missing. The following table describes the contents of each data record (Table 16):

Table 16 – HMET_SAMPSON file format

The SAMSON CD-ROM is available by contacting the NCDC at:

Climatic Services Branch

National Climatic Data Center

Federal Building

Asheville, NC 28801-2696

Phone - (704) 271-4800

Fax - (704) 271-4876

9.3.1.3 HMET_SUFAWAYS - NOAA/NCDC Surface Airways Format (Not recommended)

The Surface Airways data format is included in the list of supported formats because surface meteorological data ordered from NCDC are delivered in this format. In the Surface Airways data format each data element (wind, temperature, radiation, etc.) is written on a separate line in the file. Each line contains 24 hourly values. Missing data are not flagged.

The details of this data format are discussed in “TD-3280, Hourly Surface Airways Observations”, published by the National Oceanic and Atmospheric Administration, National Climatic Data Center. GSSHA users are strongly encouraged to reformat Surface Airways observations into the WES data format. There is a strong likelihood that data in the Surface Airways format will produce problems during the GSSHA simulation. This occurs not as consequence of faulty data, but because of the large number of data flags in the Surface Airways output format make the Surface Airways data format difficult to read in a robust fashion.

9.3.1.4 HMET_ASCII Gridded Data (Recommended format)

(Version 7 and higher) This format is for data that comes from other models or is otherwise pre-processed into gridded data. The grid shape and size of the hmet data does not need to match the GSSHA grid (and usually does not), but it does need to match itself (be consistent across data types and over time.)

The data are in individual grid files for each data type and each time step. They are usually placed in a folder by themselves as they will be quite numerous. They have a specific naming convention that GSSHA relies on: YYYYMMDDHH_TYPE.asc

The data for the ASCII Gridded format follows the HMET_WES list of data.

The type flags and units for each class of data are as follows:

The hmet file format is straightforward: it is simply the list of paths and date/time stamps for each hour. GSSHA will take the path and date/time info and append the underscore and type name to open the individual data files.

In this example, the files in the HMET_DATA folder would look like this:

To use gridded ascii data in your model, in the project file use the card HMET_ASCII "filename.hmet".

9.2.2 Missing Hydrometeorological Data

Frequently, hydrometeorological records have periods of missing data because observational instruments are inoperative, the data are not reliable, or certain observations (e.g. direct radiation, global radiation) are not available for a given station. Because hourly hydrometeorological inputs are required for the entire simulation period, missing data are synthetically produced within the GSSHA model.

9.2.2.1 Radiation Data

Many hydrometeorological stations do not record solar radiation measurements. Missing solar radiation data are simulated in GSSHA based on the location of the watershed, day of the year, and time of day, as described in Section 9.1. When only a portion of the solar radiation data is missing, it is replaced as described in the next section.

9.2.2.2 Short Duration Gaps in Hydrometeorological Measurements

When required data are missing, GSSHA creates the needed data by either persistence or persistence adjusted for the time of day. Missing variables without strong diurnal variability are filled in by “persistence” estimates. These variables include:

• barometric pressure,
• wind speed,
• dew point temperature and,
• percent cloud cover.

Persistence means that the value is held constant until a new observation becomes available. Missing variables with strong diurnal components:

• air temperature,
• global radiation,
• and direct radiation.

are replaced with the last good reading from the same time of day. It is the responsibility of the user to check for missing data. Extended periods of missing record (many days) will cause GSSHA to simulate conditions for the last day of good record over and over again. This scheme requires that the hydrometeorological data file begin with at least one day of good observations of all required variables.

9.3 Snowfall Accumulation and Melting

Snow melt is the single most important event of the water year in many parts of the world (Gray and Prowse, 1993). Snow accumulation and melt are both spatially heterogeneous, thus distributed domains often give more potential for addressing a variety of real world problems (Kirnbauer, Bloeschl et al. 1994). When GSSHA is run in the LONG_TERM simulation mode, snowfall accumulation and melt are automatically simulated, with the hydrologic response being a combination of the multiple other processes (infiltration, overland flow, groundwater, etc.) utilized within GSSHA. GSSHA employs multiple snow model options to use depending on available data and the objective of the study. These models will be discussed below.

Because of the distributed domain, the model can also account for Orographic Effects.

Within GSSHA, the snow pack consists of a single layer in each cell. While multi-layer models have certain advantages - such as time variations of liquid water content (Bloschl & Kirnbauer 1991), interflow within the snow pack layers due to ice sheets, and avalanche modeling (Colbeck 1991) - the required data to populate model such models on a watershed level is currently unrealistic in most basins. Multi-layer snow models are typically deployed at the site-scale, where spatially-close data is available. The methods used in GSSHA require only the standard precipitation, Chapter 6, and hydrometeorological data, Chapter 9.5, required for any long term continuous simulation in GSSHA, and are readily available from a variety of sources. While optional inputs are possible, the model can produce accurate results with no additional inputs other than standard precipitation and hydrometeorological inputs. The methods currently used in GSSHA to simulate snow accumulation and melt can be found here Follum and Downer 2013

Snow Related Cards for the GSSHA Project File
No additional input cards are required to simulate snow accumulation, melt, and meltwater routing. The default is for precipitation to be added to the snowpack anytime the air temperature, as specified by the hourly input values, is below 0 C. Melt occurs according to the hybrid energy balance method, as described below. Flow through the snow pack and to the overland occurs in waves. Overland routing will ignore the snow pack. All of these defaults can be altered as described below.

Snow Card Inputs - Optional

GSSHA Snow Accumulation
Precipitation accumulates as snow when Air Temperature is less than PXTEMP, where PXTEMP is the temperature (°C) at which precipitation is considered snow (Default of 0°C). PXTEMP was omitted from Follum & Downer (2013) in error, but included in Follum et al. (2015).

Because gaging systems often underestimate the amount of precipitation fallen in the form of snow, two multiplication factors are applied to the gage measurement of precipitation (P_x) to correct the precipitation of newly fallen snow (P_n), Equation 6. The snow adjustment factor (SCF) is considered a calibration parameter while the fraction of precipitation in the form of snow (f_s) is considered constant at 1.0 when temperature are at or below 0° C, and 0.0 when above 0° C.

If snowfall occurs, a warning will be printed to the screen and to the Summary file. When snow accumulation occurs the amount of snow in the watershed is reported at the beginning and end of each event summary in the Summary file.

GSSHA Melt Algorithms

 Energy Balance (EB)

 Temperature Index (TI)

 Hybrid Energy Balance (HY) Default 

 Radiation-derived Temperature Index (RTI) 

Adjustment of Atmospheric Radiation Fluxes within the Snow Model
GSSHA currently employs methods to account for both longwave and shortwave radiation in each cell. Longwave radiation is mainly a function of temperature, clouds, and atmospheric emissivity, while the shortwave radiation calculations take into consideration albedo, topographic shading, aspect of the terrain in relation to the sun, albedo of snow as it ages, atmospheric absorption and reflection, clouds, and vegetation. The calculated longwave and shortwave radiation values are then used within the EB and HY models to simulate the melting of snow.

In order to account for the Effects of Shading and Aspect on the available energy to melt snow, GSSHA employs a method that reduces the amount of energy available to melt snow when the Energy Balance or Hybrid Energy Balance snow melt algorithms are used. The Effects of Shading and Aspect are set as a default within GSSHA and are only used for snow fall, but are currently being examined in the Evapo-Transpiration methods as well. Shading and Aspect angles are not accounted for when raster-based HMET data is input into the model (http://www.gsshawiki.com/Distributed_HMET_Data).

GSSHA Melt-Water Transport Algorithms
Once melt occurs within the snow pack the subsequent melt-water is routed through the snow pack and then to other hydrologic processes within GSSHA (i.e. runoff, infiltration, evaporation, etc.) using melt-water transport (MWT) algorithms. GSSHA currently employs both a Vertical MWT and Lateral MWT algorithm.

Within a grid cell, a homogeneous snow pack assumption is utilized in GSSHA to alleviate computational and data limitation concerns associated with a heterogeneous assumption (which would include flow fingers). Equivalent properties for the homogeneous snow pack are often assumed (Colbeck 1979). In GSSHA, each cell has its own snow pack properties - namely depth, density, hydraulic conductivity, saturation, and effective porosity - derived using the SNAP model (Albert & Krajeski 1998).

Flow is considered a porous medium, therefore a form of Darcy's Equation (http://en.wikipedia.org/wiki/Darcy%27s_law) is used to determine flux rates through the snow pack. Vertical MWT through the snow pack is considered unsaturated flow, while Lateral MWT between the ground surface and the bottom of the snow pack is considered saturated flow. To simulate the flow within the snow pack, accurately capturing the saturation within the pack is vital because the saturation affects both the hydraulic conductivity and effective porosity of the snow pack. GSSHA currently uses the SNAP model to determine the saturation, saturated / unsaturated hydraulic conductivity, and effective porosity of the snow pack in each grid cell. Water is routed to the overland through the snow pack in a series of waves calculated from these parameters.

GSSHA Snow Depth and Density
Snow depth and density are simulated in GSSHA independent of what melting algorithm is used. The depth and density of snow is calculated hourly by incorporating the SNAP model (Albert and Krajeski 1998) code into GSSHA. Information related to SWE, snow depth, density, snow saturation, effective porosity, and hydraulic conductivity are exchanged between the two models. The SNAP model algorithms are used to calculate all the parameters except SWE.

According to the SNAP formulation, the snow density changes in response to snow accumulation, settlement, and melt, and is computed either with the snow depth predictions, or as updated by the user (Albert and Krajeski 1998). GSSHA currently does not allow the user to update the snow density as a calibration parameter, leaving the snow density to be calculated based on the snow depth predictions of the SNAP alorithms with no calibrated parameters. For more information on the SNAP model the interested reader is encouraged to review Albert & Krajeski (1998). For more information on the depth prediction equations used within the SNAP model the interested reader is encouraged to review Anderson (1973), Jordan (1991), and Jordan (1998).

The snow depth simulations from Follum & Downer (2012) show that when the SWE is simulated accurately the SNAP model reasonably simulates the snow depth. The results show that the SNAP model may slightly overestimate snow depth when compared to observed data Follum and Downer 2013.

Overland Routing with Snow
During overland flow routing, as described in Section 5.2, GSSHA ignores snow in an overland flow cell unless the user specfies to route the flow through the snow using Darcy's law with the project card ROUTE_LAT_SNOW. Routing through the snow as free surface flow may cause simulated flows to be higher and arrive earlier than measured flows strongly influenced by the snowpack. When this card is specifed, the flow in cells with a snowpack will be computed using Darcy's law. The default is to use the SNAP calculated vertical hydrualic conductivity (K) for computation of flow through the snow in the lateral flow computations. The alternative is to specify the lateral hydraulic conductivity with the SNOW_DARCY card, which is followed with a value of hydraulic conductivity (m/s). References report hydraulic conductivities of snow on the order of 1 cm/s (0.01 m/s) (Colbeck and Anderson, 1982 for example). The current implementation of the SNAP model produces simliar values but as of v6.2 the implementation of SNAP in GSSHA is considered experimental and is not currently reccomended and it is reccomended that the user specify the lateral hydrualic conductivity of the snow using the SNOW_DARCY card.

9.4 Sequence of Events During Long-Term Simulations

At the beginning of a long-term simulation, GSSHA opens the hydrometeorology data file, and determines the date/time of the beginning and end of the hydrometeorology data. GSSHA begins at the date/time specified as the beginning of the hydrometeorology data file. If the first rainfall event in the rainfall input file does not begin during the period of hydrometeorology data, then GSSHA aborts with a warning.

In LONG_TERM mode control of process activation switches from rainfall events to the physical conditions in the watershed. When there is no moving water on the overland flow plane or in the channels, ET demand and unsaturated water movement are updated at an hourly bases. This is possible because, altough the overall model timestep is used to loop through time in the model, each process has it's own individual update time. Therefore, computations for any process can be started and stopped at any time, and any process updates can occur at multiple overall timesteps. This allows signficant time savings by eliminating needless calculations. Altough the internal control is complicated, it is generally described below.

When water appears on the overland flow plane due to rainfall, exfiltration or channel overbank flow, overland flow routing and unsaturated zone calculations begin at the overall model time step, typically around 60s or less. Channel routing begins as water in the channel accumulates, also at the overall model time step. Overland flow calculations continue until all water on the overland flow plane runs off, infiltrates, evaporates, or stops moving. Channel routing continues as long as a minimum amount of water is in the channel and the water is moving. If groundwater/channel interactions are specified, these fluxes continue to be computed and the volume in the channel is updated each groundwater time step. If exchange with the groundwater results in the minimum amount of water required for channel routing to begin, then channel routing will resume, but at the larger groundwater time step. Anytime rainfall is occurring, overland flow, ET, channel flow, unsaturated zone water movement, and saturated groundwater flow are updated each overall model time step. Events are still recognized, but are only used for accounting purposes.

When continuous simulations are performed GSSHA also outputs a file called evap_out which list basin wide average hourly values of: soil moisture, potential evaporation (cm hr^-1), actual ET (cm hr^-1), infiltration rate (cm hr^-1), groundwater recharge (cm hr^-1), continuous frozen ground index, temperature (C), gridded temperature (C), surface soil temperature (C).

10.1 Overland Erosion Formulation

10.1.1 Overland Erosion Controlling Equation

The overland sediment routing is based on a two-dimensional mass balance equation,

(96)

where h = water depth (m); c = sediment concentration (kg m^-3); t = time (s), q_x and q_y = unit discharge in x- and y-directions, respectively (m²/s); and e(x,y,t) = sediment source/sink (kg m^-2• s^-1); q_z(x,y,t)=lateral sediment inflow to the channels (m^-2 s^-1). e(x,y,t) consists of two major sources, namely, soil detached by raindrop impact and soil detached by surface runoff (or deposition):

(97)

where d_i=soil detachment by rainfall impact (kg m^-2• s^-1) and d_r=soil detachment by surface runoff (kg m^-2• s^-1).

10.1.2 Detachment of Overland Sediment Particles

Detachment of overland sediment particles occurs due to two processes, detachment by rainfall impact, and detachment by surface runoff.

10.1.2.1 Detachment by Rainfall Impact

Following a similar approach presented by Wicks and Bathurst [1996] and Gabet and Dunne [2003], detachment by raindrops is considered to be a function of rainfall intensity. Raindrop detachment also takes into account factors such as “cushioning” by ponded surface water, ground cover, and plant interception.

(98)

where d_i = detachment capacity rate (kg m^-2 s^-1), K_I = soil erodibility factor for detachment by raindrop impact (J^-1), C_w = water depth correction factor, C_G = canopy cover factor, C_i = a cover-management factor, I = rainfall intensity (mm h^-1), and β = empirical coefficient related to rainfall intensity.

The calculation of water depth correction factor is based on the assumption that surface water depths greater than a critical depth protect the soil from raindrop impact erosion Wicks and Bathurst [1996]. The expression for the correction factor, C_w, is,

(99)

where h = water depth (m), and D_m = median raindrop diameter (m). The median raindrop diameter is determined from the Laws and Parsons [1943] equation:

(100)

10.1.2.2 Surface Runoff Detachment

Surface runoff detaches soil particles by exerting a shear stress that breaks the bonds between particles. Erosion in rills is lumped and described as gross rill erosion. Within a grid cell, rill erosion and flow are assumed to be uniformly distributed. The detachment capacity rate by surface runoff has the form:

(101)

where: d_r = detachment capacity rate (kg m^-2•s^-1), a and b are empirical coefficients, τ = the flow shear stress (Pa), τ_cr is the critical shear stress, G is the sediment load (kg m^-2 s^-1), and T_c is the sediment transport capacity of surface runoff (kg m^-2 s^-1). In the WEPP model, b is assumed equal to 1.0. Equation 101 provides an upper limit on the rate at which sediment that can be taken up by flowing water.

10.1.3 Sediment transport capacity of surface runoff

Transport of sediments on the overland flow plane is constrained by the ability of the water to carry sediments. Many researchers have investigated the process and have identified multiple equations that can be used to calculate this carrying capacity of the overland flow. Several of these have been implemented in GSSHA.

10.1.3.1 Kilinc Richardson Equation

Kilinc and Richardson studied the mechanics of overland soil erosion at Colorado State University Engineering Research Center in a flume 1.2 m deep, 1.5 m wide, and 4.9 m in length. Their investigation resulted in a sediment transport equation of uniform flow sheet and rill erosion on bare sandy soil. The original one-dimensional Kilinc-Richardson (1973) equation is:

(102)

where: the factor 25,500 is an empirical constant, is the sediment unit discharge (ton m^-1 s^-1) and q and S_o have previously defined. The form of Equation 102 is consistent with many sediment transport predictors, as reviewed by Julien and Simons (1985).

Julien (1995) modified the original Kilinc-Richardson equation to expand the applicability of the equation to non-uniform flow with consideration of soil and land-use specific factors:

(103)

where: K = soil erodibility factor, with values ranging from 0 to 1, C = soil cropping factor (0-1), and P = conservation factor (0-1).

The K, C and P factors are empirical coefficients with the same conceptual meaning as those used in the Universal Soil Loss Equation (Renard et al., 1991). The factor 0.15 in the denominator represents the maximum erosivity of the sand used in the original flume experiments. Equation 103 is applied within each model grid cell in two dimensions, which allows separate calculations of the potential sediment transport rate in the x- and y-directions:

(104)

(105)

with S_ƒx and S_ƒy defined using Equation 2, and p and q calculated using Equations 6 and 7, respectively.

The numerical approach in the current version of GSSHA is to apply the modified Kilinc and Richardson equations 104 and 105 to determine potential transport of sediment within each grid cell. The methodology is explained in detail in Johnson (1997) and Johnson et al. (2000). The potential sediment discharge in each of the two flow directions in each grid cell is calculated using the following two equations:

(106)

and,

(107)

In the current version of GSSHA, the three erosioin parameters K, C, and P, have been replaced with a single erosion parameter K, which represents the three original components, without the need to specify three different parameters for each grid cell. Reducing the number of parameters greatly facilitates calibration. The form of the equation used in GSSHA is currently:

(108)

where q_s is the sediment unit discharge (kg m^-1 s^-1), q is the unit discharge (m³ s^-1), S_f is the friction slope (-), and K is a combined factor(0-1) that consists of the combined influence of soil erodibility, vegetation, and land-use.

10.1.3.2 Engelund-Hansen Equation

The Engelund-Hansen [1967] relation is used to calculate sediment transport for each sediment size class and the resulting total transport is calculated by multiplying the proportion of the size in the parent material by the calculated rate.

(109)

where G_i is the volumetric sediment transport rate of i^th size fraction (m³ s^-1), K is the calibration coefficient (=1 for standard equation) (-), F_i is the proportion of i^-th faction in the parent material or deposited layer (-), B is the width of flow (m), V is the mean water velocity (m s^-1), h is the flow depth (m), S is the water surface slope (-), s is the specific gravity of i^th fraction (-), and g is the gravitational acceleration (m s^-2), is the sediment density (kg m^-3), is the water density (kg m^-3), D_i is the mean size of i^-th fraction (m). The factor 0.05 in the equation was developed from empirical data.

The suggested applicability of the Engelund-Hansen equation is for, where "D_N" is the grain size for which N% of sediment is finer by weight, and for sand-size sediments coarser than 0.15 mm.

10.1.3.3 Stream Power Approach

Several stream power methods [Everaert 1991] can also be chosen to calculate the transport capacity of surface runoff. The general expression has the form:

(110)

where a is the transport coefficient, b is an exponent, ω = steam power; ω_c=a critical stream power. Prior to version 5.7 these values were specified by the user. For versions 5.7 and beyond the exponents are assigned based on literature values [Everaert 1991]. These there parameters are no longer input in the cmt input table file.

10.1.3.4 Median size diameter, D₅₀ , Approach

In the versions 5.7s and beyond, the Stream Power Approach defined in 10.1.3.3, are replaced by laboratory physical model imperical equations derived by Everaert [1991] which are based on the sediment particle median size diameter, D50. The replacement is done to reduce the number of calibration parameters, a, b and critical values defined in sec. 10.1.3.3.

The followings are the list of methods and their respective equations:

a) Effective stream power method:

- (a1)

where,

- (a2)

- (a3)

and <math>q_s</math> = sediment unit discharge (kg/ m/ sec)

<math>D50</math> = Grain size parameter (um)

<math>\Omega</math> = effective stream power

<math>\omega</math> = stream power (N/m/s)

<math>R</math> = flow depth (m)

<math>q</math> = unit discharge (m²/s)

<math>S_f</math> = friction slope (m/m)

<math>\rho</math> = density of water at 20 C (kg/m³)

<math>g</math> = gravitational acceleration

This method is implemented if 4 is used for the project card SOIL_EROSION.

b) unit stream power method:

- (b)

where, <math>Vel</math> = velocity (m/s)

<math>S_f</math> = friction slope (m/m)

<math>q_s</math> = sediment unit discharge (kg/ m/ sec)

<math>D_{50}</math> = Grain size parameter (um)

This method is implemented if 3 is used for the project card SOIL_EROSION.

c) slope and unit discharge method:

if(d50<33.0)

- (c1)

else if(d50>=33.0 && d50<61.0)

- (c2)

else if(d50>=61.0 && d50<122.0)

- (c3)

else if(d50>=122.0 && d50<190.0)

- (c4)

If(d50>190.0)

- (c5)

This method is implemented if 5 is used for the project card SOIL_EROSION.

d) Transport capacity by shear velocity method:

- (d1)

u = shear velocity (m/s)

- (d2)

Τ = shear stress (N/m2)

if(d50<33) and if(u>1.4)

- (d3)

if(d50>=33.0 && d50<61.0) and if(u>1.4)

- (d4)

if(d50>=61.0 && d50<122.0) and if(u>1.45)

- (d5)

if(d50>=122.0 && d50<190.0) and if(shear_vel>1.55)

- (d6)