Dimensional analysis of numerical solution of the rigid ice model frost heave with hints on how to implement the solution in Matlab

Niklas Grip

Dimensional Analysis of Numerical Solution of the Rigid Ice Model Frost Heave

with Hints on How to Implement the Solution in Matlab

MASTER'S THESIS

Civilingenjörsprogrammet

1997:147 • ISSN: 1402-1617 • ISRN: LTU-EX--97/147--SE

the Rigid Ice Model of Frost Heave with Hints on How to Implement

the Solution in Matlab

Niklas Grip¹

Department of Mathematics Lulea University of Technology

SE-971 87 Lulea, Sweden e-mail: grip@sm.luth.se

May 15, 1997

1Supported by a grant from Coldtech.

Technical note: Typeset with L^aTEX2"and some accompanying packages (B^ibTEX, MakeIndex, ^AMS-L^aTEX version 1.2 and^AMSFonts), running on a Sun Sparcstation 4.

Each year, frost action in the ground causes very costly problems in vari- ous constructions like highways, airport pavements, railways, pipelines and building foundations in the urban cold regions of the world. This thesis contains a survey of the rigid ice model, which is the most complete model of frost heave from a microscopic and physical point of view. The physi- cal foundations of the model are explained and used in a derivation of the mathematical model.

In contrast to earlier texts on this subject, we have translated the math- ematical model to dimensionless variables and documented in much more detail how to solve it using a nite element method. (The reason for this is that the model is more accessible if not everyone has to \reinvent the wheel".) The resulting equations are presented in full detail and some of the computations are done in Maple, to simplify future changes.

We also outline how to write a Matlab program for implementation of the rigid ice model. [12] and [10] give some further advice on computational details needed in a complete program.

Frostskador i marken orsakar varje ar stora skador pa bland annat vagar, ygplatser, jarnvagar, rorledningar och husgrunder runt om i varlden. Detta

ar ett examensarbete om den modell for tjallyftning som idag ar mest fullstandig bade ur mikroskopisk och fysikalisk synvinkel: stelismodellen. De fysikaliska grunderna for modellen forklaras i en harledning av den matem- atiska modellen.

Till skillnad fran tidigare verk i detta amne har vi oversatt modellen till dimensionslosa variabler och visar mer detaljerat hur modellen kan losas med en nit element metod, och de ekvationer detta leder till redovisas in i minsta detalj. (Tanken med detta ar att amnet blir mer lattillgangligt om inte alla standigt maste "ateruppnna hjulet".) Nagra nodvandiga berakningar har gjorts i Maple for att forenkla mojliga framtida modieringar.

Vi visar aven i stora drag hur stelismodellen kan implementeras i Matlab.

Ytterligare nagra viktiga raknemassiga detaljer som kravs i ett fullstandigt program tas upp i [12] och [10].

Contents

Preface viii

1 Introduction 1

2 How to read this text 1

2.1 What you should read : : : : : : : : : : : : : : : : : : : : : : 1 2.1.1 Physical background : : : : : : : : : : : : : : : : : : : 1 2.1.2 The mathematical problem : : : : : : : : : : : : : : : 1 2.1.3 Numerical calculations : : : : : : : : : : : : : : : : : : 1 2.2 Sections that can be ignored: : : : : : : : : : : : : : : : : : : 1 2.2.1 Change to dimensionless variables : : : : : : : : : : : 1 2.2.2 Material properties : : : : : : : : : : : : : : : : : : : : 2 2.2.3 The appendices : : : : : : : : : : : : : : : : : : : : : : 2

3 Physical background 2

3.1 Frost damage : : : : : : : : : : : : : : : : : : : : : : : : : : : 2 3.2 Frost heave : : : : : : : : : : : : : : : : : : : : : : : : : : : : 2 3.3 Derivation of the rigid ice model : : : : : : : : : : : : : : : : 3 3.3.1 Restrictions : : : : : : : : : : : : : : : : : : : : : : : : 5 3.3.2 Water ow : : : : : : : : : : : : : : : : : : : : : : : : 5 3.3.3 Heat conservation : : : : : : : : : : : : : : : : : : : : 5 3.3.4 Phase equilibria and transport equations: : : : : : : : 7 3.3.5 Rate of heave : : : : : : : : : : : : : : : : : : : : : : : 9 3.3.6 Ice lens initiation : : : : : : : : : : : : : : : : : : : : : 11 3.3.7 Boundary conditions : : : : : : : : : : : : : : : : : : : 12 3.4 Restatement of the complete model : : : : : : : : : : : : : : : 12

4 Change of variables 15

4.1 The dimensionless rigid ice model: : : : : : : : : : : : : : : : 16

5 Material properties 17

5.1 Ice content I as a function of iw : : : : : : : : : : : : : : : : 17 5.1.1 Simplifying the integral for  : : : : : : : : : : : : : : 19 5.2 Thermal conductivity  : : : : : : : : : : : : : : : : : : : : : 20 5.3 Hydraulic conductivity k of frozen soils : : : : : : : : : : : : 21 5.4 Volumetric heat capacity C of soil : : : : : : : : : : : : : : : 21

6 The mathematical problem 21

7 Finite element method 23

7.1 Calculations at some xed point in time : : : : : : : : : : : : 23 7.1.1 Choosing trial functions

_i : : : : : : : : : : : : : : : 25 7.2 Boundary conditions : : : : : : : : : : : : : : : : : : : : : : : 27

7.3 Estimation of^P : : : : : : : : : : : : : : : : : : : : : : : : : 27 7.4 The course of a simulation : : : : : : : : : : : : : : : : : : : : 27 7.4.1 Initial computations : : : : : : : : : : : : : : : : : : : 28 7.4.2 First iterative step : : : : : : : : : : : : : : : : : : : : 28 7.4.3 Second iterative step : : : : : : : : : : : : : : : : : : : 28 7.4.4 Further modications below ^x_f : : : : : : : : : : : : : 28 7.5 Things we didn't have time to do : : : : : : : : : : : : : : : : 28 7.5.1 Initial values of^x_b and ^P : : : : : : : : : : : : : : : : 29 7.5.2 Not all boundary conditions are used: : : : : : : : : : 29 7.5.3 The numerical mesh : : : : : : : : : : : : : : : : : : : 29

A Dimensional analysis 30

A.1 Some necessary terminology : : : : : : : : : : : : : : : : : : : 30 A.2 Dimensional analysis of the rigid ice model : : : : : : : : : : 30 A.2.1 The dimensionless rigid ice model: : : : : : : : : : : : 33

B The Buckingham pi theorem 36

C Computation of

[k;i]

,

[k;i]

and

[i]

38

C.1 The worksheet : : : : : : : : : : : : : : : : : : : : : : : : : : 38 C.1.1 Some important notes about this worksheet : : : : : : 38 C.1.2 Computations: : : : : : : : : : : : : : : : : : : : : : : 38 C.1.3 Results : : : : : : : : : : : : : : : : : : : : : : : : : : 40 C.2 A slightly better approximation ofi;i : : : : : : : : : : : : : 41

D Boundary conditions 42

D.1 Discretized boundary conditions : : : : : : : : : : : : : : : : 42

E Basic structure of a Matlab simulation 43

E.1 Main program: : : : : : : : : : : : : : : : : : : : : : : : : : : 43 E.2 Initial conditions : : : : : : : : : : : : : : : : : : : : : : : : : 44 E.3 Computing current values of variables : : : : : : : : : : : : : 45 E.4 Integrand of  : : : : : : : : : : : : : : : : : : : : : : : : : : 49 E.5 Shifting the numerical mesh : : : : : : : : : : : : : : : : : : : 49

References 50

Index 52

List of Figures

1 A schematic diagram of the frozen fringe : : : : : : : : : : : : 4 2 A schematic view of soil, ice lenses and frozen fringe : : : : : 10 3 Figure 1 repeated : : : : : : : : : : : : : : : : : : : : : : : : : 13 4 Partition of x-axis and choice of trial functions

i : : : : : : 26

5 The block matrices of C andD : : : : : : : : : : : : : : : : : 46

List of Tables

1 All variables and notations used during the development of the rigid ice model : : : : : : : : : : : : : : : : : : : : : : : : 6 2 All variables and notations used in the nal statement of the

rigid ice model : : : : : : : : : : : : : : : : : : : : : : : : : : 14 3 Soil properties for some dierent types of soil : : : : : : : : : 15

Preface

Texts like this one never get nished, but sooner or later everyone has to reach a point where they decide that \OK, now it's good enough to publish".

I feel that I have reached that point now, even if there of course are more things that would be fun and interesting to do, but I leave that to others and go on with postgraduate studies in applied mathematics.

Due to other projects and courses, I have never worked full-time on this thesis, so it took quite a time to nish it. I have done many things that I at rst didn't know had to be done and consequently I haven't done all things that I planned to do. When I started, I knew almost nothing about geology and nite element methods, and therefore it took more time than I had thought to nd, read and understand all background material needed. I really must express my gratitude toHarry Hanse. He has patiently answered all my questions and taken time to instruct me in nite element methods and other areas that was new to me. I'm also very grateful to Johan Bystrom for oering to read through this text and for all his valuable comments.

I also want to thank my family and all other relatives and friends for everything. I've often spent much time studying and then you haven't seen much of me, but I haven't forgotten you.

Great thanks also to professorLars-Erik Persson and others at the De- partment of Mathematics for believing in me and accepting me for postgrad- uate studies.

I also want to thankCSN for student grants andColdtech for granting me nancial support for this thesis.

1 Introduction

The purpose of this thesis is to examine a mathematical model of a phe- nomenon called frost heave. This introduction to the subject contains a brief account of the mathematical model, and should give a rough under- standing of the physical background. A major drawback of this approach is a great lack of justications for involved assumptions and simplications.

We will take it for granted that the omitted material is reliable. The sceptic should consult research literature to verify all details. A thorough survey of the physics of freezing phenomena in soils can be found in [8]. This introduc- tion is based on that text, the Ph. D. thesis of Daichao Sheng [12, chapter 1{2.3], an article by O'Neill [9] and an article by O'Neill and Miller [10].

2 How to read this text

2.1 What you should read 2.1.1 Physical background

I think that any reader not already acquainted with the rigid ice model should at least browse through the physical background (section 3) and every reader should denitely take a closer look at the equations and boundary conditions that it leads to (section 3.4).

2.1.2 The mathematical problem

In section 6 all results from the preceding is collected to a mathematical model. This is an absolute minimum of what you have to read to understand what comes next.

2.1.3 Numerical calculations

Section 7 shows how the mathematical model of section 6 can be simulated using a nite element method.

2.2 Sections that can be ignored 2.2.1 Change to dimensionless variables

In section 4, we translate the problem to dimensionless variables. As far as we know, this is the rst time this is done in print. You could skip this section and jump directly to section 6.

2.2.2 Material properties

The rigid ice model contains a number of physical parameters that have to be specied. This is done in section 5. If you don't care, you can skip this section and go directly to section 6.

2.2.3 The appendices

Appendix C contains a Maple worksheet for some of the computations of the nite element method that depends on the choice of \trial functions".

Thus we can change trial functions without having to do much new work.

Appendix D explains how we have handled the boundary conditions.

Appendix E outlines how to implement the nite element method solu- tion in Matlab.

The remaining appendices contain material that many readers either are acquainted with or have no interest in.

3 Physical background

3.1 Frost damage

Frost action in the ground is a great problem in the urban cold regions of the world. It causes problems to various constructions like highways, airport pavements, railways, pipelines and building foundations. In Sweden, Finland and Norway, the costs during the frost damage period each year were 1985 estimated at between 10 and 35 million U.S. $ for road maintenance, which causes trac restrictions that cost the road users some 15 to 20 million $.

3.2 Frost heave

There are two main types of frost damages: frost heaveandthaw weakening. O'Neill and Miller [10] dene the heave of freezing soil as the upward dis- placement of the soil surface due to many dierent factors, such as expansion of the water when freezing to ice, moisture ow and the type of soil:

Coarse sands and gravels

are not susceptible to heave.

Silts

can produce spectacular heave, provided that the load to be heaved is small and the water table is high.

Clays

never produce spectacular heave, but can heave very large loads.

The soil freezes in any of three dierentfreezing modes, due to properties of the soil and boundary conditions of the narrow zone where the freezing takes place. Sometimes all three modes are transversed:

Primary heave

consists of a layer of ice growing on top of unfrozen soil, through which water is drawn to join the ice. The ice might be covered by a layer of frozen soil, but essentially no ice penetrates the unfrozen soil.

Secondary mode of heave:

Ice penetrates the soil and liquid ows to- wards it, producing layers of pure ice, called ice lenses. This heave is sometimes quite powerful. Water might ow either from or to the underlying unfrozen soil.

No heave

at all may be the result if all the ice formed penetrates the soil.

Water will then be expelled from a descending zone of freezing into the unfrozen soil.

3.3 Derivation of the rigid ice model

Over 100 dierent criteria in uence the frost heave phenomenon and simple conceptual models and rules of thumb have repeatedly failed to describe it.

Since 1961 many dierent theories have been developed and according to Sheng [12], the most complete one from a microscopic and physical point of view is the so called rigid ice model, which, however, has some drawbacks:

 It is very complex and no ecient solution has been found that is useful in practice.

 A large number of physical parameters must be specied to make the model practically applicable.

 The model has not been fully veried by experimental data.

 The non-linear convection-diusion equations of the model have been solved by the standard Galerkin nite element method, without sys- tematic investigation of computational behaviour like convergence and stability. However, the numerical solutions obtained become unstable during the terminal stage of heave. Thus deeper knowledge of appro- priate values of time step and grid size is needed for every particular application.

The rigid ice model describes the secondary mode of heave, which is the most interesting mode from an engineering point of view, since it can lift the greatest loads. Important boundary conditions are

The magnitude of the load

(including already frozen soil) to be heaved.

Heat ow

extracted from the cold (upper) boundary of the frozen zone and from unfrozen soil.

Hydrologic conditions

at greater depths in uence water pressure at the warm boundary.

An essential concept in the rigid ice model is the existence of a so called frozen fringe (gure 1) between the latest (still growing) ice lens and the unfrozen soil. In the frozen fringe, both water ow and ice movement takes place. O'Neill and Miller [10] describe in detail a process calledregelation, which assures that water is constantly melting and freezing in the frozen fringe in such a way that the ice in the pores is rigidly connected to the growing ice. Also, it moves upwards like one rigid body, with a velocity equal to the heaving rate.

An simple demonstration of mechanically induced regelation involves a thin wire, with heavy weights at each end, draped over a block of ice with the weights dangling. Ice beneath the wire is subjected to high pressure, which lowers the freezing temperature so that the ice melts. Melt water ows up around the wire and refreezes immediately, so that the block of ice shows no evidence of the passage of the wire through it. Accordingly the wire slowly sinks through the solid block of ice!

Before developing the rigid ice model of the soil pores, all physical quan- tities that will be needed are presented in table 1.

Heat flow

unfrozen water mineral particles

x b

x f

latest ice lensfrozen fringe

moisture migration

x y z

Figure 1: A schematic diagram of the frozen fringe.

3.3.1 Restrictions

The following derivation is limited to solute-free, air-free soils of negligible compressibility.

We assume that the processes involved are suciently slow and the con- tact between phases suciently pervasive so that local thermal equilibrium between the phases is achieved.

3.3.2 Water ow

Since there is no internal production or consumption of water, the increment of total water mass in a unit volume of soil must equal the mass ux out of

it. @

@t⁽wW +iI+vV) +^r
(w

v

^w+i

v

ⁱ+v

v

^v) = 0 (3.1) Microscopically, parameters such as intergranular forces vary sharply, but they may be smoothed, in the macroscopic view, since they manifest themselves simply as net intergranular force per unit area of soil. It must, however, be assumed that the frost heave process is slow enough for potential micro-variations in the ice to be damped out.

The liquid ux is assumed to follow theDarcy law, which states that the rate of water ow is proportional to the pressure gradient:

v

^w=^,k^rH =^,k^r uw

wg ⁺x^=^, k

wg^ruw^,

0

@

k0 0

1

A (3.2)

whereH is the total head of water pressure. Substitution into (3.1) yields

@t@ ⁽wW +iI+vV) +^r
(i

v

ⁱ+v

v

^v) =^r

