U.S. Department of the Interior

U.S. Geological Survey

Techniques and Methods 6–A52

Prepared in cooperation with the 

Strategic Environmental Research and Development Program

SUTRA, a Model for Saturated-Unsaturated, VariableDensity Groundwater Flow with Solute or Energy 

Transport—Documentation of Generalized Boundary 

Conditions, a Modified Implementation of Specified 

Pressures and Concentrations or Temperatures, and the 

Lake Capability

Chapter 52 of

Section A, Groundwater

Book 6, Modeling Techniques
Cover. Illustration of lakes simulated using the SUTRA lake capability. Portions of the top surface of 

the groundwater model that are covered by lake water are shaded blue.
SUTRA, a Model for SaturatedUnsaturated, Variable-Density 

Groundwater Flow with Solute or Energy 

Transport—Documentation of Generalized 

Boundary Conditions, a Modified 

Implementation of Specified Pressures and 

Concentrations or Temperatures, and the 

Lake Capability

By Alden M. Provost and Clifford I. Voss

Chapter 52 of

Section A, Groundwater

Book 6, Modeling Techniques

Prepared in cooperation with the 

Strategic Environmental Research and Development Program

Techniques and Methods 6–A52

U.S. Department of the Interior

U.S. Geological Survey
U.S. Department of the Interior

DAVID BERNHARDT, Secretary

U.S. Geological Survey

James F. Reilly II, Director

U.S. Geological Survey, Reston, Virginia: 2019

For more information on the USGS—the Federal source for science about the Earth, its natural and living 

resources, natural hazards, and the environment—visit https://www.usgs.gov or call 1–888–ASK–USGS.

For an overview of USGS information products, including maps, imagery, and publications, 

visit https://store.usgs.gov.

Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the 

U.S. Government.

Although this information product, for the most part, is in the public domain, it also may contain copyrighted materials 

as noted in the text. Permission to reproduce copyrighted items must be secured from the copyright owner.

Suggested citation:

Provost, A.M., and Voss, C.I., 2019, SUTRA, a model for saturated-unsaturated, variable-density groundwater flow 

with solute or energy transport—Documentation of generalized boundary conditions, a modified implementation of 

specified pressures and concentrations or temperatures, and the lake capability: U.S. Geological Survey Techniques 

and Methods, book 6, chap. A52, 62 p., https://doi.org/10.3133/tm6A52.

ISSN 2328-7055 (online)
iii

Contents

Abstract ..........................................................................................................................................................1

Introduction.....................................................................................................................................................1

Chapter 1. Generalized Boundary Conditions .........................................................................................2

1.1. Piecewise-Linear Form of Generalized-Flow Conditions.......................................................2

1.2. Two-Point Specification of Generalized-Flow Conditions .....................................................3

1.3. Specification of Inflow and Outflow Concentrations or Temperatures...............................4

1.4. Commonly Used Forms of the Generalized-Flow Condition ..................................................4

1.4.1. Water-Table Elevation......................................................................................................4

1.4.2. General-Pressure Condition ...........................................................................................5

1.4.3. Well or Recharge..............................................................................................................7

1.4.4. Specified Pressure...........................................................................................................7

1.4.5. River and Drain..................................................................................................................7

1.4.6. Evapotranspiration ...........................................................................................................9

1.4.7. Evapotranspiration: Alternative Formulation for Unsaturated Flow ........................9

1.4.8. Seepage Face .................................................................................................................10

1.5. Linear Form of Generalized-Transport Conditions ................................................................10

1.6. Two-Point Specification of Generalized-Transport Conditions...........................................10

1.7. Commonly Used Forms of the Generalized-Transport Condition........................................11

1.7.1. Heat Transfer Across a Thermal Boundary Layer.....................................................11

1.7.2. Solute Mass Transfer Across a Boundary Layer ......................................................12

1.7.3. Specified Source or Sink of Solute Mass or Energy ................................................12

1.7.4. Specified Concentration or Temperature...................................................................12

1.8. Budget Output and Boundary-Condition Output Files ..........................................................13

1.8.1. The “.bcopg” File for Generalized-Flow Nodes.........................................................13

1.8.2. The “.bcoug” File for Generalized-Transport Nodes ................................................13

1.9. Assumptions and Limitations....................................................................................................14

1.10. Example Problem: Seepage Face (Dam with Vertical Sides)............................................14

1.10.1. Physical Setup ..............................................................................................................14

1.10.2. Objectives ......................................................................................................................14

1.10.3. Simulation Setup...........................................................................................................15

1.10.4. Parameters ....................................................................................................................15

1.10.5. Boundary Conditions....................................................................................................16

1.10.6. Initial Conditions ...........................................................................................................17

1.10.7. Results............................................................................................................................17

Chapter 2. Modified Implementation of Specified Pressures and Concentrations or

Temperatures ..................................................................................................................................18

2.1. Specified-Pressure Conditions: Scheme Used with Iterative Solvers ..............................18

2.2. Specified-Pressure Conditions: Scheme Used with the Direct Solver .............................19

2.3. Specified-Concentration or Temperature Conditions...........................................................19
iv

Chapter 3. Lake Capability.........................................................................................................................21

3.1. Preparatory Calculations ..........................................................................................................21

3.1.1. Identification of Potential Lakes ..................................................................................21

3.1.2. Lake Stage and Area as Functions of Water Volume...............................................22

3.2. Solution Procedure for Groundwater and Lakes ..................................................................23

3.2.1. Water Balances on Lakes .............................................................................................23

3.2.2. Solute Mass or Energy Balances on Lakes ...............................................................24

3.2.3. Effect of Lakes on Groundwater-Flow and Transport Boundary Conditions ........25

3.3. Special Features of the Lake Capability .................................................................................30

3.3.1. Limiting Where Lakes Can Form ..................................................................................30

3.3.2. Lake-Bottom Elevations.................................................................................................30

3.4. Some Algorithmic Details..........................................................................................................30

3.4.1. Identifying Potential Lakes Using the Hoshen-Kopelman Algorithm.....................30

3.4.2. Lake Hierarchy as a Binary Tree..................................................................................31

3.4.3. Determining the Fate of Spillover Using Steepest Descent ....................................32

3.4.4. Drying and Rewetting ....................................................................................................33

3.5. Lake-Related Output Files .........................................................................................................33

3.6. Assumptions and Limitations....................................................................................................33

3.7. Example Problem: Lakes............................................................................................................35

3.7.1. Physical Setup ................................................................................................................35

3.7.2. Objectives ........................................................................................................................35

3.7.3. Simulation Setup.............................................................................................................35

3.7.4. Parameters ......................................................................................................................36

3.7.5. Boundary Conditions......................................................................................................36

3.7.6. Initial Conditions .............................................................................................................37

3.7.7. Results..............................................................................................................................37

Acknowledgments.......................................................................................................................................40

References Cited..........................................................................................................................................40

Appendix 1. List of Symbols.......................................................................................................................43

Appendix 2. Flow Across a Conductive Layer........................................................................................48

Appendix 3. Input Data List........................................................................................................................49
v

Figures

1. Seven commonly used forms of the generalized-flow condition .........................................3

2. For a generalized-flow condition, inflow of fluid mass, Q, is a linear function of

pressure, p, subject to optional upper and lower limits on Q and (or) p .............................4

3. Two-point representation of a general-pressure boundary condition ...............................6

4. Alternative two-point representations of a specified-pressure boundary condition .......8

5. Two-point representation of a river boundary condition.......................................................8

6. Two-point representation of an evapotranspiration boundary condition...........................9

7. Two-point representation of a seepage-face boundary condition ....................................10

8. Three forms of the generalized-transport condition ............................................................10

9. For a generalized-transport condition, the gain of solute mass or energy, QU, is a

linear function of concentration or temperature, U..............................................................11

10. Boundary conditions and finite-element mesh for the seepage face example...............15

11. Near-steady-state SUTRA results for the seepage face example.....................................16

12. Nine potential lakes identified based on a given surface topography..............................22

13. Hierarchy of nine lakes identified based on the topography in figure 12 .........................22

14. Lake cells associated with the nodes in lake 3 identified in figure 12...............................22

15. Interaction of a lake with a fluid source or sink at model nodes .......................................26

16. Interaction of a lake with specified pressures at model nodes .........................................27

17. Interaction of a lake with solute or energy sources and sinks at model nodes ..............28

18. Interaction of a lake with specified concentrations or temperatures at model nodes ..29

19. The lake capability identifies potential lakes by tracking the formation of distinct

clusters of contiguous, submerged nodes on the top surface of the model as the

water stage varies from the highest to the lowest elevation on the surface...................31

20. The middle lake spills into the right-hand lake first because the right-hand sill is

lower than the left-hand sill ......................................................................................................32

21. In reality, a sufficiently strong source to the middle lake can cause it to spill over

into the two side lakes simultaneously ...................................................................................32

22. Finite-element mesh for the lake example .............................................................................35

23. Hierarchy of five lakes identified in the lake example..........................................................37

24. Initial groundwater-flow velocity vectors in the top layer of elements for the lake

example ........................................................................................................................................38

25. Areas of the model surface covered by lake water in the lake example..........................39

26. Lake-water concentration as a function of time in the lake example ...............................40

Tables

1. Summary of input parameter values for commonly used forms of the

generalized-flow condition .........................................................................................................5

2. Summary of input parameter values for two commonly used forms of the

generalized-transport condition...............................................................................................11
vi

Conversion Factors

International System of Units to U.S. customary units

Multiply By To obtain

Length

meter (m) 3.281 foot (ft) 

kilometer (km) 0.6214 mile (mi)

Area

square meter (m2) 0.0002471 acre 

square kilometer (km2) 247.1 acre

square centimeter (cm2) 0.001076 square foot (ft2)

square meter (m2) 10.76 square foot (ft2) 

square kilometer (km2) 0.3861 square mile (mi2)

Volume

liter (L) 0.2642 gallon (gal)

cubic meter (m3) 264.2 gallon (gal) 

cubic meter (m3) 35.31 cubic foot (ft3)

cubic kilometer (km3) 0.2399 cubic mile (mi3) 

cubic meter (m3) 0.0008107 acre-foot (acre-ft) 

Flow rate

cubic meter per second (m3/s) 70.07 acre-foot per day (acre-ft/d) 

cubic meter per year (m3/yr) 0.000811 acre-foot per year (acre-ft/yr) 

meter per second (m/s) 3.281 foot per second (ft/s) 

cubic meter per second per square

 kilometer [(m3/s)/km2]

91.49 cubic foot per second per square

 mile [(ft3/s)/mi2]

cubic meter per second (m3/s) 22.83 million gallons per day (Mgal/d) 

cubic meter per day per square

kilometer [(m3/d)/km2]

0.0006844 million gallons per day per square

mile [(Mgal/d)/mi2]

millimeter per year (mm/yr) 0.03937 inch per year (in/yr) 

Mass

gram (g) 0.03527 ounce, avoirdupois (oz)

kilogram (kg) 2.205 pound avoirdupois (lb)

Pressure

kilopascal (kPa) 0.009869 atmosphere, standard (atm)

kilopascal (kPa) 0.01 bar

kilopascal (kPa) 0.2961 inch of mercury at 60°F (in Hg)

kilopascal (kPa) 0.1450 pound-force per square inch (lbf/in2) 

Density

kilogram per cubic meter (kg/m3) 0.06242 pound per cubic foot (lb/ft3) 

gram per cubic centimeter (g/cm3) 62.4220 pound per cubic foot (lb/ft3) 

Energy

joule (J) 0.738 foot-pound (ft·lb)
vii

Multiply By To obtain

Specific capacity

liter per second per meter [(L/s)/m] 4.831 gallon per minute per foot [(gal/

min)/ft]

Hydraulic conductivity

meter per day (m/d) 3.281 foot per day (ft/d) 

Hydraulic gradient

meter per kilometer (m/km) 5.27983 foot per mile (ft/mi) 

Temperature in degrees Celsius (°C) may be converted to degrees Fahrenheit (°F) as follows:

°F = (1.8 × °C) + 32.

Temperature in degrees Fahrenheit (°F) may be converted to degrees Celsius (°C) as follows:

°C = (°F – 32) / 1.8.

SUTRA, a Model for Saturated-Unsaturated, VariableDensity Groundwater Flow with Solute or Energy 

Transport—Documentation of Generalized Boundary 

Conditions, a Modified Implementation of Specified 

Pressures and Concentrations or Temperatures, and the 

Lake Capability

By Alden M. Provost and Clifford I. Voss

Abstract 

Version 3.0 of the SUTRA groundwater modeling program offers three new capabilities: generalized boundary conditions, 

a modified implementation of specified pressures and concentrations or temperatures, and lakes. Two new types of “generalized” 

boundary conditions facilitate simulation of a wide range of hydrologic processes that interact with the groundwater model, such 

as rivers, drains, and evapotranspiration. For generalized-flow boundary conditions, gain (inflow) or loss (outflow) of fluid mass 

varies linearly with pressure, subject to optional upper and lower limits on flow and (or) pressure. For generalized-transport 

boundary conditions, gain or loss of solute mass or energy varies linearly with concentration or temperature, respectively. 

Two of the original types of SUTRA boundary conditions—specified-pressure and specified-concentration or temperature—

have been modified such that user-specified, conductance-like factors (known as GNUP and GNUU in previous versions of 

SUTRA) are no longer required. The new lake capability works with all types of SUTRA boundary conditions, including the 

new generalized boundary conditions, to enable simulation of the interaction of groundwater flow and transport with lake water 

“ponded” on the surface of a three-dimensional model. SUTRA uses the topography of the top surface of the model, or, optionally, user-specified lake-bottom elevations, to identify potential lakes automatically. Increases and decreases in lake stage can 

cause lakes to coalesce and divide, respectively. The lake capability may be used with saturated or unsaturated flow and solute or 

energy transport.

Introduction

This report complements the main SUTRA documentation (Voss and Provost, 2002). The report documents three new 

capabilities in SUTRA Version 3.0: generalized boundary conditions, a modified implementation of specified pressures and 

concentrations or temperatures, and lakes. Throughout this report, the terms “fluid” and “water” are used interchangeably to 

mean “water, including any dissolved solute.” The generic term “flow” refers to movement of fluid mass within the groundwater 

model, into or out of the model from or to an external source or sink, between groundwater and lakes, or between lakes.

Two new types of “generalized” boundary conditions facilitate simulation of a wide range of hydrologic processes that 

interact with the groundwater system, such as rivers, drains, and evapotranspiration. For “generalized-flow” conditions, gain 

(inflow) or loss (outflow) of fluid mass varies linearly with pressure, subject to optional upper and lower limits on flow and 

(or) pressure. For “generalized-transport” conditions, gain or loss of solute mass or energy varies linearly with concentration or 

temperature, respectively. The two new types of conditions are specified in datasets 21A and 21B of the “.inp” (main input) file 

and datasets 7A and 7B of the optional “.bcs” (time-dependent boundary condition) files.

The implementation of specified-pressure and specified-concentration or specified-temperature boundary conditions is modified in this version of SUTRA. Previously, the user was required to provide appropriate values for two conductance-like factors—GNUP and GNUU. These factors are no longer required as input and have been removed from dataset 5 of the “.inp” file.
2 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

The lake capability introduced in this version of SUTRA simulates the evolution of lakes as they exchange fluid mass and 

solute mass or thermal energy with groundwater, with each other, and with external sources or sinks. Given the topography of 

the top surface of the model, or, optionally, user-specified lake-bottom elevations, SUTRA identifies the lakes that can form by 

“ponding” of water on the top surface. Increases and decreases in lake stage can cause lakes to coalesce and divide, leading to a 

hierarchy of “parent” and “child” lakes that is automatically determined by SUTRA at the start of the simulation. On each time 

step, SUTRA solves the groundwater-flow and transport equations and updates lake stages and lake-water concentrations or temperatures. Lakes influence the groundwater simulation by interacting with user-specified flow and transport boundary conditions 

on the top surface of a model. The lake capability may be used with saturated or unsaturated flow, solute or energy transport, and 

all types of SUTRA boundary conditions, including the new generalized boundary conditions.

A list of input datasets created or modified to support the new and modified boundary conditions and the lake capability 

is given in appendix 3. A complete, up-to-date description of all SUTRA input datasets is maintained on the U.S. Geological 

Survey (USGS) Web site at https://doi.org/10.5066/P9PPEHHM.

Chapter 1. Generalized Boundary Conditions

Prior to Version 3.0, SUTRA supported only the following four types of boundary conditions, which are applied at nodes:

1. Sources and sinks of fluid mass,

2. Sources and sinks of solute mass or thermal energy,

3. Specified pressures, and

4. Specified concentrations or temperatures.

These four types, which are still supported in SUTRA 3.0, allow users to represent a variety of hydrologic processes that

interact with the groundwater system, such as pumping or injection of water, recharge, pressure due to surface water of known 

depth and density, addition of tracer, and variations in surface temperature.

To facilitate representation of a wider range of physical process, SUTRA 3.0 offers two additional types of boundary 

conditions, which are also applied at nodes:

1. Generalized-flow conditions – sources and sinks of fluid mass that depend on pressure, and,

2. Generalized-transport conditions – sources and sinks of solute mass or energy that depend on concentration or temperature, respectively.

The two new types of boundary conditions are “generalized” in the sense that they include the four original types as special

cases while offering additional capabilities. Like the four original types, the two new types can be made time-dependent using 

“.bcs” input files.

Like sources and sinks of fluid mass and specified pressures, generalized-flow conditions require specification of 

concentrations or temperatures associated with fluid inflows to the model. In addition, generalized-flow conditions require 

specification of concentrations or temperatures associated with fluid outflows from the model. The latter specification enables 

simulation of processes such as evaporation, in which fluid may leave the model at a concentration or temperature different from 

that of the groundwater.

1.1. Piecewise-Linear Form of Generalized-Flow Conditions

Figure 1 shows schematic representations of seven idealized, pressure-dependent boundary conditions: general pressure, 

well or recharge, specified pressure, river, drain, evapotranspiration (ET), and seepage face. By specifying such functions 

using the generalized-flow boundary condition, the user can create these and other boundary-condition variations. The first 

six conditions are analogous to commonly used head-dependent conditions described by Harbaugh (2005). These conditions 

represent boundary-condition processes that may be specified at any node in the finite-element model mesh. By convention, 

fluid-mass gain by (inflow to) the model is positive and fluid-mass loss (outflow) is negative.

In each case in figure 1, the relation between flow, Q [M/s], into or out of the model at the boundary-condition node, and 

simulated pressure, p [M/(L·s2)], at the node is piecewise-linear; it is defined by a combination of horizontal, vertical, and (or) 

sloping straight-line segments. For a general-pressure condition (fig. 1A), the sloping segment of the graph characterizes a 

conductive connection between the simulated pressure at the boundary-condition node and a specified external pressure. For a 
Chapter 1. Generalized Boundary Conditions 3

well or recharge condition (fig. 1B), flow is specified independently of pressure, as indicated by the horizontal line in the graph. 

For a specified-pressure condition (fig. 1C), pressure is specified independently of flow, as indicated by the vertical line in the 

graph. For a river condition (fig. 1D), the conductive connection between the boundary-condition node and the external pressure imposed by the overlying river water is conceptualized as occurring through a river bed. The horizontal segment of the plot 

indicates that, for a given river stage, the rate of seepage from the river becomes independent of pressure when the water table 

is below the river bed (Harbaugh, 2005). (The correspondence between water-table elevation and pressure assumed in formulating some generalized-flow conditions is discussed in section 1.4.1., “Water-Table Elevation.”) A drain condition (fig. 1E) is a 

special case of a river condition in which surface water does not accumulate in the channel and recharge of surface water to the 

groundwater cannot occur. For an evapotranspiration (ET) condition (fig. 1F), the sloping segment of the graph characterizes 

the change in the rate of fluid loss by ET as the pressure (and, hence, the water-table elevation) changes. When the water table is 

above a certain depth, the rate of fluid loss by ET is at its maximum possible value, and when the water table is below a certain 

depth (the “extinction depth”), the rate of fluid loss by ET is zero, as indicated by the horizontal segments of the graph. (An 

alternative interpretation of figure 1F is presented in section 1.4.7. “Evapotranspiration: Alternative Formulation for Unsaturated 

Flow.”) For a seepage-face condition, as defined for SUTRA (fig. 1G), the vertical segment of the graph below the horizontal 

axis represents conditions along the lower portion of the seepage face, where the pressure is constrained to be zero (atmospheric 

pressure) and water discharges from the model. The vertical segment of the graph above the horizontal axis represents conditions along the seepage face where the pressure is constrained to be zero and water recharges the model. (Recharge along a 

seepage face can occur along its uppermost portion when the face is not vertical.) The horizontal segment of the graph represents 

recharge at a specified rate (p < 0) in model regions above the uppermost seepage face where unsaturated conditions prevail. If 

the face is vertical, the recharge rate is zero. Each of the conditions illustrated in figure 1 is discussed in more detail in section 

1.4. “Commonly Used Forms of the Generalized-Flow Condition.”

1.2. Two-Point Specification of Generalized-Flow Conditions

A generalized-flow condition at a node requires specification of a linear dependence of flow (fluid mass gain), Q, on 

simulated pressure at the node, p. Optional upper and (or) lower limits can be set on flow and (or) pressure. Figure 2 illustrates 

the two-point scheme used to specify both the sloping linear segment and the optional limits for a generalized-flow condition. 

User-specified points (p1, Q1) and (p2, Q2) define the linear relation between Q and p. These also serve as points at which the user 

can optionally limit either Q or p (but not both at the same point). The values p1 and p2 can be used as lower and upper limiting values of p, respectively. The values Q1

 and Q2 can be used as upper and lower limiting values of Q, respectively. All of the 

boundary conditions shown in figure 1 are created by the user setting appropriate values of (p1, Q1) and (p2, Q2) and optionally 

specifying limits at these points, thereby defining the piecewise-linear relation for each particular condition.

For most realistic hydrologic processes that might be represented by a generalized-flow boundary condition, the linear 

segment will have a negative slope as shown in figure 1; as simulated pressure at the node increases relative to the external 

pressure, the flow at the node decreases (fluid mass gain from an external source becomes less positive, or fluid mass loss to an 

external sink becomes more negative). Thus, SUTRA requires p2 > p1 and Q2 ≤ Q1, which implies a negative or zero slope. 

For the sake of notational simplicity, the symbols p1, p2, Q1, and Q2 are used throughout Chapter 1 to represent the 

