Water Flow and Heat Transport in Frozen Soil: Numerical Solution
and Freeze–Thaw Applications
Klas Hansson,* Jirka Šimůnek, Masaru Mizoguchi, Lars-Christer Lundin, and Martinus Th. van Genuchten

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACTsources in Alaska and Canada was a major driving force
for studies of freezing–thawing cycles in the 1970s and1980s (Lundin, 1989).
Many of the processes that we try to formulate mathematically and simulate numerically today were alreadybeing studied and described some 70 yr ago. For example, Beskow (1935) studied the effects of freezing andthawing on railroads and roads and observed three fundamental processes typical of freezing soils. One observation was that water flows toward, and accumulates at,freezing fronts. A second observation was that water insoil pores does not freeze at 0⬚C, but is subject to afreezing-point depression caused by interactions between water, soil particles, and solutes. Beskow also
observed that soil freezing is generally quite similar tosoil drying. Edlefsen and Andersen (1943) later tried todescribe the mutual dependence of temperature, watercontent, and solutes by means of a generalized and extended Clapeyron equation using thermodynamic equilibrium theory. The original form of the Clapeyron
equation, which relates changes in pressures and temperatures, was formulated for one-component equilibrium between two phases at the same pressure (Alberty
and Silbey, 1992), in our case pure ice and liquid water.
In soil science it is customary to modify the original
Clapeyron equation, as used by Edlefsen and Andersen(1943), by including the osmotic pressure or assuming
a difference in pressures between ice and water. Specifically, the ice pressure is sometimes assumed to equal
the zero gauge pressure, with the reference pressure being atmospheric. While this assumption has often been
debated, no consensus has yet been reached. As Spaans
and Baker (1996) wrote: “the broad assumption of zerogauge pressure in the ice phase has been questioned
under certain conditions (Miller, 1973, 1980), but thus
far there is scant evidence against it, except in obvious
cases (heaving).” In particular, if a soil is unsaturated
the potential of heaving is reduced such that the assumption of zero ice pressure is more likely to hold.
Koopmans and Miller (1966) showed experimentallythat freezing curves are similar to soil water drainage(or soil water retention) curves through a scaling relationship between ice pressure (which in turn was seen as
a function of temperature alone by using the Clapeyron
equation) and air pressure (or capillary pressure). The
theory in their article applied to saturated soils that areeither free of colloidal material (type SS [Solid-to-Solid],e.g., sand, silt, or coarse clay fractions), or soils in which
the particles are always surrounded by water and thus
separated from each other (type SLS, Solid-LiquidSolid). Miller (1978) later developed a rigid ice modelthat focused on ice lens formation and frost heave in

A new method is presented to account for phase changes in a fully
implicit numerical model for coupled heat transport and variably
saturated water flow involving conditions both above and below zerotemperature. The method is based on a mixed formulation for bothwater flow and heat transport similar to the approach commonly used
for the Richards equation. The approach enabled numerically stable,energy- and mass-conservative solutions. The model was evaluatedby comparing predictions with data from laboratory column freezing
experiments. These experiments involved 20-cm long soil columns
with an internal diameter of 8 cm that were exposed at the top to a circulating fluid with a temperature of ⫺6ⴗC. Water and soil in the columns froze from the top down during the experiment, with the freezing
process inducing significant water redistribution within the soil. A
new function is proposed to better describe the dependency of thethermal conductivity on the ice and water contents of frozen soils.
Predicted values of the total water content compared well with measured values. The model proved to be numerically stable also for ahypothetical road problem involving simultaneous heat transport andwater flow. The problem was simulated using measured values of thesurface temperature for the duration of almost 1 yr. Since the roadwas snow-plowed during winter, surface temperatures varied more
rapidly, and reached much lower values, than would have been thecase under a natural snow cover. The numerical experiments demonstrate the ability of the code to cope with rapidly changing boundary
conditions and very nonlinear water content and pressure head distributions in the soil profile.

T

he importance of freezing and thawing processes
in soils has long been recognized. Much attention
initially focused on the problems of frost heave becauseof its importance in the construction and maintenanceof roads, railroads, and oil industry pipelines. Recent
applications deal with both environmental and engineering issues. For example, the often dramatically reduced infiltration capacity of frozen soils (e.g., Seyfried
and Murdock, 1997) may increase soil erosion, and evenflooding, due to increased snow melt surface runoff. In
some engineering applications soil is artificially frozen
to prevent the migration of pollutants (Andersland etal., 1995), increase soil structural stability, or prevent
leaching of water such as during tunnel construction(Jones, 1981). While frost-related damage of roads represented an important issue in Scandinavia and North
America in the 1920s, exploitation of oil and gas reK. Hansson and L.-C. Lundin, Dep. of Earth Sciences, Uppsala Univ.,
Uppsala, Sweden; J. Šimůnek, Dep. of Environmental Science, Univ.
of California, Riverside, CA; M. Mizoguchi, Dep. of Biological andEnvironmental Engineering, Univ. of Tokyo, Japan; M. Th. van Genuchten, George E. Brown, Jr. Salinity Laboratory, USDA-ARS, Riverside, CA. Received 19 Aug. 2003. Original Research Paper. *Corresponding author (klas.hansson@hyd.uu.se).
Published in Vadose Zone Journal 3:693–704 (2004).
 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

Abbreviations: SVAT, soil–vegetation–atmosphere transfer.

693

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

694

VADOSE ZONE J., VOL. 3, MAY 2004

