Using Indicator Files as input to ParFlow

A common question is how to get spatially-varying input, such as soil types, into ParFlow. Though there are a number of ways to do this, a most common approach is to use an indicator pfb file with integer values for each land coverage. This is what we've done for the coupled model where we have different land surface types (such as from IGBP) distributed over the ground surface. An example of this type of coverage is in this paper.

What we do is write the indicator value to the corresponding spatial grid location in a pfb file, then load this file in. An indicator field is a .pfb file with integer numbers for every cell in the computational domain. An example in two-dimensions might look as follows.


These are mapped to geometry names within the parflow .tcl input file and the Tcl/TK input script would look something like:

pfset GeomInput.indinput.InputType IndicatorField
pfset GeomInput.indinput.GeomNames "field channel subsurf"
pfset Geom.indinput.FileName "my_indicator_input.pfb"

pfset GeomInput.field.Value 1
pfset GeomInput.channel.Value 2
pfset GeomInput.subsurf.Value 3
You then can use the names to assign properties, say permeability

pfset Geom.Perm.Names "field channel subsurf"

pfset Geom.field.Perm.Type Constant
pfset Geom.field.Perm.Value 0.2

pfset Geom.channel.Perm.Type Constant
pfset Geom.channel.Perm.Value 0.001

pfset Geom.subsurf.Perm.Type Constant
pfset Geom.subsurf.Perm.Value 1.0
This will assign these parameter values anywhere these indicators are located in the file. You need to distribute the file before reading it in. To write the file you can modify the FORTRAN code pf_write.f90 located in the /helper_codes directory or use another program to write an ascii file and then use pftools to convert it to a .pfb file.

Subsurface Setup Options

by Nick Engdahl & Lindsay Bearup

Keywords: Solid Files, Indicator Files, heterogeneity. 

Manual Location

Manual version: July, 2014 v.693

Defining the Problem – Chapter 3.1

Little Washita Example – Chapter 3.6.2

TFG formulation – Chapter 5.3

ParFlow Solid Files (.pfsol) – Chapter 6.5

So, you’ve decided to use ParFlow and got tired of only being able to make heterogeneities that are rectangular boxes…and now you have no idea what’s going on. Domain setup is tricky, with a lot of options. Usually, we rely on two approaches and we hope this post will help to at least conceptually clarify how those two options differ.

To start, we’re going to assume that you have a processed DEM and that your domain is already defined (slopes, limits, etc…) already as either a terrain following grid (TFG) or an orthogonal (or rectilinear) grid (OG). We’ll also assume that you’ve gotten a homogeneous version of the domain to run, which is a highly recommended first-step in the modeling process. With all that done, the task we’re going to look at is how to add some complexity within those domains with spatial heterogeneity. There are two ways to do this: Indicator Files and Solid Files, and either can be used with both grid types (TFG or OG)

Common Features

The concept here is to define different zones within your model and that’s the common framework that indicators and solid files share. The model will know that certain places in the domain belong to Category 1, for example, and you can use that reference to define any of the properties ParFlow uses. You can think of this as defining hydrofacies or hydrogeologic units. Both kinds of files will need to be created externally then brought in to ParFlow and that’s where the similarities stop.

Indicator Files

An indicator file is a .pfb file that assigns an integer to each cell. In the main ParFlow .tcl script, each integer needs to be linked to a name that is then used throughout the script to assign parameters (permeability, porosity, Van Genuchten parameters, etc). In the Little Washita example, this looks something like this:

pfset GeomInput.indi_input.InputType   IndicatorField

pfset GeomInput.indi_input.GeomNames   "s1 s2 s3 s4”

pfset Geom.indi_input.FileName         "IndicatorFile_Gleeson.50z.pfb"

pfset GeomInput.s1.Value                1

pfset GeomInput.s2.Value                2

pfset GeomInput.s3.Value                3

pfset GeomInput.s4.Value                4

…

pfset Geom.Perm.Names         "domain s1 s2 s3 s4 "

pfset Geom.domain.Perm.Type   Constant

pfset Geom.domain.Perm.Value  0.2

pfset Geom.s1.Perm.Type       Constant

pfset Geom.s1.Perm.Value      0.269022595

