Changing raster cells values inside of a polygon to the closest cells just outside the polygon!

I have been wrestling with this for a while and couldn't find out how I can do this. I want each cell of a raster inside of a polygon boundary to have the value of the closest cell just outside of the polygon.

I think the closest answer was the one at this question, but I couldn't apply it!

My plan was:

  1. Use "Polygon to Raster" to change the polygon to a raster. The cells inside the polygon become 1 and the points outside of the polygon are nodata.
  2. Use Con(IsNull(Raster),0,NULL) and assign it to Raster2
  3. Add Raster1 and Raster2 to get nodata inside the polygon
  4. It should be easy from here with nodata inside the polygon. We just use "Euclidean Allocation" and it does the job.

I'm stuck at step 2 and get an error:

[SELECT * FROM VAT_Extent_RAS WHERE IsNull( "Value" )] Failed to execute (Con).

Do you have any thoughts? Is there an easier way to do this?


I'll try to summarize my progress here. I want to change all raster cells inside the polygon to "nodata" while the raster cells outside of the polygon remain the same. I used polygon to raster function first. This gives me a raster with 0 cell inside the polygon and "nodata" outside of the polygon. What do I do next?

The Euclidean Allocation tool can accept your polygon as input, but you are right you need to set the values inside the polygon to NoData for this tool to work, otherwise it will pick up the values that are already there.

First to feature to polygon, this will give values (any value, it doesn't matter, we're only interested in value or NoData) inside the polygon and NoData outside.

Now do IsNull on the polygon raster, this will give you a raster with 1 outside the polygon and 0 inside the polygon. Hereafter called IsNullRaster. This can be done in a single statement but I want to use it twice so it's worth persisting it to disc.

Con with the IsNullRaster:Con(IsNullRaster,ValueRaster)- now you have a raster that is null inside the polygon and contains values outside the polygon. I'll call this one ValueRasterHole. If this is giving you grief do something likeCon(IsNullRaster,ValueRaster,255)and then use Set Raster Properties to enforce 255 as the nodata value.

Set your cell size, extent and snap raster to the ValueRaster and do Euclidean Allocation tool against the ValueRaster with the hole. This one is AllocatedRaster and will have values outside the polygon.

To merge the AllocatedRaster and the ValueRaster together use another Con with the IsNullRaster:Con(IsNullRaster,ValueRaster,AllocatedRaster)will fill in the hole in the polygon with the allocated raster.

Discrete natural neighbour interpolation with uncertainty using cross-validation error-distance fields

Interpolation techniques provide a method to convert point data of a geographic phenomenon into a continuous field estimate of that phenomenon, and have become a fundamental geocomputational technique of spatial and geographical analysts. Natural neighbour interpolation is one method of interpolation that has several useful properties: it is an exact interpolator, it creates a smooth surface free of any discontinuities, it is a local method, is spatially adaptive, requires no statistical assumptions, can be applied to small datasets, and is parameter free. However, as with any interpolation method, there will be uncertainty in how well the interpolated field values reflect actual phenomenon values. Using a method based on natural neighbour distance based rates of error calculated for data points via cross-validation, a cross-validation error-distance field can be produced to associate uncertainty with the interpolation. Virtual geography experiments demonstrate that given an appropriate number of data points and spatial-autocorrelation of the phenomenon being interpolated, the natural neighbour interpolation and cross-validation error-distance fields provide reliable estimates of value and error within the convex hull of the data points. While this method does not replace the need for analysts to use sound judgement in their interpolations, for those researchers for whom natural neighbour interpolation is the best interpolation option the method presented provides a way to assess the uncertainty associated with natural neighbour interpolations.

A grid algorithm suitable for line and area feature label placement

The labelling problem has been central in the framework of automated cartography. The quality and efficiency of label placement have great influences on the expression and understanding of maps. Although many algorithms have been developed to address the labelling problems of point features, very little work has been directed towards those of line or area features. Owing to the weakness of these approaches, the label quality rules of line or area features were reconsidered and strengthened based on the cognizance of cartographers. Such rules should be separate from the labelling algorithms to be appropriate for the program’s flexibility. A new grid algorithm, in contrast to traditional vector-based methods, is proposed. For the line feature, the cells passed by a line are computed, and their parallel cells are selected as the bottom of the text. For the area feature, a maximal inclusive rectangle is searched for the numerical label of its corresponding polygon (area), the midpoint of which is considered the potential position. A test program was developed and shows that the algorithm is simple and appropriate. The efficiency of the algorithm is closely related to the cell density.

This is a preview of subscription content, access via your institution.

Perceptions and expected immediate reactions to tornado warning polygons

To provide people with more specific information about tornado threats, the National Weather Service has replaced its county-wide warnings with smaller warning polygons that more specifically indicate the risk area. However, tornado warning polygons do not have a standardized definition regarding tornado strike probabilities (p s) so it is unclear how warning recipients interpret them. To better understand this issue, 155 participants responded to 15 hypothetical warning polygons. After viewing each polygon, they rated the likelihood of a tornado striking their location and the likelihood that they would take nine different response actions ranging from continuing normal activities to getting in a car and driving somewhere safer. The results showed participants inferred that the p s was highest at the polygon’s centroid, lower just inside the edges of the polygon, still lower (but not zero) just outside the edges of the polygon, and lowest in locations beyond that. Moreover, higher p s values were associated with lower expectations of continuing normal activities and higher expectations of seeking information from social sources (but not environmental cues) and higher expectations of seeking shelter (but not evacuating in their cars). These results indicate that most people make some errors in their p s judgments but are likely to respond appropriately to the p s they infer from the warning polygons. Overall, the findings from this study and other research can help meteorologists to better understand how people interpret the uncertainty associated with warning polygons and, thus, improve tornado warning systems.

This is a preview of subscription content, access via your institution.

Changing raster cells values inside of a polygon to the closest cells just outside the polygon! - Geographic Information Systems

TauDEM (Terrain Analysis Using Digital Elevation Models) is a set of tools for the analysis of terrain using digital elevation models. It incorporates programs and digital elevation model (DEM) analysis functions developed over several years of research. TauDEM is currently packaged as an extendable component (toolbar plugin) to both ESRI ArcGIS (8.x and 9.0) and MapWindow Open Source GIS.

  • Distribution and Copyright
  • Download and installation
  • * Usage Notes and Limitations *
  • Overview
  • Interface
  • A tutorial example to get started
  • The functions and what they do
  • Data formats and file naming conventions
  • * Frequently Asked Questions *
  • Updates, history and command line versions
  • Acknowledgements

Distribution and Copyright

Copyright (C) 2004 Utah State University

This program is free software you can redistribute it and/or modify it under the terms of the GNU General Public License version 2, 1991 as published by the Free Software Foundation.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

A copy of the full GNU General Public License is included in file gpl.html. This is also available at:
or from:
The Free Software Foundation, Inc., 59 Temple Place - Suite 330,
Boston, MA 02111-1307, USA.

If you wish to use or incorporate this program (or parts of it) into other software that does not meet the GNU General Public License conditions contact the author to request permission.
David G. Tarboton
Utah State University
4110 Old Main Hill
Logan, UT 84322-4110
email: david.tarboton<at symbol>

Download and Installation

System Requirements. Windows 2000 or higher. To use the ESRI ArcMAP toolbar you will need ArcGIS. The spatial analyst extension is not required for TauDEM to function, but is useful to work with the results. To use the MapWindow plugin you will need MapWindow 3.0 or higher.


  • TauDEM Setup package for ArcGIS 8.3
  • TauDEM Setup package for ArcGIS 9.0-9.3
  • TauDEM Setup package for MapWindow.
  • Source codes [2.8 Mb] This only includes the source for the TauDEM toolbar and plugin. ESRI ArcGIS required to run this software is not available for me to distribute.
  • Link to download MapWindow. MapWindow source code is available from


Execute the installation file to install the necessary libraries and components. For ArcGIS, open ArcMap. Click on tools/customize. At the Customize dialog click on add from file . Browse to c:program filesTauDEMagtaudem.dll then click Open. Click OK to added objects. The entry "Terrain Analysis using digital elevation models (TauDEM)" should appear under Toolbars. Check this (if not already checked) and close the customize dialog. The TauDEM toolbar should now be present in the ArcMap environment. If it does not appear the first time, close and reopen ArcMap. For MapWindow click on Plugins and select TauDEM. The TauDEM menu should appear. If it does not check that the folder taudem created by the install is in the correct location, which is by default c:program filesmapwindowplugins.

Usage Notes and Limitations

  1. Grid File and Path Names: Spaces are not allowed in file names or paths (e.g. "My Documents" has a space, and so will not work). File names must be 13 characters or less.
  2. Grid Size: Grids must be smaller than 7000 x 7000 cells.
  3. Grid Spatial Reference/Size: All input grids are assumed to be the same size, shape, and in the same spatial reference. This is sometimes checked by the program, but not consistantly.


TauDEM (Terrain Analysis Using Digital Elevation Models) incorporates the Digital Elevation Model (DEM) analysis tools and functions developed by David Tarboton over the years with support from a variety of sponsors, whose support is gratefully acknowledged.

