Routing Streamflow

This page describes methods for setting up the separate model to take surface runoff and baseflow outputs from each VIC model simulation cell and route it to a fixed location, typically sites with USGS gauging stations.  There are currently two standard models for routing VIC model output (runoff and baseflow) used at Purdue.

It is assumed that the reader has completed defining their study domain as described in VIC model Setup (Step 1) - Establishing the Domain before completing the steps in this tutorial.

Steps here assume that the user is working with ArcGIS and has the ArcHydro toolset installed. Most of the steps involved should translate directly to geoprocessing algorithms available in QGIS, SAGA and other commonly available, non-proprietary GIS packages.

Original Lohmann Routing Model

This the the traditional way to route streamflow from the VIC model, and it has been used for almost all VIC model applications up until Gavin Yang developed the GIS based model described below. The best documentation for setting up and running the Lohmann routing model using VIC model output can be found at the University of Washington\'s VIC model web site{.external-link rel="nofollow"} under the heading "Creating a Hydrograph".

Note

The original Lohmann model is no longer supported, and has been replaced by RVIC. Limited information is available on-line about what has changed, so the assumption is that the change was mostly in renaming the original code.

PDF copies of Dag Lohmann's manuscripts related to the routing model are available through the Wiki library.

There have been some modifications to the routing model that are not part of the standard release. Links to modified versions will be made available, in particular to the version that will read a list of grid cell coordinates rather than requiring data be provided as an ASCII ArcInfo grid file. This version allows for routing of streamflow from projected grids, such as the equal area projection used for the urban impacts study.

Two versions of the routing model exist:

• rout - Original routing model, compiled for a maximum of 500 rows and columns (NROW and NCOL in rout.f) and 200 years (NYR).
• rout_hires - Modified to resample higher resolution routing information to courser resolution VIC model runs. Compiled for a maximum of 800 rows and columns (NROW and NCOL in rout.f) and 200 years (NYR). Routing control file must now contain the VIC model resolution as the last item. The direction file read routine has been modified to use ArcGIS/ArcHydro flow direction values rather than the original routing model values. Finally, station files now provide station location as latitude and longitude (these must fall within the high-resolution grid cell representing the stream channel) as this information is easier to extract from ArcGIS, especially if you check the delineation of the watersheds in ArcHydro and use the AddXY_Management tool to obtain latitude and longitude for the WatershedPoints layer.

Advantages: Plenty of people have used the model previously so help is available; Documentation is mature; Relatively easy to setup up and setup files for many major river basins have already been developed; Obtaining flow from sub-watersheds is as easy as adding the coordinates to the station location file.

Disadvantages: A single unit hydrograph is defined for all grid cells within the watershed - no matter how large the watershed; Little information exists about calibration or how best to set many parameters, therefore most applications ignore those parameters which means the model is not being used to the best of its ability; In most cases applications use the same parameters across the watershed, even when the option is there to have them vary, so spatial variation in streamflow is lost; Actual simulation times can get very long for large basins, since the model does all calculations on the fly; The original Fortran code has limitations that can make running the model annoying, including parameters that must be set in the source code and recompiled, inflexible input file formats, and unhelpful error messages.

GIS Routing Model