0

@k

g^ruw+

0

@

kw

00

1

A 1

A

(3.3) where, the eects of the vapour phase and ice ux are usually neglected.

3.3.3 Heat conservation

The change of temperature in moist soil depends on several factors.

 The sensible heat is increasing at speedC^@T_@t W=m³.

 The net heat ux due to conduction is^r

Q

^h W=m³.

 The net heat ux due to convection associated with the movement of the three water phases is^r
(TCw

v

^w+TCi

v

ⁱ+TCv

v

^v) W=m³.

symbol meaning unit

A known constant |

B known constant J=(m³
K)

C volumetric heat capacity of soil J=(m^3
C)

Cⁱ volumetric heat capacity of ice J=(m^3
C)

C^v volumetric heat capacity of vapour J=(m^3
C)

C^w volumetric heat capacity of water J=(m^3
C)

F volumetric fraction air |

G volumetric fraction grain |

H total head of water pressure m H2O

I volumetric fraction ice |

L specic latent heat of water J=kg

P overburden pressure Pa

Qh heat ux W=m²

T temperature ^C

T0 bulk water freezing point at atmospheric pressure K

V volumetric fraction vapour |

W(^iw) volumetric fraction water |

e

x vertical unit vector in x-direction (gure 1) |

g acceleration of gravity m=s²

