Freezing and Thawing Algorithm

Introduction

A soil-water-ice system can be described by simultaneously solving the water flow equation and heat flow equation based on the principles of mass and energy conservation. In the vertical direction, the two equations can be mathematically expressed as follows:

[1]

[2]

where z is the vertical depth [L], t is time [T], θw is unfrozen water content [L3/L3], θi is ice content [L3/L3], D(θw) is water diffusivity [L2/T], K(θw) is unsaturated soil hydraulic conductivity [L/T], ρw is density of water [M/L3], ρi is density of ice [M/L3], T is temperature, L is latent heat of freezing (KJ/KG), Ch is volumetric specific heat of soil (cal L-3°C-1), and λ is thermal conductivity of soil (cal L-1T-1°C-1 ).

Equation 1 is a form of the Richards equation, which is nonlinear and usually solved by numerical methods. In Equation 2, the convective heat term was not included since the heat transport by water flow is relatively small compared with that by heat conduction (Nixon, 1975; Fuchs et. al., 1978). Although soil temperature gradients may result in soil water redistribution as a result of vapor movement, the effect is usually small and neglected for practical purposes. The soil water content, however, has a significant impact on heat transfer and thus on soil temperature due to the high heat capacity of water.

DRAINMOD predicts daily water table depth using a water balance method. The method has been found comparable to numerical solutions of the flow equation (Skaggs et. al., 1991). The soil water distribution above the water table may be estimated by assuming that it is drained to equilibrium with the water table. Ignoring the redistribution of soil water caused by soil temperature gradients, we only need to solve the heat flow equation for a soil-water-ice system. The impact of ice formation on infiltration and soil permeability may be estimated by modifying soil hydraulic conductivities according to ice content. Soil temperature (T) as a function of time (t) and depth (z) can be calculated on a daily basis by solving Equation 2 using numerical methods (Hanks et. al., 1971; Campbell, 1985).

The Crank-Nicholson finite difference scheme was programmed into DRAINMOD to solve Equation 2 to determine the soil temperature as a function of depth. The profile was discretized into 20 elements with exponential increments:

[3]

where zi is the nodal depth from soil surface, and z1=0 at the soil surface. The constants, az and bz, are selected such that the bottom of the profile, z21, is a depth where annual fluctuation of temperatures are damped out and the temperature may be assumed constant. For example, a=2.5 cm and b=1.206 can describe a 5 m soil profile with nodal depth at 2.5, 5.5, 9.2, ..., 502 cm. Each element was assigned with its own soil water and thermal properties based on daily soil moisture, as calculated in the original DRAINMOD. Applying initial and boundary conditions, a 20 * 20 tri-diagonal matrix was formed and solved using the Thomas Algorithm (Wang and Anderson, 1982). The appropriate value for the upper boundary condition is the soil surface temperature. Since a long record of measured soil surface temperatures are not usually available for most applications, air temperature was used instead. Daily maximum and minimum air temperature were fitted with a sinusoidal curve in order to obtain hourly temperature input. The lower boundary was assumed as a constant soil temperature, which can be approximated as the long term average air temperature.

The volumetric heat capacity of soil was estimated as a weighted average of the heat capacity of different soil constituents (De Vries, 1963):

[4]

here x represents the volume fraction of each constituent. The subscript s, w, i and a denote soil solids, water, ice and air respectively. The volumetric heat capacity of air is relatively small and usually negligible. When freezing conditions are not present, xi = 0.

[5]

The thermal conductivity of soil is dependant on the composition of soil particles and the soil water content. The most commonly used physical model for calculating soil thermal conductivity was given by De Vries (1963) as,