user-specified pressures and flows that define a generalized-flow boundary condition. Note that in input datasets that define 

Figure 1. Seven commonly 

used forms of the generalizedflow condition: A, general 

pressure; B, well or recharge; 

C, specified pressure; D, river; 

E, drain; F, evapotranspiration; 

and G, seepage face. Q is the 

gain of fluid mass at the node 

where this boundary condition 

is specified [M/s] (see 

appendix 2), and p is simulated 

pressure at that node [M/(L·s2)]. 

The symbol “0” represents the 

origin, (p, Q) = (0, 0).
4 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

generalized-flow conditions (“.inp” dataset 21A and “.bcs” 

dataset 7A), p1, p2, Q1, and Q2 are represented by the more 

descriptive symbols pBG1, pBG2, QPBG1, and QPBG2, respectively.

1.3. Specification of Inflow and Outflow 

Concentrations or Temperatures

For all types of flow-related SUTRA boundary 

conditions, including generalized-flow conditions, the user 

specifies the concentration or temperature of any inflow to 

the model that might occur at each boundary node. For a 

specified-source/sink or specified-pressure condition at a node, 

any outflow automatically has the groundwater concentration or temperature computed at the node. Generalized-flow 

conditions offer more flexibility by allowing the user to set the 

outflow concentration or temperature. Two options are available for specifying the value of the outflow concentration or 

temperature. The first option is to specify the value directly. 

For example, to simulate ET of pure water from solute-laden 

groundwater, the outflow concentration would be set to zero. 

The second option is to specify the outflow value relative to 

the value computed by SUTRA at the boundary node. The 

relative value can be positive, zero, or negative. A positive 

(negative) relative value indicates a concentration or temperature greater than (less than) the value computed by SUTRA at the 

boundary node. For example, to simulate groundwater discharge to a drain, the outflow concentration would be set to a relative

value of zero so that the discharging water has the same concentration as the groundwater at the boundary node.

1.4. Commonly Used Forms of the Generalized-Flow Condition

The generalized-flow condition provides a framework for defining commonly used boundary conditions on groundwater flow. 

Guidance on how to formulate a generalized-flow condition for the forms shown in figure 1 is provided below. Input parameter 

values for general-pressure conditions, rivers, drains, ET, and seepage faces are summarized in table 1. Although specified-pressure 

and pressure-independent conditions, such as recharge and wells, are most easily handled using other types of SUTRA boundary 

conditions, instructions for formulating such conditions within the generalized-flow framework are given in the main text.

1.4.1. Water-Table Elevation

In the discussion below, the river, drain, and ET conditions are conceptualized in terms of water-table elevation. In a drain, 

for example, flow is zero when the water table drops below the drain elevation. However, because a SUTRA generalized-flow 

condition expresses flow as a function of pressure at the boundary-condition node, water-table elevation must be estimated as 

a function of pressure at the boundary-condition node. For this purpose, the pressure profile below land surface is assumed to 

be hydrostatic:

p z p g z z wt wt ( ) = + r ( ) − (1)

where p(z) is the estimated pressure at a particular boundary-condition node, which is at elevation z [L], pwt is the pressure at 

the water table, which is at elevation zwt [L], ρ [M/L3] is fluid density, and g [L/s2] is the acceleration of gravity. Assuming the 

pressure at the water table is zero (pwt = 0), rearrangement of equation (1) yields an approximate relation between water-table 

elevation, zwt, and simulated pressure, p, at the boundary-condition node:

z z p

g wt = +

r

(2)

Figure 2. For a generalized-flow condition, inflow of fluid mass, 

Q, is a linear function of pressure, p, subject to optional upper and 

lower limits (dashed lines) on Q and (or) p. By convention, positive 

values of Q denote inflows and negative values denote outflows.
Chapter 1. Generalized Boundary Conditions 5

When p < 0, p = 0, or p > 0 at the boundary-condition node, the water-table elevation is below, exactly at, or above the 

boundary-condition node, respectively.

The river, drain, and ET forms of the generalized-flow condition are simplified representations of real hydrologic processes. 

Although SUTRA is capable of simulating detailed saturation profiles for saturated-unsaturated flow, the river, drain, and 

ET conditions described below do not make use of the details of unsaturated zone physics. Rather, they use equation (2) to 

formulate approximate, piecewise-linear relations between flow and pressure at boundary-condition nodes. The user should 

remain mindful of this approximation when constructing models and interpreting results. (In the case of ET, an alternative 

formulation that does take into account variable saturation is presented in section 1.4.7, “Evapotranspiration: Alternative 

Formulation for Unsaturated Flow.”)

1.4.2. General-Pressure Condition

A general-pressure condition is implemented by specifying two points, (p1, Q1) and (p2, Q2), with no limits on p or Q

(fig. 3). The concentrations or temperatures associated with possible inflow and outflow at the general-pressure node are set by 

the user according to the hydrologic-geochemical process being represented. For outflow, the concentration or temperature is 

typically chosen to be that of the groundwater, unless a chemical or physical process somehow modifies it.

The common case in which the boundary condition represents flow across a conductive layer is described below. Within 

the conductive layer, permeability, fluid density and viscosity, and gravity are assumed to be uniform. Flow is driven by the 

hydraulic head difference across the layer; under hydrostatic conditions (no head difference), no flow crosses the layer. Because 

SUTRA calculates flows in terms of pressure rather than head, this concept and the equations below are presented in terms of 

pressure. Details of the derivation are presented in appendix 2.

Table 1. Summary of input parameter values for commonly used forms of the generalized-flow condition. Each form is discussed in 

detail and symbols are defined in the main text.

Special cases1 p1 Q1 p2 Q2 Notes

General pressure pext Qfd

C

g

p p Q lay

( ) ext − 2 + fd If free drainage is negligible compared to pressuredriven flow, Qfd = 0 . See equations (3)–(8). Note that 

Clay ≥ 0 .

River, Drain pext Qmax

C

g

p p Q lay

( ) ext − 2 + max For a river, use p gd ext = r riv and Q Q max f = +d r Q iv . 

Set an upper limit on flow at Q1.

For a drain, use pext = 0 and Qmax = 0 . Set an upper 

limit on flow at Q1.

See equations (3)–(8) and (9)–(11). Note that Clay ≥ 0 .

Evapotranspiration

−rgdex 0 −rgdmax q A ET

max Set upper and lower limits on flow at Q1

 and Q2, respectively. Note that dex > 0 , dmax ≥ 0 , and qET

max < 0 > 

by convention, and d d ex > max .

See equations (13) and (14).

Seepage face < 0 Qinf 0 Qinf Set an upper limit on pressure at p2

. For a vertical seepage 

face, Qinf = 0 .

See equation (15). 1

Representation of pressure-independent sources and sinks (such as wells and recharge) and specified pressures using generalized-flow conditions is discussed in

the main text. However, the simplest way to implement pressure-independent sources and sinks is to use “.inp” dataset 17 and “.bcs” dataset 3, and the simplest 

way to implement specified pressures is to use “.inp” dataset 19 and “.bcs” dataset 5.

p > ext

p > ext
6 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

The expression for flow, Q [M/s], across the layer consists of two terms that reflect the two forces—pressure and 

gravity—that drive the flow:

Q

C

g

p p Q lay = − ( ) ext + fd (3)

The first term on the right-hand side of equation (3) represents 

flow driven by the pressure difference between the boundary 

node at pressure p and an external point (on the other side of 

the layer) at pressure pext. The conductance of the layer, Clay

[L2/s], is defined by

 (4)

where

 (5)

is the hydraulic conductivity [L/s] of the layer, L [L] is its 

thickness, k [L2] is its permeability, ρ is fluid density [M/L3], 

μ is fluid viscosity [M/(L·s)], g [L/s2] is the acceleration of 

gravity, and A [L2] is the area of the conductive layer through 

which flow passes to the boundary node (see fig. B.1 in Voss 

and Provost, 2002). The second term on the right-hand side of 

equation (3), Qfd [M/s], represents gravity-driven flow (that 

is, the rate of free drainage) that would occur if there were no 

pressure difference across the conductive layer. This is given by

 (6)

where gn [L/s2] is the component of gravity across (normal to) the layer. If the layer is inclined at an angle θ above the 

horizontal, g g n = cos  ; if the layer is horizontal, g g n = 1 ; and if the layer is vertical, g g n = 0 . If pressure within the layer 

is hydrostatic ( p p = + ext g Ln ), the driving forces due to pressure and gravity cancel each other exactly in equation (3), and no 

flow crosses the layer (Q = 0).

Any two points that lie on the line represented by equation (3) may be used to define a general-pressure condition that 

represents flow through a conductive layer as described above. One may set, for example,

p p Q Q p p Q

C

g ext fd ext p p Q lay

1 1 = = 2 2 > = ( ) ext − 2 + fd , , , and (7)

If free drainage is negligible compared to pressure-driven flow across the layer (or zero, as for a vertical conductive layer), then 

Qfd ≈ 0 , and equation (7) may be approximated by

p p Q p p Q

C

g

p p ext ext

lay

1 1 2 2 ext 2 = = , , 0 > = , and ( ) − (8)

C KA

L lay =

K k g = ρ

μ

Q KA g

g fd

n = 









 

Figure 3. Two-point representation of a general-pressure 

boundary condition. No limits are set on p or Q. The case in 

which the general-pressure condition represents flow through 

a conductive medium is described in the main text. A pressureindependent source or sink of fluid is a special case in which 

Q1 = Q2. For a pumping well, Q1 = Q2 < 0, and for an injection well or 

recharge, Q1 = Q2 > 0. Parameter values are summarized in table 1. 
Chapter 1. Generalized Boundary Conditions 7

1.4.3. Well or Recharge

The simplest way to represent a pressure-independent fluid source or sink, such as recharge or a pumping well, is through 

“.inp” dataset 17 and, for time-dependent conditions, the corresponding “.bcs” dataset 3. However, such a condition can also 

be formulated as a special case of a general-pressure condition with Q1 = Q2 so that flow does not vary with pressure. Flow 

(fluid mass gain) is negative for a pumping well and positive for an injection well or recharge. The concentration or temperature 

associated with inflow is set to that of the water being recharged or injected. The concentration or temperature associated 

with outflow is set to that of the groundwater. Evaporation of pure (solute-free) water from solute-laden groundwater can be 

simulated by setting the outflow concentration to zero.

1.4.4. Specified Pressure

The simplest way to represent a specified-pressure boundary condition is through “.inp” dataset 19 and, for time-dependent 

conditions, the corresponding “.bcs” dataset 5. However, the generalized-flow framework allows two alternative formulations 

for specified-pressure boundary conditions. The first alternative (fig. 4A) is to set p1 and p2 nearly equal so that the slope is very 

large. This is essentially how specified-pressure conditions were formulated in earlier versions of SUTRA. The disadvantage 

of this method is that the slope must be chosen by trial and error such that SUTRA computes the pressure and the corresponding flow to approximately equal precision. For details, see section 3.5 and the input data listing for “.inp” dataset 5, both in 

the main SUTRA documentation (Voss and Provost, 2002). The second alternative (fig. 4B) is to set a limit on pressure at 

p = p1, set p2 greater than p1, and set Q1 and Q2 to values that will never be attained in the simulation, so that p is always held at 

its limiting value.

1.4.5. River and Drain

A river boundary condition represents the hydraulic connection, through the river bed, between the boundary-condition 

node at the base of the river bed and the overlying river water. When the pressure computed at the boundary-condition node is 

negative, the water-table elevation is below the river bed, and the rate of leakage through the river bed no longer depends on 

pressure (Harbaugh, 2005). Thus, the river condition is a variation on the general-pressure condition (3) in which p gd ext = r riv is 

the pressure exerted by the overlying river water, where driv [L] is the depth of the river water, atmospheric pressure is zero, and 

an upper limit is set on flow (fig. 5). The maximum flow rate, Qmax [M/s], is the sum of two contributions:

Q Q max f = +d r Q iv (9)

where Qfd is the free drainage given by equation (6), and Qriv [M/s] is the flow driven by the additional head supplied by the 

overlying river water,

Q KA d

L riv

riv = 



 

 r  (10)

Thus, for a river condition, one may set

p p Q Q p p Q

C

g ext max ext p p Q lay

1 1 = = 2 1 > = ( ) ext − 2 + max , , , and (11)

with p gd ext = r riv , Qmax given by equation (9), and a limit on Q set at Q1.

The drain boundary condition at a node is a special case of the river condition in which no water collects in the channel 

( pext = 0 ) and fluid mass can only be lost from the drain node; groundwater recharge is prohibited from a drain ( Qmax = 0 ). 

A river-type boundary condition can also be used to approximate the hydraulic behavior of a perforated or porous pipe buried 

beneath the ground. In that case, the pipe wall is analogous to the river bed.
8 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

Figure 5. Two-point representation of a river boundary condition. 

An upper limit of Q1 is set on Q when the pressure becomes low 

enough that the water-table elevation is below the river bed, and 

the rate of leakage through the river bed no longer depends on 

pressure. A drain boundary condition is similar, but with p1 = 0 and 

Q1 = 0. Parameter values for rivers and drains are summarized in 

table 1.

Figure 4. Alternative two-point representations of a specified-pressure boundary condition. A, The values of p1 and p2 are nearly 

equal, so the slope is very large. B, A limit is set at p = p1, p2 is set greater than p1, and Q1 and Q2 are set to a very large value outside the 

expected range of flows, so the limit on p is always active.
Chapter 1. Generalized Boundary Conditions 9

1.4.6. Evapotranspiration

By analogy with the head-dependent evapotranspiration (ET) condition described by Harbaugh (2005), the rate of fluid 

loss, Q, due to ET at an ET boundary-condition node depends linearly on simulated pressure, p, at the node when the water-table 

elevation is between an upper and a lower bound, that is, when

(12)

where z is the elevation of the ET boundary-condition node at land surface and zwt is the water-table elevation estimated using 

equation (2). The upper bound on zwt corresponds to the depth, dmax [L], above which the rate of fluid loss is at its maximum 

(most negative) value, and the lower bound corresponds to the “extinction depth,” dex [L], below which the rate of fluid loss is 

zero. Depths are defined as positive distances from land surface, and d d ex > max . Using equation (2), equation (12) can be reexpressed in terms of pressure:

− ≤ r r gd p g ≤ − d ex max (13)

resulting in a piecewise-linear dependence of Q on p (fig. 6). Thus, the ET condition can be implemented by setting

p gd Q ex p gd Q max E q A T

max

1 1 2 2 = −r r , , = = 0 − = , and (14)

where qET

max is the maximum (most negative) ET flux [M/(L2

·s)] and A is the area [L2] associated with the ET boundary-condition 

node. In accordance with SUTRA convention, qET

max is defined as negative because it represents a loss of fluid from the model.

For solute transport, the concentration of the evapotranspired (outflowing) water is typically assumed to be zero, although 

any desired value can be specified. For energy transport, the temperature associated with the evapotranspired water depends on 

the process being simulated. If only transpiration by plants 

(not evaporation) is being simulated, setting the temperature of the transpired water to that of the groundwater may 

be the most appropriate choice. If only direct evaporation 

(not transpiration) is being simulated, the latent heat of 

vaporization, ∆Hvap [E/M], can be taken into account by 

setting the associated outflow temperature to a relative value 

of ∆H c vap w , where cw [E/(M·oC)] is the specific heat of 

water. (Relative specification of outflow temperatures is 

discussed is section 1.3, “Specification of Inflow and Outflow 

Concentrations or Temperatures.”) The concentration or temperature associated with inflow is irrelevant in the case of ET, 

since inflow cannot occur.

1.4.7. Evapotranspiration: Alternative 

Formulation for Unsaturated Flow

In the formulation of the ET condition presented in section 1.4.6, the loss of water due to ET represents the total loss 

integrated over the depth of the root zone. As an approximation, this total loss of water is assigned to the node at the top 

surface of the model. Such an approximation is especially 

useful when the details of the unsaturated zone are considered 

unimportant and the flow system is modeled as saturated. In 

that case, the simulated pressure at the top-surface node provides an estimate of the water-table elevation according to (2).

If unsaturated flow is being simulated explicitly, an ET 

condition can be applied at multiple nodes at or below the 

z d z z d − ≤ ex wt ≤ − max

Figure 6. Two-point representation of an evapotranspiration 

(ET) boundary condition. The upper limit of Q1 = 0 corresponds to 

the pressure p1 at which the water-table elevation is below the 

extinction depth. The lower limit of Q2 represents the maximum rate 

of fluid loss due to ET, which is achieved when the water-table 

elevation is above a specified depth at or below land surface. 

Parameter values are summarized in table 1. (An alternative 

interpretation is presented in section 1.4.7, “Evapotranspiration: 

Alternative Formulation for Unsaturated Flow.”)
10 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

top surface of the model, allowing the ET flux at each ET node to depend on simulated pressure at the node, which is indicative of saturation at the node. In that case, each ET node is assigned its own piecewise-linear Q versus p function (fig. 6). 

Negative pressures correspond to unsaturated conditions. Pressure p1 is the (negative) pressure below which roots are unable 

to extract groundwater. As pressure increases above p1, saturation increases, and therefore so does the rate of water loss due to 

ET. Pressure p2 is the (negative) pressure above which roots extract groundwater at the maximum possible rate. Note that in this 

alternative formulation, the loss of water due to ET can be distributed over multiple nodes at different depths, with a different 

rate of loss at each node.

1.4.8. Seepage Face

A vertical seepage face discharges to atmospheric pressure below the water table and acts as a no-flux boundary above the 

water table (Freeze and Cherry, 1979). Wherever the portion of the face above the water table is not vertical, the possibility of 

infiltration exists, and the face need not act as a no-flux boundary there.

A seepage-face boundary condition (fig. 7) can be implemented by setting

p Q 1 1 Q p inf i 2 2 Q Qnf < = 0 0 , , = = , and (15)

where Qinf is the infiltration rate [M/s] when the boundary 

node is unsaturated (p < 0). A limit is imposed on pressure at 

point 2. For a vertical seepage face, Qinf = 0 . The infiltration 

rate is commonly set to zero even when the seepage face is not 

vertical.

1.5. Linear Form of Generalized-Transport 

Conditions

Figure 8 shows schematic representations of three 

idealized concentration- or temperature-dependent boundary 

conditions: general concentration or temperature, specified 

source or sink of solute mass or energy, and specified concentration or temperature. In each case, the relation between the 

gain of solute or energy, QU [Ms/s or E/s], and concentration or 

temperature, U [Ms/ M or oC], is linear. By convention, gain of 

solute or energy by the model is positive and loss of solute or 

energy by the model is negative.

1.6. Two-Point Specification of GeneralizedTransport Conditions

A generalized-transport condition requires specification 

of a linear dependence of the solute or energy gain, QU, on 

concentration or temperature, U. Figure 9 illustrates how two 

Figure 7. Two-point representation of a seepage-face 

boundary condition. An upper limit of p2 = 0 is set on p, and Q1

and Q2 are both set to the unsaturated infiltration rate. For the 

vertical seepage face described by Freeze and Cherry (1979), 

the unsaturated infiltration rate is zero. Parameter values are 

summarized in table 1.

Figure 8. Three forms of 

the generalized-transport 

condition: A, general 

concentration or temperature; 

B, specified source or sink of 

solute mass or energy; and 

C, specified concentration or 

temperature.
Chapter 1. Generalized Boundary Conditions 11

user-specified points, (U1, QU

1 ) and (U2

, QU

2 ), define the line. 

Unlike generalized-flow conditions, generalized-transport conditions do not allow specification of bounds on QU or U.

For the sake of notational simplicity, the symbols U1, U2, 

QU

1 ,

 and QU

2 are used throughout Chapter 1 to represent the 

user-specified concentrations or temperatures and rates of gain 

or loss of solute or energy that define a generalized-transport 

boundary condition. In input datasets that define generalizedtransport conditions (“.inp” dataset 21B and “.bcs” dataset 

7B), U1, U2, QU

1 ,

 and QU

2 are represented by the more descriptive symbols UBG1, UBG2, QUBG

U

1, and QUBG

U

2, respectively.

1.7. Commonly Used Forms of the GeneralizedTransport Condition

The generalized-transport condition provides a framework for defining commonly used boundary conditions on 

solute or energy transport. Guidance on how to formulate 

a generalized-transport condition for two particular forms 

is provided below. Input parameter values are summarized 

in table 2. Although specified-concentration, specifiedtemperature, and concentration- or temperature-independent 

conditions are most easily handled using other types of 

SUTRA boundary conditions, instructions for formulating 

such conditions within the generalized-transport framework 

are given in the main text.

1.7.1. Heat Transfer Across a Thermal Boundary Layer

If a thermal boundary layer exists adjacent to the model boundary, the rate of heat transfer from an external point to a 

boundary node (for example, from the atmosphere to land surface) is proportional to the heat-transfer coefficient, the available 

area of the layer across which the heat transfer occurs, and the temperature difference across the layer:

Q h A T T U = − T ( ) ext (16)

where QU is the rate of heat transfer [E/s], hT is the heat-transfer coefficient [E/(L2·oC·s)], A is the area [L2] associated with the 

generalized-transport boundary condition node, Text is the temperature [oC] at the external point, and T is the temperature [oC] at 

Figure 9. For a generalized-transport condition, the gain of 

solute mass or energy, QU, is a linear function of concentration or 

temperature, U. By convention, positive values of QU denote gains 

and negative values denote losses. Parameter values for heat 

and solute transport across boundary layers are summarized in 

table 2.

Table 2. Summary of input parameter values for two commonly used forms of the generalized-transport condition. Each form is 

discussed in detail and symbols are defined in the main text.

Special cases1 U1 QU

1 U2 QU2 Notes

Thermal boundary 

layer

0 h AT Text Text 0 This formulation can be used if Text ≠ 0 . Otherwise, any two points that give the 

desired slope can be used to define the linear relation between QU and U. 

See equations (16) and (17).

Concentration 

boundary layer

0 h A M rCext Cext 0 This formulation can be used if Cext ≠ 0 . Otherwise, any two points that give the 

desired slope can be used to define the linear relation between QU and U. 

See equations (18) and (19). 1

Representation of concentration or temperature-independent sources and sinks and specified concentrations or temperatures using generalized-transport conditions is discussed in the main text. However, the simplest way to implement concentration or temperature-independent sources and sinks is to use “.inp” dataset 