h;h0 height of soil column and its initial value m

k(W) hydraulic conductivity m=s

n porosity of soil,n= 1^,G |

r^pc rate of ice water phase change kg=(m³
s)

t time s

uⁱ ice pressure Pa

uⁿ eective pore pressure Pa

u^w pore water pressure Pa

v^h rate of heave (=speed of ice lensing) m=s

v

i volumetric ice ux in frozen fringe, (m³ ice)=(s
m²soil) m=s

v

v volumetric vapour ux in frozen fringe, (m³ vapour)=(s
m² soil) m=s

v

w volumetric water ux in frozen fringe, (m³ water)=(s
m² soil) m=s

v^w the x-component^ex
vwof^vw m=s

x See gure 1 m

x^b Warm surface of latest ice lens (gure 1) m

x^f frost front (gure 1) m

x^w xed warm boundary (gure 2, page 10) m

 thermal conductivity of soil W=(m
^C)

^{iw iw}=u^i,u^w Pa

^{w a w a=}u^w,u^{a Pa}

ⁱ density of ice kg=^m3

^v density of vapour kg=^m3

^w density of water kg=m³

^iw ^i,^{w kg}=^m3

hⁱ the unit of variable |

r 0

@

@

@x

@

@y

@

@z 1

A |

(W) dened in equation (3.20f) |

