Theoretical modeling framework for an unsaturated freezing soil

Available online at www.sciencedirect.com

Cold Regions Science and Technology 54 (2008) 19 – 35 www.elsevier.com/locate/coldregions

Theoretical modeling framework for an unsaturated freezing soil Ning Li a,b,⁎, Feixiong Chen c , Bin Xu a , Gunter Swoboda d b

a Xi'an University of Technology, Xi'an 710048, China State Key Laboratory of Frozen Soil Engineering, CAREERI, Chinese Academy of Sciences, Lanzhou 730000, China c Chengdu metro corporation Ltd.; Chengdu, 610031, China; Xi'an University of Technology, Xi'an 710048, China d Institute of Structural Analysis, University of Innsbruck, A-6020 Innsbruck, Austria

Received 9 July 2006; accepted 3 December 2007

Abstract Based on the authors' theoretical modeling framework of the saturated freezing soil, this paper further discusses the air phase in the unsaturated freezing soil in two ideal situations: air in the soil is linked to the atmosphere (the open porous medium) and air in the soil is isolated from the atmosphere (the sealed porous medium). The corresponding theoretical modeling framework for a multi-phase porous medium with interaction of water, heat and deformation is established in this paper. The proposed theoretical model is then extended to the general case (unsaturated half-open and half-sealed porous medium). Also, a finite element numerical solution to the heat-moisture-deformation coupling behavior of the unsaturated freezing soil is obtained. The computer software for solving this problem is also developed. The in-situ measurements for Hua Shixia highway roadbed located in Qinghai–Tibetan Plateau are introduced to compare with the numerical results obtained by the proposed model. The measurement results indicate that the numerical results for the temperature field variations with time are found to be in good correlation with the roadbed depth, just as the relative freezing deformations are correlated with road surface. Numerical modeling results, when compared to field measurements of temperature and deformation, are within 10 to 20%. The comparative analysis shows that proposed theoretical modeling framework and its numerical solution for the temperature–moisture-deformation coupling behavior is suitable and acceptable. © 2007 Elsevier B.V. All rights reserved. Keywords: Freezing soil; Unsaturated; Frost heave; Thaw settlement; Modeling framework; FEM model; Quasi-confined media

1. Introduction In the U.S. Alaska, Russian Siberia and on the Qinghai Tibetan Plateau of China, there have been many engineering problems associated with frost heave and ⁎ Corresponding author. Permanent address: Institute of Geotechnical Engineering, Xi'an University of Technology, Jin Hua South Road 5 Xi'an Shannxi 710048, China. Tel./fax: +86 29 83293863. E-mail addresses: [email protected] (N. Li), [email protected] (F. Chen), [email protected] (B. Xu), [email protected] (G. Swoboda). 0165-232X/$ - see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.coldregions.2007.12.001

thawing settlement in the construction of the highway, railway and building foundations. At the early stage, the empirical formulas were adopted to predict frost heave and thawing settlements. As a result, various freezing models for coupled heat and moisture flow and their numerical solutions were developed in the 1970s. The physical and theoretical basis for heat and mass transfer in frost models has been reviewed in many papers and technical reports (Kay and Perfect, 1988; Shen and Ladanyi, 1987; Sheng, 1994; Kujala, 1997). All these models can be classified into four different groups, i.e.,

20

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

