Physical processes and properties
Q: What options does SUTRA have for modeling groundwater flow?
A: SUTRA can model saturated or unsaturated, constant-density or density-dependent groundwater
flow in two or three spatial dimensions (2D or 3D). Flow can be steady-state or transient.
Q: What options does SUTRA have for modeling transport?
A: SUTRA can model transport of either a single solute or thermal energy (but
not both at the same time). Transport can be steady-state or transient (although
it must be transient if flow is transient).
Q: What options does SUTRA have for modeling geochemistry?
A: Although SUTRA is not designed to perform complex geochemical simulations, it does
simulate equilibrium adsorption (linear, Freundlich, or Langmuir) of solute onto
the solid, as well as zeroth-order and first-order production or decay of solute
and adsorbate.
Q: How do physical properties vary spatially in a SUTRA model?
A: There are three ways in which physical properties can vary spatially in SUTRA: nodewise (a value is assigned to each node); elementwise (a value is
assigned to each element); and as a function of pressure, concentration, or
temperature. Porosity varies nodewise. Permeabilities and dispersivities vary
elementwise. Unsaturated flow properties vary both nodewise and elementwise
insofar as a "region number" is assigned to each node and element; the
dependence of the unsaturated flow properties on pressure and region number is
programmed by the user in subroutine UNSAT. (See "Q:
How can I specify the unsaturated flow properties of my system?"
below.) Fluid density varies linearly with
solute concentration or temperature. For energy transport, fluid viscosity is a function of
temperature; for solute transport viscosity is constant and uniform. All
remaining physical properties are assigned a single, constant value throughout
the model.
Q: How can I specify the unsaturated flow properties of my system?
A: Unsaturated flow properties, i.e., saturation and relative permeability as
functions of pressure, are specified by programming subroutine UNSAT. (The SUTRA code must then be
recompiled using a Fortran 90 compiler.) In datasets 14B and 15B of the main
(".inp") input file, each node and element is assigned an
"unsaturated flow property region number," which is passed to
subroutine UNSAT and can be used to assign different properties to different
groups of nodes and elements. For details, see Section 7.5 and Appendix B of the SUTRA
documentation, and
"Unsaturated flow functions"
on the "Special topics" page.
Geometry and mesh structure
Q: Are there any constraints on the geometry of the model?
A: The geometry of the model is limited only by the ability to represent it
using a valid finite-element mesh. To learn about constraints on the mesh
structure, see "Q: Are there any
constraints on the structure of the mesh?" below.
Q: Are there any constraints on the structure
of the mesh?
A: A 2D SUTRA mesh is composed of quadrilateral finite elements. The four
corners of each element are called nodes. Adjacent elements share an edge and
the two nodes that define that edge. Avoid creating elements in which nodes coincide or edges cross, or in which the interior
angle between two adjacent edges is zero or greater than or equal to 180°.
While each 2D element has exactly four nodes, a node can be shared by any number of
elements.
A 3D SUTRA mesh is composed of generalized hexahedral finite elements. The
eight corners (nodes) of each element are connected by straight edges. Adjacent
elements share a face and the four nodes that define that face. Avoid creating elements in which nodes coincide or faces intersect,
or in which the interior angle between two adjacent faces is zero or greater
than or equal to 180°.
While each 3D element has exactly eight nodes, a node can be shared by any number of
elements. In SUTRA Version 2.0 (2D3D.1), 3D meshes had to be logically
rectangular. This restriction was removed beginning with SUTRA Version 2.1.
Matrix solver selection and convergence
Q: Which matrix solver should I use?
A: The DIRECT solver is a Gaussian elimination solver; it reduces the matrix
to upper-triangular form, then back-substitutes to obtain the solution. Its main
advantage is its reliability; it is guaranteed to provide an answer, the
accuracy of which is limited only by the rounding error inherent in the
calculations. The iterative solvers, CG, GMRES, and ORTHOMIN, each rely on an iterative
process to converge to the correct solution to within a specified tolerance.
Unfortunately, the iterations are not guaranteed to converge. Still, the
iterative solvers are useful because for large problems they typically require less
storage and compute an acceptable solution in less time than the
DIRECT solver.
Which type of solver is best for a given problem depends on the size of the
problem, its numerical "difficulty" (e.g., how ill-conditioned the
matrix is), the level of accuracy required in the solution, and the modeler's
time constraints. As a rule, the DIRECT solver is preferred for relatively small
problems, while the iterative solvers are preferred for large problems unless
severe convergence problems are encountered. The CG solver generally works well
for solving the flow equation without upstream weighting. The relative performance of the GMRES and
ORTHOMIN solvers is problem-dependent, and neither appears to consistently
outperform the other. Use whichever works best for your particular problem. If
the solver iterations have difficulty converging, the trouble usually lies with
the problem formulation; switching to a different iterative solver typically
does not help, although it is easy to try. See "Q:
My matrix solver iterations are not converging. What can I do?"
below for advice on handling convergence problems.
Q: My matrix solver iterations are not
converging. What can I do?
A: If the matrix solver iterations for the flow equation and/or the transport equation are
not converging (or are converging very poorly) early in a transient simulation, the most
likely cause is the way in which the initial conditions have been specified. If
possible, try to set initial pressures (p) and concentrations or temperatures
(U) that are mutually consistent. If your initial p and U fields are seriously
out of equilibrium with each other, they will try to come to equilibrium as
rapidly as possible as soon as the simulation begins. The sudden changes that
result can cause convergence problems. One way to ensure compatibility between
the initial p and U fields is to use a two-step approach. First, specify the
desired initial U field and a good guess for the corresponding steady-state p field, and perform a steady-state
SUTRA run (with a large convergence tolerance for U, so that any U
solution is accepted immediately). Then, copy the p solution from the
resulting restart (".rst") file into the original initial conditions
(".ics") file, leaving the original U field unchanged, and perform the
'real' SUTRA run using this modified ".ics" file. Using this
method allows the simulation to start with p and U fields that are
"equilibrated" with each other. This method is especially useful if
the early transient behavior of the system is not of great concern.
Problems that involve strong anisotropy or large permeability contrasts can
be particularly challenging for iterative solvers. For such problems,
"bootstrapping" often helps; the desired steady-state (or
near-steady-state) solution is obtained by performing a series of simulations,
with the results of each run serving as the initial guess for the next run. The
idea is to work your way from an "easy" problem to the original
"hard" problem, facilitating convergence by generating good initial
guesses all along the way. For example, if permeabilities range over many orders
of magnitude and you suspect this is causing convergence problems, reduce the
contrast (e.g., by decreasing the highest permeabilities, perhaps by several
orders of magnitude), until you obtain a problem the solver can handle. Then use
this solution as the initial guess for the next problem, in which the
permeabilities are set somewhat closer to their original values. Continue the
cycle until the permeabilities are back to their original values, and your
original problem is solved.
Another common cause of convergence failure is the use of time steps that are
too large. Reducing the time step size can help the simulation successfully get
through numerically difficult time steps, such as those during which boundary
stresses undergo a sudden change. Other possible causes are improper setting of
the boundary conductances GNUP and GNUU in dataset 5; failure to set reasonable
fluid and solid compressibilities in datasets 9 and 10 (transient runs only);
and insufficiently fine spatial discretization (see Section 7.2 of the SUTRA
documentation).
If all else fails, consider switching to the DIRECT solver, which does not
depend on an iterative scheme to converge to the solution, but which typically
requires considerably more storage and time to solve large problems than do
iterative solvers.
Boundary and initial conditions
Q: What types of boundary conditions does SUTRA support?
A: SUTRA supports the following four types of boundary conditions, all of
which are specified at nodes: fluid sources or sinks (dataset 17); solute mass or energy sources or sinks
(dataset 18); specified pressures (dataset 19); and specified solute concentrations or temperatures
(dataset 20). Boundary conditions are not restricted to nodes
at the physical boundary of the model; they can be specified at any
node. By default, there is no flux of mass or energy in or out of the model domain at any
node for which a boundary condition is not explicitly specified.
Q: How can I implement time-dependent boundary conditions?
A: In SUTRA Version 2.1 and earlier, time-dependent boundary conditions,
such as recharge that changes with time, are implemented by programming subroutine BCTIME.
Beginning with Version 2.2, time-dependent boundary conditions can alternatively be specified
without programming, using one or more input files.
When time-dependent boundary conditions are specified via subroutine BCTIME, the SUTRA code
must be recompiled using a Fortran 90 compiler. Nodes at which time-dependent boundary
conditions are specified are indicated by entering a negative node number in
dataset 17, 18, 19, or 20 (depending on the type of boundary condition) of the
main (".inp") input file.
Time-dependent boundary conditions specified via the optional ".bcs" files take precedence
over corresponding specifications in the ".inp" file or subroutine BCTIME. The input format is similar to that of
".inp" datasets 17 - 20.
For details, see Section 7.5 and Appendix B of the SUTRA
documentation, and "Time-dependent boundary
conditions" on the "Special
topics" page.
Q: What happens at nodes at which no boundary conditions are specified?
A: At nodes at which no boundary conditions have been specified, there is no
flux of mass or energy in or out of the model. Thus, a model boundary at which
no conditions have been specified is automatically a no-flux boundary.
Q: What do the boundary conductances GNUP and GNUU do, and why are they
important?
A: For reasons that have to do with matrix structure and efficient
calculation of boundary flows, specified pressure boundary conditions are
implemented somewhat indirectly in SUTRA. The flow equation is formulated as
though each pressure boundary node were connected to an external node set to the
specified pressure. The connection between each pressure boundary node and its
external partner is assigned a "conductance," GNUP. If GNUP is zero,
there is no hydraulic communication between each pressure boundary node and its
external partner, and the pressures computed at the pressure boundary nodes will be
independent of the pressure specifications. On the other hand, if GNUP is large
enough, each pressure boundary node is in perfect hydraulic communication with
its external partner, and the pressures computed at the pressure boundary nodes
will match the pressure specifications perfectly (i.e., to the full precision
provided by the computer). Intermediate values of GNUP produce partial matching
of the computed pressures to the pressure specifications.
If the only purpose of this numerical device were to enforce the pressure
boundary conditions as closely as possible, one would simply set GNUP to some
huge value. However, this device is also used to compute the flow in or out of
each pressure boundary node. The flow is equal to the conductance multiplied by
the difference between the pressure at the external node (i.e., the specified
pressure) and the pressure computed at the pressure boundary node. If GNUP is
set too low, the internal and external pressures match poorly, and an inaccurate
pressure solution is computed, which leads to an inaccurate calculation of
boundary flows. If GNUP is set too high, the internal and external pressures
match perfectly, and the boundary flows are computed as zero, regardless of
their true values. Therefore, GNUP must be chosen carefully to ensure that the
pressure specification is satisfied reasonably well, while leaving enough
mismatch to allow accurate calculation of boundary flows. The recommended value
of GNUP is that which causes the internal (computed) and external (specified)
pressures to match to the first 6 or 7 digits, and differ in the rest. This
leaves about half of the digits (in double precision) for computing boundary
flows. The result is that the pressure specifications are enforced to
single-precision accuracy, and the boundary flows are computed, in effect, in
single precision.
Note that SUTRA uses only a single value of GNUP for all the pressure
boundary nodes. If the magnitudes of the boundary flows do not vary widely over the model, a
single value of GNUP should suffice to satisfy the ideal "6 or 7 digit
matching" criterion reasonably well. However, if the magnitudes of the boundary flows vary
over many orders of magnitude, the ideal matching criterion cannot be satisfied
simultaneously at all pressure boundary nodes using one value of GNUP. In such
cases, the user must exercise judgment and set GNUP to satisfy the matching criterion most
closely in parts of the model where accuracy is most important.
Finding the right value of GNUP is a matter of trial and error. The procedure
is as follows: guess a value of GNUP (say, 0.01) and run SUTRA; look at the
output pressures and compare them to the specified pressures, noting the number
of digits to which they match; if they match to more than 6 or 7 digits, reduce
GNUP; if they match to fewer than 6 or 7 digits, increase GNUP; repeat as
necessary. As a rule, each order-of-magnitude change in GNUP changes the degree
of matching by one digit.
The conductance GNUU influences the enforcement of specified concentration or
temperature boundary conditions in the same way that GNUP influences the
enforcement of specified pressure conditions. GNUU should be chosen such that
the internal (computed) and external (specified) concentrations or temperatures
match to the first 6 or 7 digits, and differ in the rest.
The utility program CheckMatchBC is designed specifically to assist in
choosing appropriate values for GNUP and GNUU. See the SutraSuite
home page for a link to the latest version.
Q: What is the difference between a "cold" and a
"warm" start, and which one should I use?
There are two methods for initializing a SUTRA run: a cold start and a
warm start. A cold start is used to initialize a simulation that is being
run for the first time, i.e., one that is not simply a continuation of an
earlier run. As a rule, a cold start should be used for a continuation run if
the input parameters or boundary conditions have changed since the previous run.
A warm start is used to continue an earlier run as though it had never been
interrupted.
The initial conditions for a SUTRA run are read from the initial
conditions (".ics") input file. This file can be generated in
one of three ways:
- manually, using a text editor;
- automatically, using preprocessing software (e.g., SutraGUI or SutraPrep);
or
- by copying a restart (".rst") file, generated previously by SUTRA,
to the ".ics" file.
Any one of these methods can be used to generate a ".ics" file for
a SUTRA run that uses a cold start, as long as datasets 1 - 3 of this
file are specified completely and correctly. For a warm start, one must
use method 3 because the ".rst" file contains additional information
(which is recorded automatically by SUTRA, and which appears after
datasets 1 - 3) that allows a simulation to be continued as though it had never
been interrupted.
To use a cold start, set CREAD='COLD' in dataset 4 of the main input (".inp")
file. To use a warm start, set CREAD='WARM'.
|
U.S. Department of the Interior | U.S. Geological Survey