Developed by Gavin Yang (see Yang et al., 2011 the GIS based routing model improves on some of the limitations of the original Lohmann routing model by introducing significant GIS pre-processing using a higher resolution DEM and resultant stream network to estimate hydrographs of streamflow arrival from each grid cell to the selected outlet. The pre-processing requires are more significant than what is required for the Lohmann model, and must be conducted for each outlet point, but once completed actual routing of VIC model output to the defined outlet is a quick process.

Advantages: Once setup is complete, the model runs very quickly; Setup files preserve spatial variations in stream networks throughout the basin; Most input parameters are derived from the DEM, so calibration is not needed; Routing model can handle inputs for any time step (though I believe changing the time step requires developing new input files).

Disadvantages: Significant setup required with GIS tools; New setup files must be developed for all outlet points; New model with minimal documentation and very few watersheds have been setup (currently limited to the White River in Indiana); Current version of the model requires pre-processing steps in ArcView, which is no longer supported by ESRI or by ECN, so these steps will have to be migrated to ArcGIS before the model can be widely applied.

Setting up the Routing Model for a new river basin

Both routing models start with the delineation of the river basin. Instructions for delineating a watershed using the ArcHydro module within ArcMap is provided here:

• Lab Document with Instructions for delineating the St. Joseph River watershed in northeastern Indiana.
• Tar file with the ArcGIS files required to run the lab.

Setup of the Lohmann Routing Model Using ArcGIS

The original Lohmann routing model has been replaced with RVIC, documentation and source code can be found at http://rvic.readthedocs.io/en/latest/.

Setup for the GIS Routing Model Using ArcGIS

Here I have attempted to compile all of the documents developed by the various students who have worked on making this routing model function for VIC model simulations (Gavin Yang, Chun-mei Chiu, Sarah Rutkowski and Wei-Chih Chen) into a single set of instructions for completing the pre-processing steps in ArcGIS and Python.

Theory behind this document please refers to Liu et al. (2003):

• A diffusive transport approach for flow routing in GIS-based flood modeling, Journal of Hydrology.
• WetSPA, http://www.vub.ac.be/WetSpa/

For all the documents of the GIS-based routing model for the VIC model, please also refer to Yang et al. (2011):

• Yang, G., L. C. Bowling; K. A Cherkauer, and B. C Pijanowski (2011), The impact of urban development on hydrologic regime from catchment to basin scales, Landscape and Urban Planning, doi: 10.1016/j.landurbplan.2011.08.003.

All tools used for model setup have been compiled into a single ArcGIS toolbox, available on RCAC systems in the directory /depot/phig/apps/ArcGIS_toolboxes (mount disk as Y: drive for most direct application of tools, other drive letters will work but will require more changes to input file paths). The tool box can also be downloaded here: GIS Routing Model Tools.tbx.

1. Input Files Required

The following files must be clipped to the watershed boundary and use the same projection. Typically this analysis is done in an equal area or UTM projection. See the VIC model pre-processing documentation for help in preparing these files.

• A high-resolution raster DEM
• A raster land use map of the same resolution as the above mentioned DEM
• A polygon file that demarks the boundaries of the VIC model solution elements (grid cells or HUCs depending on the application), where each polygon has a unique value, typically the grid cell number file. Should be in the resolution used by the VIC model, but projected to match the raster files.

For the tutorial, the following DEM and land use maps are downloaded from the USGS seamless data server. \

• A 100 meter digital elevation map, and
• A land use raster (2006 NLCD).

2. ArcHydro processing of routing model inputs

You should first complete the model domain setup and hydrologic analysis described in VIC Model Pre-Processing Instructions. This will move you through the first several steps required to setup the routing model, which requires delineated watersheds for each point at which you want output from the routing model.

Note

The GIS routing model accounts for differences in roughness between overland flow and channelized flow. When flow becomes channelized is controlled by the threshold value use to create the stream network. In the ArcHydro Stream Definition tool this is set as either \'Number of cells:\' or if the data is projected as \'Area (square km):\'. There is no standard value for this, it should be set so that stream channels develop in the GIS analysis at about the same location as in the real world. Set the number too small, and your dense channel network will drain the watershed too quickly, set it too large and water will move too slowly off of the land scape.*

3. Routing model set-up steps for entire simulation domain

The steps described here are specific to the GIS based routing model, and are not currently supported outside of ArcGIS.

WARNING: The steps here are based on the raw DEM, so the travel times may not reflect the watershed boundaries generated by ArcHydro if there was substantial tweaking of inputs to get the correct delineation (e.g., burning streams, putting up walls, filling sinks). Determining whether having some grid cells that are artificially too close to the basin outlet (low travel times, since they drain out of the watershed) makes a difference in the routing model inputs will have to be left for future assessment.

1. Download and install the routing model setup toolbox.

   • Put the '.tbx' file in an easy to find directory.
   • Open the ArcMap project where you completed the watershed delineation.
   • Open ArcToolbox from ArcMap by clicking on the red toolbox icon.
   • Right click on 'ArcToolbox' at the top of the Toolbox window, and select 'Add Toolbox...'
   • Use the browse window that opens to find the '.tbx' file, then click 'Open'.
   • This will open a Toolbox called 'GIS Routing Model Tools' (or 'Toolbox' for the original version).
   • NOTE: Some of the tools may be unable to run (shown by a red 'X' next to their name). To fix this each component with the 'X' must be opened in the model builder editor (right-click on name and select 'Edit') and then 'Saved'. This will correct problems introduced by the movement of the toolbox to a new workspace.

2. Set the environment variables

   • Whenever starting a new ArcGIS session (even opening an existing project) you should double check that the environment settings are correct for your current application.

   • • Set the Current Workspace: From the ArcMap Project select Geoprocessing -> Environment... -> Workspace, then make sure that the Workspace is set to the correct directory. This is where all output from the processing steps below will be put. It works best if this is the same directory used by ArcHydro, so that it contains the Layers folder and the archydro geodatabase.

3. Create RoutingLayers directory for storing all output from the GIS routing model setup process

   • Select GIS Routing Model Tools -> 0.

4. Compute the Shreve stream order number for the watershed.

   • Select GIS Routing Model Tools -> 1. Compute Shreve Stream Order

     • Fdr (default = Fdr) - Flow direction raster. (ArcHydro output).
     • Str (default = Str) - A Boolean raster where '1' means the cell is int he stream channel, and '0' means that it is not. (ArcHydro output).
     • ShreveStreamOrder (default = ShreveStreamOrder) -

   • This is used to scale the Manning's number for the channels, so DO NOT use the Strahler stream order it will change the channel roughness results.

5. Compute the slope of land and channel surfaces

   • Select GIS Routing Model Tools -> 2. Slope of Land-River

     • Minimum slope threshold (default = 0.01%) - Slopes below this will be reset to the minimum value
     • Stream Network Raster (default = Str) - A Boolean raster where '1' means the cell is int he stream channel, and '0' means that it is not. (ArcHydro output).
     • Filled DEM (default = 'Fil') - A DEM raster file with sinks filled. (ArcHydro output).

   • Slope is calculated differently for the land surface and along stream channels, so you should use the output from this routine NOT from the ArcHydro slope calculation.

6. Compute the Hydraulic Radius

   • Select GIS Routing Model Tools -> 3. Estimate Hydraulic Radius

     • ap (default = 0.05) - Parameter for hydraulic radius calculation. Value changes by return period of flood (T) being used to estimate travel times. Values of ap are T2 = 0.05, T10 = 0.12, and T100 = 0.18. Default uses two year return period (T2) for normal flow conditions.
     • bp (default = 0.48) - Parameter for hydraulic radius calculation. Value changes by return period of flood (T) being used to estimate travel times. Values of bp are T2 = 0.48, T10 = 0.52, and T100 = 0.55. Default uses two year return period (T2) for normal flow conditions.
     • Flow Accumulation (default = Fac) - Flow accumulation file. (ArcHydro output).
     • HydraRadius (default = hydraradius) - Output raster file.

   • This function will create a hydraulic radius grid based on the flow accumulation grid. Currently locked to a 2 year return flood, will look into activating other flood options in the future. Select a flood frequency from the pop up window e.g. 2 year, 10 year or 100 year flood. Generally, a flood frequency with 2-year return period is chosen for normal floods.

7. Compute the Manning's N for land surface and channel

   • Select GIS Routing Model Tools -> 4. Compute Manning's N

     • Shreve Stream Order (default = ShreveStreamOrder) - Raster layer generated from Step 3, above.
     • Land Use Map (default = landuse) - Raster layer classified using the GIS routing land use classes, see Creating Vegetation parameter Files Using MODIS Data{.external-link rel="nofollow"} for instructions on reclassifying the raster.
     • Land Use Remap Table (default = landuse_remap.dbf) - ArcGIS reclass table associating the GIS routing land use classes to values for Manning's N.
     • ManningGrid (default = ManningGrid) - Output raster file.
     • maxStreamOrder (default = 50) - Maximum Shreve stream order for scaling purposes. Should be at least equal to maximum value from Step 3 output grid.
     • minStreamOrder (default = 1) - Minimum Shreve stream order for scaling purposes. Should be at least equal to minimum value from Step 3 output grid, which should be 1.
     • maxManningN ( default = 0.075) - Maximum Manning's N for stream channel. This represents a sluggish, reedy natural channel with deep pools (see http://www.fsl.orst.edu/geowater/FX3/help/8_Hydraulic_Reference/Mannings_n_Tables.htm){.external-link rel="nofollow"}, and will produce the slowest flows in the smallest channels (minStreamOrder).
     • minManningN (default = 0.035) - Minimum Manning's N for stream channel. This represents a clean, straight, full flowing natural channel, adn will produce highest flow velocities in the larger channels (maxStreamOrder).

   • This script computes values of Manning's N along the stream channel and across the hillslopes. Channel values are scaled based on the Shreve Stream Order number between the provided minimum and maximum Manning's N values. Hillslope values are based on the land use class, with values stored in the table landuse_remap.dbf.

   • By defining maximum and minimum stream order independent of the range of Shreve Order values in the current domain, the resulting channel roughness values will be the same no matter whether routing files are being generated for the entire watershed or just for a sub-catchment.

8. Compute flow velocity

   • Select GIS Routing Model Tools -> 5. Calculate Flow Velocity

     • Land-River Slope (default = LandRivSlope) - Raster output from Step 4.
     • Hydraulic Radius (default = HydraRadius) - Raster output from Step 5.
     • Mannings N (default = ManningGrid) - Raster output from Step 6.
     • maxVelo (default = 3.00 m/s) - Maximum velocity of water in the biggest channels.
     • minVelo (default = 0.005 m/s) - Minimum velocity of water on the roughest slopes.

   • This function creates a flow velocity grid based on the Manning's n, hydraulic radius and slope grid.

9. Estimate travel time in hours

   • Select GIS Routing Model Tools -> 6. Compute Travel Time in Hours

     • Flow Direction (default = Fdr) - Flow direction file. (ArcHydro output).
     • Velocity (deafult = velocity) - Output from Step 7.
     • Travel Time Raster (default = t0_h) - Output: raster of travel time to the outlet in hours.
     • Travel Time ASCII Raster (default = t0_h.asc) - ASCII raster output file.

   • This function creates a flow travel time grid in hours from each cell to the catchment outlet using the weighted FLOWLENGTH routine.

4. Routing model set-up steps for each watershed

The steps described here will be repeated for each watershed representing a point where you want routed streamflow. Watershed boundaries are derived from the ArcHydro Batch Watershed Delineation Tool, where watersheds (NOT subwatersheds) are delineated for Batch Points. Each Batch Point should represent a location where routed streamflow is desired. Steps described below are specific to the GIS based routing model, and are not currently supported outside of ArcGIS.

1. Extract travel time raster statistics for each watershed.

   • Select GIS Routing Model Tools -> 7. Calculate Watershed Specific Routing Statistics

     • Travel Time Raster (default = t0_h) - A raster file of travel times in hours from each grid cell to the overall basin outlet.

     • Watershed (default = Watershed) - The multi-polygon file that is output from the ArcHydro pre-processing 'Batch Watershed Delineation' function. It contains one or more polygons representing the delineated boundaries for each defined outlet point.

     • Model Simulation Cell Numbers (default = None) - A polygon file representing each VIC model simulation cell, should be in the projection used for development of the routing model parameters.

     • Cell Number Field (default = gridval) - The field in Model Simulation Cell Numbers polygon that provides VIC model cell numbers as defined in the soil file and used for simulations.

     • Fraction Threshold (default = 0.01 or 1%) - Minimum fraction of VIC simulation cell that must be included in the watershed being processed for that cell to be included in the watershed drainage map.

     • Output file name format (default = T0_%Value%) - Name used to store the extracted files, MUST include %Value% so that new files have unique names for each watershed. Default is that %Value% is set to the 'Name' column of the Watershed attribute file.

       Using the 'Name' column is a problem, since names can result in filenames exceeding the character limit for raster files. I converted the program to use HydroIDs column, but not sure that exists everywhere, might want it to just be the ObjectID column, which will always exist but is dependent on the order of watersheds in the input file.

   • This script cycles through the Batch Watershed polygon file from ArcHydro. For each watershed is masks off the delineated domain and extracts travel times for the current watershed are adjusted by subtracting the travel time from the current outlet to the basin outlet for all pixels. This sets the outlet travel time to zero. Any values that end up negative (because they drain out of the watershed) are set to No Data.

   • VIC simulation cells (zones) with less than Fraction Threshold percent of the cell contributing to watershed streamflow are excluded. (default removes VIC simulation cells with contributing areas of less than 1%).

   • The mean and standard deviation of travel times are found for each VIC model simulation cell.

   • Output for each defined watershed is a table with complete zonal statistics, where each zone is a VIC simulation cell.

2. Export GIS routing model input files

   • Select GIS Routing Model Tools -> 8. Create Routing Input Files - Gridded / Non-Gridded
   • This script extracts the files required for the GIS routing model to work with VIC model simulation output. Two sets of instructions are provided, one for traditional gridded applications, the second for non-gridded (e.g., HUC) VIC model applications. The files required by the routing model are the same, but their format must change to deal with the configuration of the data.

   • GRIDDED: Using gridded VIC model domain and gridded inputs to the routing model (repeat for each routing watershed)

     • Requires the following parameters:
       • Batch Delineated Watersheds (default: Watershed) - Polygon output file from ArcHydro batch watershed delineation, which contains the full watersheds, not clipped sub-catchments. Used in the previous step.
       • FieldID (default: HydroID) - Field within the Batch Delineated Watersheds layer that contains unique watershed identification number. This will be used for filtering and naming all output files.
       • VIC Cell Number Raster (default: NONE) - This should be a raster file or layer that is in the correct resolution and projection for the VIC model simulation. It will be used as the template for creating all routing model input files.
       • Routing Model Output Path (default: NONE) - This is the path to which the final routing model input files should be written.
       • Output Workspace (default: RoutingLayers) - Location of all previous routing model development data.
     • Within arcinfo, relate the VIC grid cell raster to the Zonal Statistics table using the gridcode.
     • Create two new new rasters using the 'MEAN' and 'STD' columns of the statistics table.
     • Create a mask file using the Spatial Analyst 'Con' function to set grid cells with MEAN values > 0 equal to '1'. This indicates that the grid cell should be included in the routing calculation.
     • Export all three grid files to ASCII raster files, for each delineated watershed.

   • NON-GRIDDED: If you are using non-gridded inputs to the routing model, either because the VIC model is setup on non-gridded simulation cells, or because you prefer not to regrid the data (repeat for each routed watershed)

     • [The contents of the Zonal Statistics table must be converted from DBF to a CSV file. Easiest way is to open the DBF file in Excel and save to CSV.]{style="color: rgb(255,102,0);"}
     • [*** Now have to modify Fortran code to work with this file!!!! *** ]

     Note

     This program also needs to export the fraction of each VIC simulation cell that contributes runoff and baseflow to stream discharge. This then needs to be reformatted, probably with an external script, to produce the fraction file for the routing model. Really, the routing model should be able to build the file lists directly, and make use of a file with contribution fraction directly.

     Actually, need to output FractValue once it has been masked by ZonalRoutingMask in CalculateWatershedRoutingStatistics.py

5. Running the Routing Model (Fortran code)

The steps described here will be repeated for each watershed draining to a point where you want routed streamflow. Watershed boundaries are derived from the ArcHydro Batch Watershed Delineation Tool, where watersheds (NOT subwatersheds) are delineated for Batch Points. Each Batch Point should represent a location where routed streamflow is desired. Steps described below are specific to the GIS based routing model, and are not currently supported outside of ArcGIS.

1. Compute the unit hydrographs using the first and second moments of travel time (mean and variance of travel time) based on the solution from Olivera and Maidment (1999).

   • Run the Fortran program iuhcell.f90

     • Usage: iuhcell

     • This creates the unit hydrograph file using a convolution integral described in Yang (2011) and Maidment (1993).

     • Inputs include ASCII raster files for the

     • 1. watershed mask,
       2. mean travel time, and
       3. standard deviation of travel time.

     • Output is an ASCII table defining the unit hydrograph from every cell within the watershed to the outlet, based on the selected time step.

   • This program also allows the user to specify whether they would like 3 hour streamflow routing (to use with the conv3h.f90 routing program) or daily streamflow routing (for the convdaily.f90 program).

   • The choice of routing time controls the time step of the VIC model output file used for routing, thus a 3 hourly output file is required for 3 hourly routing, and a 24 hourly file is required for daily routing.

2. Route the current streamflow

   • For Daily VIC model output, use the program convdaily.f90

     • Usage: convdaily

     • Program requires:

       • The path where flux files are stored,

       • The ASCII raster Mask file from the GIS processing

       • The Unit Hydrogaph file produced by iuhcell.f90, and

       • ~~An ASCII grid file with the fractional contribution of each grid cell to overall streamflow. The mask file can be used if it is assumed that all cells still in the watershed fully contribute to watershed discharge.~~

       • This last item appears to actually be a list file with the following columns: route id, vic id, flux id, fract

         Going to have to report back.

         The route ID and VIC ID are not used in convdaily.f90. flux ID should be the name of each VIC output file, but without the path that was provided directly to the program. Final column is the fractional contribution for each grid cell.

         The current version of this code does not accept VIC model output headers, and assumes that runoff is column 6 and baseflow is column 7. The version installed on /depot/phig/apps/linux as been modified to work with IN-CCIA VIC_Drianage_Out files, which have headers and an extra column (air temperature) before runoff and baseflow. This code must be updated to be more flexible. Should it stay in fortran, or be converted to C? That will depend on how easily changes can be made in Fortran. C libraries exist to handle headers and variable column output in C, but have not been developed for Fortran - can they / should they be jointly compiled?

     • Output from the program is routed streamflow in...

Existing Routing Model Datasets

• Need to pull contents of existing routing model setup files from PBWorks site 17 Mar 2016