Table 1: All variables and notations used during the development of the rigid ice model. Pressures are gauge pressure measured against a standard atmo- sphere.

Bold face

printing denotes vectors and | denotes dimensionless (dened in section A.1) quantities. See also gure 1.

 The latent heat Lrpc W=m³ is released because of water-ice phase change.

Since energy is conserved it follows that

C@T@t ^+r

Q

^h+^r
(TCw

v

^w+TCi

v

ⁱ+TCv

v

^v) =Lrpc (3.4) The heat ux

Q

^h is assumed to follow the Fourier law

Q

^h=^,
^rT (3.5)

where the thermal conductivity (T) of the soil may depend on T. The total rate of ice-water phase change rpc depends both on the movement of the already frozen ice and the degree of melting/freezing in a point following the ice.

rpc=i^r

v

ⁱ+i@I

@t ^(3.6)

Substitution of (3.5) and (3.6) into (3.4) yields

C@T@t ^+r
(TCw

v

^w+TCi

v

ⁱ+TCv

v

^v) =^r
(^rT) +iL



r

v

ⁱ+ @I

@t



(3.7)

3.3.4 Phase equilibria and transport equations

Let us now look at those parts of the rigid ice model that concern the regelation process in the frozen fringe introduced in section 3.3, page 4.

This theory builds on the proposition that, for whatever reason, water very close to a grain surface is attracted to it by a force whose strength decays with the distance from the surface. Also, this force is greater for liquid than for air or ice.