The TauDEM plug-in currently provides the following capability:

  • Pit removal by flooding to ensure hydraulic connectivity within the watershed
  • Computation of flow directions and slopes
  • Contributing area using single and multiple flow direction methods
  • Multiple methods for the delineation of channel networks including curvature-based methods sensitive to spatially variable drainage density
  • Objective methods for determination of the channel network delineation threshold based on stream dropsStream ordering
  • Delineation of watersheds and subwatersheds draining to each stream segment and association between watershed and segment attributes for setting up hydrologic models.
  • Specialized functions for terrain analysis, including:
    • Wetness index
    • Distance to streams
    • Downslope influence function to map locations downslope that may be influenced by activities in an area
    • Upslope dependence function to map the locations upslope where activities have an effect on a downslope location
    • Decaying accumulation that evaluates upslope contribution subject to decay or attenuation
    • Concentration limited accumulation
    • Transport limited accumulation
    • Reverse accumulation


    The TauDEM plugin consists of the following three menus.

    Basic Grid Analysis contains the functions that are core to most basic Digital Elevation Model (DEM) analyses, and provide inputs to many other functions. The general order in which they should be executed is top to bottom.

    Network Analysis contains the functions required to delineate channel networks and subwatersheds. Again the sequence is top to bottom.

    Specialized Grid Analysis contains more advanced functions that can be invoked as needed.

    The section below on the functions and what they do below describes in detail how each function works. Generally each time a menu item is selected a dialog box appears showing the inputs and outputs for a particular function. Once a base DEM has been selected the file names for most inputs and outputs are populated following the default file naming convention. You are free to change these, but working with the defaults saves some typing. [A side effect of this, is that if a new base DEM is selected, all file names are reset back to their defaults.] Before running each function you should ensure that all necessary input information exists, or expect an error. Help on the individual files is available by clicking on the label to the left of the field where file names can be changed.

    A tutorial example to get started

    Download and unzip the files in into a convenient place for you to work.

    1. Use ArcToolbox | Import to Raster | ASCII to grid to import the file "demo.asc" as a grid named "demo"

      Use the float option for Grid type.
    2. Open ArcMAP and add the TauDEM toolbar. [Click on tools | Customize | Add from file and select the file c:program filesTaudemagtaudem.dll]

      You should get a toolbar that looks like

      This may be docked.
    3. Use Add Data to load the grid named "demo". [OK to the message about missing spatial reference information.]
    4. From Basic Grid Analysis "Select Base DEM grid" and click OK with the Base DEM layer as "demo".

    1. Open ArcMAP and activate the TauDEM plugin.

    2. From Basic Grid Analysis "Select Base DEM grid".

      and browse to open the file "demo.asc or demo.bgd ". This identifies this as the Base DEM and sets up default file names for all other inputs.
    3. Invoke the functions in DEM Processing Functions in sequence from top to bottom, starting with "Fill Pits".

      then "D8 flow directions"

      and so on. Examine the output at each step. You could have done this all at once with the function "Do All".
    4. The layer named "demosrc.bgd" that results from "Full River Network Raster" is a rasterized version of channel network to be mapped. This provides a background for the positioning of outlets on the streams.
    5. Use the MapWindow shape editor plugin to create a new point shapefile named "myoutlets".

    6. The shapefile when created has no data yet. Use the shapefile editor to carefully place a new point on a stream path where you want an outlet. Multiple points for multiple channel networks can be added if desired.

    7. FromTauDEM menu "Select outlets shapefile . " and select the file "myoutlets.shp". If you have trouble creating this you may use the "demooutlet.shp" file provided.

    8. From TauDEM menu invoke functions in the Network and Watershed Processing Functions group in sequence from top to bottom, or select "Do All Network and Watershed Delineation Steps" and examine the output. If outlets are not selected, "Network Delineation" delineates all streams in the domain, and all watersheds draining to these streams. The methods and parameters for stream delineation are adjusted on the "River Network Raster . " dialog.

    The functions and what they do

    Each function is described below. Optional inputs are shown in [], with file name suffixes that refer to the table in the data formats section below given in (). TauDEM assumes implicitly that all grids used have the same grid cell size and extent. This assumption is checked in some places, but not consistently so where data is obtained from diverse sources it needs to be converted to a consistent grid cell size and extent (and the same projection) for use with TauDEM.

    Fill Pits

    Method: This identifies all pits in the DEM and raises their elevation to the level of the lowest pour point around their edge. Pits are low elevation areas in digital elevation models (DEMs) that are completely surrounded by higher terrain. They are generally taken to be artifacts that interfere with the routing of flow across DEMs, so are removed by raising their elevation to the point where they drain off the edge of the DEM. The pour point is the lowest point on the boundary of the "watershed" draining to the pit. This step is not essential if you have reason to believe that the pits in your DEM are real. This step can be circumvented by copying the raw DEM source data onto the file with suffix "fel" that is the output of "Fill Pits". Also if a few isolated pits are known, but others need to be filled, the isolated pits should have "no data" elevation values inserted at their lowest point. "no data" values serve to define edges in the domain, and elevations are only raised to where flow is off an edge, so an internal "no data" value will stop a pit from being filled if necessary.

    The flow path grid to enforce drainage along existing streams is an optional input. The flow directions in the flow path grid grid take precedence over flow directions determined from the DEM and where these are uphill, the elevations along these flow paths are lowered, rather than upflow elevations raised. The verified flow path grid output when the optional flow path grid is used has loops and ambiguities present in the original flow path grid removed.

    The enforcing of flow along a flow path should be used when the stream data source is deemed to be better than the DEM. The input flow path grid uses the D8 direction encoding, i.e. 1 - East, 2 - North East, 3 - North, 4 - North West, 5 - West, 6 - South West, 7 - South, 8 - South East. No data values indicate off stream locations. The flow path grid can be created by the network editor plugin in MapWindow or in ArcGIS by burning in a stream feature dataset using the following steps.
    1. Convert features to raster retaining the same cell size and extent as the target DEM. Call the resulting grid strgrd.
    2. Use raster calculator to subtract a large number from each elevation value that corresponds to a stream. This results in a temporary DEM with deep canyons along the streams. Call the resulting grid demcanyon.
    3. Use "Fill Pits" and "D8 flow directions" to calculate flow directions on demcanyon. The flow directions calculated will be demcanyonp.
    4. Use raster calculator to evaluate demcanyonp/strgrd. This will result in no data values off the stream raster due to a divide by 0, but will retain flow directions calculated on the stream raster. The convention for naming the result is to use the suffix fdr. This is the grid input to the "Fill pits function" to enforce stream flow directions.

    D8 Flow Directions

    Method: The flow direction from each grid cell to one of its adjacent or diagonal neighbors is calculated using steepest descent. The encoding is shown
    , i.e. 1 - East, 2 - North East, 3 - North, 4 - North West , 5 - West, 6 - South West, 7 - South, 8 - South East. Slope is evaluated in the direction of steepest descent and is reported as drop/distance, i.e. tan of the angle. Flow direction is reported as "no data" for any grid cell adjacent to the edge of the DEM domain, or adjacent to a no data value in the DEM. In flat areas flow directions are assigned away from higher ground and towards lower ground using the method of Garbrecht and Martz (1997). The flow path grid to enforce drainage along existing streams is an optional input, and if input, takes precedence over elevations for the setting of flow directions. The D8 flow direction algorithm may be applied to a DEM that has not had pits filled, but will then result in no data values for flow direction and slope associated with the lowest point of the pit.

    Dinf Flow Directions

    Method: The Dinf approach assigns a flow direction based on steepest slope on a triangular facet (Tarboton, 1997).

    Flow direction defined as steepest downward slope on planar triangular facets on a block centered grid.

    Flow direction is encoded as an angle 'ang' in radians counter-clockwise from east as a continuous (floating point) quantity between 0 and 2 pi. The flow direction angle is determined as the direction of the steepest downward slope on the eight triangular facets formed in a 3 x 3 grid cell window centered on the grid cell of interest. A block-centered representation is used with each elevation value taken to represent the elevation of the center of the corresponding grid cell. Eight planar triangular facets are formed between each grid cell and its eight neighbors. Each of these has a downslope vector which when drawn outwards from the center may be at an angle that lies within or outside the 45 o (pi/4 radian) angle range of the facet at the center point. If the slope vector angle is within the facet angle, it represents the steepest flow direction on that facet. If the slope vector angle is outside a facet, the steepest flow direction associated with that facet is taken along the steepest edge. The slope and flow direction associated with the grid cell is taken as the magnitude and direction of the steepest downslope vector from all eight facets. Slope is measured as drop/distance, i.e. tan of the slope angle. In the case where no slope vectors are positive (downslope), the flow direction is set using the method of Garbrecht and Martz (1997) for the determination of flow across flat areas. This makes flat areas drain away from high ground and towards low ground. The flow path grid to enforce drainage along existing streams is an optional input and if used, takes precedence over elevations for the setting of flow directions. The Dinf flow direction algorithm may be applied to a DEM that has not had pits filled, but will then result in no data values for Dinf flow direction and slope associated with the lowest point of the pit.

    D8 Contributing Area

    Method: Contributing area counted in terms of the number of grid cells (or summation of weights) is calculated using a recursive procedure due to (Mark, 1988). The contribution at each grid cell is taken as one (or from the weight grid when the optional weight grid input is used). The contributing area for each grid cell is taken as its own contribution plus the contribution from upslope neighbors that drain in to it. This is evaluated recursively starting from points (outlets) in the outlets shapefile, or when this is not input at each point in the grid. Starting the recursive evaluation at outlet points results in only the contributing area that drains to the designated outlets being evaluated.

    The contributing area programs check for edge contamination. This is defined as the possibility that a contributing area value may be underestimated due to grid cells outside of the domain not being counted. This occurs when drainage is inwards from the boundaries or areas with no data values for elevation. The algorithm recognizes this and reports no data for the contributing area. It is common to see streaks of no data values extending inwards from boundaries along flow paths that enter the domain at a boundary. This is the desired effect and indicates that contributing area for these grid cells is unknown due to it being dependent on terrain outside of the domain of data available. Edge contamination checking may be overridden in cases where you know this is not an issue or want to ignore these problems, if for example the DEM has been clipped along a watershed outline.

    Dinf Contributing Area

    Method: Contributing area counted in terms of the number of grid cells (or summation of weights) is calculated for the multiple flow direction Dinf approach using a recursive procedure that is an extension of the very efficient recursive algorithm for single directions (Mark, 1988). The contribution at each grid cell is taken initially as one (or from the weight grid when the optional weight grid input is used). The contributing area of each grid cell is then taken as its own contribution plus the contribution from upslope neighbors that have some fraction draining to it. The flow from each cell either all drains to one neighbor, if the angle falls along a cardinal (0, p /2, p , 3 p /2) or diagonal ( p /4, 3 p /4, 5 p /4, 7 p /4) direction, or is on an angle falling between the direct angle to two adjacent neighbors. In the latter case the flow is proportioned between these two neighbor pixels according to how close the flow direction angle is to the direct angle to those pixels, as illustrated in the Dinf flow direction figure above.

    Where no weight grid is used as input, the result is reported in terms of specific catchment area, the upslope area per unit contour length, taken here as the number of cells times grid cell size (cell area divided by cell size). This assumes that grid cell size is the effective contour length, in the definition of specific catchment area and does not distinguish any difference in contour length dependent upon the flow direction. Where a weight grid is used the result is reported directly as a summation of weights, without any scaling.

    Contributing area is evaluated recursively starting from points (outlets) in the outlets shapefile, or when this is not input at each point in the grid. Starting the recursive evaluation at outlet points results in only the contributing area that drains to the designated outlets being evaluated.

    The contributing area programs check for edge contamination. This is defined as the possibility that a contributing area value may be underestimated due to grid cells outside of the domain not being counted. This occurs when drainage is inwards from the boundaries or areas with no data values for elevation. The algorithm recognizes this and reports no data for the contributing area. It is common to see streaks of no data values extending inwards from boundaries along flow paths that enter the domain at a boundary. This is the desired effect and indicates that contributing area for these grid cells is unknown due to it being dependent on terrain outside of the domain of data available. Edge contamination checking may be overridden in cases where you know this is not an issue or want to ignore these problems, if for example the DEM has been clipped along a watershed outline.

    Grid Network Order and Flow Path Lengths

    • Strahler Network order grid (gord)
    • Longest Upslope Length grid (plen).
    • Total Upslope Length grid (tlen).

    Method: The D8 flow direction grid defines a grid network that extends to each grid cell. This function orders this network according to the Strahler ordering system. Cells that don't have any other grid cells draining in to them are order 1. When two (or more) flow paths of different order join the order of the downstream flow path is the order of the highest incoming flow path. When two (or more) flow paths of equal order join the downstream flow path is increased by 1. Algorithmically this is implemented as:
    Order = Max(Highest incoming flow path order, Second highest incoming flow path order + 1). This generalizes the common definition to cases where more than two flow paths join at a point.

    The longest upslope length is the length of the flow path from the furthest cell that drains to each cell. The total upslope path length is the length of the entire grid network upslope of each grid cell. Lengths are measured between cell centers taking into account cell size and whether the direction is adjacent or diagonal.

    Where the optional raster grid and threshold value are input, the function is evaluated only considering grid cells that lie in the domain with raster grid value greater than or equal to the threshold value. Source (first order) grid cells are taken as those that do not have any other grid cells from inside the domain draining in to them, and only when two of these flow paths join is order propagated according to the ordering rules. Lengths are also only evaluated counting paths within the domain greater than or equal to the threshold.

    River Network Raster Function

    This function is invoked by both the "Full River Network Raster" menu item under "Basic Grid Analysis" and "River Network Raster Upstream of Outlets" menu item inder "Network Delineation"

    • Pit Filled Elevation grid (fel).
    • D8 Flow Direction grid (p).
    • D8 Contributing Area grid (ad8).
    • Dinf Slope Grid (slp).
    • Dinf Contributing Area grid (sca).
    • Network Order grid (gord)
    • Longest Upslope Length grid (plen)
    • Verified Flow Path grid (fdrn)
    • Outlets Shapefile

    Methods: The output is a grid indicating streams, by the grid cell value 1 on streams and 0 off streams. Six stream delineation methods are provided each with their own parameters that may be adjusted on the interface. These methods are:

    • Use existing streams. To use this method existing streams need to have been "burned in" by using an enforced flow path grid (fdr) with "Fill Pits" and "D8 Flow Directions" functions. The stream network raster is then defined from these flow directions. No parameters are required.
    • DEM curvature based. The DEM is first smoothed by a kernel with the weights at the center, sides and diagonals as specified. The Peuker and Douglas (1975) method (also explained in Band, 1986) is then used to identify upwards curved grid cells. This method flags the entire grid, then examines in a single pass each quadrant of 4 grid cells and unflags the highest. The remaining flagged cells are deemed "upwards curved" and if viewed resemble a channel network, although sometimes lacking connectivity, or requiring thinning, issues that were discussed in detail by Band (1986). The thinning and connecting of these grid cells is achieved here by computing the contributing area using only these upwards curved cells. An accumulation threshold on the number of these cells is used to map the channel network.
    • Contributing area threshold. A threshold on the contributing area (in number of cells) computed by the D8 (suffix ad8) method is used to delineate streams.
    • Grid order threshold. A threshold on the network order grid (suffix gord) is used to delineate streams. This is the network pruning by order approach suggested by Peckham (1995) and used in RiverTools.
    • Area and slope threshold. A threshold is applied to the product A S y with the threshold and exponent specified. A is the Dinf specific catchment area (suffix sca) and S is the Dinf slope (suffix slp). This method was suggested by Montgomery and Dietrich (1992). (They used the exponent y = 2 and threshold C = 200 m in their study).
    • Area and length threshold. This is an experimental method that might be justified by searching for a departure from Hack's law. Streams are mapped as initiating when A > M L y . Here A is the D8 contributing area (suffix ad8) and L the longest upstream flowpath (suffix plen). In branching systems, Hack's law suggests that L = 1/M A 1/y with 1/y = 0.6 (or 0.56) (y about 1.7). In parallel flow systems L is proportional to A (y about 1). This method tries to identify the transition by using an exponent y somewhere inbetween (y about 1.3)

    If drop analysis is checked, then the threshold is searched between the lowest and highest values given, using the number of steps given on a log scale. For the science behind the drop analysis see Tarboton et al. (1991, 1992), Tarboton and Ames (2001). The smallest threshold in the set searched with absolute value of the t statistic less than 2 is selected. This is done automatically during the River Network Raster (Upstream of Outlets) step. The threshold selected is saved in the threshold variable so may be inspected afterwards by viewing the function interface following completion. Drop analysis is only possible when outlets have been specified, because if an entire grid domain is analyzed, as the threshold varies, shorter streams draining off the edge may not meeth the threshold criterion and be excluded from the analysis. This makes defining drainage density difficult and it is somewhat inconsistent to compare statistics evaluated over differing domains.

    Stream Order Grid and Network Function

    • Pit Filled Elevation grid (fel).
    • D8 Flow Direction grid (p).
    • D8 Contributing Area grid (ad8).
    • Stream Raster Grid (src).
    • Outlets Shapefile.

    Method: This function produces a vector network from the Stream Raster grid by tracing down from each source grid cell. The network topological connectivity is stored in the Stream Network Tree file, and coordinates and attributes from each grid cell along the network in the Stream Network Coordinates file. The Stream Raster grid is used as the source for the stream network, with Flow Direction Grid used to facilitate tracing down the flow paths. Elevations and Contributing Area are used to determine the Coordinate elevation and contributing area attributes. Points in the Outlets Shapefile are used to logically split stream reaches to facilitate representing watersheds upstream and downstream of monitoring points. The program has an option to trace downslope from outlet points until the stream reach network is first encountered to accommodate situations where gage locations are not accurately registered with respect to the stream network. The program looks for an attribute field "id" in the Outlets Shapefile and if found uses the numbers therein as identifiers in the Tree file.

    Stream Shapefile and Watershed Grid Function

    Method: This function translates the text file vector network representation in the Network Tree and Coordinates files into a Shape file. Further attributed are also evaluated. The subwatershed draining to each stream segment (reach) is delineated and an labeled with the value identifier that corresponds to WSNO in the Stream Shape File. The program has an option to delineate a single watershed comprising the entire area draining to the Stream Network.

    Watershed Grid to Shapefile Function

    Method: This function translates the grid watershed representation into a polygon shapefile. The Flow Direction grid is used to resolve ambiguities and prevent slivers of missing area between watersheds. Each shape in the Watershed Shapefile is assigned an identifier the same as the value identifier in the grid from which it was created, that in turn maps back to WSNO in the Stream Shape File.

    Drop Analysis Function

    The "Drop Analysis" command does not have any formal input parameters. Rather it takes its parameters from those set by the River Network Raster function. The following is displayed:

    • The threshold used in the network delineation algorithm.
    • The drainage density (inverse length units - typically m) of the resulting network.
    • Number of first order streams (Strahler ordering) in network with specified threshold.
    • Number of higher order streams (Sequential segments of the same order are counted as one Strahler stream) in network with specified threshold.
    • Mean drop (elevation difference between start and end) of first order streams.
    • Mean drop of higher order streams.
    • Standard deviation of first order stream drops.
    • Standard deviation of higher order stream drops.
    • Students t statistic for the difference between the first order and higher order mean stream drops.

    The procedure suggested in Tarboton et al. (1991, 1992) and Tarboton and Ames (2001) is to select the smallest threshold for which the absolute value of the t statistic is less than 2. This selects the highest resolution network consistent with the "constant drop law". In the display above the threshold of 14 would be selected. It is worthwhile to view the drop analysis because the threshold steps are sometimes quite coarse, so the user may want to use the River Network Raster function to change the number of intervals and lowest and highest values used in the search. Also it can occur that a large threshold is "correct" but a small threshold results in second or higher order streams also extending into the region that should not be streams. If these dominate the sample the following sequence of t statistics might result:
    <2, <2, <2, >2, >2, <2, <2, .
    The automatic procedure in these cases would pick the lowest threshold, but it is probably better (and this is admittedly subjective, unfortunately) to pick the 6th threshold. In making this judgement I feel that it is best to consider many things, like sample sizes (is the t statistic robust), the visual impression in comparison to contour crenulations, and the drainage density that would result from a Slope versus Area plot as discussed in Tarboton et al., (1991, 1992). The automated procedure is therefore not foolproof and some degree of judgement and subjectivity is required.

    Slope/Area (Wetness Index) function

    Method: This function takes the ratio Slope/(Specific Catchment area). This is algebraically related to the more common ln(a/tan beta) wetness index, but contributing area is in the denominator to avoid divide by 0 errors when slope is 0.

    Flow Distance to Streams function

    Method: This function computes the distance from each grid cell moving downstream until a stream grid cell as defined by the Stream Raster grid is encountered.

    Downslope Influence function

    Method: The Downslope Influence function (or influence zone) of a set y within the domain is
    I(xy) = A[i(xy)]
    where A[] is the weighted accumulation operator evaluated using the Dinf Contributing Area function. I(xy) says what the contribution from the set of points y is at each point x in the map. i(xy) is an indicator (1,0) function on the set y and I is evaluated using the weighted contributing area function only on points in the set y. The disturbance grid is used to encode the indicator function and must be 1 inside the zone y and 0 over the rest of the domain.

    This may be useful for example to track where sediment or contaminant moves.

    Upslope Dependence function

    Method: This function quantifies the amount a point x contributes to the point or zone y. It is the inverse of the downslope influence function
    D(xy) = I(yx)

    This is useful for example to track where a contaminant may come from.

    Decaying Accumulation function

    • Dinf Flow Direction grid (ang).
    • Decay Multiplier grid.
    • [Weight grid].
    • [Outlets Shapefile]

    Method: A decayed accumulation operator DA[.] takes as input a mass loading field m(x) expressed at each grid location as m(i, j) that is assumed to move with the flow field but is subject to first order decay in moving from cell to cell. The output is the accumulated mass at each location DA(x). The accumulation of m at each grid cell can be numerically evaluated.

    Here d(x) = d(i ,j) is a decay multiplier giving the fractional (first order) reduction in mass in moving from grid cell x to the next downslope cell. If travel (or residence) times t(x) associated with flow between cells are available d(x) may be evaluated as exp(- l t(x)) where l is a first order decay parameter. The weight grid is used to represent the mass loading m(x). If not specified this is taken as 1. If the outlets shapefile is used the function is only evaluated on that part of the domain that contributes flow to the locations given by the shapefile.

    Useful for a tracking contaminant or compound subject to decay or attenuation.

    Concentration Limited Accumulation function

    • Dinf Flow Direction grid (ang).
    • Weight grid
    • Decay Multiplier grid.
    • Indicator grid (dg).
    • Concentration threshold
    • [Outlets Shapefile]

    Method: This function applies to the situation where an unlimited supply of a substance is loaded into flow at a concentration or solubility threshold C sol over a set of points y. The indicator grid (dg) is used to delineate the area of the substance supply denoted using the (0,1) indicator function i(xy). A[] denotes the weighted accumulation operator evaluated using the Dinf Contributing Area function. The weight grid gives the loading for the supply to the flow (e.g. the excess rainfall if this is overland flow) denoted as w(x). The specific discharge is then given by
    Over the substance supply area concentration is at the threshold (the threshold is a saturation or solubility limit).
    If i(x y) = 1
    C(x) = C sol
    L(x) = C sol Q(x)
    Where L(x) denotes the load being carried by the flow.
    At remaining locations the load is determined by load accumulation and the concentration by dilution

    Here d(x) = d(i ,j) is a decay multiplier giving the fractional (first order) reduction in mass in moving from grid cell x to the next downslope cell. If travel (or residence) times t(x) associated with flow between cells are available d(x) may be evaluated as exp(- l t(x)) where l is a first order decay parameter. The Concentration grid output is C(x). The weighted accumulation grid output is Q(x). If the outlets shapefile is used the function is only evaluated on that part of the domain that contributes flow to the locations given by the shapefile.

    Useful for a tracking a contaminant released or partitioned to flow at a fixed threshold concentration.

    Transport Limited Accumulation function

    • Dinf Flow Direction grid (ang).
    • Supply grid (sup).
    • Transport Capacity grid (tc).
    • [Concentration in supply grid (cs).]
    • [Outlets Shapefile]
    • Transport Limited Accumulation grid (tla).
    • Deposition grid (tdep).
    • [Concentration grid (ctpt).]

    Method: This function applies to the situation where there is a supply of substance (e.g. erosion) and capacity for transport of the substance (e.g. sediment transport capacity). This function accumulates the substance flux subject to the rule that the transport out of any grid cell is the minimum of the transport in to that grid cell and the transport capacity. There is then deposition in the amount of the difference.

    Here E is the supply (sup) and T cap the transport capacity (tc). T out at each grid cell becomes T in for downslope grid cells and is reported as Transport limited accumulation (tla). D is deposition (tdep). The function provides the option to evaluate concentration of a compound (contaminant) adhered to the transported substance. This is evaluated as follows

    Where L in is the total incoming compound loading and C in and T in refer to the Concentration and Transport entering from each upslope grid cell.
    If then there is no erosion from the cell, so

    where C s is the concentration supplied locally and the difference in the second term on the right represents the additional supply from the local grid cell. Then

    C out at each grid cell comprises is the concentration grid output from this function. If the outlets shapefile is used the function is only evaluated on that part of the domain that contributes flow to the locations given by the shapefile. Transport limited accumulation is useful for modeling erosion and sediment delivery, including the spatial dependence of sediment delivery ratio and contaminant that adheres to sediment.

    Reverse Accumulation function

    Method: This works in a similar way to evaluation of weighted Contributing area, except that the accumulation is by propagating the weight loadings upslope along the reverse of the flow directions to accumulate the quantity of weight loading downslope from each grid cell. The function accumulates loading from the weight grid in excess of a weight threshold. The function also reports the maximum value of the weight loading downslope from each grid cell in the Maximum Downslope grid.

    This function is designed to evaluate and map the hazard due to activities that may have an effect downslope. The example is land management activities that increase runoff. Runoff is sometimes a trigger for landslides or debris flows, so the weight grid here could be taken as a terrain stability map. Then the reverse accumulation provides a measure of the amount of unstable terrain downslope from each grid cell, as an indicator of the danger of activities that may increase runoff, even though there may be no potential for any local impact.

    Grid Converter

    The GridConverter (Utilities) Dialog helps to convert one grid format to another. Input grid format can be an ESRI, ASCII or Binary and this function converts between these formats. Also the tool allows the user to convert the grid to different grid datatypes (Short,Integer or Float).

    Data formats and file naming conventions

    The programs are written to access following grid file formats.

    • ESRI proprietary grid format. The file is a "folder" on the computer comprising multiple parts.
    • ASCII format. A file with extension .asc.
    • Binary grid format. A file with extension .bgd.

    The ASCII grid format is that used by ESRI for export of files from ArcView and Arc/Info and comprises lines of header data followed by lists of cell values. The header data includes the following keywords and values:

    • ncols - number of columns in the data set.
    • nrows - number of rows in the data set.
    • xllcenter or xllcorner - x-coordinate of the center or lower-left corner of the lower-left cell.
    • yllcenter or yllcorner - y-coordinate of the center or lower-left corner of the lower-left cell.
    • cellsize - cell size for the data set.
    • nodata_value - value in the file assigned to cells whose value is unknown. This keyword and value is optional. The nodata_value defaults to -9999.

    For example an ASCII grid file looks like,

    The first row of data is at the top of the data set, moving from left to right. Cell values should be delimited by spaces. No carriage returns are necessary at the end of each row in the data set. The number of columns in the header is used to determine when a new row begins. The number of cell values must be equal to the number of rows times the number of columns.

    The ESRI and binary grid file formats are not defined here. They are only accessed through the grid object that is part of the software. This uses the gridio application programmers interface to access ESRI binary grids, and provides a standard interface for accessing all the grids.

    7. NEARPT3 (Nearest Points in E 3 )

    Nearpt3 is a pair of subroutines to find nearest points in E 3 . Preprocess takes a list of fixed points, and preprocesses them into a data structure. Query takes a query point, and returns the closest fixed point to it.

    Nearpt3 is very fast, very small, and can process very nonhomogeneous data. The largest example run on a laptop computer was St Matthew from the Stanford Digital Michelangelo Project Archive, with 184,088,599 fixed points.

    Preprocessing the David data set, with 28,158,109 fixed points, into a 723 3 grid took 0.6 microseconds per fixed point. Each query required 6 microseconds. The platform was an IBM T43p laptop computer with a 2GHz Intel Pentium M processor and 2GB of memory.

    Finding closed equipotential surfaces on a 3D grid

    In short, I'm looking for either: (1) Publications or other sources dealing with contour/isosurface finding algorithms, so that I can write my own implementation (and parallelize as best I can), or (2) Suggestions about high-level libraries in Fortran, C/C++, or other languages, that may be useful for finding such surfaces from a scalar 3D grid.

    I am currently analyzing data in checkpoint files from a very large hydrodynamics simulation (AMR w/ Lagrangian sink particles -- data in HDF5 format), and would like to find an algorithm for efficiently computing equipotential surfaces (a.k.a. countours, or isosurfaces) from the gravitational scalar potential grid. In particular, I need to be able to find the largest closed equipotential surface around a chosen local minimum (to within some acceptable step error), that may not enclose any other minima. This surface will be used as a selection criteria for further analysis of the volume formed by the enclosed cells. This link (particularly section 3) provides a description of the procedure I'm roughly attempting to replicate, though the IDL code provided takes 2D FITS maps of column density as input, and their data is much smaller:

    This seems to be very similar to the problem of computing and rendering isosurfaces in graphics (i.e. some variant of the marching cubes algorithm seems relevant), but since I'm interested in using the surface as a boundary for selecting all interior grid cells, rather than visualizing it, most of the libraries I'm familiar with seem ill suited (e.g. VTK). I would happily be proven wrong, however, if someone knows how I might use a graphics algorithm for this purpose -- maybe there's a way to find the isosurface with a graphics library routine and then simply extract information about the surface location or boundary cells?

    Any suggestions about how to approach this problem, or ideas about potentially useful libraries are much appreciated. If I write the algorithm more or less from scratch, I will likely use Fortran (easier for me to parallelize), but I would be happy to use C/C++ or other languages if they support libraries that are better suited to the task.

    Update: After a bit more research, I'm optimistic that a contour tree will be the perfect data structure for this problem. Here is a rather informal description. It seems to provide a wonderfully abstract graph representation of how the topology of the set of isosurfaces varies with isovalue.

    Local maxima and minima are represented as leaf nodes, and are said to create or terminate a "component" (graph edge). Interior vertices/nodes represent the joining or splitting of components, or the changing of genus, at critical points -- topological events, in other words. It is also possible to sweep through the space between the isosurfaces associated with two events, and assign all data points in this region to the component connecting the two events. This more thorough publication calls this an augmented contour tree.

    With this in mind, my current strategy is pretty straightforward: once an augmented contour tree has been produced, then the cells I want to associate with a local minimum are simply those associated with the component between the local minimum leaf node and the closest topological event vertex.


    The following flow chart summarizes the methodology developed for the large scale agrichemical transport in the Midwest rivers:

    Figure 4.1 Methodology of the large scale modeling of agrichemical concentrations in the Midwest rivers.

    4.1 Representative agricultural chemicals

    • nitrate and atrazine are representative of nutrients and herbicides, respectively,
    • nitrate and atrazine are present in measurable quantities in many Midwest streams, and
    • nitrate plus nitrite as nitrogen was the only nutrient measured during the USGS reconnaissance study (Scribner et al., 1993).

    It is assumed that the nitrate plus nitrite concentrations in Midwest streams are mainly derived from chemical fertilizers.

    4.1.1 Nitrate

    Nitrogen (N) in soils occurs as organic or inorganic N. The inorganic forms, include ammonium (NH 4 - ), nitrite (NO 2 - ), nitrate (NO 3 - ), nitrous oxide (N 2 O), nitric oxide (NO), and elemental N (N 2 ). The three most important forms, NH 4 - , NO 2 - ,and NO 3 - , usually represent 2 to 5% of the total soil N. The source of NH 4 + is from mineralization of organic N and from fertilizers. Some of NH 4 + is converted to NO 2 - , which is toxic to plant roots by bacteria Nitrosomonas (2NH 4 - + 3O 2 = 2NO 2 - + 2H 2 O + 4H + ), and then oxidized to NO 3 - by Nitrobacter (2NO 2 - + O 2 = 2NO 3 - ). The NO 3 - anion is very mobile and subject to leaching losses (Tisdale et al., 1993).

    Nitrate in streams is derived from many anthropogenic and natural resources including chemical fertilizers, animal wastes, domestic sewage, legumes, mineralization of vegetation, soil organic matter, and from atmosphere through electrical, combustion and industrial processes. NO 3 - is a very soluble and mobile anion. It can be transported from agricultural fields in both overland flow and subsurface flow, and by volatilization into the atmosphere. Ammonium is adsorbed by the soil colloids and moves very little until converted to NO 3 -

    In contrast to the overland transport in which nitrate takes minutes or hours to get to a stream, downward vertical leaching and the subsequent underground travel is a long process which takes months or years. The soil system has a strong memory with respect to nitrate production and leaching. Jones and Burt (1993) presented a study in which 64% of annual nitrate concentration in streams was explained by a stepwise regression involving the year of measurement, and each of the previous two years. It may take nitrate years or even decades to appear in rivers as a base flow pollutant.

    There are some losses of nitrate due to erosion, but for humid temperate climates, erosion is generally an insignificant process compared with leaching and runoff (OECD, 1986). Other authors indicate that the adsorption has no marked influence on the rate of NO 3 - movement (Keeney, 1983, Bailey and Swank 1983). The predominant losses of nitrate in an agricultural field are due to assimilation by row crops and by other terrestrial and aquatic plants.

    • Temporal and spatial distribution of the rainfall. Because of N mobility in soils, the greater the surplus rainfall, the grater the possibility of loss of N through leaching.
    • Temporal and spatial distribution of the temperature. Since higher temperatures enhance nitrification, ammonical N applied before planting is more subject to nitrification and leaching.
    • Technical factors. In late winter the ground may be to wet for machinery to be operated and spring application is usually too late for small grain to respond in yield to the applied nitrogen fertilizer.

    The time of the nitrogen fertilizers application as well as the decrease in temperature during late fall, winter, and early spring, cause high nitrate plus nitrite concentrations in surface waters. These concentrations decrease in late spring and summer when the plant demands for nutrients are high (Goolsby and Battaglin, 1993, Davis and Keller, 1983). In addition, a significant decrease in nitrate concentration in surface water may result from assimilation of nitrate by algae and by instream riparian macrophytes (Heathwaite, 1993, Moore, 1991) as well as nitrate may be converted by the denitrification bacteria and various chemical processes into free nitrogen and nitrogen oxides which escape into the atmosphere. Since these processes are stimulated by high temperatures and low flow rates, the highest loss of nitrate in lakes and rivers occurs during the summer. Lakes and reservoirs act as a "buffer," thus they are less responsive to seasonal changes then rivers are (OECD, 1986).

    Denitrification and other processes of biochemical degradation of nitrate can be modeled as a first order reaction with an "overall" decay rate dependent on temperature and carbon content. Some losses may result from infiltration (river water seepage into groundwater).

    4.1.2 Atrazine

    Atrazine is a herbicide that controls broadleaf weeds in fields of corn and sorghum. It is one of the most widely used herbicides in the United states (Comfort and Roeth, 1996). The EPA set the drinking water health limit for atrazine at 3 µg/L (ppb). The conventional water treatment does not remove this herbicide. Recent studies conducted in 29 communities throughout the Midwest, Louisiana, and the Chesapeake Bay detected high concentrations of atrazine in tap water in months from May through July, some of them exceeding EPA Maximum Contamination Level (EWG, 1996).

    Numerous laboratory tests as well as field studies have been performed to determine the behavior of this herbicide in different chemical and physical environments for almost half of the century. Some of the published parameters are as follows: Atrazine (2-chloro-4-ethylamino-6-isopropylamino-s-triazine) is a low solubility herbicide its water solubility in typical temperature and pH varies from 30 to 35 ppm (mg/L). The volatility of this chemical is very low ( vapor pressure = 0.3*10 -6 mm Hg = 0.00004 Pa). Published values of the octanol extrability coefficient (soil sorption coefficient) Koc are from 130 to 172. The octanol-water partition coefficient Kow equals 251 (e.g., Weber, 1972, Hamaker, 1975, Wauchope, 1978, Rao, et al., 1983, Weber, 1988, work cited by Plimmer, 1988).

    There are two parameters that characterize the chemical decay in soil: half-life and chemical persistence. Kruger et al. (1993) reported half-life of atrazine in soil under unsaturated conditions ranged from 41 d to 231 d, whereas in saturated soil at the 90-120 cm depth, the half-life was 87 d. Similar values were cited by Goring et al. (1975). Wauchope (1978) determined that atrazine persistence in soil, that is the approximate time for 90% disappearance from soil, is 12 months and may vary from 6 to 18 months, that corresponds to half-life from 55 d to 165 d, depending on climate and soil.

    Atrazine decays into many degradation products which can be more persistent and mobile than their parent compound. During a USGS 1991 Midcontinent survey of near-surface aquifers, deethylatrazine, an atrazine metabolite, was the most frequently detected compound followed by atrazine, and then deisopropylatrazine, another atrazine metabolite (Kolpin et al., 1991). The half life of the deisopropylatrazine is much longer in surface water than in soil (Goolsby et al., 1993 pp. 51-62, Comfort and Roeth, 1996). However, Kruger et. al. (1993) found that deisopropylatrazine may be even less persistent under saturated conditions then in saturated soil. Their estimations of the deisopropylatrazine half-life ranged from 32 d to 173 d in unsaturated soil at top 30 cm, and from 58 d to 173 d in saturated soil at the 90 to 120 cm depth.

    Weber (1988) presented results of research, which he conducted with J. A. Best, on the effect of pH on the dissipation of atrazine applied to soil. No parent component volatilization was detected 0.1%-0.2% was found in leachate, plants used two to four percent. Ninety percent of atrazine was retained in the soil layer.

    The adsorption and movement of s-triazines in soil depends upon such factors as soil organic matter, clay minerals, pH, temperature, soil moisture, concentration and species of other ions in the system (Weber 1972, Weber 1988, Goring et al. 1975). In CREAMS model, the pesticide in runoff is partitioned between the solution phase and the sediment phase (Knisel et al., 1983).

    Study by the USGS of the occurrence of herbicides in precipitation in the Midwest and Northeastern United States showed that significant amounts of atrazine were lost through volatilization and subsequently returned to the land through precipitation washout. The amount of atrazine in precipitation washout was equal to one half the loading found in the Mississippi River (Goolsby et al., 1993 pp. 75-86). Plimmer (1988) discussed reports which described the identification of atrazine in fog. These findings were apparently in contradiction to the findings about negligible volatilization published by Weber (1972, 1988). As was pointed out by Kenaga (1975), vapor pressure is useful where molecules are volatilizing from a liquid or solid mass of the same compound. When a pesticide is applied to heterogeneous surfaces such as natural water, soil, foliage, wood, or glass, volatility is variable because sorption varies. There is a possibility that atrazine enters the atmosphere adsorbed on particulate matter through dust blowing from the land surface. The process of atrazine volatility is not well understood.

    4.2 Selection of analysis region and map coordinate system

    The Mississippi - Missouri basin above Thebes, Illinois (drainage area 2.3*10 6 km² ) together with the Ohio River basin above Grand Chain, Illinois (drainage area 0.5*10 6 km²) constitute the primary region that is used for estimation of statistical model parameters. Its extend is determined by the USGS reconnaissance study of selected herbicides and nitrate in Midwestern United States (Scribner et al., 1993). This region is one of the most extensive agricultural areas in the country, producing over 80% of all US corn and soybeans (Oberle and Burkart, 1994).

    Figure 4.2 The Mississippi-Missouri and Ohio River Basins.

    Since processing data for an area that covers almost three million square kilometers requires a computer with adequate operational and storage memory, the final model is built, and verified on the selected subregion, i.e., the Iowa-Cedar River watershed in Iowa. There are two sites for which measurement of agricultural chemicals are available: Old Man's Creek near Iowa City and Cedar River at Palisades. Flow rate is recorded in as many as 36 USGS gauging stations. Diversity of such geographic features as lakes, reservoirs, wetlands, hills, plains, and a wide range of stream sizes constitutes a very good representation of the majority of the Midwest. Figure 4.3 shows the Iowa and Cedar Rivers, and the location of gauging stations.

    Figure 4.3 The Iowa River with tributaries and the USGS gauging stations.

    All maps utilized in this research are represented in Albers Conical Equal Area Projection, generally accepted for large maps of the USA. This projection has the valuable property of equal area representation, combined with a scale error that is practically the minimum attainable in any system covering such a large area in a single sheet.

    The following parameters are applied (standard for USA): units = meters, first standard parallel = 29䣆'00'', second standard parallel = 45䣆'00'', latitude of projection's origin = 23䢨'00'', false easting = 0.000 m, false northing = 0.000 m, longitude of central meridian = -96䢨'00''. Using one common projection for the whole Midwest eliminates the problems associated with merging separately modeling regions into one unit.

    The Albers projection is of the conical type, in which the meridians are straight lines meeting in common point beyond the limits of the map, and the parallels are concentric circles, the center of which is at the point of intersection of the meridians. The meridians and the parallels intersect at right angles and the arcs of longitude along any given parallel are of equal length. The spheroid is intersected by a cone at two parallels known as the standard parallels for the area to be represented. On the two standard parallels, arcs of longitude are represented in their true lengths, or in their exact scale. Between the standard parallels, the scale along the meridians is too large and beyond them too small (Deetz and Adams 1969).

    4.3 Mathematical description

    4.3.1 Overview of transport equations

    There are two basic mechanisms that are responsible for the transport of dissolved and suspended solutes in surface waters: advection and diffusion/dispersion. These two processes are described by the advection-dispersion equation which is a fundamental equation for majority of the pollutant transport models. Equation (4.1) describes one-dimensional advection-dispersion in the reach with uniform cross-sectional area.

    c = concentration of pollutant [g/m³]
    t = time [s]
    x = distance along the river [m]
    = mass flux due to the longitudinal dispersion [g/m²s]
    E x = longitudinal dispersion coefficient [m²/s]
    - vc = mass flux due to advection [g/m²s]
    S i = i -th source/sink of the constituent [g/m³s]
    i = 1. n , n = number of sources/sinks
    K j = decay rate due to the j -th process [1/s]
    v = flow velocity [m/s]

    The source/sink term, S i and reaction term, K j c represent a wide range of features such as lateral flux, transient storage, biotic and abiotic retention, benthic flux, periphyton retention, sediment retention, reaeration, photosynthesis, and nitrification (U. of Mississippi, 1990, James and Elliot, 1993).

    • Numerical solution of the advection-dispersion equation requires subdivision of the time domain into relatively short intervals (minutes, hours), but this is less useful when calculating monthly means of chemical transport in extensive stream networks over large areas
    • It is difficult to write a procedure which solves the dispersion-advection equation using a GIS script language, so such a model has to be solved external to the GIS.

    Lagrangian type transport models, such as a Moving Segment Model (MSM), can be efficiently incorporated into GIS. In MSM (James and Elliot, 1993) the stream is subdivided into series of segments. Within each segment the variations of chemical concentration are calculated by summing, for example, hourly changes due to all the processes involved. The process within the block is described by the following equation:

    where: S i is the i -th source/sink and K j is the j -th reaction coefficient.

    • Sorption and desorption between dissolved and particulate component in the sediment and water column
    • Settling and resuspension of particles
    • Diffusive exchange between the sediment and water column
    • Loss and gain of the chemical due to the chemical and biochemical reactions such as biodegradation, volatilization, photolysis
    • Advective and diffusive transport of the chemical in water and as a bed transport
    • Net deposition and loss of chemical to deep sediments.

    Figure 4.4 shows these major reaction mechanisms and transfer routes of chemicals and solids in both river water and the river bed.

    Figure 4.4 Reactions and transfers in a natural water system.

    The processes is surface waters can be described by the mass balance equations for the dissolved and particulate components (Eq. 4.3 and Eq 4.4). For a specific chemical, e.g. nitrate or atrazine, some processes must be included in the transport model whereas some processes are not significant and therefore can be neglected.

    K v , K f , K s , K u = the bulk transfer coefficients of volatilization, dissolved
    exchange, settling and scour respectively
    c s , c p = dissolved and particulate concentrations in water, respectively
    c b , c bp = dissolved and particulate concentrations in bed
    K c , K p = decay coefficients of dissolved and particulate form
    K 1 = K 0 r c m K 0 = adsorption coefficient, r c = adsorptive capacity,
    m = concentration of solids
    K 2 = desorption coefficient
    K q = coefficient of dilution due to the groundwater inflow.

    If no dispersion is assumed ( E x = 0), no sources and sinks exist ( S i = 0), and the chemical is undergoing first order decay, i.e., the losses can be modeled by the first order reaction ( K j c = K d c ), the solution of the Eq. (4.1) is:

    This equation can be easily incorporated into GIS to model the transport agrichemicals in surface water.

    4.3.2 Regression equation development

    There is no complete and consistent data base of parameters required for comprehensive modeling the agrichemical runoff from the field and its transport in surface waters of the Midwest. The runoff and transport parameters have to be estimated from such data as: (1) the observed flow rate and the measured agrichemical concentrations in the variety of watersheds scattered over the upper Mississippi-Missouri and Ohio River basins and (2) the DEM from which a watershed morphometry can be calculated. Only a model based on the statistical analysis comes into consideration.

    Commonly used by the USGS equation for event mean concentrations of eleven water quality parameters as well as the total storm event loads has the following form (Huber, 1993):

    c = mean concentration
    X 1 . X n = explanatory variables
    ß 0 . ß n = regression coefficients
    BCF = bias correction factor.

    Cohn et al. (1992) used similar equation to describe the concentrations of several nutrient species entering Chesapeake Bay through its tributaries. Their equation (Eq. 4.7) includes periodic functions of a time variable to explain cyclic behavior:

    where Q' is the normalized flow rate and T' is the normalized time variable (yr).

    The following Section introduce a set of explanatory variables that have potential application in the large scale model of agrichemical runoff from the field and transport in the river network.

    4.3.3 Agrichemical runoff from the field

    Agrichemical application . Since it is a common practice to put more chemical on the field than the amount that can be completely utilized by the vegetation and the chemical - microbiological processes, some of the nutrients and herbicides runoff from the field. Many studies proved that the more chemical is applied on field U , the higher chemical loads are carried by a river (e.g. Battaglin et al., 1993).

    Land slope . The higher the slope of the land S L , the chemicals are more susceptible to the washout, therefore the greater losses of the chemical from the field can be expected. This complex process is mostly influenced by the gravitational forces. Since the vegetation has a smaller chance to uptake nutrient or herbicide when the land slope is higher (due to the shorter residence time), a greater portion of the mass applied on the field will reach the stream network. All models of erosion and sediment transport comprise the slope parameter.

    Distance to a stream . The influence of the average distance L L between the point of herbicide or nutrient application and the location of the stream network on the amount of mass that enters the stream is not as obvious as in the case of previously discussed watershed slope. Intuitively it may be expected that since the longer the average travel path is, the larger the chemical losses are and thus the coefficient ß 3, from equation ln(c) = . ß 3 * ln( L L ), estimated by the regression analysis should be negative. But, since the predictor variable L L represents complex watershed features, the coefficient ß 3 does not necessarily have to be negative. The average length of the land-flow path highly depends on the stream density, i.e., on the length of the streams per unit area of watershed. The more streams there are within a given watershed, the shorter the length of the land flow. The density of the stream network influences the amount of chemical that enters a unit length of the stream. Therefore, for watersheds characterized by the higher value of L L , the amount of chemical that enters a stream (mass per unit length) is higher than the amount reaching the stream network which is located in a watershed with smaller L L .

    Climate . Both, the precipitation and the temperature have an impact on the vegetation growth. Moreover, the greater the water surplus (precipitation minus the potential evaporation), the greater the possibility of loss of agrichemical through leaching if the crop is not growing vigorously or through washout if the land is not protected by a plant cover. Denitrification depends on the amount of water in soil, whereas nitrification rate is highly correlated with the temperature (Tisdale et al., 1993).

    The influence of the climate on the concentrations in surface waters can be illustrated by comparison of studies performed in different climatic regions. For example, the largest number of atrazine detections in Swedish stream waters was in July, August and September (Kreuger and Brink 1988), whereas in the Midwestern United State the major atrazine runoff occurs in May and June (Scribner et al., 1994). Swedish vegetation period is short (6-8 months) and cold (3º-17º C)

    Since the region under scrutiny extends from about 37º N to 50º N (latitude) and from 79º W to 114º W (longitude) the spatial distribution of the average temperature and the average precipitation influences not only the spatial distribution of the chemical runoff from the field but it also causes the spatially different agrichemical application time. It is believed that difference in climate conditions within studied region can be represented by a spatially distributed adjustment coefficient as well as a time shift introduced into periodic functions that describe seasonal variation of agrichemical concentration in the surface water.

    4.3.4 Transport in rivers

    The losses of agrichemical in stream can be described by an exponential function of the travel time (the traditional approach) or by an exponential function of the travel distance (Smith. et al., 1993). Although it is possible to estimate the travel time from the observed flow time series, in this research the agrichemical decay has been related to the travel distance, i.e., the amount of chemical that enters the stream in point i (cell i in raster representation of river) decays as it travels downstream according to the following equation:

    R k = total mass runoff from k -th watershed
    R i = mass that enters the stream in point i (cell i )
    n = number of all cells that constitute the stream network located within k -th watershed
    k S = overall distance decay coefficient
    L ik = length of the flow path from stream cell i in which the chemical runoff from a field enters the stream network to the watershed outlet located in cell k [102 km]
    R k = total mass runoff from k -th watershed.

    The travel distance L ik can be estimated from the DEM using GIS functions. The average regional distance decay coefficient k S may be calculated by a trial and error method which goal is to eliminate a regression coefficient ß 4 . Symbolically it can be described as follows:

    1. Assume k S = 1
    2. Using GIS calculate for each watershed k for which the concentrations were measured
    3. Estimate regression coefficient ß 4 of the concentration equation
    4. Assume different value of decay coefficient k S
    5. Repeat steps from (2) to (4) until the coefficient ß 4 is close enough to 1.

    Since the method outlined above is very computation intensive, it hasn't been fully utilized in this research. The total mass runoff R k from k -th watershed has been replaced by an explanatory variable E s defined by formula:

    The equation (4.10) assumes that the unit mass of chemical enters the stream network and it is uniformly distributed over all stream cells.

    Another variable that has potential to explain the transport of agrichemical in rivers is the stream slope S S .

    4.3.5 Seasonal variations

    The seasonal variations of the nitrate concentration as well as the atrazine concentration in surface waters are modeled by a periodic function:

    SW = seasonal component
    sin(2 n¶m /12), and cos(2 n¶m /12) = components of the Fourier series (the cycle corresponding to n = 1 has a 12-month period, n = 2, 3, 4, 5, and 6 are harmonics of period 12/ m months
    m = month of a year (1 for January)
    ß n,6 , ß n,7 , . = regression coefficients.

    Research performed by Walling and Webb (1984, work discussed by Jones and Burt, 1993) shows that in most streams the annual nitrate variations do not have a perfectly symmetrical sinusoidal form the annual minimum and maximum occurred 4-6 weeks later and 2-3 weeks earlier then the timing suggested by a single harmonic. Furthermore, no clear seasonal variations of nitrate concentration exist in catchments where groundwater with a consistently high nitrate concentration mixes with quick flow of varying origin and concentration (Jones and Burt, 1993).

    4.3.6 Exemplary equations for chemical concentration and load

    To make the presentation of methodology complete this Section presents two exemplary equations that describe the agrichemical concentration (Eq. 4.12) and load (Eq. 4.13) in the Midwest rivers:

    L = constituent load ( L = Qc ) [g/s]
    Q = flow rate at a given stream location [m³/s]
    c = constituent concentration [g/m³]
    U = annual chemical use in drainage area above the sampling point [kg/yr]
    U R = annual chemical application rate [kg/km²/yr]
    S L = average slope of the land [dimensionless]
    L L = average length of the constituent travel path, from the point of application to the stream network within a given watershed to the sampling point [km]
    A = drainage area [km²]
    E S = average of the exponent of negative flow distance in streams
    n c = number of cells that constitute the stream network within sampled watershed
    L ik = length of the flow path from the i -th stream cell to the watershed
    outlet k [102 km]
    S S = average slope of the stream network [dimensionless]
    sin(2 n¶m /12), and cos(2 n¶m /12) = components of the Fourier series (the cycle corresponding to n = 1 has a 12-month period, n = 2, 3, 4, 5, and 6 are harmonics of period 12/ m months
    m = month of a year (1 for January)
    ß 0 , ß 1 , ß 2 . = regression coefficients.

      describe changes in mass applied U (or U R ) on a field as it travels from the point of application to a stream

    4.3.7 Extracting values of explanatory variables for the regression analysis

    The estimation of values of explanatory variables such as agrichemical application rate, total application, average stream slope, average land slope, exponent of the negative flow length, and average distance from the field to the closest stream is performed in two steps: (1) Grids of spatially distributed parameters are constructed for the Upper Mississippi - Missouri River and the Ohio River basins, and (2) for each cell that represents the sampled watershed outlet the value of the explanatory variable is extracted.

    Here, the grid of spatially distributed parameters means a grid which each cell contains average or sum calculated for the total drainage area upstream to the given cell. Figure 4.5 shows selected cells that contain a value that characterize the upstream watershed.

    Figure 4.5 Example for explanation of the grid of spatially distributed values of explanatory variables.

    4.3.8 Application of the regression models

    • Estimation of monthly and annual loads in rivers of the Upper Mississippi-Missouri basin, including the Ohio River basin. The calculations are made using grid-map algebra. The loads are estimated in all cells (9.6 *10 6 cells of size 500 m ) that constitute the basin. The equations that describe chemical load in rivers require only a map of distribution of the total agricultural chemical use and the parameters of watersheds that can be easily determined from digital elevation model.
    • Calculation of agrochemical concentrations in the Iowa-Cedar River basin. The map of flow rate is required in obtaining the spatial image of atrazine and nitrate concentrations in surface waters. The estimations are performed for the watershed subdivided into subwatersheds of average area 31.6 km² utilizing maps in a vector format. The regression equations are applied to each subwatershed utilizing the data from the attribute table that characterize the upstream drainage area. A GIS approach to model the concentrations in rivers is discussed in the following Section 4.4. The procedure of the spatial distribution of the historical flow record is presented in Section 4.5.

    4.4 GIS model description

    4.4.1 GIS and cascade modeling

    The GIS technology gives the opportunity to construct versatile "cascade" models. The idea of cascade modeling incorporated here into GIS has been extracted from the methodology used in forecasting the municipal water use (Maidment and Parzen, 1984, Mizgalewicz, 1991).The GIS cascade structure is two dimensional. It can be applied in spatial dimension as well as in a temporal one. In the time domain, a general model describes spatial distribution of average annual amounts of agricultural chemicals in rivers. The annual values are then broken into seasonal or monthly values. For example, the annual average concentrations c(y) for the Midwest can be estimated by the general function:

    where: X is a vector of explanatory variables such as annual agrichemical application and watershed morphometry (area, land slope, stream slope, stream length, overland flow length) which practically do not change in time.

    The Equation (4.14) can be further extended by adding a seasonal component, a monthly fractions S(m) to brake down the annual predictions of concentration c(y) into monthly values c(m):

    Cascade modeling enables further extending of the Eq. (4.15) by adding such elements as trend and irregular component.

    Cascade modeling in spatial domain implies that the modeling process is subdivided into several levels of resolution, i.e., the Upper Mississippi-Missouri and Ohio River can be subdivided into three basins: the Missouri River Basin above junction with the Mississippi River, the Upper Mississippi River Basin, and the Ohio River Basin. The transport of agricultural pollutants in the Upper Mississippi River is estimated by utilizing results of sub-models that describe transport in individual component basins such as the Des-Moines River, the Skunk River, the Iowa-Cedar River, and Wisconsin River Basins.

    Additional spatial subdivision of the Midwest can be performed by introduction of climate zones. Each zone is defined, for instance, by specific range of temperatures and precipitation depths. Thus, the Eq. (4.15) would have the following form:

    where: c(m,g) is the monthly agrichemical concentration in rivers modified for the climate zones g( T,P,Z ) is a function of T - zonal temperature, P - zonal precipitation depth, and Z - zone location.

    The diagram shown in Figure 4.6 illustrates the example of spatio-temporal cascade modeling within GIS.

    Figure 4.6 Example of the cascade modeling within the GIS.

    Parameters of the mathematical description of the chemical applicationrunoff process and chemical losses in streams are stored in a GIS attribute table. Moreover, the equations are stored as objects in a database table. Storing equations as a database objects not only permits to fully implement the cascade modeling technique into the GIS but also it simplifies the structure of the model, allows one to test and compare unlimited number of equations, and permits changes to the mathematical description by simple record editing operations. The model prototype is constructed within ArcView, a GIS application (ESRI, 1995).

    4.4.2 Subdivision of study region into modeling units

    Elementary watersheds, i.e., modeling units that are considered lumped systems, constitute the smallest units into which the region is subdivided. Each unit is characterized by set of parameters such as area, agricultural chemical application, slope, depth of precipitation, water runoff, average elevation, and an equation that relates the mass of chemical applied and the mass of chemical runoff.

    Since the GIS offers very convenient tools for data storage and manipulation, all attributes can be stored and extracted by models that operate on different spatial scales. The order of processing the individual watersheds is the researcher's choice. Moreover, each watershed does not need to be divided into modeling units of the same size. To make the modeling process efficient, such features as density of spatial information and diversity of terrain should influence assumed size of the elementary unit.

    1. Points in which the drainage exceeds a threshold value. In these points streams originate. Here the threshold value of 25 km² has been assumed
    2. Points located immediately upstream of the stream junction and
    3. Gauging station sites.

Figure 4.7 shows example of selected watershed outlets and corresponding modeling units:

Figure 4.7 Example of watershed outlets: (1) beginning of the river, (2) stream junction, and (3) gauging station.

The watershed structure, developed by Maidment (1993) for hydrologic modeling utilizes type(2) and (3) watershed outlets which are positioned at the stream junctions and at the gauging station locations. In this study this set of watershed outlets has been extended by adding the points in which the stream network, delineated from DEM, begins (type 1 outlet). The test have been performed to determine the influence of additional unit watersheds on the uniformity of region subdivision and the control on the unit area. The Iowa-Cedar River basin has been subdivided using two sets of outlets. Full set: Median = 25.8 km² and mean = 31.6 km², reduced set (type 1 outlets not included): median = 37.9 km² and mean = 46.7 km². By including the outlets of type 1, the median of unit watershed drainage areas is very close to the used threshold. Figure 4.8 compares the frequencies of modeling unit area.

Figure 4.8 Comparison of the frequency of modeling unit area subdivided using different sets of watershed outlets.

  • The watershed is subdivided into more uniform (similar area) modeling units. one has mode control on the area of units
  • It makes possible to determine the flow and the constituent load in all nodes of the stream network. Each node has defined a contributing area and
  • The representation of all streams can be standardized. Each reach has input (inflow or load that enters the reach at the upstream end), lateral loss or gain, and output (outflow or load that leaves the reach at its end).

The threshold area of 25 km² has been selected after many tests with different values were performed. Smaller threshold area resulted in very dense stream network and thus large number of areas of size of 1-2 cells. After Arc/Info conversion of these small units from raster format (grid) into vector format many of them disappeared. In addition, smaller than 25 km² values are not justified by the data used in the research. For example, the agrichemical application rate is published with the county-size resolution. Stream network delineated from DEM, is slightly more dense that the one represented by the RF1 (1:500,000 digital map of rivers).

The larger threshold value than 25 km² resulted in very coarse subdivision of studied region and a low density stream network. Figure 4.9 presents the Iowa-Cedar River watershed divided into modeling units of different sizes.

Figure 4.9 Division of the Iowa-Cedar River into modeling units using different threshold drainage area: 25 km², 400 km², and 2500 km².

4.4.3 Unit watershed flow system

Network flow system is common function in vector GIS (Maidment 1993). Here the flow system has been extended: instead of arcs, unit watersheds (polygons) compose the flow system. The flow topology is described by two numbers: the modeling unit ID and the ID of the downstream unit, i.e., the next unit on the flow path. The ID = 0 of the next unit indicates that there are no more downstream units. Figure 4.10 shows an example of the description of modeling units flow topology.

Figure 4.10 An example of the flow topology of the unit watersheds ( x,y ), x is the unit ID, y is the ID of downstream to x unit.

The flow direction indicator is the basic concept in the raster-based hydrologic modeling. For example, Arc/Info Grid denotes the next cell in the flow path by one of eight numbers: 1 represent flow into E (East) neighbor cell, 2 in SE cell, 4 in S cell, 8 in SW cell, 16 in W cell, 32 in NW cell, 64 in N cell, and 128 in NE cell. These numbers are called "flow direction". Since modeling units are not regular spatial shapes, as the cells are, it is not possible to create an uniform numbering system of the flow direction. Proposed here method describes the flow connectivity by specifying the ID of the next unit on the flow path.

Addition of the item that describes the flow direction into the unit watershed attribute table made possible to utilize most of the concepts of hydrologic modeling, as far used in the raster GIS, into the vector environment. Such functions as flow accumulation, basin delineation, flow length can be applied for any shapes and thus for unit drainage areas. The flow system described here may be adopted for all models in which the conditions in a modeling unit does not influence the conditions in upstream unit such as kinematic wave routing and constituent decay as it flows downstream.

  • weighted average of a feature (for example agrichemical application ) for the total drainage area upstream of each (or selected) modeling unit
  • accumulated value over total upstream drainage area going along the flow path (vector version of the Grid flowaccumulation function) and
  • difference between the inputs and the output for each modeling unit, "flow-un-accumulation".

The unit watershed approach of modeling agrichemical transport can be considered as a semantic data model:

. semantic data modeling --- creating abstractions of geographic data layers which map one to one with their geographic representation but which are simplified in a functional description to the level needed for hydrologic modelling. (Maidment, 1993)

Figure 4.11 shows examples of the conceptual stream network. The links have been developed by connecting the cells immediately below each modeling unit outlet. Thus, Although the system is conceptual, its each node is located on the river represented in Grid format.

Figure 4.11 Flow system in the Iowa-Cedar River basin subdivided into modeling units of different sizes (threshold drainage area: 25 km², 400 km², and 2500 km²).

In addition, the arc attribute table can contain wide range of items that describe the links (length, slope, RF1-ID, agrichemical decay coefficient, travel time) as well as the parameters of the beginning node and the ending node (elevation, coordinates, and length from the watershed outlet). In this research only the conceptual stream network represented by unit watersheds is utilized.

4.4.4 Ordering system of the modeling units

To enhance the calculations performed by the GIS model, the following system of ordering modeling units has been utilized: All the most upstream units are assigned an order one. An interior unit has the order equal to the maximum order of the upstream units increased by one. Figure 4.12 compares this ordering method with the Strahler and the Shreve ordering systems (ESRI, 1992). The proposed ordering system allows to perform calculations in consecutive manner: initially all the first order units are processed, then the order is increased by one and all units of order two are evaluated. The routing is performed until the unit of the highest order is calculated.

Figure 4.12 Comparison of the stream ordering systems: (a) Strahler, (b) Shreve, and (c) utilized in the agrichemical transport model.

4.4.5 Enhancement of the stream delineation process

  • Since the location of a stream that is delineated from elevation data depends on the cell size, the stream networks determined from the DEMs of different resolutions are not compatible. Thus, the gauging stations linked to one gridded river system will not be in agreement with the other systems derived from DEM grids of different cell sizes.
  • The stream system constitutes the best framework for the spatial flow. It took hundreds and thousands of years for a river bed to develop to its current form. Since the river practically does not change, the other information such as location of gauging stations may be related to the location of stream reach.
  • In the flat regions, the streams delineated from the DEM tend to be straight lines, whereas, in reality the rivers have a tendency to meander. This causes an overestimation of the stream slope and an underestimation of the river length.
  • RF1 represents the true river system, whereas, the stream network delineated from the map of elevations just approximates the same system.
  • Since the elevations of the DEM are changed, the modified DEM can not be used for other tasks than the flow direction estimation.
  • The stream network has to be represented by a single lines. The double line description of rivers can be utilized as far as the distance between the lines is smaller that the cell size.
  • All existing loops have to be removed or opened to eliminate the ambiguous flow paths.
  • All lakes have to be converted into line representations or to polygons.
  • The cell width applied for adjustment process should be smaller than the half of the distance between any streams in RF1, to avoid connections of stream networks from different basins, that may be created when converting from vector format into raster one.

If the river network does not fulfill the above mentioned requirements, some editing after converting into grid format is necessary.

Incorporating the RF1 into grid has an additional advantage. By assigning the reach ID from RF1 to the raster river representation, the attribute table of RF1 can be linked with the attribute table of derived grids. Thus such information as the average flow velocity, average flow rate or stream names that are in RF1 attribute table can be used for grid models and vice versa, the parameters estimated in grid such as reach slope, drainage area and flow length can be assigned to the streams in vector RF1. The RF1 - Grid link extends the grid-network procedure for hydrologic modeling developed by Maidment (1992).

4.5 Redistribution of the flow record over ungauged rivers

The flow rate is essential to estimate the agrichemical concentration and load in all rivers of the area under investigation. In this study, the historical flow record is utilized instead of synthetic values. The high or low flow conditions are modeled by selecting a year from the past that had the high or low flow rate recorded. The flow measurements are available only in locations in which the USGS gauging stations are located. Therefore, a procedure that calculates the flow rate in ungauged stations has been developed. This section describes such a procedure.

4.5.1 GIS database of monthly flow rate and the precipitation depth

  • Compactness of the system. Spatial and temporal data are stored in one format that is specified by the GIS
  • Efficiency of the system. The spatio-temporal features can be viewed, queried, and processed directly by the procedures built into such GIS software as Arc/Info or ArcView
  • Easy to maintain. Data can be organized by spatial units (e.g., political or hydrologic units) and/or by temporal unit (e.g., year, month or decade)
  • Simple to understand. Data are stored in attribute tables, records represent the spatial domain whereas the columns (items) represent time domain.

The GIS database of monthly flow rate and the database of the average monthly precipitation depth is developed in two steps: first a map of monitoring stations is created and then, for each map, the attribute table with the measurements is build.

The maps are created utilizing the latitude and longitude published by Hydrosphere (1993 a, b, 1994). Although there are 38 USGS stations in the Iowa-Cedar River basin only the stations with the complete flow record (28) has been used for analysis. Figure 4.13 shows the map of the USGS gauging stations selected for modeling the spatial distribution of recorded flow rate.

Figure 4.13 The Iowa-Cedar River basin: subwatersheds and selected USGS gauging stations (numbers represent station ID).

The map of the weather stations contains 86 stations that are located within the Iowa-Cedar River basin and within the 50-km buffer zone outside the basin. Figure 4.14 shows the map of the weather stations. All stations are utilized to redistribute the flow rate record.

Figure 4.14 Weather stations applied to analysis of the hydrologic conditions in the Iowa-Cedar River watershed.

A fragment of the point attribute table that contains average monthly flow rates is presented in Table 4.1.

Table 4.1 An example of PAT--point attribute table of the USGS gauging station coverage. Full table contains monthly flow record for 38 stations, for period from 1940 to 1992 in m³/s

4.5.2 Average precipitation depth in modeling units

The process of spatial redistribution of the measured monthly average flow requires the average precipitation depth for each modeling unit. GIS water quality models such as SWAT-GRASS utilizes the rainfall depth observed in the closest weather station to subbasin (Ramanarayanan et al., 1996, Krysanova et al., 1996). The commonly used in hydrology methods of spatial estimation of rainfall from rain gauges are the Thiessen polygon method (Chow et al., 1988) and the inverse distance-squared method (Smith, 1993). There are two other functions available in Arc/Info GIS: kriging and trend (fitting a polynomial regression surface). In this research the Arc/Info IDW function (inverse distance-squared method) has been applied. It creates a grid of spatially distributed values extracted from the attribute table of point coverage. The calculations are very fast, thus it was feasible to create maps of spatial distribution of monthly average precipitation depth for the Iowa-Cedar River basin for the period from 1950 to 1992. Figure 4.15 presents the map of precipitation depth estimated for June 1990.

Figure 4.15 Spatial distribution of precipitation in the Iowa-Cedar River watershed in June 1990.

4.5.3 Mathematical description

There are about 38 USGS gauging stations in the Iowa-Cedar River basin (some of them are not in service). One or more gauging stations determine a drainage area (or partial drainage area), referred to as a gauging station zone or zone. The gauging station constitutes a point through which a known amount of water flows from one zone to an other zone. The runoff from each modeling unit depends on two factors: 1) the water balance calculated for the unit's respective zone and 2) the distribution of the precipitation depth over the gauging zone. Figure 4.16 illustrates an example of a gauging zone and modeling units.

Figure 4.16 Cedar River watershed above Charles City, Iowa: An example of gauging station zones and modeling units.

  1. Estimation of an average runoff coefficient relating monthly precipitation to discharge
  2. Approximation of the flow in streams
  3. Evaluation of error
  4. Correction of estimated flow rate.

Estimation of average runoff coefficient. Using recorded outflow from the first order zones, i.e. watersheds determined by the most upstream gauging stations, an average runoff coefficient C f is calculated according to the following equation:

C f (m) = average runoff coefficient for the m -th month [dimensionless]
Q j (m) = average flow recorded in j -th gauging station during m -th month [m³/s]
P j (m) = average monthly precipitation depth [mm/d]
A j (m) = j -th watershed area [km²]
m = month
j = index of the first order watersheds
a = units conversion factor.

First approximation of the flow in streams. The runoff from the modeling units is summed along the flow path, moving downstream from a first order stream toward the outlet of the basin. The mass balance for the i -th unit is described by the following formula:

Q' i (m) = estimated cumulative flow at the outlet point of the i -th unit [m³/s]
Q' k (m) = estimated cumulative flow at the outlet point of the k -th unit [m³/s]
k = index of units in the immediate upstream vicinity of the i -th unit
C f (m) , P i (m) , A i (m), = average runoff coefficient, average monthly precipitation depth in unit i , and i -th unit area, respectively

To make calculations more efficient, when the flow path crosses the border of the gauging station zone, i.e. a gauging station, the calculated cumulative flow is substituted by the measured value Q j (m) . This substitution of values ensures that the observed inflow is used to calculate accumulated flow in each zone, and the resulting error at the zone outlet is due only to inaccuracy of water balance estimated within the zone (errors do not propagate from zone to zone). Thus, the water balance in unit downstream to the gauging station is:

Q' i (m), C f (m) , P i (m) , A i (m), = same as for (Eq. 4.18)
Q j (m) = observed inflow into modeling unit i , from zone j . This value equals to
the outflow from zone j [m³/s]

Error evaluation . In this step, the difference between estimated and observed flow, , is calculated for each zone j :

Q j (m) = observed outflow from zone j
Q' i (m) = estimated cumulative flow at the outlet point of the i -th unit that
is also the outlet point of the j -th gauging station zone.

Correction of the cumulative flow. The correction of the estimated cumulative flow is the crucial process of the flow distribution method. Two weighting coefficients k 1 and k 2 are applied. The coefficient k 2 redistributes error according to the cumulative runoff calculated for each zone separately - the flow in rivers is created only by the estimated runoff C f P(m)A(m). There are no inflows from the upstream zones. Thus the coefficient k 2 redistributes error according to the drainage area and the precipitation depth within each zone separately. It is calculated for each modeling unit as a proportion of the runoff Q x i (m) from the drainage area upstream of the given, i -th unit outlet (within the zone j ) to the total runoff Q x j (m) from the zone j (at the beginning of the flow path it is very small, at the zone outlet it equals one):

To put more burden of the error correction on the major rivers rather, then small streams the coefficient k 1 has been introduced. There is higher probability that the estimated error is due to the losses/gains of the river between two gauging stations than to the losses/gains in small streams located far from measurements points. In another words, further from the gauged reach a stream is located, the higher the uncertainty of the flow is and therefore, the more appropriate is to apply the average basin runoff coefficient to estimate stream flow rate.

The coefficient k 1 is calculated using total cumulative flow, i.e., the approximated flow in the streams of the whole basin as a proportion of the estimated cumulative flow Q' i (m) at the outlet point of the i -th modeling unit to the cumulative flow Q' j (m) at the outlet of zone j in which the i -th unit is located:

The adjustment of estimated flow (Eq. 4.18 and Eq. 4.19) can be represented by the following formula:

Q'' i (m) = corrected cumulative flow at the outlet point of the i -th modeling unit
Q' i (m) = estimated cumulative flow at the outlet point of the i -th modeling unit
= error, difference between observed and estimated cumulative flow at the outlet of zone j
k 1 ,k 2 = weighting coefficients.

4.6 Exponential decay model

This section discusses a model which estimates loads in rivers assuming the chemical losses in rivers are governed by the first order reaction, i.e., the agrichemical mass exponentially decays as it travels from one modeling unit to the downstream unit. The preliminary model has been developed to calculate the concentrations in rivers utilizing the results of CEEPS (Comprehensive Environmental Economic Policy Evaluations System) metamodel developed by the Iowa State University's Center for Agricultural and Rural Development (Bouzaher and Monale, 1993).

The CEEPS applies two models of chemical runoff from the field: PRZM (Mullins et al., 1993) for pesticides and EPIC for nutrients. The results, among others, are proportions of agrichemical mass applied that leaves the field. Unfortunately, the author of this dissertation was unable to obtain the spatial distribution of these export factors over the Iowa-Cedar river neither to estimate the rate of agrichemical losses in rivers necessary to develop and test the complete GIS model. Since it is not in the scope of the research to develop such a model, only its outline for the completeness of discussion, is included below.

4.6.1 Exponential decay model overview

The spatial model frame is based on the watershed divided into modeling units described in previous sections. The mass that enters the stream in unit watershed is calculated using export factors:

M i = agrichemical mass that enters surface water in the i -th modeling unit
U i = total agrichemical application in the i -th unit and
E i = export factor that depends, e.g., on agricultural management practices.

As the agrichemical mass (load) travels along the flow path it decays. The decay process can be related to the travel time or to the travel distance. It is assumed here that the time decay coefficient or distance decay coefficient represent such processes as decay, sorption, desorption, volatilization, exchange between sediment and water column, settling, scour, and dilution. The equation (4.25) describes this process utilizing the travel time and the time decay coefficient:

M out,i = agrichemical mass that leaves i-th modeling unit
M out,k = input mass from the k-th upstream unit
M i = agrichemical runoff from i-th unit
K t,i = time loss coefficient
t i = time that takes the constituent to travel through i-th unit.

To describe the agrichemical losses by the exponential function of river reach length the K t,i t i component of the equation (4.25) must be substituted by the distance loss coefficient K K,i and the travel distance L i .

Since the database of the average monthly precipitation depth, the monthly water runoff from the field as well as the average monthly flow rate in river for all modeling units it is feasible to apply concentrations (traditional approach: ) instead of loads. The GIS model prototype is coded in ArcView script language. The space in the polygon attribute table is reserved for all required parameters.

4.6.2 Travel time approximation

Since no information about the river cross sections is available as yet, a special procedure has been developed to determine the distribution of time of travel and flow velocity. The visual analysis of daily flow record along the flow path has revealed that there is a time shift between flow time series. Figure 4.17 presents an example of the flow rate time series along the flow path from the gauging station "Winnebago River at Mason City, Iowa" to the gauging station "Iowa River at Wapello, Iowa".

Figure 4.17 Example of time series recorded in gauging stations located along the Cedar River. The data are for the water year 1990. The logarithmic scale is used to show all time series in one picture.

A cross-correlation coefficients can be applied to test the linear relationship between flow time series recorded in the gauging stations located on the same flow path. Six gauging stations along the Cedar River shown in Figure 4.18 have been selected for the preliminary analysis.

Figure 4.18 Flow path and location of gauging stations which have been used to illustrate the potential method for time of travel estimation.

The cross-correlation coefficients between the time series recorded in gauging station "Cedar River near Conesville, IA" and the time series recorded in remaining gauging stations have been calculated. Figure 4.19 shows the estimated coefficients. For example, the cross-correlation between the flow recorded in gauging station "Cedar River near Conesville" and the flow recorded in gauging station "Cedar River at Cedar Rapids" has the maximum value for a two-day lag. This suggests that it takes about two days for the water to flow from Cedar Rapids to Conesville.

Figure 4.19 Cross-correlation coefficients: the strength of the linear relationship between the flow rate in the Cedar River recorded near Conesville, IA, and the flow rate recorded at indicated locations. All gauging stations are on one flow path (shown in Figure 4.18).

The projection of the travel time over the ungauged regions is much more difficult that the projection of the flow velocity, except the case when the simple relation between the travel time and the travel distance is applied. Figures 4.20a and 4.20b show such relationship. Figure 4.20c presents the flow velocity as a function of the square root of the stream slope.

Figure 4.20 Analysis of travel time and flow velocity:
a) cumulative travel time versus cumulative flow distance
b) travel time vs. flow distance
c) velocity vs. square root of stream slope.