18 and “.bcs” dataset 4, and the simplest way to implement specified concentrations or temperatures is to use “.inp” dataset 20 and “.bcs” dataset 6.
12 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

the node. A straightforward way to implement this type of boundary condition is to specify the two intercepts of the line represented by equation (16) by setting

T Q h AT T T Q U

T ext ext

U

1 1 2 2 = = 0 0 , , = = , and (17)

However, any two distinct points that lie on the desired line may be used. The user may specify one point and choose the second 

point such that the desired slope is obtained.

1.7.2. Solute Mass Transfer Across a Boundary Layer

Solute mass transfer across a boundary layer is analogous to the case of a thermal boundary layer described above. The 

rate of mass transfer from an external point to the generalized-transport boundary condition node is proportional to the masstransfer coefficient, the available area of the layer across which the solute transfer occurs, and the concentration difference across 

the layer:

Q h A C C U = − M r( ) ext (18)

where QU is the rate of solute mass transfer [Ms/(L2·s)], hM is the mass-transfer coefficient [L/s], A is the area [L2] associated 

with the boundary-condition node, ρ is the average fluid density [M/L3] within the boundary layer, and Cext and C are the mass 

fractions [Ms/M] of solute at the external point and the boundary-condition node, respectively. 

A straightforward way to implement this type of boundary condition is to specify the two intercepts of the line represented 

by equation (18) by setting

C ρ 1 = 0 QU

1 = hM A Cext , C2 = Cext , , and QU 0 2 = 0 (19)

However, any two distinct points that lie on the desired line may be used.

Note that although no explicit bounds are set on C, concentrations generated by the boundary condition (18) are inherently 

limited when there are no other sources of solute. Consider, for example, a block of salt dissolving into the groundwater. The 

maximum solute concentration is Csat, the concentration of saturated brine, and the rate of dissolution is given by equation (18) 

with Cext =Csat. Assuming C is initially less than Csat, the rate of dissolution is initially positive; salt dissolves into the groundwater. As C increases, the rate of dissolution decreases, approaching zero as C approaches Csat. Thus, in the absence of other 

sources of solute, C will never exceed Csat, so the concentration is inherently limited by the form of the boundary condition 

(equation 18).

1.7.3. Specified Source or Sink of Solute Mass or Energy

The simplest way to represent a concentration- or temperature-independent source or sink of solute mass or energy is 

through “.inp” dataset 18 and, for time-dependent conditions, the corresponding “.bcs” dataset 4. However, such a condition can 

also be formulated as a special case of a general concentration or temperature condition with Q Q U U

1 2 = . By convention, gains of 

solute or energy are positive and losses are negative.

1.7.4. Specified Concentration or Temperature

The simplest way to represent a specified-concentration or temperature boundary condition is through “.inp” dataset 20 

and, for time-dependent conditions, the corresponding “.bcs” dataset 6. However, the generalized-transport framework allows 

an alternative formulation for specified concentration or temperature that is analogous to the formulation of specified pressure 

depicted in figure 4A: U1 is set nearly equal to U2 so that the slope is very large. This is essentially how specified-concentration 

or temperature conditions were formulated in earlier versions of SUTRA. The disadvantage of this method is that the slope must 

be chosen by trial and error such that SUTRA computes the concentration or temperature and the corresponding solute or energy 

gain to approximately equal precision. For details, see section 3.5 and the input data listing for “.inp” dataset 5, both in the main 

SUTRA documentation (Voss and Provost, 2002).
Chapter 1. Generalized Boundary Conditions 13

1.8. Budget Output and Boundary-Condition Output Files

Fluid and solute or energy budget information is calculated and reported for generalized-flow and generalized-transport 

conditions, just as it is for the original four types of SUTRA boundary conditions (Voss and Provost, 2002, section 5.8). Budget 

information is written to the “.lst” output file.

In addition, two new files are available for output of information about generalized-flow and generalized-transport conditions: the “.bcopg” and “.bcoug” files, respectively. These files are similar in form to the “.bcof”, “.bcos”, “.bcop”, and “.bcou” 

files that report information about the original four types of SUTRA boundary conditions.

1.8.1. The “.bcopg” File for Generalized-Flow Nodes

The “.bcopg” file contains the following information on generalized-flow nodes:

• Node number.

• Where the condition was defined (“.inp” file, subroutine BCTIME, or “.bcs” file).

• “.bcs” identifier (if defined in a “.bcs” file).

• Resultant fluid source/sink rate (the rate at which fluid enters or leaves the system at that node to satisfy the generalizedflow condition).

• Concentration or temperature of the resultant fluid source/sink.

• Resultant source/sink of solute or energy (the rate at which solute mass or energy enters or leaves the system as a result 

of fluid inflow/outflow at that node).

• Computed pressure.

• Pressure at the first of two points used to define the relation between flow and pressure.

• Flow at the first of two points used to define the relation between flow and pressure.

• Pressure at the second of two points used to define the relation between flow and pressure.

• Flow at the second of two points used to define the relation between flow and pressure.

• Limit on pressure or flow (if any) set at the first point, and whether the limit is active.

• Limit on pressure or flow (if any) set at the second point, and whether the limit is active.

• Solute concentration or temperature associated with fluid outflow from that node, and whether it is specified directly 

(“DIR”) or relative to the value computed at the node (“REL”).

1.8.2. The “.bcoug” File for Generalized-Transport Nodes

The “.bcoug” file contains the following information on generalized-transport nodes:

• Node number.

• Where the condition was defined (“.inp” file, subroutine BCTIME, or “.bcs” file).

• “.bcs” identifier (if defined in a “.bcs” file).

• Resultant solute or energy source/sink rate (the rate at which solute or energy enters or leaves the system at that node to 

satisfy the generalized-transport condition).

• Computed concentration or temperature.

• Concentration or temperature at the first of two points used to define the relation between solute mass gain/loss and concentration or energy gain/loss and temperature, respectively.
14 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

• Solute mass or energy gain/loss at the first of two points used to define the relation between solute mass gain/loss and 

concentration or energy gain/loss and temperature, respectively.

• Concentration or temperature at the second of two points used to define the relation between solute mass gain/loss and 

concentration or energy gain/loss and temperature, respectively.

• Solute mass or energy gain/loss at the second of two points used to define the relation between solute mass gain/loss and 

concentration or energy gain/loss and temperature, respectively.

1.9. Assumptions and Limitations

The following is a summary of the assumptions and limitations inherent in SUTRA’s generalized boundary conditions:

• Like other types of SUTRA boundary conditions, generalized boundary conditions are applied at model nodes. (See p. 2.)

• For generalized-flow conditions, the dependence of the gain or loss of fluid mass on pressure is linear or piecewise-linear. 

(See section 1.1.)

• For generalized-transport conditions, the dependence of the gain or loss of solute mass or energy on concentration or 

temperature is linear. (See section 1.5.)

• The “river,” “drain,” and “ET” forms of the generalized-flow condition described in sections 1.4.5 and 1.4.6 represent 

hydrologic processes in a simplified way that depends on the assumption that the pressure profile in the vicinity of the 

boundary-condition node is approximately hydrostatic, allowing the water-table elevation to be estimated from the pressure computed at the node. The alternative formulation of ET presented in section 1.4.7 does not depend on this assumption. (See section 1.4.1.)

• The “general-pressure” form of the generalized-flow condition is based on the assumption of uniform and isotropic 

physical properties within the conductive layer. (See appendix 2.)

• The “heat transfer” and “solute mass transfer” forms of the generalized-transport condition are based on the assumption 

that a boundary layer exists, and that the rate of heat or solute-mass transfer across the boundary layer is proportional to 

the product of a heat- or mass-transfer coefficient and the concentration or temperature difference across the layer.

1.10. Example Problem: Seepage Face (Dam with Vertical Sides)

1.10.1. Physical Setup

This problem involves two-dimensional (2D), constant-density groundwater flow through a 10-meter-wide and 10-meterhigh dam with vertical sides under steady-state conditions (fig. 10). The water level is held even with the top of the dam on the 

upstream side, where hydrostatic conditions prevail, and water does not accumulate on the downstream side, where water exits 

through a seepage face (fig. 11).

The problem is nonlinear and may be solved by approaching the steady state gradually with a series of time steps. Initially, 

the dam is fully saturated with water, and at time zero, water begins to seep out from the downstream side of the dam and flow in 

through the upstream side. As the water-table elevation within the dam decreases, inflows and outflows become nearly equal and 

steady state is approached. A total simulation time of approximately 82. hr is selected, which is sufficient time for the problem to 

essentially reach steady state at the scale simulated. Near-steady-state simulation results are compared with predictions based on 

Dupuit’s formula for discharge through the dam (Muskat, 1937) and Polubarinova-Kochina’s (1962) semianalytical solution for 

the height of the seepage face.

1.10.2. Objectives

1. Demonstrate the use of a generalized-flow boundary condition to simulate a seepage face. 

2. Compare SUTRA results with analytical and semianalytical solutions.
Chapter 1. Generalized Boundary Conditions 15

1.10.3. Simulation Setup

The mesh (fig. 10) consists of 40 by 40 elements, each of size 0.25 m by 0.25 m, resulting in a total of 1,681 nodes (NN = 

1,681) and 1,600 elements (NE = 1,600). The mesh thickness, B, represents 1.0 m of the dam thickness, which extends infinitely 

in the direction perpendicular to the mesh. The initial time step is Δto = l. s, with time steps increasing by a factor of 1.5 on each 

subsequent step up to a maximum time step of Δt = 3600. s (1. hr). The simulation consists of 100 time steps. The total simulation time is approximately 2.9 x 105

 s (82. hr). Nonlinearity iterations are used to improve convergence to steady state. Convergence of nonlinearity iterations is controlled by the tolerance for pressure, RPMAX, which is set to 100. kg/(m⋅s2

) (equivalent to 

approximately 0.01 m of freshwater head). As few as one and as many as eight nonlinearity iterations are used on each time step, 

with later time steps generally requiring fewer iterations.

The left vertical boundary, which is held at hydrostatic pressure, acts as a source of fluid. The right vertical boundary, along 

which a seepage face forms, acts as a sink of fluid. No flow occurs across the top and bottom boundaries. Formulation of the 

boundary conditions is described in the section 1.10.5, “Boundary conditions.”

1.10.4. Parameters

α = 1.0 x 10–8 m⋅s2/kg β = 4.47 x 10–10 m⋅s2/kg

(The α and β values imply that Sop = 9.0447 x 10-9 m⋅s2/kg.)

k = 1.0 x 10–12 m2 ε = 0.1

ρ = 1000. kg/m3 μ = 1.0x10–3 kg/(m⋅s)

g = 9.81 m/s2

Figure 10. Boundary 

conditions and finite-element 

mesh for the seepage face 

example (dam with vertical 

sides; m, meter).
16 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

 for 0

0.3 0.7 1 5 x 10 for 0 -5 S -0.5

p

p p w =

≥

+ + ( ) 

  <









1

2



=  −



 



 − −  −



 

























k S S r

w w 0 3

0 7

1 1 0 3

0 7

0 5 2 0 5

.

.

.

.

. .



































2

 

Unsaturated properties functions

of Van Genuchten (1980)

with 5 x 10 m s /kg

and

-5 2 S a wres = = 0 3. , ⋅ ,

2. See eqns. (2.8) and (2.21) of

Voss and Provost (200

n =

2).

Here, α and β are compressibilities of the porous matrix and the fluid, respectively, Sop is specific pressure storativity, k is permeability, e is porosity, ρ is fluid density, μ is fluid viscosity, g is the acceleration of gravity, Sw is water saturation, Swres is residual 

water saturation, a and n are Van Genuchten model parameters, and kr is relative permeability. These physical properties are 

defined and discussed in more detail in Voss and Provost (2002).

1.10.5. Boundary Conditions

No flow occurs across the top and bottom boundaries. Hydrostatic pressure is specified at all nodes along the left vertical 

boundary. At the node in the upper-left corner, pressure is held at zero.

All nodes along the right vertical boundary are generalized-flow boundary-condition nodes set up such that a seepage face 

can form to allow water to exit, and no flow crosses the vertical boundary above the seepage face. This condition is set up as 

described in the section 1.4.8, “Seepage Face,” and figure 7. The two points in the two-point specification are (p1, Q1) = (-1., 0.) 

and (p2, Q2) = (0., 0.), and an upper bound of zero is set on pressure by setting a limit at p2. Thus, the seepage face is a zeropressure boundary, and above the seepage face is a no-flux boundary.

Figure 11. Near-steadystate SUTRA results for the 

seepage face example (dam 

with vertical sides). The thick 

line indicates the numerically 

simulated position of the water 

table (p = 0). The intersection 

of the thick line with the right 

vertical boundary indicates the 

numerically simulated top of 

the seepage face. The black 

circle indicates the top of the 

seepage face predicted by 

the semianalytical solution of 

Polubarinova-Kochina (1962). 

Dots indicate element centers 

from which flow-velocity 

vectors (thin line segments) 

emanate. Vector length 

is proportional to velocity 

magnitude. Vectors are plotted 

for every second element in 

each coordinate direction. (m, 

meter)
Chapter 1. Generalized Boundary Conditions 17

1.10.6. Initial Conditions

Initially, hydrostatic pressure is set at each node. At all nodes along the top boundary, initial pressures are set at zero, so the 

simulation begins fully saturated.

1.10.7. Results

Figure 11 shows the simulated, near-steady-state water table and flow field after approximately 82. hr of simulated time. 

Along the right vertical boundary, no flow crosses the boundary above the water table, and below the water table water discharges through the seepage face. The slight decrease in the slope of the water table as it nears the right vertical boundary is 

due to the fact that at that boundary the simulated water table must coincide exactly with a node at which p = 0 is specified; the 

water table cannot fall between two nodes at that boundary. The intersection of the water table with the right vertical boundary 

marks the top of the seepage face. The semianalytical solution of Polubarinova-Kochina (1962), which does not account for flow 

through the unsaturated zone, predicts a slightly lower seepage face height (fig. 11, black circle). The semianalytical seepage 

face height (approximately 3.7 m) was determined from figure 2 of Polubarinova-Kochina (1962) for parameter values l = 10 m, 

H = 10 m, and h = 0 m, where l, H, and h are the dam width, the height of the water on the upstream side, and the height of the 

water on the downstream side, respectively.

The Dupuit formula for discharge, Qdam [m3/s per meter of thickness], through a dam with vertical sides (Muskat, 1937),

Q k g

l

H h dam = − ( ) ρ

μ

2

2 2

2 (20)

is based on the Dupuit-Forcheimer assumptions of horizontal flow and a hydraulic gradient equal to the slope of the water table. 

Nevertheless, equation (20) is an exact expression for discharge through the dam (see, for example, Polubarinova-Kochina, 

1962) if it is assumed that no flow occurs through the unsaturated zone. For this example, the discharge predicted by equation 

(20) is 4.905x10-5 m3/s per meter of dam thickness. This is 6 percent less than the discharge simulated by SUTRA, which is 

5.228 x10-5 m3/s (5.228 x10-2 kg/s) per meter of dam thickness and includes flow through the unsaturated zone.
18 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

Chapter 2. Modified Implementation of Specified Pressures and Concentrations 

or Temperatures

In versions of SUTRA prior to Version 3.0, the implementation of specified-pressure and specified-concentration or temperature boundary conditions required the user to provide values for the factors GNUP and GNUU in dataset 5 of the “.inp” 

input file. Beginning with SUTRA Version 3.0, these types of boundary conditions are implemented without the use of GNUP or 

GNUU, and the latter have been removed from dataset 5. Modified input instructions are provided in appendix 3.

The SUTRA iterative and direct solvers handle storage of the coefficient matrix differently. Therefore, different numerical 

schemes for specified pressures are used with the two types of solvers, as described below. SUTRA automatically chooses the 

appropriate numerical scheme, depending on the types of solvers used.

2.1. Specified-Pressure Conditions: Scheme Used with Iterative Solvers

Specifying pressure pi

BC at boundary node i results in an associated gain of fluid mass, Qi

PBC , at that node. Addition of 

Qi

PBC as a source/sink term to the right-hand side of the discretized flow equation for node i gives

a pii i i a pj j b Q

j i

i i

PBC + = + ≠

∑ (21)

where aii is the diagonal element of the coefficient matrix in the absence of the specified-pressure condition, aij are the offdiagonal elements, bi

 is the right-hand side in the absence of the specified-pressure condition, pi is the computed pressure at 

node i, pj is the computed pressure at neighboring node j, and the summation is over the neighbors of node i, not including itself. 

Equation (21) above is equivalent to (3.39) in the main SUTRA documentation (Voss and Provost, 2002), except that the latter is 

expressed in terms of head to simplify the presentation.

In previous versions of SUTRA, Qi

PBC was expressed as

Q p p i

PBC

p i

BC = −  ( )i (22)

which, when substituted into equation (21), yielded the discretized flow equation

a p a p b p ii p i ij j

j i

i p i

BC ( ) + + = + ≠

  ∑ (23)

at node i. The symbol  p represents the factor GNUP, which was set large enough so that equation (23) effectively reduced to 

p p i i

BC ≈ ; the user chose GNUP such that the computed value of pi matched pi

BC to about the first half of the available doubleprecision digits (making the simulated pressure match the specified pressure for the first six or seven digits). The other half of 

the available precision was used by SUTRA to compute Qi

PBC using equation (22) after the discretized flow equation was solved.

In the new scheme used with iterative solvers in SUTRA 3.0, Qi

PBC is calculated from a rearrangement of equation (21) 

with p p i i

BC = :

Q a p a p b i

PBC

ii i

BC

ij j

j i

= + − i

≠

∑ (24)

As before, Qi

PBC is calculated in a postprocessing step after the discretized flow equation has been solved for pressure.

Three alternatives for incorporating specified-pressure conditions into the discretized flow equation are:

a. A specified pressure at node i could be enforced by simply replacing equation (21) with the equation p p i i

BC = . This 

straightforward approach would, however, introduce asymmetry into the matrix and overwrite the original values of aii, aij, 

and bi, which are needed to evaluate Qi

PBC using equation (24). Therefore, this approach would be incompatible with the 

CG solver, which requires a symmetric matrix, and additional storage would be required to save the original values of aii, 

aij, and bi.

b. To instead preserve matrix symmetry if it exists, and to minimize the additional storage required, one could use an 

extremely large value of GNUP in equation (23). This approach, however, is likely to result in loss of precision when 
Chapter 2. Modified Implementation of Specified Pressures and Concentrations or Temperatures 19

recovering the original values of aii and bi from aii +  p and b p i p i

BC +  , respectively, because it would typically involve 

adding and subtracting numbers of vastly different magnitudes. 

c. The approach used with iterative solvers in SUTRA 3.0 is to enforce specified-pressure conditions by multiplying the diagonal term in equation (21) by an extremely large number, ηp

, (set by SUTRA, with no user input required) and replacing 

the right-hand side with p ii i

BC a p to give.

  p ii i ij j

j i

p ii i

BC a p + = a p a p ≠

∑ (25)

For extremely large ηp, the summation term is negligible, and equation (25) is effectively equivalent to p p i i

BC = . The 

original values of aij are retained in the matrix, and the value of aii can be recovered with negligible loss of precision by dividing 

ηp aii by ηp. Thus, only bi in equation (21) need be saved in additional storage.

The scheme described above is incompatible with the direct solver, which overwrites the coefficient matrix during the solution process. Therefore, SUTRA uses this scheme only with the iterative matrix solvers. Because this scheme does not introduce 

asymmetry into the matrix, it is compatible with the CG solver, which requires a symmetric matrix.

2.2. Specified-Pressure Conditions: Scheme Used with the Direct Solver

Equation (21) is the discretized flow equation for node i (the node of interest at which a specified-pressure condition is 

being enforced), and there exists an analogous equation (with or without a specified-pressure source term, as appropriate) for 

every node in the model. For node i and nodes that are neighbors of node i, the discretized flow equations include the simulated 

pressure at node i, pi. 

The scheme for specified-pressure nodes that is used with the SUTRA direct solver involves substituting the specified pressure at node i, pi

BC , directly for the simulated pressure at node i, pi

, wherever the latter appears in the discretized flow equation 

for every node except node i. At node i, equation (21) is formulated under the temporary assumption that Qi

PBC = 0 . Upon solving the system of discretized flow equations just described, the resulting pressure at node i is denoted by pi

( ) 0 , and equation (21) 

can be written as

a p a p b ii i ij j

j i

i

( ) 0 + = ≠

∑ (26)

Combining equation (26) with equation (24) then gives

Q a p p i

PBC

ii i

BC = − ( ) i

( ) 0 (27)

which is used in a postprocessing step to compute Qi

PBC . Because the direct solver overwrites aii, the value of aii is saved in additional storage for use in equation (27). After the matrix problem is solved and Qi

PBC is evaluated, the value of pi

( ) 0 is no longer 

needed and is replaced in the p-solution vector by pi

BC, the desired value of pi

.

SUTRA uses this scheme only with the direct solver. Because this scheme introduces asymmetry into the matrix, it is 

incompatible with the CG solver, although it could, in principle, be used with the GMRES and ORTHOMIN solvers.

2.3. Specified-Concentration or Temperature Conditions

The scheme for handling specified-concentration or temperature conditions is entirely analogous to the direct-solver 

scheme for handling specified-pressure conditions described just above. For transport, this scheme is used with all of SUTRA’s 

matrix solvers, including the direct solver, except CG (which cannot be used for transport in any case).

The specified concentration or temperature at boundary node i, Ui

BC , is substituted directly for Ui

 , wherever the latter 

appears in the matrix equations, except in the equation for node i. At node i, the concentration or temperature is solved for under 

the temporary assumption that Qi

UBC = 0 , where QiUBC is the gain of solute mass or energy at node i resulting from the specifiedconcentration or temperature condition. If the resulting concentration or temperature at node i is denoted by Ui

( ) 0 , the corresponding equation is
20 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

   a U a U b ii i ij j

j i

i

( ) 0 + = ≠

∑ (28)

where “~” indicates that the matrix and right-hand side terms are for the transport equation. Combining equation (28) with the 

analog of equation (24),

Q a U a U b i

UBC

ii i

BC

ij j

j i

= + − i

≠

  ∑  (29)

then gives

Q a i U U UBC

ii i

BC = − ( ) i  ( ) 0 (30)

which is used in a postprocessing step to compute Qi

UBC . The value of aii is saved for use in equation (30) when the direct solver 