Then, the grains in the frozen fringe will be surrounded by an adsorbed lm of mobile water. If a bubble of air approaches the immersed grain it will experience a disjoining force, balanced by the pressure wa = uw ^,

ua. The mean curvature of the water-air interface outside force eld is also determined by wa, which thus also is a measure of to which degree water replaces air in capillary space. Consequently, in the unfrozen soil the volumetric fractions of waterW, air F and grainG are complementary functions ofwa:

W(wa) +F(wa) +G= 1 (3.8)

Likewise, if surface energies of crystal-melt interfaces of ice are independent of crystallographic orientation at temperatures not too far below 0^C, then W(iw) +I(iw) +G= 1 (3.9) for air-free soil, where

iw=ui^,uw (3.10)

and an expression forI(iw) is chosen in section 5.1.

Thermodynamic reasoning, as found in standard textbooks, leads to a generalized form of the Clapeyron equation. In spite of some simplifying physical and mathematical approximations, this equation is quite accurate over the range of temperature and pressures of interest in the rigid ice model:

uw

w ^,ui

i = LT

T0 (3.11)

Substitutingui from (3.11) into (3.10) yields

iw=

i

w ^,1



uw^,iLT

T0 =Auw+BT (3.12) whereA andB are known constants, A is dimensionless¹ and

hBⁱ=

iL T0



= (kg=m³)
(J=kg)

K = Jm³
K (3.13)

Now assume that variations ofk,uw, andw as well as components of

v

ⁱ are negligible in the y- and z-direction. Also suppose that the eects of vapour phase can be neglected in the mass equation (3.3). Then

v

ⁱ =vhI

e

^x (3.14)

wherevh(the frost heave rate) is what we want to compute in our simulation later. Insertion of (3.9) into (3.3) yields

@t@ ⁽w(1^,G) +iwI) + @

@x⁽ivhI)^, @

@x

k g@uw

@x ⁺kw



= 0 (3.15) where

iw =i^,w (3.16)

Suppose that the compressibility of the soil, water and ice is negligible, i.e.,

@w

@t ⁼ @i

@t ⁼ @G

@t ⁼ @i

@x ^{= 0}

1That is, ^hAⁱ is some positive real number. See the discussion about dimensions in section A.1 at page 30.

and since the ice is moving rigidly,

@vh

@x ^{= 0 (3.17)}

Then equation (3.15) and the chain rule yields

iw dI diw@iw

@t ⁺ivh dI diw@iw

@x ^, @

@x

k g

@uw

@x ⁺gw



= 0

Insertion of (3.12) and use of the chain rule yields the rst rigid ice model dierential equation.

iwA dIdiw@uw

@t ⁺iwB dIdiw@T

@t ⁺ivh



A dIdiw@uw

@x ⁺B dIdiw@T

@x



,

@x@

k g

@uw

@x ⁺gw



= 0 It is now assumed that the processes involved are slow enough to allow neglecting the convective terms in the heat equation (3.7). Provided that variations of T and  in the y- and z- directions are negligible, insertion of (3.14) into (3.7) yields

C@T@t ^, @

@x



@T@x



,iL^@(vhI)

@x ⁺@I

@t



= 0 Since I =I(iw) usage of (3.12), (3.17) and the chain rule yields

C@T@t ^, @

@x



@T@x



,iL^vh dI diw



A@u@x^{w +}B@T@x



+ dI diw



A@u@t^{w +}B@T@t



= 0 This is the second rigid ice model dierential equation, which can be written



C^,iLB dIdiw

@T

@t ^,iLA dIdiw@uw

@t ^,iLvh dI diw



A@u@x^{w +}B@T@x



,

@x@



@T@x



= 0

3.3.5 Rate of heave