where subscript s, w and a denote for solids, water, and air respectively. The summation extends over different solid soil constituents, characterized by their thermal conductivity and shape. k represents the ratio of the space average of the temperature gradient in the soil solids, ksi, (or air, ka) to that of the temperature gradient in the continuous medium. De Vries (1963) presented in detail the procedure for using Equation 5 to calculate thermal conductivity for different soils. Skaggs and Smith (1968) tested the procedure given by De Vries and added a correction factor for dry soils. Wierenga et. al. (1969) calculated thermal conductivities of soils over a range of soil moisture contents using the same procedure mentioned above, including the correction factor suggested by Skaggs and Smith (1968), and found good agreement between field measurements and calculations. For use in DRAINMOD, Equation 5 was used to calculate thermal conductivity of different soil layers outside of the model. A range of water contents was considered and the following regression equation (Equation 6) was obtained to describe the effect of soil water content on thermal conductivity . When there is ice formation in soil, the equation has an additional term (Karvonen, 1988),

[6]

[7]

Equation 6 and 7 are used in the modified DRAINMOD with coefficient a and b as model inputs.

Soil freezing and thawing processes at any point and time are indicated by soil temperature. When predicted soil temperature is below the threshold temperature for freezing conditions, the unfrozen water content (θw) in soil can be obtained from soil freezing characteristic (θw ~ T) curve, which relates unfrozen water content to below-zero temperature. The soil freezing characteristic curve can be measured on site, calculated empirically (Dillon and Anderson, 1966; Anderson and Tice, 1972), or estimated from the generalized form of Clapeyron equation as (Miller, 1980),

[8]

where Pw and Pi are soil pore water and ice pressure in kN/m2, ρw and ρi are density of water and ice, respectively (kg/m3), L is the latent heat of fusion (kJ/kg), and 273.15 is the freezing temperature of pure water in Kelvin and T is the below-zero temperature in centigrade. Assuming that Pi is equal to the atmospheric pressure, Equation 8 can be rewritten as (Karvonen, 1988),

[9]

where h is the soil water pressure (or potential) in m. Assuming that the soil water characteristic curve for unfrozen soil is valid for a partially frozen soil, the pressure head, h at a given below-zero temperature T, calculated from Equation 9 can be used to obtain the amount of unfrozen water content from the soil water characteristic curve. A tabulated form of the soil freezing characteristic curve was used as model input in the modified DRAINMOD.

The latent heat of fusion during the freezing and thawing may be incorporated into the heat transfer equation by using the apparent heat capacity, Ca, as,

[10]

where Ch is the volumetric heat capacity as defined in Equation 2, L is the latent heat of freezing [J/M], ρw is the density of water [M/L3], dθw/dT is the slope of freezing characteristic curve, and θw is the unfrozen water content in soil.

Snow Accumulation/Melt

Precipitation is read in on an hourly basis in DRAINMOD. A rain/snow dividing base temperature is proposed to identify the precipitation as rain or snow. Precipitation is considered as snowfall if the daily average air temperature is below the base temperature. According to Corps of Engineers (1956), 90 percent of precipitation is in the form of snow when the daily average air temperature is below 0 C. A subroutine was added to DRAINMOD to calculate a running total for the cumulative depth of snow. For snow cover, soil surface temperature was modified assuming an equilibrium condition between snow and the first soil increment with the following equations (Karvonen, 1988):

[11]

[12]

where Ksnow is the thermal conductivity of snow, which is a function of snow density; Ksoil is the thermal conductivity of the top soil increment, Tair is the average air temperature, Told is the soil temperature at the first node below the soil surface on the previous day, zsnow is the snow depth, Tsurf is the soil surface temperature and Δz1 is the thickness of the top soil increment. Snow density is initially an input in DRAINMOD and then updated according to amount of old snow left, new snowfall and snowmelt each day.

When the air temperature rises above the snowmelt base temperature, snowmelt is calculated using the degree-day method (Davar, 1973):

[13]

where M is the equivalent depth of water of the daily snow melt [L], Ta is mean daily temperature, Tb is snowmelt base temperature, and C is the degree day coefficient [L/C.day].

