SUTRA— A code for simulation of saturated-unsaturated, variable-density groundwater flow with solute or energy transport—Documentation of the version 4.0 enhancements—Freeze-thaw capability, saturation and relative-permeability relations, spatially varying properties, and enhanced budget and velocity outputs

Clifford I. Voss; Alden M. Provost; Jeffrey M. McKenzie; Barret L. Kurylyk

doi:10.3133/tm6A63

SUTRA: A Code for Simulation of Saturated-Unsaturated, Variable-Density Groundwater Flow With Solute or Energy Transport—Documentation of the Version 4.0 Enhancements—Freeze-Thaw Capability, Saturation and Relative-Permeability Relations, Spatially Varying Properties, and Enhanced Budget and Velocity Outputs

Techniques and Methods 6-A63

By: Clifford I. Voss, Alden M. Provost, Jeffrey M. McKenzie, and Barret L. Kurylyk

https://doi.org/10.3133/tm6A63

Metrics

Web analytics dashboard Metrics definitions

Links

Document: Report (3.80 MB pdf) , HTML , XML
Software Release: USGS software release - SUTRA—A model for 2D or 3D saturated-unsaturated, variable-density groundwater flow with solute or energy transport
Download citation as: RIS | Dublin Core

Acknowledgments

This enhancement to the U.S. Geological Survey’s Saturated-Unsaturated Transport (SUTRA) code was partly supported by the U.S. Department of Defense Strategic Environmental Research and Development Program (SERDP). SERDP supports a broad spectrum of basic and applied research and advanced development in environmental science and technology, on projects that are planned and executed in partnership with the U.S. Department of Energy and the U.S. Environmental Protection Agency. The authors much appreciate the suggestions of technical reviewers, Audrey Sawyer (Ohio State University) and Samuel Zipper (Kansas Geological Survey, University of Kansas), which have helped the authors to improve this document.

Abstract

Version 4.0 of the Saturated-Unsaturated Transport (SUTRA) software code provides the capability to simulate the freezing and thawing of groundwater during energy transport simulations under saturated and unsaturated conditions. In addition to the types of hydrogeologic processes that SUTRA has been able to simulate in the past, this version can be used to study the effects of the freeze-thaw process on the flow and energy dynamics of hydrogeologic systems. The freeze-thaw simulation capability accounts for the latent heat of fusion and allows thermal property values to vary with changing total-water saturation, liquid-water saturation, and ice saturation. It allows the effective permeability of the porous medium to change as a result of freezing and thawing. This version also provides several user-selectable relations for the dependence of total-water saturation on fluid pressure, the dependence of liquid-water saturation on temperature during freezing and thawing, and the dependence of relative permeability on liquid saturation, as well as three user-selectable formulae for defining the bulk thermal conductivity of a mixture of solid grains, liquid water, ice, and air. For unsaturated simulations without freezing, the selectable total-water saturation relations eliminate the need for the user to program these and their associated relative-permeability functions, as had been required in previous SUTRA versions. Optional nonlinear dependence of fluid density on temperature, which covers the range from supercooled (about −50 degrees Celsius) to superheated (about 400 degrees Celsius), is also provided.

Additionally, this version makes it possible to spatially vary parameters that, in previous versions of SUTRA, were required to be spatially uniform: solid-matrix properties, adsorption parameters, and parameters for production of solute mass or energy. Spatial variation is also allowed for the newly included freeze-thaw process parameters. Additional enhancements provide (1) output of water-mass and energy budgets that include values of all component terms in the governing balance equations, and (2) output of Darcy velocities (fluid fluxes), in addition to the velocity output provided by previous SUTRA versions. These enhanced outputs allow fuller interpretation of simulation results, especially for freeze-thaw phenomena.

The set of processes simulated by this version of SUTRA are useful for studying a wide range of hydrogeologic system types, conditions, and questions. For cryohydrogeologic simulations, however, this version of the code is limited in that (1) it does not simulate thermomechanical effects of freeze-thaw, (2) pressure changes due to water density change during freezing are neglected, (3) ice saturation cannot exceed the initial porosity of the simulated medium, and (4) cryosuction, the migration of liquid water toward freezing fronts, is neglected. Furthermore, this version does not account for air flow or for water vaporization and sublimation under unsaturated conditions.

Introduction

The Saturated-Unsaturated Transport (SUTRA) version 4.0 enhancements presented in this documentation allow representation of the freezing and thawing of groundwater during energy transport simulations under saturated and unsaturated conditions. All previous versions of SUTRA and the current version (code, documentation, and user-support software) may be found in U.S. Geological Survey (2022). In addition to the processes simulated by the most recent previous versions of SUTRA (Provost and Voss, 2019; Voss and Provost, 2002), this version supports simulations of (1) the effects of freezing and thawing on the flow and energy dynamics of hydrogeologic systems, ranging from small field sites to large aquifer systems and (2) hydrologic and thermal phenomena at time scales ranging from short (for example, daily temperature or laboratory experiment freeze-thaw cycles) to intermediate (for example, seasonal freeze-thaw cycles) to long (for example, freeze-thaw during climate oscillations lasting tens of thousands of years). This version can be applied to the study of groundwater flow and the potential effects of climate change in regions with permafrost or subject to seasonal ground freezing. Earlier versions of this enhanced code have been used in studies of saturated and unsaturated freeze-thaw hydrology by McKenzie and others (2006, 2007), Ge and others (2011), McKenzie and Voss (2013), Wellman and others (2013), Kurylyk and others (2014a, 2016), Briggs and others (2014), Evans and others (2018, 2020), Lamontagne-Hallé and others (2018), Walvoord and others (2019), Rey and others (2020), Rushlow and others (2020), Chen and others (2021, 2023), and Miller and others (2022).

This current version of SUTRA further extends the SUTRA saturated porous medium freeze-thaw capability developed by McKenzie and others (2007) to allow simulation of the freeze-thaw process under both saturated and unsaturated conditions, including some generalizations of the thermal energy balance equation. Development and testing versions of SUTRA leading to the current version successfully reproduce analytical solution results that consider heat transfer through conduction, advection, and latent heat released or absorbed during pore-water phase change (McKenzie and others, 2007; Kurylyk and others, 2014b) and have been successfully compared with other codes that simulate similar saturated freeze-thaw dynamics for some benchmark problems (Rühaak and others, 2015; Grenier and others, 2018).

The freeze-thaw capability of SUTRA accounts for the large effect of the latent heat of fusion on the subsurface energy balance during ice-liquid phase change. It allows the effective permeability of the porous medium to decrease when pore water freezes and to increase when pore ice melts. This version includes user-selectable functions for the dependence of total-water saturation on fluid pressure, the dependence of liquid-water and ice saturations on temperature, and the dependence of relative permeability on liquid-water saturation. The user-selectable total-water-saturation and relative-permeability functions may also be used for simulations of unsaturated conditions without freezing, providing an enhancement to the unsaturated-zone simulation capability of SUTRA; previous SUTRA versions required user programming to define these functions. In addition, this version includes user-selectable formulae for the bulk thermal conductivity of a mixture of solid grains, liquid water, ice, and air. Also provided is nonlinear dependence of fluid density on temperature (in addition to the linear dependence included in previous versions), covering the range from freezing temperatures (about −50 degrees Celsius [°C]) to superheated fluid temperatures (about 400 °C), with maximum fluid density occurring at 4 °C.

Further, this version allows for spatially variable properties that, in previous versions of SUTRA, had spatially uniform values—solid-matrix properties, adsorption parameters, and parameters for production of solute mass or energy. Spatial variation is also allowed for the new freeze-thaw process parameters. Additional enhancements provide (1) output of water-mass and energy budgets that include values of all component terms in the governing balance equations, and (2) output of Darcy velocities of liquid water (fluid fluxes), in addition to the average linear fluid velocity output provided by previous SUTRA versions. These enhanced outputs allow fuller interpretation of simulation results, especially for freeze-thaw phenomena.

This version of SUTRA is not intended for simulation of systems that involve significant cryosuction or volumetric changes in the porous medium such as those during frost heave (see for example, Rempel, 2012). Furthermore, this version does not represent ice-rich media, in which ice is the continuous medium within which mineral grains are dispersed; total porosity has a fixed value during a simulation, and the ice content never exceeds the total porosity. The dynamics of unsaturated zone air and water-vapor flow, the ice-vapor sublimation process, and dynamic concentration effects on groundwater energy are not considered in this version.

1. Freeze-Thaw Capability

The SUTRA freeze-thaw capability can be optionally activated by the user for energy transport simulations, allowing groundwater under saturated or unsaturated conditions to gradually freeze as temperature decreases below the user-selected freezing temperature, and allowing ice in the pores to gradually melt as subfreezing temperatures increase. This enhancement accounts for the latent heat of fusion, changes in relative permeability and mobile porosity (the connected volume fraction of the porous medium containing liquid water through which water can flow) due to blockage of pores by ice, and changes in bulk thermal properties of the medium due to freeze-thaw.

Groundwater flow is driven by pressure boundary conditions, fluid sources and sinks, and fluid-density differences resulting from temperature differences, as it has been for energy-transport simulations in all previous versions of SUTRA. Also, as in previous versions of SUTRA, the temperature scale is Celsius. The set of processes simulated by this version of SUTRA is expected to be useful for studying the wide range of cryohydrogeologic systems in which thermomechanical effects can be neglected (such as most systems described by Walvoord and Kurylyk, 2016).

The migration of liquid water toward freezing fronts (cryosuction), the volumetric expansion of water when it freezes, and volumetric changes in the porous medium due to freezing (such as during frost heave) cannot be simulated in this version of SUTRA. The porous medium is assumed to have a fixed volume and a time-independent total porosity, and the ice volume cannot exceed the volume of pores.

1.1. Water-Mass Balance With Freeze-Thaw

Here, the development of the total water-mass balance (the groundwater-flow equation) Voss and Provost (2002, secs. 2.1 and 2.2) is extended to include freezing. Two types of saturation must be distinguished—liquid-water saturation (S_L) and frozen-water (ice) saturation (S_I). These sum to equal the total-water saturation (S_W) of the porous medium:

S_W

=

S_L

+

S_I

,

(1)

where

S_W(x,y[,z],t): [1] is total-water saturation (volume of frozen plus unfrozen water per volume of porous medium voids)
S_L(x,y[,z],t): [1] is liquid-water saturation (volume of unfrozen water per volume of voids)
S_I(x,y[,z],t): [1] is ice saturation (volume of ice per volume of voids).

Here, (x, y[, z], t) indicates the possibility of two-dimensional (x, y) and three-dimensional (x, y, z) spatial representations, with “t” indicating that these parameters are time-dependent. The units of the parameter appear in brackets; for the purposes of this documentation, “[1]” indicates a dimensionless quantity, “[L]” indicates length units, “[M]” indicates mass units, “[E]” indicates energy units, and “[s]” indicates time units in seconds. A complete list of units, symbols, and abbreviations used in this documentation of SUTRA version 4.0 is provided in appendix 1.

The term “water” is ambiguous when describing systems in which phase change occurs. It may be understood to imply all water (H₂O; including impurities such as solutes) irrespective of phase, or it may be understood to imply only the liquid phase. Thus, care must be taken when referring to all the water and to each of the water phases. Previous versions of SUTRA simulated variable-density, saturated-unsaturated groundwater flow and solute or energy transport without freezing. Therefore, in previous versions, the terms “water,” “liquid,” and “fluid” all meant liquid water. In this document, the terms “liquid” and “fluid” are used to refer to the liquid-water phase (including any dissolved species), and “ice” refers to the frozen-water phase. The terms “total water” and “water” refer to all the H₂O; they include both the liquid-water phase (and any dissolved species) and the ice phase.

Note: Parameters are defined as referring to “Total Water,” “Liquid Water” and “Ice” by including the letter “W,” “L,” or “I” for each parameter as a subscript or superscript in this document, and as part of the variable name in the computer code. Depending on the context, the letter “S” as a subscript or superscript in this document, or as part of a code parameter, refers to either the porous medium matrix of “Solid” mineral grains or to the chemical species being tracked by SUTRA, wherever it exists in the porous medium (that is, to “Solute,” the dissolved portion of the chemical species, or to “Sorbate,” the portion of the chemical species that is sorbed onto and within the solid mineral grains). Generic units (app. 1) that pertain to the porous medium matrix or to individual solid mineral grains are subscripted with the letter “G” to distinguish them from generic units that pertain to chemical species, which are subscripted with the letter “S.”

Note: For simulations of nonfreezing conditions, there is no ice, and the total-water and liquid parameters become equivalent to each other and to the same parameters that existed in previous versions of SUTRA. In this case, SUTRA simulates exactly the same processes as in previous versions of the code. The nonfreezing processes are backward compatible, so this version of SUTRA includes all types of simulation provided by previous versions.

The mean total-water density (ρ_W) is defined by

ρ_WS_W

=

ρ_LS_L

+

ρ_IS_I

,

(2)

where

ρ_W(x,y[,z],t): [M/L³] is mean total-water density (including both liquid water and ice);
ρ_L(x,y[,z],t): [M/L³] is liquid-water density (user-specified: may be held constant, for example, at 1,000 kilograms per cubic meter [kg/m³], which is the approximate value from 0 °C to 20 °C; set to a linear function of temperature; or set to a nonlinear function of temperature, as described in section 1.5, “Nonlinear Liquid-Water Density Function;” in mass per volume); and
ρ_I(x,y[,z],t): [M/L³] is ice density (user-specified as a constant value; about 917 kg/m³ at 0 °C; about 919 kg/m³ at –10 °C).

For the case of no freezing, there is no ice, S_I=0, and it follows from equations 1 and 2 that S_W=S_L and ρ_W=ρ_L.

Note: For most parameters that require a user-specified value, a typical value or value range is usually provided where the parameter is first defined in this documentation. Users should refer to the scientific literature to find values that are appropriate to the system that they are simulating.

The mass balance for total water (frozen plus unfrozen) is unchanged from previous versions of SUTRA, and the ice and solid matrix of minerals are assumed to be immobile:

\frac{\partial (ε S_{W} ρ_{W})}{\partial t} = - \underline{\nabla} \cdot (ε S_{L} ρ_{L} \underline{v}) + Q_{p}

,

(3)

where

t: [s] is time (in seconds);
ε(x, y[, z]): [1] is porosity (volume of voids per total volume; user-specified value at each node; typical range for geologic materials from 10⁻⁵ to 0.5);
v(x, y[, z], t): [L/s] is average fluid velocity; and
Q_p(x, y[, z], t): [M/(L³‧s)] is fluid-mass source (user-specified value at each fluid-source node).

The quantity being differentiated on the left-hand side of equation 3 is the total mass of water (frozen and unfrozen) per volume of porous medium. The quantity in parentheses on the right-hand side of equation 3 is the mass flux of liquid water.

Note: The quantity “total water content,” typically used by soil scientists, is the product of porosity and total-water saturation (εS_W). This quantity is the total water volume given as a fraction of the total porous medium volume.

Substitution of equation 2 into 3 gives

\frac{\partial (ε S_{L} ρ_{L})}{\partial t} + \frac{\partial (ε S_{I} ρ_{I})}{\partial t} = - \underline{\nabla} \cdot (ε S_{L} ρ_{L} \underline{v}) + Q_{p}

.

(4)

The first term on the left-hand side of equation 4 represents the rate of change in storage of liquid water. Assuming S_L is a function of pressure and temperature, this term can be written as

\frac{\partial (ε S_{L} ρ_{L})}{\partial t} = \frac{\partial (ε S_{L} ρ_{L})}{\partial p} \frac{\partial p}{\partial t} + \frac{\partial (ε S_{L} ρ_{L})}{\partial T} \frac{\partial T}{\partial t}

,

(5)

where

p(x, y[, z], t): [M/(L s²)] is fluid pressure and
T(x, y[, z], t): [°C] is temperature (in degrees Celsius).

Note: All water phases and minerals are assumed to have the same temperature at any point in the porous medium; local thermal equilibrium is assumed.

Using the product rule of differentiation, and introducing compressibilities and specific pressure storativity as in the SUTRA flow equation presented in Voss and Provost (2002), the derivative with respect to p in the first term of the right-hand side of equation 5 can be written as

\frac{\partial (ε S_{L} ρ_{L})}{\partial p} = S_{L} ρ_{L} \frac{\partial ε}{\partial p} + ε S_{L} \frac{\partial ρ_{L}}{\partial p} + ε ρ_{L} \frac{\partial S_{L}}{\partial p} = ρ_{L} (S_{L} S_{op}^{L} + ε \frac{\partial S_{L}}{\partial p})

,

(6)

where

S_{op}^{L}

[M/(L s²)]⁻¹ is specific pressure storativity for liquid.The specific pressure storativity for liquid depends on compressibilities defined by

\frac{\partial ε}{\partial p} = (1 - ε) α_{S}

and

(7)

\frac{\partial ρ_{L}}{\partial p} = ρ_{L} β_{L}

,

(8)

where

α_S: [M/(L‧s²)]⁻¹ is porous-matrix compressibility (user-specified value at each node; α_S ranges from about 10^–10 per kilogram per meter per second squared [(kg/m‧s²)⁻¹] for sound bedrock to about [10^–7 kg/(m‧s²)]^–1 for clay) and
β_L: [M/(L‧s²)]⁻¹ is liquid-water compressibility (user-specified constant value for the entire model; about 4.47×10^–10 [kg/(m‧s²)]⁻¹ at 20 °C to about 5.09×10^–10 [kg/(m‧s²)]⁻¹ at 0 °C).

The specific pressure storativity is then defined by

S_{op}^{L} \equiv (1 - ε) α_{S} + ε β_{L}

.

(9)

Again, using the product rule, and assuming that porosity does not vary with temperature (∂ε/∂T=0), the derivative with respect to T in the second term on the right-hand side of equation 5 can be written as

\frac{\partial (ε S_{L} ρ_{L})}{\partial t} = S_{L} ρ_{L} \frac{\partial ε}{\partial T} + ε S_{L} \frac{\partial ρ_{L}}{\partial T} + ε ρ_{L} \frac{\partial S_{L}}{\partial T} = ε (S_{L} \frac{\partial ρ_{L}}{\partial T} + ρ_{L} \frac{\partial S_{L}}{\partial T})

.

(10)

Substitution of equations 6 and 10 into equation 5 then gives the following expression for the rate of change in storage of liquid water:

\frac{\partial (ε S_{L} ρ_{L})}{\partial t} = ρ_{L} (S_{L} S_{op}^{L} + ε \frac{\partial S_{L}}{\partial p}) \frac{\partial p}{\partial t} + ε (S_{L} \frac{\partial ρ_{L}}{\partial T} + ρ_{L} \frac{\partial S_{L}}{\partial T}) \frac{\partial T}{\partial t}

.

(11)

The second term on the left-hand side of equation 4 represents the rate of change in storage of ice. An expression for this term can be derived by following a development that is analogous to the foregoing development for liquid. Assuming that S_I is a function of pressure and temperature, the second term on the left-hand side of equation 4 can be written as

\frac{\partial (ε S_{I} ρ_{I})}{\partial t} = \frac{\partial (ε S_{I} ρ_{I})}{\partial p} \frac{\partial p}{\partial t} + \frac{\partial (ε S_{I} ρ_{I})}{\partial T} \frac{\partial T}{\partial t}

.

(12)

Using the product rule of differentiation and introducing the compressibility of ice and a corresponding quantity that is mathematically analogous to the specific pressure storativity, the derivative with respect to p in the first term of the right-hand side of (12) can be expanded as

\frac{\partial (ε S_{I} ρ_{I})}{\partial p} = S_{I} ρ_{I} \frac{\partial ε}{\partial p} + ε S_{I} \frac{\partial p}{\partial t} + ε ρ_{I} \frac{\partial S_{I}}{\partial p} = ρ_{I} (S_{I} S_{op}^{I} + ε \frac{\partial S_{I}}{\partial p})

,

(13)

where

S_{op}^{I}

[M/(L‧s²)]⁻¹ is a quantity mathematically analogous to the specific pressure storativity for liquid but computed using the compressibility of ice (β_I) instead of the compressibility of liquid water (β_L). Ice compressibility and

S_{op}^{I}

are defined by

\frac{\partial ρ_{I}}{\partial p} = ρ_{I} β_{I}

and

(14)

S_{op}^{I} \equiv (1 - ε) α_{S} + ε β_{I}

,

(15)

where

β_I: [M/(L‧s²)]⁻¹ is ice compressibility (user-specified constant value for the entire model; about 1.14×10^–10 kg/(m s²)^–1 at 0 °C).

Assuming again that porosity does not vary with temperature (∂ε/∂T=0), the derivative with respect to T in the second term of the right-hand side of equation 12 can be written as

\frac{\partial (ε S_{I} ρ_{I})}{\partial T} = S_{I} ρ_{I} \frac{\partial ε}{\partial T} + ε S_{I} \frac{\partial ρ_{I}}{\partial T} + ε ρ_{I} \frac{\partial S_{I}}{\partial T} = ε (S_{I} \frac{\partial ρ_{I}}{\partial T} + ρ_{I} \frac{\partial S_{I}}{\partial T})

.

(16)

Substitution of equations 13 and 16 into equation 12 then gives the following expression for the rate of change in storage of ice:

\frac{\partial (ε S_{I} ρ_{I})}{\partial t} = ρ_{I} (S_{I} S_{op}^{I} + ε \frac{\partial S_{I}}{\partial p}) \frac{\partial p}{\partial t} + ε (S_{I} \frac{\partial ρ_{I}}{\partial T} + ρ_{I} \frac{\partial S_{I}}{\partial T}) \frac{\partial T}{\partial t}

.

(17)

Finally, substitution of equations 11 and 17 and Darcy’s Law

\underline{v} = - (\frac{\underline{\underline{k}} k_{r}}{ε S_{L} μ_{L}}) (\underline{\nabla} p - ρ \underline{g})

(18)

into equation 4, and rearrangement gives

Q_{p} = [S_{W} ρ_{W} S_{op}^{eff} + ε (ρ_{L} \frac{\partial S_{L}}{\partial p} + ρ_{I} \frac{\partial S_{I}}{\partial p})] \frac{\partial p}{\partial t} + ε (S_{L} \frac{\partial ρ_{L}}{\partial T} + S_{I} \frac{\partial ρ_{I}}{\partial T} + ρ_{L} \frac{\partial S_{L}}{\partial T} + ρ_{I} \frac{\partial S_{I}}{\partial T}) \frac{\partial T}{\partial t} - \underline{\nabla} \cdot [(\frac{\underline{\underline{k}} k_{r} ρ_{L}}{μ_{L}}) (\underline{\nabla} p - ρ_{L} \underline{g})]

,

(19)

where

µ_L(x, y[, z], t): [M/(L‧s)] is liquid viscosity (dependent on temperature for energy transport, see description of function provided in Voss and Provost (2002); user-specified constant value for the entire model for solute transport; about 1×10^–3 kg/(m‧s) at 20 °C)
$\underline{\underline{k}} (x, y [, z])$: [L²] is solid matrix permeability tensor (user-specified values for principal components of this tensor for each finite element)
k_r(x, y[, z], t): [1] is relative permeability for fluid flow (assumed to be independent of direction and dependent on saturation of liquid water; preprogrammed or user-defined function user-specified for each spatial region; see sections 2.1, “Overview of Available Functions,” and 2.4, “Preprogrammed Relative-Permeability Functions”);
$\underline{g}$: [L/s²] is gravitational acceleration (gravity vector; user-specified constant value for the entire model; for earth systems about 9.81 m/s² in downward direction); and
$S_{op}^{eff} (x, y [, z], t)$: [M/(L‧s²)]^–1 is effective specific pressure storativity.

The effective specific pressure storativity ( $S_{op}^{eff}$ ), which takes into account the compressibilities of the solid matrix, liquid water, and ice, is defined by

S_{W} ρ_{W} S_{op}^{eff} \equiv S_{L} ρ_{L} S_{op}^{L} + S_{I} ρ_{I} S_{op}^{I} = (1 - ε) S_{W} ρ_{W} α_{S} + ε (S_{L} ρ_{L} β_{L} + S_{I} ρ_{I} β_{I})

.

(20)

If ice is absent, and because, in the absence of ice, ∂S_L/∂T=0 (liquid-water saturation does not depend on temperature), equation 19 reduces to the flow equation (eq. 2.24) of Voss and Provost (2002), and thus, this total water mass balance equation is backward-compatible with previous versions of SUTRA.

1.1.1. Assumption Regarding Pressure Effect of Liquid-Ice Density Difference in Water-Mass Balance

SUTRA version 4.0 is not intended for simulation of thermomechanical effects of freeze-thaw. Cryosuction is neglected, and porosity changes associated with compressibility of the solid matrix are accounted for implicitly in the specific pressure storativity, but the value of porosity is not explicitly updated, and creation of new pore space such as occurs during frost heave is not represented. Therefore, in some situations, the volumetric expansion associated with freezing could cause anomalously high pressures to be simulated. To prevent this, the volumetric expansion associated with freezing is disregarded by setting the ice density equal to the liquid-water density in particular terms of the water mass balance (eq. 19) and by setting the derivative of ice density with respect to temperature equal to that of liquid water.

The effect of volumetric expansion associated with freezing caused by temperature change enters equation 19 through the terms that involve derivatives of saturations with respect to temperature: $ρ_{L} \frac{\partial S_{L}}{\partial T} + ρ_{I} \frac{\partial S_{I}}{\partial T}$ . Using equation 1 to express liquid-water saturation as S_L=S_W−S_I, these terms can be written as

ρ_{L} \frac{\partial S_{L}}{\partial T} + ρ_{I} \frac{\partial S_{I}}{\partial T} = ρ_{L} \frac{\partial S_{W}}{\partial T} + (ρ_{I} - ρ_{L}) \frac{\partial S_{I}}{\partial T}

.

(21)

Setting ρ_I=ρ_L in the term

ρ_{I} \frac{\partial S_{I}}{\partial T}

in equation 19 zeroes the second term on the right-hand side of equation 21, that is, the term that represents the effect of volumetric expansion associated with freezing caused by temperature change.

Analogously, the effect of volumetric expansion associated with freezing caused by pressure change enters equation 19 through the terms that involve derivatives of saturations with respect to pressure: $ρ_{L} \frac{\partial S_{L}}{\partial p} + ρ_{I} \frac{\partial S_{I}}{\partial p}$ . Again, using equation 1 to express liquid-water saturation as S_L=S_W−S_I, these terms can be written as

ρ_{L} \frac{\partial S_{L}}{\partial p} + ρ_{I} \frac{\partial S_{I}}{\partial p} = ρ_{L} \frac{\partial S_{W}}{\partial p} + (ρ_{I} - ρ_{L}) \frac{\partial S_{I}}{\partial p}

.

(22)

Setting ρ_I=ρ_L in the term

ρ_{I} \frac{\partial S_{I}}{\partial p}

in equation 19 zeroes the second term on the right-hand side of equation 22, that is, the term that represents the effect of volumetric expansion associated with freezing caused by pressure change.

With ρ_I=ρ_L in equations 21 and 22 and different ice and liquid density values in all other terms, substitution of equations 21 and 22 into equation 19, together with the assumption that total saturation (S_W) depends only on pressure such that $\frac{\partial S_{W}}{\partial T}$ =0 (as discussed later regarding equation 37 in section 1.3.2, “Freeze-Thaw Under Unsaturated Conditions”), and with the assumption that the ice density temperature derivative is equal to that of liquid water $(\frac{\partial ρ_{I}}{\partial T} = \frac{\partial ρ_{L}}{\partial T})$ , results in the form of the groundwater-flow equation that is solved by SUTRA version 4.0:

Q_{p} = (S_{W} ρ_{W} S_{op}^{eff} + ε ρ_{L} \frac{\partial S_{W}}{\partial p}) \frac{\partial p}{\partial t} + ε S_{W} \frac{\partial ρ_{L}}{\partial T} \frac{\partial T}{\partial t} - \underline{\nabla} \cdot [(\frac{\underline{\underline{k}} k_{r} ρ_{L}}{μ_{L}}) (\underline{\nabla} p - ρ_{L} \underline{g})]

.

(23)

1.2. Energy Balance With Freeze-Thaw

Here, the development of the energy-balance equation in Voss and Provost (2002, section 2.3, “Energy Transport in Ground Water”) is generalized to include the freeze-thaw process. The balance of total thermal energy for liquid water, ice, and the solid matrix is

\frac{\partial}{\partial t} [ε (S_{L} ρ_{L} e_{L} + S_{I} ρ_{I} e_{I}) + (1 - ε) ρ_{S} e_{S}] = - \underline{\nabla} \cdot (ε S_{L} ρ_{L} e_{L} \underline{v}) - \underline{\nabla} \cdot {\underline{q}}_{e} + Q_{p} e_{L}^{*} + ε (S_{L} ρ_{L} γ_{0}^{L} + S_{I} ρ_{I} γ_{0}^{I}) + (1 - ε) ρ_{S} γ_{0}^{S}

,

(24)

where

e_L(x, y[, z], t): [E/M] is energy per unit mass of liquid water;
e_I(x, y[, z], t): [E/M] is energy per unit mass of ice;
e_S(x, y[, z], t): [E/M] is energy per unit mass of solid mineral grains;
$e_{L}^{*} (x, y [, z], t)$: [E/M] is energy per unit mass of fluid source;
ρ_S(x, y[, z], t): [M/L³] is density of solid grains in the mineral matrix (user-specified value at each node; about 2,600 kg/m³ for typical mineral grains);
${\underline{q}}_{e} (x, y [, z], t)$: [E/(L²‧s)] is conductive-dispersive heat flux;
$Q_{p} e_{L}^{*}$: [E/(M‧s)] is the rate at which thermal energy is added by an external fluid source to a point in the model domain;
$γ_{0}^{L} (x, y [, z], t)$: [E/(M‧s)] is a distributed energy source in liquid (user-specified value at each node);
$γ_{0}^{I} (x, y [, z], t)$: [E/(M‧s)] is a distributed energy source in ice (user-specified value at each node); and
$γ_{0}^{S} (x, y [, z], t)$: [E/(M‧s)] is a distributed energy source in solid mineral grains (user-specified value at each node).

The term on the left-hand side of equation 24 represents the rate of change of total thermal energy in the porous medium, including energy in liquid, ice, and mineral components. The first term on the right-hand side of equation 24 represents the (negative) divergence of advective thermal energy transport (local influx of thermal energy due to groundwater flow) at a point in the modeled spatial domain. The second term represents the rate at which thermal energy enters a point in the model domain by conduction and mechanical dispersion. The last terms are distributed energy sources, which add energy to points in the model domain. These distributed energy sources may, for example, be due to radioactive decay, chemical reactions, or bacterial metabolism. These processes are not explicitly represented in SUTRA, but their total energy contribution can be included in the energy balance by user specification of values for $γ_{0}^{L}$ , $γ_{0}^{I}$ , and $γ_{0}^{S}$ .

Using the product rule of differentiation for the terms on the left-hand side of the equation, equation 24 can be written as