non-colloidal soils. Ice lenses represent situations where
the ice pressure is not equal to zero. Fuchs et al. (1978)
included a water retention function in their analysis of
the freezing process, thus extending the theory to unsaturated conditions. Using thermodynamic principles,
Kay and Groenevelt (1974) and Groenevelt and Kay(1974) developed several generalized Clapeyron equations, as well as the basic theory describing coupled waterand heat transport in frozen soils. A generalized Clapeyron equation was later derived from a more traditional physical chemistry perspective by Loch (1978).
A potential problem with the equilibrium assumptionof this theory is that equilibrium may never be reached,or may take a very long time to establish (Spaans andBaker, 1996). We note that the Clapeyron equation sometimes is referred to as the Clausius–Clapeyron equation,even though the latter is a simplified version of theClapeyron equation applicable only to equilibria between gases and liquids (e.g., Alberty and Silbey, 1992).
The introduction of computers during the 1970s enabled the development of numerical codes specificallydesigned to solve the heat and water flow equations,
whether coupled or not. Harlan (1973) presented amodel for coupled transport of heat and fluid in partiallyfrozen soils, and to our knowledge was the first to givea numerical solution to the problem using finite differences. He implemented a formulation that involved theapparent heat capacity as proposed by Lukianov andGolovko (Fuchs et al., 1978). Guymon and Luthin(1974) presented a finite element solution to the sameproblem. Neither Harlan (1973) nor Guymon andLuthin (1974) in their analyses considered the effects
of solutes, which are known to have a significant effecton the amount of unfrozen water at subzero temperatures. Cary and Mayland (1972) and Fuchs et al. (1978)were among the first to introduce models incorporating
osmotic effects on freezing processes in soils.
In the late 1970s, one-dimensional numerical soil–vegetation–atmosphere–transfer (SVAT) models emerged,
such as the SOIL model of Jansson and Halldin (1980)
and the model by Lykosov and Palagin (1980). Thesemodels combined advanced treatments of the atmospheric surface layer with soil profile simulations. They
typically included processes in the snow cover, which is
a complicated problem by itself in part because of thepresence of a continuously deforming porous medium
constituting the snow pack. Conditions in the upper soil
layers of natural field soils are more or less uniquelydetermined by external forces exerted on the soil by
prevailing climatic conditions in the atmospheric surface
layer. We note that SVAT models are now increasingly
used to define the lower boundary condition of global
circulations models (e.g., van den Hurk et al., 1995).
Several one-dimensional numerical codes currently
exist for simulating water and heat transport, including
freezing and thawing. Most of these codes are for specific applications, such as the ICM (Lytton et al., 1993)
and FrostB (Guymon et al., 1993) models designed forroads, airstrips, or similar environments. Other, moregeneral, models include the SHAW model (Flerchinger
and Saxton, 1989; Flerchinger et al., 1996) and the COUP

model (Jansson and Karlberg, 2001), the latter being an
extension of the SOIL code mentioned earlier.
The primary objective of this paper is to present a
new, numerically stable, mass-and energy-conservative
method for dealing with phase changes of water in a fully
implicit numerical solution of the coupled equations
governing heat transport and variably saturated flow,
for conditions both above and below zero temperature.
The method is as an extension of the procedure proposed by Celia et al. (1990) for variably saturated flow
to coupled water and energy transport. The numericalmodel was evaluated by comparisons with short-termlaboratory column freezing experiment data. The numerical stability of the model was further tested bymeans of a seasonal simulation of a hypothetical road
construction problem involving heat transport and water flow in response to measured surface temperatures.
In this simulation we investigate how the code predicts
the dynamics of naturally varying climatic conditions,
including thawing behavior. Since the ultimate goal is
to create a two-dimensional model that may be used to
predict road moisture and temperature distributions, it
was felt important to ensure that the code performed
well for relatively extreme upper temperature boundary
conditions. Even though the code was designed to describe processes for both positive and negative temperatures, this paper is limited to frost related processes
since other parts of the code are described in detailelsewhere (Šimůnek et al., 1998; Scanlon et al., 2003).
THEORY
Variably Saturated Flow
Variably saturated water flow for above- and subzero temperatures is described using the modified Richards equation
as follows (e.g., Fayer, 2000; Noborio et al., 1996):⳵␪u(h)
␳ ⳵␪ (T )
⫹ i i
⫽
⳵t
␳w ⳵t
⳵
⳵h
⳵T
⳵h
⳵T
KLh(h) ⫹ KLh(h) ⫹ KLT(h) ⫹ Kvh(␪) ⫹ KvT(␪)
⫺S
⳵z
⳵z
⳵z
⳵z
⳵z

冤

冥

[1]

where ␪u is the volumetric unfrozen water content (L L⫺3)(⫽␪ ⫹ ␪v), ␪ is the volumetric liquid water content (L3 L⫺3),
␪v is the volumetric vapor content expressed as an equivalentwater content (L3 L⫺3), ␪i is the volumetric ice content (L3L⫺3), t is time (T), z is the spatial coordinate positive upward(L), ␳i is the density of ice (M L⫺3) (931 kg m⫺3), ␳w is thedensity of liquid water (M L⫺3) (approximately 1000 kg m⫺3),
h is the pressure head (L), T is the temperature (K), and S isa sink term (T⫺1) usually accounting for root water uptake.
Water flow in Eq. [1] is assumed to be caused by fivedifferent processes with corresponding hydraulic conductivities (given below in parentheses). The first three terms on theright-hand side of Eq. [1] represent liquid flows due to apressure head gradient (KLh, [L T⫺1]), gravity, and a temperature gradient (KLT, [L2 T⫺1 K⫺1]), respectively. The next two
terms represent vapor flows due to pressure head (Kvh, [L T⫺1])
and temperature (KvT, [L2 T⫺1 K⫺1]) gradients, respectively.
Equation [1] is highly nonlinear, mainly due to dependencies of the water content and the hydraulic conductivity on
the pressure head, that is, ␪(h ) and KLh(h ), respectively, and3

695

www.vadosezonejournal.org

due to freezing–thawing effects that relate the ice content withthe temperature, that is, ␪i(T ).

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Soil Hydraulic Properties
While different functions for the unsaturated soil hydraulicproperties may be used, we will in this study invoke the expressions of van Genuchten (1980), Mualem (1976), and van Genuchten et al. (1991) with independent m and n parameters:

Se(h) ⫽

␪(h) ⫺ ␪r⫽ (1 ⫹ |␣h|n)⫺m␪s ⫺ ␪ r

KLh(Se) ⫽ KsSle[I␨ (p, q)]2

[2]
[3]

where Se is effective saturation, ␪r and ␪s denote the residual
and saturated water contents (L3 L⫺3), respectively, Ks is thesaturated hydraulic conductivity (L T⫺1), and ␣ (L⫺1), n, m,
and l are empirical parameters. Further, I␵ is the incomplete
beta-function, ␨ ⫽ ␪ l/m
e , p ⫽ m ⫹ 1/n, and q ⫽ 1 ⫺ 1/n. If m ⫽1 ⫺ 1/n, Eq. [3] reduces to the commonly used form:m 2
KLh(Se) ⫽ KsS le[1 ⫺ (1 ⫺ S l/m
e ) ]

