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. 

patch order and solid files...

I've gotten a few emails regarding patch order and pfsolid files so I figured it was time to post something here. Solid files provide a lot of flexibility in boundary conditions but a few rules need to be followed. The order in which the patches are specified in the solid file needs to correspond to the order in which the patch names are specified in the tcl script. This is important because ParFlow does not know the names or attach any meaning to them; a bc could be named X1, x-upper, front or Cole and ParFlow really only knows it as "patch 0" -- the names are for our benefit.

If you use the pf_solid_file_create.f90 program I wrote in the /helper_codes directory, it has a specific patch order which is as follows (these are comments from the code starting at line 114):
! notes:
! patch(1,:) = upper surface (z = zmax)
! patch(2,:) = lower surface (z = z0)
! patch(3,:) = x = x0
! patch(4,:) = x = xmax
! patch(5,:) = y = y0
! patch(6,:) = y = ymax

This patch order corresponds to an input block in your tcl script and these need to match exactly as follows:
#---------------------------------------------------------
# The Names of the GeomInputs
#---------------------------------------------------------
pfset GeomInput.Names "solidinput"

pfset GeomInput.solidinput.InputType SolidFile
pfset GeomInput.solidinput.GeomNames domain
pfset GeomInput.solidinput.FileName mysolidfile.pfsol

pfset Geom.domain.Patches "z-upper z-lower x-lower x-upper y-lower y-upper"

The last line is the important one here as the names get assigned to a numbered list, z-upper is patch(1:) in the F90 code, z-lower is patch(2:), etc. This is the only connection between the solid file and ParFlow BC's, so it needs to be set up correctly. Once it is set up these named boundary pataches can be used to define initial and boundary states for the simulation.

Happy modeling!

-Reed

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.

Using Tcl/TK to load files and generate input parameters

A question that comes up frequently in ParFlow is how to load in external parameters from a file. This could be for transient forcing for boundary conditions, for Monte Carlo simulation (in addition to changing a random seed) or even in calibration using parameters generated from another program, such as PEST. The example below was from a question recently asked by a ParFlow user and demonstrates how this can be done for a Monte Carlo simulation, though this same approach is adaptable to many other situations.

First, one would load in an array of values from a file using tcl. This is not in the ParFlow manual or addendum, but you can easily do this with native Tcl/TK commands. We use the example from the harvey_flow.tcl. Right now that script does something like this:

set fileId [open stats4.txt r 0600]
set kgu [gets $fileId]
set varu [gets $fileId]
set kgl [gets $fileId]
set varl [gets $fileId]
close $fileId

pfset Geom.upper_aquifer.Perm.GeomMean $kgu
pfset Geom.upper_aquifer.Perm.Sigma $varu
This allows us to provide parameters externally, such as from PEST for calibration (as was done in this example). To modify things, one could do something like this, say to load in file of n Kg's:

set fileId [open stats4.txt r 0600]
for {set i 1} {$i <= $n} {incr i 1} {
set kg($i) [gets $fileId]
}
close $fileId
Farther down in the harvey_flow.tcl script, where we loop over the realizations and change the random seed, we can also loop through the file of values we've read into the arrays for Kg:

#
# Loop through runs
#
for {set k 1} {$k <= $n_runs} {incr k 1} {
#
pfset Geom.upper_aquifer.Perm.Seed [ expr 33333+2*$k ]
pfset Geom.upper_aquifer.Perm.GeomMean $kg($k)...

This approach should work for a number of circumstances as long as you can generate a text file of values outside of the ParFlow simulation script (or one could even generate some relationships within the Tcl script). This can easily be done in FORTRAN, in a spreadsheet or other mathematical program. This approach can easily be modified to load data for wells, boundary conditions or any transient condition as well.