\begin{array}{l} ε S_{L} ρ_{L} \frac{\partial e_{L}}{\partial t} + ε S_{I} ρ_{I} \frac{\partial e_{I}}{\partial t} + (1 - ε) ρ_{S} \frac{\partial e_{S}}{\partial t} + e_{L} \frac{\partial}{\partial t} (ε S_{L} ρ_{L}) + e_{I} \frac{\partial}{\partial t} (ε S_{I} ρ_{I}) + e_{S} \frac{\partial}{\partial t} [(1 - ε) ρ_{S}] = \\ - \underline{\nabla} \cdot (ε S_{L} ρ_{L} e_{L} \underline{v}) - \underline{\nabla} \cdot {\underline{q}}_{e} + Q_{p} e_{L}^{*} + ε (S_{L} ρ_{L} γ_{0}^{L} + S_{I} ρ_{I} γ_{0}^{I}) + (1 - ε) ρ_{S} γ_{0}^{S} \end{array}

.

(25)

Definition and specification of some terms in the above equation are required to obtain the form of the energy balance simulated by SUTRA, as described in the following paragraphs.

The thermal-energy contents of liquid water, ice, and solid matrix, respectively, relative to the values they have at the freezing temperature of free water (T_ff) are here defined as

e_{L} \equiv e_{L}^{ff} + c_{L} (T - T_{ff})

,

(26A)

e_{I} \equiv e_{I}^{ff} + c_{I} (T - T_{ff})

, and

(26B)

e_{S} \equiv e_{S}^{ff} + c_{S} (T - T_{ff})

,

(26C)

where

T_ff: [°C] is freezing temperature of pure, free water (solute-free water not confined to pores); T_ff is set to 0 °C in SUTRA;
$e_{L}^{ff}$: [E/M] is thermal energy per unit mass of liquid water at free-water freezing temperature, T=T_ff;
$e_{I}^{ff}$: [E/M] is thermal energy per unit mass of ice at free-water freezing temperature, T=T_ff;
$e_{S}^{ff}$: [E/M] is thermal energy per unit mass of solid mineral grains at free-water freezing temperature, T=T_ff;
c_L(x, y[, z], t): [E/(M‧°C)] is specific heat of liquid (user-specified constant value for the entire model; c_L is about 4.182×10³ J/(kg‧°C) at 20 °C and about 4.218×10³ J/(kg∙°C) at 0 °C);
c_I(x, y[, z], t): [E/(M °C)] is specific heat of ice (user-specified constant value for the entire model; c_I is about 2.115×10³ J/(kg‧°C) at 0°C, about 2.037×10³ J/(kg‧°C) at –10 °C, and about 1.959×10³ J/(kg‧°C) at −20 °C); and
c_S(x, y[, z], t): [E/(M‧°C)] is specific heat of solid mineral grains (user-specified value at each node; c_S about 8.4×10² J/(kg‧°C) for quartz or calcite; typical range: 600 to 900 J/(kg‧°C) for other minerals).

The freezing temperature for pure, free water (T_ff), which is set to a value of 0 °C in SUTRA, serves as a reference temperature for a formal definition of energy content. However, the actual freezing temperature of water may be lower (for example, owing to the presence of solute and confinement to pores). To allow SUTRA to account, at least approximately, for freezing-point depression, the actual freezing temperature, T_f, (the maximum temperature at which ice can exist) can be assigned values different from 0 °C:

T_f(x, y[, z]): [°C] maximum freezing temperature of pore water (user-specified constant value in each spatial region).

Note: T_f is usually set to 0 °C for water containing only dilute solutes and for a relatively coarse mineral matrix. For groundwater with high solute concentration, for groundwater in fine-grained porous media (with large surface-energy effects), or for both, the maximum freezing temperature may be set to a different value.

The difference between the thermal-energy contents of liquid water and ice when both are at the free-water freezing temperature (T=T_ff) is the latent heat of fusion. Subtracting equation 26B from equation 26A yields

{(e_{L} - e_{I})}_{T = T_{ff}} = (e_{L}^{ff} - e_{I}^{ff}) \equiv Δ H

,

(27)

where

∆H(x, y[, z]): [E/M] is latent heat of fusion of free water (∆H>0); usually ∆H is about +3.34×10⁵ J/kg when free water has zero solute concentration (in energy per mass units); but may be user-specified to be different from this value, to represent latent heat of fusion change (usually a decrease) resulting from presence of solutes in the groundwater. (User-specified value for each spatial region.)

Note: The dependence of thermal energy content on solute concentration causes freezing-point depression and dependence of the latent heat of fusion on solute concentration. However, in equation 26A–C, it is assumed that thermal energy content depends only on temperature; the effects of pressure, solute concentration, and surface energy are neglected. As mentioned earlier, the user-specified freezing temperature (T_f) can be used to approximate the effects of freezing-point depression, and the user-specified latent heat of fusion (∆H) can be used to account approximately for changes in the latent heat due to the presence of solute. SUTRA considers both user-specified parameters to be constant in time, but the user may vary them spatially by region within the modeled area or volume to represent the effect of, for example, different salinities in different aquifers or subsurface zones. As the current version of SUTRA does not simultaneously simulate energy transport and solute transport, these solute concentration effects must be considered to be constant with time and are implied by the user’s choices for the spatial distributions of latent heat and freezing temperature.

For SUTRA, the energy contents of liquid and solid grains at the freezing temperature of free water (T_ff) are defined as zero: $e_{L}^{ff} \equiv 0$ and $e_{S}^{ff} \equiv 0$ . According to equation 27, when liquid and ice are both at the freezing temperature, T_ff=0 °C, the difference in their energy is exactly the latent heat ∆H, with the liquid having the higher energy content. Thus, $e_{I}^{ff} = - Δ H$ , and the energy definitions (eq. 26A–C) are simplified to

e_L

≡

c_LT

,

(28A)

e_I

≡ −∆

H

+

c_IT

, and

(28B)

e_S

≡

c_ST

,

(28C)

where T is the temperature, in degrees Celsius. The c_[x]T terms in each energy content equation are referred to as the “sensible heat,” the energy associated with temperature change. This “sensible heat” is included in the equation for ice (eq. 28B), but the energy content of ice also includes another quantity, the latent heat of fusion (∆H), the amount of energy released from water when it freezes. On the basis of the above formal definitions of energy content, SUTRA assigns an energy content of zero to liquid water and solid mineral grains when they have a temperature of 0 °C. At the same temperature, 0 °C, ice has an energy content equal to the negative latent heat of fusion (−∆H).

Following Voss and Provost (2002), the solid matrix is assumed to be immobile, so the mass balance for the solid matrix consisting of mineral grains is

\frac{\partial [(1 - ε) ρ_{S}]}{\partial t} = 0

.

(29)

Then noting that ∆H, c_L, c_I, and c_S are constant in time, substitution of equation 28A–C and 29 into equation 25 and rearrangement gives

\begin{array}{l} [ε (S_{L} ρ_{L} c_{L} + S_{I} ρ_{I} c_{I}) + (1 - ε) ρ_{S} c_{S}] \frac{\partial T}{\partial t} + c_{L} T \frac{\partial}{\partial t} (ε S_{L} ρ_{L}) + (- Δ H + c_{I} T) \frac{\partial}{\partial t} (ε S_{I} ρ_{I}) = \\ - \underline{\nabla} \cdot (ε S_{L} ρ_{L} c_{L} \underline{v} T) - \underline{\nabla} \cdot {\underline{q}}_{e} + Q_{p} c_{L} T^{*} + ε (S_{L} ρ_{L} γ_{0}^{L} + S_{I} ρ_{I} γ_{0}^{I}) + (1 - ε) ρ_{S} γ_{0}^{S} \end{array}

,

(30)

where

T*(x, y[, z], t): [°C] is temperature of the fluid source (user-specified value at each fluid-source node).

The conductive-dispersive heat flux is expressed in terms of an effective, or “bulk,” thermal conductivity, which can take into account contributions from the liquid water, ice, solid matrix, and air, and a mechanical dispersion tensor:

{\underline{q}}_{e} = - (λ \underline{\underline{I}} + ε S_{L} ρ_{L} c_{L} \underline{\underline{D}}) \underline{\nabla} T

,

(31)

where

λ(x, y[, z], t): [E/(s‧L‧°C)] is bulk thermal conductivity (depends on thermal conductivity values of components specified by the user as described in section 1.4, “Thermal Conductivity Functions”), and
$\underline{\underline{D}} (x, y [, z], t)$: [L²/s] is the mechanical dispersion tensor.

Equation 30 is the so-called “divergent” form of the energy balance. However, SUTRA solves the so-called “convective” form of the energy balance, obtained by expanding the divergence term that represents the contribution of fluid advection to the energy balance (the term that contains velocity) on the right-hand side of equation 30, then subtracting the water mass balance (eq. 4) multiplied by c_LT from both sides of equation 30. Using the product rule for divergence, and noting that c_L is assumed to be spatially uniform, the divergence term can be written as

\underline{\nabla} \cdot [ε S_{L} ρ_{L} c_{L} \underline{v} T] = ε S_{L} ρ_{L} c_{L} \underline{v} \cdot \underline{\nabla} T + c_{L} T \underline{\nabla} \cdot (ε S_{L} ρ_{L} \underline{v})

.

(32)

Substitution of equations 31 and 32 into equation 30, and subtraction of c_LT times (eq. 4) from both sides of the resulting equation gives

\begin{array}{l} [ε (S_{L} ρ_{L} c_{L} + S_{I} ρ_{I} c_{I}) + (1 - ε) ρ_{S} c_{S}] \frac{\partial T}{\partial t} + [(c_{I} - c_{L}) T - Δ H] \frac{\partial}{\partial t} (ε S_{I} ρ_{I}) = \\ - ε S_{L} ρ_{L} c_{L} \underline{v} \cdot \underline{\nabla} T + \underline{\nabla} \cdot [(λ \underline{\underline{I}} + ε S_{L} ρ_{L} c_{L} \underline{\underline{D}}) \cdot \underline{\nabla} T] + Q_{p} c_{L} (T^{*} - T) + ε (S_{L} ρ_{L} γ_{0}^{L} + S_{I} ρ_{I} γ_{0}^{I}) + (1 - ε) ρ_{S} γ_{0}^{S} \end{array}

.

(33)

Expansion of

\frac{\partial}{\partial t} (ε S_{I} ρ_{I})

using equation 17 and rearrangement gives

\begin{array}{l} [ε (S_{L} ρ_{L} c_{L} + S_{I} ρ_{I} c_{I}) + (1 - ε) ρ_{S} c_{S}] \frac{\partial T}{\partial t} + [(c_{I} - c_{L}) T - Δ H] [(ρ_{I} S_{I} S_{op}^{I} + ε ρ_{I} \frac{\partial S_{I}}{\partial p}) \frac{\partial p}{\partial t} + ε (S_{I} \frac{\partial ρ_{I}}{\partial T} + ρ_{I} \frac{\partial S_{I}}{\partial T}) \frac{\partial T}{\partial t}] = \\ - ε S_{L} ρ_{L} c_{L} \underline{v} \cdot \underline{\nabla} T + \underline{\nabla} \cdot [(λ \underline{\underline{I}} + ε S_{L} ρ_{L} c_{L} \underline{\underline{D}}) \cdot \underline{\nabla} T] + Q_{p} c_{L} (T^{*} - T) + ε (S_{L} ρ_{L} γ_{0}^{L} + S_{I} ρ_{I} γ_{0}^{I}) + (1 - ε) ρ_{S} γ_{0}^{S} \end{array}

.

(34)

If ice is absent, equation 34 becomes equivalent to the energy-transport variant of equation 2.52 of Voss and Provost (2002), and thus, this energy balance equation is backward-compatible with previous versions of SUTRA.

1.2.1. Assumption Regarding Pressure Effect of Liquid-Ice Density Difference in Energy Balance

As mentioned in section 1.1.1, “Assumption Regarding Pressure Effect of Liquid-Ice Density Difference in Water-Mass Balance” (see discussion related to equation 22), the volumetric expansion associated with freezing is removed by setting the ice density, ρ_I, equal to the liquid-water density (ρ_L) in two terms, $ρ_{I} \frac{\partial S_{I}}{\partial p}$ and $ρ_{I} \frac{\partial S_{I}}{\partial T}$ , of the total water mass balance (eq. 19). For consistency with the water-mass-balance equation, the value of ρ_I is also set equal to ρ_L for exactly the same two terms, where these appear in the energy balance (eq. 34). Also, as in the total water balance, it is assumed that ice density dependence on temperature, where it appears in the term $ε (S_{I} \frac{\partial ρ_{I}}{\partial T}) \frac{\partial T}{\partial t}$ , is equal to the liquid water temperature dependence: $\frac{\partial ρ_{I}}{\partial T} = \frac{\partial ρ_{L}}{\partial T}$ in this version of SUTRA. This results in the form of the energy transport equation that is solved by this version of SUTRA (version 4.0):

\begin{array}{l} [ε (S_{L} ρ_{L} c_{L} + S_{I} ρ_{I} c_{I}) + (1 - ε) ρ_{S} c_{S}] \frac{\partial T}{\partial t} + [(c_{I} - c_{L}) T - Δ H] [(ρ_{I} S_{I} S_{op}^{I} + ε ρ_{L} \frac{\partial S_{I}}{\partial p}) \frac{\partial p}{\partial t} + ε (S_{I} \frac{\partial ρ_{L}}{\partial T} + ρ_{I} \frac{\partial S_{I}}{\partial T}) \frac{\partial T}{\partial t}] = \\ - ε S_{L} ρ_{L} c_{L} \underline{v} \cdot \underline{\nabla} T + \underline{\nabla} \cdot [(λ \underline{\underline{I}} + ε S_{L} ρ_{L} c_{L} \underline{\underline{D}}) \underline{\nabla} T] + Q_{p} c_{L} (T^{*} - T) + ε (S_{L} ρ_{L} γ_{0}^{L} + S_{I} ρ_{I} γ_{0}^{I}) + (1 - ε) ρ_{S} γ_{0}^{S} \end{array}

.

(35)

1.3. Freeze-Thaw Relations

The water-mass-balance equation (eq. 23), the energy-balance equation (eq. 35), and their associated relations described in the preceding sections require specification of how total-water saturation (S_W) varies with pressure and temperature as in previous versions of SUTRA, but also for the current version require specification of “freeze-thaw” relations that describe how liquid-water saturation (S_L) and ice saturation (S_I) vary with pressure and temperature. SUTRA uses ad-hoc relations to compute S_W, S_L, and S_I as functions of pressure and temperature under saturated and unsaturated conditions. For clarity of presentation, the freeze-thaw relations are first discussed here in the context of fully saturated systems and then in the more general context of unsaturated systems.

Note: All freeze-thaw and unsaturation relations ignore the possibility of hysteresis, for which properties of the porous medium may change depending on the history of freeze-thaw and unsaturation. All relations provide a single value for saturations and porous medium properties depending on temperature and pressure and may be considered to be state functions.

1.3.1. Freeze-Thaw Under Fully Saturated Conditions

When the porous medium is fully saturated, liquid-water and ice saturations, S_L and S_I, are assumed to depend exclusively on ambient temperature (T), not on fluid pressure (p). It is also assumed that pore water cannot freeze at temperatures above the user-specified freezing temperature, T_f (usually set by the user to 0 °C), and that freezing of pore water can occur within a range of temperatures below T_f. In a sense, T_f can be considered as analogous to the user-specified air-entry pressure (p_ent) for unsaturation (see section 2.2, “Preprogrammed Total-Water Saturation Functions for Unsaturated Conditions”), because air does not enter the porous medium (thereby decreasing the liquid saturation) until the pressure drops to the air-entry pressure value, and similarly, ice does not form (thereby decreasing the liquid saturation) until the temperature drops to the freezing temperature. At any temperature below T_f, SUTRA determines the amount of ice and liquid in the porous medium by first calculating the liquid-water saturation (S_L), under fully saturated conditions, according to a user-specified function ( $S_{L}^{sat}$ ) that depends only on temperature:

S_{L} = S_{L}^{sat} (T)

,

(36)

where

S_{L}^{sat}

[1] is liquid-water saturation resulting from the freeze-thaw process in a fully saturated porous medium; a function of temperature (see section 2.3, “Preprogrammed Liquid-Water Saturation Functions for Freezing in Fully Saturated Porous Media”).Then, ice saturation under fully saturated conditions is calculated as

S_{I} = 1 - S_{L} = 1 - S_{L}^{sat} (T)

.

(37)

It is assumed that some residual liquid water always remains during freezing of groundwater, even at very low temperature, although the residual amount may be very small. Therefore, the function $S_{L}^{sat} (T)$ is formulated such that it cannot have a value less than a user-specified minimum value, $S_{Lres}^{sat} > 0$ , where $S_{Lres}^{sat}$ is a dimensionless quantity [1]. The exponential liquid-water saturation function (described in section 2.3, “Preprogrammed Liquid-Water Saturation Functions for Freezing in Fully Saturated Porous Media”) approaches $S_{Lres}^{sat}$ asymptotically. The temperature at and below which the modified power-law function has a value of $S_{Lres}^{sat}$ is determined by the particular shape of the function, which depends on the other input parameters for that function. The temperature at and below which the piecewise-linear function has a value of $S_{Lres}^{sat}$ is specified by the user:

$S_{Lres}^{sat} (x, y [, z])$: [1] is minimum liquid-water saturation that can occur as a result of freezing the porous medium under saturated conditions (user-specified value for each spatial region) and
T_Lres(x, y[, z]): [°C] is for the piecewise-linear liquid-water saturation function, the temperature at and below which the minimum liquid-water saturation under saturated conditions, $S_{Lres}^{sat}$ , occurs (user-specified value for each spatial region).

Note: Often, the values of $S_{Lres}^{sat}$ and T_Lres are not known for the porous medium being simulated. In this case, for the processes being considered, $S_{Lres}^{sat}$ may be set to a relatively low value and T_Lres (for the piecewise-linear function) may be set to an arbitrary but reasonable value. It is advisable to check the sensitivity of simulation results to the values selected, by testing with a range of values in a series of simulations.

As a freeze-thaw simulation for saturated conditions progresses, S_L and S_I, and their derivatives, determined from the user-specified relation $S_{L}^{sat} (T)$ , are used in the various terms of the water-mass-balance equation (eq. 23) and the energy-balance equation (eq. 35) and their associated formulae. Users may either program their own functions or use one of the preprogrammed $S_{L}^{sat} (T)$ relations provided (as discussed in section 2.3, “Preprogrammed Liquid-Water Saturation Functions for Freezing in Fully Saturated Porous Media”).

1.3.2. Freeze-Thaw Under Unsaturated Conditions

Regardless of whether or not freezing occurs at a given point in the model, the total-water saturation (sum of liquid and ice saturations, eq. 1) at that point is assumed to be a user-selected function of pressure alone:

S_W

=

S_W

(

p

).

(38)

This is the same water-saturation function used for unfrozen, unsaturated media, as implemented in previous versions of SUTRA (Voss and Provost, 2002). Preprogrammed S_W(p) relations are discussed in section 2.2, “Preprogrammed Total-Water Saturation Functions for Unsaturated Conditions.” It is assumed that the porous medium cannot desaturate completely, that is, become completely devoid of water, although the remaining amount of water may be very small. Therefore, the relation S_W(p) is formulated such that S_W cannot be less than a user-specified minimum value, S_Wres, with the condition that S_Wres>0, where

S_Wres(x, y[, z]): [1] is residual total-water saturation, the value below which total-water saturation cannot fall because the fluid becomes effectively immobile in the type of porous medium (for example, silt, sand, or gravel) being simulated. This parameter is specified regionwise to be S_Wres>0 (constant value in each spatial region).

For unsaturated conditions (S_W<1), it is assumed (just as it is for fully saturated conditions) that pore water does not freeze at temperatures above the freezing temperature, T_f, and that freezing of pore water can occur at temperatures below T_f. The liquid-water saturation (the liquid that remains after any freezing occurs), S_L, is computed according to equation 36, but is subject to the limitation that S_L cannot exceed the total-water saturation, S_W. Thus,

S_{L} = \min \{S_{L}^{sat} (T), S_{W} (p)\}

,

(39)

where “min” is a function that selects the lesser of its arguments. The ice saturation, S_I, is then given by

S_{I} = S_{W} - S_{L} = \max \{S_{W} (p) - S_{L}^{sat} (T), 0\}

,

(40)

where “max” is a function that selects the greater of its arguments. The SUTRA freeze-thaw relations under unsaturated conditions, equations 39 and 40, simplify to the freeze-thaw relations under saturated conditions, equations 36 and 37, when S_W=1.

The ice saturation in unsaturated conditions is determined by both the liquid water saturation and total saturation state functions (and is thus a function of both temperature and fluid pressure). The ice saturation state function represents the following behavior. For lower total water saturations in the state function due to lower fluid pressure at constant temperature, the liquid water content remains constant and ice content becomes lower. Kurylyk and Watanabe (2013) explain that, in recent laboratory measurements (Jame, 1977; Watanabe and Wake, 2009; Zhou and others, 2014), the total unsaturated water content (which is a function of pressure) does not affect the liquid water content over a range of freezing temperatures. For example, Kurylyk and Watanabe (2013) point out that this physical behavior can be observed in Watanabe and Wake (2009, fig. 2), where for a wide range of frozen soil types, the measured liquid water content depends only on temperature, and does not depend on the total water saturation. Moreover, for even lower total water saturations at which there is no longer ice in the pores, the state functions indicate the liquid saturation becomes equal to the total saturation and decreases with decreasing fluid pressure, as occurs in normal unfrozen porous media. In summary, these state functions imply that the liquid-water saturation in partially frozen unsaturated soils is independent of the total-water saturation, provided that the total water saturation is greater than the maximum possible liquid saturation (the value that occurs for fully saturated conditions as a function of temperature).

In the case where $S_{Lres}^{sat} < S_{Wres}$ (fig. 1A), the lower limit on liquid saturation due to freezing alone has been specified by the user to be lower than the total residual water saturation. This implies that when S_W=S_Wres, which occurs at relatively low pressure, some of the residual total water can be frozen, depending on the temperature, such that the liquid saturation is less than the total-water saturation. The minimum liquid saturation due to freezing alone, $S_{L} = S_{Lres}^{sat}$ , is reached at relatively low temperature. Alternatively, in the case where $S_{Lres}^{sat} > S_{Wres}$ (fig. 1B), the lower limit on liquid saturation due to freezing alone has been specified by the user to be higher than the total residual water saturation. This implies that when S_W=S_Wres, none of the residual total water can freeze, irrespective of how low the temperature may be. In all cases, when both pressure desaturation and freezing are taken into account, the minimum possible liquid saturation, denoted S_Lres, is the lesser of $S_{Lres}^{sat}$ and S_Wres. The value of S_Lres is determined by SUTRA from the user-specified values of $S_{Lres}^{sat}$ and S_Wres.

Ice-free region is in the top left side of the graphs — Figure 1.
Graphs showing the SUTRA freeze-thaw relations for saturated and unsaturated conditions. A, Minimum liquid-water saturation is less than residual water saturation ( $S_{Lres}^{sat} < S_{Wres}$ ), for which ice can exist in the minimum residual total water resulting from desaturation. B, Minimum liquid-water saturation is greater than residual water saturation ( $S_{Lres}^{sat} > S_{Wres}$ ), for which no ice can exist in the minimum residual total water resulting from desaturation. The liquid-water saturation (S_L) and the ice saturation (S_I) depend on the total-water saturation [S_W(p)] and on the liquid-water saturation under fully saturated conditions [ $S_{L}^{sat} (T)$ ]. In effect, the diagram illustrates how S_L and S_I depend on pressure and temperature. On the horizontal axis, S_W is a proxy for pressure (p), and the direction of decreasing p favors desaturation of the medium. On the vertical axis, $S_{L}^{sat}$ is a proxy for temperature (T), and the direction of decreasing T favors freezing. Lines of constant S_L are red, and lines of constant S_I are blue. Light-blue shading indicates the ice-free region of the state diagram. Grey shading indicates the region of the state diagram excluded by the lower bounds on S_W and S_L. See appendix 1 for a list of symbols.

1.3.3. Permeability for Freeze-Thaw Under Saturated and Unsaturated Conditions

Under fully saturated conditions, the relative permeability (k_r) of a partially frozen medium is assumed to depend exclusively on liquid-water saturation (S_L). This is analogous to the approach used for the relative permeability of unfrozen unsaturated media. It is typically assumed that some small residual permeability remains during freezing and unsaturation of groundwater, even for very low liquid-water saturation, because there are still connected pathways or films of liquid water on the mineral (or ice) grains through which flow can occur. Therefore, the relative permeability relation for SUTRA is formulated such that it cannot have a value less than k_rmin(x, y[, z]), a dimensionless quantity [1] that may be specified by the user for each spatial region subject to the condition k_rmin>0. Relative permeability thus ranges from the minimum value of k_rmin to a maximum value of one.

Note: Commonly, the value of k_rmin is not known for the porous medium being simulated. In this case, it may be set to a relatively low value (for example, 10^–6) for the processes being considered to essentially make liquid flow unimportant for conditions with low liquid saturations.

Note: To aid in interpreting simulation results, the value of liquid-water saturation at which k_r=k_rmin, denoted $S_{L_{krmin}}$ , may be determined mathematically by the user from the user-selected relative permeability function. SUTRA does not calculate this value of liquid-water saturation, nor is it required for input. (The exception is for the piecewise-linear relative-permeability function, for which both $S_{L_{krmin}}$ and k_rmin are values that are input by the user.) It is suggested that users calculate the $S_{L_{krmin}}$ value for their chosen function. Then, to improve interpretation of the simulation results, this saturation, $S_{L_{krmin}}$ , may be compared by the user to the value of $S_{Lres}^{sat}$ (the minimum possible liquid saturation for the frozen medium under saturated conditions) or to S_Wres (the minimum possible total-water saturation for unsaturated conditions) to understand whether or not the minimum relative permeability occurs in the model domain, and if so, where and when.

Relative permeability is calculated from an additional relation that must be specified by the user as a function of liquid-water saturation only:

k_r

=

k_r

(

S_L

),

(41)

where k_rmin ≤ k_r ≤ 1 is greater than or equal to k_rmin and less than or equal to one. When S_L=1, k_r=1, and

S_{L} \leq S_{L_{k r min}}

, k_r = k_rmin occurs. The total effective permeability (product of permeability and minimum relative permeability, kk_r) must always be greater than zero as explained in section 2, “Saturation and Relative-Permeability Functions.”

In previous versions of SUTRA, freeze-thaw was not simulated, so total-water saturation, S_W, was used in functions that provide relative permeability values. In this version of SUTRA, freeze-thaw can be simulated, so relative permeability is controlled by liquid-water saturation, S_L, which can be less than total-water saturation. This becomes equivalent to using total-water saturation for simulations in which no ice forms, providing backward compatibility with previous versions of SUTRA. Preprogrammed k_r relations are discussed in section 2, “Saturation and Relative-Permeability Functions.”

1.3.4. Behavior of the Freeze-Thaw Model

The effects of pressure and temperature on the liquid and ice saturations in the SUTRA freeze-thaw model for saturated and unsaturated conditions are illustrated in figure 1. Because the total-water saturation, S_W, is a decreasing function of pressure below atmospheric pressure or below the air-entry pressure (until the lower limit S_Wres is reached), S_W(p) can be considered as a proxy for p on the horizontal axis. Similarly, because the liquid-water saturation under fully saturated conditions, $S_{L}^{sat}$ , is a decreasing function of temperature below the freezing point, T_f (until the lower limit $S_{Lres}^{sat}$ is reached), $S_{L}^{sat} (T)$ can be considered a proxy for T on the vertical axis. The state diagram (fig. 1) can then be understood by considering two cases, as follows:

• Cooling at constant pressure.—A partially saturated, unfrozen medium is cooled while holding p constant. The total-water saturation, S_W, remains constant throughout this process, which can be tracked on the diagram by starting at any point on the top edge ( $S_{L}^{sat} = 1$ ) and proceeding vertically downward to $S_{L}^{sat} = S_{Lres}^{sat}$ . During cooling, the liquid-water saturation, S_L, remains constant and the ice saturation, S_I, remains zero until the temperature is cold enough to begin freezing the pore water. Thereafter, as T decreases below T_f, S_L decreases while S_I increases (that is, liquid becomes ice), until S_L decreases to the residual liquid-water saturation, $S_{Lres}^{sat}$ , after which further decreases in T have no effect on saturations.
• Desaturation at constant temperature.—A saturated, partially frozen medium is desaturated by decreasing p while holding T constant. As p decreases below the air-entry pressure, total-water saturation, S_W, decreases and the ice saturation, S_I, decreases. This can be tracked on the diagram by starting at any point on the right-hand edge of the state diagram (S_W=1) and proceeding horizontally leftward to S_W=S_Wres. Note that the liquid-water saturation, S_L, remains constant as long as some ice remains (S_I>0); thus, desaturation occurs by loss of ice. When pressure decreases to a certain value, all the ice disappears (S_I=0). The liquid-water saturation, S_L, then becomes equal to the total-water saturation, S_W. Thereafter, as p continues to decrease, the liquid and total-water saturations remain equal, and both decrease until their value reaches the residual total-water saturation, S_Wres, after which further decreases in p have no effect on saturations.

For a quantitative description of this behavior and for references, see section 1.3.2, “Freeze-Thaw Under Unsaturated Conditions,” and section 2.3 “Preprogrammed Liquid-Water Saturation Functions for Freezing in Fully Saturated Porous Media.”

1.3.5. General Procedure for Saturation Calculations

The procedure for computing total water, liquid, and ice saturations, their derivatives, and relative permeability has the following logical structure (fig. 2):

1. Compute total-water saturation, S_W:
1. a. If fully saturated conditions (p≥p_ent), set S_W=1.
2. b. Otherwise (p<p_ent), set S_W=S_W(p), where S_W(p) is the user-selected total-water saturation function.
2. Compute liquid-water saturation, S_L:
1. a. If freezing is not possible (T≥T_f), set S_L=S_W.
2. b. Otherwise (T<T_f):
  1. i. Compute $S_{L}^{sat} (T)$ , the user-selected liquid-water saturation function under fully saturated conditions.
  2. ii. Set S_L equal to the lesser of $S_{L}^{sat} (T)$ and S_W(p).
3. Compute ice saturation, S_I=S_W−S_L.
4. Compute derivatives:
1. a. If ice is present (S_I>0):
  1. i. ∂S_L/∂p=0
  2. ii. $\partial S_{L} / \partial T = \partial S_{L}^{sat} / \partial T$
  3. iii. ∂S_I/∂p=dS_W/dp
  4. iv. $\partial S_{I} / \partial T = - \partial S_{L}^{sat} / \partial T$
2. b. Otherwise (S_I=0):
  1. i. ∂S_L/∂p=dS_W/dp
  2. ii. ∂S_L/∂T=0
  3. iii. ∂S_I/∂p=0
  4. iv. ∂S_I/∂T=0
5. Compute relative permeability from the user-specified function k_r(S_L).

During a simulation, the above process occurs once on each time step for each required location (at nodes and Gauss points) or is repeated for each iteration during a time step (until convergence) should iterations for nonlinearity be employed by the user.

How equations in this report are related to one another — Figure 2.
Flow chart showing general procedure for computing total-water, liquid-water, and ice saturations, their derivatives, and relative permeability. See appendix 1 for a list of symbols.

1.4. Thermal Conductivity Functions

