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