The travel time shown in Figure 4.20 represents the lag time for which the cross-correlation coefficient has the maximum value. The slope for each stream reach has been calculated from the digital elevation model (DEM). The flow length has been calculated by the Arc/Info function flowlength .

The same exercise has been repeated for the Iowa River. Two river sections have been excluded from the analysis: one between Rowan and Marshalltown (maximum correlation between no-lagged series) and other between Marengo and Iowa City (the Lake Coralville caused that the maximum correlation coefficient has been estimated for nine-day lagged series). The following relationship between "travel time", the stream slope, and the stream length has been estimated for both, the Iowa River and the Cedar River:

t = travel time [d]
So = stream slope
L = stream length [km]

The Equation (4.26) gives a good approximation of the "travel time" in the Iowa River and the Cedar River. Figure 4.21 shows that the error of prediction for all stream reaches is below 0.5 day (the precision of determining a lag time of the maximum correlation is 0.5 day.)

Figure 4.21 Estimation of the "travel time" for major reaches of the Iowa River and the Cedar River, IA.

Figure 4.22 compares the flow velocity extracted from the RF1 database with the flow velocity estimated by the correlation analysis. There is no visible relation between these two data sets. The flow velocity estimated from the correlation of flow record is about three times higher then the one published with RF1. According to the RF1 velocities, it takes constituent almost a month to travel from Austin, MN to Wapello, IA. Further analysis is needed to clarify the reasons for this discrepancy.