In SUTRA version 4.0, the bulk thermal conductivity, λ_LIS, at any point in the model domain may be calculated as a weighted average of the thermal conductivities of liquid water, ice, solid mineral grains, and air. The thermal conductivities of liquid water, ice, and solid mineral grains are averaged by using one of three formulae, and the resulting value is volume-averaged with the thermal conductivity of air, as explained below.

For the combination of liquid, ice, and solid-mineral-grain components, the user-selectable bulk thermal conductivity formulae for three common volumetrically weighted means (Nield and Bejan, 1999) are as follows:

arithmetic mean:

λ_LIS

≡

ε

(

S_Lλ_L

+

S_Iλ_I

) + (1 −

ε

)

λ_S

,

(42A)

geometric mean:

λ_{LIS} \equiv {(λ_{L}^{S_{L}} λ_{I}^{S_{I}})}^{ε} λ_{S}^{(1 - ε)}

, and

(42B)

harmonic mean:

λ_{LIS} \equiv {[ε (\frac{S_{L}}{λ_{L}} + \frac{S_{I}}{λ_{I}}) + \frac{1 - ε}{λ_{S}}]}^{- 1}

,

(42C)

where

λ_L: [E/(s‧L‧°C)] is thermal conductivity of liquid water, user-specified value [about 0.60 J/(s‧m‧°C) at 20 °C; about 0.56 J/(s‧m‧°C) at 0°C] (user-specified constant value for the entire model);
λ_I: [E/(s‧L‧°C)] is thermal conductivity of ice, user-specified value [about 2.1 J/(s‧m‧°C) at 0 °C; about 2.3 J/(s‧m‧°C) at –10 °C] (user-specified constant value for the entire model);
λ_S: [E/(s‧L‧°C)] is thermal conductivity of a solid mineral grain, user-specified value [1 J/(s‧m‧°C)<λ_S<10 J/(s‧m‧°C) for most minerals] (user-specified value for each finite element); and
λ_LIS(x, y[, z], t): is the effective thermal conductivity of mixture of liquid water, ice, and solid mineral grains, calculated by SUTRA according to user-selected equation 42A, B, or C.

In previous versions of SUTRA, only the arithmetic mean (eq. 42A) value was available.

Of the three options available, the arithmetic and harmonic means give the highest and lowest estimates, respectively, of the bulk thermal conductivity. The geometric mean gives an intermediate value that is often a good estimate when experimental data on the bulk thermal conductivity are unavailable (Lunardini, 1981; Nield and Bejan, 1999).

Unsaturated media also have an air phase with influences on effective thermal conductivity defined in the following paragraphs.

λ_A(x, y[, z]): [E/(s‧L‧°C)] is the effective thermal conductivity of air (user-specified value for each finite element) [about 0.025 J/(s‧m‧°C) for static air at 0 °C; for circulating air, the effective value can be orders of magnitude greater than the thermal conductivity of liquid, ice, or solid grains].

The volumetric saturation of the air phase is defined by

S_A

=1−

S_W

=1−

S_I

−

S_L

,

(43)

where

S_A(x, y[, z], t): [1] is air saturation (volume of air per volume of voids).

The effective total thermal conductivity (bulk thermal conductivity), λ, which includes the influence of liquid water, ice, solid mineral grains, and air, is calculated as a volume average (analogous to the arithmetic mean, shown above) of the effective thermal conductivity of the mixture of liquid, ice, and solid grains, λ_LIS, with the air thermal conductivity, λ_A, as follows:

λ

=

λ_LIS

+

εS_Aλ_A

.

(44)

Note: Should the air thermal conductivity be set to zero by the user (in effect, ignored), the value of the total thermal conductivity determined by equation 44 will equal the value calculated for the liquid-ice-mineral mixture, causing the air phase to have no effect.

The air thermal conductivity for static air is practically negligible in comparison with the other thermal conductivities and was set to zero in previous versions of SUTRA. However, should air circulate among pores, thereby transferring energy convectively, its heat transfer effect may be much greater than would be indicated by its low thermal conductivity. Heat transfer by air circulation could be an important process to consider (for example, see Niu and others, 2016), especially in coarse, blocky material that is common in some permafrost environments (Wicky and Hauck, 2017). Because this is a different process than pure conduction, the thermal conductivity of air is not included in the volumetrically weighted λ_LIS values (eq. 42A–C) but is rather included later, as described just above (eq. 44).

Note: Although the physical process of air convection is not explicitly represented in the current version of SUTRA, the effective thermal conductivity of air parameter, λ_A, can be used to approximate the thermal effect of air circulation by setting it to a large value. Values of this parameter may be determined by converting published heat-transfer coefficients for circulating air into effective thermal conductivities. Note that in SUTRA version 4.0, thermal conductivities are constant for the entire simulation and cannot change with time; thus, any required seasonal changes in the effective thermal conductivity of air must be handled by running a series of consecutive simulations in which the parameter value is changed, as desired by the user, from one simulation to the next (that is, from one season to the next).

1.5. Nonlinear Liquid-Water Density Function

A nonlinear function describing the liquid density is implemented in the form of the empirical Thiesen-Scheel-Diesselhorst (TSD) equation (Tilton and Taylor, 1937; Kell, 1967), as described by McKenzie and Voss (2013), with the addition of constant-density limits for very low and very high temperatures. The liquid density, ρ_L [kg/m³], depends on temperature, T [°C], as follows:

ρ_{L} = \{\begin{matrix} 908.0 kg / m^{3} & for T \leq - 52.1784 °C \\ 1, 000.0 [\frac{(T + 288.9414) {(T - 3.9863)}^{2}}{508929.2 (T + 68.12963)}] kg / m^{3} & for - 52.1784 °C < T < 422.9571 °C \\ 500.0 kg / m^{3} & for T \geq 422.9571 °C \end{matrix}

.

(45)

The corresponding derivative of density (in kilograms per cubic meter per degree Celsius) with respect to temperature is

\frac{\partial ρ_{L}}{\partial T} = \{\begin{matrix} 0 kg / (m^{3} \cdot °C) & for T \leq - 52.1784 °C \\ \frac{- 0.00392982 (T - 3.9863) (T^{2} + 246.665 T + 20125.58)}{{(T + 68.12963)}^{2}} kg / (m^{3} \cdot °C) & for - 52.1784 °C < T < 422.9571 °C \\ 0 kg / (m^{3} \cdot °C) & for T \geq 422.9571 °C \end{matrix}

.

(46)

The function, as implemented in SUTRA, is shown in figure 3. In this function, density has a maximum value of 1,000 kg/m³ at 4 °C.

Other than two small horizontal parts, the rest of the line is curved — Figure 3.
Graph showing liquid-water density as a function of temperature (T), as implemented in the SUTRA code. The curved part of the function, from –52.1784 to 422.9571 degrees Celsius is based on the Thiesen-Scheel-Diesselhorst equation (see eq. 45). Low- and high-temperature cut offs provide constant density for temperatures lower and higher, respectively, than the cut-off temperatures. See appendix 1 for a list of symbols. kg/m³, kilogram per cubic meter; °C, degree Celsius.

Note: This nonlinear function is provided as an alternative to the linear density function implemented in all previous versions of SUTRA, which remains available in the current version. Users may either select the linear density function that was available in previous versions of SUTRA (for simulated ranges of temperature in which the nonlinear density values are approximately linear) or this nonlinear function, according to their simulation needs.

For temperatures below –40 °C, density decreases rapidly. This decrease, provided by the TSD function and the constant density below 52.1784 °C, roughly approximates the density given by the two-liquid-phase super-cooled water model of Holten and Anisimov (2012) for typical earth-system fluid hydrostatic pressures between 0.1 megapascal (MPa; atmospheric pressure at ground surface) and 20 MPa (about 2-kilometer depth; Holten and Anisimov, 2012, fig. 3). Although this constant density assumption provides only approximate liquid-density values, errors in density typically will not affect fluid flow or hydraulic behavior, because at temperatures below –40 °C there will be almost no liquid-phase water. Moreover, the relative permeability for any remaining unfrozen water at such low temperatures will be extremely low. For temperatures above 625 °C, density values provided by the TSD nonlinear function become negative. Thus, at an arbitrarily selected maximum temperature, a minimum allowed liquid density value is implemented in SUTRA (fig. 3).

Note: Use of the code for high-temperature simulations must be carried out with caution. Users must be aware that fluid density also depends on fluid pressure. For very high pressure, density may be as much as 10 percent higher than calculated by the nonlinear TSD function. Furthermore, for pressure and temperature ranges encountered in hydrogeologic systems, it is possible that phase change between liquid water and steam may occur, particularly at higher temperatures. SUTRA does not simulate phase change to/from steam; this would require a code that tracks pressure or temperature plus internal energy or enthalpy as primary variables. For high temperature and low pressure, density may be only 20 or 30 percent of the value calculated by the nonlinear density function. Thus, for simulation of such high-temperature geothermal systems, the density value calculated by equation 45 (decreasing with increasing temperature and leveling off at an arbitrary density value) is only a pragmatic proxy that should be used with caution.

For simulation of systems with small ranges of temperature just above and below 4 °C, variable-density effects on fluid flow should be small (unless permeabilities are relatively high) due to near-constant density; however, for temperature ranges that are larger and perhaps more typical of cold regions where ground ice may form (for example, –20 °C mid-winter air temperature to 20 °C summer air temperature, with even greater temperatures at depth in the subsurface), fluid density differences will be sufficient to drive flow wherever the permeability of the geologic fabric is high enough.

2. Saturation and Relative-Permeability Functions

This section describes functions that must be specified for SUTRA simulations that involve unsaturated conditions, freeze-thaw conditions, or both (for saturated simulation without freezing, none of the functions described below are required):

• Unsaturated-flow simulations without freezing.—SUTRA requires user specification of a function to obtain values of total-water saturation and its derivative with respect to pressure, S_W(p), and dS_W/dp. The functionality of total-water saturation functions is discussed in the main SUTRA documentation (Voss and Provost, 2002).
• Saturated-flow simulations with freezing.—SUTRA computes ice saturation, S_I(T), from the liquid-water saturation, $S_{L}^{sat} (T)$ . Thus, for saturated freeze-thaw simulations, SUTRA requires user specification of a saturated-conditions liquid-water saturation function, $S_{L}^{sat} (T)$ , to obtain values of liquid-water saturation and its derivative with respect to temperature, $\partial S_{L}^{sat} / \partial T$ .
• Unsaturated-flow simulations with freezing.—SUTRA determines the ice saturation, S_I(p,T), and liquid saturation, S_L(p,T), from both the liquid-water saturation, $S_{L}^{sat} (T)$ , and the total-water saturation, S_W(p) (see section 1.3.5, “General Procedure for Saturation Calculations;” fig. 2). Thus, for unsaturated freeze-thaw simulations, SUTRA requires user specification of a saturated-conditions liquid-water saturation function, $S_{L}^{sat} (T)$ , a total-water saturation function, S_W(p), and their derivatives with respect to pressure and temperature.

The three required user-selected (or user-defined) functions are listed below:

1. For all unsaturated simulations (with or without freeze-thaw): S_W(p), total-water saturation as a function of pressure (and dS_W/dp) is required.
2. For all freeze-thaw simulations (saturated or unsaturated): $S_{L}^{sat} (T)$ , liquid-water saturation as a function of temperature under fully saturated conditions (and $\partial S_{L}^{sat} / \partial T$ ) is required.
3. For all freeze-thaw simulations and for all unsaturated simulations: k_r(S_L), relative permeability as a function of liquid-water saturation is required (for nonfreezing unsaturated simulations, SUTRA automatically sets liquid saturation, S_L, equal to total-water saturation, S_L=S_W).

Two user-selected (or user-defined) functions are required for unsaturated simulations (functions 1 and 3), three for unsaturated freeze-thaw simulations (functions 1, 2, and 3), and two for saturated freeze-thaw simulations (functions 2 and 3). The user-selectable forms are described in detail in this section. The expressions for the total-water, liquid-water, and ice saturations are equations 38, 39, and 40. Relative permeability is assumed to be a function of liquid-water saturation (eq. 41).

Note: The user selects these functions through the input dataset 11B–D (see app. 2). Particular user-selectable forms of these functions have been made available in SUTRA version 4.0, as described in the following section. Additionally, there is an option for the user to program other forms in Fortran subroutines UNSAT, LIQSAT and RELPERM.

2.1. Overview of Available Functions

Before SUTRA version 4.0, the freezing capability was not available, and functions defining unsaturated properties could be specified only via user-programming of subroutine UNSAT and subsequent recompilation of the source code. The source code came with the van Genuchten (1980) unsaturated functions programmed as an example, with the model parameters set to particular values that could be changed only by reprogramming.

Beginning with SUTRA version 4.0, a number of total water saturation (for unsaturated conditions), liquid-water saturation (for freezing conditions), and relative permeability functions (for unsaturated and freezing conditions) are preprogrammed, and their parameters can be specified through SUTRA input dataset 11B–D (see app. 2). Additionally, the different functions and parameter sets can be set by the user to apply in different spatial regions of the area or volume being modeled. Kurylyk and Watanabe (2013) provide a review of these and other possible functions, with some explanations of their basis and their applicability.

The available preprogrammed total-water saturation functions (for unsaturated flow simulations) are

• van Genuchten (1980),
• Brooks and Corey (1964), and
• piecewise-linear.

The available preprogrammed liquid-water saturation functions (for freeze-thaw simulations) are

• exponential (Mottaghy and Rath, 2006),
• modified power-law (modified for SUTRA from Anderson and Tice, 1972), and
• piecewise-linear.

The available preprogrammed relative-permeability functions are

• van Genuchten (1980),
• Brooks and Corey (1964), and
• piecewise-linear.

The above-listed functions are coded in three Fortran subroutines contained within the file USUBS.f:

• UNSAT (for total-water saturation, S_W(p), and dS_W /dp),
• LIQSAT (for liquid-water saturation of a saturated porous medium, $S_{L}^{sat} (T)$ , and $\partial S_{L}^{sat} / \partial T$ ), and
• RELPERM (for relative permeability depending on liquid-water saturation, k_r(S_L)).

Subroutine ALLSAT evaluates these functions as necessary and calculates the corresponding values of liquid-water saturation, S_L, and ice saturation, S_I, each of which can vary with pressure and temperature.

Additionally, any other functions that the user wishes to employ in simulations can be user-programmed into these subroutines. Any user-programmed functions are recognized within the same framework, which allows users to input their new function parameters through dataset 11B–D. Input instructions for the total-water saturation and liquid-water saturation functions are provided in appendix 2.

All preprogrammed and user-programmed functions operate in the following manner:

• For simulations of unsaturated, nonfreezing conditions, S_W and dS_W/dp are calculated from the user-selected (or user-programmed) total-water saturation function of pressure (and S_L=S_W and S_I=0).
• For simulations of saturated, freezing conditions, $S_{L}^{sat}$ and $\partial S_{L}^{sat} / \partial T$ are calculated from the user-selected (or user-programmed) liquid-water saturation function of temperature (and S_W=1 and $S_{L} = S_{L}^{sat}$ , finally giving $S_{I} = 1 - S_{L}^{sat}$ ).
• For simulations of unsaturated, freezing conditions,
- • S_W and dS_W/dp are calculated from the user-selected (or user-programmed) total-water saturation function of pressure.
- • Then, $S_{L}^{sat}$ and $\partial S_{L}^{sat} / \partial T$ are calculated from the user-selected (or user-programmed) liquid-water saturation function of temperature, assuming fully saturated conditions.
- • Then, S_L is set to the lesser of S_W and $S_{L}^{sat}$ , such that S_L ≤ S_W (finally giving S_I=S_W−S_L).
• For simulations of unsaturated nonfreezing conditions, saturated freezing conditions and unsaturated freezing conditions, relative permeability, k_r, is always based on liquid-water saturation (S_L), as described in section 2.4, “Preprogrammed Relative-Permeability Functions.”

The procedure used in SUTRA to determine saturations and relative permeability is summarized in figure 2.

Note: Preprogrammed and user-programmed total-water saturation, relative-permeability, and liquid-water saturation functions may be used in any combination. For example, the van Genuchten total-water saturation function may be used with the piecewise-linear relative-permeability function and a user-programmed liquid-water saturation function. The function name “NONE” can be user-selected to indicate that no function will be used by SUTRA to define a saturation and relative permeability or both; in this case, SUTRA assumes that the saturation value or relative permeability value or both equal one.

2.2. Preprogrammed Total-Water Saturation Functions for Unsaturated Conditions

When pressure is nonnegative (p≥0), the porous medium is fully saturated (S_W=1). When pressure is negative (p<0), SUTRA computes the total-water saturation as a function of pressure [S_W=S_W(p)]. In some of the functions described below, two of the parameters used are:

p_ent(x, y[, z]): [M/(L‧s²)] which is air-entry pressure, that is, the (negative) pressure above which the porous medium effectively remains fully saturated; normally p_ent<0 (user-specified value for each spatial region), and
S_Wres(x, y[, z]): [1] which is residual total-water saturation, that is, the value below which total-water saturation cannot fall due to decrease in pressure because the fluid becomes effectively immobile (user-specified value for each spatial region).

Note: A value for the air-entry pressure, p_ent, may be calculated from a known hydraulic head value at air entry as p_ent=h_ent/(ρ|g|), where h_ent<0 is the (negative value of) head at air entry in length units, ρ is the fluid density, and |g| is the magnitude of gravitational acceleration.

2.2.1. Van Genuchten Total-Water Saturation Function

In the van Genuchten function, originally formulated for unfrozen groundwater (van Genuchten, 1980), total-water saturation varies with pressure according to

\begin{matrix} S_{W} = (1 - S_{Wres}) [\frac{1}{1 + {(- a_{VG} p)}^{n_{VG}}}]^{(\frac{n_{VG} - 1}{n_{VG}})} + S_{Wres} & for p \leq 0 \end{matrix}

,

(47)

with a natural minimum total water saturation value of S_Wres, where S_Wres and the following two parameters control the shape of the function:

a_VG(x, y[, z]): [(L‧s²)/M] the van Genuchten model parameter usually dependent on inverse of absolute value of the air-entry pressure (user-specified value for each spatial region) and
n_VG(x, y[, z]): [1] the van Genuchten model parameter dependent on pore-size distribution (user-specified value for each spatial region).

An air-entry pressure value cannot be specified for this function. SUTRA uses equation 47 to compute total-water saturation (liquid plus ice) whether or not freezing has occurred. Differentiation of equation 47 with respect to p yields

\frac{d S_{W}}{d p} = \frac{a_{VG} (n_{VG} - 1) (1 - S_{Wres}) {(- a_{VG} p)}^{(n_{VG} - 1)}}{{[1 + {(- a_{VG} p)}^{n_{VG}}]}^{(2 n_{VG} - 1) / n_{VG}}}

.

(48)

2.2.2. Brooks-Corey Total-Water Saturation Function

In the Brooks-Corey function, originally formulated for unfrozen groundwater (Brooks and Corey, 1964), water saturation varies with pressure according to

S_{W} = \{\begin{matrix} 1 & for p \geq p_{ent} \\ (1 - S_{Wres}) {(\frac{p}{p_{ent}})}^{- λ_{BC}} + S_{Wres} & for p < p_{ent}, p_{ent} \neq 0 \end{matrix}

,

(49)

with a natural minimum total water saturation value of S_Wres, where S_Wres, p_ent, and the following parameter control the shape of the function:

λ_BC(x, y[, z]): [1] which is the pore-size distribution index and is user-specified value for each spatial region. In equation 49, p and p_ent are negative.

Differentiation of equation 49 with respect to p yields:

\frac{d S_{W}}{d p} = \{\begin{matrix} 0 & for p \geq p_{ent} \\ - \frac{λ_{BC}}{p_{ent}} (1 - S_{Wres}) {(\frac{p}{p_{ent}})}^{- λ_{BC} - 1} & for p < p_{ent}, p_{ent} \neq 0 \end{matrix}

.

(50)

2.2.3. Piecewise-Linear Total-Water Saturation Function

In the piecewise-linear model, total-water saturation varies with pressure according to

S_{W} = \{\begin{matrix} 1 & for p > p_{ent} \\ (1 - S_{W r e s}) (\frac{p - p_{Wres}}{p_{ent} - p_{Wres}}) + S_{W r e s} & for p_{Wres} \leq p \leq p_{ent} \\ S_{W r e s} & for p < p_{Wres} \end{matrix}

,

(51)

where p_ent, S_Wres, and the p_Wres(x, y[, z]) are the parameters that control the shape of the function; p_Wres(x, y[, z]) [Μ/(L‧s²)] is the pressure at which the sloping linear segment of the function in equation 51 gives S_W=S_Wres (user-specified value for each spatial region); and all pressure values are negative. Differentiation of equation 51 with respect to p yields:

\frac{d S_{W}}{d p} = \{\begin{matrix} 0 & for p > p_{ent} \\ \frac{1 - S_{Wres}}{p_{ent} - p_{Wres}} & for p_{Wres} \leq p \leq p_{ent} \\ 0 & for p < p_{Wres} \end{matrix}

.

(52)

Note: The piecewise-linear function is often used when starting a new unsaturated simulation analysis, because the destabilizing effect of the nonlinearity of this function on the numerical solution is typically less than that of the other nonlinear functions. Once the trial simulation is working well, the user may switch to other, more appropriate functions for the particular porous-medium type being studied. This piecewise-linear saturation function is not intended to accurately represent details of physical behavior in the unsaturated zone.

Note: A useful application of the piecewise-linear function is to approximately simulate the effects of water-table movement (that is, changes in thickness of the saturated zone) on flow and drawdown, without considering behavior within the unsaturated zone in any detail. Specific yield is not explicitly included in the SUTRA equations, and this function may also be used when the only purpose of simulating unsaturated conditions is to track the movement of the water table. This function effectively provides a specific-yield value to the water balance being studied. In this case, the maximum specific yield (which occurs when the porous medium is maximally drained), S_y, is

S_y

=

ε

−

S_Wres

.

(53)

To simulate a certain thickness of unsaturated zone, the user merely needs to specify a pressure and saturation value at the top of the unsaturated zone, where the residual or lowest-desired saturation value will occur. The value of the pressure would be a negative value equivalent to the hydrostatic pressure that would occur there from a pressure of zero at the water table that linearly decreases with upward distance from the water table, and the value of total saturation would be a residual or other user-preferred value at this location. These will define one end of the sloping segment of the piecewise-linear total saturation function, while the other end of this segment is at the water table, with pressure of zero and saturation of one.

2.3. Preprogrammed Liquid-Water Saturation Functions for Freezing in Fully Saturated Porous Media

Many studies have demonstrated that the relation between liquid-water content and subfreezing temperature in fully saturated soils can be expressed using a well-defined relation, often referred to as the soil-freezing curve (for example, Williams, 1967; Anderson and Tice, 1972; Jame, 1977; Kozlowski, 2007). For certain soil types, the soil-freezing curve can be obtained from unsaturated functions for drainage of soils without freezing (soil water characteristic curves) through the Clausius-Clapeyron equation, using sorptive or capillary theory (Koopmans and Miller, 1966; Spaans and Baker, 1996; Kurylyk and Watanabe, 2013). In SUTRA version 4.0, soil-freezing curves are specified independently of the unsaturated functions to allow maximum flexibility in functionality.

The user-specified freezing temperature below which pore water begins to freeze is T_f. Surface energies in pores and the existence of solutes may lower the freezing temperature to below 0 °C, so the user-specified value of T_f may be below 0 °C (as discussed in sections 1.2, “Energy Balance With Freeze-Thaw,” and 1.3, “Freeze-Thaw Relations”). When the temperature is above T_f (that is, T>T_f), the pore water is entirely in liquid form (S_L=1). When the temperature is at or below the freezing temperature (T ≤ T_f), ice can begin to form, and SUTRA computes the liquid-water saturation as a function of temperature [S_L=S_L(T)].

The functions described below provide liquid-water saturation under fully saturated conditions, $S_{L}^{sat}$ , as a function of temperature, T, for temperatures below the freezing temperature, T_f. For unsaturated conditions, SUTRA uses $S_{L}^{sat}$ to compute the actual liquid-water saturation, S_L, as described in section 1.3.5, “General Procedure for Saturation Calculations.” Then, ice saturation is computed as depending on total-water saturation (determined as described in section 2.2, “Preprogrammed Total-Water Saturation Functions for Unsaturated Conditions”) and on the liquid saturation as S_I=S_W−S_L. Residual liquid-water saturation, $S_{Lres}^{sat}$ , is the user-specified value below which it is assumed the liquid-water saturation cannot fall due to freezing alone.

The preprogrammed functions provide the liquid-water saturation under fully saturated conditions, $S_{L}^{sat}$ . Three preprogrammed functions are currently provided in SUTRA—exponential, modified power law, and, piecewise linear (see fig. 4). Each of the three functions is expressed in terms of a relative temperature, T_∆, defined as the difference between the local porous medium temperature, T, and the user-specified maximum freezing temperature of pore water, T_f:

T_∆

≡

T

−

T_f

.

(54)

When T_f=0 °C (the most common user choice), the relative temperature becomes the actual temperature, T_∆=T.

Note: The shape of each of the three preprogrammed liquid-water saturation ( $S_{L}^{sat}$ ) functions depends on one or more parameters, the values of which are set by the user. Setting T_f<0 °C simply shifts the $S_{L}^{sat}$ curve in the negative direction along the temperature axis, without changing its shape. Thus, the function should be designed by the user for the un-shifted case and these parameter values should therefore always be set by the user for the unshifted case, T_f=0 °C.

The shapes of these functions for parameter values that were selected to highlight the differences in function shapes are shown in figure 4. The power law function tends to undergo a steep initial decrease in liquid water saturation as temperatures drop below the freezing temperature, and then a gradual decrease to the minimum liquid saturation. The exponential function tends to exhibit a delay in onset of freezing at temperatures just below the freezing temperature, then a steep drop in liquid saturation, quickly reaching and remaining at the minimum liquid saturation for lower temperatures. The linear function can be configured to drop slowly or quickly when freezing begins and then it reaches and remains at the minimum liquid saturation.

The lines in the first graph are almost entirely flat and shoot up at the end, whereas
in the second graph, climb early and even out — Figure 4.
Graphs showing examples of preprogrammed liquid-water saturation function ( $S_{L}^{sat}$ ) shapes used in SUTRA simulations for freezing in fully saturated porous media. A, Plotted as function of T. B, Plotted as function of log₁₀(T). See appendix 1 for a list of symbols. °C, degrees Celsius. Example parameter values: all functions ( $S_{Lres}^{sat}$ = 0.05, T_f=0.0 °C), modified power law (α_POW=0.1 ${(° C)}^{- β_{POW}}$ , β_POW=−0.5, for which T_POW=−0.01 °C), exponential (w_EXP=0.6 °C), and piecewise linear (T_Lres=−2.0 °C).

2.3.1. Exponential Liquid-Water Saturation Function

In the exponential function, adapted from Lunardini (1988) and Mottaghy and Rath (2006), liquid-water saturation during freezing, under fully saturated flow conditions, varies with temperature according to

S_{L}^{sat} = \{\begin{matrix} 1 & for T_{Δ} \geq 0 \\ (1 - S_{Lres}^{sat}) \exp \{- {(\frac{T_{Δ}}{w_{EXP}})}^{2}\} + S_{Lres}^{sat} & for T_{Δ} < 0 \end{matrix}

,

(55)

where exp is the exponential function and T_Δ is the relative temperature defined by equation 54. The shape of the function is controlled by

S_{Lres}^{sat}

, the asymptotic lower limit of the function, and by

w_EXP(x, y[, z]): [°C] an exponential-model parameter (specified by the user for each spatial region, subject to w_EXP>0).

Differentiation of equation 55 with respect to T yields:

\frac{d S_{L}^{sat}}{d T} = - (1 - S_{Lres}^{sat}) (\frac{2}{w_{EXP}}) (\frac{T_{Δ}}{w_{EXP}}) \exp \{- {(\frac{T_{Δ}}{w_{EXP}})}^{2}\}

.

(56)

2.3.2. Modified Power-Law Liquid-Water Saturation Function

One often-used soil freezing curve relates the liquid-water saturation to the temperature in accordance with a power law in the form, $S_{L}^{sat} = α {(- T)}^{β}$ (Anderson and Tice, 1972; Anderson and Morgenstern, 1973), where β<0. Although this function does fit a certain range of laboratory data on freezing geologic materials, it is not appropriate for general simulation over the full range of possible temperatures because it causes $S_{L}^{sat}$ to become infinite when T=0. The following modified form of the power law function (with the desired maximum value of 1 when T=T_f) is used in SUTRA:

S_{L}^{sat} = \{\begin{matrix} 1 & for T_{Δ} \geq T_{POW} \\ α_{POW} {(- T_{Δ})}^{β_{POW}} & for T_{Δ} < T_{POW} \end{matrix}

,

(57)

subject to the restriction that

S_{L}^{sat} \geq S_{Lres}^{sat}

, where T_Δ is the relative temperature defined by equation 54, and

S_{Lres}^{sat}

and the following parameters control the shape of the function:

β_POW(x, y[, z]): [1] is the modified power-law model parameter (β_POW<0; user-specified value for each spatial region) and
α_POW(x, y[, z]): [ ${(° C)}^{- β_{POW}}$ ] is the modified power-law model parameter (α_POW>0) (user-specified value for each spatial region).

T_{POW} \equiv - α_{POW}^{- (1 / β_{POW})}

is the relative temperature at and above which the power-law function has a value of 1. T_POW is not a user-defined value. It is calculated by SUTRA depending on values of α_POW and β_POW, which have user-defined values.

Differentiation of equation 56 with respect to T yields:

\frac{\partial S_{L}^{sat}}{\partial T} = \{\begin{matrix} 0 & for T_{Δ} \geq T_{POW} and S_{L}^{sat} > S_{Lres}^{sat} \\ - α_{POW} β_{POW} {(- T_{Δ})}^{β_{POW} - 1} & for T_{Δ} < T_{POW} and S_{L}^{sat} > S_{Lres}^{sat} \\ 0 & for S_{L}^{sat} = S_{Lres}^{sat} \end{matrix}

.

(58)

Many literature sources express the power-law soil-freezing curve in terms of percent moisture content by dry weight rather than, as done for SUTRA, by volumetric liquid content (liquid-water saturation):

ϕ_{L}^{sat} = α_{ϕ} {(- T)}^{β_{POW}}

,

(59)

where

ϕ_{L}^{sat}

is percent moisture content by dry weight (mass of liquid per mass of dry solid) under fully saturated conditions; α_ϕ is power-law model parameter (α_ϕ>0), as a percentage of

{(° C)}^{- β_{POW}}

, and it is assumed that T_f=0 °C (that is, that T_∆=T). In that case, the value of α_POW in equation 57 can be computed from α_ϕ by

α_{POW} = (\frac{ρ_{S}}{ρ_{L}}) (\frac{1 - ε}{ε}) \frac{α_{ϕ}}{100}

.

(60)

Andersland and Ladanyi (1994, p. 40) present values of α_ϕ and β_POW for 33 different types of soils, and Anderson and Tice (1972) demonstrate that these parameters can be obtained empirically from the specific surface area of the soil grains.

2.3.3. Piecewise-Linear Liquid-Water Saturation Function

In the piecewise-linear function, liquid-water saturation during freezing, under fully saturated flow conditions, varies with temperature according to

S_{L}^{sat} = \{\begin{matrix} 1 & for T_{Δ} \geq 0 \\ (1 - S_{Lres}^{sat}) (1 - \frac{T_{Δ}}{T_{Lres}}) + S_{Lres}^{sat} & for T_{Lres} < T_{Δ} < 0 \\ S_{Lres}^{sat} & for T_{Δ} \leq T_{Lres} \end{matrix}

,

(61)

where T_∆ is the relative temperature defined by equation 54, and

S_{Lres}^{sat}

and the relative temperature T_Lres control the shape of the function. Differentiation of equation 61 with respect to T yields

\frac{d S_{L}^{sat}}{d T} = \{\begin{matrix} 0 & for T_{Δ} \geq 0 \\ \frac{1 - S_{Lres}^{sat}}{- T_{Lres}} & for T_{L r e s} < T_{Δ} < 0 \\ 0 & for T_{Δ} \leq T_{Lres} \end{matrix}

.

(62)

2.4. Preprogrammed Relative-Permeability Functions

The relative permeability for unsaturated and freezing groundwater systems in SUTRA version 4.0 is assumed to depend only on liquid saturation. Relative-permeability functions are thus expressed below in terms of liquid-water saturation, S_L (which automatically becomes equal to the total saturation, S_W, should there be no ice at the point being considered).

It is also necessary that the total permeability value (kk_r) always be greater than zero; if not, then numerical solution becomes impossible. This requires that k_r>0. The user may specify a minimum possible value of relative permeability, k_rmin, to avoid the occurrence of k_r=0. Then, during simulation, if the user-selected relative permeability function would provide a value lower than k_rmin, the value is set by SUTRA to k_rmin.

Note: For relative permeability functions other than the piecewise linear function, if it is desired that the relative permeability be allowed to asymptotically approach zero, and not be limited by k_rmin, the user may set k_rmin=0. For the piecewise linear relative permeability function, k_rmin>0 always determines the lower limit on relative permeability, which is reached at a user-specified liquid-water saturation, $S_{L_{krmin}}$ .

Note: If no relative permeability function is specified by the user, then relative permeability will have a constant value of one, k_r=1.

For convenience and to simplify the form of most relative-permeability expressions presented below, a drainable liquid-water saturation, $S_{L}^{*}$ , (equivalent to the ratio of the volume fraction in the pores of drainable liquid, S_L−S_Lres, to the maximum possible volume fraction of drainable liquid in the pores, 1−S_Lres) is defined as

S_{L}^{*} \equiv \frac{S_{L} - S_{Lres}}{1 - S_{Lres}}

,

(63)

where S_Lres is the minimum possible (residual) liquid-water saturation, taking into account both groundwater desaturation and freezing.

Note: The value of S_Lres is not specified directly by the user. It is determined automatically by SUTRA in accordance with equation 39. When groundwater desaturation and freezing occur simultaneously, S_Lres is the lesser of the user-specified residual saturations $S_{Lres}^{sat}$ and S_Wres. In the absence of freezing, S_Lres=S_Wres. Under fully saturated conditions with freezing, $S_{Lres} = S_{Lres}^{sat}$ .

2.4.1. Van Genuchten Relative-Permeability Function

For the van Genuchten function, which was originally formulated for unfrozen porous media (van Genuchten, 1980), relative permeability varies with liquid-water saturation, S_L, according to

k_{r} = S_{L}^{*}^{1/2} {\{1 - {[1 - S_{L}^{*}^{(\frac{n_{VG}}{n_{VG} - 1})}]}^{(\frac{n_{VG} - 1}{n_{VG}})}\}}^{2} k_{r} \geq k_{rmin}

,

(64)

where

S_{L}^{*}

is the drainable liquid-water saturation defined by equation 63, and n_VG is a model parameter that appears in equation 47. Thus, k_rmin and two parameters, n_VG and S_Lres (which is the lesser of S_Wres and

S_{Lres}^{sat}

) control the shape of this function.

2.4.2. Brooks-Corey Relative-Permeability Function

For the relation presented by Brooks and Corey (1964), originally formulated for unfrozen porous media, relative permeability varies with pressure according to

k_{r} = \{\begin{matrix} 1 & for p > p_{ent} \\ {(\frac{p}{p_{ent}})}^{- 2 - 3 λ_{BC}} & for p \leq p_{ent}, p_{ent} \neq 0 \end{matrix} k_{r} \geq k_{rmin}

,

(65)

where pressures are negative, and p_ent and λ_BC are the air-entry pressure and pore size distribution index, respectively, which also appear in equation 49. Solving equation 49 for p, substituting into equation 65, and replacing S_W with S_L gives

k_{r} = {(S_{L}^{*})}^{3 + 2 / λ_{BC}} k_{r} \geq k_{rmin}

,

(66)

where

S_{L}^{*}

is the drainable liquid-water saturation defined by equation 63. The shape of this function is controlled by k_rmin and two parameters, λ_BC and S_Lres (which is the lesser of S_Wres and

S_{Lres}^{sat}

).

2.4.3. Piecewise-Linear Relative-Permeability Function

In the piecewise-linear model, relative permeability varies with liquid-water saturation, S_L, having a value of 1 when S_L=1 $(S_{L}^{†} = 1)$ and decreasing to a minimum value of k_rmin when $S_{L} = S_{L_{krmin}} (S_{L}^{†} = 0)$ , according to

k_{r} = (1 - k_{rmin}) S_{L}^{†} + k_{rmin}

,

(67A)

where

S_{L}^{†} \equiv \frac{S_{L} - S_{L_{krmin}}}{1 - S_{L_{krmin}}}

.

(67B)

Thus, user-specified parameters

S_{L_{krmin}}

and k_rmin control the shape of this function.

2.4.4. Note Regarding the “Impedance Factor” Approach

Some early studies (for example, Taylor and Luthin, 1978; Jame and Norum, 1980; Lundin, 1990) used an empirical “impedance factor” to incorporate the effect of freezing on the effective permeability (or hydraulic conductivity) of an unsaturated porous medium. The utility of the impedance factor was later questioned, however, and it is now commonly viewed as an arbitrary correction that is not physically based (Newman and Wilson, 1997; Kurylyk and Watanabe, 2013; Zhao and others, 2013). Recent studies have shown that properly parameterized classic unsaturated permeability functions can yield simulation results that concur with laboratory data of water flow in partially frozen soils without the inclusion of an impedance factor (for example, Painter, 2011; Azmatch and others, 2012). SUTRA version 4.0 therefore does not support use of an impedance function, but the following information is provided to allow users to manually program an impedance function into SUTRA, if its use is desired.

When used to describe permeability of frozen porous media, the impedance function controls effective permeability as follows:

k_{eff} = Ω k_{eff}^{uf}

,

(68)

where

k_eff: [L²] is effective permeability that takes into account both partial saturation and freezing,
$k_{eff}^{uf}$: [L²] is effective permeability that takes into account partial saturation under unfrozen conditions, and
Ω: [1] is the impedance factor (usually in the form 10^−constant).

In the SUTRA freeze-thaw model,

k_eff

=

kk_r

,

(69)

where

k: [L²] is permeability under fully saturated, unfrozen conditions (entered in dataset 15 of the “.inp” file) and
k_r: is the relative permeability discussed in previous sections.

Evaluation of equation 69 under unfrozen conditions and substitution into equation 68 gives

k_{eff} = Ω k_{r}^{uf} k

,

(70)

where

k_{r}^{uf}

[1] is the relative permeability under unfrozen conditions.

Comparison of equations 69 and 70 indicates that the relative permeability is related to the impedance factor by

k_{r} = Ω k_{r}^{uf}

.

(71)

The relative permeability under unfrozen conditions is a function of the total-water saturation,

k_{r}^{uf} = k_{r}^{uf} (S_{W})

, and the impedance factor is typically expressed as a function of the ice content and porosity, Ω=Ω(S_I, ε). A review of impedance functions, listing three separate types, is presented by Kurylyk and Watanabe (2013). Therefore, use of the typical impedance-factor approach requires the relative permeability to be a function of the total-water and liquid-water saturations and porosity, k_r=k_r(S_W, S_L, ε).

The preceding discussion was framed in terms of scalar permeability for the sake of simplicity. Analogous arguments can be made in the case of the tensor permeability allowed in SUTRA. Also note that relative permeability remains a scalar quantity.

Note: In the preprogrammed functions available in this version of SUTRA, k_r is a function only of S_L. However, subroutine RELPERM provides the infrastructure necessary for the user to program a more general relative-permeability function that is consistent with the impedance-factor approach. If the user wishes to make the impedance factor dependent on ice saturation, S_I, this can be simply accomplished in subroutine RELPERM by using the relation S_I=S_W−S_L (eq. 40).

3. Spatial Variation of Formerly Constant Parameters

In earlier versions of SUTRA, solid-matrix properties, adsorption parameters, and parameters for production of solute mass or energy were each assigned a single value that applied throughout the model domain. Beginning with this version of SUTRA (version 4.0), most material properties vary spatially. The region-wise values are input in datasets 11A through 11E (see app. 2).

All SUTRA parameters that describe material properties, indicating their allowed spatial variability, are listed in table 1. The spatial variability (that is, the discretization) provided for each parameter depends on how it appears in the numerical formulation of the flow and transport equations in SUTRA. Cellwise parameters are associated only with time-derivative terms and with terms that do not involve spatial derivatives. Nodewise parameters are those that vary spatially across finite elements, bilinearly in two dimensions (2D) and trilinearly in three dimensions (3D), and these are used for representing terms that involve spatial fluxes and spatial derivatives. Elementwise parameters are constant throughout each 2D or 3D finite element, which simplifies the numerical integration over each element of the terms in which they appear. (See Voss and Provost, 2002, for detailed information on numerical methods.) Some parameters have two types of spatial variability because they appear in two different types of terms. Detailed input instructions are listed in appendix 2.

Table 1.

Spatial variation of Saturated-Unsaturated Transport (SUTRA) code parameters.

^{[The user specifies input parameters by providing a single, constant value (CONST),
or by assigning a value to each node (NODE), element (ELEM), or region (REG). Within
SUTRA, spatial variation of parameters can be constant (CONST), for which a single
value of the parameter applies throughout the model; cellwise (CELL), for which the
parameter is constant over each cell surrounding each node; nodewise (NODE), for which
the parameter varies bilinearly for two dimensions (2D) and trilinearly for three
dimensions (3D) among nodes within each finite element; and elementwise (ELEM), for
which the parameter is constant within each finite element. Parameter superscripts/subscripts
are W for total water, L for liquid, I for ice, S for solid matrix, solid mineral
grain, or solute/sorbate. NA, not applicable]}

Table 1. Spatial variation of Saturated-Unsaturated Transport (SUTRA) code parameters.
Parameter name	Symbol or abbreviation	Fortran name	Unit	User input	Spatial variation
Liquid compressibility	β_L	COMPFL	[M/(L‧s²)]⁻¹	CONST	CONST
Liquid specific heat	c_L	CL	[E/(M‧°C)]	CONST	CONST
Solute diffusivity in liquid or liquid thermal conductivity	σ_L (D_m or λ_L)	SIGMAL	[L²/s] or [E/(L‧°C‧s)]	CONST	CONST
[for linear density variation] Density of liquid at base concentration or temperature	ρ₀	RHOLØ	[M/L³]	CONST	CONST
[for linear density variation] Base concentration or temperature	U₀	URHOLØ	[M_s/M] or [°C]	CONST	CONST
[for linear density variation] Coefficient of liquid density change with concentration or temperature	∂ρ_L/∂U	DRLDU	[M²/(L³‧M_s)] or [M/(L³‧°C)]	CONST	CONST
For solute transport: fluid viscosity	µ_L	VISCØ	[M/(L‧s)] or [units conversion factor]	CONST	for solute transport: CONST
For energy transport:					for energy transport:
[for temperature-dependent viscosity] Scale factor (positive value) that multiplies the SUTRA-calculated, T-dependent liquid viscosity (in [kg/m‧s])	µ_L	VISCØ	[M/(L‧s)] or [units conversion factor]	CONST	CONST
[for constant viscosity] Specified as a negative number. The absolute (positive) value of this number is used as the constant (T-independent) liquid viscosity (in user-preferred units that are consistent with the units used throughout the model input, for example, [kg/m‧s]). This constant value of viscosity is used throughout the entire model (instead of a SUTRA-calculated, T-dependent viscosity).	µ_L	VISCØ	[M/(L‧s)]	CONST	CONST
Effective thermal conductivity of air (for energy transport only)	λ_A	SIGMAA	[E/(s‧L‧°C)]	ELEM	ELEM
Ice compressibility	β_I	COMPI	[M/(L‧s²)]⁻¹	CONST	CONST
Ice specific heat	c_I	CI	[E/(M‧°C)]	CONST	CONST
Ice density	ρ_I	RHOI	[M/L³]	CONST	CONST
Ice thermal conductivity	λ_I	SIGMAI	[E/(s‧L‧°C)]	CONST	CONST
Latent heat of fusion	ΔH	HTLAT	[E]	REG	CELL
Freezing temperature of liquid water	T_f	TFREEZ	[°C]	REG	CELL, ELEM
Solid matrix compressibility	α_S	COMPMA	[M/(L‧s²)]⁻¹	NODE	CELL
For energy transport: solid grain specific heat	c_S	CS1	[E/(M‧°C)]	NODE	CELL
For solute transport: arbitrary value (SUTRA sets this to zero)	c_S	CS1	NA	NODE	CELL
For energy transport: thermal conductivity of a solid grain	σ_S	SIGMAS	[E/(L‧°C‧s)]	ELEM	ELEM
For solute transport: arbitrary value (SUTRA sets this to zero)	σ_S	SIGMAS	NA	ELEM	ELEM
For energy transport: value always required	ρ_S	RHOS	[M/L³]	NODE	CELL
For solute transport: density of a solid grain, value required only for sorption simulations	ρ_S	RHOS	[M/L³]	NODE	CELL
For solute transport only, name of adsorption model	Input values: NONE LINEAR FREUNDLICH LANGMUIR	ADSMOD	Unitless names for: No model Linear Freundlich Langmuir	REG	CELL
For solute transport only, linear sorption coefficient, or first Freundlich or Langmuir distribution coefficient	χ₁	CHI1	[L³_f /M_G]	REG	CELL
For solute transport only, second Freundlich or Langmuir distribution coefficient; arbitrary value for linear sorption	χ₂	CHI2	Freundlich [dimensionless quantity]	REG	CELL
For solute transport only, second Freundlich or Langmuir distribution coefficient; arbitrary value for linear sorption	χ₂	CHI2	Langmuir [L³_f /M_S]	REG	CELL
Name of total-water saturation model for unsaturated conditions	Input values: NONE VGEN BCOR PLIN UDEF	SWMOD	Unitless names for: No model Van Genuchten Brooks-Corey Piecewise-linear User-defined	REG	CELL, ELEM
Residual (irreducible) total-water saturation	S_Wres	SWRES	Dimensionless quantity	REG	CELL, ELEM
Air-entry pressure (used for some of the total saturation functions)	p_ent	PENT	[Μ/(L‧s²)]	REG	CELL, ELEM
Parameters of total-water saturation model for unsaturated conditions	VGEN: S_Wres, α_VG, n_VG	VGEN: SWRES, AA, VN	VGEN: dimensionless quantity, [(L‧s²)/M], dimensionless quantity	REG	CELL, ELEM
Parameters of total-water saturation model for unsaturated conditions	BCOR: S_Wres, p_ent, λ_BC	BCOR: SWRES, PENT, RLAMB	BCOR: dimensionless quantity, [M/(L‧s²)], dimensionless quantity	REG	CELL, ELEM
Parameters of total-water saturation model for unsaturated conditions	PLIN: S_Wres, p_ent, p_Wres	PLIN: SWRES, PENT, PSWRES	PLIN: dimensionless quantity, [M/(L‧s²)], [M/(L‧s²)]	REG	CELL, ELEM
Parameters of total-water saturation model for unsaturated conditions	UDEF: user-defined parameters	UDEF: SWPAR	UDEF: user-defined	REG	CELL, ELEM
Name of liquid-water saturation model for freeze-thaw conditions	Input values: NONE EXPO POWR PLIN UDEF	SLMOD	Unitless names for: No model Exponential Modified power-law Piecewise-linear User-defined	REG	CELL, ELEM
Parameters of liquid-water saturation model for freeze-thaw conditions	EXPO: $S_{Lres}^{sat}$ , w_EXP	EXPO: SLSATRES, W	EXPO: dimensionless quantity, [°C]	REG	CELL, ELEM
Parameters of liquid-water saturation model for freeze-thaw conditions	POWR: $S_{Lres}^{sat}$ , α_POW, β_POW	POWR: SLSATRES, ALPHA, BETA	POWR: dimensionless quantity, [(°C)⁻^β^POW], dimensionless quantity	REG	CELL, ELEM
Parameters of liquid-water saturation model for freeze-thaw conditions	PLIN: $S_{Lres}^{sat}$ , T_Lres	PLIN: SLSATRES, TLRES	PLIN: dimensionless quantity, [°C]	REG	CELL, ELEM
Parameters of liquid-water saturation model for freeze-thaw conditions	UDEF: user-defined parameters	UDEF: SLPAR	UDEF: user-defined	REG	CELL, ELEM
Name of relative-permeability model for unsaturated and freezing conditions	Input values: NONE VGEN BCOR PLIN UDEF	RKMOD	Unitless names for: No model Van Genuchten Brooks-Corey Piecewise-linear User-defined	REG	CELL, ELEM
Relative-permeability model parameters	VGEN: n_VG, k_rmin	VGEN: VN, RKMIN	VGEN: dimensionless quantity, dimensionless quantity	REG	ELEM
Relative-permeability model parameters	BCOR: λ_BC, k_rmin	BCOR: RLAMB, RKMIN	BCOR: dimensionless quantity, dimensionless quantity	REG	ELEM
Relative-permeability model parameters	PLIN: $S_{L_{krmin}}$ , k_rmin	PLIN: SLRKMIN, RKMIN	PLIN: dimensionless quantity, dimensionless quantity	REG	ELEM
Relative-permeability model parameters	UDEF: user-defined parameters	UDEF: RKPAR	UDEF: user-defined	REG	ELEM
Zero-order rate of solute mass or energy production in the liquid	$γ_{0}^{L}$	PRODL0	[(M_s/M)/s] for solute mass production; [(E/M)/s] for energy production	NODE	CELL
Zero-order rate of solute mass production in the ice (arbitrary value for solute transport and for energy transport without freezing, for which SUTRA sets this to zero)	$γ_{0}^{I}$	PRODI0	[s^-1]	NODE	CELL
Zero-order rate of solute mass or energy production in the solid grains	$γ_{0}^{S}$	PRODS0	[(Ms/M_G)/s] for adsorbate mass production; [(E/M_G)/s] for energy production	NODE	CELL
First-order rate of solute mass production in the liquid (arbitrary value for energy; SUTRA sets this to zero)	$γ_{1}^{L}$	PRODL1	[s^-1]	NODE	CELL
First-order rate of solute mass production in the solid grains (arbitrary value for energy; SUTRA sets this to zero)	$γ_{1}^{S}$	PRODS1	[s^-1]	NODE	CELL
Porous medium porosity	ε	POR	Dimensionless quantity	NODE	CELL
Maximum permeability	k_max	PMAX	[L²]	ELEM	ELEM
Middle permeability (3D only)	k_mid	PMID	[L²]	ELEM	ELEM
Minimum permeability	k_min	PMIN	[L²]	ELEM	ELEM
First permeability-direction angle	θ₁	ANGLE1	[°]	ELEM	ELEM
Second permeability-direction angle (3D only)	θ₂	ANGLE2	[°]	ELEM	ELEM
Third permeability-direction angle	θ₃	ANGLE3	[°]	ELEM	ELEM
Longitudinal dispersivity for flow in max-perm direction	α_Lmax	ALMAX	[L]	ELEM	ELEM
Longitudinal dispersivity for flow in mid-perm direction (3D only)	α_Lmid	ALMID	[L]	ELEM	ELEM
Longitudinal dispersivity for flow in min-perm direction	α_Lmin	ALMIN	[L]	ELEM	ELEM
Transverse dispersivity for flow in max-perm direction	α_Tmax	ATMAX	[L]	ELEM	ELEM
Transverse dispersivity for flow in mid-perm direction (3D only)	α_Tmid	ATMID	[L]	ELEM	ELEM
Transverse dispersivity for flow in min-perm direction	α_Tmin	ATMIN	[L]	ELEM	ELEM
Gravitational acceleration components	g_x, g_y, (2D and 3D) g_z (3D only)	GRAVX, GRAVY GRAVZ	[L/s²]	CONST	CONST

4. Enhanced Output

In SUTRA version 4.0, two output options have been enhanced. These are the budget output, which has additional values that explain the model calculations, and the velocity output, to which “Darcy velocity” has been added as a further option.

4.1. Enhanced Budget Output

Budget output, an optional user-selected SUTRA capability, is a postprocessing step that reports additional simulation results on user-selected time steps. The budget output that was available in previous versions of SUTRA is enhanced in the current version to provide additional results that may be helpful in interpreting simulated behaviors, especially for understanding complex interactions that may occur in freeze-thaw simulations. The additional results include rates of change of water mass, and either solute mass or total energy, for subcomponent parts of the values previously reported and for new terms associated with ice and freeze-thaw. Enhanced budget output is here described for all possible types of SUTRA simulations, including fluid mass, energy, and solute balances.

Most budget output values reported are representative of the model domain as a whole. These consist of the total rates of gain or loss of water mass and solute mass or energy, and in the enhanced output, component rates for each process accounted for by the fluid mass balance and solute mass or energy-balance equations being solved by SUTRA. In addition, localized information within the model domain is reported in the budget output, consisting of all water mass and solute mass or energy sources that occur at user-specified boundary condition nodes, as was done in previous SUTRA versions. Budget output is written to an output file that is typically given the filename extension “.lst.” Users may select the frequency of these outputs to occur after particular or all time steps.

Detailed information about the quantitative meaning of each value listed in the enhanced budget output is provided in this section by associating terms in the SUTRA balance equations with the values that appear in the budget output. To link the equations and printed budget output, the Fortran source code variable name is shown in the sections below next to the term (or terms) in the governing equation that the value represents and next to the line in the budget output where that value is printed.

4.1.1. Water-Mass Budget for Energy-Transport Simulations

The components of the water water-mass budget output depend on whether SUTRA is simulating energy transport or solute transport. In this section, the enhanced water-mass budget for energy transport simulations is presented. The manner in which the variables are explained here is exemplified by the Fortran source code variable STW, which is shown here in the SUTRA water-mass balance (reproduced below from equation 4):

\overset{STW ​ total water-mass change rate f (p, T)}{\overset{︷}{\frac{\partial (ε S_{L} ρ_{L})}{\partial t} + \frac{\partial (ε S_{I} ρ_{I})}{\partial t}}} = - \underline{\nabla} \cdot (ε S_{L} ρ_{L} \underline{v}) + Q_{p}

,

(72)

where f(p,T) indicates dependence on both pressure changes and temperature changes.

This value represents the total water-mass change rate, given by the first two terms of equation 72 [STW total water-mass change rate f(p,T)]. The text that accompanies the label “STW” in equation 72 provides the meaning of this value and indicates that STW depends on changes in both p and T. Each term of the enhanced water-mass budget is defined in an equation below and, with each term, an indication is given of whether each rate depends on pressure changes, temperature changes, or both.

The value of STW is printed in the water-mass-budget output at the end of each user-selected completed time step. STW may be found in the following line of the budget output (shown in its entirety after equation 74 below):

TOTAL RATE OF CHANGE IN H2O STORED IN MODEL [ S+, S-, S] —STW—

The variable name, STW, is highlighted as —STW— to indicate that this text is not actually written in the budget output and is included here for informational purposes only.

The other terms in the enhanced water-mass-balance output for energy transport are obtained by combining an expanded form of equation 4, given by equation 19, with the definition of effective storativity (eq. 20). Substitution of equation 20 into equation 19, and rearrangement of the water-mass balance two different ways to gather terms with similar physical origins, provides the definitions of the other water-mass budget for energy transport values:

\begin{array}{l} \overset{STW t o t a l w a t e r - m a s s c h a n g e r a t e f (p, T)}{\overset{︷}{(S_{L} ρ_{L} S_{op}^{L} + S_{I} ρ_{I} S_{op}^{I}) \frac{\partial p}{\partial t} + ε (ρ_{L} \frac{\partial S_{L}}{\partial p} + ρ_{I} \frac{\partial S_{I}}{\partial p}) \frac{\partial p}{\partial t} + ε (S_{L} \frac{\partial ρ_{L}}{\partial T} + S_{I} \frac{\partial ρ_{I}}{\partial T}) \frac{\partial T}{\partial t} + ε (ρ_{L} \frac{\partial S_{L}}{\partial T} + ρ_{I} \frac{\partial S_{I}}{\partial T}) \frac{\partial T}{\partial t}}} - \\ \underline{\nabla} \cdot [(\frac{\underline{\underline{k}} k_{r} ρ_{L}}{μ}) (\underline{\nabla} p - ρ_{L} \underline{g})] = Q_{p} \end{array}

and

(73)

\begin{array}{l} \overset{STP water-mass change rate f (p)}{\overset{︷}{\underset{compressive f (p)}{\underset{⎵}{\underset{SCLP}{\underset{︸}{S_{L} ρ_{L} S_{op}^{L} \frac{\partial p}{\partial t}}} + \underset{SCIP}{\underset{︸}{S_{I} ρ_{I} S_{op}^{I} \frac{\partial p}{\partial t}}}}} + \underset{saturation f (p)}{\underset{⎵}{\underset{SLP}{\underset{︸}{ε ρ_{L} \frac{\partial S_{L}}{\partial p} \frac{\partial p}{\partial t}}} + \underset{SIP}{\underset{︸}{ε ρ_{I} \frac{\partial S_{I}}{\partial p} \frac{\partial p}{\partial t}}}}}}} + \overset{STU water-mass change rate f (T)}{\overset{︷}{\underset{density f (T)}{\underset{⎵}{\underset{SDLU}{\underset{︸}{ε S_{L} \frac{\partial ρ_{L}}{\partial T} \frac{\partial T}{\partial t}}} + \underset{SDIU}{\underset{︸}{ε S_{I} \frac{\partial ρ_{I}}{\partial T} \frac{\partial T}{\partial t}}}}} + \underset{saturation f (T)}{\underset{⎵}{\underset{SLU}{\underset{︸}{ε ρ_{L} \frac{\partial S_{L}}{\partial T} \frac{\partial T}{\partial t}}} + \underset{SIU}{\underset{︸}{ε ρ_{I} \frac{\partial S_{I}}{\partial T} \frac{\partial T}{\partial t}}}}}}} - \\ \underline{\nabla} \cdot [(\frac{\underline{\underline{k}} k_{r} ρ_{L}}{µ}) (\underline{\nabla} p - ρ_{L} \underline{g})] = \underset{QIN}{\underset{︸}{Q_{p}}} \end{array}

,

(74)

where

f(p): indicates dependence on pressure changes, and
f(T): indicates dependence on temperature changes.

Note that ice density values in SIP and SIU terms of equations 73 and 74 are set equal to liquid density values, and, in the SDIU term, the ice density derivative with respect to temperature is set to the value of the liquid density derivative with respect to temperature; see section 1.1.1, “Assumption Regarding Pressure Effect of Liquid-Ice Density Difference in Water-Mass Balance.” In SUTRA, all flows in and out of the model domain are conceptualized as occurring through the boundary conditions. Therefore, the divergence term ( $\underline{\nabla} \cdot$ ) accounts only for internal redistributions of water within the model and does not contribute to the overall water budget.

The enhanced water-mass budget output has the following organization in the .lst file (each budget term in the above equations may be found below by locating its source code variable name):

H 2 O  M A S S  B U D G E T  (MASS/SECOND)
.............................................................................
RATE OF CHANGE OF H2O STORED IN MODEL DUE TO
PRESSURE CHANGE:
LIQUID-COMPRESSIVE STORAGE —SCLP—
LIQUID SATURATION —SLP—
ICE-COMPRESSIVE STORAGE —SCIP—
ICE SATURATION —SIP—
TOTAL RATE OF CHANGE OF H2O IN MODEL DUE TO PRESSURE CHANGE —STP—
.............................................................................
RATE OF CHANGE OF H2O STORED IN MODEL DUE TO
TEMPERATURE CHANGE:
LIQUID DENSITY —SDLU—
LIQUID SATURATION —SLU—
ICE DENSITY —SDIU—
ICE SATURATION —SIU—
TOTAL RATE OF CHANGE OF H2O IN MODEL DUE TO TEMPERATURE CHANGE —STU—
.............................................................................
_________________________________________________________
TOTAL RATE OF CHANGE IN H2O STORED IN MODEL [ S+, S-, S] —STW—
.............................................................................
GAIN/LOSS OF LIQUID THROUGH INFLOWS/OUTFLOWS AT:
LIQUID SOURCES AND SINKS	(sum of all —QIN— values in model)
SPECIFIED PRESSURE NODES	(sum of all —QPL— values in model)
GENERALIZED-FLOW NODES	(sum of all —QPG— values in model)
_____________________________________________________________
TOTAL RATE OF GAIN/LOSS OF LIQUID THROUGH FLOWS [ F+, F-, F] —QFF—
.............................................................................
H2O MASS BALANCE ACTIVITY [ A=((S+) - (S-)+(F+) - (F-))/2]
ABSOLUTE H2O MASS BALANCE ERROR [ S - F]
RELATIVE H2O MASS BALANCE ERROR [ 100*(S - F)/A]	(in percent)
.............................................................................
LIQUID GAINS OR LOSSES DUE TO LIQUID SOURCE NODES	(output only for time-dependent sources; list of nodal values)
LIQUID GAINS OR LOSSES DUE TO SPECIFIED PRESSURE NODES	(list of nodal values)
LIQUID GAINS OR LOSSES DUE TO GENERALIZED-FLOW NODES	(list of nodal values)
.............................................................................

The terms S+, S-, S, and F+, F-, F, are positive, negative, and net contributions to the rate of fluid-mass change in the entire model domain, respectively, with values listed in columns along the right side of the budget output in the .lst file. Water-mass gains and losses at specified pressure nodes and generalized-flow nodes are all contributions to the water-mass balance in the same manner as fluid sources and sinks (QIN) and are listed for every boundary-condition node at the end of each water-mass budget output.

4.1.2. Energy Budget

The energy budget is calculated using the divergent form of the energy-balance equation 30 (without subtraction of the fluid mass balance, which led to equation 35, the mathematically equivalent, convective form actually solved by SUTRA), because this form includes the total thermal energy of the system of most interest when interpreting results. The terms of the energy balance that are printed in the enhanced energy budget in the SUTRA .lst file are described in the following and are associated with Fortran variables used in the code for each term. Beginning with equation 30, by expanding terms using equations 11 and 17, and by gathering terms with similar physical origins, all energy budget terms are defined below. For each term in the equations, the information given provides the meaning of the term, and whether the value of the term depends on pressure changes, temperature changes, or both.

The manner in which the variables are explained here is exemplified by the Fortran source code variable STE. This value is printed in the energy budget output, representing the total energy change rate (given by the sum of terms on the left-hand side of equation 75, shown below) at the end of each user-selected completed time step. The STE value may be found in the following line of the energy budget output, which is shown in its entirety after equation 78 below:

TOTAL RATE OF CHANGE (R1+R2+R3) OF ENERGY STORED IN MODEL [S+,S-,S] —STE—

The variable name STE is highlighted as —STE— to indicate that this text is not actually written in the budget output and is included here for informational purposes only.

The equations that define the energy budget terms are as follows:

\begin{array}{l} \overset{STE total energy change rate f (p, T) ​}{\overset{︷}{[ε (S_{L} ρ_{L} c_{L} + S_{I} ρ_{I} c_{I}) + (1 - ε) ρ_{S} c_{S}] \frac{\partial T}{\partial t} + c_{L} T \frac{\partial}{\partial t} (ε S_{L} ρ_{L}) + (- Δ H + c_{I} T) \frac{\partial}{\partial t} (ε S_{I} ρ_{I})}} = \\ - \underline{\nabla} \cdot (ε S_{L} ρ_{L} c_{L} \underline{v} T) - \underline{\nabla} \cdot {\underline{q}}_{e} + Q_{p} c_{L} T^{*} + ε (S_{L} ρ_{L} γ_{0}^{L} + S_{I} ρ_{I} γ_{0}^{I}) + (1 - ε) ρ_{S} γ_{0}^{S} \end{array}

and

(75)

\begin{array}{l} \overset{STS sensible energy change f (T)}{\overset{︷}{\underset{FLD}{\underset{︸}{ε S_{L} ρ_{L} c_{L} \frac{\partial T}{\partial t}}} + \underset{FIU}{\underset{︸}{ε S_{I} ρ_{I} c_{I} \frac{\partial T}{\partial t}}} + \underset{SLD}{\underset{︸}{(1 - ε) ρ_{S} c_{S} \frac{\partial T}{\partial t}}}}} + \underset{DNS}{\underset{︸}{c_{L} T \frac{\partial}{\partial t} (ε S_{L} ρ_{L})}} + \underset{UIS}{\underset{︸}{(- Δ H + c_{I} T) \frac{\partial}{\partial t} (ε S_{I} ρ_{I})}} = \\ - \underline{\nabla} \cdot (ε S_{L} ρ_{L} c_{L} \underline{v} T) - \underline{\nabla} \cdot {\underline{q}}_{e} + \underset{QIU}{\underset{︸}{Q_{p} c_{L} T^{*}}} + \overset{PRS ​ energy production}{\overset{︷}{\underset{P0L}{\underset{︸}{ε S_{L} ρ_{L} γ_{0}^{L}}} + \underset{P0I}{\underset{︸}{ε S_{I} ρ_{I} γ_{0}^{I}}} + \underset{P0S}{\underset{︸}{(1 - ε) ρ_{S} γ_{0}^{S}}}}} \end{array}

.

(76)

In the SUTRA budget calculations, the terms DNS and UIS defined in equation 76 are further expanded using equations 11 and 17 as shown below:

\overset{\begin{array}{l} DNS ​ energy change due to \\ liquid mass change f (p, T) \end{array}}{\overset{︷}{c_{L} T \frac{\partial (ε S_{L} ρ_{L})}{\partial t}}} = \overset{DNP ​ liquid mass f (p)}{\overset{︷}{\underset{USCP compressibility}{\underset{︸}{c_{L} T ρ_{L} (S_{L} S_{op}^{L}) \frac{\partial p}{\partial t}}} + \underset{USLP ​ liquid saturation}{\underset{︸}{c_{L} T ρ_{L} (ε \frac{\partial S_{L}}{\partial p}) \frac{\partial p}{\partial t}}}}} + \overset{DNU ​ liquid mass f (T)}{\overset{︷}{\underset{USDL liquid density}{\underset{︸}{c_{L} T ε (S_{L} \frac{\partial ρ_{L}}{\partial T}) \frac{\partial T}{\partial t}}} + \underset{USLU ​ liquid saturation}{\underset{︸}{c_{L} T ε (ρ_{L} \frac{\partial S_{L}}{\partial T}) \frac{\partial T}{\partial t}}}}}

and

(77)

\begin{array}{l} \overset{UIS latent plus sensible energy change due to ice mass change f (p, T)}{\overset{︷}{(- Δ H + c_{I} T) \frac{\partial (ε S_{I} ρ_{I})}{\partial t}}} = \\ \overset{UIP ice mass f (p)}{\overset{︷}{\underset{USCIP ice compressibility}{\underset{︸}{(- Δ H + c_{I} T) (ρ_{I} S_{I} S_{op}^{I}) \frac{\partial p}{\partial t}}} + \underset{FLP ice saturation}{\underset{︸}{(- Δ H + c_{I} T) (ε ρ_{I} \frac{\partial S_{I}}{\partial p}) \frac{\partial p}{\partial t}}}}} + \overset{UIU ice mass f (T)}{\overset{︷}{\underset{USDI ice density}{\underset{︸}{(- Δ H + c_{I} T) ε (S_{I} \frac{\partial ρ_{I}}{\partial T}) \frac{\partial T}{\partial t} +}} \underset{FLU ice saturation}{\underset{︸}{(- Δ H + c_{I} T) ε (ρ_{I} \frac{\partial S_{I}}{\partial T}) \frac{\partial T}{\partial t}}}}} \end{array}

.

(78)

The ice-density values in FLP and FLU terms are set equal to liquid-density values; and, in the USDI term, the ice density derivative with respect to temperature is set to the value of the liquid density derivative with respect to temperature; see section 1.1.1, “Assumption Regarding Pressure Effect of Liquid-Ice Density Difference in Water-Mass Balance.” The quantity (−∆H+c_IT) is the energy content of ice per unit mass of ice (defined in equation 28B and in the source code it has the variable name EICE).

The enhanced energy-budget output has the following organization in the .lst file. Each budget term in the above equations may be found below by locating its variable name:

E N E R G Y  B U D G E T (ENERGY/SECOND)
.........................................................................
RATE OF CHANGE OF ENERGY STORED IN MODEL DUE TO
TEMPERATURE CHANGE IMPACTS ON THERMAL ENERGY IN:
LIQUID —FLD—
ICE —FIU—
SOLID GRAINS —SLD—
__________________________________________________________
TOTAL RATE OF ENERGY CHANGE DUE TO TEMPERATURE CHANGE (R1) —STS—
.........................................................................
RATE OF CHANGE OF ENERGY STORED IN MODEL DUE TO
CHANGE IN MASS OF LIQUID CAUSED BY PRESSURE CHANGE IMPACTS ON:
LIQUID COMPRESSIVE STORAGE —USCP—
LIQUID SATURATION —USLP—
TOTAL RATE OF ENERGY CHANGE DUE TO PRESSURE AND LIQUID MASS CHANGE —DNP—
RATE OF CHANGE OF ENERGY STORED IN MODEL DUE TO
CHANGE IN MASS OF LIQUID CAUSED BY TEMPERATURE CHANGE IMPACTS ON:
LIQUID DENSITY —USDL—
LIQUID SATURATION —USLU—
TOTAL RATE OF ENERGY CHANGE DUE TO TEMP. AND LIQUID MASS CHANGE —DNU—
_____________________________________________________________
TOTAL RATE OF CHANGE OF ENERGY DUE TO LIQUID MASS CHANGE (R2) —DNS—
.........................................................................
RATE OF CHANGE OF ENERGY STORED IN MODEL DUE TO
CHANGE IN MASS OF ICE CAUSED BY PRESSURE CHANGE IMPACTS ON:
ICE COMPRESSIVE STORAGE —USCIP—
ICE SATURATION —FLP—
TOTAL RATE OF ENERGY CHANGE DUE TO PRESSURE AND ICE MASS CHANGE —UIP—
RATE OF CHANGE OF ENERGY STORED IN MODEL DUE TO
CHANGE IN MASS OF ICE CAUSED BY TEMPERATURE CHANGE IMPACTS ON:
ICE DENSITY —USDI—
ICE SATURATION —FLU—
(INCLUDES LATENT HEAT OF FREEZE-THAW PHASE CHANGE)
TOTAL RATE OF ENERGY CHANGE DUE TO TEMP. AND ICE MASS CHANGE —UIU—
__________________________________________________________
TOTAL RATE OF CHANGE OF ENERGY DUE TO ICE MASS CHANGE (R3) —UIS—
.........................................................................
_______________________________________________________________________
TOTAL RATE OF CHANGE (R1+R2+R3) OF ENERGY STORED IN MODEL [ S+, S-, S] —STE—
.........................................................................
RATE OF CHANGE OF ENERGY STORED IN MODEL DUE TO
ZERO-ORDER PRODUCTION/DECAY OF ENERGY IN:
LIQUID —P0L—
ICE —P0I—
SOLID GRAINS —P0S—
______________________________________________________
TOTAL RATE OF PRODUCTION/DECAY OF ENERGY [ P+, P-, P] —PRS—
.........................................................................
GAIN/LOSS OF ENERGY
THROUGH INFLOWS/OUTFLOWS AT:
LIQUID SOURCES AND SINKS —QIU—
SPECIFIED PRESSURE NODES —QPU—
GENERALIZED-FLOW NODES —QPGU—
AT:
ENERGY SOURCES AND SINKS —QQU—
SPECIFIED TEMPERATURE NODES —QUL—
GENERALIZED-TRANSPORT NODES —QUG—
__________________________________________________________
TOTAL RATE OF GAIN/LOSS OF ENERGY —QFS—
THROUGH INFLOWS/OUTFLOWS & SOURCES/SINKS [ F+, F-, F]
.............................................................................
ENERGY BALANCE ACTIVITY [ A=((S+) - (S-)+(P+) - (P-)+(F+) - (F-))/2]
ABSOLUTE ENERGY BALANCE ERROR [ S - P - F]
RELATIVE ENERGY BALANCE ERROR [ 100*(S - P - F)/A]	(in percent)
.............................................................................
ENERGY SOURCES OR SINKS AT FLUID SOURCES AND SINKS	(list of nodal values)
ENERGY SOURCES OR SINKS DUE TO FLUID INFLOWS OR OUTFLOWS AT POINTS OF SPECIFIED PRESSURE	(list of nodal values)
ENERGY SOURCES OR SINKS DUE TO FLUID INFLOWS OR OUTFLOWS AT POINTS OF GENERALIZED FLOW	(list of nodal values)
TIME-DEPENDENT ENERGY SOURCES OR SINKS SET IN SUBROUTINE BCTIME	(list of nodal values)
ENERGY SOURCES OR SINKS DUE TO SPECIFIED TEMPERATURES	(list of nodal values)
ENERGY SOURCES OR SINKS DUE TO GENERALIZED-TRANSPORT CONDITIONS	(list of nodal values)
.............................................................................

The terms S+, S-, S, P+, P-, P, and F+, F-, F, are positive, negative, and net contributions to the rate of energy change in the entire model domain, respectively, with values listed in columns along the right side of the budget output in the .lst file. Energy gains and losses at liquid sources and sinks, specified pressure nodes, specified temperature nodes, and generalized-flow nodes are all contributions to the energy balance in the same manner as energy sources and sinks (QUIN) and are listed for every boundary-condition node at the end of each budget output.

4.1.3. Water-Mass Budget for Solute-Transport Simulations

The components of the water-mass budget output depend on whether SUTRA is simulating energy transport or solute transport. In this section, the enhanced water-mass budget for solute-transport simulations is presented. The manner in which the variables are explained here is exemplified by the source code variable STW, which is shown here for the SUTRA water-mass-balance equation (eq. 4) rewritten for simulation of solute-mass transport, in which the ice saturation is set to zero and temperature is replaced by solute mass fraction, C:

\overset{STW total water-mass change rate f (p, C)}{\overset{︷}{\frac{\partial (ε S_{L} ρ_{L})}{\partial t}}} = - \underline{\nabla} \cdot (ε S_{L} ρ_{L} \underline{v}) + Q_{p}

.

(79)

This value represents the total water-mass change rate, given by the first term of equation 79. The text that accompanies the label “STW” in equation 79 provides the meaning of this value and indicates that STW depends on changes in both p and C. Each term of the enhanced water-mass budget is defined in an equation below and, with each term, an indication is given of whether each rate depends on pressure changes, concentration changes, or both. The value of STW is printed in the water-mass budget output at the end of each user-selected completed time step. STW may be found in the following line of the budget output, which is shown in its entirety after equation 80 below:

TOTAL RATE OF CHANGE IN H2O STORED IN MODEL [ S+, S-, S] —STW—

The variable name STW is highlighted as —STW— to indicate that this text is not actually written in the budget output and is included here for informational purposes only.

The other terms in the enhanced water-mass balance output for solute transport are obtained from an expanded form of equation 79, through a development analogous to the one used to obtain equation 73, but for solute transport. This provides the definitions of the other water-mass budget for solute-transport values:

\overset{STW total water-mass change rate f (p, C)}{\overset{︷}{(S_{L} ρ_{L} S_{op}^{L}) \frac{\partial p}{\partial t} + ε (ρ_{L} \frac{\partial S_{L}}{\partial p}) \frac{\partial p}{\partial t} + ε (S_{L} \frac{\partial ρ_{L}}{\partial C}) \frac{\partial C}{\partial t}}} - \underline{\nabla} \cdot [(\frac{\underline{\underline{k}} k_{r} ρ_{L}}{μ}) (\underline{\nabla} p - ρ_{L} \underline{g})] = Q_{p}

and

(80)

\overset{STP water-mass change rate f (p)}{\overset{︷}{\underset{SCLP compressive f (p) ​}{\underset{︸}{S_{L} ρ_{L} S_{op}^{L} \frac{\partial p}{\partial t}}} + \underset{SLP saturation f (p)}{\underset{︸}{ε ρ_{L} \frac{\partial S_{L}}{\partial p} \frac{\partial p}{\partial t}}}}} + \overset{STU water-mass change rate f (C)}{\overset{︷}{\underset{SDLU density f (C)}{\underset{︸}{ε S_{L} \frac{\partial ρ_{L}}{\partial C} \frac{\partial C}{\partial t}}}}} - \underline{\nabla} \cdot [(\frac{\underline{\underline{k}} k_{r} ρ_{L}}{μ}) (\underline{\nabla} p - ρ_{L} \underline{g})] = \underset{QIN}{\underset{︸}{Q_{p}}}

,

(81)

where

f(C): indicates dependence on concentration changes and
f(p,C): indicates dependence on both pressure changes and concentration changes.

In SUTRA, all flows in and out of the model domain are conceptualized as occurring through the boundary conditions. Therefore, the divergence term ( $\underline{\nabla} \cdot$ ) accounts only for internal redistributions of water within the model and does not contribute to the overall water budget.

The enhanced water-mass budget output for solute transport has the following organization in the .lst file. Each budget term in the above equations may be found below by locating its variable name in the source code:

H 2 O  M A S S  B U D G E T  (MASS/SECOND)
.............................................................................
RATE OF CHANGE OF H2O STORED IN MODEL DUE TO
PRESSURE CHANGE IMPACTS ON:
LIQUID-COMPRESSIVE STORAGE —SCLP—
LIQUID SATURATION —SLP—
TOTAL RATE OF CHANGE OF H2O IN MODEL DUE TO PRESSURE CHANGE —STP—
.............................................................................
RATE OF CHANGE OF H2O STORED IN MODEL DUE TO
CONCENTRATION CHANGE IMPACTS ON:
LIQUID DENSITY —SDLU—
TOTAL RATE OF CHANGE OF H2O IN MODEL DUE TO CONCENTRATION CHANGE —STU—
.............................................................................
_________________________________________________________
TOTAL RATE OF CHANGE IN H2O STORED IN MODEL [ S+, S-, S ] —STW—
.............................................................................
GAIN/LOSS OF LIQUID THROUGH INFLOWS/OUTFLOWS AT:
LIQUID SOURCES AND SINKS	(sum of all —QIN— values in model)
SPECIFIED PRESSURE NODES	(sum of all —QPL— values in model)
GENERALIZED-FLOW NODES	(sum of all —QPG— values in model)
_____________________________________________________________
TOTAL RATE OF GAIN/LOSS OF LIQUID THROUGH FLOWS [ F+, F-, F ] —QFF—
.............................................................................
H2O MASS BALANCE ACTIVITY [ A=((S+) - (S-)+(F+) - (F-))/2 ]
ABSOLUTE H2O MASS BALANCE ERROR [ S - F ]
RELATIVE H2O MASS BALANCE ERROR [ 100*(S - F)/A ]	(in percent)
.............................................................................
LIQUID GAINS OR LOSSES DUE TO LIQUID SOURCE NODES	(output only for time-dependent sources; list of nodal values)
LIQUID GAINS OR LOSSES DUE TO SPECIFIED PRESSURE NODES	(list of nodal values)
LIQUID GAINS OR LOSSES DUE TO GENERALIZED-FLOW NODES	(list of nodal values)
.............................................................................

The terms S+, S-, S, and F+, F-, F, are positive, negative, and net contributions to the rate of fluid-mass change in the entire model domain, respectively, with values listed in columns along the right side of the budget output in the .lst file. Water-mass gains and losses at specified pressure nodes and generalized-flow nodes are all contributions to the water-mass balance in the same manner as fluid sources and sinks (QIN) and are listed for every boundary-condition node at the end of each water-mass budget output.

4.1.4. Solute-Species Mass Budget

The solute-species mass budget is calculated using the full solute-species mass-balance equation (Voss and Provost, 2002, eq. 2.3, which is equivalent to Voss and Provost, 2002, eq. 2.47), the unified transport equation for solute-species mass or energy in terms of unified transport variable, U, given in terms of concentration. This balance is informally termed “solute transport equation,” but it is actually a complete mass balance of this chemical species within the porous medium, because it also allows simulation of the same chemical species when this species is immobile, while adsorbed on the surface or interior of solid mineral grains. This version of the solute transport equation, used for the solute-species mass-balance budget output, is without subtraction of the fluid-mass balance, because this form includes the total mass of the chemical species being tracked in the system, which is of most interest when interpreting results. [Subtraction of the fluid mass balance leads to the mathematically equivalent form actually solved for solute concentration by SUTRA (Voss and Provost, 2002, eq. 2.52), as it did for energy transport, resulting in equation 34.]

The terms of the solute-species mass balance that are printed in the enhanced solute mass budget in the SUTRA .lst file are described in the following and are associated with variables used in the Fortran source code for each term. Beginning with equation 30, by expanding one term using equation 11, both written for solute-species mass in terms of concentration rather than temperature, and by gathering terms with similar physical origins, all solute-species mass-budget terms are defined below. For each term in the equations, the information given provides the meaning of the term and whether the value of the term depends on pressure changes, solute and sorbate concentration changes, or both pressure and concentration.

The SUTRA solute-species mass-balance equation is equivalent to equation 30 rewritten for simulation of solute-mass transport of a particular chemical species, in which the ice saturation is set to zero (removing several terms from the equation); temperature is replaced by solute mass fraction, C; the specific heat of liquid, c_L, is set to a value of one; the conductive-dispersive heat flux, ${\underline{q}}_{e}$ , is replaced by its solute-species mass equivalent, ${\underline{q}}_{S}$ , the diffusive-dispersive solute flux; zero-order energy production is replaced by zero-order solute-species mass production in the water and in or on the solid mineral grains (for example, in any sorbate), and additional terms for first-order production of the solute-species in the water and in the immobile (sorbed) phase are added (see Voss and Provost, 2002, eq. 2.37):

C(x, y[, z], t): [M_S/M] is solute-species concentration (mass fraction);
C_S(x, y[, z], t): [M_S/M_G] is concentration of species as sorbate on solid grains expressed as a mass fraction (mass adsorbate per mass solid grains plus adsorbate);
C*(x, y[, z], t): [M_S/M] is solute concentration of fluid sources (mass fraction);
${\underline{q}}_{S} (x, y [, z], t)$: [M/(L²‧s)] is diffusive-dispersive solute-mass flux;
$γ_{0}^{L} (x, y [, z], t)$: [s⁻¹] is zero-order solute-species mass production in liquid;
$γ_{0}^{S} (x, y [, z], t)$: [s⁻¹] is zero-order solute-species mass production in solid mineral grains;
$γ_{1}^{L} (x, y [, z], t)$: [(M_S/M)/s)] is first-order solute-species mass production in liquid; and
$γ_{1}^{S} (x, y [, z], t)$: [(M_S/M_G)/s)] is first-order solute-species mass production in solid mineral grains.

M_S, M_G, and M are the units for mass of solute species, solid grains, and water, respectively. These parameters are fully explained in Voss and Provost (2002).

The manner in which the variables of the solute-mass balance appear in the enhanced budget output are explained here, as exemplified by the Fortran source code variable STE. The STE value, which represents the total solute-species mass change rate (given by the sum of terms on the left-hand side of equation 82, shown below), is printed in the solute-species mass budget output at the end of each user-selected completed time step. Note that the terms “solute-species mass” and “solute mass” are used interchangeably in the budget output and mean the same thing. The STE value may be found in the following line of the budget output, which is shown in its entirety after equation 84 below:

TOTAL RATE OF CHANGE (R1+R2+R3) OF SOLUTE MASS STORED IN MODEL [S+,S-,S] —STE—

The variable name STE is highlighted as —STE— to indicate that this text is not actually written in the budget output and is included here for informational purposes only.

The equations that define the solute-species mass-budget terms are as follows (from Voss and Provost (2002), equations 2.31 and 2.47):

\overset{STE total solute-species mass change rate f (p, C, C_{S})}{\overset{︷}{\frac{\partial}{\partial t} (ε S_{L} ρ_{L} C) + \frac{\partial}{\partial t} [(1 - ε) ρ_{S} C_{S}]}} = - \underline{\nabla} \cdot (ε S_{L} ρ_{L} \underline{v} C) - \underline{\nabla} \cdot {\underline{q}}_{S} + Q_{P} C^{*} + ε S_{L} ρ_{L} (γ_{1}^{L} C + γ_{0}^{L}) + (1 + ε) ρ_{S} (γ_{1}^{S} C + γ_{0}^{S})

and

(82)

\begin{array}{l} \overset{STS solute-species mass change due to change in C and C_{S} f (C, C_{S}) ​}{\overset{︷}{\underset{FLD}{\underset{︸}{ε S_{L} ρ_{L} \frac{\partial C}{\partial t}}} + \underset{DNS}{\underset{︸}{C \frac{\partial}{\partial t} (ε S_{L} ρ_{L})}} + \underset{SLD}{\underset{︸}{(1 - ε) ρ_{S} \frac{\partial C_{S}}{\partial t}}}}} = \\ - \underline{\nabla} \cdot (ε S_{L} ρ_{L} \underline{v} C) - \underline{\nabla} \cdot {\underline{q}}_{S} + \underset{QIU}{\underset{︸}{Q_{p} C^{*}}} + \overset{PRS solute-species mass production}{\overset{︷}{\underset{P1L}{\underset{︸}{ε (S_{L} ρ_{L} γ_{1}^{L} C)}} + \underset{P1S}{\underset{︸}{(1 - ε) ρ_{S} γ_{1}^{S} C_{S}}} + \underset{P0L}{\underset{︸}{ε (S_{L} ρ_{L} γ_{0}^{L})}} + \underset{P0S}{\underset{︸}{(1 - ε) ρ_{S} γ_{0}^{S}}}}} \end{array}

,

(83)

where

f(C, C_S): indicates dependence on solute and sorbate concentration changes and
f(p, C, C_S): indicates dependence on both pressure changes and concentration changes.

In the SUTRA budget calculations, the DNS term defined in equation 83 is further expanded using the solute-mass equivalent of equation 11 as shown below:

\overset{DNS solute-species mass change due to liquid mass change f (p, C)}{\overset{︷}{C \frac{\partial (ε S_{L} ρ_{L})}{\partial t}}} = \overset{DNP l i q u i d m a s s c h a n g e f (p)}{\overset{︷}{\underset{USCP ​ compressibility}{\underset{︸}{C ρ_{L} (S_{L} S_{op}^{L}) \frac{\partial p}{\partial t}}} + \underset{USLP liquid saturation}{\underset{︸}{C ρ_{L} (ε \frac{\partial S_{L}}{\partial p}) \frac{\partial p}{\partial t}}}}} + \overset{DNU l i q u i d m a s s c h a n g e f (C)}{\overset{︷}{\underset{USDL liquid density}{\underset{︸}{C ε (S_{L} \frac{\partial ρ_{L}}{\partial C}) \frac{\partial C}{\partial t}}}}}

,

(84)

where f(p) indicates dependence on pressure changes.

The enhanced solute-mass-budget output has the following organization in the .lst file. Each budget term in the above equations may be found below by locating its Fortran variable name:

S O L U T E  B U D G E T  I N (SOLUTE MASS/SECOND)
...........................................................................
RATE OF CHANGE OF SOLUTE STORED IN MODEL DUE TO
CONCENTRATION CHANGE IMPACTS ON SOLUTE MASS IN:
LIQUID —FLD—
ADSORBATE —SLD—
____________________________________________________________
TOTAL RATE OF SOLUTE CHANGE DUE TO CONCENTRATION CHANGE (R1) —STS—
...........................................................................
RATE OF CHANGE OF SOLUTE STORED IN MODEL DUE TO
CHANGE IN MASS OF LIQUID CAUSED BY PRESSURE CHANGE:
LIQUID COMPRESSIVE STORAGE —USCP—
LIQUID SATURATION —USLP—
TOTAL RATE OF SOLUTE CHANGE DUE TO PRESSURE AND LIQUID MASS CHANGE —DNP—
RATE OF CHANGE OF SOLUTE STORED IN MODEL DUE TO
CHANGE IN MASS OF LIQUID CAUSED BY CONCENTRATION CHANGE:
LIQUID DENSITY —USDL—
TOTAL RATE OF SOLUTE CHANGE DUE TO CONC. AND LIQUID MASS CHANGE —DNU—
_____________________________________________________________
TOTAL RATE OF CHANGE OF SOLUTE DUE TO LIQUID MASS CHANGE (R2) —DNS—
...........................................................................
_______________________________________________________________________
TOTAL RATE OF CHANGE (R1+R2+R3) OF SOLUTE STORED IN MODEL [ S+, S-, S] —STE—
...........................................................................
RATE OF CHANGE OF SOLUTE STORED IN MODEL DUE TO
FIRST-ORDER PRODUCTION/DECAY OF SOLUTE IN:
LIQUID —P1L—
SOLID GRAINS —P1S—
ZERO-ORDER PRODUCTION/DECAY OF SOLUTE IN:
LIQUID —P0L—
SOLID GRAINS —P0S—
______________________________________________________
TOTAL RATE OF PRODUCTION/DECAY OF SOLUTE [ P+, P-, P] —PRS—
...........................................................................
GAIN/LOSS OF SOLUTE
THROUGH INFLOWS/OUTFLOWS AT:
LIQUID SOURCES AND SINKS —QIU—
SPECIFIED PRESSURE NODES —QPU—
GENERALIZED-FLOW NODES —QPGU—
THROUGH INFLUXES/OUTFLUXES AT:
SOLUTE SOURCES AND SINKS —QQU—
SPECIFIED CONCENTRATION NODES —QUL—
GENERALIZED-TRANSPORT NODES —QUG—
_______________________________________________________________
TOTAL RATE OF GAIN/LOSS OF SOLUTE
THROUGH INFLOWS/OUTFLOWS & INFLUXES/OUTFLUXES [ F+, F-, F]
...........................................................................
SOLUTE BALANCE ACTIVITY [ A=((S+) - (S-)+(P+) - (P-)+(F+) - (F-))/2
ABSOLUTE SOLUTE BALANCE ERROR [ S - P - F]
RELATIVE SOLUTE MASS BALANCE ERROR [ 100*(S - P - F)/A]
...........................................................................
SOLUTE SOURCES OR SINKS
AT FLUID SOURCES AND SINKS	(list of nodal values)
SOLUTE SOURCES OR SINKS DUE TO FLUID INFLOWS OR OUTFLOWS
AT POINTS OF SPECIFIED PRESSURE	(list of nodal values)
SOLUTE SOURCES OR SINKS DUE TO FLUID INFLOWS OR OUTFLOWS
AT POINTS OF GENERALIZED FLOW	(list of nodal values)
TIME-DEPENDENT SOLUTE SOURCES OR SINKS SET IN SUBROUTINE BCTIME	(list of nodal values)
SOLUTE SOURCES OR SINKS DUE TO SPECIFIED TEMPERATURES		(list of nodal values)
SOLUTE SOURCES OR SINKS DUE TO GENERALIZED-TRANSPORT CONDITIONS	(list of nodal values)
.............................................................................

The terms: S+, S-, S, P+, P-, P, and F+, F-, F, are positive, negative, and net contributions to the rate of solute change in the entire model domain, respectively, with values listed in columns along the right side of the budget output in the .lst file. Solute gains and losses at liquid sources and sinks, specified pressure nodes, specified temperature nodes, and generalized-flow nodes are all contributions to the solute balance in the same manner as solute sources and sinks (QUIN) and are listed for every boundary-condition node at the end of each budget output.

4.2. Enhanced Velocity Output: Darcy Velocity

The output of fluid-velocity magnitude and direction at the center of each finite element, which SUTRA has provided in all previous versions, allows visualization of the groundwater flow field. The velocity values output by SUTRA are calculated with equation 18 (repeated here):

\underline{v} = - (\frac{\underline{\underline{k}} k_{r}}{ε S_{L} μ}) (\underline{\nabla} p - ρ_{L} \underline{g})

.

(85)

This is Darcy’s Law, generalized for variably saturated, variable-density flow and written in terms of fluid pressure, as explained in Voss and Provost (2002, section 2.2, “Saturated-Unsaturated Groundwater Flow”).

The velocity in equation 85, $\underline{v}$ [L/s], is the average velocity at which a conservative tracer would move through a porous medium, also called the “tracer velocity” and “particle velocity.” Thus, $\underline{v}$ is a measure of the average bulk movement of the liquid-water molecules and nonreactive substances within the liquid, relative to the stationary solid matrix. The velocity $\underline{v}$ is an “average” because it aggregates velocities that vary from point to point within a porous medium due to variations in permeability, porosity, and connectivity at a spatial scale smaller than that at which typical laboratory and field measurements are made.

Visualizations of $\underline{v}$ provide insight into the speed and direction of water movement but not about the quantity of water that moves through the porous medium being modeled. Interpreting the quantity of water movement from visualizations of $\underline{v}$ can sometimes be misleading, as explained in the following.

Particularly for simulations of unsaturated and freezing systems, velocity magnitudes may decrease, remain unchanged, or even may increase when the porous medium desaturates or freezes. The expectation that fluid fluxes change in the same manner as fluid velocity magnitudes might contradict the intuitive expectation that total fluid flux should decrease during desaturation or freezing. This range of possible decreases or increases occurs because the change in fluid velocity resulting from a decrease in liquid saturation, S_L, depends on the ratio of relative permeability to liquid saturation, k_r/S_L, according to equation 85. If both S_L and k_r decrease at the same rate when S_L decreases due to desaturation or freezing (as defined by the user-selected liquid saturation and relative permeability relations; see section 2, “Saturation and Relative-Permeability Functions”), the velocity magnitude will remain constant as S_L decreases. Should S_L decrease faster than k_r, velocity magnitude will increase as S_L decreases, although total fluid flux always decreases as S_L decreases.

Such fluid velocity behavior may occur in a simulation, depending on the user-selected saturation and relative permeability functions, and depending on what these functions imply about the connectivity and volume of flowing liquid zones within the porous medium. Only if k_r decreases faster than S_L during desaturation or freezing will velocity magnitudes decrease. As a result, the velocity output does not provide insight into the quantity of fluid (that is, the fluid flux) that is flowing during desaturation or freezing; it only shows how fast the fluid is moving.

The rate at which fluid volume is moving is instead given by the “Darcy velocity”:

\underline{q} = ε S_{L} \underline{v} = - (\frac{\underline{\underline{k}} k_{r}}{μ}) (\underline{\nabla} p - ρ_{L} \underline{g})

,

(86)

where

\underline{q}

[L/s] is the Darcy velocity or volumetric fluid flux [(fluid volume per time per cross-sectional area of total porous medium)/s].

Because the porosity is always less than one, the Darcy velocity, $\underline{q}$ , is always less than the average fluid velocity, $\underline{v}$ , and so is not a direct indicator of the speed of water movement. Rather, Darcy velocity is a flux of fluid, representing the volume of fluid flowing across a unit cross-sectional area of the porous medium, per unit time. Thus, the Darcy velocity gives information on the quantity of fluid flow, which may be expected to decrease during desaturation or freezing. Darcy velocity (eq. 86), the volumetric fluid flux, always decreases as liquid saturation, S_L, decreases, irrespective of whether the saturation decrease is due to desaturation or freezing.

Therefore, to allow fuller interpretation of SUTRA unsaturated and freeze-thaw simulation results, the current version of SUTRA additionally provides output of Darcy velocities, $\underline{q}$ , at the centroid of each finite element, at the end of user-selected time steps. Both Darcy velocity, $\underline{q}$ , and particle velocity, $\underline{v}$ , outputs are provided in SUTRA output files (with the extensions .lst and .ele).

Note: The velocities that are output (in the .lst and .ele output files) are at the centroids of each element in the finite-element mesh. These output values are provided solely for user information and postprocessing purposes. Internally, SUTRA uses only velocities at Gauss points within each element (not at element centroids) to solve the governing differential equations (see Voss and Provost, 2002, sec. 4.3, “Gaussian Integration”). The mean value of the velocities at the Gauss points in each finite element provides the value reported at the element centroid.

5. Example Simulations

Two simple examples are described in this section. These demonstrate the main new physics functionality provided by SUTRA version 4.0, highlighting saturated and unsaturated freezing and thawing.

5.1. Saturated Groundwater Freezing Example—Containment by Frozen Groundwater Wall

This example demonstrates the effect of freezing on the flow field in a groundwater system. The frozen wall is an engineered intervention that creates a region of low permeability upstream of a contaminated region of an aquifer. Shielding the contaminated area from regional groundwater flow avoids the spread of contaminants to a larger area, thereby allowing the contaminated area to be remediated. A similar “frozen wall problem” was first proposed by McKenzie, Voss, and Siegel (2007).

The model domain considered for this analysis is a 25-by-15-meter (m) area in which a hypothetical contaminant exists within a 4-by-4-m square region centered along the southern edge of the area (fig. 5). Initially (fig. 6A) the entire shallow, uniform aquifer is at a temperature of +5 °C and regional groundwater flow with the same temperature occurs uniformly from west to east, driven by the difference in hydraulic heads along the eastern and western boundaries (implemented in the model as specified pressure conditions). An L-shaped cooling system 0.5 m wide is emplaced around the contaminated zone by burying freezing coils through the entire representative 1-m thickness of the shallow aquifer. The cooling system is installed just upstream and along the lateral side of the contaminated area, at a distance of 0.5 m from the contaminated zone, to avoid entraining contaminants within the coils during emplacement. The cooling system extends 5 m in the north-south and 5 m in the east-west directions. The cooling system is implemented in the simulation in a simple manner as a zone of specified temperature conditions held at –5 °C for the duration of the simulation once the system is activated just after 0 hours.

Contaminated is contained by cooling wall with regional flow around it — Figure 5.
Illustration showing setup of “containment by frozen groundwater wall” saturated freezing example. This is an areal groundwater flow system, viewed from above. Observation points are located at y=0 meter (m) and x=9, 13, 17, 21, 25 m. See appendix 1 for a list of symbols.

Ice saturation curvature expands around the contaminated block as time passes — Figure 6.
Illustrations showing simulated evolution of ice saturation (S_I) for “containment by frozen groundwater wall” example. Colored contours show S_I values (S_I=1, represents completely frozen groundwater). Velocity vectors that show flow direction and velocity magnitude are represented in a “windsock” style: the direction of flow is indicated by the direction of the line that emanates from the dot, the dot indicates the location of the velocity vector, and the line length is proportional to the magnitude of the velocity. Velocity vectors are plotted at centroid of every second finite element in each direction, and only vectors with the three highest orders of magnitude are shown. Panel A (0 hours) shows uniform flow field before initiation of frozen wall. Orange squares show hypothetical area of contamination; blue squares show location of cooling system (nodes with specified temperature boundary condition). Small green circles denote the observation points defined in figure 5. See appendix 1 for a list of symbols.

Figure 6.
Illustrations showing simulated evolution of ice saturation (S_I) for “containment by frozen groundwater wall” example. Colored contours show S_I values (S_I=1, represents completely frozen groundwater). Velocity vectors that show flow direction and velocity magnitude are represented in a “windsock” style: the direction of flow is indicated by the direction of the line that emanates from the dot, the dot indicates the location of the velocity vector, and the line length is proportional to the magnitude of the velocity. Velocity vectors are plotted at centroid of every second finite element in each direction, and only vectors with the three highest orders of magnitude are shown. Panel A (0 hours) shows uniform flow field before initiation of frozen wall. Orange squares show hypothetical area of contamination; blue squares show location of cooling system (nodes with specified temperature boundary condition). Small green circles denote the observation points defined in figure 5. See appendix 1 for a list of symbols.

Note that the same results as presented here for an aquifer of unit thickness are obtained for an aquifer of any thickness in which the contaminated zone and freezing coils extend from the aquifer top to the aquifer bottom. Also note that this setup may equivalently be considered as the northern half of an area with an east-west axis of symmetry. Further, note that flow field effects of varying groundwater density are neglected in this two-dimensional areal example.

The domain is discretized into 100 elements in the x-direction and 60 elements in the y-direction (∆x=∆y=0.25 m). Time steps begin at ∆t=1 hour, increasing each subsequent time step by a factor of 1.02 to a maximum of ∆t=1 month. Due to the small initial time-step size, iterations for nonlinearity are not used, so there is one iteration per time step. The initial pressure distribution is set as a simple linear function of x, causing uniform flow to occur initially from left to right before freezing begins in the wall. Five SUTRA observation points are positioned along the southern edge of the model domain. The complete model configuration and parameters are provided in figure 5 and table 2. Simulation results are shown in figures 6 and 7.

Table 2.

Parameters of “containment by frozen groundwater wall” simulation.

^{[m, meter; kg/(m s²), kilogram per meter per second squared; J/(kg °C), joule per kilogram per degree
Celsius; J/(s m °C), Joule per second per meter per degree Celsius; kg/m³, kilogram per cubic meter; °C, degree Celsius; kg/(m³ °C), kilogram per cubic meter per degree Celsius; m², square meter; m/s, meter per second; m/s², meter per second squared; s, second; m, meter; h, hour; cm, centimeter; ≥, greater
than or equal to]}

Table 2. Parameters of “containment by frozen groundwater wall” simulation.
Parameter	Value
Aquifer thickness	1 m
Model extent (east-west x, north-south y)	25 m, 15 m
Porosity	0.1
Fluid compressibility	4.47×10^-10 [kg/(m‧s²)]⁻¹
Fluid specific heat	4,182. J/(kg‧°C)
Fluid thermal conductivity	0.6 J/(s‧m‧°C)
Linear fluid density function: [here set to a constant density of 1,000 kg/m³]
Base density (RHOL0)	1,000 kg/m³
Base temperature (URHOL0)	20 °C
Coefficient of density change (DRLDU)	0 kg/(m³‧°C)
Solid matrix compressibility (COMPMA)	1×10^-8 [kg/(m‧s²)]⁻¹
Solid grain density (RHOS)	2,600 kg/m³
Solid grain specific heat	840 J/(kg‧°C)
Solid grain thermal conductivity	3.5 J/(s‧m‧°C)
Ice compressibility (COMPI)	0 [kg/(m‧s²)]⁻¹
Ice density	920 kg/m³
Ice thermal conductivity	2.14 J/(s‧m‧°C)
Ice specific heat	2,108 J/(kg‧°C)
Latent heat of fusion	3.34×10⁵ J/kg
Permeability (max and min)	1×10⁻¹⁰ m² (hydraulic conductivity about 10⁻³ m/s)
Longitudinal dispersivity (in max and min permeability directions)	0.1 m
Transverse dispersivity (in max and min permeability directions)	0.01 m
Gravity vector components (GRAVX, GRAVY)	0, 0 m/s²
Piecewise-linear liquid-water saturation function:
Minimum liquid saturation (SLSATRES)	0.01
Temperature at which minimum liquid saturation occurs (TLRES)	–2 °C
Piecewise-linear relative permeability function:
Liquid saturation at which minimum relative permeability occurs (SLRKMIN)	0.01
Minimum relative permeability (RKMIN)	1×10⁻⁶
Numerical parameters or controls
Finite elements (in x,y directions) -fishnet mesh-	100, 60
Initial time step	1 h
Time step change cycle (NTCYC)	1
Time step multiplier (TCMULT)	1.02
Maximum allowed time step (TCMAX)	730.5 h (= 1 month)
Time steps in pressure and temperature solution cycle (NPCYC, NUCYC)	1, 1
Maximum time steps in simulation (NTMAX)	1,000 (any number≥420)
Maximum absolute simulation time	1×10⁵ h
Output cycles (NPRINT, NCOLPR, LCOLPR)	10, 10, 10
Iterations for nonlinearity (ITRMAX)	1
Solver for pressure, solver for temperature	DIRECT, DIRECT
Observation locations (five selected nodes along southern edge, y=0 m; fig. 5)	x=90 m, 130 m, 170 m, 210 m, 250 m
Initial and boundary conditions
Initial temperature of aquifer	+5 °C
Temperature of cooling system	specified T=−5 °C
Initial pressure of aquifer	p=500−20x [kg/(m‧s²)]
West boundary condition for fluid mass balance	specified p=500 [kg/(m‧s²)] (about 5 cm hydraulic head) —inflowing water has T=+5 °C
East boundary condition for fluid mass balance	specified p=0 [kg/(m‧s²)] (0 cm hydraulic head) —any inflowing water has T=+9,999 °C (arbitrary value as no water ever flows in here)
East and west boundary conditions for energy balance	insulated (energy enters via the western inflow and exits via the eastern outflow)
North and south boundary conditions for fluid mass balance	no flow
North and south boundary conditions for energy balance	insulated

First set of graph lines is mostly grouped on the bottom, and second and third groups
are mostly flat at the top and decline at different slopes. — Figure 7.
Graphs showing simulated histories of pressure, temperature, and liquid saturation (S_L) at observation nodes for “containment by frozen groundwater wall” example. Colors indicate observation point location numbers (see figs. 5 and 6 for locations). See appendix 1 for a list of symbols. Liquid saturation (S_L) is dimensionless. kg/m‧s², kilogram per meter-second-squared; °C, degree Celsius.

The constant west-to-east pressure gradient initially drives a uniform eastward fluid velocity of about 1.3×10⁻⁵ meters per second (m/s; about 1.1 meters per day [m/d]). This would spread the contaminant eastward at a rate greater than 1 m/d (assuming that the contaminant does not undergo physical or chemical retardation). As soon as the groundwater within the freezing coil zone freezes (instantaneously in this simulation, immediately after the state shown in fig. 6A, because the coils are instantaneously at the specified temperature of −5 °C), a barrier is formed to flow, and fluid velocity throughout the contaminated zone decreases by about 10 times, decreasing the rate of contaminant release downstream. Thereafter, while the entire contaminated zone gradually freezes to maximum ice saturation (0.99), groundwater velocity in the growing, mostly frozen, zone decreases to less than one thousandth (1/1,000) of its initial value (fig. 6B–D), essentially immobilizing the contaminant. Finally, fluid in the entire contaminated zone becomes immobilized (fig. 6E). As the frozen wall develops with time (fig. 7A), observed pressure upstream of the wall increases, while observed pressure downstream first decreases due to the ever-increasing size of the frozen barrier blocking the flow field, then increases once the location is fully frozen. Temperatures at all observation points (fig. 7B) decrease with time, undergoing an inflection during ice formation (between 0 °C and −2 °C). Liquid saturations decrease at all observation points (fig. 7C) while ice forms. Although it takes several years for the system to reach a steady state with the coils continuously cooling the aquifer, the dynamics that cause practical immobilization of the contaminant occur within the first few months of operation. Even for this simple example of freezing in a saturated aquifer system, it can be understood that the dynamic evolution of physical and thermal conditions associated with subsurface freezing can be quite complex.

5.2. Unsaturated Groundwater Freeze-Thaw Example—Soil Column Experiment

This example demonstrates unsaturated zone freezing owing to a decrease in surface temperature and subsequent thawing owing to upward flow (water-table rise) from the unfrozen aquifer below the water table. It represents conditions of a hypothetical two-part, 1-year controlled cryohydrogeologic laboratory experiment on a 2-m-high column of soil. Part A of the experiment comprises top-down freezing of the unsaturated zone and part B comprises bottom-up thawing of the part-A column owing to liquid water inflow at the bottom.

The part-A experiment begins with an unsaturated column in steady state, in which the water table is held at the bottom, and there is a natural linear hydrostatic pressure decrease upwards. Total water saturation decreases from fully saturated conditions at the bottom edge (which represents the top of an aquifer’s saturated zone) to nearly irreducible total saturation at the top edge. Initially the entire column is at +2 °C, then at the start of the part-A experiment the top temperature is decreased to –8 °C and held at this value, while the bottom is held at +2 °C (that is, no change at the bottom). This initiates downward freezing of the unsaturated zone in the column. With these fixed top and bottom pressure and temperature settings, conditions in the column are allowed to evolve for 180 days, by which time a new near-steady-state condition has developed.

The initial condition for the part B experiment is the unsaturated, partly frozen column that developed in part A upon reaching 180 days. At the start of part B, liquid water at a temperature +2 °C is injected into the unfrozen bottom of the column at a rate of 1 centimeter (cm) of standing liquid water per day, and this injection rate is held constant for the 180-day duration of part B. The fluid injection at the bottom causes liquid to flow upward, fully saturating the empty pore space as it rises, eventually passing through the entire column and reaching the top. The top of the column is thermally insulated. Before the column becomes fully saturated to the top, there is no discharge at the top. When the column becomes fully saturated, the top pressure is then held at approximately atmospheric pressure and discharge of liquid occurs naturally at the column top. Due to supercooled initial conditions in the column (temperatures less than 0 °C, at which freezing occurs according to the freezing function specified), the upward-flowing water begins to thaw the column due to the temperature of inflowing liquid being greater than 0 °C. Eventually, from the bottom up, the column gradually becomes ice-free, and after 180 days of upward flow the temperature of the column nearly reaches the inflow temperature of +2 °C throughout.

The model configuration and parameters are shown in figure 8 and table 3. Simulation results are shown in figures 9, 10, 11, and 12.

Rectangle with squares at each corner and arrows from bottom left corner going up
along the x and y axes — Figure 8.
Illustration showing geometry of “soil column experiment” example. See appendix 1 for a list of symbols. m, meter.

Table 3.

Parameters of “soil column experiment” simulation.

^{[m, meter; kg/(m s²), kilogram per meter per second squared; J/kg, joule per kilogram; °C, degree Celsius;
J/(s m °C), joule per second per meter per degree Celsius; kg/m³, kilogram per cubic meter; kg/(m³ °C), kilogram per cubic meter per degree Celsius; m², square meter; d, day]}

Table 3. Parameters of “soil column experiment” simulation.
Parameter	Value
Column thickness and width (x-direction, x=0 m at left edge)	1 m
Column height (y-direction, y=0 m at bottom)	2 m
Porosity	0.3
Fluid compressibility	4.47×10⁻¹⁰ [kg/(m‧s²)]⁻¹
Fluid specific heat	4,182 J/(kg‧°C)
Fluid thermal conductivity	0.6 J/(s‧m‧°C)
Linear fluid density function: [here set to a constant density of 1,000 kg/m³]
Base density (RHOL0)	1,000 kg/m³
Base temperature (URHOL0)	20 °C
Coefficient of density change (DRLDU)	0 kg/(m³‧°C)
Solid matrix compressibility (COMPMA)	1×10⁻⁸ [kg/(m‧s²)]⁻¹
Solid grain density (RHOS)	2,600 kg/m³
Solid grain specific heat	840 J/(kg‧°C)
Solid grain thermal conductivity	3.5 J/(s‧m‧°C)
Ice compressibility (COMPI)	0 [kg/(m‧s²)]⁻¹
Ice density	920 kg/m³
Ice thermal conductivity	2.14 J/(s‧m‧°C)
Ice specific heat	2,108. J/(kg‧°C)
Latent heat of fusion	3.34×10⁵ J/kg
Permeability (max and min)	5×10⁻¹¹ m² (hydraulic conductivity about 5×10⁻⁴ m/s)
Longitudinal dispersivity (in max and min permeability directions)	0.001 m
Transverse dispersivity (in max and min permeability directions)	0.001 m
Gravity vector components (GRAVX, GRAVY)	0, –9.81 m/s²
Piecewise-linear total-water saturation function:
Residual total-water saturation due to pressure decrease (SWRES)	0.01
Air-entry pressure (PENT)	0 [kg/(m s²)]
Pressure at residual total-water saturation (PSWRES)	–2×10⁺⁴ [kg/(m‧s²)]
Piecewise-linear relative permeability function:
Liquid saturation at which minimum relative permeability occurs (SLRKRES)	0.01
Minimum relative permeability (RKRES)	0.01
Piecewise-linear liquid-water saturation function:
Minimum liquid saturation due to temperature decrease (freezing) (SLSATRES)	0.01
Temperature at which minimum liquid saturation occurs (TLRES)	–2 °C
Numerical parameters or controls
Finite elements (in horizontal, vertical directions)	1, 100
Time step	Part A: 0.01 d, Part B: 0.001 d
Time steps in pressure and temperature solution cycle (NPCYC, NUCYC)	1, 1
Maximum absolute simulation time	180 d
Maximum time steps in simulation (NTMAX)	Part A: 1,800; Part B: 18,000
Output cycles (NPRINT, NCOLPR, LCOLPR)	100, 100, 100
Solver for pressure, solver for temperature	DIRECT, DIRECT
Observations	at all nodes on left edge
Part A
Initial pressure of column (Hydrostatic, with p=0 at bottom where y=0; initial p values are calculated from the formula)	p=(9.81)(1,000.)(−y) [kg/(m‧s²)]
Initial temperature of column	+2 °C
Top boundary condition for fluid mass balance	No−flow (implemented via ‘drain’ condition, see Part B top boundary condition for fluid mass balance)
Top boundary condition for energy balance	Specified T=−8 °C
Bottom boundary condition for fluid mass balance	Specified p=0 [kg/(m s²)] and any inflowing water has T=+2 °C
Bottom boundary condition for energy balance	Specified T=+2 °C
Iterations for Nonlinearity (ITRMAX)	1
Part B
Initial pressure of column	Pressures taken from the end of part A
Initial temperature of column	Temperatures taken from the end of part A
Top boundary condition for fluid mass balance:	0 [kg/(m‧s²)], 0 kg/s
Generalized-flow boundary condition, a drain [see Provost and Voss (2019, section 1.4.5, “River and Drain”)]
p at first point (PBG1), rate of inflow at first point (QPBG1)	196.2 [kg/(m‧s²)]
p at second point (PBG2), rate of inflow at second point (QPBG2)	–4.96×10⁻⁴ kg/s
Limit on flow at point 1,2 (CPQL1, CPQL2)	Q (means that max limit is 0 [kg/(m‧s²)]), N (means no limit)
Temperature of any water entering top (arbitrary value as this never occurs) (temperature of any water exiting top is at temperature of ambient water at column top, changing with time)	9,999.9999 °C
Top boundary condition for energy balance:	Insulated for conduction
Bottom boundary condition for fluid mass balance:	Specified total inflow QIN=4.96032×10⁻⁴ kg/s (equivalent to 1 centimeter per day liquid water flowing into base of column
Bottom boundary condition for energy balance:	Inflowing water has T=+2 °C Specified T=+2 °C
Max iterations for nonlinearity (ITRMAX)	10
Convergence criterion for p (RPMAX)	275 [kg/(m‧s²)]
Convergence criterion for T (RUMAX) (set very large, so this does not control iterations)	999999 °C

Boxes changing colors based on level of ice or liquid saturation — Figure 9.
Illustrations of part A of experiment: freezing of unsaturated zone downward from the column top, showing ice and liquid-water saturations at indicated time steps (in days). The total-water saturation remains constant throughout part A, with values equivalent to the liquid saturation panel for 0 days. Minimum total and liquid saturation at 0 days is 0.02881 at top of column. Maximum ice saturation is 0.5988 at 180 days. (Ice saturation of zero and liquid saturation of 0.01 are shown as white.) See appendix 1 for a list of symbols.

Two of the graphs only have one line, the rest have five lines — Figure 11.
Graphs showing simulated evolution of vertical profiles along A, column of pressure, B, temperature, C, total water saturation, D, liquid saturation, and E, ice saturation for part A of the “soil-column experiment” example. Part A is simulated for 180 days. See appendix 1 for a list of symbols. kg/m‧s², kilogram per meter per second squared; °C, degree Celsius.

Two of the graphs are mostly straight lines, the rest have curves — Figure 12.
Graphs showing simulated evolution of vertical profiles along column of A, pressure, B, temperature, and C, total water saturation, D, liquid saturation, and E, ice saturation for part B of the “soil-column experiment” example. Part B is simulated for 180 days. Values for part B at 0 days are those that result from part A at 180 days. See appendix 1 for a list of symbols. m, meter; kg, kilogram; s, second; °C, degree Celsius; d, day; >, greater than.

Part A.—Following the instantaneous decrease in top temperature, simulated temperatures in the column decrease and a diffuse freezing front propagates downward through the variably saturated column, increasing ice saturations and decreasing liquid saturations (figs. 9 and 11). Nearly steady conditions in temperature and saturations develop in the soil column by about 90 days. Pressures naturally remain constant throughout the simulation, due to gravity-produced hydrostatic conditions, based on the specified pressure of 0 at the bottom. Pressures control only total-water saturation via the total-water saturation function (table 3) and do not affect apportioning of total water between liquid and ice phases. Total saturation is unchanging in part A, resulting from the steady pressure distribution, and it is equivalent to the liquid saturation at 0 days (fig. 11A). Thus, at any given location in the column, the decrease in liquid saturation is equal to the increase in ice saturation. The simulated liquid and ice saturations are controlled by temperature changes affecting the freezing function (table 3). Because column temperature increases downward and available liquid saturation decreases upward, the zone of highest ice saturation values occurs just below mid-column (figs. 9A and 11E).

Part B.—Flow of +2 °C liquid water entering the bottom of the column leads to a more complex evolution of total water, liquid, and ice saturations in part B than in part A of the experiment (figs. 10 and 12). The column resaturates from the bottom up and the newly arrived liquid begins to thaw the ice (figs. 10 and 12D–E). At about 30 days, the pressure at the top of the column exceeds zero and the boundary condition there changes from one of no flow to a generalized pressure boundary (both of which are part of the prescribed drain condition; see Provost and Voss [2019] for information on the drain condition). Thereafter the simulated top pressure remains just above 0 kg/(m‧s²). At this time (about 30 days), the column has become fully saturated with water, even though ice remains in the column until nearly 180 days.

Note that iterations for nonlinearity are required to maintain numerical stability early in the simulation when the forcing of 1 centimeter per day of liquid water into a frozen column causes numerical sensitivity of the ice-liquid transition [notice the wobbles in the ice saturation profile (fig. 12E) at an elapsed time of 3 days]. This sensitivity is also the reason that a smaller time step size is required for part B than for part A of the experiment. A finer mesh (smaller finite elements) in the vertical direction than employed here might also aid the stability of the early period of the numerical solution for part B. Another purpose of nonlinear iterations is to resolve the switchover from a no-flow condition to the generalized pressure boundary. Without such iterations, beginning at the switchover, the numerical solution for pressure at the top nodes would oscillate for the rest of the simulation with a repeating cycle of a few time steps, with the oscillations obscuring the numerical results. Nonlinearity iterations may occur only for the time step on which the switchover occurs, but to allow this to happen, such iterations must be allowed for the entire simulation (table 3). A third purpose of nonlinear iterations is to resolve the significant changes in the column that occur as soon as liquid water is injected into the column base. The convergence tolerance value for pressure, RPMAX, was selected by experimentation in a series of simulations, so that it would provide a reasonable numerical solution for these two sensitive parts of the simulated dynamics, without requiring many iterations per time step. Use of nonlinear iteration does not cause iterations on all time steps—only on those for which it is needed. The convergence tolerance for temperature, RUMAX, is set extremely high so that it is always satisfied and will never cause additional iterations.

For the entire simulated part B column experiment, thawing occurs from below, due to the warming effect of the upward-flowing liquid. This thawing occurs very slowly, taking more than 3 months, due to the relatively large latent-heat energy required to melt the ice. The energy input to the thawing process is provided by the +2 °C liquid water entering the bottom of the column. Due to the large latent-heat energy required to melt the ice, the long period required to thaw the column is a consequence of the relatively low amount of energy provided by the inflow. Almost all the ice is finally melted just after 90 days, at which time the temperature at the top of the column finally reaches 0 °C (fig. 10A [at 90 days] and 12B and E). Thereafter, the ice-free column warms to nearly +2 °C by 180 days. Nonlinear and complex interactions between groundwater flow and freeze-thaw processes occur, especially for unsaturated conditions, even for this simple one-dimensional example.

References Cited

Andersland, O.B., and Ladanyi, B., 1994, An introduction to frozen ground engineering: New York, Chapman and Hall, 352 p., accessed March 23, 2022, at https://doi.org/10.1007/978-1-4757-2290-1.

Anderson, D.M., and Morgenstern, N.R., 1973, Physics, chemistry and mechanics of frozen ground—A review, in North American Contribution, Permafrost, Second International Conference, Yakutsk, U.S.S.R., July 13–28, 1973: National Research Council of Canada and U.S. National Academy of Sciences, p. 257–288.

Anderson, D.M., and Tice, A.R., 1972, Predicting unfrozen water contents in frozen soils from surface area measurements: Highway Research Record, v. 393, p. 12–18. [Also available at https://trid.trb.org/view/126420.]

Azmatch, T.F., Sego, D.C., Arenson, L.U., and Biggar, K.W., 2012, Using soil freezing characteristic curve to estimate the hydraulic conductivity function of partially frozen soils: Cold Regions Science and Technology, v. 83–84, p. 103–109, accessed March 23, 2022, at https://doi.org/10.1016/j.coldregions.2012.07.002.

Briggs, M.A., Walvoord, M.A., McKenzie, J.M., Voss, C.I., Day-Lewis, F.D., and Lane, J.W., 2014, New permafrost is forming around shrinking Arctic lakes, but will it last?: Geophysical Research Letters, v. 41, no. 5, p. 1585–1592, accessed March 23, 2022, at https://doi.org/10.1002/2014GL059251.

Brooks, R.H., and Corey, A.T., 1964, Hydraulic properties of porous media: Colorado State University Hydrology Paper no. 3, 27 p. [Also available at https://mountainscholar.org/bitstream/handle/10217/61288/HydrologyPapers_n3.pdf.]

Chen, L., Fortier, D., McKenzie, J.M., Voss, C.I., and Lamontagne-Hallé, P., 2023, Subsurface porewater flow accelerates talik development under the Alaska highway, Yukon—A prelude to road collapse and final permafrost thaw?: Water Resources Research, v. 59, no. 4, paper e2022WR032578, 21 p., accessed June 29, 2023, at https://doi.org/10.1029/2022WR032578.

Chen, L., VossC.I., Fortier, D., and McKenzie, J.M., 2021, Surface energy balance of sub-Arctic roads with varying snow regimes and properties in permafrost regions: Permafrost and Periglacial Processes, v. 32, no. 4, p. 681-701, accessed June 29, 2023, at https://doi.org/10.1002/ppp.2129.

Evans, S.G., Ge, S., Voss, C.I., and Molotch, N.P., 2018, The role of frozen soil in groundwater discharge predictions for warming alpine watersheds: Water Resources Research, v. 54, no. 3, p. 1599–1615, accessed March 23, 2022, at https://doi.org/10.1002/2017WR022098.

Evans, S.G., Godsey, S.E., Rushlow, C.R., and Voss, C., 2020, Water tracks enhance water flow above permafrost in upland Arctic Alaska hillslopes: Journal of Geophysical Research: Earth Surface, v. 125, no. 2, paper e2019JF005256, 18 p., accessed March 23, 2022, at https://doi.org/10.1029/2019JF005256.

Ge, S., McKenzie, J.M., Voss, C.I., and Wu, Q., 2011, Exchange of groundwater and surface-water mediated by permafrost response to seasonal and long term air temperature variation: Geophysical Research Letters, v. 38, no. 14, L14402, 6 p., accessed March 23, 2022, at https://doi.org/10.1029/2011GL047911.

Grenier, C., Anbergen, H., Bense, V., Chanzy, Q., Coon, E., Collier, N., Costard, F., Ferry, M., Frampton, A., Frederick, J., Gonçalvès, J., Holmén, J., Jost, A., Kokh, S., Kurylyk, B., McKenzie, J., Molson, J., Mouche, E., Orgogozo, L., Pannetier, R., Rivière, A., Roux, N., Rühaak, W., Scheidegger, J., Selroos, J.-O., Therrien, R., Vidstrand, P., and Voss, C., 2018, Groundwater flow and heat transport for systems undergoing freeze-thaw—Intercomparison of numerical simulators for 2D test cases: Advances in Water Resources, v. 114, p. 196–218, accessed March 23, 2022, at https://doi.org/10.1016/j.advwatres.2018.02.001.

Holten, V., and Anisimov, M.A., 2012, Entropy-driven liquid–liquid separation in supercooled water: Scientific Reports, v. 2, article 713, 7 p., accessed March 23, 2022, at https://doi.org/10.1038/srep00713.

Jame, Y.-W., 1977, Heat and mass transfer in freezing soils: University of Saskatchewan dissertation, 212 p. [Also available at http://hdl.handle.net/10388/etd-11282008-091133.]

Jame, Y.-W., and Norum, D.I., 1980, Heat and mass-transfer in a freezing unsaturated porous medium: Water Resources Research, v. 16, no. 4, p. 811–819, accessed March 23, 2022, at https://doi.org/10.1029/WR016i004p00811.

Kell, G.S., 1967, Precise representation of volume properties of water at one atmosphere: Journal of Chemical & Engineering Data, v. 12, no. 1, p. 66–69, accessed March 23, 2022, at https://doi.org/10.1021/je60032a018.

Koopmans, R.W.R., and Miller, R.D., 1966, Soil freezing and soil water characteristics curves: Soil Science Society of America Journal, v. 30, no. 6, p. 680–685, accessed March 23, 2022, at https://doi.org/10.2136/sssaj1966.03615995003000060011x.

Kozlowski, T., 2007, A semi-empirical model for phase composition of water in clay-water systems: Cold Regions Science and Technology, v. 49, no. 3, p. 226–236, accessed March 23, 2022, at https://doi.org/10.1016/j.coldregions.2007.03.013.

Kurylyk, B.L., and Watanabe, K., 2013, The mathematical representation of freezing and thawing processes in variably-saturated, non-deformable soils: Advances in Water Resources, v. 60, p. 160–177, accessed March 23, 2022, at https://doi.org/10.1016/j.advwatres.2013.07.016.

Kurylyk, B.L., Hayashi, M., Quinton, W.L., McKenzie, J.M., and Voss, C.I., 2016, Influence of vertical and lateral heat transfer on permafrost thaw, peatland landscape transition, and groundwater flow: Water Resources Research, v. 52, no. 2, p. 1286–1305, accessed March 23, 2022, at https://doi.org/10.1002/2015WR018057.

Kurylyk, B.L., MacQuarrie, K.T.B., and Voss, C.I., 2014a, Climate change impacts on the temperature and magnitude of groundwater discharge from shallow, unconfined aquifers: Water Resources Research, v. 50, no. 4, p. 3253–3274, accessed March 23, 2022, at https://doi.org/10.1002/2013WR014588.

Kurylyk, B.L., McKenzie, J.M., MacQuarrie, K.T.B., and Voss, C.I., 2014b, Analytical solutions for benchmarking cold regions subsurface water flow and energy transport models—One-dimensional soil thaw with conduction and advection: Advances in Water Resources, v. 70, p. 172–184, accessed March 23, 2022, at https://doi.org/10.1016/j.advwatres.2014.05.005.

Lamontagne-Hallé, P., McKenzie, J.M., Kurylyk, B.L., and Zipper, S.C., 2018, Changing groundwater discharge dynamics in permafrost regions: Environmental Research Letters, v. 13, no. 8, article 084017, 12 p., accessed March 23, 2022, at https://doi.org/10.1088/1748-9326/aad404.

Lunardini, V.J., 1981, Heat transfer in cold climates: New York, Van Nostrand Reinhold Company, 731 p.

Lunardini, V.J., 1988, Freezing of soil with an unfrozen water content and variable thermal properties: U.S. Army Corps of Engineers, Cold Regions Research and Engineering Laboratory Report 88–2, 23 p. [Also available at http://hdl.handle.net/11681/9076.]

Lundin, L.-C., 1990, Hydraulic properties in an operational model of frozen soil: Journal of Hydrology, v. 118, no. 1–4, p. 289–310, accessed March 23, 2022, at https://doi.org/10.1016/0022-1694(90)90264-X.

McKenzie, J.M., and Voss, C.I., 2013, Permafrost thaw in a nested groundwater-flow system: Hydrogeology Journal, v. 21, no. 1, p. 299–316, accessed March 23, 2022. 10.1007/s10040-012-0942-3.

McKenzie, J.M., Siegel, D.I., Rosenberry, D.O., Glaser, P.H., and Voss, C.I., 2006, Heat transport in the Red Lake Bog, Glacial Lake Agassiz Peatlands: Hydrological Processes, v. 21, no. 3, p. 369–378, accessed March 23, 2022, at https://doi.org/10.1002/hyp.6239.

McKenzie, J.M., Voss, C.I., and Siegel, D.I., 2007, Groundwater flow with energy transport and water-ice phase change—Simulations, benchmarks, and application to freezing in peat bogs: Advances in Water Resources, v. 30, no. 4, p. 966–983, accessed March 23, 2022, at https://doi.org/10.1016/j.advwatres.2006.08.008.

Miller, O., Voss, C.I., Solomon, D.P, Miège, C., Forster, R., Schmerr, N., and Montgomery, L., 2022, Hydrologic modeling of a perennial firn aquifer in southeast Greenland: Journal of Glaciology, v. 69, no. 275, p. 607-622, accessed June 29, 2023, at https://doi.org/10.1017/jog.2022.88.

Mottaghy, D., and Rath, V., 2006, Latent heat effects in subsurface heat transport modelling and their impact on paleotemperature reconstructions: Geophysical Journal International, v. 164, no. 1, p. 236–245, accessed March 23, 2022, at https://doi.org/10.1111/j.1365-246X.2005.02843.x.

Newman, G., and Wilson, G., 1997, Heat and mass transfer in unsaturated soils during freezing: Canadian Geotechnical Journal, v. 34, no. 1, p. 63–70, accessed March 23, 2022, at https://doi.org/10.1139/t96-085.

Nield, D.A., and Bejan, A., 1999, Convection in porous media (2d ed.): New York, Springer-Verlag, 546 p. [Also available at https://doi.org/10.1007/978-1-4757-3033-3.]

Niu, F., Cheng, G., Niu, Y., Zhang, M., Luo, J., and Lin, Z., 2016, A naturally-occurring ‘cold earth’ spot in northern China: Scientific Reports, v. 6, article 34184, 7 p., accessed March 23, 2022, at https://doi.org/10.1038/srep34184.

Painter, S.L., 2011, Three-phase numerical model of water migration in partially frozen geological media—Model formulation, validation, and applications: Computational Geosciences, v. 15, no. 1, p. 69–85, accessed March 23, 2022, at https://doi.org/10.1007/s10596-010-9197-z.

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., accessed March 23, 2022, at https://doi.org/10.3133/tm6A52.

Rempel, A.W., 2012, Hydromechanical processes in freezing soils: Vadose Zone Journal, v. 11, no. 4, article vzj2012.0045, 10 p., accessed March 23, 2022, at https://doi.org/10.2136/vzj2012.0045.

Rey, D.M., Walvoord, M.A., Minsley, B.J., Ebel, B.A., Voss, C.I., and Singha, K., 2020, Wildfire‐initiated talik development exceeds current thaw projections—Observations and models from Alaska’s continuous permafrost zone: Geophysical Research Letters, v. 47, no. 15, paper 087565, 11 p., accessed March 23, 2022, at https://doi.org/10.1029/2020GL087565.

Rühaak, W.H., Anbergen, C., Grenier, J., McKenzie, J.M., Kurylyk, B., Molson, J., Roux, N., and Sass, I., 2015, Benchmarking numerical freeze/thaw models: Energy Procedia, v. 76, p. 301–310, accessed March 23, 2022, at https://doi.org/10.1016/j.egypro.2015.07.866.

Rushlow, C.R., Sawyer, A.H., Voss, C.I., and Godsey, S.E., 2020, The influence of snow cover, air temperature, and groundwater flow on the active-layer thermal regime of Arctic hillslopes drained by water tracks: Hydrogeology Journal, v. 28, p. 2057–2069, accessed March 23, 2022, at https://doi.org/10.1007/s10040-020-02166-2.

Spaans, E.J.A., and Baker, J.M., 1996, The soil freezing characteristic—Its measurement and similarity to the soil moisture characteristic: Soil Science Society of America Journal, v. 60, no. 1, p. 13–19, accessed March 23, 2022, at https://doi.org/10.2136/sssaj1996.03615995006000010005x.

Taylor, G.S., and Luthin, J.N., 1978, A model for coupled heat and moisture transfer during soil freezing: Canadian Geotechnical Journal, v. 15, no. 4, p. 548–555, accessed March 23, 2022, at https://doi.org/10.1139/t78-058.

Tilton, L.W., and Taylor, J.K., 1937, Accurate representation of the refractivity and density of distilled water as a function of temperature: Journal of Research of the National Bureau of Standards, v. 18, paper RP971, p. 205–214. [Also available at https://nvlpubs.nist.gov/nistpubs/jres/18/jresv18n2p205_A1b.pdf.]

U.S. Geological Survey, 2022, SUTRA—A model for 2D or 3D saturated-unsaturated, variable-density ground-water flow with solute or energy transport: U.S. Geological Survey software, accessed July 22, 2022, at https://doi.org/10.5066/P9PPEHHM.

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, accessed March 23, 2022, at https://doi.org/10.2136/sssaj1980.03615995004400050002x.

Voss, C.I., and Provost, A.M., 2002, SUTRA—A model for 2D or 3D saturated-unsaturated, variable-density ground-water flow with solute or energy transport [version 2D3D.1, September 22, 2010]: U.S. Geological Survey Water-Resources Investigations Report 02–4231, 291 p., accessed March 23, 2022, at https://doi.org/10.3133/wri024231.

Voss, C.I., and Provost, A.M., 2024, SUTRA—A model for 2D or 3D saturated-unsaturated, variable-density groundwater flow with solute or energy transport, version 4.0: U.S. Geological Survey software release, https://doi.org/10.5066/P9OL5IYX.

Walvoord, M.A., and Kurylyk, B.L., 2016, Hydrologic impacts of thawing permafrost—A review: Vadose Zone Journal, v. 15, no. 6, article vzj2016.01.0010, 20 p., accessed March 23, 2022, at https://doi.org/10.2136/vzj2016.01.0010.

Walvoord, M.A., Voss, C.I., Ebel, B., and Minsley, B.J., 2019, Development of perennial thaw zones in boreal hillslopes enhances potential mobilization of permafrost carbon: Environmental Research Letters, v. 14, no. 1, article 015003, 11 p., accessed March 23, 2022, at https://doi.org/10.1088/1748-9326/aaf0cc.

Watanabe, K., and Wake, T., 2009, Measurement of unfrozen water content and relative permittivity of frozen unsaturated soil using NMR and TDR: Cold Regions Science and Technology, v. 59, no. 1, p. 34–41, accessed March 23, 2022, at https://doi.org/10.1016/j.coldregions.2009.05.011.

Wellman, T.P., Voss, C.I., and Walvoord, M.A., 2013, Impacts of climate, lake size, and supra- and sub-permafrost groundwater flow on lake-talik evolution, Yukon Flats, Alaska (USA): Hydrogeology Journal, v. 21, no. 1, p. 281–298, accessed March 23, 2022, at https://doi.org/10.1007/s10040-012-0941-4.

Wicky, J., and Hauck, C., 2017, Numerical modelling of convective heat transport by air flow in permafrost talus slopes: The Cryosphere, v. 11, p. 1311–1325, accessed March 23, 2022, at https://doi.org/10.5194/tc-11-1311-2017.

Williams, P.J., 1967, Properties and behavior of freezing soils: Norwegian Geotechnical Institute Publication 72, 119 p.

Zhao, Y., Nishimura, T., Hill, R., and Miyazaki, T., 2013, Determining hydraulic conductivity for air-filled porosity in an unsaturated frozen soil by the multistep outflow method: Vadose Zone Journal, v. 12, no. 1, article vzj2012.0061, 10 p., accessed March 23, 2022, at https://doi.org/10.2136/vzj2012.0061.

Zhou, X., Zhou, J., Kinzelbach, W., and Stauffer, F., 2014, Simultaneous measurement of unfrozen water content and ice content in frozen soil using gamma ray attenuation and TDR: Water Resources Research, v. 50, no. 12, p. 9630–9655, accessed March 23, 2022, at https://doi.org/10.1002/2014WR015640.

Appendix 1. List of Units, Symbols, and Abbreviations

If a symbol appears in a numbered equation, the description of the symbol includes in parentheses the number of the equation in which the symbol first appears. If the symbol does not appear in a numbered equation, or if the first symbol definition in the text is separated significantly from the place where the symbol is first used, the location of the first appearance and first definition is noted.

Table 1.1.

Abbreviations and symbols used in SUTRA model version 4.0.

Table 1.1. Abbreviations and symbols used in SUTRA model version 4.0.
Abbreviation or symbol	Unit	Definition
Generic units
	[E]	energy, equivalent to [M∙L²/s²]
	[L]	length
	[L³]	volume (general)
	[L_G³]	solid mineral grain or solid matrix volume
	[L_I³]	ice volume
	[L_L³]	liquid-water (fluid) volume
	[M]	mass (general)
	[M_G]	solid mineral grain or solid matrix mass
	[M_I]	ice mass
	[M_L]	liquid-water (fluid) mass
	[M_S]	solute or sorbate mass
Units
	[1]	dimensionless
	[°C]	degree Celsius (unit of temperature)
	[cm]	centimeter (unit of [L])
	[d]	day (unit of time)
	[gr]	gram (unit of [M])
	[h]	hour (unit of time)
	[J]	joule (energy unit [kg∙m²/s²])
	[kg]	kilogram (unit of [M])
	[m]	meter (unit of [L])
	[s]	second (unit of time)
	[yr]	year (unit of time)
	[°]	degree (measure of an angle)
Greek lowercase
α_s	[Μ/(L‧s²)]^–1	Porous-matrix compressibility (eq. 7)
α_φ	[(percent)∙ ${(° C)}^{- β_{POW}}$ ]	Power-law model parameter (eq. 59)
α_Lmax	[L]	Scaled longitudinal dispersivity value for the maximum permeability direction (defined in dataset 15; description in app. 2)
α_Lmid	[L]	Scaled longitudinal dispersivity value for the middle permeability direction (defined in dataset 15; description in app. 2)
α_Lmin	[L]	Scaled longitudinal dispersivity value for the minimum permeability direction (defined in dataset 15; description in app. 2)
α_POW(x, y[, z])	[ ${(° C)}^{- β_{POW}}$ ]	Modified power-law model parameter (eq. 57)
α_Tmax	[L]	Scaled transverse dispersivity value for the maximum permeability direction (defined in dataset 15; description in app. 2)
α_Tmid	[L]	Scaled transverse dispersivity value for the middle permeability direction (defined in dataset 15; description in app. 2)
α_Tmin	[L]	Scaled transverse dispersivity value for the minimum permeability direction (defined in dataset 15; description in app. 2)
β_I	[Μ/(L∙s²)]^–1	Ice compressibility (eq. 14)
β_L	[Μ/(L∙s²)]^–1	Liquid-water compressibility (eq. 8)
β_POW(x, y[, z])	[1]	Modified power-law model parameter (eq. 57)
$γ_{0}^{I} (x, y [, z], t)$	[E/(M_I∙s)]	Zero-order energy source in ice (eq. 24)
$γ_{0}^{L} (x, y [, z], t)$	[E/(M_L∙s)] or [(M_S/M_L)/s)]	Zero-order energy or solute-species mass source in liquid water (for energy, eq. 24; for solute, eq. 82, definition before eq. 82)
$γ_{0}^{S} (x, y [, z], t)$	[E/(M_G∙s)] or [(M_S/M_G)/s)]	Zero-order energy or solute-species mass source in solid mineral grains (for energy, eq. 24; for solute, eq. 82, definition before eq. 82)
$γ_{1}^{L} (x, y [, z], t)$	[(M_S/M_L)/s)]	First-order solute-species mass production in liquid water (eq. 82, definition before eq. 82)
$γ_{1}^{S} (x, y [, z], t)$	[(M_S/M_G)/s)]	First-order solute-species mass production in solid mineral grains (eq. 82, definition before eq. 82)
ε(x, y[, z])	[1]	Porosity (eq. 3)
θ₁	[°]	Scaled first permeability angle value (defined in dataset 15; description in app. 2)
θ₂	[°]	Scaled second permeability angle value (defined in dataset 15; description in app. 2)
θ₃	[°]	Scaled third permeability angle value (defined in dataset 15; description in app. 2)
λ(x, y[, z], t)	[E/(s∙L∙°C)]	Bulk thermal conductivity (eq. 31)
λ_A(x, y[, z])	[E/(s∙L∙°C)]	Effective thermal conductivity of air (eq. 44, definition before eq. 44)
λ_BC(x, y[, z])	[1]	Pore-size distribution index for Brooks-Corey total-water saturation function (eq. 49)
λ_I	[E/(s∙L∙°C)]	Thermal conductivity of ice (eq. 42A)
λ_L	[E/(s∙L∙°C)]	Thermal conductivity of liquid water (eq. 42A)
λ_LIS(x, y[, z], t)	[E/(s∙L∙°C)]	Effective thermal conductivity of mixture of liquid water, ice, and solid mineral grains (eq. 42A)
λ_S	[E/(s∙L∙°C)]	Thermal conductivity of a solid mineral grain (eq. 42A)
µ_L(x, y[, z], t)	[M/(L∙s)]	Liquid water viscosity (eq. 19)
ρ_I(x, y[, z], t)	[M_I/L_I³]	Ice density (eq. 2)
ρ_L(x, y[, z], t)	[M_L/L_L³]	Liquid-water density (eq. 2)
ρ_S(x, y[, z], t)	[M_S/L_S³]	Density of solid mineral grain (eq. 24)
ρ_W(x, y[, z], t)	[M/L³]	Mean total-water density (eq. 2)
σ_A	[E/(s∙L∙°C)]	Scaled effective air thermal conductivity value, in the context of input to SUTRA (for notational consistency with σ_L, which is used for both energy transport and solute transport); identically equal to λ_A (defined in dataset 15; description in app. 2)
σ_I	[E/(s∙L∙°C)]	Thermal conductivity of ice, in the context of input to SUTRA (for notational consistency with σ_L, which is used for both energy transport and solute transport); identically equal to λ_I (defined in dataset 9; description in app. 2)
σ_L	[E/(s∙L∙°C)] or [L²/s]	For energy transport, represents liquid-water thermal conductivity, λ_L; for solute transport, represents molecular diffusivity of solute in pure liquid water, D_m (defined in dataset 9; description in app. 2)
σ_S	[E/(s∙L∙°C)]	Scaled solid mineral grain thermal conductivity value, in the context of input to SUTRA (for notational consistency with σ_L, which is used for both energy transport and solute transport); identically equal to λ_S for energy transport (defined in dataset 15; description in app. 2)
$ϕ_{L}^{sat}$	[percent]	Percent moisture content by dry weight under fully saturated conditions (eq. 59)
χ₁	[L_L³/M_G]	Linear, Freundlich or Langmuir adsorption distribution coefficient (defined in dataset 11A; description in app. 2)
χ₂	[1]	Freundlich or Langmuir adsorption coefficient (defined in dataset 11A; description in app. 2)
Greek uppercase
Ω	[1]	Impedance factor (eq. 68)
Roman lowercase
α_VG(x, y[, z])	[(L∙s²)/M]	Van Genuchten model parameter (eq. 47)
c_I(x, y[, z], t)	[E/(M_I∙°C)]	Specific heat of ice (eq. 26B)
c_L(x, y[, z], t)	[E/(M_L∙°C)]	Specific heat of liquid water (eq. 26A)
c_S(x, y[, z], t)	[E/(M_S∙°C)]	Specific heat of solid mineral grains (eq. 26C)
e_I(x, y[, z], t)	[E/M_I]	Energy per unit mass of ice (eq. 24)
$e_{I}^{ff}$	[E/M_I]	Thermal energy per unit mass of ice at free-water freezing temperature (eq. 26B)
e_L(x, y[, z], t)	[E/M_L]	Energy per unit mass of liquid water (eq. 24)
$e_{L}^{*} (x, y [, z], t)$	[E/M_L]	Energy per unit mass of fluid source (eq. 24)
$e_{L}^{ff}$	[E/M_L]	Thermal energy per unit mass of liquid water at free-water freezing temperature (eq. 26A)
e_S(x, y[, z], t)	[E/M_S]	Energy per unit mass of solid mineral grains (eq. 24)
$e_{S}^{ff}$	[E/M_S]	Thermal energy per unit mass of solid mineral grains at free-water freezing temperature (eq. 26C)
$\underline{g}$	[L/s²]	Gravitational acceleration (gravity vector) (eq. 19)
k	[L²]	Permeability under fully saturated, unfrozen conditions (eq. 69)
$\underline{\underline{k}} (x, y [, z], t)$	[L²]	Solid matrix permeability tensor (eq. 19)
k_eff	[L²]	Effective permeability that takes into account both partial saturation and freezing (eq. 68)
$k_{eff}^{uf}$	[L²]	Effective permeability that takes into account partial saturation under unfrozen conditions (eq. 68)
k_max	[L²]	Scaled maximum permeability value (defined in dataset 15; description in app. 2)
k_mid	[L²]	Scaled middle permeability value (defined in dataset 15; description in app. 2)
k_min	[L²]	Scaled minimum permeability value (defined in dataset 15; description in app. 2)
k_r(x, y[, z], t)	[1]	Relative permeability for fluid flow (eq. 19)
$k_{r}^{uf}$	[1]	Relative permeability under unfrozen conditions (eq. 70)
k_rmin(x, y[, z])	[1]	Minimum value of relative permeability (eq. 41, definition before eq. 41)
n_VG(x, y[, z], t)	[1]	Van Genuchten model parameter (eq. 47)
p(x, y[, z], t)	[M/(L∙s²)]	Fluid pressure (eq. 5)
p_ent(x, y[, z], t)	[Μ/(L∙s²)]	Air-entry pressure, first appears (eq. 49, definition before eq. 47)
p_Wres(x, y[, z])	[Μ/(L∙s²)]	Residual-saturation pressure for the piecewise-linear total-water saturation function (eq. 51)
$\underline{q}$	[L/s] or [L_L³/(L²∙s)]	Darcy velocity or volumetric fluid flux (eq. 86, definition before eq. 86)
${\underline{q}}_{e} (x, y [, z], t)$	[E/(L²∙s)]	Conductive-dispersive heat flux (eq. 24)
${\underline{q}}_{S} (x, y [, z], t)$	[M_S/(L²∙s)]	Diffusive-dispersive solute-mass flux (eq. 82, definition before eq. 82)
t	[s]	Time (eq. 3)
$\underline{v} (x, y [, z], t)$	[L/s]	Average fluid velocity (eq. 3)
w_EXP(x, y[, z])	[°C]	Exponential liquid-water saturation function parameter (eq. 55)
Roman uppercase
C(x, y[, z], t)	[M_S/M_L]	Solute-species concentration (mass fraction) (eq. 82, definition before eq. 82)
C_S(x, y[, z], t)	[M_S/M_G]	Concentration of species as sorbate on solid mineral grains expressed as a mass fraction (eq. 82, definition before eq. 82)
C*(x, y[, z], t)	[M_S M_L]	Solute concentration of fluid sources (mass fraction; eq. 82, definition before eq. 82)
$\underline{\underline{D}} (x, y [, z], t)$	[L²/s]	Mechanical dispersion tensor (eq. 31)
D_m	[L²/s]	Molecular diffusivity of solute in pure liquid water (defined in table 2 and dataset 9; description in app. 2)
∆H(x, y[, z])	[E/ M_L]	Latent heat of fusion of free water (eq. 27)
Q_p(x, y[, z], t)	[M_L/(L³∙s)]	Fluid-mass source (eq. 3)
S_A(x, y[, z], t)	[1]	Air saturation (eq. 44, definition before eq. 44)
S_I(x, y[, z], t)	[1]	Ice saturation (eq. 1)
S_L(x, y[, z], t)	[1]	Liquid-water saturation (eq. 1)
$S_{L}^{sat}$	[1]	Liquid-water saturation as a result of the freeze-thaw process in a fully saturated porous medium (eq. 36)
$S_{L}^{*}$	[1]	Drainable liquid-water saturation (defined in eq. 63)
$S_{L}^{†}$	[1]	Relative, scaled saturation for the piecewise-linear relative-permeability function (defined in eq. 67B)
$S_{L_{krmin}}$	[1]	Value of liquid-water saturation at which the relative permeability has its minimum value (k_r=k_rmin; eq. 67B, definition before eq. 41)
S_Lres	[1]	Minimum liquid-water saturation that can occur, taking into account both pressure desaturation and freezing (eqs. 55 and 63, definition after eq. 40)
$S_{Lres}^{sat} (x, y [, z])$	[1]	Minimum liquid-water saturation that can occur as a result of freezing the porous medium (eq. 55, defined after eq. 37)
$S_{op}^{eff} (x, y [, z], t)$	[M/(L∙s²)]^–1	Effective specific pressure storativity (eq. 19)
$S_{op}^{I}$	[Μ/(L∙s²)]^–1	Quantity mathematically analogous to the specific pressure storativity for liquid, but computed using the compressibility of ice (eq. 13)
$S_{op}^{L}$	[Μ/(L∙s²)]^–1	Specific pressure storativity for liquid (eq. 6)
S_W(x, y[, z], t)	[1]	Total-water saturation (eq. 1)
S_Wres(x, y[, z], t)	[1]	Residual total-water saturation (eq. 47, definition after eq. 38)
S_y	[1]	Specific yield (eq. 53)
T(x, y[, z], t)	[°C]	Temperature (eq. 5)
T*(x, y[, z], t)	[°C]	Temperature of fluid source (eq. 30)
T_∆	[°C]	Relative temperature (eq. 54)
T_f(x, y[, z], t)	[°C]	Maximum freezing temperature of pore water (eq. 54, definition before eq. 27)
T_ff	[°C]	Freezing temperature of pure, free water (eq. 26A)
T_Lres(x, y[, z], t)	[°C]	Residual-saturation temperature for the piecewise-linear liquid-water saturation function (eq. 61, definition after eq. 37)

Appendix 2. New and Modified Input Datasets

This appendix lists the U.S. Geological Survey’s Saturated-Unsaturated Transport (SUTRA) code input datasets created or modified by the freezing capability, the total-water and liquid-water saturation functions, the relative permeability function, and spatially varying properties introduced in this version of SUTRA (version 4.0) and described in this report. A complete, up-to-date description of all SUTRA input datasets is included with each SUTRA software release (Voss and Provost, 2024). The SUTRA permeability and dispersion tensors are described in Voss and Provost (2002).

Table 2.1.

List of main input file (.inp) input data created or modified by the newly introduced freezing capability, total-water, and liquid-water saturation functions, relative permeability function, and spatially varying properties for SUTRA model version 4.0.

^{[Input variables of “character” type are enclosed in single quotation marks (‘ ‘)
to provide maximum compatibility across computing platforms. A detailed discussion
of the U.S. Geological Survey Saturated-Unsaturated Transport (SUTRA) code permeability
and dispersion functions is provided at the end of the listing of dataset 15B in Voss and Provost (2002)]}

Table 2.1. List of main input file (.inp) input data created or modified by the newly introduced freezing capability, total-water, and liquid-water saturation functions, relative permeability function, and spatially varying properties for SUTRA model version 4.0.
Variable	Type	Description

Table 2.1. List of main input file (.inp) input data created or modified by the newly introduced freezing capability, total-water, and liquid-water saturation functions, relative permeability function, and spatially varying properties for SUTRA model version 4.0.
Variable	Type	Description
Dataset 2A: Simulation type (one line) [“FREEZING” has been added as an option for the type of transport]
SIMULA	Character	Four words [If the version specification is omitted, SUTRA will read all input datasets in the version 2.0 (2D3D.1) format]. The first word must be “SUTRA.” The second and third words indicate the SUTRA version number and must be “VERSION X.X, “ where X.X is “4.0,” “3.0,” “2.2,” “2.1,” “2.0,” or “2D3D.1.” The fourth word indicates the type of transport and must be either “ENERGY,” FREEZING,” or “SOLUTE.” Freezing is not supported in versions prior to 4.0. Any additional words are ignored by SUTRA. Examples: For an energy-transport simulation without freezing, one may write the following: ‘SUTRA VERSION 4.0 ENERGY TRANSPORT’ For an energy-transport simulation with freezing, one may write the following: ‘SUTRA VERSION 4.0 FREEZING TRANSPORT’ For a solute-transport simulation, one may write the following: ‘SUTRA VERSION 4.0 SOLUTE TRANSPORT’ In these examples, the word “TRANSPORT” is ignored by SUTRA but is included to make the input more readable.
Dataset 8B: Output controls and options for “.nod” file (nodewise results listed in columns; one line) [Options for outputting liquid-water and ice saturations have been added]
NCOLPR	Integer	Nodewise output cycle for a transient simulation. Output is produced in the “.nod” file on time steps numbered n \| NCOLPR \| (where n is a positive integer). In addition, 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, make NCOLPR 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 NCOLPR.
NCOL	Character	List of names of variables to be printed in columns in the “.nod” file. Up to nine columns may be specified. The ordering of columns corresponds to the ordering of variable names in the list. Names may be repeated and may appear in any order, except as noted below. Valid names are as follows: ‘N’=node number (if used, it must appear first in list) ‘X’=x-coordinate of node ‘Y’=y-coordinate of node ‘Z’=z-coordinate of node (3D only) ‘P’=pressure ‘U’=concentration or temperature ‘S’=total water saturation ‘L’=liquid water saturation ‘I’=ice saturation ‘-’=end of list (any names following ‘-’ are ignored) The symbol ‘-’ (a single dash) must be used at the end of the list. Example: To output the 3D node coordinates, pressure, and solute concentration in columns in the “.nod” file every 5 time steps, but not on the first time step, write the following: -5 ‘X’ ‘Y’ ‘Z’ ‘P’ ‘U’ ‘-’
Dataset 8C: Output controls and options for “.ele” file (velocities and Darcy velocities at element centroids listed in columns; one line) [Options for outputting Darcy velocity components have been added]
LCOLPR	Integer	Elementwise output cycle for a transient simulation. Output is produced in the “.ele” file on time steps numbered n \| LCOLPR \| (where n is a positive integer). In addition, for transient solutions, output is produced for initial conditions and on the first and last time steps. For steady-state solutions, output is produced irrespective of the value of LCOLPR, and the velocities are reported only once (for time step 1). Velocities for time step 1 are always reported.
LCOL	Character	List of names of variables to be printed in columns in the “.ele” file. Up to 10 columns may be specified. The ordering of columns corresponds to the ordering of variable names in the list. Names may be repeated and may appear in any order, except as noted below. Valid names are as follows: ‘E’=element number (if used, it must appear first in list) ‘X’=x-coordinate of element centroid ‘Y’=y-coordinate of element centroid ‘Z’=z-coordinate of element centroid (3D only) ‘VX’=x-component of fluid (liquid water) velocity ‘VY’=y-component of fluid (liquid water) velocity ‘VZ’=z-component of fluid (liquid water) velocity (3D only) ‘qX’=x-component of Darcy velocity ‘qY’=y-component of Darcy velocity ‘qZ’=z-component of Darcy velocity (3D only) ‘-’=end of list (any names following ‘-’ are ignored) The symbol ‘-’ (a single dash) must be used at the end of the list. Note: Reported velocities for time step 1 are based on initial or steady-state pressures. Reported velocities for subsequent time steps are based on pressures from the previous time step. Velocities used to formulate the transport equation within SUTRA are based on pressures from the previous nonlinearity iteration; thus, the updated velocities used internally may be different from the values reported for each time step in the “.ele” file. Example: To output the 3D element centroid coordinates and fluid (liquid water) velocity components in columns in the “.ele” file every 10 time steps, write the following: 10 ‘X’ ‘Y’ ‘Z’ ‘VX’ ‘VY’ ‘VZ’ ‘-’
Dataset 9: Liquid water and ice properties (one line) [Input variables COMPI, CI, SIGMAI, and RHOI have been added]
COMPL	Real	Liquid water compressibility, β_L=(1/ρ_L)(∂ρ_L/∂ρ). [M/(L∙s²)]⁻¹.
CL	Real	Liquid water specific heat, c_L=E/(M∙°C) Set to any arbitrary number (for example, zero) for a solute-transport simulation.
SIGMAL	Real	Liquid-water diffusivity, σ_L. For energy transport, represents liquid-water thermal conductivity, λ_L [E/(L‧°C‧s)]. For solute transport, represents molecular diffusivity of solute in pure liquid water, D_m [L²/s]. May be decreased from value in pure water to account for tortuosity of fluid paths. (User-specified constant value for the entire model.)
RHOLØ	Real	For solute transport: Density of liquid water at base concentration for linear liquid density function. [M/L³]. (User-specified constant value for the entire model.) ρ_L=RHOLØ+DRLDU (U−URHOWØ) For energy transport: either – Density of liquid water at base temperature for linear liquid density function. [M/L³], ρ_L=RHOLØ+DRLDU (U−URHOWØ), or – Set to a negative number for SUTRA to use the nonlinear liquid density function. The nonlinear liquid density function is available only for an energy transport simulation. The function calculates liquid density in units of [kg/m³] as a function of temperature in [°C]. The absolute value of RHOLØ can be used to convert units of [kg/m³] to desired units of density: [desired units]=−RHOLØ * [kg/m³] To have SUTRA use density units of [kg/m³], set RHOLØ=−1.0
URHOLØ	Real	For linear liquid density function, represents base value of solute concentration (as mass fraction) or temperature of liquid water at which base liquid density, RHOWØ, is specified. [Ms/M] or [°C]. (User-specified constant value for the entire model.) ρ_L=RHOLØ+DRLDU (U−URHOWØ) For nonlinear liquid density function, an arbitrary value (usually set to 0).
DRLDU	Real	For linear liquid density function, represents coefficient of liquid water density change with concentration (fraction) or temperature [M²/(L³‧M_s)] or [M/(L³‧°C)]. (User-specified constant value for the entire model.) ρ_L=RHOLØ+DRLDU (U−URHOWØ). For nonlinear liquid density function, an arbitrary value (usually set to 0).
VISCØ	Real	For solute transport: Liquid water viscosity, μ [M/(L‧s)]. (User-specified constant value for the entire model.) For energy transport: If set to a positive number, this value is a scale factor. (User-specified constant value for the entire model.) It multiplies the viscosity, which is calculated as a function of temperature in [°C] in units of [kg/(m‧s)]. VISCØ may be used for energy transport to convert units of [kg/(m‧s)] to desired units of viscosity. [desired units]=VISCØ * [kg/(m‧s)] To have SUTRA use viscosity units of [kg/(m‧s)], set VISCØ=1.0 If set to a negative number, the viscosity for energy transport is not calculated as a function of temperature, and the absolute value of this number becomes the user-specified constant viscosity value for the entire model. The units in which the viscosity is expressed are chosen by the user but must be consistent with the units used throughout the model input. [constant viscosity value]=AbsoluteValue(VISCØ)
COMPI	Real	Ice compressibility, β_I=(1/ρ_I)(∂ρ_I/∂p). [M/(L‧s²)]⁻¹ May be omitted for simulations that do not involve FREEZING.
CI	Real	Ice specific heat, c_I [E/(M‧°C)] May be omitted for simulations that do not involve FREEZING.
SIGMAI	Real	Ice thermal conductivity, σ_I ≡ λ_I [E/(L‧°C‧s)] May be omitted for simulations that do not involve FREEZING.
RHOI	Real	Ice density, ρ_I. [M/L³] May be omitted for simulations that do not involve FREEZING.
Dataset 10: Solid matrix properties (one line) [This dataset is no longer used. Solid matrix properties are now spatially variable and are input in datasets 14 and 15]
DATASET 11: Region-wise Parameters (1 line, plus 2 to 6 lines per parameter region, depending on the type of simulation)This dataset has been modified and expanded to include all input parameters that are now specified region-wise. These include parameters for adsorption or thermal conductivity, total water saturation, relative permeability, and liquid water saturation, and the freezing temperature and latent heat.
		Line 1:
NPMREG	Integer	Number of parameter regions, which must be greater than or equal to 1. For each of the NPMREG parameter regions, include datasets 11A–E as necessary, depending on the type of simulation. In other words, all the required datasets for region 1 are listed first, followed by all the required datasets for region 2, etc. Region [NR] numbers must appear in numerical order starting with 1. For each of the NPMREG regions:
NR	Integer	Region number. (Region numbers must appear in numerical order starting with 1.)
Dataset 11A: Adsorption parameters for solute transport simulations (one line per region) [This form of dataset 11A involves adsorption parameters and is used for solute transport simulations only. The form used for energy transport simulations is described further below]
ADSMOD(NR)	Character	Name of the adsorption model for region NR (for solute transport only). For no sorption, set to ‘NONE’. Otherwise, set to ‘LINEAR’, ‘FREUNDLICH’, or ‘LANGMUIR’: For linear sorption model, set to ‘LINEAR’. For Freundlich sorption model, set to ‘FREUNDLICH’. For Langmuir sorption model, set to ‘LANGMUIR’. If ADSMOD(NR)=‘NONE’: Omit the rest of the line. If ADSMOD(NR)=‘LINEAR’, ‘FREUNDLICH’, or ‘LANGMUIR’: (line continued):
CHI1(NR)	Real	Value of linear, Freundlich or Langmuir distribution coefficient for region NR, depending on sorption model chosen in ADSMOD(NR). χ₁. [ $L_{f}^{3} / M_{G}$ ].
CHI2(NR)	Real	Value of Freundlich or Langmuir coefficient for region (NR), depending on sorption model chosen in ADSMOD(NR). Set to any real value for linear sorption. χ₂. Dimensionless quantity for Freundlich. [ $L_{f}^{3} / M_{S}$ ] for Langmuir.
Dataset 11A: Thermal conductivity model for energy transport simulations (one line per region) [This form of dataset 11A involves the thermal conductivity model and is used for energy transport simulations (with or without freezing) only. The form used for solute transport simulations is described above]
TCMOD(NR)	Character	Name of the thermal conductivity model for region NR (for energy transport only). Set to ‘ARITHM’, ‘GEOMET, or ‘HARMON’: For arithmetic-mean bulk thermal conductivity, set to ‘ARITHM’. For geometric-mean bulk thermal conductivity, set to ‘GEOMET’. For harmonic-mean bulk thermal conductivity, set to ‘HARMON’.
Dataset 11B: Total water saturation parameters (one line per region) [This dataset is required only if the simulation is UNSATURATED or involves FREEZING. Omit otherwise]
For no specified total water saturation function (simulated total water saturation will always have a value of 1):
SWMOD(NR)	Character	Set to ‘NONE’.
For the van Genuchten total water saturation function:
SWMOD(NR)	Character	Set to ‘VGEN’.
SWRES(NR)	Real	Value of the residual total water saturation in the van Genuchten function for region NR, S_Wres. [dimensionless quantity]
AA(NR)	Real	van Genuchten function parameter α_VG for region NR. [(L‧s²)/M]
VN(NR)	Real	van Genuchten function parameter n_VG for region NR. [dimensionless quantity]
For the Brooks-Corey total water saturation function:
SWMOD(NR)	Character	Set to ‘BCOR’.
SWRES(NR)	Real	Value of the residual total water saturation in the Brooks-Corey function for region NR, S_Wres [dimensionless quantity]
PENT(NR)	Real	Value of the air-entry pressure (usually has value of zero or is a negative number) in the Brooks-Corey function for region NR, p_ent [M/(L‧s²)]
RLAMB(NR)	Real	Value of the pore size distribution index in the Brooks-Corey function for region NR, λ [dimensionless quantity]
For the piecewise-linear total water saturation function:
SWMOD(NR)	Character	Set to ‘PLIN’.
SWRES(NR)	Real	Value of the residual total water saturation in the piecewise-linear function for region NR, S_Wres. [dimensionless quantity]
PENT(NR)	Real	Value of the air-entry pressure (usually has value of zero or is a negative number) in the piecewise-linear function for region NR, p_ent. [M/(L‧s²)]
PSWRES(NR)	Real	Value of the pressure (usually is a negative number) below which S_W=S_Wres in the piecewise-linear function for region NR, $p_{S_{Wres}}$ . [M/(L‧s²)]
For a user-defined total water saturation function:
SWMOD(NR)	Character	Set to ‘UDEF’. If the user-defined function is used, it must be programmed by the user in the designated section of subroutine UNSAT.
NSWPAR(NR)	Integer	The number of user-defined function parameters specified for region NR, which must be less than or equal to 10.
SWPAR(NR)	Real	List of function parameter values for the user-defined function for region NR. [user-specified units].
Dataset 11C: Relative Permeability Parameters (one line per region) [This dataset is required only if the simulation is UNSATURATED or involves FREEZING. Omit otherwise]
For no specified relative permeability function (simulated relative permeability will always have a value of 1):
RKMOD(NR)	Character	Set to ‘NONE’.
For the van Genuchten relative permeability function:
RKMOD(NR)	Character	Set to ‘VGEN’.
VN(NR)	Real	van Genuchten function parameter n_VG for region NR. [dimensionless quantity]
RKMIN(NR)	Real	Value of the minimum relative permeability, k_rmin, allowed a for region NR. [dimensionless quantity] Setting RKMIN(NR)=0 allows the relative permeability in region NR to approach zero asymptotically as determined by the van Genuchten function.
For the Brooks-Corey relative permeability function:
RKMOD(NR)	Character	Set to ‘BCOR’.
RLAMB(NR)	Real	Value of the pore size distribution index in the Brooks-Corey function for region NR, λ. [dimensionless quantity]
RKMIN(NR)	Real	Value of the minimum relative permeability, k_rmin, allowed a for region NR. [dimensionless quantity] Setting RKMIN(NR)=0 allows the relative permeability in region NR to approach zero asymptotically as determined by the Brooks-Corey function.
For the piecewise-linear relative permeability function:
RKMOD(NR)	Character	Set to ‘PLIN’.
SLRKMIN(NR)	Real	Value of the liquid-water saturation in the piecewise-linear function at which relative permeability reaches its minimum value for region NR, $S_{L_{krmin}}$ . [dimensionless quantity]
RKMIN(NR)	Real	Value of the minimum relative permeability, k_rmin, allowed a for region NR. [dimensionless quantity]
For a user-defined relative permeability function:
RKMOD(NR)	Character	Set to ‘UDEF’. If the user-defined function is used, it must be programmed by the user in the designated section of subroutine RELPERM.
NRKPAR(NR)	Integer	The number of user-defined function parameters specified for region NR, which must be less than or equal to 10.
RKPAR(NR)	Real	List of function parameter values for the user-defined function for region NR. [user-specified units].
Dataset 11D: Liquid Water Saturation Parameters (one line per region) [This dataset is required only if the simulation involves FREEZING. Omit otherwise]
For no specified liquid water saturation function (simulated liquid water saturation at full saturation will always have a value of 1):
SLMOD(NR)	Character	Set to ‘NONE’.
For the exponential liquid water saturation function:
SLMOD(NR)	Character	Set to ‘EXPO’.
SLSATRES(NR)	Real	Value of the residual liquid water saturation in the exponential function for region NR, $S_{Lres}^{sat}$ . [dimensionless quantity]
W (NR)	Real	exponential parameter w_EXP for region NR. [°C]
For the modified power-law liquid water saturation function:
SLMOD(NR)	Character	Set to ‘POWR’.
SLSATRES(NR)	Real	Value of the residual liquid water saturation in the modified power law function for region NR, $S_{Lres}^{sat}$ . [dimensionless quantity]
ALPHA (NR)	Real	Modified power law model parameter, α_POW. [ ${(° C)}^{- β_{POW}}$ ]
BETA (NR)	Real	Modified power law model parameter, β_POW. [dimensionless quantity]
For the piecewise-linear liquid water saturation function:
SWMOD(NR)	Character	Set to ‘PLIN’.
SLSATRES(NR)	Real	Value of the residual liquid water saturation in the piecewise-linear function for region NR, $S_{Lres}^{sat}$ . [dimensionless quantity]
TLRES(NR)	Real	Value of the relative temperature, T−T_f, below which $S_{L}^{sat} = S_{Lres}^{sat}$ in the piecewise-linear function for region NR, $T_{Lres}$ . [M/(L‧s²)]
For a user-defined liquid water saturation function:
SLMOD(NR)	Character	Set to ‘UDEF’. If the user-defined function is used, it must be programmed by the user in the designated section of subroutine LIQSAT.
NSLPAR(NR)	Integer	The number of user-defined function parameters specified for region NR, which must be less than or equal to 10.
SLPAR(NR)	Real	List of function parameter values for the user-defined function for region NR. [user-specified units].
Dataset 11E: Freezing Temperature and Latent Heat (one line per region) [This dataset is required only if the simulation involves FREEZING. Omit otherwise]
TFREEZ(NR)	Real	Maximum freezing temperature of pore water, T_f. Typically, 0 °C, but may be set to, for example, a negative value that is representative of the freezing-point depression caused by solutes dissolved in the groundwater. [°C]
HTLAT(NR)	Real	Latent heat of fusion, ∆H. [E/M]
Dataset 12: Production of Energy or Solute Mass (one or two lines) [This dataset is no longer used. Production parameters are now spatially variable and are input in dataset 14]
Dataset 14A: Scale Factor for Nodewise Data (one line) [Input variables COMPMAF, CSF, RHOSF, PRODLØF, PRODSØF, PRODL1F, PRODS1F, and PRODIØF have been added]
	Character	Line must begin with the word ‘NODE’.
SCALX	Real	The scaled x-coordinates of nodes in dataset 14B are multiplied by SCALX in SUTRA. May be used to change from map to field scales, or from imperial to metric units. A value of 1.0 gives no scaling.
SCALY	Real	The scaled y-coordinates of nodes in dataset 14B are multiplied by SCALY in SUTRA. May be used to change from map to field scales or from imperial to metric units. A value of 1.0 gives no scaling.
SCALZ	Real	For 3D problems, the scaled z-coordinates of nodes in dataset 14B are multiplied by SCALZ in SUTRA. May be used to change from map to field scales or from imperial to metric units. A value of 1.0 gives no scaling. For 2D problems, the scaled element (mesh) thicknesses at nodes in dataset 14B are multiplied by SCALZ in SUTRA. May be used to easily change entire mesh thickness or to convert imperial to metric units. A value of 1.0 gives no scaling.
PORFAC	Real	The scaled nodewise porosities of dataset 14B are multiplied by PORFAC in SUTRA. May be used to easily assign a constant porosity value to all nodes by setting PORFAC=porosity, and all POR(II)=1.0 in dataset 14B.
COMPMAF	Real	The scaled nodewise solid matrix compressibilities of dataset 14B are multiplied by COMPMAF in SUTRA. May be used to easily assign a constant solid matrix compressibility value to all nodes by setting COMPMAF= solid matrix compressibility, and all COMPMA(II)=1.0 in dataset 14B.
CSF	Real	The scaled nodewise solid grain specific heats of dataset 14B are multiplied by CSF in SUTRA. May be used to easily assign a constant solid grain specific heat value to all nodes by setting CSF=solid grain specific heat, and all CS(II)=1.0 in dataset 14B. Set to any arbitrary number (for example, zero) for a solute-transport simulation.
RHOSF	Real	The scaled nodewise densities of a solid grain of dataset 14B are multiplied by RHOSF in SUTRA. May be used to easily assign a constant density of a solid grain value to all nodes by setting RHOSF=density of a solid grain, and all RHOS(II)=1.0 in dataset 14B.
PRODIN	Character	Set to ‘PROD’ if parameters for production or decay of solute or energy are specified; otherwise, set to ‘NOPROD’ to indicate that there is no production or decay.
PRODLØF	Real	The scaled nodewise values of zero-order rate of production in the liquid water (dataset 14B) are multiplied by PRODL0F in SUTRA. May be used to easily assign a constant zero-order rate of production in the liquid water value to all nodes by setting PRODL0F=zero-order rate of production in the liquid water, and all PRODL0(II)=1.0 in dataset 14B. May be omitted if there is no production or decay (PRODIN =’NOPROD’).
PRODSØF	Real	The scaled nodewise values of zero-order rate of production in the immobile phase (dataset 14B) are multiplied by PRODS0F in SUTRA. May be used to easily assign a constant zero-order rate of production in the immobile phase value to all nodes by setting PRODS0F=zero-order rate of production in the immobile phase, and all PRODS0(II)=1.0 in dataset 14B. May be omitted if there is no production or decay (PRODIN =’NOPROD’).
PRODL1F	Real	The scaled nodewise values of first-order rate of solute mass production in the liquid water (dataset 14B) are multiplied by PRODL1F in SUTRA. May be used to easily assign a constant first-order rate of solute mass production in the liquid water value to all nodes by setting PROLF1F=first-order rate of solute mass production in the liquid water, and all PRODL1(II)=1.0 in dataset 14B. Omit for energy transport simulations with or without freezing. May be omitted if there is no production or decay (PRODIN =’NOPROD’).
PRODS1F	Real	The scaled nodewise values of first-order rate of adsorbate mass production in the immobile phase (dataset 14B) are multiplied by PRODS1F in SUTRA. May be used to easily assign a constant first-order rate of adsorbate mass production in the immobile phase value to all nodes by setting PRODS1F=first-order rate of adsorbate mass production in the immobile phase, and all PRODS1(II)=1.0 in dataset 14B. Omit for energy transport simulations with or without freezing. May be omitted if there is no production or decay (PRODIN =’NOPROD’).
PRODIØF	Real	For energy transport simulations with FREEZING, the scaled nodewise values of zero-order rate of production in the ice (dataset 14B) are multiplied by PRODI0F in SUTRA. May be used to easily assign a constant zero-order rate of production in the ice value to all nodes by setting PRODI0F=zero-order rate of production in the ice, and all PRODI0(II)=1.0 in dataset 14B. Omit for energy transport simulations with or without freezing. May be omitted if there is no production or decay (PRODIN =’NOPROD’).
Dataset 14B: Nodewise Data (one line for each of NN nodes) [Input variables COMPMA, CS, RHOS, PRODLØ, PRODSØ, PRODL1, PRODS1, and PRODIØ have been added. When the DIRECT solver is used, the order in which the nodes are numbered affects the bandwidth of the global banded matrix, NBI, which in turn affects computational and storage efficiency. In this case, the user should take care to number the nodes to minimize NBI. SUTRA sets NBI equal to one plus twice the maximum difference in node numbers in the element containing the largest node number difference in the mesh. See Voss and Provost (2002, fig. 7.1) for an example. When an iterative solver is used, it is still advantageous to minimize NBI, although not as critical as in the case of the DIRECT solver]
II	Integer	Number of node, i, to which data on this line refers.
NREG (II)	Integer	Region number to which node II belongs. (1 ≤ NREG ≤ NPMREG; see NPMREG in dataset 11.) The node region number is usually the same as the region number of the elements in which it appears. When the node is at the boundary of two regions, it may be assigned to either region.
X(II)	Real	Scaled x-coordinate of node II, x_i. [L]
Y(II)	Real	Scaled y-coordinate of node II, y_i. [L]
Z(II)	Real	For 3D problems, scaled z-coordinate of node II, z_i. [L] For 2D problems, scaled thickness of mesh at node II. [L] To simulate radial cross sections, set Z(II)=(2π)(radius_i), where radius_i is the radial distance from the vertical center axis (axis of radial symmetry) to node i.
POR(II)	Real	Scaled porosity value at node II, ε_i. [dimensionless quantity]
COMPMA(II)	Real	Solid matrix compressibility at node II, $α_{S_{i}} = {(1 - ε_{i})}^{- 1} (\partial ε_{i} / \partial p)$ . [M/(L‧s²)]⁻¹.
CS(II)	Real	Solid grain specific heat value at node II, $c_{S_{i}}$ . [E/(M‧°C)] Set to any arbitrary number (for example, zero) for a solute-transport simulation.
RHOS(II)	Real	Density of a solid grain value at node II, $ρ_{S_{i}}$ . [M/L³]. Value used only for energy-transport simulations or solute-transport simulations with sorption; may otherwise be set to any arbitrary number (for example, zero).
PRODLØ(II)	Real	Zero-order rate of production in the liquid water at node II, $γ_{0_{i}}^{L}$ . [(E/M)/s] for energy production, [(Ms/M)/s] for solute mass production.
PRODSØ(II)	Real	Zero-order rate of production in the immobile phase at node II, $γ_{0_{i}}^{S}$ . [(E/MG)/s] for energy production, [(Ms/MG)/s] for adsorbate mass production.
PRODL1(II)	Real	First-order rate of solute mass production in the liquid water at node II, $γ_{1_{i}}^{L}$ . [s⁻¹]. Omit for energy transport simulations with or without freezing.
PRODS1(II)	Real	First-order rate of adsorbate mass production in the immobile phase at node II, $γ_{1_{i}}^{S}$ . [s⁻¹]. Omit for energy transport simulations with or without freezing.
PRODIØ(II)	Real	Zero-order rate of production in the ice at node II, $γ_{0_{i}}^{I}$ . [(E/M)/s]. Omit for energy transport simulations with or without freezing.
Dataset 15A: Scale Factors for Elementwise Data (one line) [Input variables SIGMASF and SIGMAAF have been added]
	Character	Line must begin with the word ‘ELEMENT’.
PMAXFA	Real	The scaled maximum permeability values, PMAX, in dataset 15B are multiplied by PMAXFA in SUTRA. May be used to convert units or to aid in assignment of maximum permeability values in elements.
PMIDFA	Real	The scaled middle permeability values, PMID, in dataset 15B are multiplied by PMIDFA in SUTRA. May be used to convert units or to aid in assignment of maximum permeability values in elements. Omit for 2D problems.
PMINFA	Real	The scaled minimum permeability values, PMIN, in dataset 15B are multiplied by PMINFA in SUTRA. May be used to convert units or to aid assignment of minimum permeability values in elements.
ANG1FA	Real	The scaled angles ANGLE1 in dataset 15B are multiplied by ANG1FA in SUTRA. For 2D problems, may be used to easily assign a uniform direction of anisotropy by setting ANG1FA=angle, and all ANGLE1(L)=1.0 in dataset 15B.
ANG2FA	Real	The scaled angles ANGLE2 in dataset 15B are multiplied by ANG2FA in SUTRA. Omit for 2D problems.
ANG3FA	Real	The scaled angles ANGLE3 in dataset 15B are multiplied by ANG3FA in SUTRA. Omit for 2D problems.
ALMAXF	Real	The scaled longitudinal dispersivities ALMAX in dataset 15B are multiplied by ALMAXF in SUTRA. May be used to convert units or to aid in assignment of dispersivities.
ALMIDF	Real	The scaled longitudinal dispersivities ALMID in dataset 15B are multiplied by ALMIDF in SUTRA. May be used to convert units or to aid in assignment of dispersivities. Omit for 2D problems.
ALMINF	Real	The scaled longitudinal dispersivities ALMIN in dataset 15B are multiplied by ALMINF in SUTRA. May be used to convert units or to aid in assignment of dispersivities.
ATMAXF	Real	The scaled transverse dispersivities ATMAX in dataset 15B are multiplied by ATMAXF in SUTRA. May be used to convert units or to aid in assignment of dispersivity.
ATMIDF	Real	The scaled transverse dispersivities ATMID in dataset 15B are multiplied by ATMIDF in SUTRA. May be used to convert units or to aid in assignment of dispersivity. Omit for 2D problems.
ATMINF	Real	The scaled transverse dispersivities ATMIN in dataset 15B are multiplied by ATMINF in SUTRA. May be used to convert units or to aid in assignment of dispersivity.
SIGMASF	Real	The scaled solid grain thermal conductivity values, SIGMAS, in dataset 15B are multiplied by SIGMASF in SUTRA. May be used to convert units or to aid in assignment of solid grain thermal conductivity values in elements. May be omitted for solute transport simulations.
SIGMAAF	Real	The scaled effective air thermal conductivity values, SIGMAA, in dataset 15B are multiplied by SIGMAAF in SUTRA. May be used to convert units or to aid in assignment of effective air thermal conductivity values in elements. May be omitted for solute transport simulations.
Dataset 15B: Elementwise Data (one line for each of NE elements) [Input variables SIGMAS and SIGMAA been added]
L	Integer	Number of element to which data on this line refers.
LREG(L)	Integer	Region number to which this element belongs. (1 ≤ LREG ≤ NPMREG; see NPMREG in dataset 11.)
PMAX(L)	Real	Scaled maximum permeability value, k_max, of element L. [L²]
PMID(L)	Real	Scaled middle permeability value, k_mid, of element L. [L²] Isotropic permeability requires: PMID(L)=PMAX(L). Omit for 2D problems.
PMIN(L)	Real	Scaled minimum permeability value, k_min, of element L. [L²] Isotropic permeability requires: PMIN(L)=PMAX(L).
ANGLE1(L)	Real	Scaled first angle, θ₁ [°], used to define the directions of maximum, middle, and minimum permeability in element L.
ANGLE2(L)	Real	Scaled second angle, θ₂ [°], used to define the directions of maximum, middle, and minimum permeability in element L. Omit for 2D problems.
ANGLE3(L)	Real	Scaled third angle, θ₃ [°], used to define the directions of maximum, middle, and minimum permeability in element L. Omit for 2D problems.
ALMAX(L)	Realsoft	Scaled longitudinal dispersivity value, α_Lmax, of element L that controls longitudinal dispersion along the maximum permeability direction when the flow direction is in the maximum permeability direction. [L]
ALMID(L)	Real	Scaled longitudinal dispersivity value, α_Lmid, of element L that controls longitudinal dispersion along the middle permeability direction when the flow direction is in the middle permeability direction. [L] Omit for 2D problems.
ALMIN(L)	Real	Scaled longitudinal dispersivity value, α_Lmin, of element L that controls longitudinal dispersion along the minimum permeability direction when the flow direction is in the minimum permeability direction. [L]
ATMAX(L)	Real	Scaled transverse dispersivity value, α_Tmax, of element L that controls transverse dispersion in the maximum permeability direction when the flow direction is perpendicular to the maximum permeability direction. [L]
ATMID(L)	Real	Scaled transverse dispersivity value, α_Tmid, of element L that controls transverse dispersion in the middle permeability direction when the flow direction is perpendicular to the middle permeability direction. [L] Omit for 2D problems.
ATMIN(L)	Real	Scaled transverse dispersivity value, α_Tmin, of element L that controls transverse dispersion in the minimum permeability direction when the flow direction is perpendicular to the minimum permeability direction. [L]
SIGMAS(L)	Real	Scaled solid grain thermal conductivity value, σ_S≡λ_S, of element L. [E/(L‧°C‧s)] May be omitted for solute transport simulations.
SIGMAA(L)	Real	Scaled effective air thermal conductivity value, σ_A≡λ_S, of element L. [E/(L‧°C‧s)] May be omitted for solute transport simulations.

References Cited

Voss, C.I., and Provost, A.M., 2002, SUTRA—A model for 2D or 3D saturated-unsaturated, variable-density ground-water flow with solute or energy transport [version 2D3D.1, September 22, 2010]: U.S. Geological Survey Water-Resources Investigations Report 02–4231, 291 p., accessed March 23, 2022, at https://doi.org/10.3133/wri024231.

Voss, C.I., and Provost, A.M., 2024, SUTRA—A model for 2D or 3D saturated-unsaturated, variable-density groundwater flow with solute or energy transport, version 4.0: U.S. Geological Survey software release, https://doi.org/10.5066/P9OL5IYX.

Conversion Factors

International System of Units to U.S. customary units


Multiply	By	To obtain
meter (m)	3.281	foot (ft)
square meter (m²)	10.76	square foot (ft²)
meter per second (m/s)	3.281	foot per second (ft/s)
kilogram (kg)	2.205	pound avoirdupois (lb)
kilopascal (kPa)	0.009869	atmosphere, standard (atm)
kilogram per cubic meter (kg/m³)	0.06242	pound per cubic foot (lb/ft³)
joule (J)	0.0000002	kilowatthour (kWh)
meter per day (m/d)	3.281	foot per day (ft/d)

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

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

Abbreviations

[Abbreviations related to equations are listed in appendix 1]

2D: two-dimensional
3D: three-dimensional
DOE: U.S. Department of Energy
EPA: U.S. Environmental Protection Agency
SERDP: Strategic Environmental Research and Development Program
SUTRA: Saturated-Unsaturated Transport [code]
TSD: Thiesen-Scheel-Diesselhorst [equation]
USGS: U.S. Geological Survey

For more information, contact

Director, Earth System Processes Division

U.S. Geological Survey

Mail Stop 411

12201 Sunrise Valley Drive

Reston, VA 20192

Publishing support provided by the Pembroke and Reston Publishing Service Centers

Disclaimers

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

Voss, C.I., Provost, A.M., McKenzie, J.M., and Kurylyk, B.L., 2024, SUTRA—A code for simulation of saturated-unsaturated, variable-density groundwater flow with solute or energy transport—Documentation of the version 4.0 enhancements—Freeze-thaw capability, saturation and relative-permeability relations, spatially varying properties, and enhanced budget and velocity outputs: U.S. Geological Survey Techniques and Methods, book 6, chap. A63, 91 p., https://doi.org/10.3133/tm6A63.

ISSN: 2328-7055 (online)

Additional publication details
Publication type	Report
Publication Subtype	USGS Numbered Series
Title	SUTRA— A code for simulation of saturated-unsaturated, variable-density groundwater flow with solute or energy transport—Documentation of the version 4.0 enhancements—Freeze-thaw capability, saturation and relative-permeability relations, spatially varying properties, and enhanced budget and velocity outputs
Series title	Techniques and Methods
Series number	6-A63
DOI	10.3133/tm6A63
Year Published	2024
Language	English
Publisher	U.S. Geological Survey
Publisher location	Reston, VA
Contributing office(s)	WMA - Earth System Processes Division
Description	Report: vii, 91 p.; Software Release
Online Only (Y/N)	Y
Additional Online Files (Y/N)	N

Table of Contents

Links

Acknowledgments

Abstract

Introduction

1. Freeze-Thaw Capability

1.1. Water-Mass Balance With Freeze-Thaw

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

1.1.1. Assumption Regarding Pressure Effect of Liquid-Ice Density Difference in Water-Mass Balance

(21)

(22)

(23)

1.2. Energy Balance With Freeze-Thaw

(24)

(25)

(26A)

(26B)

(26C)

(27)

(28A)

(28B)

(28C)

(29)

(30)

(31)

(32)

(33)

(34)

1.2.1. Assumption Regarding Pressure Effect of Liquid-Ice Density Difference in Energy Balance

(35)

1.3. Freeze-Thaw Relations

1.3.1. Freeze-Thaw Under Fully Saturated Conditions

(36)

(37)

1.3.2. Freeze-Thaw Under Unsaturated Conditions

(38)

(39)

(40)

1.3.3. Permeability for Freeze-Thaw Under Saturated and Unsaturated Conditions

(41)

1.3.4. Behavior of the Freeze-Thaw Model

1.3.5. General Procedure for Saturation Calculations

1.4. Thermal Conductivity Functions

(42A)

(42B)

(42C)

(43)

(44)

1.5. Nonlinear Liquid-Water Density Function

(45)

(46)

2. Saturation and Relative-Permeability Functions

2.1. Overview of Available Functions

2.2. Preprogrammed Total-Water Saturation Functions for Unsaturated Conditions

2.2.1. Van Genuchten Total-Water Saturation Function

(47)

(48)

2.2.2. Brooks-Corey Total-Water Saturation Function

(49)

(50)

2.2.3. Piecewise-Linear Total-Water Saturation Function