[4]

The hydraulic conductivity KLT for liquid phase fluxes due
to a gradient in T is defined as (e.g., Fayer, 2000; Noborio
et al., 1996):冢

1 d␥
KLT ⫽ KLh hGwT
␥o dT

冣

[5]

where GwT is the gain factor, ␥ is the surface tension of soil
water (M T⫺2), evaluated as ␥ ⫽ 75.6 ⫺ 0.1425T ⫺ 2.38 ⫻10⫺4T 2, and ␥0(25⬚C) ⫽ 71.89 g s⫺2. The isothermal, Kvh, andthermal, KvT, vapor hydraulic conductivities in Eq. [1] aredefined, among others, by Fayer (2000) and Scanlon et al.
(2003). Since vapor transport is not a primary objective of thispaper, we will not further discuss these two conductivities.
When the soil is frozen, the presence of ice (instead ofliquid water) in some pores may significantly increase theresistance of the porous medium to water flow and lead to an
apparent blocking effect. To account for this blocking, thehydraulic conductivity is often reduced by means of an impedance factor, ⍀ (Lundin, 1990), which in our study is multiplied
by Q, the ratio of the ice content to the total (minus theresidual) water content. The parameter Q accounts for thefact that blocking becomes more effective as the ice content
part of the total water content increases. The impedance factor
reduces the hydraulic conductivity for the liquid phase of thepartially frozen soil, KfLh, as follows

KfLh ⫽ 10⫺⍀Q KLh[6]which shows that even a small value of ⍀ can have a significant
effect on the conductivity of the liquid phase as the ice portion increases.

Heat Transport
Heat transport during transient flow in a variably saturated
porous medium is described as (e.g., Nassar and Horton,
1989, 1992):

⳵CpT
⳵␪⳵␪ (T)⫺ Lf ␳i i ⫹ L0(T) v
⫽
⳵t
⳵t
⳵t
⳵
⳵q T
⳵q T
⳵T
⳵q
␭(␪)
⫺ Cw l ⫺ Cv v ⫺ L0(T) v ⫺ CwST
⳵z
⳵z
⳵z
⳵z
⳵z
[7]

冤

冥

where the first term on the left-hand side represents changes inthe energy content, and the second and third terms representschanges in the latent heat of the frozen and vapor phases,
respectively. The terms on the right-hand side represent, respectively, soil heat flow by conduction, convection of sensible
heat with flowing water, transfer of sensible heat by diffusionof water vapor, transfer of latent heat by diffusion of water
vapor, and uptake of energy associated with root water uptake.
The volumetric heat capacity of the soil, Cp (J m⫺3 K⫺1, M L⫺1
T⫺2 K⫺1), in Eq. [7] is defined as the sum of the volumetricheat capacities of the solid (Cn), liquid (Cw), vapor (Cv), andice (Ci), phases multiplied by their respective volumetric fractions ␪:

Cp ⫽ Cn␪n ⫹ Cw␪ ⫹ Cv␪v ⫹ Ci␪i

[8]

Furthermore, in Eq. [7], L0 is the volumetric latent heat ofvaporization of water (J m⫺3, M L⫺1 T⫺2), L0 ⫽ Lw␳w, Lw is
the latent heat of vaporization of water (J kg⫺1) (⫽2.501 ⫻106 ⫺ 2369.2 T [⬚C]), and Lf is the latent heat of freezing(J kg⫺1, L2 T⫺2) (approximately 3.34 ⫻ 105 J kg⫺1).
The first two terms of Eq. [7] are often combined as follows:

⳵CpT
⳵␪
⳵C T
d␪ ⳵T
⫺ L f ␳i i ⫽ p ⫺ L f ␳i i
⳵t
⳵t
⳵t
dT ⳵t
d␪
⳵ Cp ⫺ L f ␳ i i T
⳵C T
dT
⫽ a
⫽
⳵t
⳵t

冢

冣

[9]

which, after incorporating the latent heat of fusion, leads tothe following definition of the apparent volumetric heat capacity, Ca (J m⫺3 K⫺1, M L⫺1 T⫺2 K⫺1):

Ca ⫽ C p ⫺ L f ␳ i

d␪i
dT

[10]

Assuming that the ice gauge pressure is zero, and that theosmotic pressure is zero, the relationship between capillarypressure and temperature is defined by the generalized Clapeyron equation (e.g., Williams and Smith, 1989):

L
dP
⫽ f
dT VwT[11]

where P is the pressure (Pa, M L⫺1 T⫺2) (⫽ ␳wgh), and Vw isthe specific volume of water (L3 M⫺1) (approximately 0.001 m3
kg⫺1). Using the generalized Clapeyron equation, the apparent
volumetric specific heat, Ca, can be redefined using the hydraulic (or soil moisture) capacity, C (T⫺1), which is the slope ofthe retention curve, defined as the derivative of the watercontent with respect to the pressure head:Ca ⫽ C p ⫺ L f ␳ id␪i
L2␳ d␪
L2 ␳⫽ Cp ⫹ f i⫽ Cp ⫹ f iC
dT
gT dh
gT
[12]

Equations [1] and [7] are tightly coupled due to their mutual
dependence on water contents, pressure heads, and temperatures. Equation [7] is highly nonlinear because of the dependency of the apparent volumetric heat capacity, Ca, on temperature, which is shown in Fig. 1 for three soil textural classes,that is, sand, loam, and silty clay (Carsel and Parrish, 1988).
Notice that while the apparent heat capacity increases by
about two orders of magnitude for silty clay when the freezing
point is reached, the increase for sand is almost five ordersof magnitude (Fig. 1). For sands the increase in the apparent
heat capacity due to freezing becomes negligible below about⫺0.05⬚C, which corresponds to a pressure head of ⫺6.2 m. At
that temperature almost all soil water is frozen, except for the
residual water content. For fine-textured soils the increase inReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

696

VADOSE ZONE J., VOL. 3, MAY 2004

Fig. 1. Apparent volumetric heat capacity, Ca (J m⫺3 K⫺1), for three soil textural classes (sand, loam, and silty clay) as compiled by Carsel andParrish (1988).

