冻土热流分析

VOL. 9, NO. 5 WATER RESOURCES RESEARCH OCTOBER 1973 Analysis of Coupled Heat-Fluid Transport in Partially Frozen Soil R. L. HARLAN Water Resources Branch, Department o] the Environment, Ottawa, Ontario, Canada An analogy can be made between the mechanisms of water transport in partially frozen soils and those in unsaturated soils. By use of this analogy a Darcian approach is applied to the analysis of coupled heat-fluid transport in porous media with freezing and thawing. With the aid of a numerical model, freezing-affected soil-water redistribution and infiltration to frozen soil are examined from a phenomenological.point of view, and the effects of soil type and initial conditions on the response of a hypothetical soil column are studied. In general, the model shows that the rate of upward redistribution of soil water to a freezing zone at the soil surface decreases from coarse-textured soils to fine-textured soils and decreases with in- crease in depth to the water table. Subsequent redistribution during melting of the frost wedge is shown to occur at a rate less than that associated with freezing. Infiltration to par- tially frozen soil is also shown to have a significant influence on soil-water redistribution and the response of the groundwater table. Conceptually, water in soils at moisture con- tents less than saturation can be visualized as occupying wedge-shaped volumes intercon- nected by thin liquid films existing between clay platelets and on particle surfaces [Kemper, 1960]. Liquid transfer between wedge-shaped volumes and therefore through soil in general must take place through these films. Because the properties of these adsorbed films do not change significantly with temperatures above 0øC, liquid transfer under a thermal gradient is assumed to be insignificant for most field and laboratory situations. The validity of this assumption for temperatures above 0øC has been confirmed by investigations using salt tracer techniques [Gurr et al., 1952; Kuzmak and $ereda, 1957]. For temperatures below 0øC, however, Dirksen and Miller [1966], among others, have shown the processes of water move- ment in porous materials to be altered consid- erably by the presence of an ice phase. In gen- eral, the freezing process induces both heat and mass transfer from warm regions to cold re- gions [Jumikis, 1958; Palmer, 1967]. Observations of the behavior of soil-water systems have shown liquid water to exist as films adsorbed on the surfaces of soil particles in equilibrium with ice at temperatures below 0øC, the normal freezing point of water. The thickness of these films has been shown to de- Copyright @ 1973 by the American Geophysical Union. pend mainly on temperature [Anderson and Hoekstra, 1965b] and except for very dry soils to be nearly independent of total water con- tent, that is, liquid water plus ice [Low et al., 1968a; Williams, 1968]. The thickness of the liquid films decreases from 50 A or more at 0øC to about 9 A at --5øC; from --5øC to liquid nitrogen temperatures, film thickness de- creases from about 9 A to 3 A [Anderson, 1968]. When soil and water at a constant vapor pressure and temperature are in equilibrium, the resultant energy state, or soil-water poten- tial, may be expressed by the Gibbs free energy •b, namely, • -- (Rk/M). In a, (1) where R is the ideal gas constant, k is the abso- lute temperature in. degrees Kelvin, M is the molecular weight of water in grams, and a, is the activity of liquid water. On the assumption that water vapor be- haves as an ideal gas, the activity a, of liquid water is equal to the relative vapor pressure of water in equilibrium with soil at the tem- perature and pressure of the system [Kijne and Taylor, 1964]. If the ice crystals forming in the soil pores have the properties of bulk ice (i.e., normal hexagonal ice crystals), the soil-water potential is fixed by the vapor pres- sure of the ice for each temperature at which ice and liquid water are in equilibrium and 1314 I-IARLAN: WATER TRANSPORT IN FROZEN SO•L equal to that of bulk ice at the temperature of the system [Low et eg., 1968a, b]. A unique relationship between soil-water potential and liquid water content does not exist if the me- dium exhibits hysteresis in the free energy water content relationship or if the concentra- tion of dissolved salts varies with time. As soil water segregates and freezes, the free energy of the residual unfrozen water is de- creased [Keune and Hoekstra, 1967]. Since the freezing point temperature is depressed further as the liquid water content is reduced, pre- sumably the greatest depressipn of the freezing point occurs in the water most firmly held and closest to the soil particles [Scho•eld, 1935; Williams, 1968]. Also, as ice crystals that do not contain significant quantities of dissolved salts are formed, salts are concentrated in residual liquid fraction, this concentration re- sulting in a further depression of the freezing point temperature [Ayers and Campbell, 1951; Williams, 1968]. Hoekstra and Chamberlain [1964] have shown the water present in liquid films ad- sorbed on the surfaces of soil particles and in clay platelets to be mobile and to be trans- ported under electrical as well as thermal gra- dients. X ray diffraction measurements have also shown changes in interlamellar spacings in clay-water gels before and after freezing [Ahlrichs and White, 1962; Norfish and Rattsell-Colom, 1962]. The expanded clay lattice tends to collapse on freezing but reexpands on thawing. Further, Anderson and Hoekstra [1965a] have shown the interlamellar water to behave in a liquidlike manner; only below about -80øC do the properties of liquid water in the interfacial regions approach those of the solid state [A•derson, 1968]. Since the liquid water content in a partially frozen soil is mainly a function of temperature, a thermal gradient in a partially frozen soil is analogous to a water content gradient in an unsaturated homogeneous soil [Hoekstra, 1967]. Provided continuity of the unfrozen water films is maintained, liquid transfer under a tempera- ture gradient should occur to the cold side at a rate dependent on the magnitude of the thermal gradient, the temperature of the system, and the surface or conducting area. If the continuity is broken, for example, through the formation 1315 of an ice lens, migration toward the cold side will be affected significantly. For a partially frozen soil therefore the tem- perature of the system places a 'physical re- straint' on the mobility of water. Since the thickness of adsorbed films decreases directly with temperature, the hydraulic conductivity should also decrease with temperature below 0øC and should be nearly independent of total water content and porosity [Hoekstra, 19.66]. It follows that the hydraulic conductivity of a partially frozen soil should be related to the energy state of the soil-water-ice system in much the same manner as the hydraulic con- ductivity of an unsaturated soil is related to the soil-water pressure head. Although the hy- draulic conductivity relationships should be similar for frozen and unfrozen materials, they are not necessarily identical. From the foregoing discussion an analogy can be drawn between the mechanisms of water transport in unsaturated soil and those in par- tially frozen soil. By use of this analogy a Daretan approach is applied from a hydrologic point of view to the mathematical description of the hydrodynamics of fluid transport in partially frozen soils. This approach is in sharp contrast, at least in the conceptual aspects, to the capillary model, widely used in studies of frost heaving phenomena. Specifi- cally, the model as presented in this paper can be used to examine freezing-affected soil-water redistribution as well as infiltration into frozen soil. Although the main concern in this paper is with soil-water redistribution within the un- frozen zone in close proximity to the freezing front, the model is equally applicable to anal- ysis of water movement within the frozen zone itself. MAThEMATiCAL FORMVLAT•ON The mathematical description of simultaneous heat and fluid transport in a partially frozen porous medium requires two systems of equa- tions expressing the interrelationships among the laws of fluid and heat flow, the equations of continuity for mass and energy, and the characteristics of the fluids and medium in- volved. Since analytic solutions to these equa- tions are unavailable except for special cases, numerical or approximate techniques must be employed. The mathematical model described 1316 HARLAN' WATER TRANSPORT IN FROZEN SOIL provides a finite difference solution to the one- dimensional coupled heat-fluid flow problem with freezing and thawing in a homogeneous rigid porous medium. For generality, the mass and heat transport equations can also be formu- lated in deforming coordinates, in, which the velocity fields are defined in relation to soil par- ticles rather than to a fixed point in space. Mass transport equation. As was considered previously, at temperatures less than 0øC the energy state of liquid water in equilibrium with ice is a function of temperature, and except for very dry conditions it is independent of tothl water content. Consequently, the ice phase can be conceptualized as behaving as a sink, or source, whereby water is added to or removed from storage on the basis of predefined tem- perature and energy balance criteria. On the assumption that the air phase and vapor transfer have negligible effects on net water transfer, the generalized mass transport equation for one-dimensional steady or un- steady flow in a saturated or partially saturated heterogeneous porous medium with freezing or thawing can be written as 0 I 0• 1 0(D/0/) OX DIK(X' T, 7')•xx - Ot + AS (2) where X, Pl, K, T, position c6ordinate, centimeters; time, minutes; density of liquid fraction, gm cm-•; volumetric liquid fraction, cm • cm-•; 'effective' hydraulic conductivity, cm min-1; temperature, øC; total head, centimeters; matrix or capillary pressure head, centimeters; change in ice per unit volume per unit time, gm cm -• min -1. The total head • is defined by [Rose, 1966] = z (3) where ½ is the Gibbs free energy or soil-water potential and combines the osmotic and matrix or capillary pressure heads, G is the pneumatic pressure head, and Z is the elevation head. For this study the effect of variation in pneu- matic and osmotic pressure heads on mass transfer is assumed to be negligible. The effective hydraulic conductivity K is de- fined as a function of position x, temperature T, and the matrix or capillary pressure head r. Implicit in this definition are the assumptions that the osmotic pressure head is more or less constant over the region of interest and negli- gible in relation to the matrix or capillary pressure head. Hea& transport equation. For the case in which convective heat transport associated with movement of the gaseous phase is small and its effect on net heat transfer is negligible, the one-dimensional steady or nonsteady convec- tion-conduction heat transport equation may be written as [Harlan, 1971] 0-• X(x, T, t)•xx • el Dl ot (4) where x, thermal conductivity, cal cm -1 øC -1 min-1; T, temperature, øC; Cl, bulk specific heat of water, cal gm -1 øC-l; v•, fluid flow velocity in x direction, cm min-•; c-•, 'apparent' volumetric specific heat, cal cm -• The apparent volumetric specific heat, in- corporating the latent heat of fusion, is defined by = .(x, T, t) - (5) where H, volumetric specific heat capacity, cal cm L0, latent heat of fusion, cal gm-1; 0,, volumetric ice fraction, cm • cm-a; o,, ice density, gm cm -a. The volumetric specific heat capacity H as de- fined here incorporates the specific heat of the soil material, the specific heat of the liquid water fraction, and the specific heat of the ice fraction. The assumptions in the development of the heat transfer equation are summarized as fol- lows' (1) The temperature of the fluid or fluids entering the medium is equal to the tem- perature along the appropriate boundary or boundaries. (2) The temperature distribution is smooth and continuous in both time and space. (3) Thermal resistance between fluids and the soil matrix is small, so that local fiuid and matrix temperatures are equal or approximately equal. (4) Within any volume element the HARLAN' WATER TRANSPORT IN FROZEN SOIL medium is homogeneous. (5) Free thermal con- vection is negligible. NUMERICAL SOLUTION When a change of phase, i.e., freezing or thawing, is involved, the heat and fluid flow fields are interdependent, and a change in direc- tion or intensity of one field effects a change in the other. As a result of this interdependency, the heat flow and fluid flow equations are coupled mathematically through the 'phase change component.' An optimization procedure therefore must be incorporated into the com- putational scheme, and a solution must be iterated between the heat. transfer equation and the mass transfer equation. To implement the solution of these equations, a fully implicit finite difference scheme is employed. For the region of interest in the time-space domain, a rectangular grid system with variable mesh spacing is used to discretize the boundary value problem. The solution in terms of the temperature, soil-water potential, ice content, and mass flux velocity distributions is advanced from known values of these distributions at t -- n -- I to unknown values at t -- n. For an interior node (j, n), where j is the spatial node and n is the time increment, the appropriate finite difference approximation for the mass transport equation (2) is At n K(•i_(1/• ) •i_1 n • o Axi + Axi+• + [C(•i n-(1/2)) _ At • Axi • AXial n--(1/2)) K(fi_(•/• At • Axi Axi • Axi.• K , At • -- 1/2)) 1 __ = C(•f ( •i •- Ax• • AXial K(•i+(1/2) (•i n-1 -- •i_l n-1 + Ax• n-- (1/2) At n K(•i_(1/2 ) ) Ax i + AXi+l AXi+l 1317 -- •i+i n-I - 2az) __ P•s (osin __ osin--1) (6) Pt where/Xx is the positional node spacing in centi- meters, /xt is the time step interval in minutes, /xZ is the change in elevation head (positive upward) in centimeters, and C -- 00•/0½ is the specific water capacity. The finite difference approximation for the heat transport equation (4) is (fl + o)T•_? + (-• - fl -, -- • n + • - c•)T, + (• - ,)T•+• --- (--f• -- (•)Ti_l n-1 -']- (Or -1- /• -']-• where • -- c•n-1)Ti n-1 + (--o• + ""y)Ti+l n-1 (7) Atnxi + ( 1/2)n- (1/2) Og = (Axi + AXi+l)(AXi+l ) A t,•hi_ ( •/•)•-(1/2) ctp• At•+(•/•) * = 2(Sx• + 5x•+•) C•p• •tnvi_(•/•) n-(1/•) •= 2(Sx• + 5Xi+l) and v is the mass flux in centimeters per min- ute. All numerical calculations were carried out on a CDC 6400 digital computer at the Computer Science Center, Department of Energy, Mines and Resources, Ottawa, Ontario. SoIL PROPERTIES In general, available data on the physical, thermal, and hydraulic transmission properties of frozen and partially frozen soils are inade- quate to permit the quantitive comparison of field or laboratow observations with the pre- dicted theoretical soil pro•e response. In most field or laboratory studies involving freezing, for example, the soil profile is described only qualitatively, if it is described at all, and essen- tia•y no data on the physical or thereal prop- erties are reposed. Consequently, for this study the storage and 1318 transmission properties of frozen and partially frozen soils are assumed to be directly analo- gous to those for unfrozen soil at a similar energy state as represented by the soil-water pressure head. That is, for a given soil the soil-water potential-unfrozen water content and the soil- water potential-hydraulic conductivity func- tional relationships are assumed to be equiva- lent to the pressure head-water content and the pressure head-hydraulic conductivity functional relationships. The results as presented therefore are indicative only of the relative magnitude of the response of the soil profile for varying initial conditions and soil types. For comparison three soil types are used in this study. These are (1) I)el Monte sand [Liakopoulos, 1965], (2) Geary silt loam [Hanks a•d Bowers, 1962], and (3) Yolo clay [Philip, 1969]. The pressure head-water con- H•tRL^•' W^TER TR^•sro•T • F•ozE• Sore tent and the pressure head-hydraulic conduc- tivity functional relationships for these soils are shown in Figure 1. Hysteresis in the functional relationship is not considered, but average rela- tionships are used in the calculations. Thermal conductivities are calculated by the method of de Vries [1952] from the physical properties of the soil matrix and water and ice contents according to the formula i --0 / i --0 where x• is the volume fraction of each soil constituent, X• is its thermal conductivity, and n is the number of soil constituents. The value of the weighting factor k, depends on the ratio of the average temperature gradient in the ith kind of particle to the average temperature gradient in the continuous medium. For calcu- Fig. 1. -400 -300 -200 -lOO o o.o (a) DEL MONTE SAND -400 I -300 -200 1• - 100 I • 0 0.10 0.20 0.30 0.0 0.010 0.020 0.(•30 - 10,000 -1,000 -100 -10 -1 0.0 - 10,000 -1,000 -100 -10 (b) GEARY SILT LOAM -10,000 -1,000 -100 -10 i i '011 0.2 0.3 0.4 0.5 (c) YOLO CLAY - 10,000 0.0' 011' 0'.2' 0'.3' 0.4 0'.5 WATER CONTENT cm 3 cm '3 -1 • 0.0 0.000001 \ \ \ 0.0001 -1 ,ooo -lOO I 0.0 0.0002 0.0004 0.0006 HYDRAULIC CONDUCTIVITY - cm min -1 i O.Ol Pressure head-water content and pressure head-hydraulic conductivity functional relationships for (a) Del Monte sand, (b) Geary silt. loam, and (c) Yolo clay. HARLAN' WATER TRANSPORT IN FROZEN SOIL lations presented here, k• was taken as 1 for all i. RESULTS Freezing-afjected soil-water redistribution. The effect of the soil type and the initial depth to a free water surface on upward redistribu- tion of soil water to a freezing front is shown in Figure 2. For each soil the initial soil-water content distribution at t -- 0 was taken as the equilibrium no-flux condition for the specified initial position of the groundwater table, that is, 100 or 200 cm below the soil surface. The total depth of the profile for all runs was taken as 500 cm. The initial temperature at all depths in the profile and the basal temperature bound- ary condition for all times were -F1.0øC. The surface temperature boundary condition was varied as is shown in Figure 2 (inset). To limit changes within the soil profile to those asso- ciated directly with freezing and subsequent thawing, the soil profile was treated as a closed system, and thus mass transfer across the soil boundaries was excluded. Although the general pattern of soil-water redistribution as affected by freezing is similar for all soils and for both initial water table positions (Figure 2), there are significant dif- 1319 ferences in both the rate and the magnitude of the response. First, the magnitude of the effect of upward soil-water migration to the freezing front on groundwater levels decreases from coarse-textured soils to fine-textured soils and with increase in depth to the water table. An increase in initial depth to the water table from 100 to 200 cm, for example, decreases the re- cession in water levels, as is indicated by the position of the free water surface from 25.5 to 5.1 cm and from 0.1 to (0.005 cm for