O'Neill and Miller [10] have investigated two methods of calculating vh, which appears in both rigid ice model dierential equations. Both methods are based on mass conservation, and redundancy with the general mass

soil ice lens

warmest ice lens

h

x

x

x b

f

fixed warm boundary w

x

frozen fringe

unfrozen soil

Figure 2: A schematic view of soil, ice lenses and frozen fringe.

balance equation (3.3) is avoided in two dierent ways: Balances are drawn up either ^➀ only over the lens-fringe interface at xb or ^➁ over the entire frozen fringe and unfrozen soil down to a xed warm boundary at x =xw

(gure 2). We will choose x-axis such that xw = 0. In both cases, we take advantage of our knowledge that the warmest ice lens consists of pure ice, moving at a rate ofvh.

➀: At the frozen fringe side of xb, unfrozen water ows into the lens and freezes at a rate of wvw while the soil ice ow into the lens is iIvh. The mass ux at the lens side ofxb is ivh. Mass conservation at the plane throughxb yields the equation

ivh =iIvh+wvw at x=xb

which together with the onedimensional form of the Darcy law (3.2) yields an expression forvh:

vh=^,k^{@u@xw +}wg

ig(1^,I) atx=xb (3.18) Here, all quantities at the right-hand side are evaluated at x =xb and are determined fromuw and T (suitable expression for I is discussed in section 5.1). However, since the ice is assumed to be rigidly interconnected, the values ofvh apply everywhere in the ice.

➁: Another expression forvh derives from a mass balance over the space between xf and xw The total mass ux through the plane at xb is given byivh, while the mass ux through the plane at xw is given by wvw(xw).

The dierence between these two uxes must equal the rate of change of the aggregate mass content inbetween:

ivh^,wvw(xf) = @

@t

Z x^b

x^w(wW +iI)dx

The warm boundary xw is xed over time. The value value of xb is changed instantaneously each time a new lens starts to grow, and is constant the remaining time. Thus, except for a nite number of points in time,xw

andxb are constant and after eliminatingW using equation (3.9), the above expression may be rewritten as

vh= w

ivw(xf) +iw

i @

@t

Z x^b x^w Idx

or, after applying the one-dimensional form of Darcy's law (3.2), vh =^, k

ig @uw

@x







x=x^w,

wk(xw)

i +iw

i @

@t

Z x^b

x^w Idx (3.19)

3.3.6 Ice lens initiation

The pressure in the frozen fringe depends both on the ice pressure ui, the water pressure uw and the pressure of everything above, the overburden pressure P. If water is owing from the unfrozen soil thenuw <0, but the eective pore pressure² un, which is a weighted sum ofui and uw, might be positive. A necessary condition for initiation of an ice lens is that un =P, so that the regelation process can separate \ oating" solid particles.

undef= uw+ (1^,)ui =P

where  depends on W.  is called stress partitioning factor. Given a relation between iw and W, the relation iw(I) can be determined from laboratory experiments (see section 5.1). Then a theoretically based sugges- tion [13] yields an integral that can be evaluated by e.g., quadrature:

(W) = 0:5 n



W ^{, 0}:3

iw(W)

Z n

W iwI⁰(iw)diw(W)



where nis the porosity of the soil.

2Also callednet pressure.

3.3.7 Boundary conditions

There are two boundary conditions at the warm surface of the latest (grow- ing) ice lens. One states that the overburden pressureP at xb (see gure 1 at page 4) is fully carried by the ice lens. Insertion of P = ui into the Clapeyron equation (3.11) yields

(1 +A)uw+BT =P at x=xb

Also, the rate of heat removal must equal the rate of liberation of latent heat due to freezing atxb.

,



@T@x



x^b++



@T@x



x^b,=Lwvw(xb)

Wherexb+ denotes the limit whenx approachesxb from the cold side and xb^,is the corresponding limit from the warm side.

The total heighth(t) of the soil column is h(t) =h0+

Z t 0 vhdt

3.4 Restatement of the complete model

Finally, the complete rigid ice model is summarized. All symbols needed are dened in table 2. The rigid ice model dierential equations are

iwA dIdiw@uw

@t ⁺iwB dIdiw@T

@t ⁺ivh



A dIdiw@uw

@x ⁺B dIdiw@T

@x



,

@x@

k g

@uw

@x ⁺gw



= 0 (3.20a) and



C^,iLB dIdiw

@T

@t ^,iLA dIdiw@uw

@t ^,iLvh dI diw



A@u@x^{w +}B@T@x



,

@x@



@T@x



= 0 (3.20b) We want to compute the frost heave rate vh. One way to do this is to use the formula

vh=^,k^{@u@xw +}wg

ig(1^,I) atx=xb (3.20c)

Heat flow

unfrozen water mineral particles

x b

x f

latest ice lensfrozen fringe

moisture migration

x y z

Figure 3: A schematic diagram of the frozen fringe. (Same as gure 1.) but we will instead use the expression

vh=^, k

ig @uw

