Water Resources of the United States
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.
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.
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:
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'.