Water and heat flow in a directionally frozen silty soil
K. Watanabe
Graduate School of Bioresources, Mie University,1577 Kurima-Machiya, Tsu 514-8507, Japan

ABSTRACT: Directional freezing experiments on silty soil were carried out. The water andheat flows were calculated using the modified version of the HYDRUS-1D code, which includes a soil freezing model. In this model, the liquid water pressure at subzero temperatures
was determined using temperatures, and the liquid (unfrozen) water content was estimated fromsoil water characteristic (retention curve) at room temperature. Unfrozen water content profiles
can be simulated when proper temperature profiles are calculated. The model can also simulate
water flow from the unfrozen region to the freezing front and the moisture profile in the unfrozen region. The calculated ice content roughly agreed with the column experiment when an impedance factor for the hydraulic conductivity was adjusted. However, water flow in the frozen
region was not obtained since the impedance factor reduced the hydraulic conductivity toomuch. Better estimation of the hydraulic conductivity of frozen soils will be needed in future.
1 INTRODUCTION
In soil, some water remains unfrozen at subzero temperatures and the amount of unfrozen water
decreases with the temperature. The relationship between the amount of unfrozen water, θl, andtemperature, T, is called the soil freezing curve (SFC). Understanding how unfrozen water flows
through frozen ground is important in investigations of water and solute redistribution (Baker &Spaans, 1997), soil microbial activity (Watanabe & Ito, 2008), mechanical stability and frost
heaving (Wettlaufer & Worster, 2006), waste disposal (McCauley et al., 2002), and climatechange in permafrost areas (Lopez, 2007). To simulate the unfrozen water flow in unsaturated
frozen soils, it is necessary to know not only how to express the hydraulic and thermal conductivities of the frozen soil but also how to determine its soil retention curve (soil water characteristics SWC; relationship between θl and the pressure head, h) and the SFC.
Williams (1964) and Koopmans & Miller (1966) measured the SFC and SWC of the same
soils under freezing and drying processes and found a unique relationship between the negative
temperature at which a given unfrozen water content occurs and the suction, corresponding to asimilar moisture content at room temperature. Harlan (1973) derived the SFC from SWC using
this relationship and analyzed the coupled heat and water flow in partially frozen soil numerically. In the Harlan’s simulation, the unsaturated hydraulic conductivity of soil at room temperature was also applied to that of frozen soil, assuming the same pore water geometry for frozen
and unfrozen soils. However, numerical simulations have suggested that this assumption overestimates water flow near the freezing front (Harlan, 1973; Taylor & Luthin, 1974; Jame &Norum, 1980). When the soil is frozen, the presence of ice in some pores may block water flow.
To account for this blocking, several impedance factors have been introduced (e.g., Jame &Norum, 1980; Lundin, 1990; Smirnova et al., 2000). However, Black & Hardenberg (1991)
criticized the use of an impedance factor, stating that it is a potent and wholly arbitrary correction function for determining the hydraulic conductivity of frozen soils. Newman & Wilson(1997) also concluded that an impedance factor is unnecessary when an accurate SWC and therelationship between hydraulic conductivity and water pressure are defined.

15

Harlan’s concept and the impedance factor have been improved in several numerical studies(e.g., Flerchinger & Saxton, 1989; Zhao et al., 1997; Stähli et al., 1999). Hansson et al. (2004)
also included these models in the HYDRUS-1D code and analyzed both laboratory and fieldsoil-freezing experiments. In this study, we performed a laboratory directional freezing experiment of unsaturated silty soil and simulated movement of water and heat in the soil using themodified HYDRUS-1D code to verify the models with the impedance factor and to estimate thethermal and hydraulic conductivity of the frozen soil.
2 SOIL FREEZING MODEL
2.1 Water and heat flow equations
Assuming that vapor and ice flows are negligible, variably saturated water flow in abovefreezing and subzero soil is described using a modified Richards’ equation as follows (e.g., Noborio et al., 1996; Hansson et al., 2004):.
∂θl (h) ρi ∂θi (T ) ∂ ⎛ ∂h∂T ⎞+= ⎜ Kh+ K h + KT
⎟
ρl ∂t
∂t
∂z ⎝
∂z
∂z ⎠

(1)

where θi is the volumetric ice content, t is time, z is the spatial coordinate, ρi is the density of ice,and Kh and KT are the hydraulic conductivities of the flow due to a pressure head gradient anddue to a temperature gradient, respectively. When ice is formed in soil pores, it releases latentheat, Lf, and the heat transport is described as follows:
Cp

∂θ
∂T
∂ ⎛ ∂T ⎞∂T− L f ρi i = ⎜ λ⎟ − Cl ql∂t∂t ∂z ⎝ ∂z ⎠∂z(2)where Cp and Cl are the volumetric heat capacity of the soil particles and liquid water, respectively, λ is the thermal conductivity, and ql is the liquid water flux. The left-hand side of equation (2) can be rewritten using the apparent volumetric heat capacity, Ca.
∂θi ⎞ ∂T ⎛∂θl ⎞ ∂T∂T⎛
⎜ C p − L f ρi ∂T ⎟ ∂t = ⎜ C p + L f ρi ∂T ⎟ ∂t = Ca ∂t
⎝
⎠
⎝
⎠

(3)

Equations (1) and (2) are a tightly coupled duet due to their mutual dependence ofn the water content, pressure head, and temperature, and can be solved when the SWC and SFC are available.

2.2 Soil water pressure in a frozen soil
When ice and liquid water coexist, a state equation of phase equilibrium, known as a generalized
form of the Clausius-Clapeyron equation (GCCE), arises.
vl

dPl
dP L f− vi i =
dT
dT T

(4)

where Pl and Pi are the liquid water and ice pressures, and vl and vi are the specific volumes ofliquid water and ice, respectively. Assuming that GCCE is also valid in frozen soil with Pl = ρghand Pi = 0, the matric potential of unfrozen water in frozen soil at an equilibrium state can be estimated from the temperature:
h=

Lf
g

ln

T
Tm

(5)

where Tm is the freezing temperature of bulk water in Kelvin.
When the soil pores illustrated in Figure 1a are filled with solute-free water and cooled below
0°C, water near the centers of the pores freezes easily, whereas water near the soil particles andat the corners among particles tends to remain in a liquid state due to the decrease in free energy
resulting from surface and capillary forces. Further lowering of the temperature induces more
ice formation, resulting in a decrease in the unfrozen water thickness with decreasing tempera-

16

ture. Williams (1964) and Koopmans & Miller (1966) regarded unfrozen water in freezing soil
as having the same geometry as water in drying unsaturated soils (Fig. 1b), and assumed thesame pressure difference between unfrozen water-ice interfaces and water-air interfaces. Under
these assumptions, frozen soil at an h corresponding to the T from equation (5) contains thesame amount of liquid water as unfrozen unsaturated soil at h; that is the SFC can be estimatedfrom the SWC. Furthermore, the slope of the SFC appearing in equation (3) is derived from theslope of SWC through GCCE:dθl dθl dh ρ gL f dθl
=
=
dT
dh dT
vwT dhLiquid water(6)

Liquid water

Ice

Soil
particleAir

Soilparticle(a)(b)

Fig. 1 Schematic illustration of liquid water geometry in soil pores: (a) freezing with the ice-liquid waterinterface in a saturated soil; (b) drying with the air-liquid water interface under room temperature.

2.3 Hydraulic and thermal properties
The water retention curve (SWC) and hydraulic conductivity, Kh, of unsaturated soil at room
temperature are sometimes expressed using a formula proposed by Mualem (1976) and van
Genuchten (1980):
Se =

(

)

θ ( h) − θ rn −m= 1+ αhθs − θrm
K h = K s Sel ⎡⎢1 − (1 − Se1/ m ) ⎤⎥
⎣
⎦

(7)

2

(8)

where Se is the effective saturation, θs and θr are the saturated and residual water content, respectively, Ks is the saturated hydraulic conductivity, and α, n, m, and l are empirical parameters.
The soil water pressure at –1°C is estimated as –12,500 cm from equation (5), indicating that a
SWC model that can express a relatively low pressure region is preferable for simulating frozen
soil. For this purpose, in this study, we use the following equation, derived by Durner (1994),
which combines two van Genuchten equations (7) using a weighting factor w:

Se = w1 ⎡1 + (α1h ) 1 ⎤
⎣
⎦
n

− m1

+ w2 ⎡1 + (α 2 h ) 2 ⎤
⎣
⎦
n

(9)

( w1Se1 + w2 Se 2 ) ⎛⎜ w1α1 ⎡⎣⎢1 − (1 − Sel1m ) ⎤⎦⎥ + w2α 2 ⎡⎣⎢1 − (1 − Sel 2m ) ⎤⎦⎥ ⎞⎟l

Kh = Ks

− m2

⎝

1

m1

( w1α1 + w2α 2 )

2

2

m2

⎠

2

(10)

The hydraulic conductivity KT for liquid water flow due to a temperature gradient is defined as(e.g., Hansson et al., 2004):⎛
1 dγ ⎞
KlT = Klh ⎜ hG⎟
γ 0 dT ⎠
⎝

(11)

17

where G is the enhanced factor, γ is the surface tension, and γ0 is the tension at 25°C. If we assume the same liquid water geometry as shown in Figure 1, it is thought that the decrease in hydraulic conductivity of frozen soil with decreasing unfrozen water is also estimated by equation(8). However, several studies have reported that the use of the unsaturated hydraulic conductivity for unfrozen soil for frozen soil overestimates the water flow in frozen soil (e.g., Jame &Norum, 1980; Lundin, 1990). Therefore, in this study, we invoke a modification of equation (8)using an impedance factor, Ω:K fh = 10−Ωθi / φ K h(12)where φ is the porosity and θi/φT is the degree of ice saturation of the soil.
The soil heat capacity, Cp, can be estimated by summing the heat capacity, C, multiplied by
the volumetric fraction, θ, of each soil element. With subscripts n, o, a , l, and i, for soil particles, soil organic matter, air, liquid water, and ice, respectively, and assuming that unfrozen water has the same heat capacity as liquid water at room temperature:C p = θ n Cn + θ o Co + θ a Cair + θl Cl + θi Ci(13)Campbell (1985) introduced the relationship between the amount of liquid water and the thermal
conductivity of soils, and Hansson et al. (2004) expanded this model to frozen soils by using theice fraction parameter, F,{

λ = C1 + C2 (θ + Fθi ) − (C1 − C4 ) exp − ⎡⎣C3 (θ + Fθi )⎤⎦F = 1 + F1θi F2

C5

}

(14)
(15)

Liquid water content, θ (cm3cm-3)

where C1…5, F1, and F2 are empirical parameters.
0.6

Hanging water
Pressure plate
Vapore pressure
Dew point
NMR
Eq .(7)
Eq. (9)

0.4

0.2

0

Fig. 2

10

2

3

10

4

10

5

10

|h | (cm)Soil water characteristics of Fujinomori silt.

3 COLUMN FREEZING EXPERIMENT
3.1 Material and Methods
The samples used in this study consisted of Fujinomori silt, which is highly susceptible to frost
and retains much liquid water, even when T < –1°C. Figure 2 shows the SWC measured using
several physical methods. Silt was mixed with water at θ = 0.4 and packed at a bulk density of
1.18 into an acrylic column with an internal diameter of 7.8 cm and a height of 35 cm. Fifteen
copper–constantan thermocouples and seven time domain reflectometry (TDR) probes were inserted into each column, and the side wall of the column was insulated. The TDR probes were
initially calibrated for the measured unfrozen water content by comparison with a pulsed nuclear
magnetic resonance (NMR) measurement. The column was allowed to settle at an ambient temperature of 2°C for 24 h to establish the initial water and temperature profiles and was then frozen from the upper end by controlling the temperature at both ends of the column (TL = –8°C

18and TH = 2°C). During the experiment, no water flux was allowed from either end, and the profiles of temperature and water content were monitored using the thermocouples and TDR
probes. A series of experiments with different durations of freezing was then performed for eachfreezing condition. At the end of the experimental series, the sample was cut into 2.5-cm sections to measure the total water content using the dry-oven method. The thermocouple and TDR
readings confirmed that each column had the same temperature and water profiles during the series of experiments.
0

28h

50h

10

6h

6h

20

30

0h

0h
Depth (cm)Depth (cm)

0

(a)

10

28h

30

-6

-4

-2

0

50h

20

0.1

2

Temperature (oC)

(b)
0.2

0.3

0.4

0.5

Liquid and total water content (cm3cm-3)

Liquid water content, θ (cm3cm-3)

Liquid water content, θ (cm3cm-3)

Fig 3. (a) Temperature and (b) moisture profiles measured in the frozen silt column (0, 6, 28, 50 h after
freezing started). The solid line and dashed line in moisture profiles represent total water and unfrozen
water contents, respectively.
0.6

0.6

(a)
Column experiment
NMR measurement

0.4

0.4

0.2

0

(b)

0.2

-6

-4

-2

0

Temperature (oC)Column experimentEq. (9)
0

10

2

10

3

4

10

5

10

|h | (cm)Fig. 4 (a) Soil freezing curve observed during the freezing experiment and measured by pulsed NMR
methods under thermal equilibrium condition. (b) Soil water characteristics of Fujinomori silt estimatedfrom (a).

3.2 Experimental results
Figure 3a shows the temperature profile of the freezing silt. When both ends of the column wereset at different temperatures, the soil near the column ends reached the required temperatures
quickly. The 0°C isotherm advanced at 1.57, 0.34, and 0.16 cm h-1 for 0–6, 6–24, and 24–48 h,
respectively, and lowering of the freezing point by approximately 0.5°C was observed. Thechanges in the temperature profiles were smaller than expected from the thermal conductivity,
implying heat flow from the side wall that prevented soil freezing. Although complete insulation
was difficult in the laboratory experiment, the differences in the temperature and the location of
the freezing front in soils at the center and near the wall of the column can be estimated within
0.5°C and to less than 1 cm, so we still regard it as directional freezing.
Figure 3b presents water profiles in silt at the same freezing times as shown in Figure 3a. The
solid line indicates total water content, θT, measured using the dry-oven method, and the dashedline indicates the unfrozen water content, θl, measured using TDR. The ice content, θi, was obtained by subtracting the unfrozen water from the total water content. The silt had a relatively

19

vertical initial θl = θT profile, having similar θl values for h < 100 cm (Fig. 2). An increase in θT,decrease in θl in the frozen area, and decrease in θl in the unfrozen area with the advancing
freezing front were observed, implying that the soil water flowed not only through the unfrozen
area but also through the frozen area.
Figure 4a shows the amount of unfrozen water measured with TDR during the freezing experiment (SFC). The amount of unfrozen water decreased sharply with temperature, althoughover 0.1 cm3 cm-3 of water remained as a liquid at –8°C. Figure 4b compares the unfrozen watercontent as a function of the corresponding pressure based on equation (5) with the fitted SWC as
shown in Figure 2. Since the frozen soil characteristics reasonably agreed with the unfrozen
SWC well, we confirmed SWC can be applied as the frozen soil characteristics in the numericalsimulation.
4 CALCULATIONS
The water and heat flows in the freezing experiment were simulated using a modified version of
the HYDRUS-1D code. The measured temperature and water content (0 h in Fig. 3) in the 35cm vertical silt column were given as the initial conditions. No water flux and a constant temperature (Ttop = –8°C and Tbottom = 2°C) were applied at both ends of the silt column for 48 h. No
solute effect was considered in this simulation. Table 1 lists the hydraulic and thermal parameters used. From the fitted SWC (Figs. 2 and 4b) and measured saturated hydraulic conductivity, lin equation (10) was estimated with comparison to data from Watanabe & Wake (2008). The
thermal conductivity of the frozen silt at different temperatures was first measured in the laboratory and the thermal parameters were estimated.
Table
1. Hydraulic and thermal parameter values for silt_______________________________________________________________________________________________________θsθr
α1 n1 α2n2 w2
Ks lθn θo C1 C2 C3 C4
F1 F2
3 -3
-1
-1
-1
3 -3
m
m
m
m
m
d
m
m
_______________________________________________________________________________________________________
0.57
0.06 0.35 3.1 0.011 1.7 0.461 0.06 -0.08 0.55 00.72 0.84 8.38 0.093 13 1_______________________________________________________________________________________________________
0

0

50h

10

20

30

(a)
-6

0h

28h
Depth (cm)Depth (cm)

6h

0h

28h
20

30
-4

-2

0

2

0.1

Temperature (oC)

6h

10

50h

(b)
0.2

0.3

0.4

0.5

Liquid and total water content (cm3cm-3)

Fig. 5 Profiles of (a) temperature and (b) moisture in a directionally frozen silt calculated by using HYDRUS-1D code. The solid line and dashed line in moisture profiles represent total water and unfrozen water contents, respectively.

The calculated freezing rate underestimated the measured rate, since the heat flow from theside wall was not negligible in the laboratory experiment. Heat outflow through the wall causes
quicker freezing than under fully insulated conditions, while inflow results in slower freezing.
The calculated temperature profiles became congruent with the measured profiles (Fig. 5a) when
the apparent thermal conductivity slightly larger than that obtained from the parameters in Table
1 was applied.
Figure 5b shows the water profile at this temperature profile. This model simulated theamounts of liquid water in both the frozen and unfrozen regions well, although the total amountof water in the frozen region was highly dependent on the impedance factor Ω (Fig. 6a). There20

was a one-to-one relationship between temperature and the pressure of the unfrozen water in thefrozen region. Therefore, the profile of the unfrozen water can be determined from the temperature profile directly, if the SFC is equivalent to the SWC. In other words, a reasonable unfrozen
water profile can be obtained when proper temperature profiles are calculated. The water flow
into and through the frozen soil appeared as a change in the total amount of water (or ice). When
Ω = 0, a huge pressure difference between the frozen and unfrozen regions induced water toflow to the freezing front, and the soil near the freezing front quickly reached ice saturation, sothat water could no longer pass through it. Decreasing the hydraulic conductivity of frozen soil
using Ω to reduce the water flow near the frozen front resulted in a decrease in the amount of ice
near the frozen front; however, Ω cannot be evaluated from any soil properties and needs to becalibrated from the ice profile data itself. Furthermore, equation (12) provides an extremely
small hydraulic conductivity for frozen soil according to Ω and the amount of ice (Fig. 6b) regardless of whether the amount of unfrozen water corresponds to the water path, and allows virtually no water flow in the frozen region (Fig. 5b, 6c, d). Further study of the hydraulic conductivity of frozen soil is needed to predict the total amount of water and water flow in frozen soil,
which is important for estimating water balance, solute redistribution, gas emission, and snow
water infiltration in cold regions.
0

0

5

2

0
Depth (cm)Depth (cm)

Ω = 20
10

20

(a)

30

θs
0.2

0.1

0.3

0.4

0.5

0

2

10

5
20

Ω = 20

30

(b)

10

12

10

10

10

Liquid and total water content (cm3cm-3)

8

6

4

10

10

Klh (cm/h)

0

0

(c)

(d)

10

Depth (cm)Depth (cm)

2

10

20

Ω=0 2
5 20
30

10

Ω=0
20

2

5

20
30

-10

-8

-6

-4

-2

0

0

h ( x104 cm)

0.01

0.02

0.03

Flux (cm/h)

Fig.6 Profiles of (a) moisture, (b) hydraulic conductivity (c) pressure head, and (d) liquid water flux in asilt column which was directionally frozen 50 h, calculating with different impedance factor Ω.

5 SUMMARY
A freezing model for unsaturated frozen soil implemented in HYDRUS-1D could simulate waterflow from the unfrozen region to the freezing front and the unfrozen water profile in the frozen
region. It is very useful to consider the detailed mechanisms of water and heat flow in unsaturated soil under freezing conditions.
Since the impedance factor for the hydraulic conductivity is given as a function of ice water,
it may underestimate the hydraulic conductivity as the ice water content increases. It would benecessary to describe the impedance factor in accordance with the unfrozen liquid water content
instead of the ice content. Although the model assumed to reach the phase equilibrium instanta-

21neously, phase transition from liquid water to ice may take time to reach the equilibrium, andthe transition rate would be proportional to the supercooling degree. If these time-dependent ice
formations cannot be ignored, it is necessary to take into account for the non-equilibrium effects. As similar to the non-equilibrium flow and transport model, dual porosity or dual permeability formulation, for instance, may be also useful for the freezing model. Furthermore, soluteconcentration increases near the frozen front because of the solute exclusion, resulting in depression of freezing point of soil water. Future verification of these issues from both experimentand numerical simulation is important to evaluate water balance, solute redistribution, and snow
melt infiltration in the frozen soil.
REFERENCES
Baker, J. M. & Spaans, E. J. A. 1997. Mechanics of meltwater movement above and within frozen soil,
Proc. Int. Sym. Physics, Chemistry, and Ecology of Seasonally Frozen Soils. CRREL Special Report,
97–10, pp. 31–36, US Army Cold Regions Research and Engineering Lab (CRREL), New Hampshire.
Black, P. B. & Hardenberg, M. J. 1991. Historical perspectives in frost heave research, the eary works of
S. Taber and G. Beskow, CRREL Special Report, 91–23.
Campbell, G.S. 1985. Soil physics with BASIC, Elsevier, New York.
Dash, J. G., Fu, H. & Wettlaufer, J. S. The premelting of ice and its environmental consequences, Rep.
Prog. Phys., 58, 115–167.
Durner, W. 1994. Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water
Resour. Res., 30(2), 211-223.
Flerchinger, G. N. & Saxton, K. E. 1989. Simultaneous heat and water model of a freezing snow-residuesoil system: I. Theory and development, Trans. ASAE, 32, 565–571.
Hansson, K., Šimůnek, J., Mizoguchi, M., Lundin, L.–C. & van Genuchten, M. Th. 2004. Water flow andheat transport in frozen soil: Numerical solution and freeze-thaw applications, Vadose Zone J., 3, 693704.
Harlan, R. L. 1973. Analysis of coupled heat-fluid transport in partially frozen soil, Water Resour. Res., 9,1314–1323.
Jame, Y. W. & Norum, D. I. 1980. Heat and mass transfer in freezing unsaturated porous media, Water
Resour. Res., 16, 811–819.
Koopmans, R. W. R., $ Miller, R. D. 1966. Soil freezing and soil water characteristic curves, Soil Sci.
Soc. Am. Proc., 30, 680–685.
Lopez, B. 2007. Cold Scapes, National Geographic, 212, 137–151.
Lundin, L.-C. 1990. Hydraulic properties in an operational model of frozen soil, J. Hydrolo., 118, 289–
310.
McCauley, C. A., White, D. M., Lilly, M. R. & Nyman, D. M. 2002. A comparison of hydraulic conductivities, permeabilities and infiltration rates in frozen and unfrozen soils, Cold Regions Sci. Tech., 34,
117–125.
Newman, G. P. & Wilson, G. W. 2004. Heat and mass transfer in unsaturated soils during freezing, Can.
Geotech. J., 34, 63–70.
Smirnova, T. G., Brown, J. M., Benjamin, S. G. & Kim, D. 2000. Parameterization of cold-season processes in the MAPS land-surface scheme, J. Geophys. Res., 105, 4077–4086.
Stähli, M., Jansson, P.–E. & Lundin, L–C. 1999. Soil moisture redistribution and infiltration in frozen
sandy soils, Water Resour. Res., 35, 95–103.
Taylor, G. S. & Luthin, J. N. 1978. A model for coupled heat and moisture transfere during soil freezing.
Can Geotech. J., 15, 548–555.
Watanabe, K & Ito, M. 2008. In situ observation of the distribution and activity of microorganisms in frozen soil, Cold Regions Sci. Tech, doi:10.1016/j.coldregions.2007.12.004.
Wettlaufer, J. S., & Worster, M. G. 2006. Premelting dynamics, Annu. Rev. Fluid. Mech., 38, 427–452.
Williams, P. 1964. Unfrozen water content of frozen soil and soil moisture suction, Geotechnique, 14,
231–246.
Zhao, L., Gray, D. M. & Male, D. H. 1997. Numerical analysis of simultaneous heat and mass transfer
during infiltration into frozen ground, J. Hydrol., 200, 345–363.

22