is used. After the matrix problem is solved and Qi

UBC

 is evaluated, the value of Ui

( ) 0 is no longer needed and is replaced in the 

U-solution vector by Ui

BC , the desired value of Ui

.
Chapter 3. Lake Capability 21

Chapter 3. Lake Capability

The lake capability introduced in SUTRA Version 3.0 simulates the evolution of lakes on the surface of a three-dimensional 

(3D) model as the lakes exchange fluid mass and solute mass or thermal energy with groundwater, with each other, and with 

external sources or sinks. Lakes influence the groundwater simulation by interacting with user-specified flow and transport 

boundary conditions on the top surface of the model. Lakes can increase boundary-condition pressures, intercept or supply 

recharge, and receive groundwater discharge and runoff. As lakes expand and contract during the simulation, the presence or 

absence of lake water can control whether a given boundary condition is applied or not applied, giving the user considerable 

control over how fluid and solute or energy enters or exits the model in the presence or absence of lakes. Exchanges of water 

between lakes, with groundwater, and with external sources and sinks can represent a variety of physical processes that include 

spillover from an adjacent lake, groundwater base flow, direct precipitation into lakes, evaporation of lake water, and surfacewater runoff into lakes.

Given the topography of the top surface of the model, or, optionally, user-specified lake-bottom elevations, SUTRA identifies the lakes that can form by “ponding” of water on the top surface. Increases and decreases in lake stage can cause lakes to 

coalesce and divide, leading to a hierarchy of possible “parent” and “child” lakes that is automatically determined by SUTRA

at the start of the simulation. On each time step of a transient simulation, SUTRA solves the groundwater flow and transport 

equations and updates lake stages and lake-water concentrations or temperatures. The calculations, which are described in more 

detail below, conserve water mass and solute mass or thermal energy within each lake, the groundwater system, and the SUTRA 

model as a whole. 

3.1. Preparatory Calculations

At the start of a simulation that involves lakes, SUTRA performs a number of preparatory calculations that pertain to lakes. 

These calculations include identification of potential lakes and determination of stage and area as functions of water volume for 

each lake.

3.1.1. Identification of Potential Lakes

The surface topography on which lake water can “pond” is determined by the elevations of the nodes on the top surface of 

the SUTRA model, or, optionally, by user-specified lake-bottom elevations. The lake capability can be used with 3D SUTRA 

meshes that are REGULAR, BLOCKWISE, or LAYERED. Each of these three types of mesh consists of layers of nodes (Voss 

and Provost, 2002). LAYERED meshes need not be logically rectangular but must have the same (regular or irregular) twodimensional structure within each layer. SUTRA infers the location of the top layer from the node numbering, which is assumed 

to conform to the convention used in SutraGUI (Winston and Voss, 2004). The lake capability cannot presently be used with 3D 

IRREGULAR meshes, nor with 2D meshes of any type.

To illustrate how SUTRA determines the set of potential lakes, consider the example shown in figure 12. To simplify the 

explanation, the topography in this example varies along only one grid direction, effectively reducing the problem to two dimensions, and the node numbering includes only surface nodes.

SUTRA identifies potential lakes by starting with the surface node that has the highest elevation and working down through 

surface nodes in order of decreasing elevation, finding the topographic ridges, or “sills,” that separate lakes. When the lake stage 

drops below a sill, a lake splits into two lakes.

In this example, the search begins with the water stage set at 1.0 (the highest elevation, which is found at nodes 1 and 31). 

For this lake stage, all the lake water is coalesced in one lake, which is designated as lake 1. When the stage drops below an 

elevation of 0.9 (node 6), lake 1 divides into two lakes: lake 2 on the right (nodes 7–30) and lake 3 on the left (nodes 2–5).

The next division occurs when the stage drops below 0.8 (node 26). Then lake 2 divides into lake 4 (nodes 8–25) and lake 5 

(nodes 27–30). Similar divisions occur as the stage decreases through elevations of 0.7 (node 12) and 0.6 (node 20).

The nine potential lakes identified in figure 12 form a hierarchy of “parents” and “children” that is represented by an 

inverted tree in figure 13. As the stage decreases, parent lakes split into children. Children of the same parent lake are called 

“siblings.” As the stage increases, sibling lakes coalesce into their parents. SUTRA assumes that the hierarchy of lakes is represented by a binary tree; each parent has exactly two children. Implications of this assumption are discussed in section 3.4.2, 

“Lake Hierarchy as a Binary Tree.” The algorithm used to identify potential lakes is discussed in more detail in section 3.4.1, 

“Identifying Potential Lakes Using the Hoshen-Kopelman Algorithm.” The “basin” of a lake is defined as the collection of all 

surface nodes that belong to that lake or any of its descendants in the binary tree.

Water balances are maintained on “active” lakes. A lake is considered active if it is not coalesced with its sibling or 

split into its children. If a lake is coalesced with its sibling, the lake and its sibling are both full, and a single water balance is 
22 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

maintained on the resulting parent lake. Thus, the original lake and its sibling become inactive, and the parent lake becomes 

active. If a lake is split into its two children, a separate water balance is maintained on each child lake. Thus, the parent lake 

becomes inactive and the child lakes become active.

3.1.2. Lake Stage and Area as Functions of Water Volume

Each surface node is overlain by a “lake cell” that can hold lake water. The volume of water in a lake is the sum of the 

volumes of water in all the lake cells within that lake. For the purpose of volume calculations, the lake-bottom elevation is 

assumed to be uniform within each lake cell, as shown in figure 14.

Figure 13. Hierarchy of nine lakes identified based on the 

topography in figure 12. Moving down the tree as stage 

decreases, lake 1 divides into lakes 2 and 3, lake 2 divides into 

lakes 4 and 5, and so forth. Moving up the tree as stage increases, 

lakes 8 and 9 coalesce into lake 6, lakes 6 and 7 coalesce into lake 

4, and so forth. Lake 6 is the “parent” of lakes 8 and 9; lakes 8 and 

9 are “siblings,” as they are both “children” of lake 6.

Figure 14. Lake cells (outlined in red and filled with pink) 

associated with the nodes (black dots) in lake 3 identified in 

figure 12. For the purpose of volume calculations, the lake-bottom 

elevation is assumed to be uniform within each lake cell. Straightline segments connecting nodes are included only to guide the 

eye in viewing the surface topography and are not used in volume 

calculations.

Figure 12. Nine potential lakes identified based on a given surface topography. Although the lake capability is 

designed to work only with three-dimensional SUTRA models, this schematic representation is presented in two 

dimensions to clearly illustrate how lakes are identified. The node numbering includes only surface nodes (black 

dots). Blue and green horizontal lines show the stages at which the various lakes come into existence by splitting 

of larger lakes as the stage decreases. Two different colors are used to help distinguish adjacent (“sibling”) lakes.
Chapter 3. Lake Capability 23

The rate at which the lake stage increases as the volume of water in the lake increases is inversely proportional to the 

surface area of the lake. For the purpose of lake-stage calculations, the surface area of the lake is the sum of the areas associated with all the submerged lake cells in the lake. Whenever an additional lake cell is submerged, the lake area increases by the 

area of that lake cell. Thus, lake stage is a continuous, piecewise-linear function of volume, with a decrease in slope occurring 

each time an additional lake cell is submerged. Each lake has its own stage-volume function, which is computed automatically 

by SUTRA.

3.2. Solution Procedure for Groundwater and Lakes

When the lake capability is used in a transient-flow, transient-transport SUTRA simulation, the following sequence of 

calculations is performed on each time step or nonlinearity iteration:

• Solve the groundwater-flow equation (Voss and Provost, 2002, equations (2.24) and (4.112)), taking into account all 

sources and sinks of groundwater, including exchanges of water with lakes via groundwater-flow boundary conditions. 

Lake stages are from the previous time step or nonlinearity iteration.

• Solve the groundwater-transport equation (Voss and Provost, 2002, equations (2.52) and (4.145)), taking into account all 

sources and sinks of solute mass or thermal energy, including exchanges of solute mass or energy with lakes via groundwater-flow and transport boundary conditions. Lake concentrations or temperatures are from the previous time step or 

nonlinearity iteration.

• Solve the set of water-mass and solute or energy balance equations for the lakes (described in section 3.2.1, “Water 

Balances on Lakes,” and section 3.2.2, “Solute Mass or Energy Balances on Lakes”), taking into account all sources and 

sinks of lake water and solute or energy, including exchanges with the groundwater-flow system. To conserve water mass 

and solute mass or energy in the simulation, exchanges are consistent with those used in solving the groundwater-flow 

and groundwater-transport equations. The result of solving the balance equations is an update of lake volumes, stages, 

and concentrations or temperatures.

When the lake capability is used in a steady-flow, steady-transport simulation, SUTRA uses the lake stages and concentrations or temperatures specified by the user or saved in the restart information from a previous simulation to formulate boundary 

conditions for the groundwater flow and transport equations. Lakes are not updated. This simulation mode is useful for computing initial conditions for transient simulations. The lake capability is not programmed for use in the steady-flow, transienttransport simulation mode in this version of SUTRA.

3.2.1. Water Balances on Lakes

After the groundwater flow and transport equations have been solved, exchanges of fluid between lakes and groundwater 

are computed at all boundary-condition nodes that interact with lakes. Exchanges of fluid with groundwater, other lakes, and 

external sources and sinks are then summed for each active lake (that is, each lake that is not coalesced with its sibling or split 

into its children). Rates of fluid exchange with groundwater, other lakes, and external sources and sinks are assumed to remain 

constant throughout each time step. The net inflows or outflows to or from lakes over a time step can cause splitting or coalescence of lakes.

3.2.1.1. Spillover, Coalescence, and Splitting of Lakes

If a lake receives a net inflow of water, its stage increases. If the lake fills completely, that is, if its stage rises to the elevation of the topographic sill that separates the lake from its sibling lake, any further net inflow spills over into the sibling lake (or 

an active descendant thereof, as described in section 3.4.3, “Determining the Fate of Spillover Using Steepest Descent”). Spillover occurs by a surface-water process that is not explicitly modeled and is assumed to occur instantaneously; excess water is 

simply reallocated from the lake to its sibling. If both the lake and its sibling fill completely, they coalesce into their parent lake.

If a lake experiences a net outflow of water, its stage decreases. If the stage of a lake decreases to the elevation of the topographic sill that separates its child lakes, any further net outflow causes the lake to split into its children. If the stage decreases 

below the elevation of the lowest surface node in the lake, the lake dries out completely. Complete drying and subsequent rewetting are described in section 3.4.4, “Drying and Rewetting.”
24 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

3.2.1.2. Equations Governing the Mass of Water in a Lake

During a given time step, the mass of water in each active lake varies linearly with time:

M M F t t l l l

net = + ( ) − 0

0 (31)

where

Ml = mass of water in lake l at time t [M], 

Ml

0

= mass of water in lake l at time t

0

 [M],

 t = current time [s],

 t0 = time at which the last breakpoint occurred [s], and

Fl

net

 = net rate of inflow of water to or outflow of water from lake l [M/s].

The net rate of inflow or outflow of water, Fl

net

, is the sum of the sources and sinks associated with spillover to or from other 

lakes, groundwater-flow boundary conditions (which can act as external sources and sinks when intercepted by a lake), and 

removal of lake water due to lake-area limitations (see section 3.3.1, “Limiting Where Lakes Can Form”):

Fl FFF F F net

k l

spill

k

NAL

l l

spill

l

bc in

l

bc out

l

lim = − + + + =

∑ , , ′

, ,

1

(32)

where

NAL = number of active lakes,

Fk l

spill

, = rate of spillover from lake k to lake l [M/s],

Fl l

spill

, ′ = rate of spillover from lake l to its sibling, lake l’ [M/s],

Fl

bc,in = rate of inflow to lake l from groundwater-flow boundary conditions [M/s],

Fl

bc,out = rate of outflow from lake l to groundwater-flow boundary conditions [M/s], and

Fl

lim = rate of outflow from lake l due to lake-area limitations [M/s].

Rates of inflow to lake l are positive and rates of outflow from lake l are negative by definition in equation (32), except rates of 

spillover. Spillover represents outflow from one lake and inflow into another. Therefore, rates of spillover are defined as positive, 

and Fl l

spill

, ′ is included in equation (32) with a negative sign to indicate an outflow from lake l. Note that a lake can receive spillover from its own sibling or from siblings of inactive lakes further up the hierarchy (as described in section 3.4.3, “Determining 

the Fate of Spillover Using Steepest Descent”).

3.2.1.3. Calculation of Water Volume from Water Mass in a Lake

The volume of water in each lake, which in turn determines the stage, is computed from the mass of water under the 

approximation that lake water has the density of water at the base concentration or temperature, RHOW0 (specified in dataset 9 

of the “.inp” input file). This greatly simplifies the lake calculations by making water volume a linear function of water mass and 

decoupling the water balances from the solute-mass or thermal-energy balances. In systems in which lake-water density departs 

significantly from RHOW0, the error introduced in the volumetric calculations may not be negligible. However, this approximation does not compromise the conservation of mass in each lake, the groundwater-flow system, or the model as a whole.

3.2.2. Solute Mass or Energy Balances on Lakes

Lakes are treated as well-mixed tanks; water entering a lake from various sources instantaneously mixes with the water 

already in the lake. The equation that governs the evolution of concentration or temperature in lake l under these conditions is

M u M u G t t l l l l l

net = + − 0 0

0 ( ) (33a)

and
Chapter 3. Lake Capability 25

G F u F u F u F u l

net

k

NAL

k l

spill

k l l

spill

l l l l

bc in

l

bc i = − = ′ 1 , , + + lim , , n

l

bc out

l l

bc in

l

bc out ∑ + + F u G G+ , , , (33b)

where

ul = concentration or thermal energy content in lake l [Ms

/M or E/M],

uk = concentration or thermal energy content in lake k [Ms

/M or E/M],

ul

bc,in = concentration or thermal energy content of inflows to lake l from groundwater-flow boundary conditions [Ms

/M 

or E/M],

Gl

bc,in = rate of inflow of solute mass or energy to lake l from groundwater-transport boundary conditions [Ms

/s or E/s], 

and

Gl

bc,out = rate of outflow of solute mass or energy from lake l to groundwater-transport boundary conditions [Ms

/s or 

E/s].

For solute-transport simulations, u represents solute concentration [Ms/M]. For energy-transport simulations, u represents 

thermal energy content [E/M], which is defined as u c = wT , where cw is the specific heat [E/(M∙°C)] of the water and T is its 

temperature [°C]. Rates of inflow to lake l are positive and rates of outflow from lake l are negative by definition in equation 

(33), except rates of spillover, as discussed in section 3.2.1.2, “Equations Governing the Mass of Water in a Lake.”

3.2.3. Effect of Lakes on Groundwater-Flow and Transport Boundary Conditions

Lakes influence the groundwater-flow and transport simulation by modifying the flows, pressures, and concentrations or 

temperatures associated with boundary-condition nodes on the top surface of the model. The way a lake interacts with a boundary-condition node depends on the type of boundary condition.

When lake water overlies a node at which a source or sink of fluid mass, solute mass, or energy is specified, the lake 

intercepts the source or sink; fluid, solute, or energy is added to or removed from the lake water instead of the groundwater. 

The details of this type of interaction are described in sections 3.2.3.1, “Interaction with Fluid Sources and Sinks,” and 3.2.3.5, 

“Interaction with Sources and Sinks of Solute Mass or Energy.”

When lake water overlies a node at which pressure, concentration, or temperature is specified, the lake water modifies 

the boundary pressure, concentration, or temperature that controls the flow of water, solute, or energy in or out of the groundwater system, and the lake water supplies or receives the resultant inflow or outflow. The details of this type of interaction are 

described in sections 3.2.3.2, “Interaction with Specified Pressures,” and 3.2.3.6, “Interaction with Specified Concentrations or 

Temperatures.”

Note that in the case of a specified source or sink, the lake exchanges fluid, solute, or energy directly with the source or sink 

at a rate specified by the user. By contrast, in the case of a specified pressure, concentration, or temperature, the lake exchanges 

fluid, solute, or energy with the groundwater at a rate computed by SUTRA.

Generalized boundary conditions (described earlier in this report) include specified-source/sink and specified-pressure/concentration/temperature conditions as special cases. Therefore, in specifying how each generalized boundary condition behaves 

when interacting with a lake, the user has two options. One option is to set the generalized boundary condition to behave like a 

specified source or sink by allowing the lake to intercept the source or sink of fluid, solute, or energy. This option can be used to 

simulate, for example, evaporation of pure water from solute-laden lake water at a user-specified rate. The other option is to set 

the generalized boundary condition to behave like a specified pressure, concentration or temperature by exchanging fluid, solute, 

or energy with the lake. This option can be used to simulate, for example, groundwater discharge across a conductive lake bed 

into a lake at a rate controlled by the pressure exerted by the overlying lake water. The two options are described in detail in 

section 3.2.3.3, “Interaction with Generalized-Flow Conditions.”

Boundary-condition nodes on the top surface of the model that are not overlain by lake water interact with the groundwater 

model in the same way as in previous versions of SUTRA (Voss and Provost, 2002). However, such nodes can also contribute 

water to lakes via runoff. The rate of runoff is related to the rate of groundwater recharge or discharge at the node via a userspecified ratio, as described in more detail in section 3.2.3.4, “Runoff to Lakes.”

The presence or absence of lake water can also optionally control whether a given boundary condition is applied or not 

applied by SUTRA. Only boundary-condition nodes at which lakes are allowed to form (see section 3.3.1) can be affected by 

the presence or absence of lake water. This feature, together with the interactions between lakes and active boundary-condition 

nodes described above, allows considerable flexibility in formulating lake-related boundary conditions of practical interest. In 

the discussion that follows, each node in question is assumed to be an active boundary-condition node on the top surface of the 

model, unless specified otherwise.
26 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

3.2.3.1. Interaction with Fluid Sources and Sinks

Associated with each fluid source or sink node (“.inp” dataset 17 and “.bcs” dataset 3) is a user-specified rate of inflow or 

outflow, QIN [M/s]. By definition, the value of QIN is positive for inflow to and negative for outflow from the model. Figure 15

shows how QIN enters or exits the model, depending on whether or not the node is submerged by a lake.

If the node is submerged by a lake (figs. 16A and 16B), QIN is added to or withdrawn from the lake. This allows direct addition of water to lakes, as by precipitation, or direct removal of water from lakes at the lake-water concentration or temperature, 

as by pumpage directly from the lake.

If the node is not submerged by a lake (figs. 16C and 16D), QIN is added to or withdrawn from the groundwater system. 

In addition, runoff at a rate of φ|QIN| is added to the lake in whose basin the node resides, where φ is the user-specified runoff 

expressed as a fraction of recharge or discharge. In the case of recharge (QIN>0), the magnitude of the total inflow to the groundwater-lake system is (1+ φ)| QIN|, and in the case of discharge (QIN<0), the magnitude of the total outflow from the groundwaterlake system is (1-φ)| QIN|. Runoff is described in more detail in section 3.2.3.4, “Runoff to Lakes.”

Water added to the groundwater flow system or lake from outside the model has the user-specified concentration or temperature assigned to the boundary node. Water withdrawn from a lake has the concentration or temperature of the lake water. 

Groundwater discharge and any associated runoff has the concentration or temperature computed at the boundary node.

3.2.3.2. Interaction with Specified Pressures

Associated with each specified-pressure node (“.inp” dataset 19 and “.bcs” dataset 5) is a user-specified pressure, pBC

[M/(L·s2)], and a resultant rate of fluid mass gain, QPBC [M/s]. By definition, the value of QPBC is positive for inflow to and negative for outflow from the model. Figure 16 shows how QPBC enters or exits the model, depending on whether or not the node is 

submerged by a lake.

If the node is submerged by a lake (figs. 17A and 17B), the user-specified pressure pBC is augmented by adding the pressure 

associated with the overburden of lake water. The overburden pressure is the product of the depth and density of the lake water 

and the acceleration of gravity. The resultant groundwater recharge or discharge, QPBC, is supplied or received by the lake.

If the node is not submerged by a lake (figs. 17C and 17D), the boundary pressure remains at its user-specified value, pBC. 

The resultant groundwater recharge or discharge, QPBC, is supplied or received by an external source or sink. In addition, runoff 

at a rate of φ|QPBC| is added to the lake in whose basin the node resides, where φ is the user-specified runoff expressed as a fraction of recharge or discharge. In the case of recharge (QPBC>0), the magnitude of the total inflow to the groundwater-lake system 

is (1+φ)|QPBC |, and in the case of discharge (QPBC<0), the magnitude of the total outflow from the groundwater-lake system is 

(1-φ)|QPBC |. Runoff is described in more detail in section 3.2.3.4, “Runoff to Lakes.”

Water added to the groundwater flow system or lake from outside the model has the user-specified concentration or 

temperature for inflow assigned to the boundary node. Groundwater recharge from the lake has the concentration or temperature of the lake water. Groundwater discharge and any associated runoff has the concentration or temperature computed at the 

boundary node.

Figure 15. Interaction of a lake (light blue) with a fluid source (inflow) or sink (outflow) at model nodes (black circles). 

Solid arrows represent the specified inflow or outflow, QIN, to or from the aquifer or lake at each node. Four cases 

are illustrated: A, inflow at a submerged node; B, outflow at a submerged node; C, inflow at a nonsubmerged node; 

and D, outflow at a nonsubmerged node. At submerged nodes A and B, QIN (solid arrow) is added to or withdrawn 

from the lake. At nonsubmerged nodes C and D, QIN (solid arrow) is added to or withdrawn from the aquifer. Also, the 

nonsubmerged nodes contribute runoff (outlined arrow) to the lake at a non-negative rate that is determined by the 

user-specified ratio, φ, of runoff to specified inflow or outflow: runoff = φ|QIN|.
Chapter 3. Lake Capability 27

3.2.3.3. Interaction with Generalized-Flow Conditions

The interaction of lakes with generalized-flow conditions is analogous to their interaction with either fluid sources and 

sinks or specified pressures. The type of interaction is chosen by the user in the “.lkbc” input file, which is summarized in section 3.2.3.8, "Application of Boundary Conditions: The “.lkbc” Input File." The input format for the “.lkbc” file is described in 

detail in appendix 3.

Associated with each generalized-flow node (“.inp” dataset 21A and “.bcs” dataset 7A) are user-specified pressures, pBG1