@x







x=x^w,

wk(xw)

i + iw

i @

@t

Z x^b

x^w Idx (3.20d) A necessary condition for initiation of a new ice lens is that un=P, so that the regelation process can separate \ oating" solid particles:

un def= uw+ (1^,)ui =P (3.20e) where depends on W. Given a relation betweeniw and W, the relation

iw(I) can be determined from laboratory experiments. Then a theoreti- cally based suggestion [13] yields an integral that can be evaluated by e.g., quadrature:

(W) = 0:5 n



W ^{, 0}:3

iw(W)

Z n

W iwI⁰(iw)diw(W)



(3.20f) wherenis the porosity of the soil.

The rst boundary condition states that the overburden pressure P at xb is fully carried by the ice lens.

(1 +A)uw+BT =P at x=xb (3.20g)

The second one says that the rate of heat removal must equal the rate of liberation of latent heat due to freezing atxb.

,



@T@x



x^b++



@T@x



x^b,=Lwvw(xb) (3.20h) Wherexb+ denotes the limit whenx approachesxb from the cold side and xb^,is the corresponding limit from the warm side.

symbol meaning (value and) unit

A known constant ^,₉₉₇⁸⁰ |^a

B known constant ^,6:7858₆₀₇
¹⁰⁸ J=(m³
K)^a C volumetric heat capacity of soil (section 5.4) 1:76
10⁶ J=(m^3
C)

I volumetric fraction ice |

L specic latent heat of water 333
10³ J=kg

P overburden pressure^bc Pa

T temperature ^C

T0 bulk water freezing point at atmospheric pressure 273.15 K

W(^iw) volumetric fraction water |

e

x vertical unit vector in x-direction (gure 3) | g acceleration of gravity (in Lulea, Sweden) 9.823 m=s²

n porosity of soil^b,n^{= 1,}G^{(table 1) |}

h;h0 height of soil column and its initial value m k(W) hydraulic conductivity (=permeability) See sec-tion 5.3. m=s

t time s

uⁱ ice pressure^c Pa

uⁿ eective pore pressure^c Pa

u^w pore water pressure^c Pa

v^h rate of heave (=speed of ice lensing) m=s v^w frozen fringe water ux in the x-direction,

(m³ water)=(s
m² soil) m=s

x See gure 3 m

x^b warm surface of latest ice lens (gure 3) m

x^w xed warm boundary (gure 2, page 10) m

⁽T⁾ thermal conductivity of soil^b W=^(m
C)

^{iw iw}=u^i,u^wc Pa

ⁱ density of ice (at -4^C) 917 kg=m³

^w density of water 997 kg=m³

^iw ^i,^w -80 kg=m³

(W) dened in equation (3.20f) |

Table 2: All variables and notations used in the nal statement of the rigid ice model.

Bold face

printing denotes vectors and | denotes dimensionless (dened in appendix A.1) quantities. Se also gure 3.

aCalculated from equation (3.12) without introducing any round-o errors.

bThe value diers for dierent types of soil, as can be seen in table 3.

cGauge pressure measured against a standard atmosphere.

Finally, the total height h(t) of the soil column is h(t) =h0+

Z t

0 vhdt (3.20i)

4 Change of variables

It is shown in appendix A how a change into dimensionless variables trans- forms the problem of section 3.4 into a form better suited for our purposes.

(Note that from now on ^{boldface italic} font denotes dimensionless vari- ables, not vectors. This will simplify the notation.)

symbol parameter unit 1: sand 2: silt 3: clay

^d dry density of soil kg=m³ 2:0
10³ 1:4
10³ 1:1
10³

n porosity | 0.25 0.005 0.006

k^{s saturated}permeability m=^{s 0.1 10,5 10,8}

^s thermal conductivity

of soil solids m^WC 2.3 2.3 2.3

T^t average top

temperature ^C -5.0 -5.0 -5.0

T^b average bottom

temperature ^C 5.0 5.0 5.0

h0 initial heigth (h, gure 2, page 10) of

soil column m 0.20 0.20 0.20

W0

initial volumetric fraction unfrozen water at -1.0

C(equation (5.3)) | 0.25 0.005 0.006

 parameter in

equation (5.3) ^C -1.19 -1.79 -2.65

 parameter in

equation (5.3) | 1.6 2.85 1.99

Table 3: Soil properties for some dierent types of soil. This is an arbitrary choice of values that is supposed to be representative for the three types of soil. Note that the rigid ice model is derived for silt (see also section 3.2).

 and  varies for dierent samples of the same type and here we use the values of three of Kujala's [6, table 3] samples. The remaining parameters are collected from an example by Sheng [12, table 3.1].

Kujala's samples had other W0-values, but we want to start without any ice, i.e., with W0 =n and hope that this still proves to be a realistic set of parameters.

It is found in appendix A that the variables

C = C

B ^{P =} PA

Li ^ui= uiA Li

un= unA