pfset Geom.s2.Perm.Type       Constant

pfset Geom.s2.Perm.Value      0.043630356

pfset Geom.s3.Perm.Type       Constant

pfset Geom.s3.Perm.Value      0.015841225

pfset Geom.s4.Perm.Type       Constant

pfset Geom.s4.Perm.Value      0.007582087

Anytime you have specifications like this, there must be complete entries for every entry in the “Names” list or your model won’t run. If you wanted to fill-in the space represented by any of the zones, like “s3” with a distribution of permeability instead of a single value, you can swap “Constant” for turning bands or another option (see the manual for those).

Like other input files, the indicator file needs to be distributed and undistributed to run on multiple processors:

pfdist IndicatorFile_Gleeson.50z.pfb

set runname "LW"

puts $runname

pfrun    $runname

pfundist IndicatorFile_Gleeson.50z.pfb

For more information on indicator files, check out an earlier blog post, here.

Indicator files can be used for OG and TFG domains but it is important to remember that all positions in a TFG are relative to the model top. Numerically, the grid is still just a rectangle but the internal geometric transform using the slopes makes the model think it is tilted. The variable thickness is the same way; it’s a multiplier within the solver, which is unknown at the time the program does its domain setup tasks.

Many people would like a step-by-step guide describing how to make 3-D shapes in an indicator file but there are far too many ways for us to cover. The simplest way is to use some the common indicator geostatistical methods (like sequential indicator simulation, or transition probability based indicators) that produce a 3-D array of indicators. If you convert those outputs in to a PFB using pftools, you’ll be ready to go.

Solid Files

Solid files are a very different approach and are basically exactly what you would expect: solid objects. Instead of defining zones at every cell in the domain, we define the volumetric extents of solids and let ParFlow decide what portion of the domain is inside each. This has the useful property that the solids are independent of any particular grid as long as the bounds of the domain completely contain the solids. For example, if your initial grid uses a dx=5.0 with nx=10, the same solid file can be used for dx=1.0 and nx=50, where dx and nx are the cell spacing and number of cells, respectively.

The most common use of solid files is to define the extent of the domain (see the crater2d.tcl example in the /test directory and this blog post), but they can also be used to identify zones within the domain, just like indicators. In our obviously biased opinion, a great example of this is The East Inlet model described in Engdahl & Maxwell (2015). The different zones corresponding to different geologic environments are clearly visible in Figure 2 of that article and each has been filled in with a Gaussian random permeability field.

Solid files are powerful tools for representing complex shapes but they are not easy to create. This isn’t a ParFlow specific issue; no matter what your application, it can be tricky defining 3-D objects that are complete, enclosed volumes. On the plus side, the pfsol file is plain text, so you can read it. Its contents are a collection of vertices, connections and object id’s that make up a collection of small, planar objects and the outward normal of these triangular planes define the geometry. An earlier blog post explains how a utility in the pftools directory can help with the task. The simplest approach is to use objects with vertical edges along every side but the top and bottom surfaces can vary significantly. Again, the East Inlet model is a good example of the complexity the surfaces can take.

In prepping this post, we realized that none of the ParFlow test cases to date (up to release 710) contain examples of assigning properties using solid files for zones; only defining domain limits. We’ll address this in a future post and add an example that walks you through the creation of the East Inlet domain using solid files for both OG and TFG equivalents. For now, we’ll just look at the necessary keys in crater2d.tcl:

pfset GeomInput.Names        "solidinput $Zones background"

pfset GeomInput.solidinput.InputType  SolidFile

pfset GeomInput.solidinput.GeomNames  domain

pfset GeomInput.solidinput.FileName   crater2D.pfsol

This is used to define the active domain, from which the mask is computed. If you plot the mask, you’ll see only the portions of the domain that reside within the pfsol object. If another solid file was read in and given different GeomName, it could be referred to just like “domain” for any property in the domain. The same procedure shown here is merely repeated to create other zones as geometric objects.

An important point is that when you are defining zones within the domain, and not just the domain limits, the order that the solid files are read in specifies their hierarchical relationships. In other words, if two overlap, the last one that is read in wins. Another point is that strange things will happen if the solid file that defines your domain (if you are using one) doesn’t fit within the computational grid.