and pBG2 [M/(L·s2)], and rates of inflow or outflow, QPBG1 and QPBG2 [M/s], that define the piecewise-linear relationship between 

computed pressure and inflow or outflow of fluid at the node. (Input parameters pBG1, pBG2, QPBG1, and QPBG2 represent the same 

quantities as the symbols p1, p2, Q1, and Q2, respectively, which are used throughout Chapter 1 for the sake of notational simplicity.) For a given pressure p computed by SUTRA at the generalized-flow node, QPBG [M/s] is the resultant rate of inflow or 

outflow at the node. By definition, the value of QPBG is positive for inflow to and negative for outflow from the model. The type 

of interaction chosen by the user affects how QPBG is apportioned between the lake and the groundwater and whether pBG1 and 

pBG2 are modified by the presence of the lake, as described below.

3.2.3.3.1. Interaction Analogous to a Fluid Source or Sink

If the interaction of lakes with a generalized-flow condition is specified to be analogous to the interaction of lakes with 

a fluid source or sink condition, QPBG is apportioned between the lake (as runoff) and the groundwater recharge to or the discharge from the model. This is done in the same manner that QIN is apportioned between the lake and the groundwater recharge/

discharge for a fluid source or sink (fig. 16). Pressures pBG1 and pBG2 are unaffected by the lake. This facilitates simulation of 

processes such as evaporation from the lake.

Water added to the groundwater flow system or lake from outside the model has the user-specified concentration or temperature for inflow assigned to the boundary node. Groundwater discharge and any associated runoff, or water withdrawn from 

a lake, has the user-specified concentration or temperature for outflow assigned to the boundary node. For example, the outflow 

concentration can be set to zero to simulate evaporation of pure water from solute-laden lake water.

3.2.3.3.2. Interaction Analogous to a Specified Pressure

If the interaction of lakes with a generalized-flow condition is specified to be analogous to their interaction with a specifiedpressure condition, pressures pBG1 and pBG2 are increased by the pressure associated with the overburden of lake water. The resultant flow, QPBG, enters or exits the groundwater-flow system in the same manner as QPBC does for a specified pressure (fig. 17). 

This facilitates simulation of processes such as discharge of groundwater across a conductive lake bed and into a lake.

Water added to the groundwater flow system or lake from outside the model has the user-specified concentration or temperature for inflow assigned to the boundary node. Groundwater recharge from the lake has the concentration or temperature of 

the lake water. Groundwater discharge and any associated runoff has the user-specified concentration or temperature for outflow 

assigned to the boundary node.

Figure 16. Interaction of a lake (light blue) with specified pressures at model nodes (black circles). Solid arrows 

represent the resultant inflow or outflow, QPBC, to or from the aquifer at each node. Four cases are illustrated: A, 

inflow at a submerged node; B, outflow at a submerged node; C, inflow at a nonsubmerged node; and D, outflow at a 

nonsubmerged node. At submerged nodes A and B, QPBC (solid arrow) is supplied or received by the lake, and the lake 

modifies the user-specified boundary pressure, pBC, by adding the pressure pLAKE associated with the overburden of 

lake water. At nonsubmerged nodes C and D, QPBC (solid arrow) is controlled by the user-specified boundary pressure 

pBC. Also, the nonsubmerged nodes contribute runoff (outlined arrow) to the lake at a non-negative rate that is 

determined by the user-specified ratio, φ, of runoff to inflow or outflow to or from the aquifer: runoff = φ|QPBC|.
28 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

3.2.3.4. Runoff to Lakes

Although the SUTRA lake capability does not simulate the dynamics of stream and overland flow explicitly, the effects of 

such runoff to lakes in SUTRA are represented in an approximate way. At each active boundary node on the top surface of the 

model, runoff to lakes is related to the rate of groundwater recharge or discharge via a user-specified factor, φ, which is runoff 

expressed as a fraction of recharge or discharge. Thus, if the gain of fluid mass is Q (which is positive for recharge and negative for discharge), the rate of runoff is φ|Q|. Runoff is assigned to the lake in whose basin the boundary node resides. If the lake 

to which runoff is assigned is inactive (that is, if it is split into its children), runoff is assigned to a descendant of that lake that 

is identified using a discrete steepest-descent algorithm described in section 3.4.3, “Determining the Fate of Spillover Using 

Steepest Descent.” Thus, runoff follows the steepest downhill direction until it encounters either lake water or a completely dry 

lake bottom. Runoff is assumed to reach the lake instantaneously, which is usually a reasonable assumption for a groundwater 

model, as surface water typically redistributes more quickly than groundwater.

In the case of recharge (for example, fig. 16C), runoff can be conceptualized as originating from an external fluid source 

that is divided into recharge and runoff. In the case of discharge (for example, fig. 16D), a portion of the discharge can be conceptualized as being diverted to runoff, with the remainder going to an external sink. The user may specify a different fraction 

φ for recharge than for discharge at a given node. Although SUTRA does not restrict the allowable values of φ, a non-negative 

value is typically desired so that runoff is non-negative (into, not out of, the lake). If runoff is being diverted from groundwater discharge, a value of φ less than or equal to 1 is typically desired so that the diverted runoff does not exceed 100% of the 

discharge.

3.2.3.5. Interaction with Sources and Sinks of Solute Mass or Energy

Associated with each solute-mass or energy source or sink node (“.inp” dataset 18 and “.bcs” dataset 4) is a user-specified 

rate of inflow or outflow of solute mass or energy, QUINC [Ms/s or E/s]. The value of QUINC is positive for inflow to and negative for outflow from the model. Figure 17 shows how QUINC enters or exits the model, depending on whether or not the node is 

submerged by a lake.

If the node is submerged by a lake (figs. 18A and 18B), QUINC is added to or removed from the lake. This provides a mechanism for adding or removing solute mass or energy directly to or from lakes. If the node is not submerged by a lake (figs. 18C

and 18D), QUINC is added to or removed from the groundwater-flow system at that node.

3.2.3.6. Interaction with Specified Concentrations or Temperatures

Associated with each specified-concentration or temperature node (“.inp” dataset 20 and “.bcs” dataset 6) is a user-specified concentration or temperature, UBC [Ms

/M or °C] and a resultant gain of solute mass or energy QUBC

U [Ms

/s or E/s]. The value 

of QUBC

U is positive for inflow to and negative for outflow from the model. Figure 18 shows how QUBCU enters or exits the model, 

depending on whether or not the node is submerged by a lake.

If the node is submerged by a lake (figs. 19A and 19B), the boundary concentration or temperature is changed from UBC to 

the lake concentration or temperature, and the resultant flow into or out of the groundwater-flow system, QUBC

U , is supplied or 

received by the lake. If the node is not submerged by a lake (figs. 19C and 19D), the boundary concentration or temperature 

remains at UBC, and QUBC

U is added to or removed from the groundwater-flow system at that node.

Figure 17. Interaction of a lake (light blue) with solute or energy sources and sinks at model nodes (black circles). 

Arrows represent the specified source or sink of solute or energy, QUINC, to or from the aquifer or lake at each 

node. Four cases are illustrated: A, source at a submerged node; B, sink at a submerged node; C, source at a 

nonsubmerged node; and D, sink at a nonsubmerged node. At submerged nodes A and B, QUINC (arrow) is added to or 

withdrawn from the lake. At nonsubmerged nodes C and D, QUINC (arrow) is added to or removed from the aquifer.
Chapter 3. Lake Capability 29

3.2.3.7. Interaction with Generalized-Transport Conditions

The interaction of lakes with generalized-transport conditions is analogous to their interaction with either sources and sinks 

of solute or energy or specified concentrations or temperatures. The type of interaction is chosen by the user in the “.lkbc” input 

file, which is summarized in section 3.2.3.8, “Application of Boundary Conditions: The “.lkbc” Input File.” The input format for 

the “.lkbc” file is described in detail in appendix 3. 

Associated with each generalized-transport node (“.inp” dataset 21B and “.bcs” dataset 7B) are user-specified 

concentrations or temperatures, UBG1 and UBG2 [Ms/M or °C], and rates of gain or loss of solute mass or energy, QUBG

U

1 and 

QUBG

U

2 , that define the linear relationship between computed concentration or temperature and the rate of gain or loss of solute 

or energy at the node. (Input parameters UBG1, UBG2, QUBG

U

1 , and QUBG

U

2 represent the same quantities as the symbols U1

, U2, 

QU

1 , and QU

2 , respectively, which are used throughout Chapter 1 for the sake of notational simplicity.) For a given concentration or temperature U computed by SUTRA at the generalized-flow node, QUBG

U [Ms

/s or E/s] is the resultant rate of gain or 

loss of solute mass or energy. The value of QUBG

U is positive for inflow to and negative for outflow from the model. The type of 

interaction chosen by the user affects how QUBG

U is apportioned between the lake and the groundwater and whether UBG1, UBG2, 

QUBG

U

1, and QUBG

U

2 are modified by the presence of the lake, as described below.

3.2.3.7.1. Interaction Analogous to a Source or Sink of Solute or Energy

If the interaction of lakes with a generalized-transport condition is specified to be analogous to their interaction with a 

source or sink of solute or energy, QUBG

U enters or exits the model in the same manner that QUINC does for a source or sink of 

solute or energy (fig. 18). Concentrations or temperatures UBG1 and UBG2 are unaffected by the lake. This allows simulation of 

heat transfer between a lake and the atmosphere across a thermal boundary layer, for example, with the heat transfer being 

driven by the temperature difference between the lake water and the atmosphere.

3.2.3.7.2. Interaction Analogous to a Specified Concentration or Temperature

If the interaction of lakes with a generalized-transport condition is specified to be analogous to their interaction with a 

specified concentration or temperature, UBG1 and UBG2 depend on whether or not lake water is present. The resultant flow, QUBG

U , 

enters or exits the groundwater-flow system in the same manner that QUBC

U does for a specified concentration or temperature 

(fig. 19).

If the node is submerged by a lake, SUTRA automatically sets the point (UBG1, QUBG

U

1 ) to (ULAKE, 0), and adjusts the point 

(UBG2, QUBG

U

2 ) such that the original slope of the QU -versus-U line (fig. 9) is preserved. Thus, the resultant gain of solute or 

energy by the groundwater-flow system, QUBG

U , which is supplied or received by the lake, becomes zero when U = ULAKE. This 

allows simulation of heat transfer between a lake and groundwater across a thermally conductive lake bed and (or) a thermal 

boundary layer, for example, with the heat transfer being driven by the temperature difference between the lake water and the 

groundwater at the boundary node.

3.2.3.8. Application of Boundary Conditions: the “.lkbc” Input File

The presence or absence of lake water can optionally control whether a given boundary condition is applied or not applied. 

Consider, for example, a model in which the recharge rate due to precipitation is specified across the top of the model, but the 

recharge rate is different for precipitation that falls on a lake than it is for precipitation that falls on the land surface. One set of 

fluid-source boundary conditions specified at all the top-surface nodes provides the recharge to the land surface. This set is not 

Figure 18. Interaction of a lake (light blue) with specified concentrations or temperatures at model nodes (black 

circles). Arrows represent the resultant source or sink of solute or energy, QUBC

U , to or from the aquifer at each 

node. Four cases are illustrated: A, source at a submerged node; B, sink at a submerged node; C, source at a 

nonsubmerged node; and D, sink at a nonsubmerged node. At submerged nodes A and B, QUBC

U (arrow) is supplied 

or received by the lake and is controlled by the lake concentration or temperature, ULAKE. At nonsubmerged nodes C

and D, QUBC

U (arrow) is controlled by the user-specified boundary concentration or temperature, UBC.
30 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

applied when lake water is present and is applied when lake water is absent; it is applied only when the land surface is exposed. 

A second set of fluid-source boundary conditions specified at all the top-surface nodes provides the recharge to lakes. This 

second set is applied when lake water is present and is not applied when lake water is absent; it is applied only when the land 

surface is inundated. Thus, at any given fluid-source node at any given time, the recharge rate is either the value for lakes or the 

value for the land surface, depending on whether or not lake water is present.

SUTRA boundary-condition nodes are grouped by type and given default parameter values in “.inp” datasets 17–21. Timevarying parameter values that override the default values are specified in “.bcs” datasets 2–7, which can be distributed over 

multiple “.bcs” files. Specification of the way lakes control the application of boundary conditions follows a parallel structure in 

the “.lkbc” input file. For each of the six types of boundary condition, default specifications that apply to all boundary conditions of that type are listed first. These are followed by optional specifications for individual “.bcs” files that override the default 

specifications. Thus, by organizing boundary conditions in multiple “.bcs” files, one can effectively construct boundary conditions that change type and (or) parameter values depending on whether or not lake water is present.

Specifications for each of the six types of boundary condition include how the presence or absence of lake water affects the 

application status of boundary conditions of that type. For the generalized-flow and generalized-transport types, the user also 

specifies whether lakes interact with that type of boundary condition analogously to a source or sink, or analogously to a specified pressure, concentration, or temperature. The input format for the “.lkbc” file is described in detail in appendix 3.

3.3. Special Features of the Lake Capability

This section describes features that give the user additional control over the formation of lakes.

3.3.1. Limiting Where Lakes Can Form

In some applications, it may be desirable to simulate lakes over only a portion or portions of the surface of the model. 

To allow this, the subset of top-surface nodes at which lakes may form are listed in the optional “.lkar” input file. Input data 

formats for the “.lkar” file are listed in appendix 3. If the optional “.lkar” file is omitted, all nodes on the top surface of the 

model are eligible to form lakes. If the optional “.lkar” file is specified, then if the lowest-elevation node in the basin of a lake is 

disqualified from lake formation, lake water cannot accumulate in that basin, so that entire basin is disqualified from lake formation. Top-surface nodes at which lakes may form must form one contiguous region; they cannot be separated into two or more 

separate regions in which lakes may form. The contiguous region may contain “holes” in which lakes are not allowed to form, 

but there must be only one contiguous region so that potential lakes are all contained within a single hierarchy.

When a lake expands such that it reaches the outer boundary of the model or a node at which lakes are not allowed to 

form, the lake volume (stage) is not allowed to increase, and any excess inflow is assumed to “spill over” out of the model. 

The amount of spillover out of the model is reported in the budget information, as summarized in section 3.5, “Lake-Related 

Output Files.”

3.3.2. Lake-Bottom Elevations

The lake capability identifies potential lakes and their hierarchy automatically, as described in section 3.1.1, “Identification 

of Potential Lakes.” If the surface topography includes numerous small undulations, a very large number of potential lakes can 

result. To keep the number of lakes and their hierarchy manageable, lake-bottom elevations can be specified by the user in the 

optional “.lkar” file. This feature allows the user to supply lake-bottom elevations that have been smoothed in a preprocessing 

step, if desired, to eliminate small undulations. Optionally, the user can choose to have SUTRA use the top-surface node elevations of the groundwater model as the lake-bottom elevations by default.

3.4. Some Algorithmic Details

This section describes in greater detail some of the procedures used by the lake capability and their ramifications for 

SUTRA models involving lakes.

3.4.1. Identifying Potential Lakes Using the Hoshen-Kopelman Algorithm

The SUTRA lake capability identifies potential lakes by starting with the highest elevation in the grid and working its way 

down, finding the “sills” at which lakes divide. On a 2D elevation grid that represents the top surface of a 3D model, the lakes 
Chapter 3. Lake Capability 31

that form at a given water stage correspond to contiguous “clusters” of submerged lake cells on the “lattice” formed by the grid. 

When the stage is lowered and a lake splits, the number of clusters increases.

Consider, for example, the regular, rectangular grid shown in figure 19. The nodes in the grid vary in elevation (arbitrary 

units), and numerals indicate the cluster (lake) number with which each node is identified. In figure 20A, the water stage is at 

1.0—the elevation of the highest node in the grid. Thus, all nodes are submerged, and all belong to one large cluster (cluster 1). 

When the stage is lowered below 0.7 (the highest sill elevation), cluster 1 divides into two distinct clusters, 2 and 3 (fig. 20B). As 

the water stage is lowered further, clusters continue to divide until finally all possible clusters (lakes) have been identified. The 

final array (fig. 20C) corresponds to a lake hierarchy that is identical to the one in figure 13.

The lake capability uses a generalization of the Hoshen-Kopelman algorithm (Hoshen and Kopelman, 1976; Al-Futaisi and 

Patzek, 2002) to identify clusters on regular or irregular grids. Thus, the lake capability can be used with REGULAR, BLOCKWISE, or LAYERED SUTRA meshes. LAYERED meshes need not be logically rectangular but must have the same (regular or 

irregular) two-dimensional structure within each layer (Voss and Provost, 2002).

3.4.2. Lake Hierarchy as a Binary Tree

The lake capability assumes that when a lake divides, it splits into exactly two children. In other words, as the lake stage 

falls, it cannot encounter two sills within the same lake simultaneously; it must encounter one of them before the other. This 

poses no dilemma for division of lakes: one division can be followed immediately by a second, effectively splitting one lake 

into three. If the two sills are distinct but close in elevation, it generally does not make much practical difference which one is 

encountered first.

When the lake stage is rising, however, which sill is encountered first can have a significant effect on the exchange of water 

between lakes. Consider the situation shown in figure 20, in which the sill on the right is slightly lower than the sill on the left 

and is therefore encountered first by the rising stage in the middle lake.

When the stage rises to the elevation of the right-hand sill, the middle lake begins to spill into the right-hand lake, and it 

continues to do so until the right-hand lake is filled. The middle and right-hand lakes then coalesce and rise together. When their 

stage rises to the elevation of the left-hand sill, they begin to spill into the left-hand lake. If, instead, the left-hand sill had been 

slightly lower than the right-hand sill, the sequence would have been reversed: the left-hand lake would have filled first, then the 

right-hand lake. Even if the elevations of the two sills had been equal, the lake capability would have encountered one of them 

first in its search algorithm, and the corresponding lake would have been designated to receive the spillover.

In a real system, the lake stage necessarily rises above the sill elevation when one lake is spilling into another. All else 

being equal, the greater the rate of spillover and the narrower the channel between the lakes, the higher the stage rises above the 

sill. In fact, it is possible for the lake stage to be above two or more sills simultaneously, and thus for a lake to spill into two or 

more lakes at the same time (fig. 21). One of the sills may begin spilling over first, but that does not necessarily preclude others 

from spilling over as well, as the stage continues to rise.

This points out a limitation of the binary lake hierarchy used in the lake capability. Of the three physically possible 

scenarios—middle lake spills into right lake first, middle lake spills into left lake first, and middle lake spills into right and left 

Figure 19. The lake capability identifies potential lakes by tracking the formation of distinct clusters of contiguous, submerged 

nodes on the top surface of the model as the water stage varies from the highest to the lowest elevation on the surface. This 

example shows a regular, rectangular grid of surface nodes. Numerals indicate the cluster number with which each node is 

identified, and colors aid in distinguishing the clusters: A, With the water stage at the maximum elevation, all nodes belong to one 

large cluster (cluster 1); B, At a lower water stage, only the colored nodes are submerged, and cluster 1 splits into two distinct 

clusters (2 and 3); and C, When the water stage reaches the surface node that has the lowest elevation, all possible clusters 

(lakes) have been identified.
32 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

lakes simultaneously—only the first two can be simulated by the lake capability; the third cannot. In the present lake capability 

paradigm, each lake has but one sibling, and spillover is routed only to that one sibling (or an active descendant thereof). Users 

should remain mindful of this limitation when formulating their models and interpreting results. If a lake has two sills of equal 

elevation, the sill chosen by SUTRA as the sill that separates the lake from its sibling might or might not be the sill that the user 

would choose to begin spilling over first. In cases in which this limitation has a significant effect on the outcome of the simulation, small adjustments to the surface topography by the user may help to force a model to choose the most important exchanges 

of water between lakes.

3.4.3. Determining the Fate of Spillover Using Steepest Descent

When a lake fills to capacity, any additional net inflow of water spills over into the lake’s sibling. If the sibling lake is not 

active at the time (because it is split into its children), the spillover is rerouted to one of the sibling’s children. This rerouting 

continues down the lake hierarchy until an active descendant of the sibling lake is found.

Because the details of the surface-water flow that distributes the spillover are not modeled, a simple search algorithm is 

used to decide how spillover is routed down the lake hierarchy. The underlying assumption is that the surface flow follows a 

path of steepest descent from node to node.

The search begins at the node that is the sill between the spilling lake and its sibling. It then moves to the neighboring node 

(on the sibling’s side of the sill) that provides the steepest drop in elevation. The search continues in a like manner, moving from 

node to neighboring node according to the steepest-descent criterion, until a childless lake (a lake at the bottom of the hierarchy) 

is reached. The path traversed by the search determines the sequence of lakes to which spillover from the original lake is routed 

until an active descendant lake is encountered.

A topographic depression can be optionally eliminated from the lake hierarchy because it is too small to warrant treatment as an independent lake, as described in section 3.3.2, “Lake-Bottom Elevations.” Any spillover that is routed into such 

a depression by the steepest-descent algorithm is rerouted by treating the depression (for this purpose only) as a lake that is 

full and spilling over its sill. Spillover that enters the depression immediately exits at the sill and resumes following a path of 

steepest descent.

Figure 20. The middle lake spills into the right-hand lake first because the right-hand sill is lower than the left-hand 

sill. The spilling process (overland flow) is not modeled explicitly, and the stage in the middle lake remains fixed at the 

level of the right-hand sill until the right-hand lake is filled. Arrows indicate inflow and spillover of water.

Figure 21. In reality, a sufficiently strong source to the middle lake can cause it to spill over into the two side lakes 

simultaneously. This cannot be simulated by the SUTRA lake capability. Arrows indicate inflow and spillover of water.
Chapter 3. Lake Capability 33

3.4.4. Drying and Rewetting

The groundwater-flow equation is solved and exchanges of water between groundwater and lakes are computed before

mass balance calculations are performed on the lakes. Thus, the net outflow (total of all outflows from less the total of all inflows 

to a lake) of water from a lake on a given time step can exceed the amount of water in the lake. In that case, the water in the lake 

is “overdrawn” on that time step, and the overdraft is reported as a water mass balance error. Similarly, the groundwater-transport equation is solved and exchanges of solute mass or energy between groundwater and lakes are computed before solute-mass 

or energy balance calculations are performed on the lakes. If the net outflow of solute mass from a lake on a given time step 

exceeds the solute mass in the lake, the solute in the lake is “overdrawn” on that time step, and the overdraft is reported as a 