hydrodynamic models (Harlan, 1973; Sheppard et al., 1978; Jansson and Halldin, 1979; Fukuda, 1982; Fukuda and Nakagawa, 1985; Guymon et al., 1980, 1993), rigidice models (Konrad and Morgenstern, 1980; Gilpin, 1980; O'Neill and Miller, 1985; Nixon, 1991), thermomechanical models (Duquennoi and Fremond, 1989; Fremond and Mikkola, 1991; Miao et al., 1999) and semi-empirical models, e.g., segregation potential models (Kujala, 1997). Based on force equilibrium and the principles of mass and energy conservation for a three phase (soil particle, ice and water) saturated porous medium, Li et al. (2000a,b) and Li et al. (2002) developed a theoretical modeling framework for heat-moisture-deformation coupling behaviors of freezing soil and jointed rock. In the proposed theoretical model, the different force interaction behaviors between the soil skeleton and ice particles in the frozen zone are taken into consideration. However, most frozen soils are unsaturated; and the soil thermal/mechanical properties as a whole vary greatly with the effective stresses which may change significantly with the pore air pressure inside the unsaturated soil. Particularly, the freezing or thawing situation might result in great changes with pore water pressure. The thickness of ice lenses varies from several millimeters to meters with the temperature change from − 1 °C to − 30 °C. When the ice lenses are very thick, the air in the thick ice can be considered to isolate from the atmosphere. In thawing soils, the air in the soil may be open to the atmosphere for shallow sands, while in deep clay, the air can be assumed as closed to atmosphere approximately. As far as the representative unit discussed in this paper is concerned, when the ice lenses are thin and homogeneously distributed in the representative unit, some air may be isolated from the atmosphere while other air may be open to the atmosphere, so this kind of soil can be assumed as a type of half-open and half-sealed porous medium. This paper continues our earlier research work and focuses on a four-phase unsaturated porous media consisting of soil particle, ice, water and air under freezing and thawing conditions. Two ideal situations of the porous features of freezing soil are considered: one is that the air in the soil is open to the atmosphere (e.g. air in the shallow sand), while the other is that air in the soil is isolated from the atmosphere (e.g. for the most deep clay). In reality, air that exists in the unsaturated frozen soil is between the two ideal conditions described above; and it is neither completely open to the atmosphere nor completely isolated from the atmosphere. Accordingly, it is very difficult to determine the real field conditions for the air in the soils, so we can assume

Fig. 1. A representative unit cross-section of an unsaturated, frozen soil.

that the soil condition is a mixture body between the two cases and therefore a half-opened and half-sealed frozen soil is examined here. At the same time, we assume that air flow satisfies Darcy's law. 2. Effective stress principle for the unsaturated freezing soil The stresses acting on all components of the unsaturated freezing soil unit are: soil particles with contacting stress σs, s is the superscript convenient for expressing tensor format, ice particles stress σi, pore water pressure of unfrozen water pw , and air pressure in soil voids pa. Considering a freezing soil unit with a cross-sectional area A, the actual soil particles contacting area is As, the unfrozen water contacting area is Aw, the ice particles contacting area is Ai and the air contacting area is Aa, as shown in Fig. 1. Accordingly, we obtain the following formula: A ¼ As þ Aw þ Ai þ Aa

ð1Þ

The total stress can be expressed as: rij ¼

As s Ai i Aw Aa rij þ rij þ pw dij þ pa dij A A A A

ð2Þ

where δij is Kronecker's delta = 0 for i≠ j and = 1 for i = j. In Eq. (2), the relationship between the average effective stress on the soil skeleton σij′ and the actual contacting stress on the soil particles σijs can be written as: rijV ¼

As s r A ij

ð3Þ

Usually, the actual ratio between the soil particle contacting area and the cross-sectional area is very small in comparison with other three items so that it can be neglected for simplicity. Therefore, we obtain: Aw Ai Aa þ þ ¼1 A A A

ð4Þ

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

where Aw/A is defined as the pore water pressure coefficient of the unfrozen water in the freezing soil unit, denoted as χ1; Aa/A is defined as the air pressure coefficient, denoted as χ2. And then Eq. (4) can be rewritten as: Ai ¼ 1  v1  v2 A

ð5Þ

21

For the opened porous soil, the air pressure decreases when water moves in. The equilibrium equation for the 2-D situation has the following form: ! 8 > ∂rV ∂sV ∂rix ∂sixy ∂Pw x xy > > ð  v Þ þ þ 1  v þ þ v1 þ qXx ¼ 0 > 1 2 < ∂x ∂y ∂x ∂y ∂x ! > ∂sixy ∂riy ∂sV ∂ryV ∂Pw > xy > > þ þ ð1  v1  v2 Þ þ þ v1 þ qXy ¼ 0 : ∂x ∂y ∂x ∂y ∂y

ð10Þ

Substituting Eqs. (3) and (5) into Eq. (2), the effective stress principle for unsaturated soil can be obtained as:

4. Continuity equation for unsaturated freezing soil

rij ¼ rijV þ ð1  v1  v2 Þriij þ v1 dij pw þ v2 dij pa

4.1. Continuity equation for open porous unsaturated freezing soil

ð6Þ

The precondition for the above effective stress principle of the unsaturated freezing soil is that na N 0, in which na is the volumetric air content (i.e., there is air in the soil). Liquid water in an unsaturated freezing soil, however, will migrate to another region when induced by a temperature gradient. As more water flows from other regions into the soil unit, air inside the soil unit is pressurized (in sealed porous medium) or transported away (in open porous media) and the air volume becomes less and less. If the water content or ice content increases to a critical value (85%), the effective stress principle for the saturated soil can be used.

In open porous unsaturated freezing soil, the air content in the soil decreases due to the deformation of the soil skeleton or volume expansion increases (1.09 times) when water becomes ice. Only when the air volume in the vicinity of the freezing front decreases to zero, can the freezing soil be expanded or heaved. This causes the soil particles to separate from each other. Soil particles are always in contact with each other when the air volume has not been eliminated by increasing ice volume. Accordingly, the continuity of the freezing soil includes two stages: the first is that the soil changes from the unsaturated (na N 0) to the saturated state, and the second is the development stage of the saturated state (na = 0).

3. Equilibrium equations for unsaturated freezing soil Both saturated and unsaturated freezing soils satisfy the following equilibrium equation: rij; j þ qXi ¼ 0

ð7Þ

where ρ is the equivalent density and X is volume force. Substituting Eq. (6) into Eq. (7), we can obtain:

 riVj; j þ ð1  v1  v2 Þriij; j þ v1 dij pw ; j þv2 dij pa ; j þqXi ¼ 0

ð8Þ For sealed porous unsaturated soil, the air pressure changes with the migrating water, Eq. (8) for the 2-D situation in this case can be written as: ! 8 i ∂si > > > ∂rxV þ ∂s xVy þ ð1  v1  v2 Þ ∂rx þ xy þ v1 ∂Pw þ v2 ∂Pa þ qXx ¼ 0 > < ∂x ∂y ∂x ∂y ∂x ∂x ! i i > ∂s ∂r ∂s V ∂rV ∂P ∂P > xy y xy y w a > > þ þ ð1  v1  v2 Þ þ þ v2 þ qXy ¼ 0 þ v1 : ∂x ∂y ∂x ∂y ∂y ∂y

ð9Þ

4.1.1. Unsaturated stage (na N 0) When the content of each component (soil particle, ice particle, unfrozen water and air) is defined as ns, ni, nw and na respectively, the volume change in each component can be obtained as (Fig. 2): q  qi : : vw  w nsnwi þ ns j  t va ns ¼ ns j  t qi

ð11Þ

  qw qw  qi n: i ¼ ni j  t vw þ  ni n: wi þ ni j  t va qi qi ð12Þ   q  qi : : vw þ 1  w nw nwi þ nw j  t va nw ¼ ðnw  1Þj  t qi

ð13Þ q  qi : : na ¼ na j  t m wa  w nanwi þ ðna  1Þj  t m ð14Þ qi

22

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

in a unit volume with the change velocity in content of dn wi (Fig. 2b). The third items on the right side from Eqs. (11)–(14) show the volume changes in each component due to the air moving out of a unit volume of soil with the velocity of t va (Fig. 2c). All the above equations from Eqs. (11)–(14) are the compatibility conditions for all the components of the unsaturated freezing soil with open pores. 4.1.1.1. Continuity equation for soil particles. The unit volume of soil changes from 1 at time t to 1 + Δv at time t + Δt due to the water entering or leaving the unit volume during the soil-freezing process. Mass conservation for soil particles can be expressed by the following equations: qs  1  ns ¼ qs ð1 þ DV Þ  ns2 DV ¼

ns ns2  ns 1¼ ns2 ns2

ð15Þ

ð16Þ

: ∂e DV ns2  ns 1 ns ¼ lim ¼  lim  ¼ Dt Y0 ∂t DtY0 Dt ns2 Dt ns

ð17Þ

where ρs is the soil density, ns2 is the soil particle content after change in volume, ε is volume strain. 4.1.1.2. Continuity equation for ice particles. Soil variation in the unit volume due to the freezing soil absorbing water isj  Y vw Dt  1. The variation due to absorbing air is j  t va Dt  1. The variation due to phase change from water into ice is qwqqi nwi ¼ qwqqi n: wi  Dt, i i where nwi is the water content caused by phase changes to form ice in the unit volume. Mass conservation for the ice skeleton can be expressed by the following equations: Fig. 2. Change in Volume of an unsaturated frozen soil per unit volume. (a) Change in volume due to water migration; (b) Change in volume due to ice formation upon freezing; (c) Change in volume due to the flow of air.

where dn s, dn i, dn w, dn a are the content change velocity of soil particle, ice particle, unfrozen water and air respectively; t vw is the water seepage velocity; ρw, ρi are the water density and ice density respectively; dn wi is the water content change velocity which takes place during phase change; and t va is the air moving velocity out of the soil. The first items on the right side from Eqs. (11)–(14) express the volume changes in each component due to the water migrating into, through and out of the soil volume (Fig. 2a). The second items on the right side from Eqs. (11)–(14) describe the volume changes in each component due to the phase change of the unfrozen water

  qw : qi  1  ni ¼ qi ð1 þ DV Þ  ni2  nwi  Dt qi

DV ¼

ni þ qqw n: wi  Dt i

ni2

1¼

ni þ 1:09n: wi  Dt 1 ni2

ð18Þ

ð19Þ

  : ∂e DV ns  ns2 qw nwi  Dt ¼ lim ¼  lim þ ð20Þ Dt Y0 ns2  Dt ∂t DtY0 Dt qi ns2  Dt : : ni q nwi ¼ þ w ni qi ni where ni2 is the ice particle content after change in volume.

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

4.1.1.3. Continuity equation for unfrozen water. Mass conservation for unfrozen water can be written as: : qw  1  nw ¼ qw ½ð1 þ DV Þ  nw2 þ j  t vw  Dt þ nwi  Dt  ð21Þ ð22Þ

  : ∂e DV nw  nw2 j  t vw Dt þ nwi  Dt ¼ lim ¼ lim  Dt Y0 nw2  Dt ∂t DtY0 Dt nw2  Dt : : t nw  j  vw  nwi ¼ ð23Þ nw where nw2 is the unfrozen water content after change in volume. 4.1.1.4. Continuity equation for air. tion for air can be written as:

Mass conserva-

qa  1  na ¼ qa ½ð1 þ DV Þ  na2 þ j  t va  Dt  na  j  t va  Dt DV ¼ 1 na2   ∂e DV na  na2 j  t va Dt ¼ lim ¼ lim  DtY0 na2  Dt ∂t DtY0 Dt na2  Dt : na  j  tva ¼ na

The air density in a sealed porous soil unit changes as the soil unit is compressed or expanded: na0 q na a0

∂ e qw  qi :  vw þ j  t va ¼ 0 nwi þ j  t ∂t qi

ð29Þ

where na0 is the initial air contents (at atmosphere pressure); ρa0 is the initial air density. The content variation rate of each component in the soil unit due to water migration can be written as: : vw ns ¼ ns j  t

ð30Þ

n: i ¼ ni j  t vw

ð31Þ

: nw ¼ ðnw  1Þj  t vw

ð32Þ

: na ¼ na j  t vw

ð33Þ

ð24Þ

The content variation rate due to phase change from water to ice can be expressed as:

ð25Þ

q  qi : : nsnwi ns ¼  w qi : ni ¼

ð26Þ

where na2 is the air content after change in volume. Substituting continuity equation of each component into Eqs. (11)–(14) respectively, the continuity equation for open porous unsaturated freezing soil can be written as: ð27Þ

The physical meaning of Eq. (27) is that the deformation of a soil unit volume in unit time equals to the sum of water and air moving into or out of the unit, plus the volume expansion due to the phase change from water to ice. 4.1.2. Development state of saturated soil (na = 0) If na = 0, the continuity equation for saturated freezing soil should be used: ∂e qw  qi :  vw ¼ 0 nwi þ j  t ∂t qi

4.2. Continuity equation for sealed porous unsaturated freezing soil

qa ¼

: nw  j  t vw  Dt  nwi  Dt DV ¼ 1 nw2

23

ð28Þ



n: w ¼

 qw qw  qi :  ni nwi qi qi



 qw  qi 1  nw n: wi qi

q  qi : nanwi n: a ¼  w qi

ð34Þ

ð35Þ

ð36Þ

ð37Þ

Owing to compression or expansion of the air in the : sealed, unsaturated frozen soil, air deformation is n ac Dt : (where n ac is the change in volumetric flux of air due to compression or expansion at the unit time) over the time period Δt. The content variation rates for each component are shown in Fig. 3. n: ¼ n n: ð38Þ s

s ac

: : ni ¼ ninac

ð39Þ

: : nw ¼ nwnac

ð40Þ

: : na ¼ ð1  na Þnac

ð41Þ

24

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

5. Energy conservation equations for unsaturated freezing soil In the same way, the energy conservation equation for unsaturated freezing soil can be written as: 

 ∂T t þ vs  jT ∂t  ∂e t ¼ j  ns J s  ns bs T Qs þ Q1 ∂t

ns qs Cs

Fig. 3. Soil per unit volume.

In contrast to open porous medium, the air density varies with the change in volumetric air content in a sealed porous medium. The deformation continuous conditions for sealed porous unsaturated freezing soil can thus be expressed as: q  qi : : : ns ¼ nsnac  w nsnwi þ ns j  t vw qi

ð42Þ

  :n ¼ n n: þ qw  qw  qi n n: þ n j  t vw i i ac i wi i qi qi

  qw  qi : : : nw nwi þ ðnw  1Þj  t vw nw ¼ nwnac þ 1  qi ð44Þ

ð45Þ

Based on Eqs. (42)–(45), the continuity equation for sealed porous unsaturated freezing soil can be written as: ∂e qw  qi : :  vw  nac ¼ 0 nwi þ j  t ∂t qi

 ∂T t þ vw  jT nw qw Cw ∂t  ∂Pw t jt ¼ j  nw J w  nw T vw þ Q 3 ∂T

ð48Þ



  ∂T t þ va  jT ∂t  ∂Pa t jt va þ Q4 ¼ j  na J a  na T ∂T

ð49Þ

na qa Ca ð43Þ

q  qi : : : nanwi þ na j  t vw na ¼ ð1  na Þnac  w qi

  ∂T t þ vi  jT ni qi Ci ∂t  ∂e t ¼ j  ni J i ni bi T Qi þ Q2 ∂t

ð47Þ

ð46Þ

Eq. (46) indicates that the freezing soil deformation equals to the total changes in water volume due to water movement and phase change, plus the compression volume of air in soil.

ð50Þ

where C s, C i, C w , C a are the thermal capacity of soil skeleton, ice skeleton, unfrozen water and air t t t t respectively; J s , J i , J w , J a are the heat flux of soil, ice, water and air respectively; βs, β i are the thermal expansion coefficients of soil particles and ice particles respectively; Q S, Q i are the energy of soil skeleton and ice skeleton required by per unit volume change due to dilation or shrinkage; Q1, Q2, Q3, Q4 are energy absorbed from other system of soil skeleton, ice skeleton, unfrozen water and air respectively. In Eqs. (47)–(50), the first term on the left side is the energy change induced by thermal capacity of each component; the second term by soil skeleton deformation, ice deformation, unfrozen water flow and migrating air respectively. On the right side, the first term is the heat conduction of soil skeleton, ice and unfrozen water in Eqs. (47)–(49), and the heat convection of air inside the soil unit in Eq. (50); the second term is the energy change by heat expansion and cold contraction of soil skeleton and ice in Eqs. (47)–(48); the pore water pressure of unfrozen water changes with

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

temperature in Eq. (49), and air associated with changes in pressure in Eq. (50); and the third term is the energy of soil skeleton, ice, unfrozen water and air conducted from other components of the soil unit from Eqs. (47)–(50). The energy conservation condition for energy exchange inside the soil unit is: Q1 þ Q2 þ Q3 þ Q4 ¼ Lqi

∂ nwi ∂t

ð51Þ

where L is the latent heat of freezing/thawing of water. Eqs. (47)–(51) can be rewritten as:
P P P  ∂T P P P þ C1t C1 þ C2 þ C3 v s  jT þ C 2 t vw  jT þ C 3 t va jT ∂t
P  P ∂ e P ∂Pw t Tj  vw ¼ j  kjT  bT Q  nw ∂t ∂T ∂Pa ∂ nwi Tj  t ð52Þ na va þ Lqi ∂T ∂t P

P

P

where C1 ¼ ns qs Cs þ ni qi Ci ; C2 ¼ nw qw Cw ; C3 ¼ na qa Ca ; P k ¼ ns ks þ ni ki þ nw kw þ na ka . ks , kw , ki , ka are the thermal conductivity of soil, water, ice and air resP pectively; b ¼ ns bs þ ni bi is the equivalent thermal P expansion coefficient; Q ¼ ns Qs þ ni Qi is the equivalent energy required for unit volume change due to dilation or shrinkage. 5.1. Energy conservation equation for open porous unsaturated freezing soil In the deformation process of an open porous saturated freezing soil, the air density and the pore air pressure do not change, that is ∂Pa t na T j va ¼ 0 ∂T

ð53Þ

expansion of air and changes in air density and pressure. The air migrating rate is zero. t va ¼ 0


P P P  ∂T P P þ C1 t C1 þ C2 þ C3 vs  jT þ C2 t vw  jT ∂t
P  P ∂e P ∂Pw ¼ j  kjT  bT Q  nw Tj  t vw ∂t ∂T ∂Pa ∂nwi ð56Þ Tj  t na va þ Lqi ∂T ∂t Since the air is completely sealed inside the soil unit, we have Eq. (29). 6. Governing equations for unsaturated freezing soil 6.1. Governing equations for half-sealed porous unsaturated freezing soil In the above context, we present the governing equations for both open and sealed porous, unsaturated, freezing soil. In fact, the most freezing soils are among the two extreme situations. Thus, the governing equations for half-sealed porous unsaturated freezing soil are needed to combine the two extreme cases. The governing equation for unsaturated freezing soil can be obtained from: 6.1.1. Static equilibrium equation The static equilibrium equation is the same as Eq. (8) mentioned above. 6.1.2. Continuity equation


P P P  ∂T P P P C1 þ C2 þ C3 vs  jT þ C2 t vw  jT þ C3 t va  jT þ C1 t ∂t
P  P ∂e P ∂Pw ∂nwi Tj  t vw þ Lqi ¼ j  kjT  bT Q  nw ∂T ∂t ∂t

∂e q  qi : : þjt vw þ j  t va  nac  w nwi ¼ 0 ∂t qi

5.2. Energy conservation equation for sealed porous unsaturated freezing soil For the sealed porous unsaturated freezing soil, the deformation of soil mainly results from compression or

ð55Þ

Therefore, the corresponding energy conservation equation becomes:

Eq. (52) can be rewritten as:

ð54Þ

25

ð57Þ

The physical meaning is that the unit soil deformation equals to the sum of in-flowing water, the air compression (or expansion) and the increment volume due to water changing into ice. 6.1.3. Energy conservation equation The energy conservation equation is the same as Eq. (52) mentioned above.

26

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

6.1.4. Supplementary equation

6.2. Governing equation for open porous unsaturated freezing soil

6.1.4.1. Compatibility conditions for content continuity deformation q  qi : : : v w þ ns j  t va  nsnac  w nsnwi ns ¼ n s j  t qi

i rV ij;j þ ð1  v1  v2 Þrij; j þ v1 dij pw; j þ qXi ¼ 0

ð58Þ   qw qw  qi : : : t t  ni nwi ni ¼ ni j  vw þ ni j  va  ninac þ qi qi ð59Þ : : t v w þ nw j nw ¼ ðnw 1Þj  t   va  nwnac q  qi : þ 1  w nw nwi qi vw þ ðna  1Þj  t va þ ð1  na Þn: ac n: a ¼ na j  t qw  qi :  nanwi qi

ð60Þ

ð61Þ

6.1.4.2. Seepage equation for unfrozen water (Darcy flow) kw t vw ¼  jðpw þ qw t gÞ gw

ð62Þ

6.1.4.3. Seepage equation for air (Darcy flow) ð63Þ

where ka and γa are the permeability and bulk density of water respectively. 6.1.4.4. Thermo-mechanical equation for air   Z nac t na0 þ 1n  j  v dt a ac pa0 na0 ð64Þ ¼ T T0 where nac is the compressed or expanded air volumetric content. 6.1.4.5. State equation for air density   Z nac qa0 na0 ¼ qa na0 þ  jt va dt 1  nac

ð66Þ

6.2.2. Continuity equation The continuity equation is the same as Eq. (27) mentioned above. 6.2.3. Energy conservation equation The energy conservation equation is the same as Eq. (54) mentioned above. 6.2.4. Supplementary equations The four supplementary equations are the same as : Eqs. (11)–(14). And at the same time, n a ¼ jt va þ ∂e na ∂ t which can be obtained from Eq. (27). If water does not freeze to ice and water does not migrate, we obtain ∂∂ et ¼ j  t va . Accordingly, it can be inferred that soil deformation is due to air migration in this case. 6.3. Governing equation for sealed porous unsaturated freezing soil

where kw and γw are the permeability and bulk density of water respectively, t g is the gravitational acceleration

ka t va ¼  jpa ga

6.2.1. Equilibrium equation

ð65Þ

The governing equation for sealed porous unsaturated freezing soil can be obtained in the same way: 6.3.1. Equilibrium equation The static equilibrium equation is the same as Eq. (8) mentioned above. 6.3.2. Continuity equation The continuity equation is the same as Eq. (46) mentioned above. 6.3.3. Energy conservation equation The energy conservation equation is the same as Eq. (56) mentioned above. 6.3.4. Supplementary equations 6.3.4.1. Compatibility conditions for content continuity deformation. The four supplementary equations for content continuity conditions are the same as Eqs. (42)–(45). And another equation can be: : : ∂e nac  na ¼ ∂t na

ð67Þ

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

When there is no water flow and no water freezing : : : into ice, na ¼ ð1  na Þnac , then ∂∂ et ¼ nac . This means that there is only deformation due to air pressure change inside the soil; also, the air compression velocity equals to the volume deformation velocity of the soil. 6.3.4.2. Thermo-mechanical condition pa ¼

nao T p0 na T0

ð68Þ

Kij11 ¼

 Z  ∂Ni ∂Nj ∂Ni ∂Nj þ E3 E1 dxdy ∂x ∂x ∂y ∂y v

ð71Þ

Kij12 ¼

 Z  ∂Ni ∂Nj ∂Ni ∂Nj þ E3 E2 dxdy ∂x ∂y ∂y ∂x v

ð72Þ

Z Kij13 ¼

v

∂Ni P N j dxdy ∂x

ð73Þ

∂Nj dxdy ∂x

ð74Þ

Z

6.3.4.3. Seepage for unfrozen water kw t v w ¼  j ð pw þ q w t gÞ gw

where

ð69Þ

6.3.4.4. State equation for air density. This State equation for air density is the same as Eq. (29) mentioned above.

Kij14

¼

Kij21

 Z  ∂Ni ∂Nj ∂Ni ∂Nj þ E3 ¼ E2 dxdy ∂y ∂x ∂x ∂y v

ð75Þ

 Z  ∂Ni ∂Nj ∂Ni ∂Nj þ E1 E3 dxdy ∂x ∂x ∂y ∂y v

ð76Þ

Ni v

Kij22 ¼

7. FEM format of governing equations The compressibility of unfrozen water and the effect of thermal dilation or shrinkage of soil skeleton and ice skeleton are ignored in preparing the governing equations for FEM. Also, a linear elastic constitutive model for soil is assumed. For frozen region, the governing equations are distributed in space by Weighted Residual Method, and then the equation of FEM for temperature field is distributed by Crank–Nicolson method and the equation of FEM for seepage field by a finite difference backwards method. At last, the differential equations are expressed by increment of displacement and pore pressure in stress and moisture field.

27

Z Kij23 ¼

v

∂Ni P N j dxdy ∂y

ð77Þ

∂Nj dxdy ∂y

ð78Þ

Z Kij24

¼

Ni v

Z F1 ¼

Ni qXx dxdy  Z

F2 ¼

∮ N f˜ dC

ð79Þ

∮ N f˜ dC

ð80Þ

C

v

Ni qXy dxdy 

C

v

j x

i y

pw ¼ pwf þ p˜ w

7.1. Dispersed in space The total displacements (u; v ) are divided into two part: additional displacements ua and va, which are caused by moisture migration and expansion because of water phase change; and mechanical displacements u and v, which are caused by mechanical effects and equal to the total displacement minus the additional displacements. The governing equations are distributed in space by Weighted Residual Method. The equilibrium equation can be expressed as: " #  " # P

Kij11 Kij12 Kij13 u p˜ w þ pwf P  v Kij21 Kij22 Kij23 " #

 Kij14 F1 þ Fxea ¼ þ ð1  vÞb 24 fT g ð70Þ F2 þ Fyea Kij P P



Fxea

ð81Þ

Z ½BT ½Dfeax gdxdy

¼

ð82Þ

v

n o Z ½BT ½D eay dxdy Fyea ¼

ð83Þ

v

b¼

Lqw T0

ð84Þ

Eð1AÞ E where E2 ¼ ð1þAÞEAð12AÞ, E3 ¼ 2ð1þA Þ ¼ G, E1 ¼ ð1þAÞð12AÞ ¼ E2 þ 2E3 , E is the Young's modulus; μ is Poisson's ratio; Ni is the shape function; F1 and F2 are nodal load vectors in X and Y directions respectively; Xx and Xy are volume force of X direction and Y direction; f˜ x and f˜ y can be obtained through stress boundary condition, that is,

28

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

˜ f˜ x ðtÞ ¼ rV x lx þ sxy ly , fy ðt Þ ¼ sxy lx þ rV y ly ; pwf is water pressure balance, that is, the soil moisture driving force; p˜ w is the induced excess pore water pressure caused by the moisture when the soil occurs to deform; Fxea and Fyea are additional node load vectors caused by εax and ε ay which are the additional volume strain amounts Rin X and Y directions respectively, that is fF ea g ¼ v ½BT ½Dfea gdxdy. The continuity equation can be expressed as:

where

P

  q2i ni C5 ¼ L P 1 þ qs ns

ð93Þ

P q2 ns0 C6 ¼ L Pi 2 qs ns

ð94Þ

h i : h i :

1 kf h 33 i  Kij31 fP u g  Kij32 fP vg  P Kij p˜ w þ pwf C 4 gw 1 1 Dq ∂nw h 34 i : Kij fT g ¼  P fF3 g þ P C4 C 4 qw ∂T ð85Þ

KijT

where

KijT 1

Z Kij31 ¼

P

Ni v

Z Kij32

¼ v

¼ v

P

∂T qXy  b ∂y Z

P

ð88Þ

KijT 3 ¼ rx

Ni v

N i Nj dxdy

ð89Þ

KijT 4 ¼ sxy

N i q˜ n dC

ð90Þ

FT ¼ 

P

v

∮

P

C

  Dq ni n2s C4 ¼ 1  1þ qw ns n2s0 P

∂T qXx  b ∂x

Z

¼

F3 ¼ 

¼

ð87Þ

Z Kij34

¼



∂ Ni ∂ Nj ∂ Ni ∂ Nj þ dxdy ∂x ∂x ∂y ∂y P

v



KijT 2

∂Nj dxdy Ni ∂y P

¼

ð95Þ

Z Ni Nj dxdy

ð96Þ

Ni Nj dxdy

ð97Þ

v

ð86Þ

P

Z  Kij33

∂Nj dxdy ∂x

  ∂Ni ∂Nj ∂Ni ∂Nj þ Ni dxdy ∂x ∂x ∂y ∂y

Z

v

v

∂Nj dxdy þ sxy ∂x

∂Nj dxdy þ ry ∂x

Z Ni

∂Nj dxdy ∂y

ð98Þ

Ni

∂Nj dxdy ∂y

ð99Þ

v

Z v

∮ N q˜ dC C

ð91Þ

Ni

Z

i T

ð100Þ

where q˜T is the thermal flux boundary condition. 7.2. Dispersed in time

where q˜ n is the water flux boundary; kf is the hydraulic conductivity of frozen soils. The energy conservation equation can be expressed as:  h i   P P ∂nw P P bkf : KijT fT g C3 þ C6 ½CT fT g þ k þ C 5 ∂T gw h i h i h i h i : : ¼ fFT g KijT1 þ KijT3 fu g KijT 2 þ KijT4 fvg

ð92Þ

In the Eq. (85), the equation of FEM for seepage field is distributed by finite difference backwards method and can be expressed as: h i h i

Dt kf b h 33 i nþ1 nþ1 p˜ w þ pwf  Kij31 fuT g  Kij32 fvT g P Kij C 4 gw h i h i Dt n n ¼  P fF3 g  Kij31 uT þ Kij32 uT ð101Þ C4

1 DP ∂nw h 34 i nþ1 Kij T  Tn þP C 4 Pw ∂T

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

Eq. (70) after being distributed with respect to time, can be rewritten together with Eq. (101 as: 9 8 3 8 9nþ1 > F1 þ Fxea > Kij11 Kij12 Kij13 P > u = < F þ F ea > = 22 23 7< 6 K 21 K K 2 P y ij ij 7 6 ij v ¼ 5 4 Dt Dt k > > : ; f > p˜ w þ pwf ; :  P F3 > Kij31 Kij32  P Kij33 C g C 4 w 4 h i 9 8 nþ1 14 > > ð1  vÞb Kij fT g 8 9 > > > > > > 0 h i = < < = nþ1 24 T ð 1  v Þb K f g 0  h i h i þ þ ij n 32 P n ; > > : K 31 fP h i  > > u g þ K f v g > > Dq ∂nw 34 ij ij > ; : Kij fT gnþ1  fT gn > qw ∂T

29

of governing equations for the unfrozen region can be obtained.

2

ð102Þ Then, in the Eq. (92), the FEM equation for temperature field is distributed by Crank–Nicolson method as:      P P ∂nw Dt P P kf b h T i C3 þ C6 ½CT  þ k þ C5 Kij 2 gw ∂T nþ1 Dt
nþ1 n  ¼ þ FT FT  T 2     P P ∂nw Dt P P kf b ½CT   k þ C5 fT n g þ C3 þ C6 2 gw ∂T

7.3. Design and software development On the basis of the geotechnical engineering simulation software FINAL (which was developed by Prof. G. Swoboda of University of Innsbruck of Austria and introduced by Prof. Li N), a computer program focusing on frozen soil engineering is developed to analyze the temperature–moisture-deformation-stress of an unsaturated freezing and thawing soil. The logic of program is illustrated in Fig. 4, and further information can be found in Swoboda (1994) and Chen (2001). 8. Application of the proposed model to Hua Shixia highway roadbed

i h i i h i h

h

 KijT 1 þ KijT 3 unþ1  un þ KijT 2 þ KijT 4 vnþ1  vn

ð103Þ Eqs. (102) and (103) are the governing equations of heat and mass transfer for the entire modeling domain. P P P When C 4 ¼ 1, C 5 ¼ C 6 ¼ 0, v ¼ 1, pwf ¼ 0, the FEM

Experimental results for Hua Shixia highway roadbed site are introduced to check the numerical results of the model described here. The field site is located in Qinghai–Tibetan Plateau with an altitude of over 4200 m above seal level, where the permafrosts are widely distributed. The engineering and technical personnel of the National State Key Laboratory of Frozen Soil Engineering have established a test station on the Hua Shixia highway and have carried out the measurements of the ground temperatures and

Fig. 4. The logic used in the formulation of the temperature–moisture-deformation model of an unsaturated frozen soil.

30

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

Table 1 Parameters of the highway roadbed

Water content by weight Unit weight(g/cm3) Dry unit weight(g/cm3) Heat capacity of unfrozen soil(J/m3.k) Heat capacity of frozen soil (J/m3.k) Thermal conductivity of unfrozen soil(W/m.k) Thermal conductivity of frozen soil(W/m.k) Unfrozen water content

Water transmissibility k(m/s)

nw0 A B Tf Af Bf ku

Char-bearing subclay

Peat-bearing subclay

Ballast-bearing subclay

Weathered bedrock

Sand-gravel filling

0.40 0.46 0.64 1.17 × 106 1.34 × 106 0.388 0.464 0.40 0.106700 − 0.5675 − 0.09 3.5 × 10− 9 13.438 3.5 × 10− 9

1.785 0.31 0.10 1.03 × 106 0.89 × 106 0.732 2.150 0.40 0.106700 − 0.5675 − 0.09 1.0 × 10− 9 13.438 1.0 × 10− 9

0.245 1.52 1.02 2.82 × 106 2.23 × 106 1.128 1.580 0.245 0.117509 − 0.5541 − 0.24 1.2 × 10− 9 13.438 1.2 × 10− 9

0.04 2.5 2.5 2.29 × 106 2.29 × 106 2.700 2.700 0.04 0.144099 − 0.5602 − 7.09 1.0 × 10− 13 0.0 1.0 × 10− 13

0.08 ~ 0.10 2.1 1.8 2.82 × 106 2.23 × 106 1.128 1.580 0.092 0.138692 − 0.5673 − 2.04 1.2 × 10− 9 13.438 1.2 × 10− 9

Note: The selection of water transmissibility has been considered moisture supply in actual frozen soil roadbed and other factor for suitability discount.

deformations under different environmental conditions with the times.

temperature T fitted by a curve according to experimental results by Zhu and Carbee (1984) is presented as:

8.1. Parameters of the highway roadbed

E ¼ 4  102 jT j0:636

The thermal, moisture and deformation parameters (Table 1) of the highway roadbed at the measured section are taken from the in-site measurements by Li (1999). The relationship between unfrozen water content nw and temperature can be expressed as:

nw0 ; N T zTf nw ¼ ð104Þ AðT ÞB ; N T bTf

T bTf

ð106Þ

8.2. Boundary conditions of the roadbed section for numerical model The numerical model domain of the Hua Shixia roadbed is shown in Fig. 5. The bottom boundary is fixed, with ∂T/∂y = 0.018 ( °C/m); the surface is under the atmospheric pressure and goes through seasonal

where A, B are the constants determined from the experiments; and Tf is freezing temperature. The parameters and relationship between unfrozen water and temperature are taken from Xu and Deng (1991). Using the experimental data of Horiguchi and Miller (1983), the hydraulic conductivity of water k (m/s) can be taken as:

ku ; N T zTf k¼ ð105Þ Af eBf T ; N T bTf where Af and Bf are again the experimental constants; and ku is the hydraulic conductivity of water in the unfrozen soil. Poisson's ratio of frozen and unfrozen soil can be taken as 0.3. The Young's modulus of unfrozen soil is given as 11.2 MPa (Lambe and Whitman, 1969); the relation between equivalent Young's modulus E and

Fig. 5. Numerical model node spacing for Hua Shixia highway roadbed cross-section.

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

31

Fig. 8. The sketch of key points of roadbed. Fig. 6. The measured temperature profiles centerline for the Hua Shixia highway at various times.

cycles of T. The surface temperature can be expressed as trigonometric function shown by Eq. (107). Also, there is no heat exchange at the side boundaries and water could enter the top boundary during the freezing. The road surface of Hua Shixia roadbed is covered by asphalt. According to adherent layer theory, the ground surface temperature with the times can be obtained through regression analysis. The measured surface temperature can be expressed (Zhu 1988) as:  T ¼ 4:5 þ 12:2 sin

2p 4p tþ 8640 3

 ð107Þ

8.3. Comparison of numerical results with field measurements 8.3.1. Temperature distribution The temperature variations along the roadbed depth at different times are shown in Fig. 6 and Fig. 7. Nu-

Fig. 7. The calculated temperature profiles centerline for the Hua Shixia highway at various times.

merical results are found to be in good correlation with the test data obtained. Seasonal change of ground surface temperature affects the underground temperature to a depth of 8.0 m, while the in-situ field measurements showed an impact to only 6.0 m. A ground temperature pattern for the numerical results is similar to the in-situ field measurements. The measured ground temperatures are slightly higher than those by numerical results, mainly because the regression equation for temperature of ground surface for a multi-years climate might differ from the real change of the atmosphere temperature in 1999. The key points are shown in Fig. 8 for FEM analysis of the roadbed; and Fig. 9 indicates that the differences and the correlations of the periodic change of temperature along the selected depths of roadbed. 8.3.2. Freezing soil deformation In order to contrast the numerical results with those obtained from the measurements, we define the measured deformation relative to its value on January 1 1999 as the basic value. The total displacements of the soil from the simulation beginning (January 1 of the 1st year) to April 15 of the 3rd year show the general magnitude of heaving deformation as Fig. 10 (a). The net deformation from January 1 of the 3rd year to April 15 of the 3rd year indicates the heaving deformation in this period as shown in Fig. 10(b). The total soil displacements from the simulation beginning (January 1 of the 1st year) to October 15 of the 3rd year still show the magnitude of heaving deformation, as shown in Fig. 10(c). But the net deformation from January 1 of the 3rd year to October 15 of the 3rd year indicates the magnitude of thawing deformation during the summer as indicated in Fig. 10(d). The displacement patterns also illustrate that the soils under natural ground surface heave from January to July and thaw from June to December while the roadbed under the asphalt road surface heave from January to April and thaw from May to December.

32

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

Fig. 9. Temperature comparison of key points between author's result and in-site measurement. (a) temperature curve of key point 1; (b) temperature curve of key point 4; (c) temperature curve of key point 7; (d) temperature curve of key point 14.

The roadbed deformation pattern can be divided into two phases: one is the frost heaving phase from November to next April; and the other is the thawing phase from May to November. These phenomena are not very distinct for the in-situ measurement cited in this paper. However, it is verified by the field measurements by Liu et al. (2002). In Fig. 11, the measured curves are plotted as relative values which should be minus the value from starting measurement to January 1, 1999. The deformations values of both the modeling and field results are shown in Fig. 11. The curves with hollow points represent the numerical results; the heaving and thawing deformations are basically consistent with site experiments denoted with black solid points. There are deviations between the measurements and the numerical results because the actual site environmental conditions of frozen soil are so complicated that it is very difficult to simulate them through numerical methods. The deformations of road surface are nonuniform in general. The similar results of the relative freezing deformations of road surface by both measurements and numerical results are shown in Fig. 12. The maximum

relative deformation appears at the side of roadbed instead of the center. This accords with another field measurement by Liu et al. (2002). This graph indicates the model's ability to predict the frozen soil roadbed deformation rational. 8.3.3. Stress distribution The stress distribution in the middle of December is shown in Fig. 13. The higher magnitude stresses mainly occur from 0.0 m to 2.0 m depth under the natural ground surface and from 0.0 m to 5.0 m depth under the asphalt roadbed. The principal stresses, presented in both the horizontal and vertical directions under the natural ground surface, are uniformly distributed. However, the complex temperature field and the existing thaw region in the roadbed result in the stresses not being distributed uniformly. The stresses are comparatively small in thawed region. 9. Conclusion Soil freezing involves the complex processes of a multi-phase porous media and interactions with water,

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

33

Fig. 10. The total displacement field at different time. (a) total displacements from the beginning (January 1st) to 3rd year April 15th; (b) increased displacements from 3rd year January 1st to 3rd year April 15th; (c) total displacement from beginning (January 1st) to 3rd year October 15th; (d) increased displacements from 3rd year January 1st to 3rd year October 15th.

heat and deformation. The basic theoretical modeling framework of the multi-phase porous medium is established in this paper in order to predict the complex freezing and thawing processes for cold region engineering applications mathematically. Further improvements made on the basis of experiments and in-situ tests during

construction of the Qinghai–Tibetan Railway will be presented in subsequent papers. The examples presented in this paper demonstrate that the proposed theoretical modeling framework and its numerical solution for the coupling behavior between the temperature–moisturedeformation regimes are suitable and acceptable.

Fig. 11. The relative displacement comparisons at Hua shixia highway roadbed.

Fig. 12. The roadbed surface deformation differences.

34

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35

Fig. 13. The modeled stress distribution in and adjacent to a crosssection of the Hua Shixia highway from coupled analyses at December middle (minus the initial stress).

Nomenclature σs soil particles with contacting stress, KPa σj ice particles stress, KPa pw pore water pressure of unfrozen water, KPa pa air pressure in soil voids, KPa A a cross-sectional area of a freezing soil unit, m2 As the actual soil particles contacting area, m2 Aw the unfrozen water contacting area, m2 Ai the ice particles contacting area, m2 Aa the air contacting area, m2 σij the total stress, KPa σijs σiji the actual contacting stress on the soil particles and ice particles respectively, KPa δij Kronecker's delta σV the average effective stress on the soil skeleton, ij KPa χ1 the pore water pressure coefficient of the unfrozen water in the freezing soil unit, Aw/A χ2 the air pressure coefficient, Aa/A X volume force, KN ρ the equivalent density, ρ = nsρs + niρi + nwρw + naρa, kg/m3 qs qw the density of soil, water, ice and air respecqi qa tively, kg/m3 ns ni the content of soil particle, ice particle, nw na unfrozen water and air respectively : n wi The content change velocity of water which phase change ice in volume, 1/s : : ns ni the content change velocity of soil particle, ice : : nw na particle, unfrozen water and air respectively, 1/s nwi The content of water which phase change ice in volume Y va the velocity of air moving out of the soil, m/s Y vw the velocity of water seepage, m/s

ε volume strain ns2 ni2 the content of soil particle, ice particle, nw2 na2 unfrozen water and air after change in volume respectively na0 nw0 the initial air content (equals outside air pressure) and initial unfrozen water content ρa0 the initial air density, kg/m3 Cs Ci the thermal capacity of soil skeleton, ice, Cw Ca unfrozen water and air respectively, J/kg °C βsβi the thermal expansion coefficient of soil particle and ice particle respectively, 1/°C t t Js Ji the heat flux of soil, ice, water and air t t JwJa respectively, J/m2 : the content variation rate of the air compresn ac sion or expansion in volume, 1/s t v st vi the volume change velocity of soil skeleton, ice skeleton respectively, m3/s QS Qi the energy of soil skeleton and ice skeleton required for unit volume change due to dilatancy or shrinkage, J P b the equivalent thermal expansion coefficient of P soil, b ¼ ns bs þ ni bi , 1/°C Q1 Q2 the energy absorbed from other system of soil Q3 Q4 skeleton, ice skeleton, unfrozen water and air respectively, J L latent thermal of phase change, J/m3 t g gravity acceleration, m/s2 γs γw the thermal conductivity of soil, water, ice and air respectively, W/m °C γi γa γw γi bulk density of water, ice and air respectively, KN/m3 γa kw ka the permeability of water and air respectively, m/s P Q the equivalent energy required for unit volume P change, Q ¼ ns Qs þ ni Qi , J nac the air volumetric content which is compressed or expanded na the air content in volume

Acknowledgements This project is supported by the National Natural Science Foundation of China, No. 19772036 and by the Knowledge Innovation Program of Chinese Academy of Sciences, No. KZCX1-SW-04. References Chen, F.X., 2001. The fully coupled modeling of the Thermal–moisturedeformation behavior for the saturated freezing soils. Ph.D. Thesis, Xi'an university of technology, Xi'an, China.

N. Li et al. / Cold Regions Science and Technology 54 (2008) 19–35 Duquennoi, C., Fremond, M., 1989. Modeling of thermal soil behavior. VTT Symposium 94, 2, 895–915. Fremond, M., Mikkola, M., 1991. Thermo mechanical modeling of freezing soil. Proceedings of the Sixth International Symposium on Ground Freezing. A.A. Balkema, Rotterdam, pp. 17–24. Fukuda, M., 1982. Experimental studies of coupled heat and moisture transfer in soils during freezing. Contribution No. 2528 from the Institute of the Low Temperature Science. Hokkaido University, Sapporo, Japan, pp. 35–91. S. Fukuda, M., Nakagawa, S., 1985. Numerical analysis of frost heaving based upon the coupled heat and water flow model. 4th International Symposium on Ground Freezing. A.A. Balkema, Boston, Massachusettes, pp. 109–117. 5–7 August 1985. Gilpin, R.R., 1980. A model for the prediction of ice lensing and frost heave in soils. Water Resources Research 16, 918–930. Guymon, G.L., Berg, R.L., Hromadka, T.V., 1993. Mathematical model of frost heave and thaw settlement in pavements. USA Cold Regions Research and Engineering Lab. Report. CRREL. 93–102. Guymon, G.L., Hromadka, T.V., Berg, R.L., 1980. A one-dimensional frost heave model based upon simultaneous heat and water flux. Cold Regions Science and Technology 3, 253–262. Harlan, R.L., 1973. Analysis of coupled heat–fluid transport in partially frozen soil. Water Resource Research 9, 1314–1323. Horiguchi, K.L., Miller, R.D., 1983. Hydraulic conductivity functions of frozen materials. Proc. 4th Int. Conf. on permafrost, Fairbanks, Alaska, July 1983. National Academy Press, Washington, D.C., pp. 504–508. Jansson, P.E., Halldin, S., 1979. Model for the annual water and energy flow in a layered soil. In: Halldin, S. (Ed.), Comparison of Forest and Energy Exchange Models, ISEM. Copenhagen/Elsevier, Amsterdam, pp. 145–163. Kay, B.D., Perfect, E., 1988. State of the art: heat and mass transfer in freezing soils. Proc. Of 5th International Symposium on Ground Freezing, pp. 3–21. Konrad, J.M., Morgenstern, N.R., 1980. A mechanistic theory of ice lens formation in fine grained soils. Canadian Geotechnical Journal 17, 476–486. Kujala, K., 1997. Estimation of frost heave and thaw weakening by statistical analyses and physical models. In: Knutsson (Ed.), Ground Freezing 97. Balkema, Rotterdam, pp. 31–41. Lambe, T.W., Whitman, R.D., 1969. Soil Mechanics. John Wiley & Sons, New York.

35

Li, D.Q., 1999. The analysis of roadbed stability on permafrost and its degradation of national highway in Qinghai province. Ph.D. Thesis, State Key Laboratory of Frozen Soil Engineering, CAS, China. Li, N., Chen, B., Chen, F.X., Xu, X.Z., 2000a. The coupled heatmoisture-mechanic model of the frozen soil. Cold Region Science and Technology 31 (3), 199–205. Li, N., Chen, B., Dang, F.N., 2000b. Coupled Temperature-seepagedeformation model of the jointed rock masses. Progress in Natural Science 10 (12), 911–918. Li, N., Chen, F.X., Su, B., Chen, G.D., 2002. Theoretical frame of the saturated freezing soil. Cold Regions Science and Technology 35, 72–80. Liu, Y.Z., Wu, Q.B., Zhang, J.M., Sheng, Y., 2002. Deformation of Highway Roadbed in Permafrost Regions of the Tibetan Plateau. Journal of Glaciology and Geocryology 24 (1), 14–19. Miao, T.D., Guo, L., Niu, Y.H., Zhang, Ch.Q., 1999. Modeling on coupled heat and moisture transfer in freezing soil using mixture theory. Science in China 42 (S), 9–16 Series D. Nixon, J.F., 1991. Discrete ice lens theory for frost heave in soils. Canadian Geotechnical Journal 28, 843–859. O'Nell, K., Miller, R.D., 1985. Exploration of a Rigid Ice Model of Frost Heave. Water Resources Research 21, 281–296. Sheng, D.C., 1994. Thermodynamics of freezing soils. Theory and application. Ph.D. Thesis, Lulea University of Technology, Sweden. Shen, M., Ladanyi, B., 1987. Modelling of coupled heat, moisture and stress filed in freezing soil. Cold Regions Science and Technology l4, 237–246. Sheppard, M.T., Kay, B.D., Loch, J.P.G., 1978. Development and testing of a computer model for heat and mass flow in freezing soils. Proc. 3rd Int. conf. on permafrost, pp. 76–81. Swoboda, G., 1994. FINAL user's manual Version6.6.Innsbruck. University of Innsbruck. Xu, X.Z., Deng, Y.S., 1991. The test study on moisture transfer in frozen soil. Science Press, Beijing. Zhu, Y.L., Carbee, D.L., 1984. Uniaxial compressive strength of frozen silt under constant deformation rates. Cold Regions Science and Technology 9 (1), 3–15. Zhu, L.N., 1988. Study of the adherent layer on different types of ground in permafrost regions on the Qinghai–Tibet Plateau. Journal of Glaciology and Geocryology 13 (2), 34–40.