Again, we realize this post might not answer all your questions just yet, but stay tuned for a post with new examples soon!

References & Links

Engdahl, N.B. and Maxwell, R.M. (2015) Quantifying changes in age distributions and the hydrologic balance of a high-mountain watershed from climate induced variations in recharge, Journal of Hydrology, 522, p 152-162, doi: http://dx.doi.org/10.1016/j.jhydrol.2014.12.032

Workflow for building a model of a real domain

by Laura Condon

Keywords: Domain setup, watershed model, getting started

Manual Location

Manual version: July, 2014 v.693

Defining the Problem – Chapter 3.1

TFG example (Little Washita) – Chapter 3.6.2

TFG formulation – Chapter 5.3

ParFlow Solid Files (.pfsol) – Chapter 6.5

Basic Workflow for Setting up a Real Domain
This section is copied from the ParFlow manual section 3.1. This workflow is for setting up a real domain for ParFlow CLM from a DEM using the terrain following grid with an indicator file to define the subsurface. Refer to this blog post for a discussion of the differences between terrain following grid and orthogonal grids. If you are setting up an orthogonal grid this general workflow is still applicable but steps 4-6 will be slightly modified. 

The general approach is as follows:

1. Gather input datasets to define the domain. First decide the resolution that you would like to simulate at. Then gather the following datasets at the appropriate resolution for yourdomain:(a) Elevation (DEM)(b) Soil data for the near surface layers (c) Geologic maps for the deeper subsurface (d) Land Cover

2. Create consistent gridded layers that are all clipped to your domain and have the same number of grid cells

3. Convert gridded files to .pfb (see manual section 6.2). One way to accomplish this is by reformatting the gridded outputs to the correct ParFlow .sa order (§ 6.7) and to convert the .sa file to .pfb using the conversion tools (see manual section 4.3 Example 1)

4. Calculate slopes in the x and y directions from the elevation dataset. This can be done with the built in tools as shown in the manual section 4.3, example 5. In most cases some additional processing of the DEM will be required to ensure that the drainage patterns are correct. To check this you can run a “parking lot test” by setting the permeability of surface to almost zero and adding a flux to the top surface. If the results from this test don’t look right (i.e. your runoff patterns don’t match what you expect) you will need to go back and modify your DEM. The built in ParFlow tools pitfill and flatfill can be used to address some issues. These tools are also shown in the manual section 4.3 Example 5. For a more thorough workflow refer to this blog post on DEM processing.

5. Create an indicator file for the subsurface. The indicator file is a 3D .pfb file with the same dimensions as your domain that has an integer for every cell designating which unit it belongs to. The units you define will correspond to the soil types and geologic units from your input datasets. This blog post contains more information on indicator files and the Little Washita annotated input script in the manual (section 3.6.2) uses an indicator file to define the subsurface. 

6. Determine the hydrologic properties for each of the subsurface units defined in the indicator file. You will need: Permeability, specific storage, porosity and vanGenuchten parameters.

7. At this point you are ready to run a ParFlow model without CLM and if you don’t need to include the land surface model in your simulations you can ignore the following steps. Either way, at this point it is advisable to run a “spinup” simulation to initialize the water table. There are several ways to approach this. One way is to start with the water table at a constant depth and run for a long time with a constant recharge forcing until the water table reaches a steady state. There are some additional key for spinup runs that are provided in the manual section 6.1.34. This blog post has some addition tips on spinning up a model. 

8. Convert land cover classifications to the IGBP1 land cover classes that are used in CLM.

9. Create a CLM vegm file that designates the land cover fractions for every cell (Refer to the clm input directory in the Washita Example for an sample of what a vegm file should look like).

10. Create a CLM driver file to set the parameters for the CLM model (Refer to the clm input directory in the Washita Example for a sample of a CLM driver file).

11. Assemble meteorological forcing data for your domain. CLM requires the following variables (also described on p. 136 of the manual):

• DLWR: Long wave radiation [W/m2]
• DSWR: Visible or short-wave radiation [W/m2].
• APCP: Precipitation [mm/s]
• Temp: Air Temperature [K]
• UGRD: East-west wind speed [m/s]
• SPFH: Specific humidity [kg/kg]
• Press: Atmospheric pressure [pa]
• VGRD: South-to-North wind speed [m/s]