solute mass balance error. Thermal energy is not subject to overdraft, and no limit is imposed on the lake-water temperature.

A lake from which all water has been removed is said to be dry. A dry lake rewets if there is net inflow into the lowest node 

in the lake. Note that inflows can include runoff of external sources and groundwater discharge, as illustrated, for example, in 

figure 15.

3.5. Lake-Related Output Files

The lake number for each surface node is output in the “.lkn” file. The binary tree that represents the lake hierarchy is 

written to the “.lkh” file.

Fluid and solute or energy budget information is reported for lake-related nodes in the SUTRA “.lst” output file. It follows 

a format similar to that used for other budget items in the “.lst” file.

Budget information for individual lakes is reported in the “.lkbu” lake-budget output file:

• Status (coalesced, split, full, and so forth).

• Stage and its rate of change with time.

• Water volume and mass and their rates of change with time.

• Solute mass and concentration and their rates of change with time (solute transport only).

• Thermal energy content and temperature and their rates of change with time (energy transport only).

• Rates of water mass and solute mass or energy loss to lake-area limits.

• Rates of water mass gain, loss, and net exchange with groundwater, external sources and sinks, runoff, and other lakes.

• Rates of solute mass gain, loss, and net exchange with groundwater, external sources and sinks, runoff, and other lakes 

(solute transport only).

• Rates of thermal energy gain, loss, and net exchange with groundwater, external sources and sinks, runoff, and other 

lakes (energy transport only).

Lake stage and depth information is output to the “.lkst” file. This file has a format similar to that of the SUTRA “.nod” file. For 

inactive lakes and top-surface nodes at which lake water is absent, a user-specified marker value (RNOLK in dataset 1 of the 

“.lkin” input file) is printed. The lake capability is compatible with the SUTRA restart capability. Lake-related restart information is written to the “.lkrs” file.

3.6. Assumptions and Limitations

The following is a summary of the assumptions and limitations inherent in SUTRA’s lake capability:

• Lakes exchange fluid mass and solute mass or energy with groundwater, with each other, and with external sources or 

sinks. Lakes influence the groundwater simulation by interacting with user-specified flow and transport boundary conditions on the top surface of the model. (See p. 25.)

• Lakes form on top of the surface topography defined by the top layer of nodes in the SUTRA model, or, optionally, by 

user-specified lake-bottom elevations. (See section 3.1.1.)

• Because the lake capability requires an identifiable “top layer” of nodes, its use is restricted to 3D REGULAR, BLOCKWISE, and LAYERED meshes. LAYERED meshes need not be logically rectangular but must have the same (regular or 
34 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

irregular) two-dimensional structure within each layer (Voss and Provost, 2002). The lake capability cannot presently be 

used with 3D IRREGULAR meshes, nor with 2D meshes of any type. (See section 3.1.1 .)

• SUTRA automatically identifies the lakes that can form on the model surface and arranges them in a numbered hierarchy 

that represents the splitting and coalescence of lakes as the surface-water stage changes. Because the lake numbering 

is generally not known by the user a priori, formulation of lake specifications, such as initial lake stages, can involve 

trial and error over multiple model runs. An option to associate lake specifications with surface-node numbers instead of 

lake numbers is available to facilitate the process in the early stages of modeling. (See section 3.1.1 and dataset 3 of the 

“.lkin” input file in appendix 3.)

• Each surface node is overlain by a “lake cell” that can hold lake water. The total volume of water in a lake is the sum of 

the volumes of water held in the lake cells associated with that lake. For the purpose of volume calculations, the lakebottom elevation is assumed to be uniform within each lake cell. (See section 3.1.2.)

• The rate at which the lake stage increases as the volume of water in the lake increases is inversely proportional to the 

surface area of the lake. For the purpose of lake-stage calculations, the surface area of the lake is the sum of the areas 

associated with all the submerged lake cells in the lake. Thus, lake stage is a continuous, piecewise-linear function of 

volume. (See section 3.1.2.)

• The groundwater-flow equation includes exchanges of water with lakes via groundwater-flow boundary conditions. 

Exchanges are based on lake stages from the previous time step or nonlinearity iteration.

• The groundwater-transport equation includes exchanges of solute mass or energy with lakes via groundwater-flow and 

transport boundary conditions. Exchanges are based on lake concentrations or temperatures from the previous time step 

or nonlinearity iteration.

• The water-mass balance equations for the lakes include exchanges of water with the groundwater-flow system. To 

conserve water mass in the simulation, exchanges are consistent with those used in solving the groundwater-flow 

equation provided there is sufficient water in each lake to satisfy the computed exchanges. If water is overdrawn from a 

lake on a given time step, the overdraft is reported as a lake water mass-balance error.

• The solute-mass or energy balance equations for the lakes include exchanges of solute mass or energy with the groundwater-flow system. To conserve solute mass or energy in the simulation, exchanges are consistent with those used in 

solving the groundwater-transport equation provided there is sufficient solute mass in each lake to satisfy the computed 

exchanges. If solute mass is overdrawn from a lake on a given time step, the overdraft is reported as a lake solute massbalance error.

• Rates of fluid exchange with groundwater and external sources and sinks are assumed to remain constant throughout 

each time step.

• The volume of water in each lake, which in turn determines the stage, is computed from the mass of water under the 

approximation of constant density. This introduces error into the volumetric calculations but does not compromise the 

conservation of mass. (See section 3.2.1.3.)

• Lakes act as well-mixed tanks. There is no stratification of concentration or temperature within a lake. (See section 3.2.2.)

• The user can specify the set of surface nodes at which lakes may form. This set of nodes can comprise only a portion of 

the top surface of the model but must form a single, contiguous region such that potential lakes are all contained within a 

single hierarchy. (See section 3.3.1.)

• The lake hierarchy is a binary tree. When a lake divides, it splits into exactly two children. As a result, a lake that is filled 

to capacity cannot spill over into two or more lakes simultaneously. (See section 3.4.2.)

• Spillover from one lake to another occurs by a surface-water process that is not explicitly modeled and is assumed to 

occur instantaneously; excess water is simply reallocated from the lake to its sibling. (See section 3.2.1.1.)

• If the sibling of a lake that is spilling over is not active, the spillover is routed to an active lake lower in the hierarchy that 

is determined using a discrete steepest-decent algorithm. (See section 3.4.3.)
Chapter 3. Lake Capability 35

3.7. Example Problem: Lakes

3.7.1. Physical Setup

This problem involves interaction of 3D groundwater flow with lakes and transport of solute between groundwater and 

lakes. The model domain covers a 10-km-by-20-km rectangular area and is approximately 100 m deep, with a variable surface 

topography that corresponds to the top surface of the groundwater model, which includes three topographic depressions (fig. 22). 

Initially, no lake water is on the surface of the model. As topographically driven groundwater flow recharges at higher elevations and discharges at lower elevations, lakes form in the three topographic depressions. A source of solute at the model surface 

enters the groundwater and discharges into one of the lakes. As the lake stages rise, the lakes coalesce, and their waters mix. 

Eventually, one large lake forms. Subsequent withdrawal of lake water causes the lake stage to fall, and the large lake splits into 

the three original lakes. At the end of the simulation, the groundwater-lake flow system is near steady state.

3.7.2. Objectives

1. Demonstrate the interaction of groundwater with lakes simulated using the SUTRA lake capability. 

2. Demonstrate solute transport in groundwater and lakes simulated using the SUTRA lake capability.

3.7.3. Simulation Setup

The mesh (fig. 23) consists of 40 by 20 elements horizontally and 10 elements vertically, resulting in a total of 9,471 nodes 

(NN = 9,471) and 8,000 elements (NE = 8,000). Each element is 500. m by 500. m horizontally and between approximately 7. 

and 10. m vertically. The bottom of the model is at 0 m elevation, and the top-surface elevation, ztop [m], is given by

z r r r top = − + ( ) − + ( ) − + ( )100 15

1 1500

20

1 2000

30

1 3000 3

4

4

4

5

4 (34)

Figure 22. Finite-element 

mesh for the lake example. 

Green cubes denote specifiedpressure nodes that act as 

solute sources. Red cube 

denotes a fluid sink node that 

withdraws lake water. (50X 

vertical exaggeration; m, 

meter; km, kilometer)
36 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

where

r x y 3

2 2

= − ( ) 16000 + − ( ) 7000 (35)

r x y 4

2 2

= − ( ) 4000 + − ( ) 5000 (36)

and

r x y 5

2 2

= − ( ) 10000 + − ( ) 5000 (37)

Based on this surface topography, the SUTRA lake capability identified and reported the lake hierarchy shown in figure 23. 

Lakes 3, 4, and 5 correspond to the right, left, and middle topographic depressions in figure 23. Coalescence of lakes 4 and 5 

forms lake 2, and coalescence of lakes 2 and 3 forms lake 1.

Groundwater flow is fully saturated and of constant density. Concentration [dimensionless] is expressed on a relative scale 

of 0. to 1. The run is transient in both pressure and concentration. The time step size is Δt = 6.048 x 105 s (1. week). The simulation consists of 500 time steps. The total simulation time is 3.024 x 108

 s (500. weeks). The CG solver is used for groundwaterflow solutions and the GMRES solver is used for groundwater-transport solutions. Groundwater-flow solutions required no more 

than five solver iterations to converge to a tolerance of 1 x 10-8. Groundwater-transport solutions required no more than 13 solver 

iterations to converge to a tolerance of 1 x 10-10.

3.7.4. Parameters

a = 1.0 x 10–8 m⋅s2/kg b = 4.47 x 10–10 m⋅s2/kg

(The α and β values imply that Sop = 9.0447 x 10-9 m⋅s2/kg.)

k = 1.0 x 10–10 m2 e = 0.1

r = 1000. kg/m3 m = 1.0x10–3 kg/(m⋅s)

g = 9.81 m/s2 αL = 





Dm = 1.0 x 10-9 m2/s αT = 10. m

Here, a and b are compressibilities of the porous matrix and the fluid, respectively, Sop is specific pressure storativity, k is permeability, e is porosity, r is fluid density, m is fluid viscosity, g is the acceleration of gravity, Dm is apparent molecular diffusivity, 

αL is longitudinal dispersivity, and αT is transverse dispersivity. These physical properties are defined and discussed in more 

detail in Voss and Provost (2002).

3.7.5. Boundary Conditions

No flow occurs across the vertical and bottom boundaries. Atmospheric (zero) pressure is specified at all nodes on the top 

surface, and all nodes on the top surface are eligible to form lakes. The fraction of discharge that runs off from inactive lake 

nodes to lakes is set to zero. At time t = 30 weeks, the concentration associated with inflow to the model is changed from 0. to 1. 

(expressed in arbitrary units) at 36 specified-pressure nodes in the upper-right corner of the model area (fig. 23). At 100 weeks, 

lake-water withdrawal from a single surface node (fig. 23) is begun at a constant rate of 1.6 x 104 kg/s.

4000. m for horizontal flow

20. m for vertical flow
Chapter 3. Lake Capability 37

3.7.6. Initial Conditions

Initially, the groundwater flow is steady state, with no lakes on the surface of the model (fig. 24). Steady-state pressures are 

obtained through an extra initial simulation. Initial concentrations are zero throughout the model.

3.7.7. Results

Figure 25 shows areas of the model surface covered by lake water at a sequence of times during the simulation. Figure 26

shows how lake-water concentrations change with time.

At the start of the simulation, no lakes are on the surface of the model. Topographically driven, fresh groundwater discharges into the topographic depressions and starts to fill them with fresh lake water. By 30 weeks (fig. 26A), lakes 3 (right) 

and 4 (left) are both full and are spilling over into lake 5 (middle). At this time, the concentration associated with recharge from 

nodes in the upper-right corner of the model area (fig. 23) changes from 0. to 1. [dimensionless], creating a source of solute to 

the model. Solute-containing recharge discharges into lake 3, which in turn spills over into lake 5, and the concentration in both 

lakes increases (yellow and red curves in fig. 27).

At 47 weeks, lakes 4 and 5 coalesce to form lake 2 (convergence of blue and red curves to form the purple curve in fig. 27). 

The water in lake 2 is a mixture of the fresh water from lake 4 and the solute-containing water from lake 5. Because the volume 

of lake 5 is much larger than that of lake 4, the resulting concentration in lake 2 is close to that of lake 5. Figure 26B shows 

the areas covered by lakes 2 and 3 at 60 weeks. At 92 weeks, lakes 2 and 3 coalesce to form lake 1 (convergence of purple and 

yellow lines to form the black curve in fig. 27). The concentration in lake 1 results from mixing of the water in lake 2 with the 

smaller but more concentrated volume of water in lake 3.

Figure 26C shows the area of the model surface covered by lake 1 at 100 weeks. At this time, lake-water withdrawal 

begins, and the lake stage starts to decrease. At 119 weeks, lake 1 splits into its children—lakes 2 and 3 (splitting of the black 

curve into the yellow and purple curves in fig. 27). Lake-water withdrawal continues from lake 2. At 346 weeks, lake 2 splits 

into its children—lakes 4 and 5 (splitting of the purple curve into the blue and red curves in fig. 27). Thereafter, the concentration in lake 4 decreases as the water in lake 4 is flushed out by discharge of fresh groundwater into the lake. Lakes 3 and 4 are 

full and are spilling over into lake 5 as the stage in lake 5 declines due to continued lake-water withdrawal. As the stage in lake 

5 decreases, the head gradient driving groundwater discharge into the lake increases, and steady-state flow is approached as 

inflows into the lake balance withdrawals. Figure 26D shows the areas covered by lakes 3, 4, and 5 at 500 weeks, when the flow 

system is near steady state.

Figure 23. Hierarchy of five lakes identified in the lake example. 

Lakes 3, 4, and 5 correspond to the right, left, and middle 

topographic depressions (fig. 23).
38 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

Figure 24. Initial groundwater-flow velocity vectors in the top layer of elements for the lake example. Grid lines show the top 

surface of the model. Dots indicate element centroids. Velocity vectors point away from element centroids. The direction of each 

vector indicates the direction of flow at the element centroid, and the length of each vector is proportional to the magnitude of the 

velocity. (50X vertical exaggeration)
Chapter 3. Lake Capability 39

Figure 25. Areas of the model surface covered by lake water (shaded blue) in the lake example. Groundwater flow recharges at 

higher elevations and discharges at lower elevations (fig. 25), supplying water to form lakes in the topographic depressions. A, By 

30 weeks, lakes 3 (right) and 4 (left) are full and are spilling over onto lake 5 (middle). B, By 60 weeks, lakes 4 and 5 are coalesced 

into their parent, lake 2. C, By 100 weeks, all lakes are coalesced into one large lake (lake 1), and lake-water withdrawal begins. D, 

By 500 weeks, lake stages have declined such that there are once again three distinct lakes, and inflows nearly balance outflows in 

each lake, so the flow system is near steady state. (50X vertical exaggeration)
40 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

Acknowledgments

This work was funded in part by Strategic Environmental Research and Development Program (SERDP) grant RC2111. 

The authors are grateful to the technical reviewers of this report, which include Delwyn Oki and Richard Yager of the U.S. 

Geological Survey. 

References Cited

Al-Futaisi, A., and Patzek,T.W., 2002, Extension of the Hoshen-Kopelman algorithm to non-lattice environments, Physica A 

321, p. 665–678.

Freeze, R.A., and Cherry, J.A., 1979, Groundwater: Englewood Cliffs, N.J., Prentice Hall.

Harbaugh, A.W., 2005, MODFLOW-2005, The U.S. Geological Survey modular ground-water model—the Ground-Water Flow 

Process: U.S. Geological Survey Techniques and Methods, book 6, chap. A16, variously paginated, at http://pubs.usgs.gov/

tm/2005/tm6A16/.

Hoshen, J., and Kopelman, R., 1976, Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm: Physical Review B, v. 14, no. 8, p. 3438–3445.

Muskat, M., 1937, The flow of homogeneous fluids through porous media: New York, McGraw Hill.

Figure 26. Lake-water concentration as a function of time in the lake example. Concentration is expressed in arbitrary units in which 

the maximum source concentration is 1. When two lakes coalesce (lakes 4 and 5 at 47 weeks, and lakes 2 and 3 at 92 weeks), the 

water in the parent lake is a mixture of the concentrations in the two child lakes. When a lake splits into its children (lake 1 at 119 

weeks, and lake 2 at 346 weeks), the concentration in each child lake is initially the concentration from the parent lake. Over time, the 

concentrations in the two child lakes diverge. Lake-water concentration (arbitrary units)

Time, in weeks

10−6

0 50 100 150 200 250 300 350 400 450 500

10−5

10−4

10−3

10−2

10−1

100
References Cited 41

Polubarinova-Kochina, P., 1962, Theory of ground water movement: Princeton, N.J., Princeton Univ. Press. (Translated into 

English from the Russian by J.M. Roger de Wiest.)

Van Genuchten, M.Th., 1980, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils: Soil Science 

Society of America Journal, v. 44, no. 5, p. 892–898.

Voss, C.I., and Provost, A.M., 2002 (Version of September 22, 2010), SUTRA, A model for saturated-unsaturated variabledensity ground-water flow with solute or energy transport: U.S. Geological Survey Water-Resources Investigations 

Report 02–4231, 291 p., at http://pubs.er.usgs.gov/publication/wri024231.

Winston, R.B., and Voss, C.I., 2004, SutraGUI, a graphical user interface for SUTRA, a model for ground-water flow with solute 

or energy transport: U.S. Geological Survey Open-File Report 03–285, 114 p., at http://pubs.er.usgs.gov/publication/ofr03285.

Appendix 1. List of Symbols 43

Appendix 1. List of Symbols

The symbols used in this report are listed below. For each mathematical symbol, the listing includes the equation or section 

of the text in which the symbol was introduced, the units associated with the symbol, and a brief description of the quantity 

represented by the symbol.

Generic Units

[1] dimensionless

[E] energy units, or [M∙L2/s2]

[L] length units

[M] fluid mass units

[Ms] solute mass units

Units

[°C] degrees Celsius

[hr] hours

[kg] kilograms

[m] meters

[s] seconds

Special Notation

j i ≠

∑ (21) summation over neighbors of node i, not including node i itself

Greek Lowercase

α §1.10.4 [M/(L∙s2

)]-1 Porous matrix compressibility

αL §3.7.4 [L] Longitudinal dispersivity

αT §3.7.4 [L] Transverse dispersivity

β §1.10.4 [M/(L∙s2

)]-1 Fluid compressibility

ε §1.10.4 [1] Porosity

ηp (25) [1] Extremely large number used in implementing specified-pressure conditions in SUTRA 3.0

θ §1.4.2 [°] Angle above horizontal at which a conductive layer is inclined

μ (5) [M/(L∙s)] Fluid viscosity

νp (22) [L∙s] Symbol for the parameter GNUP used to control enforcement of 

specified-pressure conditions in previous versions of SUTRA

ρ (1) [M/L ] Fluid density

φ §3.2.3.1 [1] User-specified ratio of runoff to flow in or out of the groundwater model at a boundary node

Greek Uppercase

ΔHvap §1.4.6 [E/M] Latent heat of vaporization of water

Roman Lowercase

a §1.10.4 [L∙s2

/M] Van Genuchten model parameter; see equations (2.8) and (2.21) 

of Voss and Provost (2002)

aii (21) [L∙s] Diagonal element in row i of the coefficient matrix of the discretized flow equation

f

3
44 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

aii (28) [L∙s] Diagonal element in row i of the coefficient matrix of the discretized transport equation

aij (21) [L∙s] Element in row i and column j of the coefficient matrix of the 

discretized flow equation

aij (28) [M/s] or

[E/(°C∙s)]

Element in row i and column j of the coefficient matrix of the 

discretized transport equation

bi (21) [M/s] ith element of the right-hand side vector of the discretized flow 

equation in the absence of a specified-pressure condition at 

node i



bi

(28) [Ms/s] or

[E/s]

ith element of the right-hand side vector of the discretized transport equation in the absence of a specified-concentration or 

specified-temperature condition at node i

cw §1.4.6 [E/(M∙°C)] Specific heat of water

dex (12) [L] Extinction depth, below which the rate of fluid loss by evapotranspiration is zero

dmax (12) [L] Depth above which the rate of fluid loss by evapotranspiration is 

at its maximum (most negative) value

driv §1.4.5 [L] Depth of river water

g (1) [L/s2

] Gravitational acceleration (magnitude of the gravity vector)

g (A1) [L/s2

] Gravity vector

gn (6) [L/s2

] Gravitational acceleration normal to a conductive layer (normal 

component of the gravity vector)

h §1.10.7 [L] Height of water on the downstream side of the dam in the 

seepage-face example problem

hM (18) [Ms/(L2

∙s)] Mass-transfer coefficient

hT (16) [E/(L2

∙oC∙s)] Heat-transfer coefficient

k (5) [L2

] Permeability

k (A1) [L2

] Permeability tensor

kr §1.10.4 [1] Relative permeability

l §1.10.7 [L] Dam width in the seepage-face example problem

n §1.10.4 [1] Van Genuchten model parameter; see equations (2.8) and (2.21) 

of Voss and Provost (2002)

n (A2) [L] Coordinate normal to the surface of a conductive layer

p §1.1 [M/(L∙s2

)] Fluid pressure

p1 §1.2 [M/(L∙s2

)] First of two pressure values that define a generalized-flow condition

p2 §1.2 [M/(L∙s2

)] Second of two pressure values that define a generalized-flow 

condition

pBC §3.2.3.2 [M/(L∙s2

)] Specified pressure

pi (21) [M/(L∙s2

)] Computed pressure at node i

pj (21) [M/(L∙s2

)] Computed pressure at node j

pext (3) [M/(L∙s2

)] Pressure at an external point that is separated from the model 

boundary by a conductive layer

pwt (1) [M/(L∙s2

)] Pressure at the water table (assumed to be zero)

p(z) (1) [M/(L∙s2

)] Estimated pressure at a boundary-condition node at elevation z

pi

( ) 0 (26) [M/(L∙s2

)] Pressure computed at node i under the temporary assumption 

that Qi

PBC = 0
Appendix 1. List of Symbols 45

pi

BC (21) [M/(L∙s2

)] Specified pressure at boundary node i

q (A1) [M/(L∙s2

)] Mass flux vector

qn (A2) [M/(L2

∙s)] Component of the mass flux vector normal to the surface of a 

conductive layer

qET

max (14) [M/(L2

∙s)] Maximum (most negative) evapotranspiration (ET) flux