Kuz’min (1961, pp 107-110) suggested a value of Tb = +2 oC and C = 5 mm/C.day for plains. If field data on snowmelt events are available, both Tb and C can be locally determined by fitting the data.

In the modified DRAINMOD, Tb and C are constants that are model inputs. The total snow depth is calculated by tracking the amount of snow water equivalent and considering the daily densification effect. Air temperature was considered as the only driving force for snow melt. The amount of infiltration water from snowmelt is limited by soil freezing, or the soil ice content. When soil is frozen, most of the snowmelt may leave the field as surface runoff. A critical ice content was set up as model input. When the ice content exceeds this critical value, infiltration ceases and snowmelt leaves surface as runoff. The value of the critical ice content can be estimated as about 60% of field capacity (Willis et. al., 1960), where, for this purpose, field capacity may be defined as the soil water content at a pressure head of -300 cm of water. When ice is formed, soil hydraulic conductivities are modified using the following equation (Motivilov, 1978):

[14]

where θice and θw are volumetric content of ice and unfrozen water in soil, and K(θice) and K(θw) are hydraulic conductivities with and without ice.

References

Anderson, D. M. and A. R. Tice. 1972. Predicting unfrozen water contents in frozen soils from surface area measurements. Highway Research Record, No. 393, page 13-18. Frost Actions in Soils, Highway Research Board.

Campbell, G. S. 1985. Soil temperature and heat flow. Chapter 4 in Soil Physics with BASIC, Transport Models for Soil-Plant System. Elsevier Science Publishers B. V.

De Vries, D. A. 1963. Thermal properties of soils. Chapter 7 in Physics of Plant Environment Edited by W. R. Van Wijk

Dillon, H. B. and O. B. Andersland. 1966. Predicting unfrozen water contents in frozen soils. Canadian Geotech. J., Vol. III, No. 2, p53-60.

Fuchs, M., G. S. Campbell and R. I. Papendick, 1978. An analysis of sensible and latent heat flow in a partially frozen unsaturated soil. Soil Sci. Soc. Am. J. 42:379-385.

Hanks, R. J., D. D. Austin and W. T. Ondrechen. 1971. Soil temperature estimation by a numeric method. Soil Sci. Soc. Am. J., 35(5):665-667.

Karvonen, T. 1988. A model for predicting the effect of drainage on soil moisture, temperature and crop yield. phD Diss., Helsinki Univ. of Technol., Helsinki, Finland.

Kuz'min P. P. 1961. Melting of snow cover. Translated from Russian, Israel Program for Scientific Translations, Jerusalem, 1972.

Miller, R. D. 1980. Freezing phenomenon in soils. Chapter 11, Application of Soil Physics. Edited by Daniel Hillel, 1980, Academic Press, Inc..

Motovilov Y. G. 1978. Mathematical model of water infiltration into frozen soil. Soviet Hydrology, 17: 62-66.

Nixon, J. F. 1975. The role of convective heat transport in the thawing of frozen soils. Can. Geotech. J. 12:425-429.

Skaggs, R. W. and E. M. Smith. 1968. Apparent thermal conductivity of soil as related to soil porosity. Trans. ASAE, 11(4):504-507

Skaggs, R. W., T. Karvonen and H. Kandil. 1991. Predicting soil water fluxes in drained lands. Paper presented at the 1991 International Summer Meeting, Am. Soc. Of Agri. Eng., Albuquerque, N. M.

U. S. Army Corps of Engineers. 1956. Snow hydrology, Portland, Oregon; North Pacific Division, Corps of Engineers.

Wierenga, P. J., D. R. Nielsen and R. M. Hagan. 1969. Thermal properties of a soil based upon field and laboratory measurements. Soil Sci. Sco. Amer. Proc., 33: 354-360.

Willis, W. O., C. W. Carlson, J. Alessi and H. J. Haas. 1960. Depth of freezing and spring runoff as related to fall soil-moisture level. Canadian J. Soil Sci. 41:115-123.