If you choose to do spatially heterogeneous forcings you will need to generate separate files for each variable. The files should be formatted in the standard ParFlow format with the third (i.e. z dimension) as time. If you are doing hourly simulations it is standard practice to put 24 hours in one file, but you can change this as needed. For an example of hetergeneous forcing files refer to the NLDAS directory in the Washita Example. Alternatively, if you would like to force the model with spatially homogenous forcings, then a single file can be provided where each variable is a column and rows designate time steps.

12. Run your simulation! Section 3.6.2 of the manual has an annotated input script for the Little Washita example which is a good place to start when setting up your run script. 

If your model doesn't run the first time don't worry there are some troubleshooting blog posts here and here and you can also check the archives of the ParFlow users mailing list here. 

Terrain Processing

by James Gilbert & Christine Pribulick 

Keywords: DEM, GRASS

Manual Location

Manual version:July, 2014 v.693

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 (http://grass.osgeo.org/download/software/mac-osx/).

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.

Domain Options: Terrain Following vs. Orthogonal Grids

by Lindsay Bearup

Keywords: Terrain Following Grid, Solid Files, Domain Geometry

Manual Location

Manual version:July, 2014 v.693

Defining the Problem – Chapter 3.1

TFG example (Little Washita) – Chapter 3.6.2

TFG formulation – Chapter 5.3

ParFlow Solid Files (.pfsol) – Chapter 6.5

Two primary methods are available for defining the gridded domain of a watershed in ParFlow: orthogonal methods (implemented using solid files) and the terrain following method (TFG). The strengths and limitations of each approach are discussed below.

Depiction of orthogonal (upper) and the terrain following (lower) grid formulations for a topographic land-surface and (left) schematics of the associated finite difference stencils. Note in this figure that i, j, and k are the x, y, and z cell indices, respectively, and that the limits of the domain would be at i = nx, j = ny and k = nz. Note also that nz is different (greater) in the orthogonal formulations (upper) than in the terrain formulation (lower) due to the difference in inactive cells between formulations. Figure from: Maxwell 2013 AWR.

Terrain Following Grid (TFG)

The terrain following grid formulation transforms the ParFlow grid to conform to topography. This formulation can be very useful for coupled surface-subsurface flow problems where groundwater flow follows the topography. This method also requires fewer inputs because it uses the provided surface slope files to determine the geometry. The following key is used to run with TFG and can only be used with Solver Richards (currently not available with Solver Impes): 

pfset Solver.TerrainFollowingGrid           True

TFG is typically applied using a "Box" input type to define the discretization. See the Little Washita example in Chapter 3.6.2 of the ParFlow Manual and the LW_var_dz.tcl script in the ParFlow test directory for an example of the TFG formulation applied to a watershed with a variable vertical discretization. 

pfset GeomInput.domaininput.GeomName  domain

pfset GeomInput.domaininput.InputType  Box 

Solid Files

Solid files rely on traditional orthogonal formulations of the geometry and discretization of the model domain. This approach provides more freedom to define irregular geologic layers and nonuniform watershed depths but requires the user to generate a pfsol file to define active/inactive cells. A solid file is a triangulated information network (TIN) file that defines the domain using a system of points, triangles, and patches. Patches are critical for applying boundary conditions to the domain and are listed in a specific order. For more information regarding solid files and patches, check out the blog entry on patch order here.  A fortran file named pf_solid_file_create.f90 is shipped with ParFlow (v.693 at the time of writing) and is located at $PARFLOW_DIR/pftools/prepostproc. This script can be modified to generate a solid file from a DEM. As written, this script assumes that the rest of the domain is rectangularly shaped and is the same size as the extents of the DEM. The keys required to use a solid file are:

pfset GeomInput.solidinput.InputType  SolidFile

pfset GeomInput.solidinput.GeomNames  domain

pfset GeomInput.solidinput.FileName   name.pfsol

References 
Maxwell, Reed M. A terrain-following grid transform and preconditioned for parallel, large-scale, integrated hydrologic modeling. Advances in Water Resources 53, doi:10.1016/j.advwatres.2012.10.001, 2013.