the apparent heat capacity extends to much lower temperatures, which reflects the fact that for these soils a significant
amount of water remains unfrozen at slightly subzero temperatures. The apparent heat capacities, Ca, shown in Fig. 1 were
calculated for fully saturated soils. For unsaturated soils, Cavalues are more or less constant above a small negative temperature that can be calculated using the generalized Clapeyron Eq. [11] from the pressure head corresponding to theactual water content. Thus, when the actual soil water content
corresponds to a pressure head of about ⫺120 m, soil water
starts freezing only at ⫺1⬚C, which circumvents the enormous
increase in the apparent heat capacity and eliminates the corresponding nonlinearity in Eq. [7].
According to Campbell (1985), the apparent thermal conductivity for unfrozen conditions, ␭(␪) (W m⫺1 K⫺1, M L T⫺3
K⫺1), can be described by␭(␪) ⫽ ␭0(␪) ⫹ ␤tCw|qw|冦␭ (␪) ⫽ C ⫹ C ␪ ⫺ (C ⫺ C )exp[⫺(C ␪) ]
0

1

2

1

4

3

C5

thermal dispersivity (L), and Ci (i ⫽ 1, . . ., 5) are constants
that can be estimated experimentally or derived from material
properties such as the volume fraction of solids. Campbell(1985) suggested that C5 ⫽ 4. Equation [13b] cannot be used
for subzero temperatures since ice conducts heat much better
than water; hence, only the liquid water content is includedin the present equation. Three models for the thermal conductivity of a frozen soil were discussed in Kennedy and Sharratt(1998). In the first model the thermal conductivity increases
linearly with ice content, in the second model it increases
nonlinearly with ice content, and in the third model it is constant. Hence, consensus is yet to be reached. Given the nonlinear nature of the thermal conductivity of the soil used in ourstudy (Fig. 2), we modified Eq. [13b] by replacing ␪, with ␪ ⫹F␪i, where F is given byF ⫽ 1 ⫹ F1␪iF2

[13a]
[13b]

where ␭0(␪) is the thermal conductivity of the porous medium(solid plus water) in the absence of flow, ␤t is the longitudinal[14]

so that␭0(␪) ⫽ C1 ⫹ C2(␪ ⫹ F␪i) ⫺ (C1 ⫺ C4)
exp{⫺[C3(␪ ⫹ F␪i)]C5}.

[15]

697

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

Fig. 2. Measured thermal conductivity, ␭0, of Kanagawa sandy loam soil (symbols) as a function of temperature and water content (Mizoguchi,1990), as well as its parameterizations used in numerical simulations (lines) using Eq. [15]. The upper graph shows the thermal conductivity
for frozen conditions at total volumetric water contents of 0.20, 0.30, and 0.40, and the lower graph the thermal conductivity of unfrozen soil
at various water contents.

Thus, F accounts for the difference between the thermal conductivities of ice and water in soils. In this way the ice cancontribute to the soil thermal conductivity for frozen conditions, which makes a significant difference for relatively low
temperatures when the liquid water content is small and theice content relatively high.

Numerical Implementation
The coupled water, heat, and vapor transport processes described above were incorporated in the numerical HYDRUS-1D code (Šimůnek et al., 1998). Equation[1], subject to appropriate initial and boundary conditions, was solved numerically using finite differences forboth spatial and temporal discretization. The employedsolution scheme is an extension of the mass–conservative iterative numerical scheme used by Celia et al. (1990).
Finite element and finite difference methods were used

for the spatial and temporal discretizations, respectively,of the heat transport Eq. [7]. We used Picard iteration
to linearize both the water flow and heat transport equations. Numerical discretization of Eq. [1] and [7] ultimately leads to a system of linear equations:Ahh ⫽ Rh

[16]

ATT ⫽ RT[17]

where Ah and AT are triagonal matrices, Rh and RT are
the known right-hand side vectors, and h and T are vectors of unknown pressure heads and temperatures at
the new time level.
The water flow and heat transport equations, as given
above, each have two unknowns: the water content ␪and the pressure head h in the modified Richards Eq.
[1], and the ice content ␪i and the temperature T in theReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

698

VADOSE ZONE J., VOL. 3, MAY 2004

heat transport Eq. [7]. We used an approach similar tothe one first used by Celia et al. (1990) to eliminate
one unknown from each equation. For the Richards
equation we first replaced the partial derivative of thewater content with respect to time with a finite difference approximation that involved the water content
␪j⫹1, k⫹1 at the new time level and the last iteration, andthe water content ␪j at the old time level, where index
j represents the time level and index k the iteration
number. We next added and subtracted the water content ␪j⫹1, k at the new time level and the previous iteration, and split this term into two parts (see also Eq. [18]below). Finally, the first difference in water contents atthe new time level and the last two iterations was replaced by the product of the hydraulic water capacity,
C, and the difference in pressure heads at the same timelevel and iterations to give:⳵␪ ⌬␪ ␪j⫹1, k⫹1 ⫺ ␪j ␪j⫹1, k⫹1 ⫺ ␪j⫹1, k ⫹ ␪j⫹1, k ⫺ ␪j
≈
⫽
⫽
⳵t ⌬t
⌬t
⌬t
⫽

␪j⫹1, k⫹1 ⫺ ␪j⫹1, k ␪j⫹1, k ⫺ ␪j⫹⌬t⌬t

⫽ C j⫹1, khj⫹1, k⫹1 ⫺ hj⫹1, k ␪j⫹1, k ⫺ ␪j⫹⌬t⌬td␪j⫹1, k[18]dh
The only unknown in Eq. [18] is now the pressure headhj⫹1, k⫹1 at the new time level and the last iteration. Allother variables are known from previous calculations
and are incorporated into the right-hand side of Eq. [16].
A similar procedure was used for the heat transport
Eq. [7] to eliminate the ice content ␪i. The partial derivative on the left side of Eq. [7] was first replaced by afinite difference approximation, which was expanded
using the ice content ␪i at the last time level and previous
iteration. Splitting this term also into two parts leads
then to⳵␪␪j⫹1, k⫹1 ⫺ ␪ ji
Lf ␳i i ⫽ Lf ␳i i
⳵t⌬t
C j⫹1, k ⫽⫽ L f ␳i

␪ ij⫹1, k⫹1 ⫺ ␪ij⫹1, k ⫹ ␪ij⫹1, k ⫺ ␪ ji⌬t

␪ ij⫹1, k⫹1 ⫺ ␪ ij⫹1, k␪ j⫹1, k ⫺ ␪ ji
⫹ L f ␳i i
[19]
⌬t
⌬t
The finite difference term of the ice content (secondterm) was next replaced by a finite difference term involving temperatures multiplied by the derivative of theice content with respect to the temperature. Finally, we
replaced the temperature in this derivative with thepressure head using the generalized Clapeyron equation, and substituted the change in water content by achange in the ice content (negative), leading to:␪ j⫹1, k⫹1 ⫺ ␪ ij⫹1, kd␪ j⫹1, k T j⫹1, k⫹1 ⫺ T j⫹1, k⫽ L f ␳i i
Lf ␳i i⌬t
dT⌬t
⫽ L f ␳i冢冣j⫹1, kT j⫹1, k⫹1 ⫺ T j⫹1, k⌬tj⫹1, kT j⫹1, k⫹1 ⫺ T j⫹1, k
⌬t
[20]

d␪
冢LgT␳ dh
冣
L␳
⫽ ⫺冢
C
gT 冣
⫽⫺

2
f i

2
f i

Again, the only unknown in Eq. [19] and [20] is thetemperature T j⫹1,k⫹1 at the new time level and the last
iteration, while all other variables are known from previous calculations and can be incorporated into the righthand side of Eq. [17].
The second term of Eq. [7] represents the change
in the amount of latent energy stored in the ice. As
demonstrated by the form of the apparent heat capacity(Eq. [12] and Fig. 1), this term can lead to significant
numerical oscillations. Let us momentarily assume asituation without flow, take a soil slightly above zerotemperature, and remove a certain amount of energy.
During the first numerical iteration at a particular time
step, the decrease in temperature of the soil is evaluatedusing the apparent heat capacity (a relatively small number) for a given above-zero temperature (Fig. 1). Thiscan lead to a considerable decrease in temperature below zero and corresponding freezing of soil water, which
can be significantly overestimated since the apparent
heat capacity increases several orders of magnitude at
the point where soil water starts freezing. When correcting for this error in the second iteration, it is possiblethat the new temperature will result in an apparent heat
capacity that will now be too low, thus returning thetemperature back to a positive value. The solution canthen easily oscillate between these two values without
converging to the right temperature. To avoid such oscillations, whenever the temperature decreases from above
to below freezing (the critical temperature), we reset
the new temperature at this critical value and restart
iterations with the maximum apparent heat capacity.
This described reset of the temperature was found to
produce relatively smooth calculations without the undesired oscillations.

APPLICATION TO A LABORATORYFREEZING EXPERIMENT
Laboratory Experiments
Mizoguchi (1990) performed freezing experiments inwhich he first packed four identical cylinders with Kanagawa sandy loam. Each cylinder was 20 cm long andhad an internal diameter of 8 cm. The samples wereprepared for the freezing test by bringing them to thesame initial state involving a uniform temperature of
6.7⬚C and a close to uniform volumetric water content
of 0.33 throughout the cylinders. Three cylinders wereused for the freezing tests, while the fourth was used toprecisely measure the initial conditions.
The sides and bottoms of three cylinders were thermally insulated, and their tops exposed to a circulating
fluid with a temperature of ⫺6⬚C. Thus, water and soil
in each cylinder was subjected to freezing from the top
down. The three cylinders were taken from the freezing
apparatus after 12, 24, and 50 h, respectively. Afterbeing removed from the cylinders, the soil samples weredivided into 1-cm-thick slices and dried in an oven to
obtain total water content (liquid water plus ice) distributions. The soil samples had a concentration of approximately 1 mg NaCl g⫺1 soil, which is sufficiently low for699

www.vadosezonejournal.org

osmotic effects to be negligible for the studied flow andtransport processes (Mizoguchi, 1990).

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Soil Hydraulic and Thermal Properties
The soil water retention curve was measured by Ishida(1983) using a combination of the following methods:a hanging water column, a pressure membrane, andequilibration over a salt solution (e.g., Jury et al., 1991).
The saturated hydraulic conductivity, Ks, was measured
directly. Van Genuchten’s (1980) analytical model withindependent m and n parameters was subsequently fitted to the retention data using the RETC model (vanGenuchten et al., 1991). From the fitted soil water retention curve and the measured Ks value, the unsaturated
hydraulic conductivity function was estimated using Eq.
[3]. The final fit resulted in the following soil hydraulicparameters: ␪r ⫽ 0.05, ␪s ⫽ 0.535, Ks ⫽ 3.2 ⫻ 10⫺6 ms⫺1,
␣ ⫽ 1.11 m⫺1, n ⫽ 1.48, m ⫽ 0.2, and l ⫽ 0.5. Figure 3
shows a very good fit of the measured retention data,
especially in the pressure head interval of interest in
this study, bounded by the two vertical dotted lines.
Simulation of the laboratory freezing experiment also
requires estimates of the thermal conductivity, both asa function of the water content for positive temperatures, and indirectly as a function of temperature forfreezing conditions. Equations [13b] and [15] were fitted
to measured thermal conductivity data (Mizoguchi,1990) as separate functions of water content and temperature. Measured negative temperatures were for thispurpose transformed into equivalent water contents using the Clapeyron Eq. [11] and the measured retention
curve. Values of the thermal conductivity close to 0⬚Cwere given extra weight in the optimization processsince most of the phase changes occur close to this temperature. Repeated nonlinear least-squares optimization
on Eq. [13b] with C5 ⫽ 4, combined with unconstrained
nonlinear optimization of Eq. [15], produced the following parameter values: C1 ⫽ 0.55 W m⫺1 K⫺1, C2 ⫽ 0.80 W
m⫺1 K⫺1, C3 ⫽ 3.07, C4 ⫽ 0.13 W m⫺1 K⫺1, C5 ⫽ 4, F1 ⫽13.05, and F2 ⫽ 1.06. It was possible to obtain excellent
fits for both Eq. [13b] and [15] when they were fitted
independently to the data. We also optimized the two
equations simultaneously, using more weight for thefrozen curve in Eq. [15], which still resulted in a good
fit of the unfrozen thermal conductivity (Fig. 2).
Heat transport was simulated using a variable heat
flux upper boundary condition:
qh ⫽ ⫺hc (TTop ⫺ TCoolant)is the heat flux (W m⫺2), h[21]

where qhc is the convective
heat transfer coefficient (W m⫺2 K⫺1), a constant representing the inverse of the surface resistance to heat exchange, and TTop and TCoolant are the temperatures at thesoil surface and of the circulating fluid, respectively.
Since there was some uncertainty as to the nature of
the coolant used in the experiment, a trial and error
procedure was used to obtain a value of 28 W m⫺2 K⫺1for hc. A zero heat-flux boundary condition was usedat the bottom of the soil column. We additionally applied zero flux boundary condition at both ends of thecolumn for the water flow calculations. As initial condi-

Fig. 3. Measured hydraulic properties of Kanagawa sandy loam (circles) and the fitted analytical model of van Genuchten (1980). The
dotted vertical lines represent pressure heads corresponding withthe highest and lowest temperatures used in the experiment, andthus specify the experimental interval of interest.

tions, we used water content and temperatures measured on the fourth experimental column. Finally, thevalue of the impedance factor, ⍀, regarded as a calibration variable, was set equal to 7. The value of ⍀ cannot
be evaluated directly when measuring soil hydraulicproperties but needs to be calibrated to the data.

Simulation Results
As mentioned in the introduction, many of the processes typical in freezing soils were already well documented some 80 yr ago. One of these is that water flows
toward freezing fronts where it changes phase from liquid to solid. This process is clearly evident in Fig. 4
where the total water content in the upper half of thecylinder increases as the column freezes. The analogy
between soil drying and freezing implies that removing

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

700

VADOSE ZONE J., VOL. 3, MAY 2004

Fig. 4. Simulated and measured values of the total volumetric water content 0, 12, 24, and 50 h after freezing started. The simulated values wereaveraged over 1-cm intervals. |␪tot|-error represents the mean absolute difference between simulated and measured water contents.

water from a soil by freezing has the same effect asdrying. Since freezing is a much quicker process thandrying, extremely high hydraulic gradients emerge andcan lead to sometimes very rapid upward flow of water.
This upward flow toward the freezing front reduces thewater content below the front, with concomitant verysignificant reductions in the unsaturated hydraulic conductivity. The freezing front is clearly visible in Fig. 4
as the depth interval where the total water content decreases rapidly. The simulation results in Fig. 4 represent
averages over 1-cm depth intervals to facilitate comparisons with the measured total water contents.
The calculated results in Fig. 4 are fairly close to themeasured values. Specifically, the rapid decrease in the
total water content at or immediately below the freezing
front and the gradual recovery deeper in the columns
are predicted well. We obtained an average absolutedifference of about 0.01 to 0.02 between simulated andmeasured total water contents, although the error varied
substantially with depth. Notice in particular that after
50 h, the calculated freezing front penetrated deeper
into the columns as compared with the experimental
data. Reasons for these differences are unclear, but
could be due to Eq. [21] providing an incomplete or
inaccurate approximation of the heat flux across theupper boundary. This heat flux should strongly affectthe position of the freezing front. The differences could
be caused also by having inaccurate estimates of theunsaturated hydraulic conductivity, which should mostly
affect the rate at which water redistributes within thesoil during freezing. Unfortunately, no measured valuesof the unsaturated hydraulic conductivity were available
to test this assumption.
The applicability of Eq. [21] in particular may needfurther study. This equation is relatively standard for
calculating heat fluxes across a surface separating a solidand a fluid, and involves a heat transfer coefficient
multiplied with the temperature difference between thesolid surface and the fluid. If, as in our case, the fluid
is pumped along the surface of the solid, the coefficientis called the convective heat transfer coefficient, hc, anddepends on the Nusselt number, which in turn depends
on the appropriate Reynolds and Prandtl numbers (e.g.,Monteith and Unsworth, 1990). Therefore, hc is usually
not a constant. A well-known example is that hc increases when the velocity of the passing fluid increases
such that a ⫺20⬚C day is far less pleasant when it is also
very windy (often referred to as wind chill). This reflects
the fact that as the wind speed increases, heat is carriedaway from a body at an accelerated rate, thus driving

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

701

Fig. 5. Simulated (symbols) and measured values (horizontal bars) of the total volumetric water content 0, 12, 24, and 50 h after freezing started.
A variable convective heat transfer coefficient, hc, was used for the first simulation, and a heat leakage bottom boundary for the second.
Simulated values were averaged over 1-cm intervals.

down body temperature. We obtained better agreement
between the simulated and measured outputs when hc
was allowed to decrease nonlinearly as a function of thesurface temperature squared, from 40 W m⫺2 K⫺1 for
0⬚C and above, to approximately 10 W m⫺2 K⫺1 for⫺4⬚C (Fig. 5). While the Nusselt number does dependon temperature, primarily because of changes in theviscosity with temperature, this effect is probably notlarge enough to motivate such a strong change in hc as
we found.
A possible reason for the simulated freezing front
penetrating deeper into the column as compared withthe experiments may have been poor insulation of thebottoms of the freezing cylinders. To test this hypothesis, and to simulate the effects of poor insulation, we
introduced the possibility of having a heat leak throughthe bottom boundary. Since no details of the leakage,if present, were available, we used Eq. [21] for this purpose. The value of TCoolant was assumed to be 6.7⬚C, equal
to the temperature of the surrounding air. Using a trialand-error procedure we tested hc values in the range
from 0.2 to 6 W m⫺1 K⫺1. These trials produced muchcloser agreement between predicted and measured total
water contents. One example is demonstrated in Fig. 5

for hc ⫽ 3 W m⫺1 K⫺1. These additional simulations
hence provide some support for the leakage hypothesis,
even if many assumptions were required to improve thematch with the data.

APPLICATION TO A SWEDISH ROAD
Problem Description
To perform a long-term test of the model, evaluate
its numerical stability, and illustrate a road application, surface temperatures gathered during nearly 1 yr
for a road in northern Sweden (the E4 highway 40 km
north of Luleå) were used as the upper temperatureboundary condition. Near-surface temperatures werecollected hourly using sensors placed about 1 cm below
the road surface.
The main differences between this and the laboratory
column experiment is that the amplitude of the upperboundary temperatures was much higher and that the soil
was subject to thawing as well as freezing. Since the surface of the road was snow-plowed during winter, the surface temperatures varied much more than would be thecase for a natural surface where temperature fluctua-

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

702

VADOSE ZONE J., VOL. 3, MAY 2004

Fig. 6. Temperature data used as the upper boundary condition, and simulated ⫺4, ⫺1, ⫺0.3, 0, 1, and 4ⴗC isotherms for a hypothetical highway
40 km north of Luleå in northern Sweden.

tions are normally dampened by the insulating snow
pack. Also, the surface temperature during summer wasmuch higher than would have been the case for a baresoil because of the ability of asphalt to absorb radiation.
We were also interested in this problem from a purely
conceptual-numerical point of view to investigate thawing behavior.
Since no soil hydraulic information was available forthe road and its subsurface at the time this manuscript
was prepared, the simulation was performed rather arbitrarily using a homogenous 3-m thick loam profile. The
following van Genuchten parameter values for loam aspresented by Carsel and Parrish (1988) were used forthis purpose: ␪r ⫽ 0.078, ␪s ⫽ 0.43, Ks ⫽ 2. 9 ⫻ 10⫺6
ms⫺1, ␣ ⫽ 3.6 m⫺1, n ⫽ 1.56, l ⫽ 0.5, and m ⫽ 0.36.
Measured near-surface temperatures were used as the
Dirichlet upper boundary condition, while the temperature at the lower end of the soil profile was kept constantat 3⬚C. A no-flow boundary condition was used for theflow calculations for both the upper and lower boundaries. We used the same thermal properties as those for
Kanagawa sandy loam in the first example.

Freezing–Thawing Simulations
Simulated temperatures below the road surface (Fig. 6)
followed distributions that are typical of relatively se-vere freezing conditions in that the freezing front moved
relatively deep into the profile and then leveled off. The
frost depth in this figure is located somewhere between
the 0 and ⫺0.3⬚C isotherms. Frost depth is usually defined as the depth where the temperature is zero, eventhough formally it should be related to the ice content.
Figure 6 shows that later in the year thawing starts
from the top, while leaving some parts of the subsurface
frozen. Notice that the soil is not completely thawed
even in early July. The progress of freezing and thawing
is also clearly visible in Fig. 7, which shows distributions
versus depth of the total water content and the temperature during the year. The water content distributions
generally show very steep gradients, with concomitant
extremely nonlinear pressure head distributions at andnear the freezing fronts. The mass- and energy-conservative iterative Picard schemes implemented in our numerical code were found to perform very well for theextreme situations of this example. We note that while
the results of this hypothetical application appear reasonable, they eventually should be compared against
measured data. We also emphasize that a highway of
the type simulated here in general has much lower watercontents in the upper part of its layered structure. Measured distributions hence may look somewhat different
for real roads.

703

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

Fig. 7. Simulated water content (top) and temperature (bottom) distributions for a hypothetical highway 40 km north of Luleå in northern Sweden.

CONCLUSIONSWe presented numerical solution schemes for coupledwater flow and heat transport problems involving aboveand subzero temperature conditions, and hence freezing–
thawing cycles. The numerical solutions performed wellfor two test problems, including one application thatconsiders rapidly varying surface temperatures. Results
for the laboratory freezing experiment were satisfactory
inasmuch as the response of the numerical solution wasvery close to the measured values of the total water content. The long-term road simulation test produced seemingly reasonable results, with no evidence of numericalinstabilities or problems with mass or heat conservation.
Still, we would like to see the model being tested eventually against a more realistic and complete field data set
with detailed information of the unsaturated soil hydraulic and thermal properties. The new expression weused for the thermal conductivity as a function of thewater content and temperature, applicable to bothfrozen and unfrozen conditions, fitted the experimental
data well and is easy to add to existing models.
ACKNOWLEDGMENTS
The first author thanks the Swedish National Road Administration and Hans Werthén Fonden for financially supporting

his visit to the George E. Brown, Jr. Salinity Laboratory inRiverside, CA, where most of the research leading to thispaper started. The authors thank Dr. Åke Hermansson of theSwedish National Road and Transport Research Institute forproviding the road surface temperature data. The authors
are grateful for the insightful comments presented by threereviewers. The work was also supported by the Swedish National Road Administration (Agreement No. AL 90 A 99:
8683) and by SAHRA (Sustainability of semi-Arid Hydrology
and Riparian Areas) under the STC Program of the National
Science Foundation (Agreement No. EAR-9876800).

REFERENCES
Alberty, R.A., and R.J. Silbey. 1992. Physical chemistry. John Wiley &Sons, Singapore.
Andersland, O.B., D.C. Wiggert, and S.H. Davies. 1995. Frozen soil
subsurface barriers: Formation and ice erosion. J. Contam. Hydrol. 23:133–147.
Beskow, G. 1935. Soil freezing and frost heaving. (In Swedish, withEnglish abstract.) SGU, Årsbok 26, Stockholm, Sweden.
Cary, J.W., and H.F. Mayland. 1972. Salt and water movement in
unsaturated frozen soil. Soil Sci. Soc. Am. Proc. 36:549–555.
Campbell, G.S. 1985. Soil physics with BASIC. Elsevier, New York.
Carsel, R.F., and R.S. Parrish. 1988. Developing joint probability
distributions of soil water retention characteristics. Water Resour.
Res. 24:755–769.
Celia, M.A., E.T. Bouloutas, and R.L. Zarba. 1990. A general massconservative numerical solution for the unsaturated flow equation.
Water Resour. Res. 26:1483–1496.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

704

VADOSE ZONE J., VOL. 3, MAY 2004

Edlefsen, N.E., and A.B.C. Andersen. 1943. Thermodynamics of soil
moisture. Hilgardia 15:31–298.
Fayer, M.J. 2000. UNSAT-H Version 3.0: Unsaturated soil waterand heat flow model. Theory, user manual, and examples. PacificNorthwest National Laboratory, Richland, WA.
Flerchinger, G.N., and K.E. Saxton. 1989. Simultaneous heat andwater model of a freezing snow-residue-soil system I. Theory anddevelopment. Trans. ASAE 32:565–571.
Flerchinger, G.N., C.L. Hanson, and J.R. Wight. 1996. Modelingevapotranspiration and surface energy budgets across a watershed.
Water Resour. Res. 32:2539–2548.
Fuchs, M., G.S. Campbell, and R.I. Papendick. 1978. An analysis ofsensible and latent heat flow in a partially frozen unsaturated soil.
Soil Sci. Soc. Am. J. 42:379–385.
Groenevelt, P.H., and B.D. Kay. 1974. On the interaction of waterand heat transport in frozen and unfrozen soils: II. The liquid
phase. Soil Sci. Soc. Am. Proc. 38:400–404.
Guymon, G.L., R.L. Berg, and T.V. Hromadka. 1993. Mathematicalmodel of frost heave and thaw settlement in pavements. CRREL
Report 93-2. CRREL, Hanover, NH.
Guymon, G.L., and J.N. Luthin. 1974. A coupled heat and moisturetransport model for arctic soils. Water Resour. Res. 10:995–1001.
Harlan, R.L. 1973. Analysis of coupled heat-fluid transport in partiallyfrozen soil. Water Resour. Res. 9:1314–1323.
Ishida, T. 1983. Water transport through soil-plant system. (In Japanese.) Ph.D. diss. University of Tokyo.
Jansson, P.-E., and S. Halldin. 1980. SOIL water and heat model:
Technical description. Swedish Coniferous Forest Proj. Tech. Rep.
26. Swedish University of Agricultural Sciences, Uppsala, Sweden.
Jansson, P.-E., and L. Karlberg. 2001. Coupled heat and mass transfer
model for soil-plant-atmosphere systems. Royal Institute of Technology, Dept of Civil and Environmental Engineering, Stockholm.
Jones, J.S. 1981. State-of-the-art report—Engineering practice in artificial ground freezing. Eng. Geol. 18:313–326.
Jury, W.A., W.R. Gardner, and W.H. Gardner. 1991. Soil physics. 5th
ed. John Wiley & Sons, New York.
Kay, B.D., and P.H. Groenevelt. 1974. On the interaction of waterand heat transport in frozen and unfrozen soils: I. Basic theory;the vapor phase. Soil Sci. Soc. Am. Proc. 38:395–400.
Kennedy, I., and B. Sharratt. 1998. Model comparisons to simulate
soil frost depth. Soil Sci. 163:636–645.
Koopmans, R.W.R., and R.D. Miller. 1966. Soil freezing and soil
water characteristic curves. Soil Sci. Soc. Am. J. 30:680–685.
Loch, J.P.G. 1978. Thermodynamic equilibrium between ice and waterin porous media. Soil Sci. 126:77–80.
Lundin, L.-C. 1989. Water and heat flows in frozen soils. Basic theory
and operational modeling. Ph.D. thesis. Uppsala University, Uppsala, Sweden.
Lundin, L.-C. 1990. Hydraulic properties in an operational model offrozen soil. J. Hydrol. (Amsterdam) 118:289–310.
Lykosov, V.N., and E.G. Palagin. 1980. Method and example of computation of heat and moisture transport in freezing soil in thepresence of a snow cover. Sov. Hydrol.: Sel. Pap. 19:69–73.
Lytton, R.L., D.E. Pufahl, C.H. Michalak, H.S. Liang, and B.J. Demp-sey. 1993. An integrated model of the climatic effects in pavements. FHWA Rep. RD-90–033. Federal Highway Administration,
McLean, VA.
Miller, R.D. 1973. Soil freezing in relation to pore water pressure andtemperature. p. 344–352. In 2nd Int. Conf. on Permafrost, Yakutsk,
Siberia, Natl. Acad. Sci., Washington.
Miller, R.D. 1978. Frost heaving in non-colloidal soils. p. 708–713. In
Proc. 3rd International Conf. on Permafrost, Edmonton, Alberta.
Natl. Res. Council of Canada, Ottawa, ON.
Miller, R.D. 1980. Freezing phenomena in soils. p. 254–299. In D. Hillel(ed.) Applications of soil physics. Academic Press, New York.
Mizoguchi, M. 1990. Water, heat and salt transport in freezing soil.
Ph.D. thesis. (In Japanese.) University of Tokyo.
Monteith, J.L., and M. Unsworth, 1990. Principles of environmental
physics. 2nd ed. Butterworth-Heinemann, Burlington, MA.
Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.
Nassar, I.N., and R. Horton. 1989. Water transport in unsaturated
nonisothermal salty soil: II. Theoretical development. Soil Sci. Soc.
Am. J. 53:1330–1337.
Nassar, I.N., and R. Horton. 1992. Simultaneous transfer of heat,water, and solute in porous media: I. Theoretical development.
Soil Sci. Soc. Am. J. 56:1350–1356.
Noborio, K., K.J. McInnes, and J.L. Heilman. 1996. Two-dimensionalmodel for water, heat and solute transport in furrow-irrigated soil:
I. Theory. Soil Sci. Soc. Am. J. 60:1001–1009.
Scanlon, B., K. Keese, R.C. Reedy, J. Šimůnek, and B. Andraski,2003. Variations in flow and transport in thick desert vadose zones
in response to paleoclimatic forcing (0– 90 kyr): Monitoring, modeling, and uncertainties, Water Resour. Res., No. 7, 1179, doi:10.1029/
2002WR001604, 13.1–13.7.
Seyfried, M.S., and M.D. Murdock. 1997. Use of air permeability toestimate infiltrability of frozen soil. J. Hydrol. (Amsterdam) 202:
95–107.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten. 1998. The HYDRUS-1D software package for simulating water flow and solutetransport in two-dimensional variably saturated media. Version 2.0.
IGWMC-TPS-70. International Ground Water Modeling Center,
Colorado School of Mines, Golden, CO.
Spaans, E.J.A., and J.M. Baker. 1996. The soil freezing characteristic:
Its measurement and similarity to the soil moisture characteristic.
Soil Sci. Soc. Am. J. 60:13–19.
van den Hurk, B.J.J.M., A. Verhoef, A.R. van den Berg, and H.A.R.
de Bruin. 1995. An intercomparison of three vegetation/soil models
for a sparse vineyard canopy. Q. J. R. Meteorol. Soc. 121:1867–1889.
van Genuchten, M.Th. 1980. A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am.
J. 44:892–898.
van Genuchten, M.Th., F.J. Leij, and S.R. Yates. 1991. The RETC
code for quantifying the hydraulic functions of unsaturated Soils.
Ver. 1.0 EPA Rep. 600/2-91/065. U.S. Salinity Lab., USDA-ARS,
Riverside, CA.
Williams, P.J., and M.W. Smith. 1989. The frozen earth—Fundamentals of geocryology. Cambridge University Press, Oxford.