r3 (35) [L] Distance (radius) defined for the lakes example problem

r4 (36) [L] Distance (radius) defined for the lakes example problem

r5 (37) [L] Distance (radius) defined for the lakes example problem

t (31) [s] Time

t0 (31) [s] Time at which the current time step began

u § 3.2.2 [Ms/M] or

[°C]

Concentration or thermal energy content

ul (33) [Ms/M] or

[°C]

Concentration or thermal energy content in lake l

uk (33) [Ms/M] or

[°C]

Concentration or thermal energy content in lake k

ul

bc,in (33) [Ms/M] or

[°C]

Concentration or thermal energy content of inflow to lake l from 

groundwater-flow boundary conditions

x (35) [L] x coordinate

y (35) [L] y coordinate

z (1) [L] z coordinate (elevation)

ztop (34) [L] Top-surface elevation in the lakes example problem

zwt (1) [L] Water-table elevation

Roman Uppercase

A (4) [L2

] Area associated with a boundary-condition node

Clay (3) [L2

/s] Conductance of a conductive layer that separates an external 

point from the model boundary

C (18) [Ms/M] Solute concentration (mass fraction) in the fluid (mass of solute 

per total mass of fluid)

C1 (19) [Ms/M] First of two concentrations that define a generalized-transport 

condition for solute transport

C2 (19) [Ms/M] Second of two concentrations that define a generalized-transport 

condition for solute transport

Cext (18) [Ms/M] Concentration at an external point that is separated from the 

model boundary by a mass-transfer boundary layer

Csat §1.7.2 [Ms/M] Concentration of saturated brine

Dm §3.7.4 [L2

/s] Apparent molecular diffusivity of solute in solution in a porous 

medium, including tortuosity effects

Fl l

spill

, ′ (32) [M/s] Rate of spillover from lake l to its sibling, lake l’

Fk l

spill

, (32) [M/s] Rate of spillover from lake k to lake l

Fl

bc,in (32) [M/s] Rate of inflow to lake l from groundwater-flow boundary conditions

Fl

bc,out (32) [M/s] Rate of outflow from lake l to groundwater-flow boundary conditions
46 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

Fl

lim (32) [M/s] Rate of outflow from lake l due to lake-area limitations

Fl

net (31) [M/s] Net rate of inflow of water to or outflow of water from lake l

Gl

bc,in (33) [Ms/M]

[E/s]

Rate of inflow of solute mass or energy to lake l from groundwater-transport boundary conditions

Gl

bc,out (33) [Ms/M]

[E/s]

Rate of outflow of solute mass or energy from lake l to groundwater-transport boundary conditions

H §1.10.7 [L] Height of water on the upstream side of the dam in the seepageface example problem

K (4) [L/s] Hydraulic conductivity

L (4) [L] Thickness of a conductive layer

Ml (31) [M] Mass of water in lake l at time t

Ml

0 (31) [M] Mass of water in lake l at time t0

NAL (32) [1] Number of active lakes

Q §1.1 [M/s] Rate of fluid-mass gain or loss (flow)

Q1 §1.2 [M/s] First of two rates of fluid-mass gain or loss (flow) that define a 

generalized-flow condition

Q2 §1.2 [M/s] Second of two rates of fluid-mass gain or loss (flow) that define a 

generalized-flow condition

Qdam §1.10.7 [L3

/s per unit thickness] Discharge through the dam in the seepage-face example problem

Qfd (3) [M/s] Rate of free (gravity-driven) drainage of fluid across a conductive layer

QIN §3.2.3.1 [M/s] Rate of fluid-mass gain or loss (flow) at a fluid source or sink 

node

Qinf (15) [M/s] Fluid-mass infiltration rate associated with a seepage-face 

boundary condition

Qmax (9) [M/s] Maximum rate of fluid-mass gain or loss (flow) associated with a 

river boundary condition

QPBC §3.2.3.2 [M/s] Rate of fluid-mass gain or loss (flow) at a specified-pressure 

node

QPBG §3.2.3.3 [M/s] Rate of fluid-mass gain or loss (flow) at a generalized-flow node

QPBG1 §1.2 [M/s] Symbol used to represent Q1 in input datasets that define 

generalized-flow conditions

QPBG2 §1.2 [M/s] Symbol used to represent Q2 in input datasets that define 

generalized-flow conditions

Qriv (9) [M/s] Rate of fluid-mass gain or loss (flow) driven by the head supplied by a column of river water

QUINC §3.2.3.5 [Ms/s] or

[E/s]

Rate of gain or loss of solute mass or energy at a solute mass or 

energy source/sink node

QU §1.6 [Ms/s] or

[E/s]

Rate of gain or loss of solute mass or energy

QU

1 §1.6 [Ms/s] or

[E/s]

First of two rates of gain or loss of solute or energy that define a 

generalized-transport condition

QU

2 §1.6 [Ms/s] or

[E/s]

Second of two rates of gain or loss of solute or energy that define 

a generalized-transport condition

Qi

PBC (21) [M/s] Rate of gain or loss of fluid mass associated with a specifiedpressure condition at node i

Qi

UBC §2.3 [Ms/s] or

[E/s]

Rate of gain or loss of solute or energy associated with a 

specified-concentration or specified-temperature condition at 

node i
Appendix 1. List of Symbols 47

QUB

C

U §3.2.3.6 [M

s/s] or

[E/s]

Rate of gain or loss of solute mass or energy at a specifiedconcentration or specified-temperature node

QUB

G

U §3.2.3.7 [M

s/s] or

[E/s]

Rate of gain or loss of solute mass or energy at a generalizedtransport node

QUB

G

U

1 §1.6 [M

s/s] or

[E/s]

Symbol used to represent 

Q

U1 in input datasets that define 

generalized-transport conditions

QUB

G

U

2 §1.6 [M

s/s] or

[E/s]

Symbol used to represent 

Q

U2 in input datasets that define 

generalized-transport conditions

Sop §1.10.4 [M/(L∙s

2

)]-1 Specific pressure storativity

Sw §1.10.4 [1] Water saturation (saturation) 

(volume of water per volume of voids)

Swres §1.10.4 [1] Residual water saturation T (16) [°C] Fluid temperature T1 (17) [°C] First of two temperatures that define a generalized-transport 

condition for energy transport

T2 (17) [°C] Second of two temperatures that define a generalized-transport 

condition for energy transport

Text (16) [°C] Fluid temperature at an external point that is separated from the 

model boundary by a thermal boundary layer

U (1.6) [°C] or

[M

s/M]

Solute concentration or temperature

U1 §1.6 [°C] or

[M

s/M]

First of two concentrations or temperatures that define a 

generalized-transport condition

U2 §1.6 [°C] or

[M

s/M]

Second of two concentrations or temperatures that define a 

generalized-transport condition

UBC §3.2.3.6 [°C] or

[M

s/M]

Specified concentration or temperature

UBG1 §1.6 [°C] or

[M

s/M]

Symbol used to represent U1 in input datasets that define 

generalized-transport conditions

UBG2 §1.6 [°C] or

[M

s/M]

Symbol used to represent U2 in input datasets that define 

generalized-transport conditions

Ui §2.3 [°C] or

[M

s/M]

Computed concentration or temperature at node i

Uj (28) [°C] or

[M

s/M]

Computed concentration or temperature at node j

Ui

( ) 0 (28) [°C] or

[M

s/M]

Concentration or temperature computed at node i under the 

temporary assumption that 

Qi

UBC

=0

Ui

BC §2.3 [°C] or

[M

s/M]

Specified concentration or temperature at boundary node i
48 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

Appendix 2. Flow Across a Conductive Layer

The mass flux vector for groundwater flow within the conductive layer, q [M/(L2·s)], is given by

q

k = − ( ) ∇ + p gρ

μ

ρ (A1)

where ρ is fluid density [M/L3], μ is fluid viscosity [M/(L·s)], k is the permeability tensor [L2] within the layer, g is the gravity 

vector [L/s2], and p is pressure [M/(L·s2)]. If the layer is assumed to be isotropic, and permeability, density, viscosity, and gravity 

can be considered uniform throughout the layer, the component of the flux normal to the surface of the layer, qn, can be written 

as

q k p

n

g n n = − ∂

∂

+ 



 



 ρ

μ

ρ (A2)

and integrated across the thickness of the layer to give

q L k n = − ( ) p p ext − + g Ln

ρ

μ

ρ (A3)

Here, n is the coordinate [L] in the normal direction, L is the thickness [L] of the layer, k is the isotropic permeability, gn is the 

component of gravity normal to the surface of the layer, p [M/(L·s2)] is the pressure on the side of the layer adjacent to the 

model boundary, and pext [M/(L·s2)] is the “external” pressure on the other side of the layer. Multiplication of both sides of (A3) 

by the area, A [L2], associated with the boundary node and rearrangement yields an expression for the total mass flow rate of 

groundwater through the layer, Q [M/s]:

Q kA

L p p

kA g ext

n = − ( ) + ρ

μ

ρ

μ

2

(A4)

which can be rewritten as (3)-(6). If permeability, density, viscosity, or gravity varies spatially within the intervening medium, 

choosing a representative, constant value preserves the linear form of equation (A4).
Appendix 3. Input Data List 49

Appendix 3. Input Data List

This section describes SUTRA input datasets that are new or have been modified in Version 3.0. The datasets described 

here follow the general rules for input described in the SUTRA documentation (Voss and Provost, 2002). A complete, up-to-date 

description of all SUTRA input datasets is maintained on the USGS Web site at https://doi.org/10.5066/P9PPEHHM.

List of New or Modified Input Data for the Main Input File (.inp)

Model Version: SUTRA 3.0

NOTE: “.bcopg”, “.bcoug”, “.lkin”, “.lkbc”, “.lkar”, “.lkbu”, “.lkst”, “.lkrs”, “.lkn”, and “.lkh” files have been added.

The file “SUTRA.FIL” contains file assignments (one line for each assignment) in the following format:

Variable Type Description

FTYPE Character File type (within single quotes ‘ ‘ as shown below). Valid values are as follows:

‘INP’ = “.inp” input file (main input)

‘BCS’ = “.bcs” input file (time-dependent sources and boundary conditions)

‘ICS’ = “.ics” input file (initial conditions)

‘LKIN’ = “.lkin” input file (lake main input)

‘LKBC’ = “.lkbc” input file (lake-boundary condition interactions)

‘LKAR’ = “.lkar” input file (lake-area nodes)

‘LST’ = “.lst” output file (main output listing)

‘RST’ = “.rst” output file (restart conditions)

‘NOD’ = “.nod” output file (nodewise results)

‘ELE’ = “.ele” output file (elementwise results)

‘OBS’ = “.obs” output file (observations)

‘OBC’ = “.obc” output file (observations)

‘BCOF’ = “.bcof” output file (specifications and results at fluid-source/sink nodes)

‘BCOS’ = “.bcos” output file (specifications and results at solute/energy-source/sink nodes)

‘BCOP’ = “.bcop” output file (specifications and results at specified-pressure nodes)

‘BCOU’ = “.bcou” output file (specifications and results at specified-concentration/temperature 

nodes)

‘BCOPG’ = “.bcopg” output file (specifications and results at generalized-flow nodes)

‘BCOUG’ = “.bcoug” output file (specifications and results at generalized-transport nodes)

‘LKBU’ = “.lkbu” output file (lake budgets)

‘LKST’ = “.lkst” output file (lake stages)

‘LKRS’ = “.lkrs” output file (lake restart) 

‘LKN’ = “.lkn” output file (lake numbers for surface nodes) 

‘LKH’ = “.lkh” output file (lake hierarchy)

‘SMY’ = “.smy” output file (simulation summary)

IUNIT Integer FORTRAN unit number to be assigned to the file. If IUNIT is not a valid FORTRAN unit number 

or if it is already assigned to another file, SUTRA will assign the next available unit number after 

IUNIT. (Unit numbers less than 11 are assumed to be unavailable.)

FNAME Character Full name of the file (within single quotes ‘ ‘ as shown below).
50 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

Notes:

Assignments for the “.bcs”, “.nod”, “.ele”, “.obs”, “.obc”, “.bcof”, “.bcos”, “.bcop”, “.bcou”, and “.lkst” files are optional. If 

any of these assignments are omitted, the corresponding output files will not be created by SUTRA. Assignment for the “.smy” 

file is also optional – if not assigned, it will receive the file name “SUTRA.SMY” and an automatically generated unit number, 

by default. Assignment for the “.rst” file is required if ISTORE ≠ 0 in dataset 4, and assignment for the “.lkrs” file is required 

in addition if the lake capability is used. Assignments for the “.inp”, “.ics”, and “.lst” files are always required. Assignments 

for the “.lkin”, “.lkbc”, “.lkbu”, “.lkn”, and “.lkh” files are always required when the lake capability is used. Assignment for 

the “.lkar” file is optional – if the lake capability is used and the “.lkar” file is not assigned, all nodes on the top surface of the 

model are eligible to form lakes. Assignments may be listed in any order. Assignment of unit numbers is performed by SUTRA 

in the order in which the files are listed, except for the “.obs” and “.obc” files. The latter are always assigned last by SUTRA.

For the observation output files (“.obs” and “.obc”), FNAME is a base filename from which the actual filenames are automatically derived by SUTRA. SUTRA generates one observation output file for each combination of schedule and output format that 

appears in the observation specifications in dataset 8D of the main input (“.inp”) file. For an example of this, see the description 

of dataset 8D.

Multiple “.bcs” files may be assigned. At the beginning of each time step, SUTRA checks each “.bcs” file and reads in any 

boundary condition specifications that pertain to that time step. The “.bcs” files are checked in the order in which they are listed 

in “SUTRA.FIL”. If the “.bcs” files contain multiple specifications of the same kind, at the same node, on the same time step, 

the specification read in last overrides any previous specifications.

Example:

‘INP’ 50 ‘project.inp’

‘BCS’ 52 ‘project.bcs’

‘ICS’ 55 ‘project.ics’

‘LST’ 60 ‘project.lst’

‘RST’ 66 ‘project.rst’

‘NOD’ 70 ‘project.nod’

‘ELE’ 75 ‘project.ele’

‘OBS’ 80 ‘project.obs’

‘OBC’ 85 ‘project.obc’

‘BCOF’ 91 ‘project.bcof’

‘BCOS’ 93 ‘project.bcos’

‘BCOP’ 92 ‘project.bcop’

‘BCOU’ 94 ‘project.bcou’

‘BCOPG’ 95 ‘project.bcopg’

‘BCOUG’ 96 ‘project.bcoug’

‘SMY’ 98 ‘project.smy’

‘LKIN’ 40 ‘project.lkin’

‘LKBC’ 41 ‘project.lkbc’

‘LKAR’ 42 ‘project.lkar’

‘LKBU’ 43 ‘project.lkbu’

‘LKST’ 44 ‘project.lkst’

‘LKRS’ 45 ‘project.lkrs’

‘LKN’ 46 ‘project.lkn’

‘LKH’ 47 ‘project.lkh’
Appendix 3. Input Data List 51

DATASET 3: Simulation Control Numbers (one line)

NOTE: Input variables NPBG and NUBG have been added.

Variable Type Description

NN Integer Exact number of nodes in finite element mesh.

NE Integer Exact number of elements in finite element mesh.

NPBC Integer Exact number of nodes at which pressure is a specified constant value or function of time.

NUBC Integer Exact number of nodes at which temperature or concentration is a specified constant value or 

function of time.

NSOP Integer Exact number of nodes at which a fluid source/sink is a specified constant value or function of time.

NSOU Integer Exact number of nodes at which an energy or solute mass source/sink is a specified constant value or 

function of time.

NPBG Integer Exact number of nodes at which a generalized-flow condition is specified.

NUBG Integer Exact number of nodes at which a generalized-transport condition is specified.

NOBS Integer Exact number of points at which observations will be made. Set to zero for no observations.

Note:

The order in which NPBC, NUBC, NSOP, and NSOU are listed above is different from the order in which the corresponding 

datasets (19, 20, 17, and 18) appear in the “.inp” file.

DATASET 5: Upstream Weighting Parameter (one line)

NOTE: Parameters GNUP and GNUU are no longer needed and have been removed.

Variable Type Description 

UP Real Fractional upstream weight for stabilization of oscillations in results due to highly advective 

transport or unsaturated flow. UP may be given any value from 0.0 to +1.0. UP = 0.0 implies no 

upstream weighting (Galerkin method). UP = 0.5 implies 50% upstream weighting. UP = 1.0 

implies full (100%) upstream weighting. Recommended value is zero. 

WARNING: Upstream weighting increases the local effective longitudinal dispersivity of the simulation by approximately UP* (ΔLL)/2 where ΔLL is the local distance between element sides along 

the direction of flow. (See section 7.2 of the main SUTRA documentation (Voss and Provost, 

2002)). Note that the amount of this increase varies from place to place, depending on flow direction and element size. Thus, a nonzero value for UP actually changes the value of longitudinal 

dispersivity used by the simulation and broadens otherwise sharp concentration, temperature, or 

saturation fronts.

DATASET 6: Temporal Control and Solution Cycling Data

NOTE: The list of input parameters remains unchanged from SUTRA Version 2.2. However, schedule “TIME_STEPS” 

must now be defined before any other schedules are defined, and time-step numbers and times must now be listed in 

ascending order within each ‘STEP LIST’ or ‘TIME LIST’ schedule.
52 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

DATASET 8E: Output Controls and Options for “.bcof”, “.bcos”, “.bcop”, “.bcou”, “.bcopg”, and “.bcoug” Files (Source and 

Boundary Condition Specifications and Results Listed in Columns) (one line)

NOTE: Input variables NBGPPR and NBGUPR have been added.

Variable Type Description

NBCFPR Integer Fluid source node output cycle for transient simulation. Output is produced in the “.bcof” file on 

time steps numbered n│NBCFPR│ (where n is a positive integer).

NBCSPR Integer Solute or energy source node output cycle for transient simulation. Output is produced in the “.bcos” 

file on time steps numbered n│NBCSPR│ (where n is a positive integer).

NBCPPR Integer Specified pressure node output cycle for transient simulation. Output is produced in the “.bcop” file 

on time steps numbered n│NBCPPR│ (where n is a positive integer).

NBCUPR Integer Specified concentration or temperature node output cycle for transient simulation. Output is 

produced in the “.bcou” file on time steps numbered n│NBCUPR│ (where n is a positive 

integer).

NBGPPR Integer Generalized-flow node output cycle for transient simulation. Output is produced in the “.bcopg” file 

on time steps numbered n│NBGPPR│ (where n is a positive integer).

NBGUPR Integer Generalized-transport node output cycle for transient simulation. Output is produced in the “.bcoug” 

file on time steps numbered n│NBGUPR│ (where n is a positive integer).

CINACT Character Set to a value of ‘Y’ to list all sources and boundary conditions, including those that have been 

rendered inactive (using negative node numbers in a “.bcs” file; see the description of “.bcs” 

datasets 3–7). Set to ‘N’ to suppress output for inactive sources and boundary conditions.

Notes:

For transient simulations, output is also produced on the first and last time steps. For steady-state simulations, output is produced 

irrespective of the values of the output control variables.

The information that is output to the “.bcof”, “.bcos”, “.bcop”, and “.bcou” files is described in section 5.9 of the main SUTRA 

documentation (Voss and Provost, 2002). The information that is output to the “.bcopg” and “.bcoug” files is described in section 1.8 of this document, “Budget Output and Boundary-Condition Output Files.”

Example:

To output source and boundary condition information to the “.bcof”, “.bcos”, “.bcop”, “.bcou”, “.bcopg”, and“.bcoug” files 

every 10 time steps, suppressing output for inactive sources and boundary conditions, write the following:

10 10 10 10 10 10 ‘N’
Appendix 3. Input Data List 53

DATASET 21A: Constant Values for Generalized-Flow Nodes (one line for each of NPBG generalized pressure nodes as indicated in dataset 3, plus one line)

NOTE: This is an entirely new dataset.

O M I T when there are no generalized-flow nodes (NPBG = 0).

Specifications made in this dataset are constant values that are used by SUTRA unless/until they are superseded by timedependent specifications in dataset 7A of an optional “.bcs” file.

Variable Type Description

Lines 1 to NPBG

IPBG Integer Number of node to which generalized-flow data on this line refers. Specifying the node number with 

a negative sign indicates to SUTRA that the manner in which the input variable values at this 

node vary with time has been programmed by the user in subroutine BCTIME.

PBG1 Real Pressure value at the first of two points that define the linear relation between fluid inflow/outflow 

and pressure at node IPBG. [M/(L∙s2)]. 

May be omitted if IPBG<0.

QPBG1 Real Rate of inflow/outflow of fluid mass when the computed pressure equals PBG1 at node IPBG. [M/s].

May be omitted if IPBG < 0.

PBG2 Real Pressure value at the second of two points that define the linear relation between fluid inflow/outflow 

and pressure at node IPBG. [M/(L∙s2)]. 

May be omitted if IPBG < 0.

QPBG2 Real Rate of inflow/outflow of fluid mass when the computed pressure equals PBG2 at node IPBG. [M/s]. 

May be omitted if IPBG < 0.

CPQL1 Character Determines whether a limit is set on fluid inflow/outflow, pressure, or neither at the first of two 

points that define the linear relation between fluid inflow/outflow and pressure at node IPBC. 

[M/(L∙s2)]. Valid values are as follows:

‘Q’ = set limit of Q ≤ QPBG1

‘P’ = set limit of p ≥ PBG1

‘N’ = set no limit at point 1

May be omitted if IPBG < 0.

CPQL2 Character Determines whether a limit is set on fluid inflow/outflow, pressure, or neither at the second of two 

points that define the linear relation between fluid inflow/outflow and pressure at node IPBC. 

[M/(L∙s2)]. Valid values are as follows:

‘Q’ = set limit of Q ≥ QPBG2

‘P’ = set limit of p ≤ PBG2

‘N’ = set no limit at point 2

May be omitted if IPBG < 0.

UPBGI Real Temperature or solute concentration of any external fluid that enters the model at node IPBC. 

UPBGI is a specified constant value. [°C] or [Ms/M].

May be omitted if IPBG < 0.

CUPBGO Character One word. Must be either ’DIR‘ or ’REL‘. Indicates whether the temperature or solute concentration, UPBGO, of any fluid that exits the model at node IPBG is being specified directly (‘DIR’) or 

relative to the computed concentration or temperature at the node (‘REL’).