Figure 4.22 Comparison of the velocity from the RF1 database with the velocity estimated from a correlation analysis of the flow record in the Iowa River and tributaries
Go to: Chapter 1 Chapter 2 Chapter 3 Chapter 5 Chapter 6 Chapter 7 Home Page

VoronoiMesh as a TogglerBar

Is it possible to use a VoronoiMesh to define a TogglerBar or SetterBar type Control ?

For example, I can customise the looks of a TogglerBar

But this doesn't change the rectangle-shaped buttons and I don't think this is the way to go about if I want to define a VoronoiMesh with clickable and "toggleable" cells.

where each cell is selected/unselected whenever I click it, adding/removing a correspondent number to a list, for example, in the case of a TogglerBar . I would like this to also work as a SetterBar .

Edit 1: Thank you all for your answers. As a follow up, I'm now interested in developing a TogglerBar -type object that allows users to hold and drag the mouse to select/deselect several cells. If you have time, please take a look at it, I'm a bit a clueless on how to do this, so any hint or idea is welcome.

Edit 2: Following Lukas Lang's answer below, I also tried to vary the grid size in Manipulate

However, this doesn't seem to behave as expected. Instead, I get

Any idea why, and how to fix this? I tried Dynamic , but didn't work.

Edit 3: As a third and (hopefully) final edit, thanks to Lukas Lang's answer, I was able to solve the original question. Now I just need to define several toggler-type meshes of the same shape. One naive attempt is simply

Which naturally doesn't yield meshes with the same shape, due to the randomness in defining the points. How can I solve this? I have tried to define the mesh outside, then I lose the dynamic update of the mesh-dependent control. I would like something like the following, where I'm able to independently update similarly shaped meshes


We thank the United Kingdom Foreign and Commonwealth Office and the British Indian Ocean Territory Administration for granting us permission to undertake the research and the Marine Resources Assessment Group (MRAG) for the collection and provision of the risk data. Funding for this project was provided by the Bertarelli Foundation and contributed to the Bertarelli Programme in Marine Science. Dr Barb Block was supported by funds from Stanford University. We also thank John Pearce and James Moir Clark at MRAG for helpful discussion about BIOT enforcement data and risk.

Watch the video: How to create Grid in ArcGIS with required Dimension (October 2021).