6. ParFlow Files
In this chapter, we discuss the various file formats used in ParFlow. To help simplify the description of these formats, we use a pseudocode notation composed of fields and control constructs.
A field is a piece of data having one of the field types listed in Table 6.1 (note that field types may have one meaning in ASCII files and another meaning in binary files).
field type |
ASCII |
binary |
---|---|---|
|
||
|
||
|
||
|
||
|
Fields are denoted by enclosing the field name with a <
on the left
and a >
on the right. The field name is composed of alphanumeric
characters and underscores (_
). In the defining entry of a field,
the field name is also prepended by its field type and a :
.
The control constructs used in our pseudocode have the keyword
names FOR
, IF
, and LINE
, and the beginning and end of
each of these constructs is delimited by the keywords BEGIN
and END
.
The FOR
construct is used to describe repeated input format
patterns. For example, consider the following file format:
<integer : num_coordinates>
FOR coordinate = 0 TO <num_coordinates> - 1
BEGIN
<real : x> <real : y> <real : z>
END
The field <num_coordinates>
is an integer specifying the number of coordinates to
follow. The FOR
construct indicates that <num_coordinates>
entries follow,
and each entry is composed of the three real fields, <x>
, <y>
,
and <z>
. Here is an example of a file with this format:
3
2.0 1.0 -3.5
1.0 1.1 -3.1
2.5 3.0 -3.7
The IF
construct is actually an IF/ELSE
construct, and is used to describe input
format patterns that appear only under certain circumstances. For example, consider
the following file format:
<integer : type>
IF (<type> = 0)
BEGIN
<real : x> <real : y> <real : z>
END
ELSE IF (<type> = 1)
BEGIN
<integer : i> <integer : j> <integer : k>
END
The field <type>
is an integer specifying the “type” of input to
follow. The IF
construct indicates that if <type>
has value 0,
then the three real fields, <x>
, <y>
, and <z>
, follow.
If <type>
has value 1, then the three integer
fields, <i>
, <j>
, and <k>
, follow. Here is an example
of a file with this format:
0
2.0 1.0 -3.5
The LINE
construct indicates fields that are on the same line of
a file. Since input files in ParFlow are all in “free format”, it is
used only to describe some output file formats. For example, consider
the following file format:
LINE
BEGIN
<real : x>
<real : y>
<real : z>
END
The LINE
construct indicates that the three real
fields, <x>
, <y>
, and <z>
, are all on the same line. Here is an example of a file
with this format:
2.0 1.0 -3.5
Comment lines may also appear in our file format pseudocode. All text
following a #
character is a comment, and is not part of the file format.
6.1. Main Input File (.tcl)
The main ParFlow input file is a TCL script. This might seem overly combersome at first but the basic input file structure is not very complicated (although it is somewhat verbose). For more advanced users, the TCL scripting means you can very easily create programs to run ParFlow. A simple example is creating a loop to run several hundred different simulations using different seeds to the random field generators. This can be automated from within the ParFlow input file.
The basic idea behind ParFlow input is a simple database. The database contains entries which have a key and a value associated with that key. This is very similiar in nature to the Windows XP/Vista registry and several other systems. When ParFlow runs, it queries the database you have created by key names to get the values you have specified.
The command pfset
is used to create the database entries.
A simple ParFlow input script contains a long list of pfset
commands.
It should be noted that the keys are “dynamic” in that many are built up from values of other keys. For example if you have two wells named northwell and southwell then you will have to set some keys which specify the parameters for each well. The keys are built up in a simple sort of heirarchy.
The following sections contain a description of all of the keys used by
ParFlow. For an example of input files you can look at the test
subdirectory
of the ParFlow distribution. Looking over some examples should give you
a good feel for how the file scripts are put together.
Each key’s entry has the form:
type KeyName default value Description
The “type” is one of integer, double, string, list. Integer and double are IEEE numbers. String is a text string (for example, a filename). Strings can contain spaces if you use the proper TCL syntax (i.e. using double quotes). These types are standard TCL types. Lists are strings but they indicate the names of a series of items. For example you might need to specify the names of the geometries. You would do this using space seperated names (what we are calling a list) “layer1 layer2 layer3”.
The descriptions that follow are organized into functional areas. An example for each database entry is given.
Note that units used for each physical quantity specified in the input file must be consistent with units used for all other quantities. The exact units used can be any consistent set as ParFlow does not assume any specific set of units. However, it is up to the user to make sure all specifications are indeed consistent.
6.1.1. Input File Format Number
integer FileVersion no default This gives the value of the input file version number that this file fits.
pfset FileVersion 4
As development of the ParFlow code continues, the input file format will
vary. We have thus included an input file format number as a way of
verifying that the correct format type is being used. The user can check
in the parflow/config/file_versions.h
file to verify that the format
number specified in the input file matches the defined value
of PFIN_VERSION
.
6.1.2. Computing Topology
This section describes how processors are assigned in order to solve the domain in parallel. “P” allocates the number of processes to the grid-cells in x. “Q” allocates the number of processes to the grid-cells in y. “R” allocates the number of processes to the grid-cells in z. Please note “R” should always be 1 if you are running with Solver Richards [Jones-Woodward01] unless you’re running a totally saturated domain (solver IMPES).
integer Process.Topology.P no default This assigns the process splits in the x direction.
pfset Process.Topology.P 2
integer Process.Topology.Q no default This assigns the process splits in the y direction.
pfset Process.Topology.Q 1
integer Process.Topology.P no default This assigns the process splits in the z direction.
pfset Process.Topology.R 1
In addition, you can assign the computing topology when you initiate your parflow script using tcl. You must include the topology allocation when using tclsh and the parflow script.
Example Usage:
[from Terminal] tclsh default_single.tcl 2 1 1
[At the top of default_single.tcl you must include the following]
set NP [lindex $argv 0]
set NQ [lindex $argv 1]
pfset Process.Topology.P $NP
pfset Process.Topology.Q $NQ
pfset Process.Topology.R 1
6.1.3. Computational Grid
The computational grid is briefly described in §3.1 Defining the Problem. The computational grid keys set the bottom left corner of the domain to a specific point in space. If using a .pfsol file, the bottom left corner location of the .pfsol file must be the points designated in the computational grid. The user can also assign the x, y and z location to correspond to a specific coordinate system (i.e. UTM).
double ComputationalGrid.Lower.X no default This assigns the lower x coordinate location for the computational grid.
pfset ComputationalGrid.Lower.X 0.0
double ComputationalGrid.Lower.Y no default This assigns the lower y coordinate location for the computational grid.
pfset ComputationalGrid.Lower.Y 0.0
double ComputationalGrid.Lower.Z no default This assigns the lower z coordinate location for the computational grid.
pfset ComputationalGrid.Lower.Z 0.0
integer ComputationalGrid.NX no default This assigns the number of grid cells in the x direction for the computational grid.
pfset ComputationalGrid.NX 10
integer ComputationalGrid.NY no default This assigns the number of grid cells in the y direction for the computational grid.
pfset ComputationalGrid.NY 10
integer ComputationalGrid.NZ no default This assigns the number of grid cells in the z direction for the computational grid.
pfset ComputationalGrid.NZ 10
real ComputationalGrid.DX no default This defines the size of grid cells in the x direction. Units are L and are defined by the units of the hydraulic conductivity used in the problem.
pfset ComputationalGrid.DX 10.0
real ComputationalGrid.DY no default This defines the size of grid cells in the y direction. Units are L and are defined by the units of the hydraulic conductivity used in the problem.
pfset ComputationalGrid.DY 10.0
real ComputationalGrid.DZ no default This defines the size of grid cells in the z direction. Units are L and are defined by the units of the hydraulic conductivity used in the problem.
pfset ComputationalGrid.DZ 1.0
Example Usage:
#---------------------------------------------------------
# Computational Grid
#---------------------------------------------------------
pfset ComputationalGrid.Lower.X -10.0
pfset ComputationalGrid.Lower.Y 10.0
pfset ComputationalGrid.Lower.Z 1.0
pfset ComputationalGrid.NX 18
pfset ComputationalGrid.NY 18
pfset ComputationalGrid.NZ 8
pfset ComputationalGrid.DX 8.0
pfset ComputationalGrid.DY 10.0
pfset ComputationalGrid.DZ 1.0
6.1.4. Geometries
Here we define all “geometrical” information needed by ParFlow. For example, the domain (and patches on the domain where boundary conditions are to be imposed), lithology or hydrostratigraphic units, faults, initial plume shapes, and so on, are considered geometries.
This input section is a little confusing. Two items are being specified, geometry inputs and geometries. A geometry input is a type of geometry input (for example a box or an input file). A geometry input can contain more than one geometry. A geometry input of type Box has a single geometry (the square box defined by the extants of the two points). A SolidFile input type can contain several geometries.
list GeomInput.Names no default This is a list of the geometry input names which define the containers for all of the geometries defined for this problem.
pfset GeomInput.Names "solidinput indinput boxinput"
string GeomInput.*geom_input_name*.InputType no default This defines the input type for the geometry input with geom_input_name. This key must be one of: SolidFile, IndicatorField, Box.
pfset GeomInput.solidinput.InputType SolidFile
list GeomInput.*geom_input_name*.GeomNames no default This is a list of the names of the geometries defined by the geometry input. For a geometry input type of Box, the list should contain a single geometry name. For the SolidFile geometry type this should contain a list with the same number of gemetries as were defined using GMS. The order of geometries in the SolidFile should match the names. For IndicatorField types you need to specify the value in the input field which matches the name using GeomInput.geom_input_name.Value.
pfset GeomInput.solidinput.GeomNames "domain bottomlayer \
middlelayer toplayer"
string GeomInput.*geom_input_name*.Filename no default For IndicatorField and SolidFile geometry inputs this key specifies the input filename which contains the field or solid information.
pfset GeomInput.solidinput.FileName ocwd.pfsol
integer GeomInput.*geometry_input_name*.Value no default For IndicatorField geometry inputs you need to specify the mapping between values in the input file and the geometry names. The named geometry will be defined whereever the input file is equal to the specifed value.
pfset GeomInput.sourceregion.Value 11
For box geometries you need to specify the location of the box. This is done by defining two corners of the the box.
double Geom.*box_geom_name*.Lower.X no default This gives the lower X real space coordinate value of the previously specified box geometry of name box_geom_name.
pfset Geom.background.Lower.X -1.0
double Geom.*box_geom_name*.Lower.Y no default This gives the lower Y real space coordinate value of the previously specified box geometry of name box_geom_name.
pfset Geom.background.Lower.Y -1.0
double Geom.*box_geom_name*.Lower.Z no default This gives the lower Z real space coordinate value of the previously specified box geometry of name box_geom_name.
pfset Geom.background.Lower.Z -1.0
double Geom.*box_geom_name*.Upper.X no default This gives the upper X real space coordinate value of the previously specified box geometry of name box_geom_name.
pfset Geom.background.Upper.X 151.0
double Geom.*box_geom_name*.Upper.Y no default This gives the upper Y real space coordinate value of the previously specified box geometry of name box_geom_name.
pfset Geom.background.Upper.Y 171.0
double Geom.*box_geom_name*.Upper.Z no default This gives the upper Z real space coordinate value of the previously specified box geometry of name box_geom_name.
pfset Geom.background.Upper.Z 11.0
list Geom.*geom_name*.Patches no default Patches are defined on the surfaces of geometries. Currently you can only define patches on Box geometries and on the the first geometry in a SolidFile. For a Box the order is fixed (left right front back bottom top) but you can name the sides anything you want.
For SolidFiles the order is printed by the conversion routine that converts GMS to SolidFile format.
pfset Geom.background.Patches "left right front back bottom top"
Here is an example geometry input section which has three geometry inputs.
#---------------------------------------------------------
# The Names of the GeomInputs
#---------------------------------------------------------
pfset GeomInput.Names "solidinput indinput boxinput"
#
# For a solid file geometry input type you need to specify the names
# of the gemetries and the filename
#
pfset GeomInput.solidinput.InputType SolidFile
# The names of the geometries contained in the solid file. Order is
# important and defines the mapping. First geometry gets the first name.
pfset GeomInput.solidinput.GeomNames "domain"
#
# Filename that contains the geometry
#
pfset GeomInput.solidinput.FileName ocwd.pfsol
#
# An indicator field is a 3D field of values.
# The values within the field can be mapped
# to ParFlow geometries. Indicator fields must match the
# computation grid exactly!
#
pfset GeomInput.indinput.InputType IndicatorField
pfset GeomInput.indinput.GeomNames “sourceregion concenregion”
pfset GeomInput.indinput.FileName ocwd.pfb
#
# Within the indicator.pfb file, assign the values to each GeomNames
#
pfset GeomInput.sourceregion.Value 11
pfset GeomInput.concenregion.Value 12
#
# A box is just a box defined by two points.
#
pfset GeomInput.boxinput.InputType Box
pfset GeomInput.boxinput.GeomName background
pfset Geom.background.Lower.X -1.0
pfset Geom.background.Lower.Y -1.0
pfset Geom.background.Lower.Z -1.0
pfset Geom.background.Upper.X 151.0
pfset Geom.background.Upper.Y 171.0
pfset Geom.background.Upper.Z 11.0
#
# The patch order is fixed in the .pfsol file, but you
# can call the patch name anything you
# want (i.e. left right front back bottom top)
#
pfset Geom.domain.Patches " z-upper x-lower y-lower \
x-upper y-upper z-lower"
6.1.5. Timing Information
The data given in the timing section describe all the “temporal” information needed by ParFlow. The data items are used to describe time units for later sections, sequence iterations in time, indicate actual starting and stopping values and give instructions on when data is printed out.
double TimingInfo.BaseUnit no default This key is used to indicate the base unit of time for entering time values. All time should be expressed as a multiple of this value. This should be set to the smallest interval of time to be used in the problem. For example, a base unit of “1” means that all times will be integer valued. A base unit of “0.5” would allow integers and fractions of 0.5 to be used for time input values.
The rationale behind this restriction is to allow time to be discretized on some interval to enable integer arithmetic to be used when computing/comparing times. This avoids the problems associated with real value comparisons which can lead to events occurring at different timesteps on different architectures or compilers.
This value is also used when describing “time cycling data” in, currently, the well and boundary condition sections. The lengths of the cycles in those sections will be integer multiples of this value, therefore it needs to be the smallest divisor which produces an integral result for every “real time” cycle interval length needed.
pfset TimingInfo.BaseUnit 1.0
integer TimingInfo.StartCount no default This key is used to indicate the time step number that will be associated with the first advection cycle in a transient problem. The value -1 indicates that advection is not to be done. The value 0 indicates that advection should begin with the given initial conditions. Values greater than 0 are intended to mean “restart” from some previous “checkpoint” time-step, but this has not yet been implemented.
pfset TimingInfo.StartCount 0
double TimingInfo.StartTime no default This key is used to indicate the starting time for the simulation.
pfset TimingInfo.StartTime 0.0
double TimingInfo.StopTime no default This key is used to indicate the stopping time for the simulation.
pfset TimingInfo.StopTime 100.0
double TimingInfo.DumpInterval no default This key is the real time interval at which time-dependent output should be written. A value of 0 will produce undefined behavior. If the value is negative, output will be dumped out every \(n\) time steps, where \(n\) is the absolute value of the integer part of the value.
pfset TimingInfo.DumpInterval 10.0
integer TimingInfo.DumpIntervalExecutionTimeLimit 0 This key is used to indicate a wall clock time to halt the execution of a run. At the end of each dump interval the time remaining in the batch job is compared with the user supplied value, if remaining time is less than or equal to the supplied value the execution is halted. Typically used when running on batch systems with time limits to force a clean shutdown near the end of the batch job. Time units is seconds, a value of 0 (the default) disables the check.
Currently only supported on SLURM based systems, “–with-slurm” must be specified at configure time to enable.
pfset TimingInfo.DumpIntervalExecutionTimeLimit 360
For Richards’ equation cases only input is collected for time step selection. Input for this section is given as follows:
list TimeStep.Type no default This key must be one of: Constant or Growth. The value Constant defines a constant time step. The value Growth defines a time step that starts as \(dt_0\) and is defined for other steps as \(dt^{new} = \gamma dt^{old}\) such that \(dt^{new} \leq dt_{max}\) and \(dt^{new} \geq dt_{min}\).
pfset TimeStep.Type Constant
double TimeStep.Value no default This key is used only if a constant time step is selected and indicates the value of the time step for all steps taken.
pfset TimeStep.Value 0.001
double TimeStep.InitialStep no default This key specifies the initial time step \(dt_0\) if the Growth type time step is selected.
pfset TimeStep.InitialStep 0.001
double TimeStep.GrowthFactor no default This key specifies the growth factor \(\gamma\) by which a time step will be multiplied to get the new time step when the Growth type time step is selected.
pfset TimeStep.GrowthFactor 1.5
double TimeStep.MaxStep no default This key specifies the maximum time step allowed, \(dt_{max}\), when the Growth type time step is selected.
pfset TimeStep.MaxStep 86400
double TimeStep.MinStep no default This key specifies the minimum time step allowed, \(dt_{min}\), when the Growth type time step is selected.
pfset TimeStep.MinStep 1.0e-3
Here is a detailed example of how timing keys might be used in a simualtion.
#-----------------------------------------------------------------------------
# Setup timing info [hr]
# 8760 hours in a year. Dumping files every 24 hours. Hourly timestep
#-----------------------------------------------------------------------------
pfset TimingInfo.BaseUnit 1.0
pfset TimingInfo.StartCount 0
pfset TimingInfo.StartTime 0.0
pfset TimingInfo.StopTime 8760.0
pfset TimingInfo.DumpInterval -24
## Timing constant example
pfset TimeStep.Type Constant
pfset TimeStep.Value 1.0
## Timing growth example
pfset TimeStep.Type Growth
pfset TimeStep.InitialStep 0.0001
TimeStep.GrowthFactor 1.4
TimeStep.MaxStep 1.0
TimeStep.MinStep 0.0001
6.1.6. Time Cycles
The data given in the time cycle section describe how time intervals are created and named to be used for time-dependent boundary and well information needed by ParFlow. All the time cycles are synched to the TimingInfo.BaseUnit key described above and are integer multipliers of that value.
list CycleNames no default This key is used to specify the named time cycles to be used in a simulation. It is a list of names and each name defines a time cycle and the number of items determines the total number of time cycles specified. Each named cycle is described using a number of keys defined below.
pfset Cycle.Names constant onoff
list Cycle.*cycle_name*.Names no default This key is used to specify the named time intervals for each cycle. It is a list of names and each name defines a time interval when a specific boundary condition is applied and the number of items determines the total number of intervals in that time cycle.
pfset Cycle.onoff.Names "on off"
integer Cycle.*cycle_name.interval_name*.Length no default This key is used to specify the length of a named time intervals. It is an integer multiplier of the value set for the TimingInfo.BaseUnit key described above. The total length of a given time cycle is the sum of all the intervals multiplied by the base unit.
pfset Cycle.onoff.on.Length 10
integer Cycle.*cycle_name*.Repeat no default This key is used to specify the how many times a named time interval repeats. A positive value specifies a number of repeat cycles a value of -1 specifies that the cycle repeat for the entire simulation.
pfset Cycle.onoff.Repeat -1
Here is a detailed example of how time cycles might be used in a simualtion.
#-----------------------------------------------------------------------------
# Time Cycles
#-----------------------------------------------------------------------------
pfset Cycle.Names "constant rainrec"
pfset Cycle.constant.Names "alltime"
pfset Cycle.constant.alltime.Length 8760
pfset Cycle.constant.Repeat -1
# Creating a rain and recession period for the rest of year
pfset Cycle.rainrec.Names "rain rec"
pfset Cycle.rainrec.rain.Length 10
pfset Cycle.rainrec.rec.Length 8750
pfset Cycle.rainrec.Repeat -1
6.1.7. Domain
The domain may be represented by any of the solid types in §6.1.4 Geometries above that allow the definition of surface patches. These surface patches are used to define boundary conditions in §6.1.24 Boundary Conditions: Pressure and §6.1.25 Boundary Conditions: Saturation below. Subsequently, it is required that the union (or combination) of the defined surface patches equal the entire domain surface. NOTE: This requirement is NOT checked in the code.
string Domain.GeomName no default This key specifies which of the named geometries is the problem domain.
pfset Domain.GeomName domain
6.1.8. Phases and Contaminants
list Phase.Names no default This specifies the names of phases to be modeled. Currently only 1 or 2 phases may be modeled.
pfset Phase.Names "water"
list Contaminant.Names no default This specifies the names of contaminants to be advected.
pfset Contaminants.Names "tce"
6.1.9. Gravity, Phase Density and Phase Viscosity
double Gravity no default Specifies the gravity constant to be used.
pfset Gravity 1.0
string Phase.*phase_name*.Density.Type no default This key specifies whether density will be a constant value or if it will be given by an equation of state of the form \((rd)exp(cP)\), where \(P\) is pressure, \(rd\) is the density at atmospheric pressure, and \(c\) is the phase compressibility constant. This key must be either Constant or EquationOfState.
pfset Phase.water.Density.Type Constant
double Phase.*phase_name*.Density.Value no default This specifies the value of density if this phase was specified to have a constant density value for the phase phase_name.
pfset Phase.water.Density.Value 1.0
double Phase.*phase_name*.Density.ReferenceDensity no default This key specifies the reference density if an equation of state density function is specified for the phase phase_name.
pfset Phase.water.Density.ReferenceDensity 1.0
double Phase.*phase_name*.Density.CompressibilityConstant no default This key specifies the phase compressibility constant if an equation of state density function is specified for the phase phase|-name.
pfset Phase.water.Density.CompressibilityConstant 1.0
string Phase.*phase_name*.Viscosity.Type Constant This key specifies whether viscosity will be a constant value. Currently, the only choice for this key is Constant.
pfset Phase.water.Viscosity.Type Constant
double Phase.*phase_name*.Viscosity.Value no default This specifies the value of viscosity if this phase was specified to have a constant viscosity value.
pfset Phase.water.Viscosity.Value 1.0
6.1.10. Chemical Reactions
double Contaminants.*contaminant_name*.Degradation.Value no default This key specifies the half-life decay rate of the named contaminant, contaminant_name. At present only first order decay reactions are implemented and it is assumed that one contaminant cannot decay into another.
pfset Contaminants.tce.Degradation.Value 0.0
6.1.11. Permeability
In this section, permeability property values are assigned to grid points within geometries (specified in §6.1.4 Geometries above) using one of the methods described below. Permeabilities are assumed to be a diagonal tensor with entries given as,
where \(K({\bf x})\) is the permeability field given below. Specification of the tensor entries (\(k_x, k_y\) and \(k_z\)) will be given at the end of this section.
The random field routines (turning bands and pgs) can use conditioning data if the user so desires. It is not necessary to use conditioning as ParFlow automatically defaults to not use conditioning data, but if conditioning is desired, the following key should be set:
string Perm.Conditioning.FileName “NA” This key specifies the name of the file that contains the conditioning data. The default string NA indicates that conditioning data is not applicable.
pfset Perm.Conditioning.FileName "well_cond.txt"
The file that contains the conditioning data is a simple ascii file containing points and values. The format is:
nlines
x1 y1 z1 value1
x2 y2 z2 value2
. . . .
. . . .
. . . .
xn yn zn valuen
The value of nlines is just the number of lines to follow in the file, which is equal to the number of data points.
The variables xi,yi,zi are the real space coordinates (in the units used for the given parflow run) of a point at which a fixed permeability value is to be assigned. The variable valuei is the actual permeability value that is known.
Note that the coordinates are not related to the grid in any way. Conditioning does not require that fixed values be on a grid. The PGS algorithm will map the given value to the closest grid point and that will be fixed. This is done for speed reasons. The conditioned turning bands algorithm does not do this; conditioning is done for every grid point using the given conditioning data at the location given. Mapping to grid points for that algorithm does not give any speedup, so there is no need to do it.
NOTE: The given values should be the actual measured values - adjustment in the conditioning for the lognormal distribution that is assumed is taken care of in the algorithms.
The general format for the permeability input is as follows:
list Geom.Perm.Names no default This key specifies all of the geometries to which a permeability field will be assigned. These geometries must cover the entire computational domain.
pfset GeomInput.Names "background domain concen_region"
string Geom.geometry_name.Perm.Type no default This key specifies which method is to be used to assign permeability data to the named geometry, geometry_name. It must be either Constant, TurnBands, ParGuass, or PFBFile. The Constant value indicates that a constant is to be assigned to all grid cells within a geometry. The TurnBand value indicates that Tompson’s Turning Bands method is to be used to assign permeability data to all grid cells within a geometry [TAG89]. The ParGauss value indicates that a Parallel Gaussian Simulator method is to be used to assign permeability data to all grid cells within a geometry. The PFBFile value indicates that premeabilities are to be read from the “ParFlow Binary” file. Both the Turning Bands and Parallel Gaussian Simulators generate a random field with correlation lengths in the \(3\) spatial directions given by \(\lambda_x\), \(\lambda_y\), and \(\lambda_z\) with the geometric mean of the log normal field given by \(\mu\) and the standard deviation of the normal field given by \(\sigma\). In generating the field both of these methods can be made to stratify the data, that is follow the top or bottom surface. The generated field can also be made so that the data is normal or log normal, with or without bounds truncation. Turning Bands uses a line process, the number of lines used and the resolution of the process can be changed as well as the maximum normalized frequency \(K_{\rm max}\) and the normalized frequency increment \(\delta K\). The Parallel Gaussian Simulator uses a search neighborhood, the number of simulated points and the number of conditioning points can be changed.
pfset Geom.background.Perm.Type Constant
double Geom.*geometry_name*.Perm.Value no default This key specifies the value assigned to all points in the named geometry, geometry_name, if the type was set to constant.
pfset Geom.domain.Perm.Value 1.0
double Geom.*geometry_name*.Perm.LambdaX no default This key specifies the x correlation length, \(\lambda_x\), of the field generated for the named geometry, geometry_name, if either the Turning Bands or Parallel Gaussian Simulator are chosen.
pfset Geom.domain.Perm.LambdaX 200.0
double Geom.*geometry_name*.Perm.LambdaY no default This key specifies the y correlation length, \(\lambda_y\), of the field generated for the named geometry, geometry_name, if either the Turning Bands or Parallel Gaussian Simulator are chosen.
pfset Geom.domain.Perm.LambdaY 200.0
double Geom.*geometry_name*.Perm.LambdaZ no default This key specifies the z correlation length, \(\lambda_z\), of the field generated for the named geometry, geometry_name, if either the Turning Bands or Parallel Gaussian Simulator are chosen.
pfset Geom.domain.Perm.LambdaZ 10.0
double Geom.*geometry_name*.Perm.GeomMean no default This key specifies the geometric mean, \(\mu\), of the log normal field generated for the named geometry, geometry_name, if either the Turning Bands or Parallel Gaussian Simulator are chosen.
pfset Geom.domain.Perm.GeomMean 4.56
double Geom.*geometry_name*.Perm.Sigma no default This key specifies the standard deviation, \(\sigma\), of the normal field generated for the named geometry, geometry_name, if either the Turning Bands or Parallel Gaussian Simulator are chosen.
pfset Geom.domain.Perm.Sigma 2.08
integer Geom.*geometry_name*.Perm.Seed 1 This key specifies the initial seed for the random number generator used to generate the field for the named geometry, geometry_name, if either the Turning Bands or Parallel Gaussian Simulator are chosen. This number must be positive.
pfset Geom.domain.Perm.Seed 1
integer Geom.*geometry_name*.Perm.NumLines 100 This key specifies the number of lines to be used in the Turning Bands algorithm for the named geometry, geometry_name.
pfset Geom.domain.Perm.NumLines 100
double Geom.*geometry_name*.Perm.RZeta 5.0 This key specifies the resolution of the line processes, in terms of the minimum grid spacing, to be used in the Turning Bands algorithm for the named geometry, geometry_name. Large values imply high resolution.
pfset Geom.domain.Perm.RZeta 5.0
double Geom.*geometry_name*.Perm.KMax 100.0 This key specifies the the maximum normalized frequency, \(K_{\rm max}\), to be used in the Turning Bands algorithm for the named geometry, geometry_name.
pfset Geom.domain.Perm.KMax 100.0
double Geom.*geometry_name*.Perm.DelK 0.2 This key specifies the normalized frequency increment, \(\delta K\), to be used in the Turning Bands algorithm for the named geometry, geometry_name.
pfset Geom.domain.Perm.DelK 0.2
integer Geom.*geometry_name*.Perm.MaxNPts no default This key sets limits on the number of simulated points in the search neighborhood to be used in the Parallel Gaussian Simulator for the named geometry, geometry_name.
pfset Geom.domain.Perm.MaxNPts 5
integer Geom.*geometry_name*.Perm.MaxCpts no default This key sets limits on the number of external conditioning points in the search neighborhood to be used in the Parallel Gaussian Simulator for the named geometry, geometry_name.
pfset Geom.domain.Perm.MaxCpts 200
string Geom.*geometry_name*.Perm.LogNormal “LogTruncated” The key specifies when a normal, log normal, truncated normal or truncated log normal field is to be generated by the method for the named geometry, geometry_name. This value must be one of Normal, Log, NormalTruncated or LogTruncate and can be used with either Turning Bands or the Parallel Gaussian Simulator.
pfset Geom.domain.Perm.LogNormal "LogTruncated"
string Geom.*geometry_name*.Perm.StratType “Bottom” This key specifies the stratification of the permeability field generated by the method for the named geometry, geometry_name. The value must be one of Horizontal, Bottom or Top and can be used with either the Turning Bands or the Parallel Gaussian Simulator.
pfset Geom.domain.Perm.StratType "Bottom"
double Geom.*geometry_name*.Perm.LowCutoff no default This key specifies the low cutoff value for truncating the generated field for the named geometry, geometry_name, when either the NormalTruncated or LogTruncated values are chosen.
pfset Geom.domain.Perm.LowCutoff 0.0
double Geom.*geometry_name*.Perm.HighCutoff no default This key specifies the high cutoff value for truncating the generated field for the named geometry, geometry_name, when either the NormalTruncated or LogTruncated values are chosen.
pfset Geom.domain.Perm.HighCutoff 100.0
string Geom.*geometry_name*.Perm.FileName no default This key specifies that permeability values for the specified geometry, geometry_name, are given according to a user-supplied description in the “ParFlow Binary” file whose filename is given as the value. For a description of the ParFlow Binary file format, see §6.3 ParFlow Binary Files (.pfb). The ParFlow Binary file associated with the named geometry must contain a collection of permeability values corresponding in a one-to-one manner to the entire computational grid. That is to say, when the contents of the file are read into the simulator, a complete permeability description for the entire domain is supplied. Only those values associated with computational cells residing within the geometry (as it is represented on the computational grid) will be copied into data structures used during the course of a simulation. Thus, the values associated with cells outside of the geounit are irrelevant. For clarity, consider a couple of different scenarios. For example, the user may create a file for each geometry such that appropriate permeability values are given for the geometry and “garbage” values (e.g., some flag value) are given for the rest of the computational domain. In this case, a separate binary file is specified for each geometry. Alternatively, one may place all values representing the permeability field on the union of the geometries into a single binary file. Note that the permeability values must be represented in precisely the same configuration as the computational grid. Then, the same file could be specified for each geounit in the input file. Or, the computational domain could be described as a single geouint (in the ParFlow input file) in which case the permeability values would be read in only once.
pfset Geom.domain.Perm.FileName "domain_perm.pfb"
string Perm.TensorType no default This key specifies whether the permeability tensor entries \(k_x, k_y\) and \(k_z\) will be specified as three constants within a set of regions covering the domain or whether the entries will be specified cell-wise by files. The choices for this key are TensorByGeom and TensorByFile.
pfset Perm.TensorType TensorByGeom
string Geom.Perm.TensorByGeom.Names no default This key specifies all of the geometries to which permeability tensor entries will be assigned. These geometries must cover the entire computational domain.
pfset Geom.Perm.TensorByGeom.Names "background domain"
double Geom.*geometry_name*.Perm.TensorValX no default This key specifies the value of \(k_x\) for the geometry given by geometry_name.
pfset Geom.domain.Perm.TensorValX 1.0
double Geom.*geometry_name*.Perm.TensorValY no default This key specifies the value of \(k_y\) for the geometry given by geom_name.
pfset Geom.domain.Perm.TensorValY 1.0
double Geom.*geometry_name*.Perm.TensorValZ no default This key specifies the value of \(k_z\) for the geometry given by geom_name.
pfset Geom.domain.Perm.TensorValZ 1.0
string Geom.*geometry_name*.Perm.TensorFileX no default This key specifies that \(k_x\) values for the specified geometry, geometry_name, are given according to a user-supplied description in the “ParFlow Binary” file whose filename is given as the value. The only choice for the value of geometry_name is “domain”.
pfset Geom.domain.Perm.TensorByFileX "perm_x.pfb"
string Geom.*geometry_name*.Perm.TensorFileY no default This key specifies that \(k_y\) values for the specified geometry, geometry_name, are given according to a user-supplied description in the “ParFlow Binary” file whose filename is given as the value. The only choice for the value of geometry_name is “domain”.
pfset Geom.domain.Perm.TensorByFileY "perm_y.pfb"
string Geom.*geometry_name*.Perm.TensorFileZ no default This key specifies that \(k_z\) values for the specified geometry, geometry_name, are given according to a user-supplied description in the “ParFlow Binary” file whose filename is given as the value. The only choice for the value of geometry_name is “domain”.
pfset Geom.domain.Perm.TensorByFileZ "perm_z.pfb"
6.1.12. Porosity
Here, porosity values are assigned within geounits (specified in §6.1.4 Geometries above) using one of the methods described below.
The format for this section of input is:
list Geom.Porosity.GeomNames no default This key specifies all of the geometries on which a porosity will be assigned. These geometries must cover the entire computational domain.
pfset Geom.Porosity.GeomNames "background"
string Geom.*geometry_name*.Porosity.Type no default This key specifies which method is to be used to assign porosity data to the named geometry, geometry_name. The only choice currently available is Constant which indicates that a constant is to be assigned to all grid cells within a geometry.
pfset Geom.background.Porosity.Type Constant
double Geom.*geometry_name*.Porosity.Value no default This key specifies the value assigned to all points in the named geometry, geometry_name, if the type was set to constant.
pfset Geom.domain.Porosity.Value 1.0
6.1.13. Specific Storage
Here, specific storage (\(S_s\) in Equation [eq:richard]) values are assigned within geounits (specified in §6.1.4 Geometries above) using one of the methods described below.
The format for this section of input is:
list Specific Storage.GeomNames no default This key specifies all of the geometries on which a different specific storage value will be assigned. These geometries must cover the entire computational domain.
pfset SpecificStorage.GeomNames "domain"
string SpecificStorage.Type no default This key specifies which method is to be used to assign specific storage data. The only choice currently available is Constant which indicates that a constant is to be assigned to all grid cells within a geometry.
pfset SpecificStorage.Type Constant
double Geom.*geometry_name*.SpecificStorage.Value no default This key specifies the value assigned to all points in the named geometry, geometry_name, if the type was set to constant.
pfset Geom.domain.SpecificStorage.Value 1.0e-4
6.1.14. dZMultipliers
Here, dZ multipliers (\(\delta Z * m\)) values are assigned within geounits (specified in §6.1.4 Geometries above) using one of the methods described below.
The format for this section of input is:
string Solver.Nonlinear.VariableDz False This key specifies whether dZ multipliers are to be used, the default is False. The default indicates a false or non-active variable dz and each layer thickness is 1.0 [L].
pfset Solver.Nonlinear.VariableDz True
list dzScale.GeomNames no default This key specifies which problem domain is being applied a variable dz subsurface. These geometries must cover the entire computational domain.
pfset dzScale.GeomNames domain
string dzScale.Type no default This key specifies which method is to be used to assign variable vertical grid spacing. The choices currently available are Constant which indicates that a constant is to be assigned to all grid cells within a geometry, nzList which assigns all layers of a given model to a list value, and PFBFile which reads in values from a distributed pfb file.
pfset dzScale.Type Constant
list Specific dzScale.GeomNames no default This key specifies all of the geometries on which a different dz scaling value will be assigned. These geometries must cover the entire computational domain.
pfset dzScale.GeomNames "domain"
double Geom.*geometry_name*.dzScale.Value no default This key specifies the value assigned to all points in the named geometry, geometry_name, if the type was set to constant.
pfset Geom.domain.dzScale.Value 1.0
string Geom.*geometry_name*.dzScale.FileName no default This key specifies file to be read in for variable dz values for the given geometry, geometry_name, if the type was set to PFBFile.
pfset Geom.domain.dzScale.FileName vardz.pfb
integer dzScale.nzListNumber no default This key indicates the number of layers with variable dz in the subsurface. This value is the same as the ComputationalGrid.NZ key.
pfset dzScale.nzListNumber 10
double Cell.*nzListNumber*.dzScale.Value no default This key assigns the thickness of each layer defined by nzListNumber. ParFlow assigns the layers from the bottom-up (i.e. the bottom of the domain is layer 0, the top is layer NZ-1). The total domain depth (Geom.domain.Upper.Z) does not change with variable dz. The layer thickness is calculated by ComputationalGrid.DZ *dZScale.
pfset Cell.0.dzScale.Value 1.0
Example Usage:
#--------------------------------------------
# Variable dz Assignments
#------------------------------------------
# Set VariableDz to be true
# Indicate number of layers (nzlistnumber), which is the same as nz
# (1) There is nz*dz = total depth to allocate,
# (2) Each layer’s thickness is dz*dzScale, and
# (3) Assign the layer thickness from the bottom up.
# In this example nz = 5; dz = 10; total depth 40;
# Layers Thickness [m]
# 0 15 Bottom layer
# 1 15
# 2 5
# 3 4.5
# 4 0.5 Top layer
pfset Solver.Nonlinear.VariableDz True
pfset dzScale.GeomNames domain
pfset dzScale.Type nzList
pfset dzScale.nzListNumber 5
pfset Cell.0.dzScale.Value 1.5
pfset Cell.1.dzScale.Value 1.5
pfset Cell.2.dzScale.Value 0.5
pfset Cell.3.dzScale.Value 0.45
pfset Cell.4.dzScale.Value 0.05
6.1.15. Manning’s Roughness Values
Here, Manning’s roughness values (\(n\) in Equations [eq:manningsx] and [eq:manningsy]) are assigned to the upper boundary of the domain using one of the methods described below.
The format for this section of input is:
list Mannings.GeomNames no default This key specifies all of the geometries on which a different Mannings roughness value will be assigned. Mannings values may be assigned by PFBFile or as Constant by geometry. These geometries must cover the entire upper surface of the computational domain.
pfset Mannings.GeomNames "domain"
string Mannings.Type no default This key specifies which method is to be used to assign Mannings roughness data. The choices currently available are Constant which indicates that a constant is to be assigned to all grid cells within a geometry and PFBFile which indicates that all values are read in from a distributed, grid-based ParFlow binary file.
pfset Mannings.Type "Constant"
double Mannings.Geom.*geometry_name*.Value no default This key specifies the value assigned to all points in the named geometry, geometry_name, if the type was set to constant.
pfset Mannings.Geom.domain.Value 5.52e-6
double Mannings.FileName no default This key specifies the value assigned to all points be read in from a ParFlow binary file.
pfset Mannings.FileName roughness.pfb
Complete example of setting Mannings roughness \(n\) values by geometry:
pfset Mannings.Type "Constant"
pfset Mannings.GeomNames "domain"
pfset Mannings.Geom.domain.Value 5.52e-6
6.1.16. Topographical Slopes
Here, topographical slope values (\(S_{f,x}\) and \(S_{f,y}\) in Equations [eq:manningsx] and [eq:manningsy]) are assigned to the upper boundary of the domain using one of the methods described below. Note that due to the negative sign in these equations \(S_{f,x}\) and \(S_{f,y}\) take a sign in the direction opposite of the direction of the slope. That is, negative slopes point “downhill” and positive slopes “uphill”.
The format for this section of input is:
list ToposlopesX.GeomNames no default This key specifies all of the geometries on which a different \(x\) topographic slope values will be assigned. Topographic slopes may be assigned by PFBFile or as Constant by geometry. These geometries must cover the entire upper surface of the computational domain.
pfset ToposlopesX.GeomNames "domain"
list ToposlopesY.GeomNames no default This key specifies all of the geometries on which a different \(y\) topographic slope values will be assigned. Topographic slopes may be assigned by PFBFile or as Constant by geometry. These geometries must cover the entire upper surface of the computational domain.
pfset ToposlopesY.GeomNames "domain"
string ToposlopesX.Type no default This key specifies which method is to be used to assign topographic slopes. The choices currently available are Constant which indicates that a constant is to be assigned to all grid cells within a geometry and PFBFile which indicates that all values are read in from a distributed, grid-based ParFlow binary file.
pfset ToposlopesX.Type "Constant"
double ToposlopeX.Geom.*geometry_name*.Value no default This key specifies the value assigned to all points in the named geometry, geometry_name, if the type was set to constant.
pfset ToposlopeX.Geom.domain.Value 0.001
double ToposlopesX.FileName no default This key specifies the value assigned to all points be read in from a ParFlow binary file.
pfset TopoSlopesX.FileName lw.1km.slope_x.pfb
double ToposlopesY.FileName no default This key specifies the value assigned to all points be read in from a ParFlow binary file.
pfset TopoSlopesY.FileName lw.1km.slope_y.pfb
Example of setting \(x\) and \(y\) slopes by geometry:
pfset TopoSlopesX.Type "Constant"
pfset TopoSlopesX.GeomNames "domain"
pfset TopoSlopesX.Geom.domain.Value 0.001
pfset TopoSlopesY.Type "Constant"
pfset TopoSlopesY.GeomNames "domain"
pfset TopoSlopesY.Geom.domain.Value -0.001
Example of setting \(x\) and \(y\) slopes by file:
pfset TopoSlopesX.Type "PFBFile"
pfset TopoSlopesX.GeomNames "domain"
pfset TopoSlopesX.FileName lw.1km.slope_x.pfb
pfset TopoSlopesY.Type "PFBFile"
pfset TopoSlopesY.GeomNames "domain"
pfset TopoSlopesY.FileName lw.1km.slope_y.pfb
6.1.17. Retardation
Here, retardation values are assigned for contaminants within geounits (specified in §6.1.4 Geometries above) using one of the functions described below. The format for this section of input is:
list Geom.Retardation.GeomNames no default This key specifies all of the geometries to which the contaminants will have a retardation function applied.
pfset GeomInput.Names "background"
string Geom.*geometry_name*.*contaminant_name*.Retardation.Type no default This key specifies which function is to be used to compute the retardation for the named contaminant, contaminant_name, in the named geometry, geometry_name. The only choice currently available is Linear which indicates that a simple linear retardation function is to be used to compute the retardation.
pfset Geom.background.tce.Retardation.Type Linear
double Geom.*geometry_name*.*contaminant_name*.Retardation.Value no default This key specifies the distribution coefficient for the linear function used to compute the retardation of the named contaminant, contaminant_name, in the named geometry, geometry_name. The value should be scaled by the density of the material in the geometry.
pfset Geom.domain.Retardation.Value 0.2
6.1.18. Full Multiphase Mobilities
Here we define phase mobilities by specifying the relative permeability function. Input is specified differently depending on what problem is being specified. For full multi-phase problems, the following input keys are used. See the next section for the correct Richards’ equation input format.
string Phase.*phase_name*.Mobility.Type no default This key specifies whether the mobility for phase_name will be a given constant or a polynomial of the form, \((S - S_0)^{a}\), where \(S\) is saturation, \(S_0\) is irreducible saturation, and \(a\) is some exponent. The possibilities for this key are Constant and Polynomial.
pfset Phase.water.Mobility.Type Constant
double Phase.*phase_name*.Mobility.Value no default This key specifies the constant mobility value for phase phase_name.
pfset Phase.water.Mobility.Value 1.0
double Phase.*phase_name*.Mobility.Exponent 2.0 This key specifies the exponent used in a polynomial representation of the relative permeability. Currently, only a value of \(2.0\) is allowed for this key.
pfset Phase.water.Mobility.Exponent 2.0
double Phase.*phase_name*.Mobility.IrreducibleSaturation 0.0 This key specifies the irreducible saturation used in a polynomial representation of the relative permeability. Currently, only a value of 0.0 is allowed for this key.
pfset Phase.water.Mobility.IrreducibleSaturation 0.0
6.1.19. Richards’ Equation Relative Permeabilities
The following keys are used to describe relative permeability input for the Richards’ equation implementation. They will be ignored if a full two-phase formulation is used.
string Phase.RelPerm.Type no default This key specifies the type of relative permeability function that will be used on all specified geometries. Note that only one type of relative permeability may be used for the entire problem. However, parameters may be different for that type in different geometries. For instance, if the problem consists of three geometries, then VanGenuchten may be specified with three different sets of parameters for the three different goemetries. However, once VanGenuchten is specified, one geometry cannot later be specified to have Data as its relative permeability. The possible values for this key are Constant, VanGenuchten, Haverkamp, Data, and Polynomial.
pfset Phase.RelPerm.Type Constant
The various possible functions are defined as follows. The Constant specification means that the relative permeability will be constant on the specified geounit. The VanGenuchten specification means that the relative permeability will be given as a Van Genuchten function [VanGenuchten80] with the form,
where \(\alpha\) and \(n\) are soil parameters and \(m = 1 - 1/n\), on each region. The Haverkamp specification means that the relative permeability will be given in the following form [Haverkamp-Vauclin81],
where \(A\) and \(\gamma\) are soil parameters, on each region. The Data specification is currently unsupported but will later mean that data points for the relative permeability curve will be given and ParFlow will set up the proper interpolation coefficients to get values between the given data points. The Polynomial specification defines a polynomial relative permeability function for each region of the form,
list Phase.RelPerm.GeomNames no default This key specifies the geometries on which relative permeability will be given. The union of these geometries must cover the entire computational domain.
pfset Phase.RelPerm.Geonames domain
double Geom.*geom_name*.RelPerm.Value no default This key specifies the constant relative permeability value on the specified geometry.
pfset Geom.domain.RelPerm.Value 0.5
integer Phase.RelPerm.VanGenuchten.File 0 This key specifies whether soil parameters for the VanGenuchten function are specified in a pfb file or by region. The options are either 0 for specification by region, or 1 for specification in a file. Note that either all parameters are specified in files (each has their own input file) or none are specified by files. Parameters specified by files are: \(\alpha\) and N.
pfset Phase.RelPerm.VanGenuchten.File 1
string Geom.*geom_name*.RelPerm.Alpha.Filename no default This key specifies a pfb filename containing the alpha parameters for the VanGenuchten function cell-by-cell. The ONLY option for geom_name is “domain”.
pfset Geom.domain.RelPerm.Alpha.Filename alphas.pfb
string Geom.*geom_name*.RelPerm.N.Filename no default This key specifies a pfb filename containing the N parameters for the VanGenuchten function cell-by-cell. The ONLY option for geom_name is “domain”.
pfset Geom.domain.RelPerm.N.Filename Ns.pfb
double Geom.*geom_name*.RelPerm.Alpha no default This key specifies the \(\alpha\) parameter for the Van Genuchten function specified on geom_name.
pfset Geom.domain.RelPerm.Alpha 0.005
double Geom.*geom_name*.RelPerm.N no default This key specifies the \(N\) parameter for the Van Genuchten function specified on geom_name.
pfset Geom.domain.RelPerm.N 2.0
int Geom.*geom_name*.RelPerm.NumSamplePoints 0 This key specifies the number of sample points for a spline base interpolation table for the Van Genuchten function specified on geom_name. If this number is 0 (the default) then the function is evaluated directly. Using the interpolation table is faster but is less accurate.
pfset Geom.domain.RelPerm.NumSamplePoints 20000
int Geom.*geom_name*.RelPerm.MinPressureHead no default This key specifies the lower value for a spline base interpolation table for the Van Genuchten function specified on geom_name. The upper value of the range is 0. This value is used only when the table lookup method is used (NumSamplePoints is greater than 0).
pfset Geom.domain.RelPerm.MinPressureHead -300
double Geom.*geom_name*.RelPerm.A no default This key specifies the \(A\) parameter for the Haverkamp relative permeability on geom_name.
pfset Geom.domain.RelPerm.A 1.0
double Geom.*geom_name*.RelPerm.Gamma no default This key specifies the the \(\gamma\) parameter for the Haverkamp relative permeability on geom_name.
pfset Geom.domain.RelPerm.Gamma 1.0
integer Geom.*geom_name*.RelPerm.Degree no default This key specifies the degree of the polynomial for the Polynomial relative permeability given on geom_name.
pfset Geom.domain.RelPerm.Degree 1
double Geom.*geom_name*.RelPerm.Coeff.*coeff_number* no default This key specifies the coeff_numberth coefficient of the Polynomial relative permeability given on geom_name.
pfset Geom.domain.RelPerm.Coeff.0 0.5
pfset Geom.domain.RelPerm.Coeff.1 1.0
NOTE: For all these cases, if only one region is to be used (the domain), the background region should NOT be set as that single region. Using the background will prevent the upstream weighting from being correct near Dirichlet boundaries.
6.1.20. Phase Sources
The following keys are used to specify phase source terms. The units of the source term are \(1/T\). So, for example, to specify a region with constant flux rate of \(L^3/T\), one must be careful to convert this rate to the proper units by dividing by the volume of the enclosing region. For Richards’ equation input, the source term must be given as a flux multiplied by density.
string PhaseSources.*phase_name*.Type no default This key specifies the type of source to use for phase phase_name. Possible values for this key are Constant and PredefinedFunction. Constant type phase sources specify a constant phase source value for a given set of regions. PredefinedFunction type phase sources use a preset function (choices are listed below) to specify the source. Note that the PredefinedFunction type can only be used to set a single source over the entire domain and not separate sources over different regions.
pfset PhaseSources.water.Type Constant
list PhaseSources.*phase_name*.GeomNames no default This key specifies the names of the geometries on which source terms will be specified. This is used only for Constant type phase sources. Regions listed later “overlay” regions listed earlier.
pfset PhaseSources.water.GeomNames "bottomlayer middlelayer toplayer"
double PhaseSources.*phase_name*.Geom.*geom_name*.Value no default This key specifies the value of a constant source term applied to phase phase _name on geometry geom_name.
pfset PhaseSources.water.Geom.toplayer.Value 1.0
string PhaseSources.*phase_name*.PredefinedFunction no default This key specifies which of the predefined functions will be used for the source. Possible values for this key are X, XPlusYPlusZ, X3Y2PlusSinXYPlus1, and XYZTPlus1PermTensor.
pfset PhaseSources.water.PredefinedFunction XPlusYPlusZ
The choices for this key correspond to sources as follows:
- X:
\({\rm source}\; = 0.0\)
- XPlusYPlusX:
\({\rm source}\; = 0.0\)
- X3Y2PlusSinXYPlus1:
- \({\rm source}\; = -(3x^2 y^2 + y\cos(xy))^2 - (2x^3 y + x\cos(xy))^2 - (x^3 y^2 + \sin(xy) + 1) (6x y^2 + 2x^3 -(x^2 +y^2) \sin(xy))\)This function type specifies that the source applied over the entire domain is as noted above. This corresponds to \(p=x^{3}y^{2}+\sin(xy)+1\) in the problem \(-\nabla\cdot (p\nabla p)=f\).
- X3Y4PlusX2PlusSinXYCosYPlus1:
- \({\rm source}\; = -(3x^22 y^4 + 2x + y\cos(xy)\cos(y))^2 - (4x^3 y^3 + x\cos(xy)\cos(y) - \sin(xy)\sin(y))^2 - (x^3 y^4 + x^2 + \sin(xy)\cos(y) + 1) (6xy^4 + 2 - (x^2 + y^2 + 1)\sin(xy)\cos(y) + 12x^3 y^2 - 2x\cos(xy)\sin(y))\)This function type specifies that the source applied over the entire domain is as noted above. This corresponds to \(p=x^{3}y^{4}+x^{2}+\sin (xy)\cos(y) +1\) in the problem \(-\nabla\cdot (p\nabla p)=f\).
- XYZTPlus1:
- \({\rm source}\; = xyz - t^2 (x^2 y^2 +x^2 z^2 +y^2 z^2)\)This function type specifies that the source applied over the entire domain is as noted above. This corresponds to \(p = xyzt + 1\) in the problem \(\frac{\partial p}{\partial t}-\nabla\cdot (p\nabla p)=f\).
- XYZTPlus1PermTensor:
- \({\rm source}\; = xyz - t^2 (x^2 y^2 3 + x^2 z^2 2 + y^2 z^2)\)This function type specifies that the source applied over the entire domain is as noted above. This corresponds to \(p = xyzt + 1\) in the problem \(\frac{\partial p}{\partial t}-\nabla\cdot (Kp\nabla p)=f\), where \(K = diag(1 \;\; 2 \;\; 3)\).
6.1.21. Capillary Pressures
Here we define capillary pressure. Note: this section needs to be defined only for multi-phase flow and should not be defined for single phase and Richards’ equation cases. The format for this section of input is:
string CapPressure.*phase_name*.Type “Constant” This key specifies the capillary pressure between phase \(0\) and the named phase, phase_name. The only choice available is Constant which indicates that a constant capillary pressure exists between the phases.
pfset CapPressure.water.Type Constant
list CapPressure.*phase_name*.GeomNames no default This key specifies the geometries that capillary pressures will be computed for in the named phase, phase_name. Regions listed later “overlay” regions listed earlier. Any geometries not listed will be assigned \(0.0\) capillary pressure by ParFlow.
pfset CapPressure.water.GeomNames "domain"
double Geom.*geometry_name*.CapPressure.*phase_name*.Value 0.0 This key specifies the value of the capillary pressure in the named geometry, geometry_name, for the named phase, phase_name.
pfset Geom.domain.CapPressure.water.Value 0.0
Important note: the code currently works only for capillary pressure equal zero.
6.1.22. Saturation
This section is only relevant to the Richards’ equation cases. All keys relating to this section will be ignored for other cases. The following keys are used to define the saturation-pressure curve.
string Phase.Saturation.Type no default This key specifies the type of saturation function that will be used on all specified geometries. Note that only one type of saturation may be used for the entire problem. However, parameters may be different for that type in different geometries. For instance, if the problem consists of three geometries, then VanGenuchten may be specified with three different sets of parameters for the three different goemetries. However, once VanGenuchten is specified, one geometry cannot later be specified to have Data as its saturation. The possible values for this key are Constant, VanGenuchten, Haverkamp, Data, Polynomial and PFBFile.
pfset Phase.Saturation.Type Constant
The various possible functions are defined as follows. The Constant specification means that the saturation will be constant on the specified geounit. The VanGenuchten specification means that the saturation will be given as a Van Genuchten function [VanGenuchten80] with the form,
where \(s_{sat}\) is the saturation at saturated conditions, \(s_{res}\) is the residual saturation, and \(\alpha\) and \(n\) are soil parameters with \(m = 1 - 1/n\), on each region. The Haverkamp specification means that the saturation will be given in the following form [Haverkamp-Vauclin81],
where \(A\) and \(\gamma\) are soil parameters, on each region. The Data specification is currently unsupported but will later mean that data points for the saturation curve will be given and ParFlow will set up the proper interpolation coefficients to get values between the given data points. The Polynomial specification defines a polynomial saturation function for each region of the form,
The PFBFile specification means that the saturation will be taken as a spatially varying but constant in pressure function given by data in a ParFlow binary (.pfb) file.
list Phase.Saturation.GeomNames no default This key specifies the geometries on which saturation will be given. The union of these geometries must cover the entire computational domain.
pfset Phase.Saturation.Geonames domain
double Geom.*geom_name*.Saturation.Value no default This key specifies the constant saturation value on the geom_name region.
pfset Geom.domain.Saturation.Value 0.5
integer Phase.Saturation.VanGenuchten.File 0 This key specifies whether soil parameters for the VanGenuchten function are specified in a pfb file or by region. The options are either 0 for specification by region, or 1 for specification in a file. Note that either all parameters are specified in files (each has their own input file) or none are specified by files. Parameters specified by files are \(\alpha\), N, SRes, and SSat.
pfset Phase.Saturation.VanGenuchten.File 1
string Geom.*geom_name*.Saturation.Alpha.Filename no default This key specifies a pfb filename containing the alpha parameters for the VanGenuchten function cell-by-cell. The ONLY option for geom_name is “domain”.
pfset Geom.domain.Saturation.Filename alphas.pfb
string Geom.*geom_name*.Saturation.N.Filename no default This key specifies a pfb filename containing the N parameters for the VanGenuchten function cell-by-cell. The ONLY option for geom_name is “domain”.
pfset Geom.domain.Saturation.N.Filename Ns.pfb
string Geom.*geom_name*.Saturation.SRes.Filename no default This key specifies a pfb filename containing the SRes parameters for the VanGenuchten function cell-by-cell. The ONLY option for geom_name is “domain”.
pfset Geom.domain.Saturation.SRes.Filename SRess.pfb
string Geom.*geom_name*.Saturation.SSat.Filename no default This key specifies a pfb filename containing the SSat parameters for the VanGenuchten function cell-by-cell. The ONLY option for geom_name is “domain”.
pfset Geom.domain.Saturation.SSat.Filename SSats.pfb
double Geom.*geom_name*.Saturation.Alpha no default This key specifies the \(\alpha\) parameter for the Van Genuchten function specified on geom_name.
pfset Geom.domain.Saturation.Alpha 0.005
double Geom.*geom_name*.Saturation.N no default This key specifies the \(N\) parameter for the Van Genuchten function specified on geom_name.
pfset Geom.domain.Saturation.N 2.0
Note that if both a Van Genuchten saturation and relative permeability are specified, then the soil parameters should be the same for each in order to have a consistent problem.
double Geom.*geom_name*.Saturation.SRes no default This key specifies the residual saturation on geom_name.
pfset Geom.domain.Saturation.SRes 0.0
double Geom.*geom_name*.Saturation.SSat no default This key specifies the saturation at saturated conditions on geom_name.
pfset Geom.domain.Saturation.SSat 1.0
double Geom.*geom_name*.Saturation.A no default This key specifies the \(A\) parameter for the Haverkamp saturation on geom_name.
pfset Geom.domain.Saturation.A 1.0
double Geom.*geom_name*.Saturation.Gamma no default This key specifies the the \(\gamma\) parameter for the Haverkamp saturation on geom_name.
pfset Geom.domain.Saturation.Gamma 1.0
integer Geom.*geom_name*.Saturation.Degree no default This key specifies the degree of the polynomial for the Polynomial saturation given on geom_name.
pfset Geom.domain.Saturation.Degree 1
double Geom.*geom_name*.Saturation.Coeff.*coeff_number* no default This key specifies the coeff_numberth coefficient of the Polynomial saturation given on geom_name.
pfset Geom.domain.Saturation.Coeff.0 0.5
pfset Geom.domain.Saturation.Coeff.1 1.0
string Geom.*geom_name*.Saturation.FileName no default This key specifies the name of the file containing saturation values for the domain. It is assumed that geom_name is “domain” for this key.
pfset Geom.domain.Saturation.FileName "domain_sats.pfb"
6.1.23. Internal Boundary Conditions
In this section, we define internal Dirichlet boundary conditions by setting the pressure at points in the domain. The format for this section of input is:
string InternalBC.Names no default This key specifies the names for the internal boundary conditions. At each named point, \({\rm x}\), \({\rm y}\) and \({\rm z}\) will specify the coordinate locations and \({\rm h}\) will specify the hydraulic head value of the condition. This real location is “snapped” to the nearest gridpoint in ParFlow.
NOTE: Currently, ParFlow assumes that internal boundary conditions and pressure wells are separated by at least one cell from any external boundary. The user should be careful of this when defining the input file and grid.
pfset InternalBC.Names "fixedvalue"
double InternalBC.*internal_bc_name*.X no default This key specifies the x-coordinate, \({\rm x}\), of the named, internal_bc_name, condition.
pfset InternalBC.fixedheadvalue.X 40.0
double InternalBC.*internal_bc_name*.Y no default This key specifies the y-coordinate, \({\rm y}\), of the named, internal_bc_name, condition.
pfset InternalBC.fixedheadvalue.Y 65.2
double InternalBC.*internal_bc_name*.Z no default This key specifies the z-coordinate, \({\rm z}\), of the named, internal_bc_name, condition.
pfset InternalBC.fixedheadvalue.Z 12.1
double InternalBC.*internal_bc_name*.Value no default This key specifies the value of the named, internal_bc_name, condition.
pfset InternalBC.fixedheadvalue.Value 100.0
6.1.24. Boundary Conditions: Pressure
Here we define the pressure boundary conditions. The Dirichlet conditions below are hydrostatic conditions, and it is assumed that at each phase interface the pressure is constant. It is also assumed here that all phases are distributed within the domain at all times such that the lighter phases are vertically higher than the heavier phases.
Boundary condition input is associated with domain patches (see §6.1.7 Domain). Note that different patches may have different types of boundary conditions on them.
list BCPressure.PatchNames no default This key specifies the names of patches on which pressure boundary conditions will be specified. Note that these must all be patches on the external boundary of the domain and these patches must “cover” that external boundary.
pfset BCPressure.PatchNames "left right front back top bottom"
string Patch.*patch_name*.BCPressure.Type no default This key
specifies the type of boundary condition data given for patch
patch_name. Possible values for this key are DirEquilRefPatch,
DirEquilPLinear, FluxConst, FluxVolumetric, PressureFile, FluxFile,
OverlandFow, OverlandFlowPFB and ExactSolution. The choice
DirEquilRefPatch specifies that the pressure on the specified patch
will be in hydrostatic equilibrium with a constant reference pressure
given on a reference patch. The choice DirEquilPLinear specifies
that the pressure on the specified patch will be in hydrostatic
equilibrium with pressure given along a piecewise line at elevation
\(z=0\). The choice FluxConst defines a constant normal flux
boundary condition through the domain patch. This flux must be specified
in units of \([L]/[T]\). For Richards’ equation, fluxes must be
specified as a mass flux and given as the above flux multiplied by the
density. Thus, this choice of input type for a Richards’ equation
problem has units of \(([L]/[T])([M]/[L]^3)\). The choice
FluxVolumetric defines a volumetric flux boundary condition through
the domain patch. The units should be consistent with all other user
input for the problem. For Richards’ equation fluxes must be specified
as a mass flux and given as the above flux multiplied by the density.
The choice PressureFile defines a hydraulic head boundary condition
that is read from a properly distributed .pfb file. Only the values
needed for the patch are used. The choice FluxFile defines a flux
boundary condition that is read form a properly distributed .pfb file
defined on a grid consistent with the pressure field grid. Only the
values needed for the patch are used. The choices OverlandFlow and
OverlandFlowPFB both turn on fully-coupled overland flow routing as
described in [KM06] and in §5.5 Overland Flow.
The key OverlandFlow corresponds to a Value key with a positive
or negative value, to indicate uniform fluxes (such as rainfall or
evapotranspiration) over the entire domain while the key
OverlandFlowPFB allows a .pfb
file to contain grid-based,
spatially-variable fluxes. The choice ExactSolution specifies that
an exact known solution is to be applied as a Dirichlet boundary
condition on the respective patch. Note that this does not change
according to any cycle. Instead, time dependence is handled by evaluating
at the time the boundary condition value is desired. The solution is specified
by using a predefined function (choices are described below). NOTE: These last
three types of boundary condition input is for Richards’ equation cases only!
pfset Patch.top.BCPressure.Type DirEquilRefPatch
string Patch.*patch_name*.BCPressure.Cycle no default This key specifies the time cycle to which boundary condition data for patch patch_name corresponds.
pfset Patch.top.BCPressure.Cycle Constant
string Patch.*patch_name*.BCPressure.RefGeom no default This key specifies the name of the solid on which the reference patch for the DirEquilRefPatch boundary condition data is given. Care should be taken to make sure the correct solid is specified in cases of layered domains.
pfset Patch.top.BCPressure.RefGeom domain
string Patch.*patch_name*.BCPressure.RefPatch no default This key specifies the reference patch on which the DirEquilRefPatch boundary condition data is given. This patch must be on the reference solid specified by the Patch.patch_name.BCPressure.RefGeom key.
pfset Patch.top.BCPressure.RefPatch bottom
double Patch.*patch_name*.BCPressure.*interval_name*.Value no default This key specifies the reference pressure value for the DirEquilRefPatch boundary condition or the constant flux value for the FluxConst boundary condition, or the constant volumetric flux for the FluxVolumetric boundary condition.
pfset Patch.top.BCPressure.alltime.Value -14.0
double Patch.*patch_name*.BCPressure.*interval_name*.*phase_name*.IntValue no default Note that the reference conditions for types DirEquilPLinear and DirEquilRefPatch boundary conditions are for phase 0 only. This key specifies the constant pressure value along the interface with phase phase_name for cases with two phases present.
pfset Patch.top.BCPressure.alltime.water.IntValue -13.0
double Patch.*patch_name*.BCPressure.*interval_name*.XLower no default This key specifies the lower \(x\) coordinate of a line in the xy-plane.
pfset Patch.top.BCPressure.alltime.XLower 0.0
double Patch.*patch_name*.BCPressure.*interval_name*.YLower no default This key specifies the lower \(y\) coordinate of a line in the xy-plane.
pfset Patch.top.BCPressure.alltime.YLower 0.0
double Patch.*patch_name*.BCPressure.*interval_name*.XUpper no default This key specifies the upper \(x\) coordinate of a line in the xy-plane.
pfset Patch.top.BCPressure.alltime.XUpper 1.0
double Patch.*patch_name*.BCPressure.*interval_name*.YUpper no default This key specifies the upper \(y\) coordinate of a line in the xy-plane.
pfset Patch.top.BCPressure.alltime.YUpper 1.0
integer Patch.*patch_name*.BCPressure.*interval_name*.NumPoints no default This key specifies the number of points on which pressure data is given along the line used in the type DirEquilPLinear boundary conditions.
pfset Patch.top.BCPressure.alltime.NumPoints 2
double Patch.*patch_name*.BCPressure.*interval_name*.*point_number*.Location no default This key specifies a number between 0 and 1 which represents the location of a point on the line on which data is given for type DirEquilPLinear boundary conditions. Here 0 corresponds to the lower end of the line, and 1 corresponds to the upper end.
pfset Patch.top.BCPressure.alltime.0.Location 0.0
double Patch.*patch_name*.BCPressure.*interval_name*.*point_number*.Value no default This key specifies the pressure value for phase 0 at point number point_number and \(z=0\) for type DirEquilPLinear boundary conditions. All pressure values on the patch are determined by first projecting the boundary condition coordinate onto the line, then linearly interpolating between the neighboring point pressure values on the line.
pfset Patch.top.BCPressure.alltime.0.Value 14.0
string Patch.*patch_name*.BCPressure.*interval_name*.FileName
no default This key specifies the name of a properly distributed .pfb
file
that contains boundary data to be read for types PressureFile and FluxFile.
For flux data, the data must be defined over a grid consistent with the
pressure field. In both cases, only the values needed for the patch will
be used. The rest of the data is ignored.
pfset Patch.top.BCPressure.alltime.FileName ocwd_bc.pfb
string Patch.*patch_name*.BCPressure.*interval_name*.PredefinedFunction no default This key specifies the predefined function that will be used to specify Dirichlet boundary conditions on patch patch_name. Note that this does not change according to any cycle. Instead, time dependence is handled by evaluating at the time the boundary condition value is desired. Choices for this key include X, XPlusYPlusZ, X3Y2PlusSinXYPlus1, X3Y4PlusX2PlusSinXYCosYPlus1, XYZTPlus1 and XYZTPlus1PermTensor.
pfset Patch.top.BCPressure.alltime.PredefinedFunction XPlusYPlusZ
The choices for this key correspond to pressures as follows.
- X:
\(p = x\)
- XPlusYPlusZ:
\(p = x + y + z\)
- X3Y2PlusSinXYPlus1:
\(p = x^3 y^2 + \sin(xy) + 1\)
- X3Y4PlusX2PlusSinXYCosYPlus1:
\(p = x^3 y^4 + x^2 + \sin(xy)\cos y + 1\)
- XYZTPlus1:
\(p = xyzt + 1\)
- XYZTPlus1PermTensor:
\(p = xyzt + 1\)
Example Script:
#---------------------------------------------------------
# Initial conditions: water pressure [m]
#---------------------------------------------------------
# Using a patch is great when you are not using a box domain
# If using a box domain HydroStaticDepth is fine
# If your RefPatch is z-lower (bottom of domain), the pressure is positive.
# If your RefPatch is z-upper (top of domain), the pressure is negative.
### Set water table to be at the bottom of the domain, the top layer is initially dry
pfset ICPressure.Type HydroStaticPatch
pfset ICPressure.GeomNames domain
pfset Geom.domain.ICPressure.Value 2.2
pfset Geom.domain.ICPressure.RefGeom domain
pfset Geom.domain.ICPressure.RefPatch z-lower
### Using a .pfb to initialize
pfset ICPressure.Type PFBFile
pfset ICPressure.GeomNames "domain"
pfset Geom.domain.ICPressure.FileName press.00090.pfb
pfset Geom.domain.ICPressure.RefGeom domain
pfset Geom.domain.ICPressure.RefPatch z-upper
6.1.25. Boundary Conditions: Saturation
Note: this section needs to be defined only for multi-phase flow and should not be defined for the single phase and Richards’ equation cases.
Here we define the boundary conditions for the saturations. Boundary condition input is associated with domain patches (see §6.1.7 Domain). Note that different patches may have different types of boundary conditions on them.
list BCSaturation.PatchNames no default This key specifies the names of patches on which saturation boundary conditions will be specified. Note that these must all be patches on the external boundary of the domain and these patches must “cover” that external boundary.
pfset BCSaturation.PatchNames "left right front back top bottom"
string Patch.*patch_name*.BCSaturation.*phase_name*.Type no default This key specifies the type of boundary condition data given for the given phase, phase_name, on the given patch patch_name. Possible values for this key are DirConstant, ConstantWTHeight and PLinearWTHeight. The choice DirConstant specifies that the saturation is constant on the whole patch. The choice ConstantWTHeight specifies a constant height of the water-table on the whole patch. The choice PLinearWTHeight specifies that the height of the water-table on the patch will be given by a piecewise linear function.
Note: the types ConstantWTHeight and PLinearWTHeight assume we are running a 2-phase problem where phase 0 is the water phase.
pfset Patch.left.BCSaturation.water.Type ConstantWTHeight
double Patch.*patch_name*.BCSaturation.*phase_name*.Value no default This key specifies either the constant saturation value if DirConstant is selected or the constant water-table height if ConstantWTHeight is selected.
pfset Patch.top.BCSaturation.air.Value 1.0
double Patch.*patch_name*.BCSaturation.*phase_name*.XLower no default This key specifies the lower \(x\) coordinate of a line in the xy-plane if type PLinearWTHeight boundary conditions are specified.
pfset Patch.left.BCSaturation.water.XLower -10.0
double Patch.*patch_name*.BCSaturation.*phase_name*.YLower no default This key specifies the lower \(y\) coordinate of a line in the xy-plane if type PLinearWTHeight boundary conditions are specified.
pfset Patch.left.BCSaturation.water.YLower 5.0
double Patch.*patch_name*.BCSaturation.*phase_name*.XUpper no default This key specifies the upper \(x\) coordinate of a line in the xy-plane if type PLinearWTHeight boundary conditions are specified.
pfset Patch.left.BCSaturation.water.XUpper 125.0
double Patch.*patch_name*.BCSaturation.*phase_name*.YUpper no default This key specifies the upper \(y\) coordinate of a line in the xy-plane if type PLinearWTHeight boundary conditions are specified.
pfset Patch.left.BCSaturation.water.YUpper 82.0
integer Patch.*patch_name*.BCPressure.*phase_name*.NumPoints no default This key specifies the number of points on which saturation data is given along the line used for type DirEquilPLinear boundary conditions.
pfset Patch.left.BCPressure.water.NumPoints 2
double Patch.*patch_name*.BCPressure.*phase_name*.*point_number*.Location no default This key specifies a number between 0 and 1 which represents the location of a point on the line for which data is given in type DirEquilPLinear boundary conditions. The line is parameterized so that 0 corresponds to the lower end of the line, and 1 corresponds to the upper end.
pfset Patch.left.BCPressure.water.0.Location 0.333
double Patch.*patch_name*.BCPressure.*phase_name*.*point_number*.Value no default This key specifies the water-table height for the given point if type DirEquilPLinear boundary conditions are selected. All saturation values on the patch are determined by first projecting the water-table height value onto the line, then linearly interpolating between the neighboring water-table height values onto the line.
pfset Patch.left.BCPressure.water.0.Value 4.5
6.1.26. Initial Conditions: Phase Saturations
Note: this section needs to be defined only for multi-phase flow and should not be defined for single phase and Richards’ equation cases.
Here we define initial phase saturation conditions. The format for this section of input is:
string ICSaturation.*phase_name*.Type no default This key specifies the type of initial condition that will be applied to different geometries for given phase, phase_name. The only key currently available is Constant. The choice Constant will apply constants values within geometries for the phase.
ICSaturation.water.Type Constant
string ICSaturation.*phase_name*.GeomNames no default This key specifies the geometries on which an initial condition will be given if the type is set to Constant.
Note that geometries listed later “overlay” geometries listed earlier.
ICSaturation.water.GeomNames "domain"
double Geom.*geom_input_name*.ICSaturation.*phase_name*.Value no default This key specifies the initial condition value assigned to all points in the named geometry, geom_input_name, if the type was set to Constant.
Geom.domain.ICSaturation.water.Value 1.0
6.1.27. Initial Conditions: Pressure
The keys in this section are used to specify pressure initial conditions for Richards’ equation cases only. These keys will be ignored if any other case is run.
string ICPressure.Type no default This key specifies the type of initial condition given. The choices for this key are Constant, HydroStaticDepth, HydroStaticPatch and PFBFile. The choice Constant specifies that the initial pressure will be constant over the regions given. The choice HydroStaticDepth specifies that the initial pressure within a region will be in hydrostatic equilibrium with a given pressure specified at a given depth. The choice HydroStaticPatch specifies that the initial pressure within a region will be in hydrostatic equilibrium with a given pressure on a specified patch. Note that all regions must have the same type of initial data - different regions cannot have different types of initial data. However, the parameters for the type may be different. The PFBFile specification means that the initial pressure will be taken as a spatially varying function given by data in a ParFlow binary (.pfb) file.
pfset ICPressure.Type Constant
list ICPressure.GeomNames no default This key specifies the geometry names on which the initial pressure data will be given. These geometries must comprise the entire domain. Note that conditions for regions that overlap other regions will have unpredictable results. The regions given must be disjoint.
pfset ICPressure.GeomNames "toplayer middlelayer bottomlayer"
double Geom.*geom_name*.ICPressure.Value no default This key specifies the initial pressure value for type Constant initial pressures and the reference pressure value for types HydroStaticDepth and HydroStaticPatch.
pfset Geom.toplayer.ICPressure.Value -734.0
double Geom.*geom_name*.ICPressure.RefElevation no default This key specifies the reference elevation on which the reference pressure is given for type HydroStaticDepth initial pressures.
pfset Geom.toplayer.ICPressure.RefElevation 0.0
double Geom.*geom_name*.ICPressure.RefGeom no default This key specifies the geometry on which the reference patch resides for type HydroStaticPatch initial pressures.
pfset Geom.toplayer.ICPressure.RefGeom bottomlayer
double Geom.*geom_name*.ICPressure.RefPatch no default This key specifies the patch on which the reference pressure is given for type HydorStaticPatch initial pressures.
pfset Geom.toplayer.ICPressure.RefPatch bottom
string Geom.*geom_name*.ICPressure.FileName no default This key specifies the name of the file containing pressure values for the domain. It is assumed that geom_name is “domain” for this key.
pfset Geom.domain.ICPressure.FileName "ic_pressure.pfb"
6.1.28. Initial Conditions: Phase Concentrations
Here we define initial concentration conditions for contaminants. The format for this section of input is:
string PhaseConcen.*phase_name*.*contaminant_name*.Type no default This key specifies the type of initial condition that will be applied to different geometries for given phase, phase_name, and the given contaminant, contaminant_name. The choices for this key are Constant or PFBFile. The choice Constant will apply constants values to different geometries. The choice PFBFile will read values from a “ParFlow Binary” file (see §6.3 ParFlow Binary Files (.pfb)).
PhaseConcen.water.tce.Type Constant
string PhaseConcen.*phase_name*.GeomNames no default This key specifies the geometries on which an initial condition will be given, if the type was set to Constant.
Note that geometries listed later “overlay” geometries listed earlier.
PhaseConcen.water.GeomNames "ic_concen_region"
double PhaseConcen.*phase_name*.*contaminant_name*.*geom_input_name*.Value no default This key specifies the initial condition value assigned to all points in the named geometry, geom_input_name, if the type was set to Constant.
PhaseConcen.water.tce.ic_concen_region.Value 0.001
string PhaseConcen.*phase_name*.*contaminant_name*.FileName no default This key specifies the name of the “ParFlow Binary” file which contains the initial condition values if the type was set to PFBFile.
PhaseConcen.water.tce.FileName "initial_concen_tce.pfb"
6.1.29. Known Exact Solution
For Richards equation cases only we allow specification of an exact solution to be used for testing the code. Only types that have been coded and predefined are allowed. Note that if this is speccified as something other than no known solution, corresponding boundary conditions and phase sources should also be specified.
string KnownSolution no default This specifies the predefined function that will be used as the known solution. Possible choices for this key are NoKnownSolution, Constant, X, XPlusYPlusZ, X3Y2PlusSinXYPlus1, X3Y4PlusX2PlusSinXYCosYPlus1, XYZTPlus1 and XYZTPlus1PermTensor.
pfset KnownSolution XPlusYPlusZ
Choices for this key correspond to solutions as follows.
- NoKnownSolution:
No solution is known for this problem.
- Constant:
\(p = {\rm constant}\)
- X:
\(p = x\)
- XPlusYPlusZ:
\(p = x + y + z\)
- X3Y2PlusSinXYPlus1:
\(p = x^3 y^2 + sin(xy) + 1\)
- X3Y4PlusX2PlusSinXYCosYPlus1:
\(p = x^3 y^4 + x^2 + \sin(xy)\cos y + 1\)
- XYZTPlus1:
\(p = xyzt + 1\)
- XYZTPlus1:
\(p = xyzt + 1\)
double KnownSolution.Value no default This key specifies the constant value of the known solution for type Constant known solutions.
pfset KnownSolution.Value 1.0
Only for known solution test cases will information on the \(L^2\)-norm of the pressure error be printed.
6.1.30. Wells
Here we define wells for the model. The format for this section of input is:
string Wells.Names no default This key specifies the names of the wells for which input data will be given.
Wells.Names "test_well inj_well ext_well"
string Wells.*well_name*.InputType no default This key specifies the type of well to be defined for the given well, well_name. This key can be either Vertical or Recirc. The value Vertical indicates that this is a single segmented well whose action will be specified by the user. The value Recirc indicates that this is a dual segmented, recirculating, well with one segment being an extraction well and another being an injection well. The extraction well filters out a specified fraction of each contaminant and recirculates the remainder to the injection well where the diluted fluid is injected back in. The phase saturations at the extraction well are passed without modification to the injection well.
Note with the recirculating well, several input options are not needed as the extraction well will provide these values to the injection well.
Wells.test_well.InputType Vertical
string Wells.*well_name*.Action no default This key specifies the pumping action of the well. This key can be either Injection or Extraction. A value of Injection indicates that this is an injection well. A value of Extraction indicates that this is an extraction well.
Wells.test_well.Action Injection
double Wells.*well_name*.Type no default This key specfies the mechanism by which the well works (how ParFlow works with the well data) if the input type key is set to Vectical. This key can be either Pressure or Flux. A value of Pressure indicates that the data provided for the well is in terms of hydrostatic pressure and ParFlow will ensure that the computed pressure field satisfies this condition in the computational cells which define the well. A value of Flux indicates that the data provided is in terms of volumetric flux rates and ParFlow will ensure that the flux field satisfies this condition in the computational cells which define the well.
Wells.test_well.Type Flux
string Wells.*well_name*.ExtractionType no default This key specfies the mechanism by which the extraction well works (how ParFlow works with the well data) if the input type key is set to Recirc. This key can be either Pressure or Flux. A value of Pressure indicates that the data provided for the well is in terms of hydrostatic pressure and ParFlow will ensure that the computed pressure field satisfies this condition in the computational cells which define the well. A value of Flux indicates that the data provided is in terms of volumetric flux rates and ParFlow will ensure that the flux field satisfies this condition in the computational cells which define the well.
Wells.ext_well.ExtractionType Pressure
string Wells.*well_name*.InjectionType no default This key specfies the mechanism by which the injection well works (how ParFlow works with the well data) if the input type key is set to Recirc. This key can be either Pressure or Flux. A value of Pressure indicates that the data provided for the well is in terms of hydrostatic pressure and ParFlow will ensure that the computed pressure field satisfies this condition in the computational cells which define the well. A value of Flux indicates that the data provided is in terms of volumetric flux rates and ParFlow will ensure that the flux field satisfies this condition in the computational cells which define the well.
Wells.inj_well.InjectionType Flux
double Wells.*well_name*.X no default This key specifies the x location of the vectical well if the input type is set to Vectical or of both the extraction and injection wells if the input type is set to Recirc.
Wells.test_well.X 20.0
double Wells.*well_name*.Y no default This key specifies the y location of the vectical well if the input type is set to Vectical or of both the extraction and injection wells if the input type is set to Recirc.
Wells.test_well.Y 36.5
double Wells.*well_name*.ZUpper no default This key specifies the z location of the upper extent of a vectical well if the input type is set to Vectical.
Wells.test_well.ZUpper 8.0
double Wells.*well_name*.ExtractionZUpper no default This key specifies the z location of the upper extent of a extraction well if the input type is set to Recirc.
Wells.ext_well.ExtractionZUpper 3.0
double Wells.*well_name*.InjectionZUpper no default This key specifies the z location of the upper extent of a injection well if the input type is set to Recirc.
Wells.inj_well.InjectionZUpper 6.0
double Wells.*well_name*.ZLower no default This key specifies the z location of the lower extent of a vectical well if the input type is set to Vectical.
Wells.test_well.ZLower 2.0
double Wells.*well_name*.ExtractionZLower no default This key specifies the z location of the lower extent of a extraction well if the input type is set to Recirc.
Wells.ext_well.ExtractionZLower 1.0
double Wells.*well_name*.InjectionZLower no default This key specifies the z location of the lower extent of a injection well if the input type is set to Recirc.
Wells.inj_well.InjectionZLower 4.0
string Wells.*well_name*.Method no default This key specifies a method by which pressure or flux for a vertical well will be weighted before assignment to computational cells. This key can only be Standard if the type key is set to Pressure; or this key can be either Standard, Weighted or Patterned if the type key is set to Flux. A value of Standard indicates that the pressure or flux data will be used as is. A value of Weighted indicates that the flux data is to be weighted by the cells permeability divided by the sum of all cell permeabilities which define the well. The value of Patterned is not implemented.
Wells.test_well.Method Weighted
string Wells.*well_name*.ExtractionMethod no default This key specifies a method by which pressure or flux for an extraction well will be weighted before assignment to computational cells. This key can only be Standard if the type key is set to Pressure; or this key can be either Standard, Weighted or Patterned if the type key is set to Flux. A value of Standard indicates that the pressure or flux data will be used as is. A value of Weighted indicates that the flux data is to be weighted by the cells permeability divided by the sum of all cell permeabilities which define the well. The value of Patterned is not implemented.
Wells.ext_well.ExtractionMethod Standard
string Wells.*well_name*.InjectionMethod no default This key specifies a method by which pressure or flux for an injection well will be weighted before assignment to computational cells. This key can only be Standard if the type key is set to Pressure; or this key can be either Standard, Weighted or Patterned if the type key is set to Flux. A value of Standard indicates that the pressure or flux data will be used as is. A value of Weighted indicates that the flux data is to be weighted by the cells permeability divided by the sum of all cell permeabilities which define the well. The value of Patterned is not implemented.
Wells.inj_well.InjectionMethod Standard
string Wells.*well_name*.Cycle no default This key specifies the time cycles to which data for the well well_name corresponds.
Wells.test_well.Cycle "all_time"
double Wells.*well_name*.*interval_name*.Pressure.Value no default This key specifies the hydrostatic pressure value for a vectical well if the type key is set to Pressure.
Note This value gives the pressure of the primary phase (water) at \(z=0\). The other phase pressures (if any) are computed from the physical relationships that exist between the phases.
Wells.test_well.all_time.Pressure.Value 6.0
double Wells.*well_name*.*interval_name*.Extraction.Pressure.Value no default This key specifies the hydrostatic pressure value for an extraction well if the extraction type key is set to Pressure.
Note This value gives the pressure of the primary phase (water) at \(z=0\). The other phase pressures (if any) are computed from the physical relationships that exist between the phases.
Wells.ext_well.all_time.Extraction.Pressure.Value 4.5
double Wells.*well_name*.*interval_name*.Injection.Pressure.Value no default This key specifies the hydrostatic pressure value for an injection well if the injection type key is set to Pressure.
Note This value gives the pressure of the primary phase (water) at \(z=0\). The other phase pressures (if any) are computed from the physical relationships that exist between the phases.
Wells.inj_well.all_time.Injection.Pressure.Value 10.2
double Wells.*well_name*.*interval_name*.Flux.*phase_name*.Value no default This key specifies the volumetric flux for a vectical well if the type key is set to Flux.
Note only a positive number should be entered, ParFlow assignes the correct sign based on the chosen action for the well.
Wells.test_well.all_time.Flux.water.Value 250.0
double Wells.*well_name*.*interval_name*.Extraction.Flux.*phase_name*.Value no default This key specifies the volumetric flux for an extraction well if the extraction type key is set to Flux.
Note only a positive number should be entered, ParFlow assignes the correct sign based on the chosen action for the well.
Wells.ext_well.all_time.Extraction.Flux.water.Value 125.0
double Wells.*well_name*.*interval_name*.Injection.Flux.*phase_name*.Value no default This key specifies the volumetric flux for an injection well if the injection type key is set to Flux.
Note only a positive number should be entered, ParFlow assignes the correct sign based on the chosen action for the well.
Wells.inj_well.all_time.Injection.Flux.water.Value 80.0
double Wells.*well_name*.*interval_name*.Saturation.*phase_name*.Value no default This key specifies the saturation value of a vertical well.
Wells.test_well.all_time.Saturation.water.Value 1.0
double Wells.*well_name*.*interval_name*.Concentration.*phase_name*.*contaminant_name*.Value no default This key specifies the contaminant value of a vertical well.
Wells.test_well.all_time.Concentration.water.tce.Value 0.0005
double Wells.*well_name*.*interval_name*.Injection.Concentration.*phase_name*.*contaminant_name*.Fraction no default This key specifies the fraction of the extracted contaminant which gets resupplied to the injection well.
Wells.inj_well.all_time.Injection.Concentration.water.tce.Fraction 0.01
Multiple wells assigned to one grid location can occur in several instances. The current actions taken by the code are as follows:
If multiple pressure wells are assigned to one grid cell, the code retains only the last set of overlapping well values entered.
If multiple flux wells are assigned to one grid cell, the code sums the contributions of all overlapping wells to get one effective well flux.
If multiple pressure and flux wells are assigned to one grid cell, the code retains the last set of overlapping hydrostatic pressure values entered and sums all the overlapping flux well values to get an effective pressure/flux well value.
6.1.31. Code Parameters
In addition to input keys related to the physics capabilities and modeling specifics there are some key values used by various algorithms and general control flags for ParFlow. These are described next :
string Solver.Linear PCG This key specifies the linear solver used for solver IMPES. Choices for this key are MGSemi, PPCG, PCG and CGHS. The choice MGSemi is an algebraic mulitgrid linear solver (not a preconditioned conjugate gradient) which may be less robust than PCG as described in [Ashby-Falgout90]. The choice PPCG is a preconditioned conjugate gradient solver. The choice PCG is a conjugate gradient solver with a multigrid preconditioner. The choice CGHS is a conjugate gradient solver.
pfset Solver.Linear MGSemi
integer Solver.SadvectOrder 2 This key controls the order of the explicit method used in advancing the saturations. This value can be either 1 for a standard upwind first order or 2 for a second order Godunov method.
pfset Solver.SadvectOrder 1
integer Solver.AdvectOrder 2 This key controls the order of the explicit method used in advancing the concentrations. This value can be either 1 for a standard upwind first order or 2 for a second order Godunov method.
pfset Solver.AdvectOrder 2
double Solver.CFL 0.7 This key gives the value of the weight put on the computed CFL limit before computing a global timestep value. Values greater than 1 are not suggested and in fact because this is an approximation, values slightly less than 1 can also produce instabilities.
pfset Solver.CFL 0.7
integer Solver.MaxIter 1000000 This key gives the maximum number of iterations that will be allowed for time-stepping. This is to prevent a run-away simulation.
pfset Solver.MaxIter 100
double Solver.RelTol 1.0 This value gives the relative tolerance for the linear solve algorithm.
pfset Solver.RelTol 1.0
double Solver.AbsTol 1E-9 This value gives the absolute tolerance for the linear solve algorithm.
pfset Solver.AbsTol 1E-8
double Solver.Drop 1E-8 This key gives a clipping value for data written to PFSB files. Data values greater than the negative of this value and less than the value itself are treated as zero and not written to PFSB files.
pfset Solver.Drop 1E-6
string Solver.PrintSubsurf True This key is used to turn on printing of the subsurface data, Permeability and Porosity. The data is printed after it is generated and before the main time stepping loop - only once during the run. The data is written as a PFB file.
pfset Solver.PrintSubsurf False
string Solver.PrintPressure True This key is used to turn on printing of the pressure data. The printing of the data is controlled by values in the timing information section. The data is written as a PFB file.
pfset Solver.PrintPressure False
string Solver.PrintVelocities False This key is used to turn on printing of the x, y and z velocity data. The printing of the data is controlled by values in the timing information section. The data is written as a PFB file.
pfset Solver.PrintVelocities True
string Solver.PrintSaturation True This key is used to turn on printing of the saturation data. The printing of the data is controlled by values in the timing information section. The data is written as a PFB file.
pfset Solver.PrintSaturation False
string Solver.PrintConcentration True This key is used to turn on printing of the concentration data. The printing of the data is controlled by values in the timing information section. The data is written as a PFSB file.
pfset Solver.PrintConcentration False
string Solver.PrintWells True This key is used to turn on collection and printing of the well data. The data is collected at intervals given by values in the timing information section. Printing occurs at the end of the run when all collected data is written.
pfset Solver.PrintWells False
string Solver.PrintLSMSink False This key is used to turn on
printing of the flux array passed from CLM
to ParFlow.
Printing occurs at each DumpInterval time.
pfset Solver.PrintLSMSink True
string Solver.WriteSiloSubsurfData False This key is used to specify printing of the subsurface data, Permeability and Porosity in silo binary file format. The data is printed after it is generated and before the main time stepping loop - only once during the run. This data may be read in by VisIT and other visualization packages.
pfset Solver.WriteSiloSubsurfData True
string Solver.WriteSiloPressure False This key is used to specify printing of the saturation data in silo binary format. The printing of the data is controlled by values in the timing information section. This data may be read in by VisIT and other visualization packages.
pfset Solver.WriteSiloPressure True
string Solver.WriteSiloSaturation False This key is used to specify printing of the saturation data using silo binary format. The printing of the data is controlled by values in the timing information section.
pfset Solver.WriteSiloSaturation True
string Solver.WriteSiloConcentration False This key is used to specify printing of the concentration data in silo binary format. The printing of the data is controlled by values in the timing information section.
pfset Solver.WriteSiloConcentration True
string Solver.WriteSiloVelocities False This key is used to specify printing of the x, y and z velocity data in silo binary format. The printing of the data is controlled by values in the timing information section.
pfset Solver.WriteSiloVelocities True
string Solver.WriteSiloSlopes False This key is used to specify printing of the x and y slope data using silo binary format. The printing of the data is controlled by values in the timing information section.
pfset Solver.WriteSiloSlopes True
string Solver.WriteSiloMannings False This key is used to specify printing of the Manning’s roughness data in silo binary format. The printing of the data is controlled by values in the timing information section.
pfset Solver.WriteSiloMannings True
string Solver.WriteSiloSpecificStorage False This key is used to specify printing of the specific storage data in silo binary format. The printing of the data is controlled by values in the timing information section.
pfset Solver.WriteSiloSpecificStorage True
string Solver.WriteSiloMask False This key is used to specify printing of the mask data using silo binary format. The mask contains values equal to one for active cells and zero for inactive cells. The printing of the data is controlled by values in the timing information section.
pfset Solver.WriteSiloMask True
string Solver.WriteSiloEvapTrans False This key is used to specify
printing of the evaporation and rainfall flux data using silo binary
format. This data comes from either clm
or from external calls to
ParFlow such as WRF. This data is in units of [L^3 T^{-1}]. The printing
of the data is controlled by values in the timing information section.
pfset Solver.WriteSiloEvapTrans True
string Solver.WriteSiloEvapTransSum False This key is used to
specify printing of the evaporation and rainfall flux data using silo
binary format as a running, cumulative amount. This data comes from
either clm
or from external calls to ParFlow such as WRF. This
data is in units of [L^3]. The printing of the data is controlled by
values in the timing information section.
pfset Solver.WriteSiloEvapTransSum True
string Solver.WriteSiloOverlandSum False This key is used to specify calculation and printing of the total overland outflow from the domain using silo binary format as a running cumulative amount. This is integrated along all domain boundaries and is calculated any location that slopes at the edge of the domain point outward. This data is in units of \([L^3]\). The printing of the data is controlled by values in the timing information section.
pfset Solver.WriteSiloOverlandSum True
string Solver.TerrainFollowingGrid False This key specifies that a terrain-following coordinate transform is used for solver Richards. This key sets x and y subsurface slopes to be the same as the Topographic slopes (a value of False sets these subsurface slopes to zero). These slopes are used in the Darcy fluxes to add a density, gravity -dependent term. This key will not change the output files (that is the output is still orthogonal) or the geometries (they will still follow the computational grid)– these two things are both to do items. This key only changes solver Richards, not solver Impes.
pfset Solver.TerrainFollowingGrid True
6.1.32. SILO Options
The following keys are used to control how SILO writes data. SILO allows writing to PDB and HDF5 file formats. SILO also allows data compression to be used, which can save signicant amounts of disk space for some problems.
string SILO.Filetype PDB This key is used to specify the SILO filetype. Allowed values are PDB and HDF5. Note that you must have configured SILO with HDF5 in order to use that option.
pfset SILO.Filetype PDB
string SILO.CompressionOptions This key is used to specify the SILO compression options. See the SILO manual for the DB_SetCompression command for information on available options. NOTE: the options avaialable are highly dependent on the configure options when building SILO.
pfset SILO.CompressionOptions ``METHOD=GZIP''
6.1.33. Richards’ Equation Solver Parameters
The following keys are used to specify various parameters used by the linear and nonlinear solvers in the Richards’ equation implementation. For information about these solvers, see [Woodward98] and [Ashby-Falgout90].
double Solver.Nonlinear.ResidualTol 1e-7 This key specifies the tolerance that measures how much the relative reduction in the nonlinear residual should be before nonlinear iterations stop. The magnitude of the residual is measured with the \(l^1\) (max) norm.
pfset Solver.Nonlinear.ResidualTol 1e-4
double Solver.Nonlinear.StepTol 1e-7 This key specifies the tolerance that measures how small the difference between two consecutive nonlinear steps can be before nonlinear iterations stop.
pfset Solver.Nonlinear.StepTol 1e-4
integer Solver.Nonlinear.MaxIter 15 This key specifies the maximum number of nonlinear iterations allowed before iterations stop with a convergence failure.
pfset Solver.Nonlinear.MaxIter 50
integer Solver.Linear.KrylovDimension 10 This key specifies the maximum number of vectors to be used in setting up the Krylov subspace in the GMRES iterative solver. These vectors are of problem size and it should be noted that large increases in this parameter can limit problem sizes. However, increasing this parameter can sometimes help nonlinear solver convergence.
pfset Solver.Linear.KrylovDimension 15
integer Solver.Linear.MaxRestarts 0 This key specifies the number of restarts allowed to the GMRES solver. Restarts start the development of the Krylov subspace over using the current iterate as the initial iterate for the next pass.
pfset Solver.Linear.MaxRestarts 2
integer Solver.MaxConvergencFailures 3 This key gives the maximum number of convergence failures allowed. Each convergence failure cuts the timestep in half and the solver tries to advance the solution with the reduced timestep.
The default value is 3.
Note that setting this value to a value greater than 9 may result in errors in how time cycles are calculated. Time is discretized in terms of the base time unit and if the solver begins to take very small timesteps \(smaller than base time unit \/ 1000\) the values based on time cycles will be change at slightly incorrect times. If the problem is failing converge so poorly that a large number of restarts are required, consider setting the timestep to a smaller value.
pfset Solver.MaxConvergenceFailures 4
string Solver.Nonlinear.PrintFlag HighVerbosity This key specifies
the amount of informational data that is printed to the *.out.kinsol.log
file. Choices for this key are NoVerbosity, LowVerbosity, NormalVerbosity
and HighVerbosity. The choice NoVerbosity prints no statistics about the
nonlinear convergence process. The choice LowVerbosity outputs the nonlinear
iteration count, the scaled norm of the nonlinear function, and the number of
function calls. The choice NormalVerbosity prints the same as for LowVerbosity
and also the global strategy statistics. The choice HighVerbosity prints the
same as for NormalVerbosity with the addition of further Krylov iteration
statistics.
pfset Solver.Nonlinear.PrintFlag NormalVerbosity
string Solver.Nonlinear.EtaChoice Walker2 This key specifies how the linear system tolerance will be selected. The linear system is solved until a relative residual reduction of \(\eta\) is achieved. Linear residuall norms are measured in the \(l^2\) norm. Choices for this key include EtaConstant, Walker1 and Walker2. If the choice EtaConstant is specified, then \(\eta\) will be taken as constant. The choices Walker1 and Walker2 specify choices for \(\eta\) developed by Eisenstat and Walker [EW96]. The choice Walker1 specifies that \(\eta\) will be given by \(| \|F(u^k)\| - \|F(u^{k-1}) + J(u^{k-1})*p \| | / \|F(u^{k-1})\|\). The choice Walker2 specifies that \(\eta\) will be given by \(\gamma \|F(u^k)\| / \|F(u^{k-1})\|^{\alpha}\). For both of the last two choices, \(\eta\) is never allowed to be less than 1e-4.
pfset Solver.Nonlinear.EtaChoice EtaConstant
double Solver.Nonlinear.EtaValue 1e-4 This key specifies the constant value of \(\eta\) for the EtaChoice key EtaConstant.
pfset Solver.Nonlinear.EtaValue 1e-7
double Solver.Nonlinear.EtaAlpha 2.0 This key specifies the value of \(\alpha\) for the case of EtaChoice being Walker2.
pfset Solver.Nonlinear.EtaAlpha 1.0
double Solver.Nonlinear.EtaGamma 0.9 This key specifies the value of \(\gamma\) for the case of EtaChoice being Walker2.
pfset Solver.Nonlinear.EtaGamma 0.7
string Solver.Nonlinear.UseJacobian False This key specifies whether the Jacobian will be used in matrix-vector products or whether a matrix-free version of the code will run. Choices for this key are False and True. Using the Jacobian will most likely decrease the number of nonlinear iterations but require more memory to run.
pfset Solver.Nonlinear.UseJacobian True
double Solver.Nonlinear.DerivativeEpsilon 1e-7 This key specifies the value of \(\epsilon\) used in approximating the action of the Jacobian on a vector with approximate directional derivatives of the nonlinear function. This parameter is only used when the UseJacobian key is False.
pfset Solver.Nonlinear.DerivativeEpsilon 1e-8
string Solver.Nonlinear.Globalization LineSearch This key specifies the type of global strategy to use. Possible choices for this key are InexactNewton and LineSearch. The choice InexactNewton specifies no global strategy, and the choice LineSearch specifies that a line search strategy should be used where the nonlinear step can be lengthened or decreased to satisfy certain criteria.
pfset Solver.Nonlinear.Globalization LineSearch
string Solver.Linear.Preconditioner MGSemi This key specifies which preconditioner to use. Currently, the three choices are NoPC, MGSemi, PFMG, PFMGOctree and SMG. The choice NoPC specifies that no preconditioner should be used. The choice MGSemi specifies a semi-coarsening multigrid algorithm which uses a point relaxation method. The choice SMG specifies a semi-coarsening multigrid algorithm which uses plane relaxations. This method is more robust than MGSemi, but generally requires more memory and compute time. The choice PFMGOctree can be more efficient for problems with large numbers of inactive cells.
pfset Solver.Linear.Preconditioner MGSemi
string Solver.Linear.Preconditioner.SymmetricMat Symmetric This key specifies whether the preconditioning matrix is symmetric. Choices fo rthis key are Symmetric and Nonsymmetric. The choice Symmetric specifies that the symmetric part of the Jacobian will be used as the preconditioning matrix. The choice Nonsymmetric specifies that the full Jacobian will be used as the preconditioning matrix. NOTE: ONLY Symmetric CAN BE USED IF MGSemi IS THE SPECIFIED PRECONDITIONER!
pfset Solver.Linear.Preconditioner.SymmetricMat Symmetric
integer Solver.Linear.Preconditioner.*precond_method*.MaxIter 1 This key specifies the maximum number of iterations to take in solving the preconditioner system with precond_method solver.
pfset Solver.Linear.Preconditioner.SMG.MaxIter 2
integer Solver.Linear.Preconditioner.SMG.NumPreRelax 1 This key specifies the number of relaxations to take before coarsening in the specified preconditioner method. Note that this key is only relevant to the SMG multigrid preconditioner.
pfset Solver.Linear.Preconditioner.SMG.NumPreRelax 2
integer Solver.Linear.Preconditioner.SMG.NumPostRelax 1 This key specifies the number of relaxations to take after coarsening in the specified preconditioner method. Note that this key is only relevant to the SMG multigrid preconditioner.
pfset Solver.Linear.Preconditioner.SMG.NumPostRelax 0
string Solver.Linear.Preconditioner.PFMG.RAPType NonGalerkin For the PFMG solver, this key specifies the Hypre RAP type. Valid values are Galerkin or NonGalerkin
pfset Solver.Linear.Preconditioner.PFMG.RAPType Galerkin
logical Solver.EvapTransFile False This key specifies specifies
that the Flux terms for Richards’ equation are read in from a .pfb
file. This file has [T^-1] units. Note this key is for a steady-state
flux and should not be used in conjunction with the transient key below.
pfset Solver.EvapTransFile True
logical Solver.EvapTransFileTransient False This key specifies
specifies that the Flux terms for Richards’ equation are read in from a
series of flux .pfb
file. Each file has [T^-1] units. Note this key s
hould not be used with the key above, only one of these keys should be set
to True
at a time, not both.
pfset Solver.EvapTransFileTransient True
string Solver.EvapTrans.FileName no default This key specifies
specifies filename for the distributed .pfb
file that contains the
flux values for Richards’ equation. This file has [T^-1] units.
For the steady-state option (Solver.EvapTransFile=True) this key
should be the complete filename. For the transient option
(Solver.EvapTransFileTransient=True then the filename is a header and
ParFlow will load one file per timestep, with the form filename.00000.pfb
.
pfset Solver.EvapTrans.FileName evap.trans.test.pfb
string Solver.LSM none This key specifies whether a land surface
model, such as CLM
, will be called each solver timestep. Choices
for this key include none and CLM. Note that CLM
must be compiled
and linked at runtime for this option to be active.
pfset Solver.LSM CLM
6.1.34. Spinup Options
These keys allow for reduced or dampened physics during model spinup or initialization. They are only intended for these initialization periods, not for regular runtime.
integer OverlandFlowSpinUp 0 This key specifies that a simplified form of the overland flow boundary condition (Equation [eq:overland_bc]) be used in place of the full equation. This formulation removes lateral flow and drives and ponded water pressures to zero. While this can be helpful in spinning up the subsurface, this is no longer coupled subsurface-surface flow. If set to zero (the default) this key behaves normally.
pfset OverlandFlowSpinUp 1
double OverlandFlowSpinUpDampP1 0.0 This key sets \(P_1\) and provides exponential dampening to the pressure relationship in the overland flow equation by adding the following term: \(P_2*exp(\psi*P_2)\)
pfset OverlandSpinupDampP1 10.0
double OverlandFlowSpinUpDampP2 0.0 This key sets \(P_2\) and provides exponential dampening to the pressure relationship in the overland flow equation adding the following term: \(P_2*exp(\psi*P_2)\)
pfset OverlandSpinupDampP2 0.1
6.1.35. CLM Solver Parameters
string Solver.CLM.Print1dOut False This key specifies
whether the CLM
one dimensional (averaged over each processor)
output file is written or not. Choices for this key include True and
False. Note that CLM
must be compiled and linked at runtime
for this option to be active.
pfset Solver.CLM.Print1dOut False
integer Solver.CLM.IstepStart 1 This key specifies the value of
the counter, istep in CLM
. This key primarily determines the start
of the output counter for CLM
. It is used to restart a run by setting
the key to the ending step of the previous run plus one. Note
that CLM
must be compiled and linked at runtime for this option to
be active.
pfset Solver.CLM.IstepStart 8761
String Solver.CLM.MetForcing no default This key specifies defines
whether 1D (uniform over the domain), 2D (spatially distributed) or 3D
(spatially distributed with multiple timesteps per .pfb
forcing file)
forcing data is used. Choices for this key are 1D, 2D and 3D. This key
has no default so the user must set it to 1D, 2D or 3D. Failure to set
this key will cause CLM
to still be run but with unpredictable values
causing CLM
to eventually crash. 1D meteorological forcing files
are text files with single columns for each variable and each timestep
per row, while 2D forcing files are distributed ParFlow binary files, one
for each variable and timestep. File names are specified in the
Solver.CLM.MetFileName variable below. Note that CLM
must be compiled
and linked at runtime for this option to be active.
pfset Solver.CLM.MetForcing 2D
String Solver.CLM.MetFileName no default This key specifies
defines the file name for 1D, 2D or 3D forcing data. 1D meteorological
forcing files are text files with single columns for each variable and
each timestep per row, while 2D and 3D forcing files are distributed
ParFlow binary files, one for each variable and timestep (2D) or one for
each variable and multiple timesteps (3D). Behavior of this key is
different for 1D and 2D and 3D cases, as sepcified by the
Solver.CLM.MetForcing key above. For 1D cases, it is the FULL FILE
NAME. Note that in this configuration, this forcing file is not
distributed, the user does not provide copies such
as narr.1hr.txt.0
, narr.1hr.txt.1
for each processor. ParFlow only
needs the single original file (e.g. narr.1hr.txt). For 2D cases, this key
is the BASE FILE NAME for the 2D forcing files, currently set to NLDAS,
with individual files determined as follows NLDAS.<variable>.<time step>.pfb
.
Where the <variable>
is the forcing variable and <timestep>
is the
integer file counter corresponding to istep above. Forcing is needed
for following variables:
- DSWR:
Downward Visible or Short-Wave radiation \([W/m^2]\).
- DLWR:
Downward Infa-Red or Long-Wave radiation \([W/m^2]\)
- APCP:
Precipitation rate \([mm/s]\)
- Temp:
Air temperature \([K]\)
- UGRD:
West-to-East or U-component of wind \([m/s]\)
- VGRD:
South-to-North or V-component of wind \([m/s]\)
- Press:
Atmospheric Pressure \([pa]\)
- SPFH:
Water-vapor specific humidity \([kg/kg]\) [clm_forcing]
Note that CLM
must be compiled and linked at runtime for this option to be active.
pfset Solver.CLM.MetFileName narr.1hr.txt
String Solver.CLM.MetFilePath no default This key specifies
defines the location of 1D, 2D or 3D forcing data. For 1D cases, this is
the path to a single forcing file (e.g. narr.1hr.txt
). For 2D and
3D cases, this is the path to the directory containing all forcing files.
Note that CLM
must be compiled and linked at runtime for this
option to be active.
pfset Solver.CLM.MetFilePath "path/to/met/forcing/data/"
integer Solver.CLM.MetFileNT no default This key specifies the number of timesteps per file for 3D forcing data.
pfset Solver.CLM.MetFileNT 24
string Solver.CLM.ForceVegetation False This key specifies whether
vegetation should be forced in CLM
. Currently this option only works
for 1D and 3D forcings, as specified by the key Solver.CLM.MetForcing
.
Choices for this key include True and False. Forced vegetation variables
are :
- LAI:
Leaf Area Index \([-]\)
- SAI:
Stem Area Index \([-]\)
- Z0M:
Aerodynamic roughness length \([m]\)
- DISPLA:
Displacement height \([m]\) [clm_forcing]
In the case of 1D meteorological forcings, CLM
requires four files
for vegetation time series and one vegetation map. The four files should
be named respectively lai.dat
, sai.dat
, z0m.dat
, displa.dat
.
They are ASCII files and contain 18 time-series columns (one per IGBP
vegetation class, and each timestep per row). The vegetation map should
be a properly distributed 2D ParFlow binary file (.pfb
) which contains
vegetation indices (from 1 to 18). The vegetation map filename is veg_map.pfb
.
ParFlow uses the vegetation map to pass to CLM
a 2D map for each
vegetation variable at each time step. In the case of 3D meteorological
forcings, ParFlow expects four distincts properly distributed ParFlow binary
file (.pfb
), the third dimension being the timesteps. The files should
be named LAI.pfb
, SAI.pfb
, Z0M.pfb
, DISPLA.pfb
. No
vegetation map is needed in this case.
pfset Solver.CLM.ForceVegetation True
string Solver.WriteSiloCLM False This key specifies whether the CLM
writes two dimensional binary output files to a silo binary format. This data
may be read in by VisIT and other visualization packages. Note that CLM
and silo must be compiled and linked at runtime for this option to be active.
These files are all written according to the standard format used for all ParFlow
variables, using the runname, and istep. Variables are either two-dimensional
or over the number of CLM
layers (default of ten).
pfset Solver.WriteSiloCLM True
The output variables are:
eflx_lh_tot
for latent heat flux total [W/m^2] using the silo variable LatentHeat;
eflx_lwrad_out
for outgoing long-wave radiation [W/m^2] using the silo variable LongWave;
eflx_sh_tot
for sensible heat flux total [W/m^2] using the silo variable SensibleHeat;
eflx_soil_grnd
for ground heat flux [W/m^2] using the silo variable GroundHeat;
qflx_evap_tot
for total evaporation [mm/s] using the silo variable EvaporationTotal;
qflx_evap_grnd
for ground evaporation without condensation [mm/s] using the silo
variable EvaporationGroundNoSublimation;
qflx_evap_soi
for soil evaporation [mm/s] using the silo variable EvaporationGround;
qflx_evap_veg
for vegetation evaporation [mm/s] using the silo variable EvaporationCanopy;
qflx_tran_veg
for vegetation transpiration [mm/s] using the silo variable Transpiration;
qflx_infl
for soil infiltration [mm/s] using the silo variable Infiltration;
swe_out
for snow water equivalent [mm] using the silo variable SWE;
t_grnd
for ground surface temperature [K] using the silo variable TemperatureGround; and
t_soil
for soil temperature over all layers [K] using the silo variable TemperatureSoil.
string Solver.PrintCLM False This key specifies whether the CLM
writes two dimensional
binary output files to a PFB
binary format. Note that CLM
must be compiled and linked
at runtime for this option to be active. These files are all written according to the
standard format used for all ParFlow variables, using the runname, and istep. Variables
are either two-dimensional or over the number of CLM
layers (default of ten).
pfset Solver.PrintCLM True
The output variables are:
eflx_lh_tot
for latent heat flux total [W/m^2] using the silo variable LatentHeat;
eflx_lwrad_out
for outgoing long-wave radiation [W/m^2] using the silo variable LongWave;
eflx_sh_tot
for sensible heat flux total [W/m^2] using the silo variable SensibleHeat;
eflx_soil_grnd
for ground heat flux [W/m^2] using the silo variable GroundHeat;
qflx_evap_tot
for total evaporation [mm/s] using the silo variable EvaporationTotal;
qflx_evap_grnd
for ground evaporation without sublimation [mm/s] using the silo
variable EvaporationGroundNoSublimation;
qflx_evap_soi
for soil evaporation [mm/s] using the silo variable EvaporationGround;
qflx_evap_veg
for vegetation evaporation [mm/s] using the silo variable EvaporationCanopy;
qflx_tran_veg
for vegetation transpiration [mm/s] using the silo variable Transpiration;
qflx_infl
for soil infiltration [mm/s] using the silo variable Infiltration;
swe_out
for snow water equivalent [mm] using the silo variable SWE;
t_grnd
for ground surface temperature [K] using the silo variable TemperatureGround; and
t_soil
for soil temperature over all layers [K] using the silo variable TemperatureSoil.
string Solver.WriteCLMBinary True This key specifies whether the CLM
writes two dimensional
binary output files in a generic binary format. Note that CLM
must be compiled and linked at
runtime for this option to be active.
pfset Solver.WriteCLMBinary False
string Solver.CLM.BinaryOutDir True This key specifies whether the CLM
writes
each set of two dimensional binary output files to a corresponding directory. These
directories my be created before ParFlow is run (using the tcl script, for example).
Choices for this key include True and False. Note that CLM
must be compiled and
linked at runtime for this option to be active.
pfset Solver.CLM.BinaryOutDir True
These directories are:
/qflx_top_soil
for soil flux;
/qflx_infl
for infiltration;
/qflx_evap_grnd
for ground evaporation;
/eflx_soil_grnd
for ground heat flux;
/qflx_evap_veg
for vegetation evaporation;
/eflx_sh_tot
for sensible heat flux;
/eflx_lh_tot
for latent heat flux;
/qflx_evap_tot
for total evaporation;
/t_grnd
for ground surface temperature;
/qflx_evap_soi
for soil evaporation;
/qflx_tran_veg
for vegetation transpiration;
/eflx_lwrad_out
for outgoing long-wave radiation;
/swe_out
for snow water equivalent; and
/diag_out
for diagnostics.
string Solver.CLM.CLMFileDir no default This key specifies what
directory all output from the CLM
is written to. This key may be
set to "./"
or ""
to write output to the ParFlow run directory.
This directory must be created before ParFlow is run. Note
that CLM
must be compiled and linked at runtime for this option
to be active.
pfset Solver.CLM.CLMFileDir "CLM_Output/"
integer Solver.CLM.CLMDumpInterval 1 This key specifies how often
output from the CLM
is written. This key is in integer multipliers
of the CLM
timestep. Note that CLM
must be compiled and linked
at runtime for this option to be active.
pfset Solver.CLM.CLMDumpInterval 2
string Solver.CLM.EvapBeta Linear This key specifies the form of
the bare soil evaporation \(\beta\) parameter in CLM
. The
valid types for this key are None, Linear, Cosine.
- None:
No beta formulation, \(\beta=1\).
- Linear:
\(\beta=\frac{\phi S-\phi S_{res}}{\phi-\phi S_{res}}\)
- Cosine:
\(\beta=\frac{1}{2}(1-\cos(\frac{(\phi -\phi S_{res})}{(\phi S-\phi S_{res})}\pi)\)
Note that \(S_{res}\) is specified by the key Solver.CLM.ResSat
below,
that beta is limited between zero and one and also that CLM
must
be compiled and linked at runtime for this option to be active.
pfset Solver.CLM.EvapBeta Linear
double Solver.CLM.ResSat 0.1 This key specifies the residual
saturation for the \(\beta\) function in CLM
specified above.
Note that CLM
must be compiled and linked at runtime for this
option to be active.
pfset Solver.CLM.ResSat 0.15
string Solver.CLM.VegWaterStress Saturation This key specifies the
form of the plant water stress function \(\beta_t\) parameter in CLM
.
The valid types for this key are None, Saturation, Pressure.
- None:
No transpiration water stress formulation, \(\beta_t=1\).
- Saturation:
\(\beta_t=\frac{\phi S -\phi S_{wp}}{\phi S_{fc}-\phi S_{wp}}\)
- Pressure:
\(\beta_t=\frac{P - P_{wp}}{P_{fc}-P_{wp}}\)
Note that the wilting point, \(S_{wp}\) or \(p_{wp}\), is
specified by the key Solver.CLM.WiltingPoint
below, that the
field capacity, S_{fc} or p_{fc}, is specified by the
key Solver.CLM.FieldCapacity
below, that beta_t is limited
between zero and one and also that CLM
must be compiled and
linked at runtime for this option to be active.
pfset Solver.CLM.VegWaterStress Pressure
double Solver.CLM.WiltingPoint 0.1 This key specifies the wilting
point for the \(\beta_t\) function in CLM
specified above. Note
that the units for this function are pressure [m] for a Pressure
formulation and saturation [-] for a Saturation formulation. Note
that CLM
must be compiled and linked at runtime for this option
to be active.
pfset Solver.CLM.WiltingPoint 0.15
double Solver.CLM.FieldCapacity 1.0 This key specifies the field
capacity for the \(\beta_t\) function in CLM
specified above.
Note that the units for this function are pressure [m] for a Pressure
formulation and saturation [-] for a Saturation formulation. Note
that CLM
must be compiled and linked at runtime for this option
to be active.
pfset Solver.CLM.FieldCapacity 0.95
string Solver.CLM.IrrigationTypes none This key specifies the form
of the irrigation in CLM
. The valid types for this key are none,
Spray, Drip, Instant.
pfset Solver.CLM.IrrigationTypes Drip
string Solver.CLM.IrrigationCycle Constant This key specifies the
cycle of the irrigation in CLM
. The valid types for this key are
Constant, Deficit. Note only Constant is currently implemented. Constant
cycle applies irrigation each day from IrrigationStartTime to
IrrigationStopTime in GMT.
pfset Solver.CLM.IrrigationCycle Constant
double Solver.CLM.IrrigationRate no default This key specifies the
rate of the irrigation in CLM
in [mm/s].
pfset Solver.CLM.IrrigationRate 10.
double Solver.CLM.IrrigationStartTime no default This key
specifies the start time of the irrigation in CLM
GMT.
pfset Solver.CLM.IrrigationStartTime 0800
double Solver.CLM.IrrigationStopTime no default This key specifies
the stop time of the irrigation in CLM
GMT.
pfset Solver.CLM.IrrigationStopTime 1200
double Solver.CLM.IrrigationThreshold 0.5 This key specifies the
threshold value for the irrigation in CLM
.
pfset Solver.CLM.IrrigationThreshold 0.2
integer Solver.CLM.ReuseCount 1 How many times to reuse a CLM
atmospheric forcing file input. For example timestep=1, reuse =1 is
normal behavior but reuse=2 and timestep=0.5 subdivides the time step
using the same CLM
input for both halves instead of needing two files.
This is particually useful for large, distributed runs when the user
wants to run ParFlow at a smaller timestep than the CLM
forcing. Forcing files will be re-used and total fluxes adjusted
accordingly without needing duplicate files.
pfset Solver.CLM.ReuseCount 5
string Solver.CLM.WriteLogs True When False, this disables writing of the CLM output log files for each processor. For example, in the clm.tcl test case, if this flag is added False, washita.output.txt.p and washita.para.out.dat.p (were p is the processor #) are not created, assuming washita is the run name.
pfset Solver.CLM.WriteLogs False
string Solver.CLM.WriteLastRST False Controls whether CLM restart files are sequentially written or whether a single file restart file name.00000.p is overwritten each time the restart file is output, where p is the processor number. If “True” only one file is written/overwritten and if “False” outputs are written more frequently. Compatible with DailyRST and ReuseCount; for the latter, outputs are written every n steps where n is the value of ReuseCount.
pfset Solver.CLM.WriteLastRST True
string Solver.CLM.DailyRST True Controls whether CLM writes daily restart files (default) or at every time step when set to False; outputs are numbered according to the istep from ParFlow. If ReuseCount=n, with n greater than 1, the output will be written every n steps (i.e. it still writes hourly restart files if your time step is 0.5 or 0.25, etc…). Fully compatible with WriteLastRST=False so that each daily output is overwritten to time 00000 in restart file name.00000.p where p is the processor number.
pfset Solver.CLM.DailyRST False
string Solver.CLM.SingleFile False Controls whether ParFlow writes
all CLM
output variables as a single file per time step. When “True”,
this combines the output of all the CLM output variables into a special
multi-layer PFB with the file extension “.C.pfb”. The first 13 layers
correspond to the 2-D CLM outputs and the remaining layers are the soil
temperatures in each layer. For example, a model with 4 soil layers will
create a SingleFile CLM output with 17 layers at each time step. The file
pseudo code is given below in ParFlow CLM Single Output Binary Files (.c.pfb) and
the variables and units are as specified in the multiple PFB
and SILO
formats as above.
pfset Solver.CLM.SingleFile True
integer Solver.CLM.RootZoneNZ 10 This key sets the number of soil
layers the ParFlow expects from CLM
. It will allocate and format all
the arrays for passing variables to and from CLM
accordingly. Note
that this does not set the soil layers in CLM
to do that the user
needs to change the value of the parameter nlevsoi
in the
file clm_varpar.F90
in the PARFLOW_DIR\pfsimulator\clm
directory to reflect the desired numnber of soil layers and recompile.
Most likely the key Solver.CLM.SoiLayer
, described below, will
also need to be changed.
pfset Solver.CLM.RootZoneNZ 4
integer Solver.CLM.SoiLayer 7 This key sets the soil layer, and
thus the soil depth, that CLM
uses for the seasonal temperature
adjustment for all leaf and stem area indices.
pfset Solver.CLM.SoiLayer 4
integer Solver.CLM.SoilLevels 10 This key sets the number of soil
levels for CLM
.
pfset Solver.CLM.SoilLevels 4
integer Solver.CLM.LakeLevels 10 This key sets the number of lake levels for CLM.
pfset Solver.CLM.LakeLevels 4
6.2. ParFlow NetCDF4 Parallel I/O
NetCDF4 parallel I/O is being implemented in ParFlow. As of now only
output capability is implemented. Input functionality will be added in
later version. Currently user has option of printing 3-D time varying
pressure or saturation or both in a single NetCDF file containing
multiple time steps. User should configure ParFlow(pfsimulatior
part) --with-netcdf
option and link the appropriate NetCDF4 library. Naming
convention of output files is analogues to binary file names. Following
options are available for NetCDF4 output along with various performance
tuning options. User is advised to explore NetCDF4 chunking and ROMIO
hints option for better I/O performance.
HDF5 Library version 1.8.16 or higher is required for NetCDF4 parallel I/O
integer NetCDF.NumStepsPerFile This key sets number of time steps user wishes to output in a NetCDF4 file. Once the time step count increases beyond this number, a new file is automatically created.
pfset NetCDF.NumStepsPerFile 5
string NetCDF.WritePressure False This key sets pressure variable to be written in NetCDF4 file.
pfset NetCDF.WritePressure True
string NetCDF.WriteSaturation False This key sets saturation variable to be written in NetCDF4 file.
pfset NetCDF.WriteSaturation True
string NetCDF.WriteMannings False This key sets Mannings coefficients to be written in NetCDF4 file.
pfset NetCDF.WriteMannings True
string NetCDF.WriteSubsurface False This key sets subsurface data(permeabilities, porosity, specific storage) to be written in NetCDF4 file.
pfset NetCDF.WriteSubsurface True
string NetCDF.WriteSlopes False This key sets x and y slopes to be written in NetCDF4 file.
pfset NetCDF.WriteSlopes True
string NetCDF.WriteMask False This key sets mask to be written in NetCDF4 file.
pfset NetCDF.WriteMask True
string NetCDF.WriteDZMultiplier False This key sets DZ multipliers to be written in NetCDF4 file.
pfset NetCDF.WriteDZMultiplier True
string NetCDF.WriteEvapTrans False This key sets Evaptrans to be written in NetCDF4 file.
pfset NetCDF.WriteEvapTrans True
string NetCDF.WriteEvapTransSum False This key sets Evaptrans sum to be written in NetCDF4 file.
pfset NetCDF.WriteEvapTransSum True
string NetCDF.WriteOverlandSum False This key sets overland sum to be written in NetCDF4 file.
pfset NetCDF.WriteOverlandSum True
string NetCDF.WriteOverlandBCFlux False This key sets overland bc flux to be written in NetCDF4 file.
pfset NetCDF.WriteOverlandBCFlux True
6.2.1. NetCDF4 Chunking
Chunking may have significant impact on I/O. If this key is not set, default chunking scheme will be used by NetCDF library. Chunks are hypercube(hyperslab) of any dimension. When chunking is used, chunks are written in single write operation which can reduce access times. For more information on chunking, refer to NetCDF4 user guide.
string NetCDF.Chunking False This key sets chunking for each time varying 3-D variable in NetCDF4 file.
pfset NetCDF.Chunking True
Following keys are used only when NetCDF.Chunking is set to true. These keys are used to set chunk sizes in x, y and z direction. A typical size of chunk in each direction should be equal to number of grid points in each direction for each processor. e.g. If we are using a grid of 400(x)X400(y)X30(z) with 2-D domain decomposition of 8X8, then each core has 50(x)X50(y)X30(z) grid points. These values can be used to set chunk sizes each direction. For unequal distribution, chunk sizes should as large as largest value of grid points on the processor. e.g. If one processor has grid distribution of 40(x)X40(y)X30(z) and another has 50(x)X50(y)X30(z), the later values should be used to set chunk sizes in each direction.
integer NetCDF.ChunkX None This key sets chunking size in x-direction.
pfset NetCDF.ChunkX 50
integer NetCDF.ChunkY None This key sets chunking size in y-direction.
pfset NetCDF.ChunkY 50
integer NetCDF.ChunkZ None This key sets chunking size in z-direction.
pfset NetCDF.ChunkZ 30
6.2.2. ROMIO Hints
ROMIO is a poratable MPI-IO implementation developed at Argonne National Laboratory, USA. Currently it is released as a part of MPICH. ROMIO sets hints to optimize I/O operations for MPI-IO layer through MPI_Info object. This object is passed on to NetCDF4 while creating a file. ROMIO hints are set in a text file in “key” and “value” pair. For correct settings contact your HPC site administrator. As in chunking, ROMIO hints can have significant performance impact on I/O.
string NetCDF.ROMIOhints None This key sets ROMIO hints file to be passed on to NetCDF4 interface.If this key is set, the file must be present and readable in experiment directory.
pfset NetCDF.ROMIOhints romio.hints
An example ROMIO hints file looks as follows.
romio_ds_write disable
romio_ds_read disable
romio_cb_write enable
romio_cb_read enable
cb_buffer_size 33554432
6.2.3. Node Level Collective I/O
A node level collective strategy has been implemented for I/O. One process on each compute node gathers the data, indices and counts from the participating processes on same compute node. All the root processes from each compute node open a parallel NetCDF4 file and write the data. e.g. If ParFlow is running on 3 compute nodes where each node consists of 24 processors(cores); only 3 I/O streams to filesystem would be opened by each root processor each compute node. This strategy could be particularly useful when ParFlow is running on large number of processors and every processor participating in I/O may create a bottleneck. Node level collective I/O is currently implemented for 2-D domain decomposition and variables Pressure and Saturation only. All the other ParFlow NetCDF output Tcl flags should be set to false(default value). CLM output is independently handled and not affected by this key. Moreover on speciality architectures, this may not be a portable feature. Users are advised to test this feature on their machine before putting into production.
string NetCDF.NodeLevelIO False This key sets flag for node level collective I/O.
pfset NetCDF.NodeLevelIO True
6.2.4. NetCDF4 Initial Conditions: Pressure
Analogues to ParFlow binary files, NetCDF4 based option can be used to set the initial conditions for pressure to be read from an “nc” file containing single time step of pressure. The name of the variable in “nc” file should be “pressure”. A sample NetCDF header of an initial condition file looks as follows. The names of the dimensions are not important. The order of dimensions is important e.g. (time, lev, lat, lon) or (time,z, y, x)
netcdf initial_condition {
dimensions:
x = 200 ;
y = 200 ;
z = 40 ;
time = UNLIMITED ; // (1 currently)
variables:
double time(time) ;
double pressure(time, z, y, x) ;
}
Node level collective I/O is currently not implemented for setting initial conditions.
string ICPressure.Type no default This key sets flag for initial conditions to be read from a NetCDF file.
pfset ICPressure.Type NCFile
pfset Geom.domain.ICPressure.FileName "initial_condition.nc"
6.2.5. NetCDF4 Slopes
NetCDF4 based option can be used slopes to be read from an “nc” file containing single time step of slope values. The name of the variable in “nc” file should be “slopex” and “slopey” A sample NetCDF header of slope file looks as follows. The names of the dimensions are not important. The order of dimensions is important e.g. (time, lat, lon) or (time, y, x)
netcdf slopex {
dimensions:
time = UNLIMITED ; // (1 currently)
lon = 41 ;
lat = 41 ;
variables:
double time(time) ;
double slopex(time, lat, lon) ;
}
netcdf slopey {
dimensions:
time = UNLIMITED ; // (1 currently)
lon = 41 ;
lat = 41 ;
variables:
double time(time) ;
double slopey(time, lat, lon) ;
}
The two NetCDF files can be merged into one single file and can be used with tcl flags. The variable names should be exactly as mentioned above. Please refer to “slopes.nc” under Little Washita test case. Node level collective I/O is currently not implemented for setting initial conditions.
string TopoSlopesX.Type no default This key sets flag for slopes in x direction to be read from a NetCDF file.
pfset TopoSlopesX.Type NCFile
pfset TopoSlopesX.FileName "slopex.nc"
string TopoSlopesY.Type no default This key sets flag for slopes in y direction to be read from a NetCDF file.
pfset TopoSlopesY.Type NCFile
pfset TopoSlopesy.FileName "slopey.nc"
6.2.6. NetCDF4 Transient EvapTrans Forcing
Following keys can be used for NetCDF4 based transient evaptrans forcing. The file should contain forcing for all time steps. For a given time step, if the forcing is null, zero values could be filled for the given time step in the “.nc” file. The format of the sample file looks as follows. The names of the dimensions are not important. The order of dimensions is important e.g. (time, lev, lat, lon) or (time,z, y, x)
netcdf evap_trans {
dimensions:
time = UNLIMITED ; // (1000 currently)
x = 72 ;
y = 72 ;
z = 3 ;
variables:
double evaptrans(time, z, y, x) ;
}
Node level collective I/O is currently not implemented for transient evaptrans forcing.
string NetCDF.EvapTransFileTransient False This key sets flag for transient evaptrans forcing to be read from a NetCDF file.
pfset NetCDF.EvapTransFileTransient True
string NetCDF.EvapTrans.FileName no default This key sets the name of the NetCDF transient evaptrans forcing file.
pfset NetCDF.EvapTrans.FileName "evap_trans.nc"
6.2.7. NetCDF4 CLM Output
Similar to ParFlow binary and silo, following keys can be used to write output CLM variables in a single NetCDF file containing multiple time steps.
integer NetCDF.CLMNumStepsPerFile None This key sets number of time steps to be written to a single NetCDF file.
pfset NetCDF.CLMNumStepsPerFile 24
string NetCDF.WriteCLM False This key sets CLM variables to be written in a NetCDF file.
pfset NetCDF.WriteCLM True
The output variables are:
eflx_lh_tot
for latent heat flux total [W/m^2] using the silo variable LatentHeat;
eflx_lwrad_out
for outgoing long-wave radiation [W/m^2] using the silo variable LongWave;
eflx_sh_tot
for sensible heat flux total [W/m^2] using the silo variable SensibleHeat;
eflx_soil_grnd
for ground heat flux [W/m^2] using the silo variable GroundHeat;
qflx_evap_tot
for total evaporation [mm/s] using the silo variable EvaporationTotal;
qflx_evap_grnd
for ground evaporation without condensation [mm/s] using the silo
variable EvaporationGroundNoSublimation;
qflx_evap_soi
for soil evaporation [mm/s] using the silo variable EvaporationGround;
qflx_evap_veg
for vegetation evaporation [mm/s] using the silo variable EvaporationCanopy;
qflx_tran_veg
for vegetation transpiration [mm/s] using the silo variable Transpiration;
qflx_infl
for soil infiltration [mm/s] using the silo variable Infiltration;
swe_out
for snow water equivalent [mm] using the silo variable SWE;
t_grnd
for ground surface temperature [K] using the silo variable TemperatureGround; and
t_soil
for soil temperature over all layers [K] using the silo variable TemperatureSoil.
6.2.8. NetCDF4 CLM Input/Forcing
Solver.CLM.IstepStart
netcdf metForcing {
dimensions:
lon = 41 ;
lat = 41 ;
time = UNLIMITED ; // (72 currently)
variables:
double time(time) ;
double APCP(time, lat, lon) ;
double DLWR(time, lat, lon) ;
double DSWR(time, lat, lon) ;
double Press(time, lat, lon) ;
double SPFH(time, lat, lon) ;
double Temp(time, lat, lon) ;
double UGRD(time, lat, lon) ;
double VGRD(time, lat, lon) ;
Note: While using NetCDF based CLM forcing, ``Solver.CLM.MetFileNT`` should be set to its default value of 1
string Solver.CLM.MetForcing no default This key sets meteorological forcing to be read from NetCDF file.
pfset Solver.CLM.MetForcing NC
Set the name of the input/forcing file as follows.
pfset Solver.CLM.MetFileName "metForcing.nc"
This file should be present in experiment directory. User may create soft links in experiment directory in case where data can not be moved.
6.2.9. NetCDF Testing Little Washita Test Case
The basic NetCDF functionality of output (pressure and saturation) and initial conditions (pressure) can be tested with following tcl script. CLM input/output functionality can also be tested with this case.
parflow/test/washita/tcl_scripts/LW_NetCDF_Test.tcl
This test case will be initialized with following initial condition file, slopes and meteorological forcing.
parflow/test/washita/parflow_input/press.init.nc
parflow/test/washita/parflow_input/slopes.nc
parflow/test/washita/clm_input/metForcing.nc
6.3. ParFlow Binary Files (.pfb)
The .pfb
file format is a binary file format which is used to store ParFlow
grid data. It is written as BIG ENDIAN binary bit ordering. The format
for the file is:
<double : X> <double : Y> <double : Z>
<integer : NX> <integer : NY> <integer : NZ>
<double : DX> <double : DY> <double : DZ>
<integer : num_subgrids>
FOR subgrid = 0 TO <num_subgrids> - 1
BEGIN
<integer : ix> <integer : iy> <integer : iz>
<integer : nx> <integer : ny> <integer : nz>
<integer : rx> <integer : ry> <integer : rz>
FOR k = iz TO iz + <nz> - 1
BEGIN
FOR j = iy TO iy + <ny> - 1
BEGIN
FOR i = ix TO ix + <nx> - 1
BEGIN
<double : data_ijk>
END
END
END
END
6.4. ParFlow CLM Single Output Binary Files (.c.pfb)
The .pfb
file format is a binary file format which is used to
store CLM
output data in a single file. It is written as
BIG ENDIAN binary bit ordering. The format for the file is:
<double : X> <double : Y> <double : Z>
<integer : NX> <integer : NY> <integer : NZ>
<double : DX> <double : DY> <double : DZ>
<integer : num_subgrids>
FOR subgrid = 0 TO <num_subgrids> - 1
BEGIN
<integer : ix> <integer : iy> <integer : iz>
<integer : nx> <integer : ny> <integer : nz>
<integer : rx> <integer : ry> <integer : rz>
FOR j = iy TO iy + <ny> - 1
BEGIN
FOR i = ix TO ix + <nx> - 1
BEGIN
eflx_lh_tot_ij
eflx_lwrad_out_ij
eflx_sh_tot_ij
eflx_soil_grnd_ij
qflx_evap_tot_ij
qflx_evap_grnd_ij
qflx_evap_soi_ij
qflx_evap_veg_ij
qflx_infl_ij
swe_out_ij
t_grnd_ij
IF (clm_irr_type == 1) qflx_qirr_ij
ELSE IF (clm_irr_type == 3) qflx_qirr_inst_ij
ELSE NULL
FOR k = 1 TO clm_nz
tsoil_ijk
END
END
END
END
6.5. ParFlow Scattered Binary Files (.pfsb)
The .pfsb
file format is a binary file format which is used to
store ParFlow grid data. This format is used when the grid data
is “scattered”, that is, when most of the data is 0. For data of
this type, the .pfsb
file format can reduce storage requirements
considerably. The format for the file is:
<double : X> <double : Y> <double : Z>
<integer : NX> <integer : NY> <integer : NZ>
<double : DX> <double : DY> <double : DZ>
<integer : num_subgrids>
FOR subgrid = 0 TO <num_subgrids> - 1
BEGIN
<integer : ix> <integer : iy> <integer : iz>
<integer : nx> <integer : ny> <integer : nz>
<integer : rx> <integer : ry> <integer : rz>
<integer : num_nonzero_data>
FOR k = iz TO iz + <nz> - 1
BEGIN
FOR j = iy TO iy + <ny> - 1
BEGIN
FOR i = ix TO ix + <nx> - 1
BEGIN
IF (<data_ijk> > tolerance)
BEGIN
<integer : i> <integer : j> <integer : k>
<double : data_ijk>
END
END
END
END
END
6.6. ParFlow Solid Files (.pfsol)
The .pfsol
file format is an ASCII file format which is used
to define 3D solids. The solids are represented by closed
triangulated surfaces, and surface “patches” may be associated
with each solid.
Note that unlike the user input files, the solid file cannot contain comment lines.
The format for the file is:
<integer : file_version_number>
<integer : num_vertices>
# Vertices
FOR vertex = 0 TO <num_vertices> - 1
BEGIN
<real : x> <real : y> <real : z>
END
# Solids
<integer : num_solids>
FOR solid = 0 TO <num_solids> - 1
BEGIN
#Triangles
<integer : num_triangles>
FOR triangle = 0 TO <num_triangles> - 1
BEGIN
<integer : v0> <integer : v1> <integer : v2>
END
# Patches
<integer : num_patches>
FOR patch = 0 TO <num_patches> - 1
BEGIN
<integer : num_patch_triangles>
FOR patch_triangle = 0 TO <num_patch_triangles> - 1
BEGIN
<integer : t>
END
END
END
The field <file_version_number>
is used to make file format
changes more manageable. The field <num_vertices>
specifies
the number of vertices to follow. The fields <x>
, <y>
,
and <z>
define the coordinate of a triangle vertex. The
field <num_solids>
specifies the number of solids to follow.
The field <num_triangles>
specifies the number of triangles
to follow. The fields <v0>
, <v1>
, and <v2>
are
vertex indexes that specify the 3 vertices of a triangle.
Note that the vertices for each triangle MUST be specified in
an order that makes the normal vector point outward from the
domain. The field <num_patches>
specifies the number of
surface patches to follow. The field num_patch_triangles
specifies the number of triangles indices to follow (these
triangles make up the surface patch). The field <t>
is
an index of a triangle on the solid solid
.
ParFlow .pfsol
files can be created from GMS .sol
files
using the utility gmssol2pfsol
located in the $PARFLOW_DIR/bin
directory. This conversion routine takes any number of
GMS .sol
files, concatenates the vertices of the solids defined
in the files, throws away duplicate vertices, then prints out
the .pfsol
file. Information relating the solid index in the
resulting .pfsol
file with the GMS names and material IDs are
printed to stdout.
6.7. ParFlow Well Output File (.wells)
A well output file is produced by ParFlow when wells are defined. The well output file contains information about the well data being used in the internal computations and accumulated statistics about the functioning of the wells.
The header section has the following format:
LINE
BEGIN
<real : BackgroundX>
<real : BackgroundY>
<real : BackgroundZ>
<integer : BackgroundNX>
<integer : BackgroundNY>
<integer : BackgroundNZ>
<real : BackgroundDX>
<real : BackgroundDY>
<real : BackgroundDZ>
END
LINE
BEGIN
<integer : number_of_phases>
<integer : number_of_components>
<integer : number_of_wells>
END
FOR well = 0 TO <number_of_wells> - 1
BEGIN
LINE
BEGIN
<integer : sequence_number>
END
LINE
BEGIN
<string : well_name>
END
LINE
BEGIN
<real : well_x_lower>
<real : well_y_lower>
<real : well_z_lower>
<real : well_x_upper>
<real : well_y_upper>
<real : well_z_upper>
<real : well_diameter>
END
LINE
BEGIN
<integer : well_type>
<integer : well_action>
END
END
The data section has the following format:
FOR time = 1 TO <number_of_time_intervals>
BEGIN
LINE
BEGIN
<real : time>
END
FOR well = 0 TO <number_of_wells> - 1
BEGIN
LINE
BEGIN
<integer : sequence_number>
END
LINE
BEGIN
<integer : SubgridIX>
<integer : SubgridIY>
<integer : SubgridIZ>
<integer : SubgridNX>
<integer : SubgridNY>
<integer : SubgridNZ>
<integer : SubgridRX>
<integer : SubgridRY>
<integer : SubgridRZ>
END
FOR well = 0 TO <number_of_wells> - 1
BEGIN
LINE
BEGIN
FOR phase = 0 TO <number_of_phases> - 1
BEGIN
<real : phase_value>
END
END
IF injection well
BEGIN
LINE
BEGIN
FOR phase = 0 TO <number_of_phases> - 1
BEGIN
<real : saturation_value>
END
END
LINE
BEGIN
FOR phase = 0 TO <number_of_phases> - 1
BEGIN
FOR component = 0 TO <number_of_components> - 1
BEGIN
<real : component_value>
END
END
END
END
LINE
BEGIN
FOR phase = 0 TO <number_of_phases> - 1
BEGIN
FOR component = 0 TO <number_of_components> - 1
BEGIN
<real : component_fraction>
END
END
END
LINE
BEGIN
FOR phase = 0 TO <number_of_phases> - 1
BEGIN
<real : phase_statistic>
END
END
LINE
BEGIN
FOR phase = 0 TO <number_of_phases> - 1
BEGIN
<real : saturation_statistic>
END
END
LINE
BEGIN
FOR phase = 0 TO <number_of_phases> - 1
BEGIN
FOR component = 0 TO <number_of_components> - 1
BEGIN
<real : component_statistic>
END
END
END
LINE
BEGIN
FOR phase = 0 TO <number_of_phases> - 1
BEGIN
FOR component = 0 TO <number_of_components> - 1
BEGIN
<real : concentration_data>
END
END
END
END
END
END
6.8. ParFlow Simple ASCII and Simple Binary Files (.sa and .sb)
The simple binary, .sa
, file format is an ASCII file format
which is used by pftools
to write out ParFlow grid data.
The simple binary, .sb
, file format is exactly the same,
just written as BIG ENDIAN binary bit ordering. The format
for the file is:
<integer : NX> <integer : NY> <integer : NZ>
FOR k = 0 TO <nz> - 1
BEGIN
FOR j = 0 TO <ny> - 1
BEGIN
FOR i = 0 TO <nx> - 1
BEGIN
<double : data_ijk>
END
END
END