Thursday, August 13, 2015

Terrain Processing

by James Gilbert & Christine Pribulick 
Keywords: DEM, GRASS

Manual Location
Setting Up a Real Domain – Chapter 3.1.2

If you would like to simulate a real location, you need to develop model inputs that reflect hydrologically-relevant properties of your area of interest. One such property is slope – intuitively we know that it controls how water flows in real life and in ParFlow it is a key variable in the simulation of overland and subsurface flow and defines the domain geometry when using the terrain-following grid methodology.
There are numerous methods for deriving slopes for a real landscape surface. We’ll discuss one method here, but all methods have some basic similarities: they rely on an elevation dataset (often digital elevation models, or DEM’s), and some methods or algorithms that try to find hydrologically-consistent flow paths across the surface represented by the DEM. An important factor to consider is that, at least with kinematic wave routing in ParFlow, ponding and backwater effects are not simulated. The slope inputs therefore need to represent flow paths that are continuous and monotonic downhill – slopes that locally reverse direction within a larger downslope region will tend to cause an unrealistic stagnation of surface water. (Incidentally, cells with anomalously high surface pressures are helpful in diagnosing potential problems with your slopes).
Some of the ParFlow tools listed in the manual (Section 4.2) were originally developed to help do this. In recent years, however, other methods have proven quick, useful, and robust to a wide range of ParFlow applications. In particular, the use of the A* algorithm in the GRASS GIS software has produced good results. Here we provide some of the basic steps that will get you from a DEM to ParFlow-ready slopes:
1. Geographic extent of your domain  - the following steps all involve GIS operations, so it’s helpful to have at least a shapefile (or geographic data file of your choice) of your domain extent.
2. Elevation file – you’ll need an elevation dataset for your domain. Such data are available from many sources and come in a variety of formats. You’ll likely want to project it to be consistent with other spatial datasets related to your project (or at the very least the domain dataset from step 1) and then clip it to your domain extent.
3. Download and install GRASS GIS – if you have this installed already, you’re set. If not, it can be downloaded and installed on Mac operating systems (
4. Add files and set up GRASS region for processing your DEM. If you’ve worked in ArcGIS but not GRASS, this may be a little different. GRASS expects all datasets to be in the same projection and will only perform operations on the region extent you define. If you set your region extent to the same as your domain extent (from Step 1) and load in your DEM in the same projection, you’ll be all ready to go. Here are some example steps for getting there:
  • Open new directory
  • Next select 'Read projection and data terms from georeferenced data file' (GIS DEM ASCII file previously made from GIS Domain Setup)
  • Setup Grid: Toolbar-> Settings->Region-> Set Region-> Resolution->Grid Res 2D =#### ( where #### is the resolution of your domain in the units of your projection; if #### was 1000, this would correspond to a region with resolution of 1000 m or 1km)
  • Add raster: Select ‘Add Raster Map Layer’ icon from the GRASS GIS Layer Manager window -> Select DEM raster created in GIS.
5. Run watershed analysis tool to get datasets you’ll use for calculating slopes.
  • Select Raster -> Hydrologic Modeling -> Watershed Analysis from the toolbar In the ‘Required’ tab, select your DEM raster dataset.
  • You can make changes to values in the ‘Input Options’ tab as desired. Depending on the size, resolution, and variability of the topography you’re working with, some values may produce more desireable results than others. It may be useful to experiment with values to find acceptable results.
  • ‘Output Options’: You’ll want the following outputs: Accumulation, flow direction, watershed basin, stream segments
  • ‘Optional’ tab: It is advisable, but not necessary, to select the following: "Enable MFD flow", "allow only horizontal and vertical flow of water", "Use Positive flow accumulation even for likely underestimates", and "Allow output files to overwrite existing files".
  • Run analysis tool – results should display in the map window.
  • Repeat/redo as needed to get basins, stream segments, and accumulation patterns that make sense for your domain.
  • NOTE: You can perform these steps with a DEM that is much higher resolution than your model grid will be. Simply ensure that your GRASS region is set to the correct resolution (e.g. 1000 m) and the result of the Watershed Analysis tool will be aggregated to that resolution as well.
6. Once you have datasets you feel good about, export them to ascii (text) raster files using File -> Export Raster Map -> ASCII grid export

7. Up to this point we’ve just created datasets that are hydrologically consistent with the DEM data we started with, but not specifically slopes. The next step assumes we want to treat convergent drainage zones (stream channels, or the data represented in the ‘stream segments’ result from the Watershed Analysis) or different subwatersheds differently than other parts of the landscape. For example, we may want to ensure that stream channel slopes are constrained within a certain range that is different from the rest of the model domain. Regardless of the options you’d like to apply, to get from the GRASS output datasets and ParFlow slope input files you’ll need to do the following in your favorite programming language or GIS software:
  • Calculate the x-direction and y-direction slope for each cell. This can be as simple as the elevation of one cell minus the elevation of an adjacent cell divided by the distance between cell centers. This will give you an array of x and y direction slopes where the sign of the slope indicates a direction
  • Adjust these slopes by the ‘flow direction’ dataset from GRASS. This will ensure a single continuous flow path from anywhere in the domain. If flow direction = 2, then the x slope should be zero and the y slope should be negative (sign change will depend on how you calculated slopes to begin with). If flow direction = 6, x slope is also zero but y  slope should be positive. Similarly for flow direction 4 and 8, y slopes are 0 and x slopes are positive and negative, respectively.
  • You can now use the ‘watershed basin’ and ‘stream segments’ GRASS output to apply different constraints to slopes selectively by watershed or by channel reach.
  • Once you’re satisfied, the final task is to convert separate x and y slope datasets to files ParFlow can read. You can follow the pseudocode in the ParFlow manual to write out your slope datasets to ParFlow binary (Section 6.2) or ascii formats (Section 6.7).
8. As a final step, you’ll want to check your slopes in ParFlow to make sure they yield reasonable routing patterns across the domain. One way to do this is through a “parking lot” test. In this test, you set up your domain with very low hydraulic conductivity values for the subsurface (to effectively make it impermeable) and apply one or many pulses of water. The simulated routing response should show streams and rivers forming and moving water (pressure) downstream. If rivers don’t form in the right place (or ever) or get stuck at a bad slope point, you may need to redo the steps above to correct whatever slope errors are causing the unwanted behavior.

No comments: