VZJ sign up for citetrack
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (10)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Hansson, K.
Right arrow Articles by van Genuchten, M. Th.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Hansson, K.
Right arrow Articles by van Genuchten, M. Th.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Hansson, K.
Right arrow Articles by van Genuchten, M. Th.
Related Collections
Right arrow Heat Transport
Right arrow Numerical Solutions
Right arrow Coupled Flow/Transport Models
Published in Vadose Zone Journal 3:693-704 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

ORIGINAL RESEARCH

Water Flow and Heat Transport in Frozen Soil

Numerical Solution and Freeze–Thaw Applications

Klas Hansson*,a, Jirka Simunekb, Masaru Mizoguchic, Lars-Christer Lundina and Martinus Th. van Genuchtend

a Dep. of Earth Sciences, Uppsala Univ., Uppsala, Sweden
b Dep. of Environmental Science, Univ. of California, Riverside, CA
c Dep. of Biological and Environmental Engineering, Univ. of Tokyo, Japan
d George E. Brown, Jr. Salinity Laboratory, USDA-ARS, Riverside, CA
* Corresponding author (klas.hansson{at}hyd.uu.se).

Received 19 August 2003.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 THEORY
 APPLICATION TO A LABORATORY...
 Simulation Results
 APPLICATION TO A SWEDISH...
 CONCLUSIONS
 REFERENCES
 
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 zero temperature. The method is based on a mixed formulation for both water 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 evaluated by 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 the thermal 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 a hypothetical road problem involving simultaneous heat transport and water flow. The problem was simulated using measured values of the surface temperature for the duration of almost 1 yr. Since the road was snow-plowed during winter, surface temperatures varied more rapidly, and reached much lower values, than would have been the case 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.

Abbreviations: SVAT, soil–vegetation–atmosphere transfer


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
THEORY
 APPLICATION TO A LABORATORY...
 Simulation Results
 APPLICATION TO A SWEDISH...
 CONCLUSIONS
 REFERENCES
 
THE IMPORTANCE OF freezing and thawing processes in soils has long been recognized. Much attention initially focused on the problems of frost heave because of its importance in the construction and maintenance of 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 even flooding, due to increased snow melt surface runoff. In some engineering applications soil is artificially frozen to prevent the migration of pollutants (Andersland et al., 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 resources in Alaska and Canada was a major driving force for studies of freezing–thawing cycles in the 1970s and 1980s (Lundin, 1989).

Many of the processes that we try to formulate mathematically and simulate numerically today were already being studied and described some 70 yr ago. For example, Beskow (1935) studied the effects of freezing and thawing 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 in soil pores does not freeze at 0°C, but is subject to a freezing-point depression caused by interactions between water, soil particles, and solutes. Beskow also observed that soil freezing is generally quite similar to soil drying. Edlefsen and Andersen (1943) later tried to describe the mutual dependence of temperature, water content, 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 zero gauge 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 experimentally that 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 are either 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-Liquid-Solid). Miller (1978) later developed a rigid ice model that focused on ice lens formation and frost heave in 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 water and 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 assumption of this theory is that equilibrium may never be reached, or may take a very long time to establish (Spaans and Baker, 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 the Clapeyron 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 specifically designed to solve the heat and water flow equations, whether coupled or not. Harlan (1973) presented a model for coupled transport of heat and fluid in partially frozen soils, and to our knowledge was the first to give a numerical solution to the problem using finite differences. He implemented a formulation that involved the apparent heat capacity as proposed by Lukianov and Golovko (Fuchs et al., 1978). Guymon and Luthin (1974) presented a finite element solution to the same problem. Neither Harlan (1973) nor Guymon and Luthin (1974) in their analyses considered the effects of solutes, which are known to have a significant effect on 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). These models 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 the presence of a continuously deforming porous medium constituting the snow pack. Conditions in the upper soil layers of natural field soils are more or less uniquely determined 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 for roads, airstrips, or similar environments. Other, more general, 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 numerical model was evaluated by comparisons with short-term laboratory column freezing experiment data. The numerical stability of the model was further tested by means 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 detail elsewhere (Simunek et al., 1998; Scanlon et al., 2003).


   THEORY
TOP
 ABSTRACT
 INTRODUCTION
 THEORY
APPLICATION TO A LABORATORY...
 Simulation Results
 APPLICATION TO A SWEDISH...
 CONCLUSIONS
 REFERENCES
 
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):

[1]
where {theta}u is the volumetric unfrozen water content (L3 L–3) (={theta} + {theta}v), {theta} is the volumetric liquid water content (L3 L–3), {theta}v is the volumetric vapor content expressed as an equivalent water content (L3 L–3), {theta}i is the volumetric ice content (L3 L–3), t is time (T), z is the spatial coordinate positive upward (L), {rho}i is the density of ice (M L–3) (931 kg m–3), {rho}w is the density of liquid water (M L–3) (approximately 1000 kg m–3), h is the pressure head (L), T is the temperature (K), and S is a sink term (T–1) usually accounting for root water uptake.

Water flow in Eq. [1] is assumed to be caused by five different processes with corresponding hydraulic conductivities (given below in parentheses). The first three terms on the right-hand side of Eq. [1] represent liquid flows due to a pressure 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, {theta}(h) and KLh(h), respectively, and due to freezing–thawing effects that relate the ice content with the temperature, that is, {theta}i(T).

Soil Hydraulic Properties
While different functions for the unsaturated soil hydraulic properties 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:

[2]

[3]
where Se is effective saturation, {theta}r and {theta}s denote the residual and saturated water contents (L3 L–3), respectively, Ks is the saturated hydraulic conductivity (L T–1), and {alpha} (L–1), n, m, and l are empirical parameters. Further, I{varsigma} is the incomplete beta-function, {zeta} = {theta}l/me, p = m + 1/n, and q = 1 – 1/n. If m = 1 – 1/n, Eq. [3] reduces to the commonly used form:

[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):

[5]
where GwT is the gain factor, {gamma} is the surface tension of soil water (M T–2), evaluated as {gamma} = 75.6 – 0.1425T – 2.38 x 10–4T2, and {gamma}0(25°C) = 71.89 g s–2. The isothermal, Kvh, and thermal, KvT, vapor hydraulic conductivities in Eq. [1] are defined, among others, by Fayer (2000) and Scanlon et al. (2003). Since vapor transport is not a primary objective of this paper, we will not further discuss these two conductivities.

When the soil is frozen, the presence of ice (instead of liquid water) in some pores may significantly increase the resistance of the porous medium to water flow and lead to an apparent blocking effect. To account for this blocking, the hydraulic conductivity is often reduced by means of an impedance factor, {Omega} (Lundin, 1990), which in our study is multiplied by Q, the ratio of the ice content to the total (minus the residual) water content. The parameter Q accounts for the fact 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 the partially frozen soil, KfLh, as follows

[6]
which shows that even a small value of {Omega} 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):

[7]
where the first term on the left-hand side represents changes in the energy content, and the second and third terms represents changes 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 diffusion of 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 volumetric heat capacities of the solid (Cn), liquid (Cw), vapor (Cv), and ice (Ci), phases multiplied by their respective volumetric fractions {theta}:

[8]

Furthermore, in Eq. [7], L0 is the volumetric latent heat of vaporization of water (J m–3, M L–1 T–2), L0 = Lw{rho}w, Lw is the latent heat of vaporization of water (J kg–1) (=2.501 x 106 – 2369.2 T [°C]), and Lf is the latent heat of freezing (J kg–1, L2 T–2) (approximately 3.34 x 105 J kg–1).

The first two terms of Eq. [7] are often combined as follows:

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

[10]

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

[11]
where P is the pressure (Pa, M L–1 T–2) (= {rho}wgh), and Vw is the 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 of the retention curve, defined as the derivative of the water content with respect to the pressure head:

[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 orders of 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 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, Ca values 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 the actual 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].



View larger version (21K):
[in this window]
[in a new window]
  		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 and Parrish (1988).

 
According to Campbell (1985), the apparent thermal conductivity for unfrozen conditions, {lambda}({theta}) (W m–1 K–1, M L T–3 K–1), can be described by

[13a]

[13b]
where {lambda}0({theta}) is the thermal conductivity of the porous medium (solid plus water) in the absence of flow, ßt is the longitudinal 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 included in 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 our study (Fig. 2) , we modified Eq. [13b] by replacing {theta}, with {theta} + F{theta}i, where F is given by

[14]
so that

[15]



View larger version (17K):
[in this window]
[in a new window]
  		Fig. 2. Measured thermal conductivity, {lambda}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 can contribute 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 the ice content relatively high.

Numerical Implementation
The coupled water, heat, and vapor transport processes described above were incorporated in the numerical HYDRUS-1D code (Simunek et al., 1998). Equation [1], subject to appropriate initial and boundary conditions, was solved numerically using finite differences for both spatial and temporal discretization. The employed solution 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:

[16]

[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 {theta} and the pressure head h in the modified Richards Eq. [1], and the ice content {theta}i and the temperature T in the heat transport Eq. [7]. We used an approach similar to the 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 the water content with respect to time with a finite difference approximation that involved the water content {theta}j+1,k+1 at the new time level and the last iteration, and the water content {theta}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 {theta}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 at the 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 time level and iterations to give:

[18]

The only unknown in Eq. [18] is now the pressure head hj+1,k+1 at the new time level and the last iteration. All other 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 {theta}i. The partial derivative on the left side of Eq. [7] was first replaced by a finite difference approximation, which was expanded using the ice content {theta}i at the last time level and previous iteration. Splitting this term also into two parts leads then to

[19]

The finite difference term of the ice content (second term) was next replaced by a finite difference term involving temperatures multiplied by the derivative of the ice content with respect to the temperature. Finally, we replaced the temperature in this derivative with the pressure head using the generalized Clapeyron equation, and substituted the change in water content by a change in the ice content (negative), leading to:

[20]

Again, the only unknown in Eq. [19] and [20] is the temperature Tj+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 right-hand 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 a situation without flow, take a soil slightly above zero temperature, 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 evaluated using the apparent heat capacity (a relatively small number) for a given above-zero temperature (Fig. 1). This can 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 possible that the new temperature will result in an apparent heat capacity that will now be too low, thus returning the temperature back to a positive value. The solution can then 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 LABORATORY FREEZING EXPERIMENT
TOP
 ABSTRACT
 INTRODUCTION
 THEORY
 APPLICATION TO A LABORATORY...
Simulation Results
 APPLICATION TO A SWEDISH...
 CONCLUSIONS
 REFERENCES
 
Laboratory Experiments
Mizoguchi (1990) performed freezing experiments in which he first packed four identical cylinders with Kanagawa sandy loam. Each cylinder was 20 cm long and had an internal diameter of 8 cm. The samples were prepared for the freezing test by bringing them to the same 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 were used for the freezing tests, while the fourth was used to precisely 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. After being removed from the cylinders, the soil samples were divided 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 osmotic effects to be negligible for the studied flow and transport processes (Mizoguchi, 1990).

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, and equilibration over a salt solution (e.g., Jury et al., 1991). The saturated hydraulic conductivity, Ks, was measured directly. Van Genuchten's (1980) analytical model with independent m and n parameters was subsequently fitted to the retention data using the RETC model (van Genuchten 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 hydraulic parameters: {theta}r = 0.05, {theta}s = 0.535, Ks = 3.2 x 10–6 ms–1, {alpha} = 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.



View larger version (14K):
[in this window]
[in a new window]
  		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 with the highest and lowest temperatures used in the experiment, and thus specify the experimental interval of interest.

 
Simulation of the laboratory freezing experiment also requires estimates of the thermal conductivity, both as a function of the water content for positive temperatures, and indirectly as a function of temperature for freezing 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 this purpose transformed into equivalent water contents using the Clapeyron Eq. [11] and the measured retention curve. Values of the thermal conductivity close to 0°C were given extra weight in the optimization process since 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 the frozen 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:

[21]
where qh is the heat flux (W m–2), hc 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 the soil 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–1 for hc. A zero heat-flux boundary condition was used at the bottom of the soil column. We additionally applied zero flux boundary condition at both ends of the column for the water flow calculations. As initial conditions, we used water content and temperatures measured on the fourth experimental column. Finally, the value of the impedance factor, {Omega}, regarded as a calibration variable, was set equal to 7. The value of {Omega} cannot be evaluated directly when measuring soil hydraulic properties but needs to be calibrated to the data.


   Simulation Results
TOP
 ABSTRACT
 INTRODUCTION
 THEORY
 APPLICATION TO A LABORATORY...
 Simulation Results
APPLICATION TO A SWEDISH...
 CONCLUSIONS
 REFERENCES
 
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 the cylinder increases as the column freezes. The analogy between soil drying and freezing implies that removing water from a soil by freezing has the same effect as drying. Since freezing is a much quicker process than drying, extremely high hydraulic gradients emerge and can lead to sometimes very rapid upward flow of water. This upward flow toward the freezing front reduces the water content below the front, with concomitant very significant 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.



View larger version (36K):
[in this window]
[in a new window]
  		Fig. 4. Simulated and measured values of the total volumetric water content 0, 12, 24, and 50 h after freezing started. The simulated values were averaged over 1-cm intervals. |{theta}tot|-error represents the mean absolute difference between simulated and measured water contents.

 
The calculated results in Fig. 4 are fairly close to the measured 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 absolute difference of about 0.01 to 0.02 between simulated and measured 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 the upper boundary. This heat flux should strongly affect the position of the freezing front. The differences could be caused also by having inaccurate estimates of the unsaturated hydraulic conductivity, which should mostly affect the rate at which water redistributes within the soil during freezing. Unfortunately, no measured values of the unsaturated hydraulic conductivity were available to test this assumption.

The applicability of Eq. [21] in particular may need further study. This equation is relatively standard for calculating heat fluxes across a surface separating a solid and a fluid, and involves a heat transfer coefficient multiplied with the temperature difference between the solid surface and the fluid. If, as in our case, the fluid is pumped along the surface of the solid, the coefficient is called the convective heat transfer coefficient, hc, and depends 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 carried away from a body at an accelerated rate, thus driving down body temperature. We obtained better agreement between the simulated and measured outputs when hc was allowed to decrease nonlinearly as a function of the surface 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 depend on temperature, primarily because of changes in the viscosity with temperature, this effect is probably not large enough to motivate such a strong change in hc as we found.



View larger version (39K):
[in this window]
[in a new window]
  		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.

 
A possible reason for the simulated freezing front penetrating deeper into the column as compared with the experiments may have been poor insulation of the bottoms 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 through the 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 trial-and-error procedure we tested hc values in the range from 0.2 to 6 W m–1 K–1. These trials produced much closer 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 the match with the data.


   APPLICATION TO A SWEDISH ROAD
TOP
 ABSTRACT
 INTRODUCTION
 THEORY
 APPLICATION TO A LABORATORY...
 Simulation Results
 APPLICATION TO A SWEDISH...
CONCLUSIONS
 REFERENCES
 
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 Luleå) were used as the upper temperature boundary condition. Near-surface temperatures were collected 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 upper boundary 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 the case for a natural surface where temperature fluctuations are normally dampened by the insulating snow pack. Also, the surface temperature during summer was much higher than would have been the case for a bare soil 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 for the 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 as presented by Carsel and Parrish (1988) were used for this purpose: {theta}r = 0.078, {theta}s = 0.43, Ks = 2. 9 x 10–6 ms–1, {alpha} = 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 constant at 3°C. A no-flow boundary condition was used for the flow 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 severe 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, even though 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 and near the freezing fronts. The mass- and energy-conservative iterative Picard schemes implemented in our numerical code were found to perform very well for the extreme 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 water contents in the upper part of its layered structure. Measured distributions hence may look somewhat different for real roads.



View larger version (44K):
[in this window]
[in a new window]
  		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 Luleå in northern Sweden.

 


View larger version (28K):
[in this window]
[in a new window]
  		Fig. 7. Simulated water content (top) and temperature (bottom) distributions for a hypothetical highway 40 km north of Luleå in northern Sweden.

 

   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 THEORY
 APPLICATION TO A LABORATORY...
 Simulation Results
 APPLICATION TO A SWEDISH...
 CONCLUSIONS
REFERENCES
 
We presented numerical solution schemes for coupled water flow and heat transport problems involving above- and subzero temperature conditions, and hence freezing–thawing cycles. The numerical solutions performed well for two test problems, including one application that considers rapidly varying surface temperatures. Results for the laboratory freezing experiment were satisfactory inasmuch as the response of the numerical solution was very close to the measured values of the total water content. The long-term road simulation test produced seemingly reasonable results, with no evidence of numerical instabilities 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 we used for the thermal conductivity as a function of the water content and temperature, applicable to both frozen 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 Werthén Fonden for financially supporting his visit to the George E. Brown, Jr. Salinity Laboratory in Riverside, CA, where most of the research leading to this paper started. The authors thank Dr. Åke Hermansson of the Swedish National Road and Transport Research Institute for providing the road surface temperature data. The authors are grateful for the insightful comments presented by three reviewers. 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
TOP
 ABSTRACT
 INTRODUCTION
 THEORY
 APPLICATION TO A LABORATORY...
 Simulation Results
 APPLICATION TO A SWEDISH...
 CONCLUSIONS
 REFERENCES
 




This article has been cited by other articles:


Home page
 Vadose Zone JHome page
J. Simunek, M. Th. van Genuchten, and M. Sejna
Development and Applications of the HYDRUS and STANMOD Software Packages and Related Codes
Vadose Zone J., May 27, 2008; 7(2): 587 - 600.
[Abstract] [Full Text] [PDF]

Home page
 Vadose Zone JHome page
R. P. Daanen, D. Misra, and H. Epstein
Active-Layer Hydrology in Nonsorted Circle Ecosystems of the Arctic Tundra
Vadose Zone J., October 8, 2007; 6(4): 694 - 704.
[Abstract] [Full Text] [PDF]

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (10)
Right arrow Citing Articles via Google Scholar