May be omitted if IPBG < 0.

UPBGO Real Temperature or solute concentration of any fluid that exits the model at node IPBG. If CUPBGO = 

’DIR’, the exit temperature or concentration is set directly to the value UPBGO. If CUPBGO = 

’REL’, the exit temperature or concentration is set to the value U+UBPGO, where U is the temperature or concentration computed at the node. [°C] or [Ms/M].

May be omitted if IPBG < 0.

Last line

Integer Placed immediately following all NPBG generalized-flow lines. Line must begin with the integer 0.
54 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

DATASET 21B: Constant Values for Generalized-Transport Nodes (one line for each of NUBG generalized-transport nodes 

indicated in dataset 3, plus one line)

NOTE: This is an entirely new dataset.

O M I T when there are no generalized-transport nodes (NUBG=0).

Specifications made in this dataset are constant values that are used by SUTRA unless/until they are superseded by timedependent specifications in dataset 7B of an optional “.bcs” file.

Variable Type Description

Lines 1 to NUBG

IUBG Integer Number of node to which generalized-transport data on this line refers. Specifying the node number 

with a negative sign indicates to SUTRA that the manner in which the input variable values at 

this node vary with time has been programmed by the user in subroutine BCTIME.

UBG1 Real Concentration or temperature value at the first of two points that define the linear relation between solute mass or energy inflow/outflow and concentration or temperature at node IUBG. 

[Ms/M or C°].

May be omitted if IUBG < 0.

QUBG1 Real Rate of inflow/outflow of solute mass or energy when the computed concentration or temperature 

equals UBG1 at node IUBG. [Ms/s or E/s].

May be omitted if IUBG < 0.

UBG2 Real Concentration or temperature value at the second of two points that define the linear relation 

between solute mass or energy inflow/outflow and concentration or temperature at node IUBG. 

[Ms/M or C°].

May be omitted if IUBG < 0.

QUBG2 Real Rate of inflow/outflow of solute mass or energy when the computed concentration or temperature 

equals UBG2 at node IUBG. [Ms/s or E/s].

May be omitted if IUBG < 0.

Last line

Integer Placed immediately following all NUBG specified temperature or concentration lines. Line must 

begin with the integer 0.
Appendix 3. Input Data List 55

List of New or Modified Input Data for the Time-Dependent Sources and Boundary Conditions Files 

(.bcs)

Model Version: SUTRA 3.0

DATASET 2: Identifier and Time-Dependent Source and Boundary Condition Counts (one line) 

NOTE: Input variables NPBG1 and NUBG1 have been added.

Required for each time step on which there is a source or boundary condition specification in the “.bcs” file (that is, for 

each time step in the controlling schedule specified in dataset 1).

Variable Type Description

BCSID Character Name that identifies the set of sources and (or) boundary conditions that follows in datasets 3– 7. 

May be up to 40 characters long and may include spaces.

NSOP1 Integer Exact number of nodes at which a fluid source/sink is specified for the current time step.

NSOU1 Integer Exact number of nodes at which an energy or solute mass source/sink is specified for the current 

time step.

NPBC1 Integer Exact number of nodes at which pressure is specified for the current time step.

NUBC1 Integer Exact number of nodes at which temperature or concentration is specified for the current time step.

NPBG1 Integer Exact number of nodes at which a generalized-flow condition is specified for the current time step.

NUBG1 Integer Exact number of nodes at which a generalized-transport condition is specified for the current time 

step.

Note:

The source and boundary condition counts NSOP1, NSOU1, NPBC1, NUBC1, NPBG1, and NUBG1 correspond to “.bcs” datasets 3, 4, 5, 6, 7A, and 7B, respectively, and must be listed in that order. The counts refer to specifications made on the current 

time step in the current “.bcs” file; specifications made on other time steps and (or) in other files are not included in the counts.

Specifications need be included in “.bcs” files only when they are changing on a given time step. Specifications are persistent – 

once set, they remain in force unless/until superseded by another specification.

If the controlling schedule includes time steps for which no “.bcs” specification is desired, one can simply set NSOP1=NSOU1=

NPBC1=NUBC1=NPBG1=NUBG1=0 on those time steps.

To assist the user in keeping track of where various boundary conditions are set, SUTRA includes the identifier BCSID with 

each source or boundary condition listed in the “.bcof”, “.bcos”, “.bcop”, “.bcou”, “.bcopg”, and “.bcoug” output files.
56 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

DATASET 7A: Time-Dependent Generalized Flows (one line for each of NPBG1 generalized-flow nodes as indicated in “.bcs” 

dataset 2, plus one line)

NOTE: This is an entirely new dataset.

O M I T when there are no generalized flows to be specified for the current time step (NPBG1 = 0 in dataset 2).

Specifications made in this dataset take precedence over those made in “.inp” dataset 21A (and subroutine BCTIME). A 

specification remains in force unless/until superseded by a corresponding specification on a later time step (or in another 

“.bcs” file read later on the same time step).

Variable Type Description

Lines 1 to NPBG1

IPBG1 Integer Number of node to which generalized-flow data on this line refers. Specifying the node number with 

a negative sign renders the generalized-flow condition inactive (it is not applied) at node IPBG1, 

which gives the user the option of suppressing output for the inactive boundary condition by setting CINACT = ’N’ in “.inp” dataset 8E.

PBG11 Real Pressure value at the first of two points that define the linear relation between fluid inflow/outflow 

and pressure at node IPBG1. [M/(L⋅s2)]. 

May be omitted if the boundary condition is inactive (IPBG1 < 0).

QPBG11 Real Rate of inflow/outflow of fluid mass when the computed pressure equals PBG11 at node IPBG1. 

[M/s]. 

May be omitted if the boundary condition is inactive (IPBG1 < 0).

PBG21 Real Pressure value at the second of two points that define the linear relation between fluid inflow/outflow 

and pressure at node IPBG1. [M/(L⋅s2)].

May be omitted if the boundary condition is inactive (IPBG1 < 0).

QPBG21 Real Rate of inflow/outflow of fluid mass when the computed pressure equals PBG21 at node IPBG1. 

[M/s].

May be omitted if the boundary condition is inactive (IPBG1 < 0).

CPQL11 Character Determines whether a limit is set on fluid inflow/outflow, pressure, or neither at the first of two points 

that define the linear relation between fluid inflow/outflow and pressure at node IPBG1. [M/(L∙s2)]. 

Valid values are as follows:

‘Q’ = set limit of Q ≤ QPBG11

‘P’ = set limit of p ≥ PBG1

‘N’ = set no limit at point 1

May be omitted if the boundary condition is inactive (IPBG1 < 0).

CPQL21 Character Determines whether a limit is set on fluid inflow/outflow, pressure, or neither at the second of two 

points that define the linear relation between fluid inflow/outflow and pressure at node IPBG1. 

[M/(L∙s2)]. Valid values are as follows:

‘Q’ = set limit of Q ≥ QPBG21

‘P’ = set limit of p ≤ PBG21

‘N’ = set no limit at point 2

May be omitted if the boundary condition is inactive (IPBG1 < 0).

UPBGI1 Real Temperature or solute concentration of any external fluid that enters the model at node IPBG1. 

UPBGI1 is a specified constant value. [°C] or [Ms/M].

CUPBGO1 Character One word. Must be either ‘DIR’ or ‘REL’. Indicates whether the temperature or solute concentration, 

UPBGO1, of any fluid that exits the model at node IPBG1 is being specified directly (‘DIR’) or 

relative to the computed concentration or temperature at the node (‘REL’).

May be omitted if the boundary condition is inactive (IPBG1 < 0).

UPBGO1 Real Temperature or solute concentration of any fluid that exits the model at node IPBG1. If CUPBGO1 

= ‘DIR’, the exit temperature or concentration is set directly to the value UPBGO1. If CUPBGO1 

= ‘REL’, the exit temperature or concentration is set to the value U+UBPGO1, where U is the 

temperature or concentration computed at the node. [°C] or [Ms/M].

May be omitted if the boundary condition is inactive (IPBG1 < 0).
Appendix 3. Input Data List 57

Last line

Integer Placed immediately following all NPBG1 specified pressure lines. Line must begin with the integer 

0.

Note:

It is not sufficient to specify only the input variable that is changing on the current time step. For example, if PBG11 is changing 

but the other input variables (for example, QPBG11) are not, one must still provide a value for each input variable.

DATASET 7B: Time-Dependent Generalized Transport (one line for each of NUBG1 generalized-transport nodes indicated in 

“.bcs” dataset 2, plus one line)

NOTE: This is an entirely new dataset.

O M I T when there are no generalized-transport nodes to be specified for the current time step (NUBG1 = 0 in dataset 

2).

Specifications made in this dataset take precedence over those made in “.inp” dataset 21B (and subroutine BCTIME). A 

specification remains in force unless/until superseded by a corresponding specification on a later time step (or in another 

“.bcs” file read later on the same time step).

Variable Type Description

Lines 1 to NUBG1

IUBG1 Integer Number of node to which generalized-transport data on this line refers. Specifying the node number 

with a negative sign renders the generalized-transport condition inactive (it is not applied) at 

node IUBG1, which gives the user the option of suppressing output for the inactive boundary 

condition by setting CINACT = ’N’ in “.inp” dataset 8E.

Any node listed in this dataset must also be listed in “.inp” dataset 21B.

UBG11 Real Concentration or temperature value at the first of two points that define the linear relation between 

solute mass or energy inflow/outflow and concentration or temperature at node IUBG1. [Ms/M or 

C°]. 

May be omitted if the boundary condition is inactive (IUBG1 < 0).

QUBG11 Real Rate of inflow/outflow of solute mass or energy when the computed concentration or temperature 

equals UBG11 at node IUBG1. [Ms/s or E/s].

May be omitted if the boundary condition is inactive (IUBG1 < 0).

UBG21 Real Concentration or temperature value at the second of two points that define the linear relation 

between solute mass or energy inflow/outflow and concentration or temperature at node IUBG1. 

[Ms/M or C°].

May be omitted if the boundary condition is inactive (IUBG1 < 0).

QUBG21 Real Rate of inflow/outflow of solute mass or energy when the computed concentration or temperature 

equals UBG21 at node IUBG1. [Ms/s or E/s].

May be omitted if the boundary condition is inactive (IUBG1 < 0).

Last line

Integer Placed immediately following all NUBG1 specified temperature or concentration lines. Line must 

begin with the integer 0.

Note:

It is not sufficient to specify only the input variable that is changing on the current time step. For example, if UBG11 is changing 

but the other input variables (for example, QUBG11) are not, one must still provide a value for each input variable.
58 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

Input Data for the Main Input File for the Lake Capability (.lkin)

Model Version: SUTRA 3.0

NOTE: This is an entirely new input file.

This input file is optional in general but is required for the lake capability. If this input file is omitted, the lake capability 

is not used.

DATASET 1: Lake Output Control (one line)

Variable Type Description

NLAKPR Integer Output cycle for lake information during a transient simulation. Output is produced in the lake output files on time steps numbered n│NLAKPR│(where n is a positive integer). Also, for transient 

solutions, output is produced for initial conditions and on the first and last time steps. To cancel 

printed output for the first time step of a transient simulation, set NLAKPR to a negative number 

(that is, place a minus sign before the desired output cycle). For steady-state solutions, output is 

produced irrespective of the value of NLAKPR.

DATASET 2: Number of lake specifications and default values (one line)

Variable Type Description

NLSPEC Integer Number of lake specifications. Set NLSPEC = 0 if there are no lake specifications.

FRROD Real Default value of the fraction of groundwater recharge diverted as runoff from generalized-flow nodes. (See 

section 3.2.3.4, “Runoff to Lakes.”) This default value is applied to the basins of all lakes not identified 

in the lake specifications in DATASET 3 below.

FDROD Real Default value of the fraction of groundwater discharge diverted as runoff from generalized-flow nodes. 

(See section 3.2.3.4, “Runoff to Lakes.”) This default value is applied to the basins of all lakes not identified in the lake specifications in DATASET 3 below.

RNOLK Real Value output (in place of lake stage) to the “.lkst” output file at nodes not submerged by lake water. May be 

set to any real value.

DATASET 3: Lake specifications (one line for each of NLSPEC optional lake specifications)

Variable Type Description

CTYPE Character One word. Must be either ‘LAKE’ or ‘NODE’.

ILON Integer Number of the lake or node to which the specification applies. If CTYPE = ‘LAKE’, ILON is the number of 

the lake to which the specification applies. If CTYPE = ‘NODE’, ILON is the number of a node within 

the basin of the lake to which the specification applies. Node ILON must be on the top surface of the 

model.

STGI Real Initial stage in the lake to which the specification applies. If STGI is below the lowest point in the lake 

bathymetry, the lake is initially dry. SUTRA resolves any conflicting initial stage specifications in favor 

of the highest stage specified.

UWI Real Initial lake-water concentration or temperature in the lake to which the specification applies. [Ms/M or C°]

FRRO Real Value of the fraction of groundwater recharge diverted as runoff from generalized-flow nodes within the 

basin of the lake to which the specification applies. (See section 3.2.3.4, “Runoff to Lakes.”)

FDRO Real Value of the fraction of groundwater discharge diverted as runoff from generalized-flow nodes within the 

basin of the lake to which the specification applies. (See section 3.2.3.4, “Runoff to Lakes.”)
Appendix 3. Input Data List 59

Note:

Lakes are identified and numbered automatically by SUTRA. Therefore, the user generally will not know the lake numbering 

scheme until a preliminary run is performed. Once the lake numbering is known, CTYPE = ‘LAKE’ can be used to allow lake 

specifications to be made conveniently using lake numbers. Until then, the user has the option to make specifications using node 

numbers by setting CTYPE = ‘NODE’. Setting the desired specifications can involve trial and error, especially for complex lake 

systems.

Input Data for the Boundary-Condition Interaction Input File for the Lake Capability (.lkbc)

Model Version: SUTRA 3.0

NOTE: This is an entirely new input file.

DATASET 1: Boundary-condition interaction specification counts (one line)

Variable Type Description

NBCIF Integer Number of interaction specifications for fluid sources/sinks.

NBCIS Integer Number of interaction specifications for sources/sinks of solute mass or energy.

NBCIP Integer Number of interaction specifications for specified pressures.

NBCIU Integer Number of interaction specifications for specified concentrations or temperatures.

NBCIPG Integer Number of interaction specifications for generalized-flow conditions.

NBCIUG Integer Number of interaction specifications for generalized-transport conditions.

DATASET 2: Boundary-condition interaction specifications for fluid sources/sinks (one line for each of NBCIF specifications)

Variable Type Description

BCSFNM Character On the first line, set BCSFNM = ‘DEFAULT’. SUTRA applies DEFAULT specifications to all fluidsource/sink nodes that can form lakes but are not identified as fluid sources/sinks in the optional list of 

“.bcs” files that follows.

On each subsequent line (optional), set BCSFNM to the name of a “.bcs” file that contains fluid sources/

sinks with which lakes will interact.

ILKF Integer Effect that the presence or absence of lake water has on the pplication of fluid sources/sinks identified by 

BCSFNM. When ILKF = -1, a fluid source/sink is applied to a node only when lake water is absent at 

that node. When ILKF=0, a fluid source/sink is a applied to a node whether lake water is present or absent at that node. When ILKF=1, a fluid source/sink is applied to a node only when lake water is present 

at that node.
60 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

DATASET 3: Boundary-condition interaction specifications for sources/sinks of solute mass or energy (one line for each of 

NBCIS specifications)

Variable Type Description

BCSFNM Character On the first line, set BCSFNM = ‘DEFAULT’. SUTRA applies DEFAULT specifications to all solute mass 

or energy-source/sink nodes that can form lakes but are not identified as sources/sinks of solute mass or 

energy in the optional list of “.bcs” files that follows.

On each subsequent line (optional), set BCSFNM to the name of a “.bcs” file that contains sources/sinks 

of solute mass or energy with which lakes will interact.

ILKS Integer Effect that the presence or absence of lake water has on the application of sources/sinks of solute mass or 

energy identified by BCSFNM. When ILKS = -1, a source/sink of solute or energy is applied to a node 

only when lake water is absent at that node. When ILKS = 0, a source/sink of solute or energy is applied 

to a node whether lake water is present or absent at that node. When ILKS = 1, a source/sink of solute 

or energy is applied to a node only when lake water is present at that node.

DATASET 4: Boundary-condition interaction specifications for specified pressures (one line for each of NBCIP specifications)

Variable Type Description

BCSFNM Character On the first line, set BCSFNM = ‘DEFAULT’. SUTRA applies DEFAULT specifications to all specifiedpressure nodes that can form lakes but are not identified as specified-pressure nodes in the optional list 

of “.bcs” files that follows.

On each subsequent line (optional), set BCSFNM to the name of a “.bcs” file that contains specified 

pressures with which lakes will interact.

ILKP Integer Effect that the presence or absence of lake water has on the application of specified pressures identified 

by BCSFNM. When ILKP = -1, a specified pressure is applied to a node only when lake water is absent 

at that node. When ILKP = 0, a specified pressure is applied to a node whether lake water is present or 

absent at that node. When ILKP = 1, a specified pressure is applied to a node only when lake water is 

present at that node.

DATASET 5: Boundary-condition interaction specifications for specified concentrations/temperatures (one line for each of 

NBCIU specifications)

Variable Type Description

BCSFNM Character On the first line, set BCSFNM = ‘DEFAULT’. SUTRA applies DEFAULT specifications to all specifiedconcentration/temperature nodes that can form lakes but are not identified as specified-concentration/

temperature nodes in the optional list of “.bcs” files that follows.

On each subsequent line (optional), set BCSFNM to the name of a “.bcs” file that contains specified concentrations/temperatures with which lakes will interact.

ILKU Integer Effect that the presence or absence of lake water has on the application of specified concentrations or temperatures identified by BCSFNM. When ILKU = -1, a specified concentration or temperature is applied 

to a node only when lake water is absent at that node. When ILKU = 0, a specified concentration or 

temperature is applied to a node whether lake water is present or absent at that node. When ILKU = 1, a 

specified concentration or temperature is applied to a node only when lake water is present at that node.
Appendix 3. Input Data List 61

DATASET 6A: Boundary-condition interaction specifications for generalized-flow conditions (one line for each of NBCIPG 

specifications)

Variable Type Description

BCSFNM Character On the first line, set BCSFNM = ‘DEFAULT’. SUTRA applies DEFAULT specifications to all 

generalized-flow nodes that can form lakes but are not identified as generalized-flow nodes in 

the optional list of “.bcs” files that follows.

On each subsequent line (optional), set BCSFNM to the name of a “.bcs” file that contains 

generalized-flow conditions with which lakes will interact.

ILKPG Integer Effect that the presence or absence of lake water has on the application of generalized-flow conditions identified by BCSFNM. When ILKPG = -1, a generalized-flow condition is applied 

to a node only when lake water is absent at that node. When ILKPG = 0, a generalized-flow 

condition is applied to a node whether lake water is present or absent at that node. When 

ILKPG = 1, a generalized-flow condition is applied to a node only when lake water is present 

at that node.

CTIPG Character Type of interaction between the generalized-flow condition and lakes. Select one of the 

following:

‘F’ = Interaction analogous to that between a fluid source/sink node and a lake

‘P’ = Interaction analogous to that between a specified-pressure node and a lake

(See section 3.2.3.3, “Interaction with Generalized-Flow Conditions.”)

DATASET 6B: Boundary-condition interaction specifications for generalized-transport conditions (one line for each of 

NBCIUG specifications)

Variable Type Description

BCSFNM Character On the first line, set BCSFNM = ‘DEFAULT’. SUTRA applies DEFAULT specifications to all generalized-transport nodes that can form lakes but are not identified as generalized-transport nodes in 

the optional list of “.bcs” files that follows.

On each subsequent line (optional), set BCSFNM to the name of a “.bcs” file that contains generalized-transport conditions with which lakes will interact.

ILKUG Integer Effect that the presence or absence of lake water has on the application of generalized-transport conditions identified by BCSFNM. When ILKUG = -1, a generalized-transport condition is applied to 

a node only when lake water is absent at that node. When ILKUG = 0, a generalized-transport condition is applied to a node whether lake water is present or absent at that node. When ILKUG = 1, a 

generalized-transport condition is applied to a node only when lake water is present at that node.

CTIUG Character Type of interaction between the generalized-flow condition and lakes. Select one of the following:

‘S’ = Interaction analogous to that between a solute or energy source/sink node and a lake

‘U’ = Interaction analogous to that between a specified-concentration or specified-temperature node 

and a lake

(See section 3.2.3.7, “Interaction with Generalized-Transport Conditions.”)
62 SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport

Input Data for the Lake-Area Input File for the Lake Capability (.lkar)

Model Version: SUTRA 3.0

NOTE: This is an entirely new input file.

This input file is optional. If the lake capability is used and this input file is omitted, all nodes on the top surface of the 

model are eligible to form lakes and lake-bottom elevations are based on the top-surface node elevations of the groundwater model.

DATASET 1A: Number of lake nodes and method of determining lake-bottom elevations (one line)

Variable Type Description

Character Line must begin with the word ‘LAKE’.

NLAN Integer Exact number of nodes on the top surface of the model that are eligible to form lakes.

CBOT Character One word. Set to ‘DEFAULT’ to have SUTRA use the top of the groundwater model to set lake-bottom 

elevations. Set to ‘SPECIFIED’ to have lake-bottom elevations be user-specified in dataset 1B.

DATASET 1B: List of lake nodes (one line for each of NLAN nodes, plus one line)

Variable Type Description

Lines 1 to NLAN

IL Integer Number of node that is eligible to form lakes. Node IL must be on the top surface of the model.

ELVLB Real User-specified lake-bottom elevation at node IL. May be omitted if lake-bottom elevations are set by default 

(CBOT = ‘DEFAULT’).

Last line

Integer Placed immediately following all NLAN lake-node lines. Line must begin with the integer 0.
Publishing support provided by the U.S. Geological Survey Science

Publishing Network, Reston Publishing Service Center

For more information, please contact:

Earth System Processes Division

U.S. Geological Survey

Mail Stop 411

12201 Sunrise Valley Drive

Reston, VA 20192

703–648–5005
Provost and Voss—SUTRA, a Model for Saturated-Unsaturated, Variable-Density Groundwater Flow with Solute or Energy Transport—TM 6–A52

ISSN 2328-7055 (online)

https://doi.org/10.3133/tm6A52