Li ^uw= uwA

Li iw= iw

Li

T =^,TB Li = T

T0 ^Tb= Tb

T0 ^Tt = Tt

T0

k= k

A^pL ^ks= ks

A^pL ^{vh =} vh

pL

vw= vw

pL ^h= hg

L ^x= xg

L

xb= xbg

L ^xw= xwg

L ^t= tg

pL

=^, g

BL^{32 d=} d

i ^w= w

i

iw= iw

i = ^A=A

I =I ^W =W ⁿ=n

transforms the problem into the form described by equations (A.1)-(A.9).

Equation (A.1) and (A.2) can be written in a more convenient form, using the vector

f =



uw

T



as shown below. Also, equations (A.3)-(A.9) are repeated.

4.1 The dimensionless rigid ice model

We are interested in how ^u_w and ^T depends on ^x and ^t. The dierential equations to solve are



iw ^,iw

,1 1

 d^I d_iw@^f

@^{t +}

0 0 0 ^,C

@^f

@^{t +}

 1 ^,1

,1 1



vh d^I d_iw@^f

@^x + @

@^x



,k 0

0 ^

@^f

@^x



=

 @@^x(^k_w) 0



(4.1) where two dierent expressions for^v_h are

vh=^,k@@uxw +^_w^A

1^,I (4.2)

and

vh =^,k @^u_w

@^x









x=^xw^,w^Ak(^x_w) +^_iw @

@^t

Z x^b

x^{f I}d^x (4.3)

(We will, however, only use the last one.) The boundary conditions at ^x_b are

(^A+ 1)^u_w^,AT =^P at ^x=^x_b (4.4) and

,





@^T

@^x



xb⁺+





@^T

@^x



xb^,=^_w^v_w(^x_b) (4.5) Wherexb+ denotes the limit whenx approachesxb from the cold side and xb^, is the corresponding limit from the warm side. With

(^W) = 0:5

n



W ,

0:3

iw(^W)

Z n

W iw^I0(_iw)d_iw(^W)



(4.6) a necessary condition for initiation of a new ice lens is that

un^def= ^u_w+ (1^,)^u_i=^P (4.7) Finallly,

h(^t) =^h0+

Z t

0 ^vh(^t)d^t (4.8)

5 Material properties

The rigid ice model contains a number of physical parameters that have to be specied. Especially parameters as the stress partitioning factor (W), the unfrozen water content W(iw) and the hydraulic conductivityk(W) of the frozen fringe are dicult to determine. Sheng [12] describes many dierent suggestions, obtained by regression of experimental data. We choose some of them here, but without more experimental data, no suggestion stands out as superior to the others.

5.1 Ice content

^I

as a function of

iw

It follows from equation (3.9) that it suces to nd eitherI(iw) orW(iw).

Since the sixties, many researchers have examined the volumetric fraction W of water³ in the frozen fringe. A common approach is to experimentally

3Koopmans and Miller [4] made the useless suggestion W = ^P5k₌₀A^kln^!iwiw



k. Since ^!iwiw is of dimension^pressure_tension =_length¹ , this yields something horric with \dimension"

P5

k=0dim(A^k)^ln^_length^{1 k}!

determineW as a function of the temperatureT. This is suitable for many applications, but we need W (orI) as a function of iw.

Either ui or uw can be eliminated from Clapeyrons equation (3.11), by insertingiw=ui^,uw:

iw=

i

w ^,1



uw^,iL TT0 (5.1a)

iw=



1^,w

i



ui^,wL TT0 (5.1b) Black [2] explains how the work of Colbeck [3] yields that

uw ^0^ui

(remember that pressures are gauge pressures measured against a standard atmosphere). Insertion of this andw ^i into (5.1) leads to the inequality

,iL TT0 ^iw^,wL TT0 (5.2) Of course, only one suchiw is correct. However, Black determines a rela- tionship of the type

W =a
(iw)^b

from measurements ofW againstT, using two dierentiwthat full (5.2), chosen to maximize the dierence between the results. The dierence is, however not visible in Black's graph of the two achieved curves.

Here we will use another model, based on more experimental data, but we will assume that the change of independent variable fromT to iw can be done as above. Furthermore we make the reasonable assumption that this holds even for thederivative _d^dWiw.

Kujala [5, 6] has measured the unfrozen water contents of 68 soil samples from clay to sand, by the method of nuclear magnetic resonance, and found that the unfrozen water content can be approximated well by the so-called Weibull survival function

n,I =^W =^W0e^,(T) (5.3) where is the temperature at which 63.2% of the water is frozen and ^W0 is the initial volumetric water content⁴. Table 3 (page 15) shows the values of and fore some dierent soil samples. Note that this function only is dened for _T ^0, i.e., for T ^0. (For T >0 we should have ^W = 0, but of courseT ^0 in the frozen fringe.)

4As before, we use ^{boldfaceitalic}font to indicate that^W0 andare dimensionless.