CRREL REPORT 95-8

Permafrost Formation Time

Virgil J. Lunardini April 1995 PERMAFROST SURFACE TEMPERATURE (

°C)

0 160,000

TIME BEFORE PRESENT (years)

–8

–16
Abstract: The age of permafrost is closely linked to

the time required for soil systems to freeze, since the

permafrost must be at least as old as the formation

time. Cycles of freeze–thaw will complicate the relation between the freeze rate and the age. A model

based on pure conduction heat transfer with freeze–

thaw is used to predict the time required for a given

thickness of permafrost to develop, either heterogenetically or syngenetically. The formation time is a

function of the long-term geothermal gradient (initial

temperature of the thawed soil), the ratios of the froHow to get copies of CRREL technical publications:

Department of Defense personnel and contractors may order reports through the Defense Technical Information Center:

DTIC-BR SUITE 0944

8725 JOHN J KINGMAN RD

FT BELVOIR VA 22060-6218

Telephone 1 800 225 3842

E-mail help@dtic.mil

msorders@dtic.mil

WWW http://www.dtic.dla.mil/

All others may order reports through the National Technical Information Service:

NTIS

5285 PORT ROYAL RD

SPRINGFIELD VA 22161

Telephone 1 703 487 4650

1 703 487 4639 (TDD for the hearing-impaired)

E-mail orders@ntis.fedworld.gov

WWW http://www.fedworld.gov/ntis/ntishome.html

A complete list of all CRREL technical publications is available from:

USACRREL (CECRL-LP)

72 LYME RD

HANOVER NH 03755-1290

Telephone 1 603 646 4338

E-mail techpubs@crrel.usace.army.mil

For information on all aspects of the Cold Regions Research and Engineering Laboratory, visit our World Wide Web

site:

http://www.crrel.usace.army.mil

zen to thawed thermal properties, and the temperature

history of the upper surface of the permafrost (higher

than the air temperature). The simple theory allows

universal graphs to be produced that predict the formation time for a given thickness of permafrost. Realistic soil property ratios and paleotemperature scenarios will then lead to estimates of the formation time of

permafrost for a specific site. The model indicates that

deep permafrost (more than 1500 m) requires formation times on the order of the complete Quaternary

Period.

Cover: Reconstruction of paleotemperature history of Prudhoe Bay, Alaska

(after Osterkamp and Gosink 1991).
CRREL Report 95-8

Permafrost Formation Time

Virgil J. Lunardini April 1995

Prepared for

OFFICE OF THE CHIEF OF ENGINEERS

Approved for public release; distribution is unlimited.

US Army Corps

of Engineers

Cold Regions Research &

Engineering Laboratory
ii

PREFACE

This report was prepared by Dr. Virgil J. Lunardini, Mechanical Engineer, Applied Research Division, Research and Engineering Directorate, U.S. Army Cold Regions Research and Engineering

Laboratory. This study was primarily funded by U.S. Army Strategic Environmental Research and

Development Program, Deep Permafrost Borehole Sites in Alaska, 2 6.00 78 203, and DA Project

4A161102AT24, Task BS/0045, Climate Impacts on Permafrost Distribution and Integrity.

The author thanks Dr. Gary Phetteplace and Michael Ferrick of CRREL for their technical review of this report.

The contents of this report are not to be used for advertising or promotional purposes. Citation of

brand names does not constitute an official endorsement or approval of the use of such commercial

products.
CONTENTS

Preface ............................................................................................................................................. ii

Nomenclature .................................................................................................................................. v

Introduction ..................................................................................................................................... 1

Origin and existence of permafrost............................................................................................. 2

Paleotemperatures ....................................................................................................................... 2

Theory.............................................................................................................................................. 6

Heterogenetic freeze relations .................................................................................................... 7

Heterogenetic model verification ............................................................................................. 10

Sensible and latent heats ........................................................................................................... 10

Neumann solution ..................................................................................................................... 11

Steady-state solution ................................................................................................................. 11

Effect of quaternary freeze–thaw cycles................................................................................... 11

Syngenetic growth of permafrost.............................................................................................. 12

Syngenetic model verification .................................................................................................. 14

Discussion...................................................................................................................................... 14

Prudhoe Bay, Alaska ................................................................................................................ 18

Deep permafrost........................................................................................................................ 19

Conclusions ................................................................................................................................... 20

Literature cited............................................................................................................................... 20

Appendix A: Soil properties and ratios ......................................................................................... 23

Appendix B: Quaternary cyclic thermal modulation..................................................................... 25

Appendix C: Heat balance integral equations for syngenetic permafrost growth ......................... 29

Appendix D: Energy flows at the permafrost base ........................................................................ 35

Appendix E: FORTRAN programs for numerical quadrature of energy equation ....................... 39

Abstract.......................................................................................................................................... 45

ILLUSTRATIONS

Figure

1. Paleotemperature variations ..................................................................................................... 3

2. Pleistocene temperature fluctuations ........................................................................................ 3

3. Vostok ice core temperature inferences ................................................................................... 4

4. Permafrost surface paleotemperature model for Prudhoe Bay, Alaska.................................... 4

5. Permafrost surface paleotemperature model for East Siberian location .................................. 5

6. Reconstructed temperature for northern China during last ice age .......................................... 5

7. Mean paleotemperature departure, Prudhoe Bay, Alaska, during one glacial cycle ................ 5

8. Paleotemperature history at Mackenzie Delta, Canada ............................................................ 6

9. Freeze of a semi-infinite region with linear initial temperature ............................................... 7

10. Temperature distribution after initial freeze of soil ............................................................... 11

11. Bottom thaw of permafrost .................................................................................................... 12

12. Thermal state at end of melting and long-term equilibrium................................................... 12

13. Geometry and coordinate system for freeze of a semi-inifinite medium with moving

upper surface ................................................................................................................... 12

iii
iv

Figure

14. Coordinate system for stationary upper surface ..................................................................... 13

15. Formation time of permafrost, saturated mineral soils, φ = 0, ε = 0.4 ................................... 15

16. Formation time of permafrost, saturated mineral soils, φ = 0, ε = 0.379 ............................... 15

17. Formation time of permafrost, saturated mineral soils, φ = 0, ε = 0.3 ................................... 16

18. Formation time for permafrost to reach 90% equilibrium thickness, saturated mineral soils 16

19. Effect of soil porosity on permafrost thickness after 15,000 years, saturated soils,

G = 0.0286°C/m .............................................................................................................. 16

20. Syngenetic growth of permafrost, saturated mineral soils, φ = 0, ε = 0.379, ST = 0.15........ 17

21. Syngenetic growth of permafrost, saturated mineral soils, φ = 0, ε = 0.379, ψ = 0.15 ......... 17

22. Growth of permafrost at Prudhoe Bay, Alaska ...................................................................... 19

TABLES

Table

1. Movement of Tf isotherm, homogeneous soil for heterogenetic freezing .............................. 10

2. Movement of Tf isotherm, homogeneous soil for syngenetic freezing .................................. 14

3. Paleotemperature scenarios, Prudhoe Bay, Alaska ................................................................ 17

4. Effect of previous cooling on permafrost growth time .......................................................... 19

5. Formation time of deep permafrost ........................................................................................ 19

6. Extreme predicted permafrost thickness ................................................................................ 20
v

β δ/X

β′ d

d

β

σ

δ depth of temperature disturbance

ε soil porosity

γN X t / 2 α1

∆T To – Tf

∆T1 Tf – Ts

ξ integrated temperature

ρ density

σ

G

T

X

∆ 1

σN dimensionless freeze depth for Neumann

solution

τ α1

1

2

G

T

t

∆





 





φ T T

T T

o f

f s

−

−

ψ

U T

G

∆ 1

α1

Subscripts

f solidification or freeze value

e equilibrium value

i initial value

o initial surface value

s surface

u thawed

1,2 frozen and thawed, respectively

∞ steady-state value

NOMENCLATURE

c specific heat

c12 c1 /c2

C ρc, volumetric specific heat

f fraction of geothermal energy used for

melt

g

α σ

ββ σ

21 1

2 2

1 ( ) +

[ ] ( ) + +

+

g′ dg/dσ

G geothermal gradient (slope of initial temperature distribution)

k thermal conductivity

kg thermal conductivity of soil solids

k12 k1/k2

l latent heat of solidification

L ρl, volumetric latent heat

M m2π2t/4tc

qg geothermal energy flow

QS,QL sensible and latent heats

S soil saturation parameter

ST 

c T T 1

l ( ) 1 s − , Stefan Number

t time

tc Xo

2/4αf, characteristic time for permafrost temperature changes

tm time to complete melt of permafrost

T temperature

U deposition rate at surface of permafrost

x depth coordinate

xw volumetric water content

X permafrost thickness

α thermal diffusivity

α12 α1/α2
INTRODUCTION

The age of permafrost is of interest to biologists, geophysicists and engineers. Clearly, permafrost must

be at least as old as the time it took for it to form; thus, the formation time of permafrost can be considered

its minimum age. A volume of permafrost can be much older than this since it may exist for many years after formation. This report will examine the formation time of permafrost using a pure heat conduction model. As we shall see, the surface temperature history of the soil mass is critical for any prediction of the permafrost formation time. Since the formation time of permafrost is expected to be on the order of millennia,

it is necessary to examine the geophysical record to obtain some bounds on realistic surface temperatures

that the Earth has experienced during the time when permafrost was growing. First, we will discuss permafrost and paleotemperature scenarios, then we will formulate a mathematical model of permafrost growth,

and, finally, we will examine some predictions by the model of permafrost formation times.

Permafrost is a widespread phenomenon that has been and still is greatly misunderstood. The term

“permafrost” is generally attributed to S.W. Muller (1945), who apparently coined the name in place of the

more awkward terms: permanently frozen ground or permanent frost. Bryan (1946) suggested the term

“pergelisol,” but this has not been adopted except in the French literature. In order to understand the concept, let us look at a general definition given in Lunardini (1981a):

Permafrost describes the thermal condition of earth materials (sand, glacial till, organic matter, etc.)

when their temperature remains at or below 32°F (0°C) continuously for a significantly long time,

but not necessarily for an entire geological period. It does not include earth materials that drop below 32°F during one winter and remain below 32°F through the following summer and into the next

winter, although for practical engineering purposes such materials may be included.

Clearly, permafrost is not so much a material as it is the thermal state of ordinary soil systems. It does

not include systems that are at or below 0°C, but contain no earth materials, e.g., ice caps, glaciers and icebergs. There is no agreement on the minimum time during which the material must remain below 0°C to

qualify as permafrost. Soils that freeze during an exceptionally severe winter and survive for 1 or 2 more

years are called “pereletoks,” and often are not classified as permafrost (Swinzow 1969).

The existence of permafrost is a result of the history and the present state of the energy balance at the

Earth’s surface—measured by the surface temperature—and the deep Earth heat flow. If permafrost exists

and the net yearly gain of energy by the entire permafrost volume is equal to the net loss of energy, then the

permafrost will remain stationary, while an excess heat gain over heat loss will result in a net loss of frozen

material. Given the same energy balances, however, i.e., net gain of energy over the year, one region may

have permafrost (albeit degrading) while another will not. This is ascribable to the thermal history of the

frozen ground in the two areas. Though both are losing or have lost permafrost, one region may have started with a larger volume of permafrost than the other. Thus, the present energy flow conditions may be such

that permafrost cannot exist in one region, whereas it will subsist in another area, although in a receding

form often referred to as “relic permafrost.”

Permafrost Formation Time

VIRGIL J. LUNARDINI
In this sense previous glaciation has very likely played an important role in the present existence of permafrost in marginal areas. It is safe to conclude that little, if any, permafrost exists beneath nonpolar glaciers, but once they withdraw, permafrost may rapidly form and grow. So, previously glaciated regions will

show a lesser volume of frozen ground than unglaciated regions with similar climatic histories. In this regard

it is significant that Canada was heavily glaciated while Russia (Siberia) had little permanent glaciation.

Thus, the permafrost thickness in Siberia is much greater than in Canada, although the climates are similar.

Since the energy balance involves meteorological conditions, surface vegetation, topography and soil

conditions, we can anticipate that we will find no simple correlation for the existence of permafrost in the

marginal or discontinuous regions. Permafrost does exist far south of the usually accepted limits in scattered

patches almost always associated with high altitudes and, thus, microclimates similar to the usual permafrost

regions.

Origin and existence of permafrost

In discussing the present distribution of permafrost, questions often arise concerning its origin and age.

These two concepts should be clearly differentiated, since they deal with two separate phenomena. The age

of a particular deposit of permafrost is the time that has elapsed since the freezing of the soil system. Actually, it may be very difficult or impossible to determine this age because thawing and freezing may cycle at

long intervals and different frequencies in different regions of the Earth. Thus, the ages of two “similar” deposits of permafrost may be quite different. In this regard, the presence of preserved animal remains may be

a reliable clue to the age of a deposit of frozen ground. The age of permafrost is a question of significance

and may be useful to paleontologists, paleobotanists, etc. The present thermal state of the permafrost—temperature, degradation, aggradation, etc.—is of interest to engineers.

The origin of permafrost involves the question of the conditions under which it can form and grow. These

same conditions will explain the present existence of permafrost at a given location. As the conditions for

the origin of permafrost are dynamic, it is certainly possible that areas now lacking permafrost once had

these underlying frozen strata, and that the present regions of permafrost could once have been thawed. In

other words, the present existence of permafrost depends upon two things: the proper energy exchange conditions and the thermodynamic state of the permafrost mass itself. The first of these conditions has little to

do with past climatic conditions, but the second is a function of the complete thermal history of the permafrost and is thus related to past climate. In this sense, it is incorrect to describe permafrost simply as a legacy

of the last great ice age. It is possible that some, perhaps most, permafrost had its origin at the beginning of

the Pleistocene era (Brown 1964), but this should not imply that the intervening thermal conditions were

without significance.

From the above discussion, it should be clear that the formation and existence of permafrost are related to

the present and historical conditions of energy exchange between the soil and the atmosphere. Nevertheless,

it is not possible to state these conditions in a simple, precise manner that will allow us to define unique permafrost indicators.

Over a sufficiently long time span, the energy exchange will be periodic, and, averaged over a number of

periods, the net energy flow for a given soil volume will determine its thermodynamic state. Dynamic equilibrium of the energy flows may exist such that the soil is perennially frozen or the thermodynamic state of

the soil may be varying because of an imbalance of the cyclic energy flows. The original formation of permafrost depends upon a net periodic (yearly) loss of energy from the soil volume that must persist for many,

many years for the permafrost to attain great thicknesses. The maintenance of the present thermodynamic

state of permafrost requires only that the net energy flow averaged over a number of years be zero.

Paleotemperatures

The thermal history of permafrost is greatly influenced by the long-term temperature variations experienced at its upper surface. The relative mean global temperature deduced from oxygen isotope data is shown

for the past 180 million years in Figure 1 (Eddy and Bradley 1991). There was probably little or no permafrost prior to the late Tertiary Period and certainly none for 100 million years prior to the long-term cooling

2
that began 35 million years before present. During the Pliocene, some 2–5 million years before present, the

temperature oscillated 0–4°C above the present values. Permafrost was probably present at locations that

now have mean temperatures less than –14°C but with greatly reduced thickness and quite variable temporal

existence.

The Quaternary Period, which includes the Pleistocene (about 1 million years) and the Holocene (present

to about 10,000 years), remains a time of greatly reduced temperatures and massive glaciations. The period

is marked by ice ages of 100,000 to 120,000 years duration, interrupted by interglacials, and with increasingly severe minimum temperatures and temperature drops. We are now in the Holocene interglacial (which

is actually cooler than the previous interglacial, the Sangamon or Eem) and the record would seem to indicate the next significant temperature move should be downward with a new ice age evolving. Figure 2

shows the temperature record (again from isotope analysis) for the past million years or so (Folland et al.

1990). We note graphically the record of at least four or five glacial periods, with the last one ending some

12,000 years ago in North America—the Wisconsin. The temperatures variations from present values swung

from highs of +3°C during interglacials to lows of about –10°C during the glacial maximums of the last million years. The overall trends reveal very rapid temperature rises over time spans of approximately 12,000–

13,000 years, followed by less rapid temperature drops over 20,000 years, followed by a period of about

100,000 years of gradually decreasing temperatures with interglacial temperature rises of 4°C.

The temperature history of the past 160,000 years has been quantified using the deep glacial ice cores

taken in Greenland and Antarctica. Figure 3 shows the temperature variations from present values, for the

time of the last great ice age, taken from the isotope analysis of the Vostok (Antarctic) core (Jouzel et al.

1987). The maximum temperature excess of the past 160,000 years was about 3°C and thus only the discontinuous permafrost zones would have been in danger of disappearing, although there was continual change

in the permafrost thickness. If we assume that the temperature variations obtained from the ocean isotope

and ice core analyses are global, then the paleotemperature history at a specific site can be reconstructed using the present temperature data. This has been done for Prudhoe Bay by Osterkamp and Gosink (1991) as

shown in Figure 4 and for East Siberia by Maximova and Romanovsky (1988) (Fig. 5). Sun and Li (1988)

also presented a quantitative model of temperature fluctuations during the last ice age in northern China

(Fig. 6).

The isotopic data all indicate that the temperature at a given site could have dropped as much as 10 to Mean Global Temperature Warm Cold 180 150 135 100 65 58 38 25 13 0 1.8

Time Before Present (millions of years)

Mesozoic Tertiary Quaternary

Figure 1. Paleotemperature variations (after Eddy and Bradley 1991).

Figure 2. Pleistocene temperature fluctuations (after Folland et al. 1990).

3
b. Based on East Siberia model.

c. Based on Brigham and Miller (1983) data

for Barrow, Alaska.

a. Based on Vostok (Antarctica) ice core

temperature inferences.

Figure 4. Permafrost surface paleotemperature model for Prudhoe Bay, Alaska (after Osterkamp

and Gosink 1991).

Figure 3. Vostok (Antarctica) ice core temperature inferences (after Jouzel et al. 1987)

4
Time Before Present (thousands of years) 

0

2

– 2

– 4

– 6

– 8

– 10

SurfaceTemperature (°C)

0 20 40 60 80 100 120 140 160

Figure 6. Reconstructed temperature

for northern China during last ice age

(after Sun and Li 1988).

Glacial

0 20 40 60 80 100 120

Years Before Present (thousands of years)

2

0

– 2

– 4

 Surface Temperature

Below Present Value (°C)

1.0

– 2.76

Inter

Glacial

Glacial

Inter

Glacial

0 20 40 60 80 100 120

Time Before Present (thousands of years)

0

– 2

– 4

– 6

Surface Temperature Below Present (°C)

– 1.2

– 6.3

a. Based on Vostok (Antarctica) ice core data,

∆T = –5.29°C.

b. Based on Brigham and Miller (1983) data,

∆T = –2.54°C.

Figure 7. Mean paleotemperature departure, Prudhoe Bay, Alaska, during one glacial cycle.

12°C below the present values for varying periods of time. While extreme temperature drops at some sites

may have been significantly greater than these values, there is no convincing evidence of this. Temperature

variations of 8 to 12°C during the last glaciation have also been reported from the Greenland ice cores

(Dansgaard and Oeschger 1989). Folland et al. (1990) note that global temperatures underwent 5–7°C variations, with changes as great as 10–15°C at middle and high latitude regions of the Northern Hemisphere.

During the Eemian interglacial, temperatures in Siberia, Canada and Greenland may have increased by 4–

8°C. Thus, it would seem prudent to use extremal paleotemperature excursions on this order of magnitude

(10–12°C) for the development of permafrost formation models. We cannot ascribe the rapid growth of

deep permafrost to extraordinarily low temperatures that are beyond the ranges we have mentioned above.

Figure 7 shows two examples of average paleotemperatures for a glacial cycle made up of glacial and inter5

12 8 4

Lower Than Present

Temperature (°C)

0

20

40

60

80

100

120

140

Time Before Present (thousands of years)

Figure 5. Permafrost surface paleotemperature model for East Siberian location

(after Maximova and Romanovsky 1988).
glacial intervals. Figure 7a is based on the Vostok ice core while Figure 7b is calculated from information

given by Brigham and Miller (1983) for Prudhoe Bay, Alaska.

There have been surprisingly few systematic thermal studies of the origin of permafrost and the total

time required for its formation. Osterkamp and Gosink (1991) studied the response of permafrost thickness

to surface temperature variations. They used quasi-steady and numerical models to predict the position of

the permafrost bottom, with arbitrary initial permafrost thicknesses, but were interested in the inverse problem of deducing paleotemperatures from the present permafrost data at Prudhoe Bay, Alaska. Allen et al.

(1988) used a quasi-steady model (which will inherently underpredict the growth time) for the same purpose in the Mackenzie Delta region of Canada, using still another paleotemperature history (Fig. 8). They

assumed that Illinois permafrost was completely melted during the Sangamon interglacial, although this is

highly unlikely given the time available, unless Illinoisian glaciation limited the permafrost thickness to

modest values. Nevertheless, they used an initial permafrost thickness and predicted about 800 m of growth

in 75,000 years. Their thermal model is based on the work of Lachenbruch et al. (1982) dealing with the

rate of thaw of thick permafrost zones during several millennia. Romanovsky et al. (1988) noted qualitative

aspects of the origin and disappearance of permafrost in the Transbaikal area of Russia. Katasonov (1988)

used “cryogenic structures” to argue for the origin of permafrost early in the Quaternary and its persistence

up to the present.

In this report an attempt is made to predict the rate of permafrost formation, starting with no permafrost,

i.e., its origin, using a simple conduction model. If the soil forming the permafrost exists before freezing

starts, the growth is heterogenetic. When the permafrost forms as the soil material is gradually deposited at

the surface, the permafrost is said to have a syngenetic origin. The thermal conditions for each type of

growth will be examined.

THEORY

The solution to conductive heat transfer problems, with solidification phase change, has interested engineers and mathematicians for over a century. These problems (often referred to as Stefan problems) are inherently nonlinear and solution methods are very restricted. A classical solution for the case of a constant

temperature, semi-infinite medium that undergoes a step change of surface temperature was given by Neumann (1860) and expanded upon by Carslaw and Jaeger (1959); it is called the Neumann solution. Tao

(1978) extended the similarity technique of Neumann to the semi-infinite slab with arbitrary initial temperature. This is precisely the solution we need, but unfortunately this exact solution is such that numerical

computations are extremely difficult because of transient functions that require an increasing number of series terms as time increases. For permafrost formation the time scales are so huge that Tao’s solution is impractical and cannot be used. Like the exact solution of Lozano and Reemsten (1981), for flux boundary

conditions, Tao’s solution is perhaps best used numerically to verify the accuracy of approximate and numerical solutions or for short-time solutions.

The search for practical solutions for engineering design has led to some convenient approximate solution methods for Stefan problems. The heat balance integral technique solves the energy equation on averFigure 8. Paleotemperature history at Mackenzie Delta, Canada (after Allen et al. 1988).

6
age over a space volume, instead of at each point of space (Goodman 1958). The concept is identical to the

well-known momentum integral technique of fluid mechanics, sometimes referred as the Karman-Pohlhausen method (Schlichting 1968). The method has been applied successfully to constant initial temperature

problems of the semi-infinite slab (Lunardini and Varotta 1981) as well as the cylindrical geometry (Lunardini 1980). A modification of the integral method utilizing a single integration over an entire nonconstant

property volume has yielded accurate solutions (Yuen 1980, Lunardini 1981b,1982).

The integral solution has been used for a problem with variable initial temperature distributions, but the

results were limited to shallow freeze depths (Lunardini 1984). This report presents an approximate solution

to the modified Neumann problem for which a linear initial temperature distribution exists. Such an initial

temperature is common for soil systems with a geothermal temperature gradient and is directly applicable to

the question of permafrost formation rates.

Heterogenetic freeze relations

Figure 9 shows the case of an infinite layer of soil with a linear initial temperature distribution (G represents a geothermal gradient). The soil is assumed to be homogeneous and conduction is the only mode of

energy transfer. At zero time the surface temperature drops to Ts and is held constant while freezing commences. At any time t > 0, there is a frozen layer, called layer 1 (0 < x < X) and a thawed layer (x > X). The

thawed layer is further divided into layer 2 (X < x < X + δ) where temperature changes occur and layer 3 (x

> X + δ) where thermal effects are not discernible. We ignore the finite time it takes for the surface temperature to drop to the freezing point, Tf. This time will be small compared to the formation time and a realistic

scenario prior to the onset of freezing is To = Tf .

The governing equations are the conduction energy equations with appropriate boundary and initial conditions; see the Nomenclature for definition of symbols not defined in the text.

For the frozen zone

α

∂

∂

∂

∂ 1

2

1 1

2 0 T

x

T

t

= ≤≤ x X (1)

T1 (X, t) = Tf (1a)

T1 (0, t) = Ts . (1b)

For the thawed zone

α

∂

∂

∂

∂ δ 2

2

2 2

2 0 T

x

T

t

= ≤≤ + x X (2)

Figure 9. Freeze of a semi-infinite region with linear initial temperature.

7

Frozen Thawed

12 3

Thawed

X δ

x

To

Tf

Ts

Ti = Gx + To
8

T2 (X, t) = Tf

∂ δ

∂

TX t

x

G 2 ( ) + = , . (2b)

The initial temperature at the beginning of freeze is

Ti = To + Gx . (2c)

The maximum depth at any time to which the temperature disturbance will be felt is X + δ. Then

T2 (X + δ, t) = (X + δ)G + To. (2d)

The energy balance at the phase change interface for the freeze process is

 

k T

x Xt k T

x

X t dX

dt 1

1

2

2

2

∂ ∂

∂ ρ ∂ ( ) , ,. − ( ) = l (2e)

The energy balance at the freezing front can also be written as two equations (Lunardini 1981b)

 

−  ( )



 



 + () () = ( ) k T Xt

x

k T Xt

x

T Xt

x

T Xt

x

1

1

2

2

2 1

1 1

2

1

2

∂

∂

∂

∂

∂

∂ ρ α

∂

∂

, ,, , l (3)

 

− () () +

 ( )



 



 = ( ) k T Xt

x

T Xt

x

k T Xt

x

T Xt

x

1

1 2

2

2

2

2 2

2

2

2

∂

∂

∂

∂

∂

∂ ρ α

∂

∂

,, , , l . (4)

Because of the initial temperature distribution, during freeze the heat flow to the interface from the

thawed region will exceed the geothermal heat flow until equilibrium is established. Likewise, during a

thaw period the heat flow from the thawed zone will be less than the deep geothermal heat flow.

An approximate solution to this problem will be obtained using the heat balance integral technique (see

Lunardini 1991). In this method, the differential equations are solved on average over a finite volume of

material rather than at each point of the region. The integration of the energy equations over the regions

where temperature changes are occurring, 0 ≤ x ≤ X + δ, detailed by Lunardini (1981b) is

 

d

dt

c T x t dx c T x t dx X c c T X

cX T G X k T t

x

k G

X

X

X

ρ ρ ρ ρρ

ρδ δ ∂

∂

δ

11 1 2 2 2 1 2 2 11

22 1

1

2

0

2

0

∫ ∫ ( ) + ( ) −+ − ( )











− + ( ) + + ( ) 



 











= − ( ) +

+

, ,

, .

l f

o (5)

Quadratic temperature profiles in regions 1 and 2 that satisfy the boundary conditions are chosen as

T T aX x X

X

aX T x X

X 1f1 1 1 = +  −





 + − ( ) −





 ∆

2

(6)

TT G X T x X GX T x X

2 f =+ + [ ] ( ) + − − + ( ) ( ) − δ

δ δ

2 2 2

2

∆ ∆ (7)

where
9

a X T

g

g

T GX X

GXT 1

21 = = ( ) +

[ ] ( ) + +

+

∆ ∆

∆

1

2 2 , . 1 αδ δ

In general, the simplest temperature profiles that will satisfy the boundary conditions should be chosen.

The accuracy of the method increases as the order of a polynomial temperature choice increases; however,

the use of high-order polynomials (third and higher) is often not justified since a small increase in accuracy

requires significantly more computational effort. Equation 4 can be used to find a relation between X and d.

In nondimensional form this is

 

β σβ φ ρ β

g

k g

S

− + [ ] ( ) + = ( ) −

21

21 2 2

2 1

T

. (8)

Equation 5, the energy integral equation, can now be written nondimensionally, using eq 6 and 7 as

τ σ

σ

= ∫ Kd

0

(9)

K

b b

g

g

g

C C

g

k

=

+− −  ′



 



 − + 





 − + ( ) ′

− 



 



 +

1 2 21

21

21

1

6

1 2

3

1

3

1 1 2

β σ σ β σ φ σβ

σ

(10)

where

 

b

S

C

T

1 21

1

3

1

=− + + 



 



φ

b C 2 21

1

3 =− φ .

The derivatives of β and g can be found from the following equations

d

d

a aa

a aa

β

σ = β′ = ++

5 1 4

3 24

. (11)

dg

d gaaσ = ′ = − β′ 1 2 . (12)

where

a

m m

1

21 1

2 = −  ( ) + ( ) +



 



 α σ φ ββ

a m 2

21

2 = 2 12 ( ) + [ ] ( ) + + α σφ σβ φ

a

g

g

S

k 3

21

21

1 2 1 = − ( ) − − ρ σ

T

a

S g

4

21

2

2 1

= + 



 



 ρ β

T
a5 = k21(β + 2)

m = β[σ(β + 2) + 2φ].

The problem has now been reduced to a simple numerical quadrature of eq 9 using the auxiliary relations

of eq 10–12. The numerical solution of eq 9 can be obtained quite easily with a personal computer and any

standard numerical integration routine. (A FORTRAN program, PERM.FOR, to carry out the integration is

listed in Appendix E.)

The model requires only the ratios of the thawed to frozen values of thermal conductivity, specific heat capacity and density for the permafrost soils. These property ratios can be estimated with good accuracy for soil

systems as noted in Appendix A. The absolute values of the frozen and thawed soil properties are not needed

to carry out the solution of eq 9–12.

Heterogenetic model verification

It is possible to check the solution for a special case. Although there is no exact solution for the phase-change case, a

solution was found for the transient location of the Tf isotherm

for a homogeneous soil with zero latent heat, i.e, infinite Stefan number (Lunardini, in prep.). The relation is

erf

f s

X

t

GX

2 T T

1

α = − − . (13)

If we let the Stefan number be large and hold the property

ratios to unity, the approximate solution can be compared to this exact relation. Table 1 notes the results for a

typical case.

The comparison indicates that the approximate technique gives good results, especially as the time increases. The results also show that the Tf isotherm requires surprisingly long times to penetrate deeply, even

without phase change. This is explainable by the sensible-to-latent heat ratios to be examined below.

Sensible and latent heats

The total energy extracted from a unit area of soil is the sum of the latent and sensible energies.

QT = QL + QS . (14)

The latent energy is

QL = XL (15)

while the sensible heat flow is

Q C T T x t dx C T T dx C T T x t

X X

X

X

S f f1 u if u i2 = − [ ] ( ) + − ( ) + − [ ] ( ) ∫ ∫∫+

, ,

0 0

δ

dx (16)

where Ti (x) is the original temperature before the freeze starts. Using the temperature relations (eq 6 and 7)

leads to

QS = Cf ∆T1 X{C21 [(σ/2) + φ + β(σ+ φ) / 3] + (0.5/g + 1) / 3} . (17)

The ratio of the sensible to the latent heat is

QS /QL = ST {C21 [(σ/2) + φ + β(σ + φ) / 3] + (0.5/g + 1)/3}. (18)

This ratio is quite large even for small Stefan numbers and tends to increase as the freeze depth increases.

Table 1. Movement of Tf isotherm, homogeneous soil for heterogenetic freezing.

ST =1000, l = 0, Tf – Ts = 10°C,

φ = 0, G = 0.0286°C/m.

X Time (years)

(m) Eq 9 Exact eq 13 % difference

27.97 2.08 2.15 3.6

69.93 22.63 25.30 10.6

139.86 204.48 234.73 12.9

314.69 53,053.2 48,141.3 –10.2

332.17 242,852.4 238,469.1 –1.8

10
Figure 10. Temperature distribution after

initial freeze of soil.

11

Neumann solution

The Neumann solution is a special case that occurs when the initial temperature gradient G is zero.

Thus, it is the conduction solution that will predict the minimum time needed to form a given thickness of

frozen material. The exact solution is well known as

X t = 2γ α N 1 . (19)

γN is a function of ST , φ and the property ratios and is known (Lunardini 1981). Equation 19 can also be

written as

σ γτ N N = 2 . (20)

Steady-state solution

Unlike the Neumann solution, the permafrost zone for the general case with a non-zero geothermal gradient reaches an equilibrium value. The equilibrium value is the thickness of permafrost that will form if infinite time is available for growth. At steady state, the net heat flux at the phase change interface will be

zero so that (dX/dt)∞ = 0. Then the temperature in region 1 is

T T x

X 1 T 1 ∞ = + s ∞

∆ . (21)

At the solidification interface

k

T

x

1 2 k G ∂

∂

1∞ = . (22)

Thus

X T

k G ∞ = ∆ 1

21

(23)

or

σ∞ = 1

21 k . (24)

Effect of quaternary freeze–thaw cycles

The previous discussion assumed that the soil was initially unfrozen, with a temperature distribution given by eq 2c, as shown in Figure 9. However, the cyclic warming and cooling that has taken place over the

past million years will tend to cyclically lower the ground temperature as compared to the previous initial

temperature distribution. If the temperature is lower than that of eq 2c, the permafrost growth will be accelerated. Figure 10 shows the thermal state of permafrost at equilibrium with a geothermal heat flow—the

solid line. The approach to the equilibrium state of the permafrost is affected by the initial temperature of the soil in

that the time to reach equilibrium will be greatly changed

but the final equilibrium thickness is not affected by the

initial assumption (see eq 23).

Let us assume that surface warming to Tf occurs and

the permafrost eventually reaches a stage where the entire

frozen zone is at Tf . The system starts to thaw in this condition and we ignore the time needed to reach this state

(small compared to melt time scales) and also neglect the

bottom freeze–thaw during this initial warming time. Figure 11 shows the thermal state as the bottom thaw commences, at x = Xo. The geothermal heat flow is used to

melt permafrost and also to warm the soil in the region X
Figure 11. Bottom thaw of permafrost. Figure 12. Thermal state at end of melting (a) and

long-term equilibrium (b).

Figure 13. Geometry and coordinate system for freeze of a

semi-infinite medium with moving upper surface.

+ δ. The maximum thaw rate occurs if all of the geothermal heat flow goes into melting, but this is physically impossible.

Given sufficient time, the entire permafrost volume melts and the soil temperature is as shown in Figure

12—curve a. A layer of soil of thickness Xo + δo will be thermally modulated. The temperature gradient for

depths greater than Xo+ δo is G, the geothermal gradient. Also shown in Figure 12 is the dashed line (curve

b) denoting the equilibrium temperature distribution with the same geothermal gradient but no cooling.

This state would be reached if the surface remained at Tf for a long time after thaw was completed. Clearly,

the sensible heat that must be removed in the modulated case (curve a) will be less than that for the original

freeze discussed earlier (curve b). This means that cooled soil (after cyclic melting or warming) will freeze

much faster with the same mean surface temperature. This temperature modulation will only be significant

for zones of relatively shallow permafrost, for which appreciable thaw can take place during the interglacials of about 15,000 years. The equations and solution are discussed in Appendix B.

Syngenetic growth of permafrost

Syngenetic growth of permafrost occurs when material is deposited at the

surface while freezing is in progress.

This is inherently a much more efficient

freeze process and the growth of frozen

layers can be greatly accelerated. Figure

13 shows a sketch of the process with the

surface deposit laid down such that the

surface, held at a constant temperature

Ts, is moving at a constant velocity U.

The total frozen zone at any time is equal

to the material frozen at the interface

X(t), plus the depositional layer Ut. The

syngenetic system will be inherently unstable with no equilibrium solution. If

the frozen zone should equal the heterogenetic steady-state value (U = 0), the

motion of the upper surface will cause

melting at the base of the permafrost.

12
13

Thus, at some time during the formation process, the bottom must be melting to compensate for the depositional growth but the process cannot be steady.

Let us fix the upper surface so that it remains stationary, as shown in Figure 14. Then it would appear

that a steady flow of material is moving at a constant velocity U and the original soil surface seems to be

moving downward at a steady velocity. The energy equation is transformed such that a convective term exists. The equations for regions 1 and 2 are

α

∂

∂

∂

∂

∂

∂ i

i ii

2

2 0 12 T

x

U T

x

T

t − −= =i , . (25)

The boundary and initial conditions are exactly the same as those of the heterogenetic case (see Appendix

C for details of the syngenetic equations).

The heat balance integral form of the energy equations is

 

d

dt

c T x t dx c T x t dx c c T X

cX T G X k T t

x

k G cU T c

X

X

X

ρ ρ ρ ρρ

ρδ δ ∂

∂ ρ ρ

δ

11 1 2 2 2 1 2 2 11 f

22 o 1

1

2 11 1 2 2

{ ( ) + ( ) − − X + ( )

− + ( ) + + ( ) 



 











= − ( ) +− −

∫ ∫

+

0

2

0

, ,

,

l

∆ U T GX [ ] ∆+ + ( ) δ

(26)

We note that this is identical to the relation for heterogenetic growth, eq 5, except for the two additional

terms on the right-hand side of eq 26. Carrying out the integrations and making eq 26 nondimensional,

leads to the following result.

τ σ

σ

= ∫ K d1

0

(27)

K

b b

g

g

g

C C

g

k SC

1

1 2 21

21

21 T 21

3 =

+− −  ′



 



 − + 





 − + ( ) ′

− 



 



 +− + + + + { } [ ] ( )

β σ σ β σ φ σβ

σ

ψ φ σβ

1

6

1 2

3

1

1 1 2 11 1 /

(28)

where

.

Figure 14. Coordinate system for stationary upper surface.
14

ψ

α = U TG

∆ 1

1

.

Now the integration of eq 27 follows exactly as was done previously (see FORTRAN program

PFTSYNB.FOR in Appendix E). The model requires only the ratios of the thawed to frozen values of thermal conductivity, specific heat capacity, and density for the permafrost soils, as was noted earlier.

Syngenetic model verification

It is possible to check the solution for a special case as was done for the heterogenetic growth. Although

there is no exact solution for the phase-change case, there is an exact solution for the transient location of

the Tf isotherm for the same problem with a homogeneous soil with zero latent heat, i.e, infinite Stefan

number (Lunardini, in prep.). The relation is

e B B A A A ψσ

τ τ

f erfc

2

erfc

2 ( ) −1 1 20 − + ( ) + = . (29)

where A = σf – ψτ

B = σf + ψτ

σf = location of Tf isotherm.

If we let the Stefan number be large and hold the property ratios to unity, the heat balance integral solution

for syngenetic growth can be compared to this exact relation. Table 2 notes the results for typical cases.

The results indicate that the approximate technique gives excellent results, especially as the time increases. Thus, the Heat Balance Integral method and the numerical quadrature are robust even for very long

time spans. See also eq C17 for further verification of the solution method with phase change.

DISCUSSION

Equation 9 was solved numerically using Simpson’s rule. This resulted in values of the permafrost depth

versus time as a function of ST, ε and the thermal property ratios of the frozen and the thawed zones. The

results are presented in Figures 15–17.

These graphs depend only upon the quantities ST, φ and ε. In Appendix A it is shown that the soil porosity, ε, determines the saturated soil property ratios. The thermal property ratios used for the graphs are

listed in Tables A1 and A2. The graphs are only valid for the particular soil ratios given. However, this

does not affect the validity of the model itself. Any specific site can be modeled by using site-specific property ratios in eq 9. Figures 15–17 can be used to estimate permafrost formation times for a wide range of

surface temperatures and geothermal gradients.

The graphs can also be applied to variable surface temperatures with a bit of manipulation. Figure 18

Table 2. Movement of Tf isotherm, homogeneous soil, syngenetic freezing.

ST =1000, l = 0, Tf – Ts = 10°C, φ = 0, G =

0.0286°C/m.

τ σf exact σf approximate % difference

a. U = 10 mm/yr

0.0014 0.0087 0.090 –1.5

0.0069 0.1634 0.1700 –4.0

2.033 0.7739 0.7700 0.5

2.8232 0.8206 0.8100 1.3

b. U = 1 mm/yr

0.005 0.1448 0.15 –3.6

0.099 0.3874 0.4006 –3.42

1.0689 0.6572 0.6662 –1.37

11.4705 0.8795 0.8720 0.85
15

Figure 15. Formation time of permafrost, saturated mineral soils, φ = 0, ε = 0.4.

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0

10– 2 10– 1 100

τ

σ

ST = 0.60

0.40

0.20

0.10

0.06

0.04

0.02

0.01

ε = 0.379

φ = 0

Figure 16. Formation time of permafrost, saturated mineral soils, φ = 0, ε = 0.379 (Prudhoe Bay, Alaska).

1.6

1.4

1.2

1.0

0.8

0.6

0.4

100 101 102 103 104

τ

ε = 0.379

φ = 0

ST = 0.60

0.40

0.30

0.20

0.15

0.10

0.08

0.06

0.04

0.02

0.01

plots the time needed to reach 90% of the equilibrium permafrost thickness at a site. Figure 19 shows the

thickness of permafrost formed after 15,000 years as a function of average surface temperature for different

soils as characterized by the soil porosity. Figures 20 and 21 show that syngenetic growth greatly reduces

the formation time during the latter stages of permafrost growth. However, there may be significant questions as to whether surface deposition has continued for long periods of time at any given location.

ST = 0.68

0.40

0.30

0.23

0.20

0.15

0.10

0.08

0.04

0.02

0.01

ε = 0.4

φ = 0

1.6

1.4

1.2

1.0

0.8

0.6

0.4

100 101 102 103 104

τ
16

σ

10– 1 100

τ

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0

10– 2

ST = 0.60

0.40

0.20

0.10

0.06

0.04

0.02

0.01

ε = 0.3

φ = 0

Figure 17. Formation time of permafrost, saturated mineral soils, φ = 0, ε = 0.3.

Time τ = α1

G

2

t ——

∆T1

80

120

160

200

240

280

0.1 0.2 0.3 0.4 0.5 0.6 0.7

Stefan Number 

C1 ∆T1 —

 L

ε = 0.400

0.379

0.200

0.300

1.6

1.4

1.2

1.0

0.8

0.6

0.4

0.2

100 101 102 103 104

τ

σ ST = 0.60

0.30

0.20

0.10

0.06

0.04

0.02

0.01

ε = 0.3

φ = 0

0 10 20 30

600

400

200

0

Permafrost Thickness (m)

∆T1, Temperature (°C)

ε = 0.300

0.379

0.400

Figure 18. Formation time for permafrost to

reach 90% equilibrium thickness, saturated

mineral soils.

Figure 19. Effect of soil porosity on permafrost

thickness after 15,000 years, saturated mineral

soils, G = 0.0286°C/m.
17

Figure 21. Syngenetic growth of permafrost, saturated mineral

soils, φ = 0, ε = 0.379, ψ = 0.01 (Prudhoe Bay, Alaska).

1.6

1.4

1.2

1.0

0.8

0.6

0.4

1 10 100 1000

τ

0.01

0.15

0.05

ST = 0.50

Syngenetic Growth

 = 0.01

 = 0.379

 = 0

ψ

ε

φ

Figure 20. Syngenetic growth of permafrost, saturated mineral

soils, φ = 0, ε = 0.379, ST = 0.15 (Prudhoe Bay, Alaska).

1.6

1.4

1.2

1.0

0.8

0.6

0.4

1 10 100 1000

τ

σ

0 0.01

0.005 ψ = 0.02

Syngenetic Growth

 = 0.15

 = 0.379

 = 0

ε

φ

ST

Table 3. Paleotemperature scenarios, Prudhoe Bay, Alaska

(after Osterkamp and Gosink 1991).

Duration before

Ts* present time Stefan

Scenario (°C) (years) number

Present (Lachenbruch et al. 1982) –10.99 10,000–15,000 0.1438

Brigham and Miller (1983) –13.69 240,000 0.1827

Matteucci (1989) –11 300,000 0.1440

Vostok, Robin (1983)† –16 160,000 0.2159

E. Siberia (Maximova and

Romanovsky 1988)† –11.3 170,000 0.1483

* Average value of the fluctuating surface temperature over the indicated

 time span.

† Fitted to present Prudhoe Bay surface temperature.
Prudhoe Bay, Alaska

Considerable information on the permafrost is available from oil wells in the Prudhoe Bay, Alaska, area

(Lachenbruch et al. 1982). Using the actual permafrost data, we note that the property ratios for ε = 0.379

(Tables A1 and A2) are very close to measured and estimated values (Lachenbruch et al. 1982). Some temperature possibilities for Prudhoe are listed in Table 3. Figures 16a and b can be used to estimate the permafrost formation time, depending upon the temperature chosen.

Example 1

At Prudhoe Bay, Alaska, the permafrost has the following present conditions:

Ts = –10.99°C

G = 0.0286°C/m

ε = 0.379

ST = 0.1440

α1 = 58.89 m2/yr

kg = 4.34 W/m K

measured k21 = 0.5795

measured present permafrost thickness, Xp = 599.3 m.

Then, the equilibrium permafrost thickness is

X T

k G e 1

21

= = m

∆ 602 8. .

The calculated equilibrium thickness is essentially the same as the measured value. How long would it take

to reach this depth if the surface temperature had been Ts = –10.99°C for an indefinite period? The nature

of the solution is such that the final equilibrium values will be reached only after very long times. From

Figure 18 we find the time to reach 90% of equilibrium (542.0 m) is

τ α = 



 

 1  =

1

G

T

t

∆

2

241 2.

thus t = 500,740 years. This time is obviously quite long and suggests that the present climate of Prudhoe

Bay is probably considerably warmer than it has been on average over the past glacial cycles. Such warming over the past 15,000 years is widely accepted.

Example 2

Prudhoe Bay, with the Brigham and Miller (1983) paleotemperature scenario, has the following: Ts = –

13.69°C for 225,000 years before warming to –11°C in the last 15,000 years. For this case Xe = 763.5 m.

For t = 225,000 years, τ = 67.3. From Figure 16, with ST = 0.1827, we read σ = 1.412. Thus, in 225,000

years the permafrost will grow to X = 626.5 m. This value will then slowly decay to the new equilibrium of

602.8 m over 15,000 years. This requires a melt rate of 1.58 mm/yr. This value is well within the estimated

average thaw rate of 2.5 mm/year for this case (see Appendix D). Note that a lower surface temperature

will greatly accelerate growth since the new equilibrium depth will be greater than before. Hence, a much

larger fraction of the growth will be during the early, rapid growth stage.

Figure 22 shows the permafrost thickness at Prudhoe Bay after six glacial cycles with some typical features of permafrost growth demonstrated. First, the initial permafrost growth is quite rapid, reaching a

thickness of 570 m after 120,000 years, with the paleotemperature model of Figure 7a. The thickness then

slowly approaches an equilibrium value of 739 m but it will surpass the present thickness of 600 m after

about 185,000 years. Thus, a paleotemperature model as cold as that of Figure 7b will yield permafrost that

is too thick. Also shown is the finite difference prediction of Osterkamp and Gosink (1991), using Figure 6

(which has the same mean temperature as Fig. 7a) and starting with 600 m of permafrost, that indicates

much faster permafrost growth and thicker permafrost. Their quasi-steady model neglects sensible heat, as18
Figure 22. Growth of permafrost at Prudoe Bay, Alaska; effect of paleotemperature scenarios.

Time (yr)

800

600

400

200

0

Frozen Depth (m)

0 1 2 3 4 5 6 7 × 105

I II III IV V VI

Prudhoe Bay, Alaska

Permafrost Growth

Ts

 = –13.76 Glac.

 –10.00 IG

Ts

 = –11.00

Osterkamp and

Gosink (1991)

Glacial Cycle Ts (°C)

 I, II G –11.00

 I –11.00

 III, IV G –13.76

 I –10.00

 V G –9.00

 I –8.00

 VI G –13.76

 I –11.00

Present 

Permafrost 

Thickness

Table 4. Effect of previous cooling on permafrost growth time.* X = 541.0 m.

Time

Case G (years) Comment

No pre-cool 0.0286 486,100 uses eq 10

Pre-cool 0.0220 78,900 uses eq 10

Pre-cool 0.0286 105,400 uses eq B15

Syngenetic 0.0286 172,500 eq 27, U = 1 mm/yr

* Soil properties for Prudhoe Bay, ∆T1 = 10.0°C.

sumes fixed geothermal heat flow to the freezing interface, and starts with permafrost depth far less than equilibrium. All of these approximations overestimate the

freeze rate. The second curve of Figure 22 shows the significant effects if the paleotemperature model is modified

by only a small amount. The predicted permafrost thickness will reach present values after about 640,000 years

and will tend to oscillate about this value.

Example 3

Consider the effect of previous cooling at

Prudhoe Bay. If the equivalent geothermal gradient is G = 0.0220°C, then the time to reach

541 m is 78,940 years. The results of three possibilities are shown in Table 4. Note the very

large effect of previous cooling or cyclic thermal modulation.

Deep permafrost

Example 1

Consider the case of very thick permafrost,

on the order of 1600 m. Let the properties be those of Prudhoe Bay but ε = 0.4, ∆T1 = 29.27°C, G = 0.0286,

α1 = 57.99 m2/yr, then Xe =1813 m. The value of the surface temperature chosen is on the order of 12°C

less than present winter temperatures experienced in parts of Canada, Russia and Greenland, although it is

doubtful that such temperatures could have persisted for 1 million years. The value used illustrates the long

time needed to form deep permafrost by conduction alone. We will find the time required to form 90% of

the permafrost or 1632 m. The calculations are as before, with ST = 0.4 (see Fig. 15 for the case without

previous cooling). The results are given in Table 5. The formation time is very long, even with previous

cooling. This example used Prudhoe Bay properties and geothermal gradient, which could be significantly

different at a site in Siberia with deep permafrost.

Example 2

Consider the question of the maximum permafrost thickness that is probable. This will occur if the frozen thermal conductivity is large, the geothermal heat flow and latent heat are low, and the surface temperature is minimal (within the constraints discussed earlier). Let ε = 0.2, Ts = –23.5°C, qg = 0.042 W/m2, and

assume coarse-grained soil with kg = 5.86 W/m K. The results of the calculation are shown in Table 6.

Table 5. Formation time of deep permafrost.*

Xe = 1813 m.

Time

Case G (years) Comment

No pre-cool 0.0286 4,190,600 uses eq 10

Pre-cool 0.0220 488,900 uses eq 10

Pre-cool 0.0286 358,900 uses eq B15

Syngenetic 0.0286 784,100 eq 27, U = 1mm/yr

Neumann (φ = 0) 0.0 64,800 absolute minimum time

* Soil properties for Prudhoe Bay, ∆T1 = 29.3°C.

19
20

Table 6. Extreme predicted permafrost thickness.

X (m)

Time X (m) syngenetic

(years) heterogenetic U = 1 mm/yr

100,000 1488.4 1539.9

532,257 2000.0 2158.6

1,000,000 2132.0 2362.9

2,000,000 2254.5 2529.5

Infinity 2600.2 2600.2

The permafrost reaches a thickness of 2132 m after 1 million

years and a value of 2255 m after 2 million years. Balobaev et al.

(1978) note that the greatest permafrost thicknesses recorded are on

the East Siberian platform and present graphs with maximum permafrost thicknesses of 1500 m. These thicknesses are on the order of the

value predicted here and it is not likely permafrost much thicker than

this has ever existed, since the required time exceeds the plausible

time available. These extreme thicknesses are not in thermal equilibrium with the present surface temperatures and are slowly thawing.

CONCLUSIONS

The calculations and examples indicate that the growth of permafrost, with pure conduction heat transfer, is governed by the transient surface temperature, the geothermal heat flow and the soil thermal properties. Permafrost grows very rapidly for an initial phase and then asymptotically approaches a steady-state

value after time spans of immense length. Very thick permafrost may have required the total Quaternary

Period to form. It is likely that permafrost is not in equilibrium at most sites. The bottom growth and decay

of permafrost are so slow that accurate methods of detecting which is occurring (or if equilibrium exists)

are not available for the field. Permafrost less than 600 m can grow within 50,000 years, with surface temperatures only slightly lower than present values, but deeper permafrost depths require time scales of several ice ages and quite low temperatures to form.

LITERATURE CITED

Allen, D., F. Michel and A. Judge (1988) Paleoclimate and permafrost in the Mackenzie Delta. In Proceedings, 5th International Conference on Permafrost, 2–5 August, Trondheim, Norway. Trondheim, Norway:

Tapir Publishers, vol. 1, p. 33–38.

Balobaev, V.T., V.N. Devyatkin and I.M. Kutasov (1978) Contemporary geothermal conditions of the existence and development of permafrost. In Permafrost: The North American Contribution to the 2nd International Conference on Permafrost, 13–28 July, Yakutsk. Washington, D.C.: National Academy of Sciences,

p. 8–12.

Brigham, J.K. and G.H. Miller (1983) Paleotemperature estimates of the Alaskan Arctic Coastal Plain during the last 125,000 years. In Proceedings, 4th International Conference on Permafrost, 17–22 July, Fairbanks, Alaska. Washington, D.C.: National Academy Press, p. 80–85.

Brown, R.J.E. (1964) Permafrost investigations on the Mackenzie Highway in Alberta and Mackenzie District. Ottawa, Canada: National Research Council, NRC 7885.

Bryan, K. (1946) Cryopedology—The study of frozen ground and intensive frost action with suggestions

on nomenclature. American Journal of Science, 244: 622–642.

Carslaw, H.S. and C.J. Jaeger (1959) Conduction of Heat in Solids, 2nd edition. Oxford: Clarendon Press.

COHMAP Members (1988) Climatic changes of the last 18,000 Years: Observations and model simulations. Science, 241: 1043–1052.

Dansgaard, W. and H. Oeschger (1989) Past environmental long-term records from the Arctic. In The Environmental Record in Glaciers and Ice Sheets (H. Oeschger and C.C. Langway, Jr., Ed). New York: John

Wiley and Sons, p. 288–317.

Eddy, J.A. and R.S. Bradley (1991) Changes in time of the temperature of the earth. EarthQuest, 5(1).

Folland, C.K., T.R. Karl and K.YA. Vinnikov (1990) Observed climate variations and change. In Climate

Change, The IPCC Scientific Assessment, chapter 7 (J.T. Houghton, G.J. Jenkins and J.J.Ephraums, Ed.).

Cambridge: Cambridge University Press.

Gold, L.W. and A.H. Lachenbruch (1973) Thermal conditions in permafrost—A review of North America

literature. In Permafrost: The North American Contribution to the 2nd International Conference on Permafrost, 13–28 July, Yakutsk. Washington, D.C.: National Academy of Sciences, p. 3–25.

Goodman, T.R. (1958) The heat-balance integral and its application to problems involving a change of

phase. ASME Transactions, 80: 335–342.
21

Jouzel, J., C. Lorius, J.R. Petit, C. Genthon, N.I. Barkov, V.M. Kotlyakov and V.N. Petrov (1987) Vostok

ice core: A continuous isotope temperature record over the last climatic cycle (160,000 years). Nature, 329:

403–408.

Katasonov, E.M. (1988) Continuous persistence of the permafrost zone during the Quaternary Period. In

Proceedings, 5th International Conference on Permafrost, 2–5 August, Trondheim, Norway. Trondheim,

Norway: Tapir Publishers, vol. 1, p. 801–804.

Kersten, M.S. (1949) Thermal properties of soils. University of Minnesota Experimental Station, Bulletin

no. 28.

Lachenbruch, A.H., J.H. Sass, B.V. Marshall and T.H. Moses, Jr. (1982) Permafrost, heat flow and the

geothermal regime at Prudhoe Bay, Alaska. Journal of Geophysical Research, 87(B11): 9301–9316.

Lozano, C.J. and R. Reemsten (1981) On a Stefan problem with an emerging free boundary. Numerical

Heat Transfer, 4: 39–245.

Lunardini, V.J. (1980) Phase change around a circular pipe. USA Cold Regions Research and Engineering

Laboratory, CRREL Report 80–27.

Lunardini, V.J. (1981a) Heat Transfer in Cold Climates. New York: Van Nostrand Reinhold.

Lunardini, V.J. (1981b) Application of the heat balance integral to conduction phase change problems.

USA Cold Regions Research and Engineering Laboratory, CRREL Report 81–25.

Lunardini, V.J. (1982) Freezing of soil with surface convection. In Proceedings, 3rd International Symposium on Ground Freezing, 22–24 June, Hanover, N.H. USA Cold Regions Research and Engineering Laboratory, p. 205–213.

Lunardini, V.J. (1984) Freezing of a semi-infinite medium with initial temperature gradient. Journal of Energy Resources Technology, 106: 103–106.

Lunardini, V.J. (1991) Heat Transfer With Freezing and Thawing. Amsterdam: Elsevier (Amsterdam, Oxford, New York, Tokyo).

Lunardini, V.J. (in prep.) Permafrost and ground temperature calculations. USA Cold Regions Research

and Engineering Laboratory Report.

Lunardini, V.J. and R. Varotta (1981) Approximate solution to Neumann problem for soils systems.

ASME Journal of Energy Resources Technology, 103(1): 76–81.

Matteucci, G. (1989) Orbital forcing in a stochastic resonance model of the late-Pleistocene climatic variations. Climate Dynamics, 3: 179–190.

Maximova, L.N. and V.Y. Romanovsky (1988) A hypothesis for the Holocene permafrost evolution. In

Proceedings, 5th International Conference on Permafrost, 2–5 August, Trondheim, Norway. Trondheim,

Norway: Tapir Publishers, vol. 1, p. 102–106.

Muller, S.W. (1945) Permafrost or permanently frozen ground and related engineering problems. Special

Report, Strategic Engineering Study No. 62, U.S. Army Engineers, Military Intelligence Division Office.

Neumann, F. (c. 1860) Lectures given in 1860s, cf. Reimann-Weber, Die partiellen Differentialgleichungen. Physik, Edition 5, vol. 2, 1912, p. 121.

Osterkamp, T.E. and J.P. Gosink (1991) Variations in permafrost thickness in response to changes in paleoclimate. Journal of Geophysical Research, 96(B3): 4423–4434.

Robin, G. de Q. (1983) The Climate Record in Polar Ice Sheets. New York: Cambridge University Press.

Romanovsky, N.N., V.N. Zaitsev, S.Yu. Volchenkoc, V.P. Volkova and O.M. Lisitsina (1988) New data on

permafrost of Kodar-Chara-Udokan Region. In Proceedings, 5th International Conference on Permafrost,

2–5 August, Trondheim, Norway. Trondheim, Norway: Tapir Publishers, vol. 1, p. 233–236.

Schlichting, H. (1968) Boundary Layer Theory, 6th edition. New York: McGraw-Hill Book Co.

Sun, Jianzhong and Li, Xinguo (1988) Division and temperature condition of the last glaciation in northern

China. In Proceedings, 5th International Conference on Permafrost, 2–5 August, Trondheim, Norway.

Trondheim, Norway: Tapir Publishers, vol. 1, p. 107–112.

Swinzow, G.K. (1969) Certain aspects of engineering geology in permafrost. Journal of Engineering Geology, 3: 177–215.

Tao, L.N. (1978) The Stefan problem with arbitrary initial and boundary conditions. Quarterly of Applied

Mathematics, 36(3): 223–233.

Yuen, W.W. (1980) Application of the heat balance integral to melting problems with initial subcooling.

International Journal of Heat and Mass Transfer, 23: 1157–1160.
23

APPENDIX A: SOIL PROPERTIES AND RATIOS

The thermal conductivity of a mixture such as a soil can be estimated by using the weighted geometric

mean (Lachenbruch et al. 1982, Gold and Lachenbruch 1973, Lunardini 1981a). This can be written for a

general soil as

kk k k k x x x x = ( ) gw ai ( ) () ( ) g w a i (A1)

where ka, kg, ki and kw are the thermal conductivities of air, soil solids, ice and water; xa, xg, xw and xw are

the volumetric fractions of air, soil solids, ice and water. The geometric mean is usually better than the assumption of parallel geometry (weighted arithmetic mean), which is often used for simplicity.

Saturated soil

Many assumptions can be made concerning the soil saturation and porosity but simple approximations

will be used here. If the soil is always saturated, has a constant void ratio ε, and all of the water freezes,

then xw = ε and the conductivity ratio can be expressed as

k kk u wg = ( )( ) − ε ε γ

1

(A2)

k kk f ig = ( )( ) − ε 1 ε (A3)

kk k k k u f 21 w i / / = = ( )ε ε γ ( ) 1− (A4)

where γ = 0.9825 is a temperature correction for kg (Lachenbruch et al. 1982) and kw/ki = 1.34/5.45 =

0.2459.

The volumetric specific heat for the system may be expressed as follows, for the thawed and frozen

states

Cu = Csu (1 – xw) + Cw xw (A5)

Cf = Csf (1 – xi) + Ci xi (A6)

where Csu and Csf are the volumetric specific heats of unfrozen and frozen solids, and Cw and Ci are the

volumetric specific heats of water and ice.

It is fortunate that the volumetric specific heats of soil solids and ice are all about the same. For example, the specific heat of organic solids is 0.461 cal/cm3 °C, for mineral solids it is 0.420, and for ice it is

0.459 (Lunardini 1981a). If one assumes that the values for the solids, except for ice, change little through

the phase change then

Cf = 0.4202 + 0.0388ε (A7)

C

C

C u

f

= = 21

+

+

0 4296 0 5708

0 4202 0 0388

. .

. . . εε

(A8)

The density ratio is

ρ

ρ

ρ

ε

ε

u

f

= = 21

−

−

1 0 6154

1 0 650

.

. . (A9)
Finally, the ratio of thermal diffusivities is

α

α

α u

f

21

21

21

= = k

C . (A10)

The latent heat is

L = 79.71 ε (A11)

The Stefan Number is

 

S C T T C T

L T T 1

1 s

1 = − ( ) = = +

l

∆ ∆ 1

1

0 4202 0 0388

79 71

. .

. . εε

(A12)

It is possible to present the results for soil systems, quite efficiently, since the property ratios can be described as functions of the soil void ratio ε (Lunardini and Varotta 1981). Using the thermal conductivities

of Table A1, the property ratios used in the calculations are given in Table A2. The thermal conductivity

ratio will be representative of soil that is not too dry. Thus, eq A4 and A10 should be acceptable if ε ≥ 0.2

(Kersten 1949).

Balobaev et al. (1978) note that for limestone and dolomite, k = 2.44 – 3.37 W/m °C and qg = 0.017

– 0.021 W/m2 at 60–800 m, anomalously low heat flow values.

Nonsaturated soil

For the nonsaturated soil, assuming that the porosity does not change during phase changes, the ratio of

thawed to frozen thermal conductivity

k k kk S S

u g wa = ( ) − − ( ) ( ) γ ε ε 1 1 (A13)

 

kk k k S S

f g i a wi wi = ( )( ) − ( ) − 1 1 ε ε ρ ε ρ (A14)

k

k

kk k

S u

f

wi a wi wi S = [ ] ( ) ( ) − ( ) − γ ε ρ ε 1 ε ρ 1 / (A15)

where S is the thawed soil saturation level, ρwi = ρw /ρi; ρi/ρw = 0.91. Interestingly, the ratio k21 can have

the same values as for the saturated case if the saturation has certain values, e.g., S = 0.756, ε = 0.379.

24

Table A2. Calculated saturated granular

soil thaw–freeze property ratios.

Soil ku/kf Cu/Cf αu/αf ρu/ρf

porosity ε eq 21 eq 24 eq 26 eq 25

0.2 0.7448 1.2706 0.5862 1.008

0.3 0.6484 1.3909 0.4662 1.0129

0.379 0.5812 1.4847 0.3915 1.0174

0.4 0.5645 1.5094 0.3740 1.0187

0.5 0.4915 1.6265 0.3022 1.0256

Table A1. Thermal conductivity of

materials.

Thermal conductivity

Substance W/(m°C)

Water 0.561

ice 2.281

air 0.0237

silicaceous soil solids* 4.29–5.87

*Lachenbruch et al. (1982)
25

Figure B1. Bottom melt process.





x = 0

Changing

Temp.

T2 (x, t)

q1 qg

(Xo – X) X δ

x

Xo

(Xo + δ)

Tf

T °

s

Frozen Thawed

1 2

APPENDIX B: QUATERNARY CYCLIC THERMAL MODULATION

Thaw process

Assuming that melting occurs with a fraction f of the geothermal heat flow going into increasing the sensible heat allows an estimate to be made of the thaw time and the final temperature, noted in Figure 12,

which will then be the initial temperature distribution for the next freeze cycle. Referring to Figure 11, we

assume that the fraction of geothermal energy available for melt at any time is

f f XX

X = ( ) m o − +

o

1

(B1)

where fm is the melt fraction at the conclusion of melt and Xo is the initial frozen thickness. The temperature in the region of changing temperature is

T2 = Tf + b(x – ∆1) + c(x –∆1)2 (B2)

where ∆1 = Xo – X and b,c are functions of time only. The initial temperature of the soil is

T

T

T Gx X

x X

x X i

f

f o

o

o

= + −







≤

( ) ≥ . (B3)

The sensible heat added to the soil is

Q C T T dx C T T dx

X X

X

X

X

s u 2i u 2i

o

o

o

o

= − ( ) + − ( ) −

+

∫ ∫

δ

. (B4)

The sensible heat addition is shown as the shaded region of Figure

B1. The sensible heat at any instant is then

Q GC X X s u = − 







1

3 2

δ . (B5)

The change in sensible heat must equal the fraction of the geothermal energy not used for melting, or

dQ

dt q f S = − g ( ) 1 . (B6)

The energy balance at the melt interface is

L dX

dt = fqg . (B7)

The solutions to these equations lead to

1

3 2

1 GC X 1 X

q t M eMt

u g δ − 





 =+ − ( ) 



 



 (B8)

X

X e

f

Mt

= ( ) −

( ) −

o

m

1

1 (B9)
26

where

M f q

X L = ( ) m g −

o

1 .

The time to complete melt is

t XL f

q f m o m

g m

= ( ) −

ln . 1 (B10)

This leads to the value of δ when melt is completed.

δo o

u

m

m

= + − − 



 



 X L

GC

f

2 f

3

1

1 ln . (B11)

We define a linear temperature distribution that will have the same sensible heat when thaw is completed.

T T Gx ′ = + ′ i f (B12)

′

=

−

+

G

G

X

X

δ

δ

o o

o o

3 . (B13)

where G’ is the equivalent geothermal gradient. Finally, the new temperature distribution at the beginning

of freeze is

T T bx cx ifo o =+ + 2 (B14)

where

 

b G X

X

c

GX

X o

o o

o o

o o

o o

= ( ) −

+ =

( ) +

δ

δ δ , . 2

This initial temperature distribution is shown as curve a

in Figure 12. Table B1 shows some results for Prudhoe

Bay. Note the long melt times even if f is as high as

90%.

Freeze of cooled soil

The freezing process is as discussed earlier except

that the initial soil temperature is lowered as noted in

Figure 10. The basic equations used earlier are still valid

except that the coefficients of eq 6 and 7 change, owing to the new initial temperature given by eq B14.

The basic equation, replacing eq 9, is

dF

d

k

g

σ

τ σ =− − 



 

 21  1 2 1 (B15)

F

g S

C m S=− + + 



 



 − − −− 









 − − ( )( ) +

















1

6

1

3 6 63

1

2

1 1 21

2 2 ρ 2

σ β ββ σ β β

T

21 o o

3 (B16)

β σβ β σ β ρ β

g

km S g

S

− −+ + [ ] ( ) + + ( ) = ( ) − 21 o o

T

2 12 1 2 2 1 21 (B17)

Table B1. Melt relations for Prudhoe Bay.

δο Melt time

fm  ̃

f (m) (years)

0.1 0.55 1.5584 7948.1 0.9064 108,926

0.2 0.60 1.0118 5265.6 0.8636 85,654

0.3 0.65 0 .720 3833.5 0.8196 73,230

0.4 0.70 0.5272 2887.3 0.7706 65,022

0.5 0.75 0.3863 2195.8 0.7139 59,023

0.6 0.80 0.2771 1659.9 0.646 54,373

0.8 0.90 0.1157 867.8 0.455 47,502

Xo = 600 m, L = 30.21 cal/cm3, Cu = 0.6457 cal/cm3 °C.

qg = 1.35 × 10–6 cal/cm2 s, G = 0.0286°C.

ln m

m

f

f − −1

1

G′

G
27

g

Q

Q = − ( ) −

( ) − 1

2

α β 21

β β (B18)

where Q = m0 (β+1) + σSo (β+1)2

m0 = b0/G

S0 = c0∆T1/G2.

The solution to these equations follows in exactly the same way as previously discussed.
29

APPENDIX C: HEAT BALANCE INTEGRAL EQUATIONS

FOR SYNGENETIC PERMAFROST GROWTH

Method 1

We will formulate the basic equations in terms of a convective system with mass flowing through the

stationary upper surface at constant velocity U, as shown in Figure 14. The governing equations are

α

∂

∂

∂

∂

∂

∂ 1

2

1

2

1 1 0 0 T

x

U T

x

T

t

− − = ≤≤ x X (C1)

T Xt T 1 f ( ,) = (C1a)

TtT 1 s (,) 0 = (C1b)

α

∂

∂

∂

∂

∂

∂

δ 2

2

2

2

2 2 0 T

x

U T

x

T

t

− − = ≤≤ + XxX (C2)

T Xt T 2 f ( ,) = (C2a)

∂ δ

∂

TX t

x

G 2 ( ,) + = (C2b)

TX t X G T 2 ( ,) ( ) . + =+ + δ δ o (C2c)

The initial temperature at the beginning of freeze is

T T Gx i o = + . (C2d)

The energy balance at the phase change interface for the freeze process is

 

k T

x

Xt k T

x

Xt U dX

dt 1

1

2

2 ∂

∂

∂

∂ ( ,) ( ,) . − += ρ ρ 2 2 l l (C3)

Because of the initial temperature distribution, during freeze, the heat flow to the interface from the thawed

region will exceed the geothermal heat flow until equilibrium is established. Likewise, during a thaw period, the heat flow from the thawed zone will be less than the deep geothermal heat flow. The energy balance at the freezing front can also be written as two equations (Lunardini 1981b)

 

− 



 



 k + =

T Xt

x

k T Xt

x

T Xt

x

T Xt

x

1

1

2

21 1 ∂

∂

∂

∂

∂

∂ ρ α

∂

∂

( ,) ( ,) ( ,) ( ,) 2

1 1

2

2 l (C4)

 

− + 



 



 k = T Xt

x

T Xt

x

k T Xt

x

T Xt

x

1

1 2

2

2

2 2

2 ∂

∂

∂

∂

∂

∂ ρ α

∂

∂

( ,) ( ,) ( ,) ( ,) .

2 2

2 l (C5)

Quadratic temperature profiles in regions 1 and 2 that satisfy the boundary conditions are

T T aX x X

X

aX T x X

X 1 f = +  −





 + −  −





 1 11

2

( ) ∆ (C6)
TT G X T x X GX T x X

2 f =+ + [ ] ( ) +  −





 − + − δ

δ δ

2 2

2

2 ∆ ∆ ( ) ( ) (C7)

where

a X T

g

g

T GX X

GXT 1 = = ( ) +

[ ] ( ) + +

+

∆ ∆

∆

1 21

2 2 , . 1 αδ δ

In general, the simplest temperature profiles that will satisfy the boundary conditions should be chosen.

The accuracy of the method increases as the order of a polynomial temperature choice increases; however,

the use of high-order polynomials (third and higher) is often not justified since a small increase in accuracy

requires significantly more computational effort. Equation C5 can be used to find a relation between X and

δ. In nondimensional form this is

β σβ φ ρ β

g

k g

S

− + [ ] ( ) + = ( ) − 21

21 2 2

2 1

T

. (C8)

The heat balance integral forms for eq C1 and C2 are as follows.

α

∂

∂

∂

∂ 1

1 1

11 1

0 0 0

0

T Xt

x

T t

x

UT Xt T t d

dt

T x t dx T dx

dt

X

( ,) (,) − ,, , 



 



 − [ ] ( ) − ( ) − ( ) + = ∫ f (C9)

α

∂ δ

∂

∂

∂

δ

δ δ δ

2

2 2

2 2

2 0

TX t

x

T Xt

x

UT X t T Xt

d

dt

T x t dx T X

d X

dt

T dX

dt X

X

( ,) ( ,) , ,

, .

+ − 



 



 − + [ ] ( ) − ( ) −

( ) + + ( ) ( ) + − =

+

∫ 2 f

(C10)

These two energy equations are summed, along with eq C3, to yield an integrated equation for the entire

region undergoing temperature changes. The result is

 

d

dt

c T x t dx c T x t dx X c c T X c X T G X

k T t

x

k G cU T

X

X X

ρ ρ ρ ρρ ρ δ δ

∂

∂ ρ ρ

δ

11 1 2 2 2 1 2 2 11 f 2 2 o

2 11 2

, ,

,

( ) + ( ) −+ − ( ) − + ( ) + + ( ) 



 

















= − ( ) +− −

+

∫ ∫ l

2

0

0

1

1 ∆ 1 cU T GX 2 [ ] ∆+ + ( ) δ . (C11)

Equation C11, the energy integral equation, can now be written nondimensionally as

τ σ

σ

= ∫ K d1

0

(C12)

K

b b

g

g

g

C C

g

k SC

1

1 2 21

21

21 21

1

6

1 2

3

1

3

1 1 2 11 1

=

+− −  ′



 



 − + 





 − + ( ) ′

− 



 



 +− + + + + { } [ ] ( )

β σ σ β σ φ σβ

σ

ψ φ σβ / T

(C13)

where

ψ

α = U TG

∆ 1

1

30
b

S 1 21 C 1

3

1

=− + + 



 





T

φ

b C 2 21

1

3 =− φ.

The derivatives of β and g can be found from the following equations

d

d

a aa

a aa

β

σ = β′ = ++

5 1 4

3 24

(C14)

dg

d gaaσ = ′ = − β′ 1 2 (C15)

where

a m m 1

21 1

2 = −  ( ) + ( ) +



 



 α σ φββ

a m 2

21

2 = 2 12 ( ) + [ ] ( ) + + α σφ σβ φ

a g

g

S

k 3

21

21

1 2 1 = − ( ) − − ρ σ

T

a

S g

4

21

2

2 1

= + 









 ρ β

T

a k 5 = + 21( ) β 2

m = + βσβ φ [ ] ( ) 2 2. +

The problem has now been reduced to a simple numerical quadrature of eq C12 using the auxiliary relations of eq C13–C15. A FORTRAN program to carry out the integration is listed as PFTSYNB.FOR in Appendix E.

Phase change model verification

A simplification of this problem can be solved in a closed form. Consider the case of a soil initially

thawed at Tf and with a zero geothermal gradient G. The problem is then one of a single phase only with eq

C1, C1a,b, C3, C4, and C9 governing the freeze process. The temperature is chosen as

 

TT P x X

X

c P x X

X

= +  −





 −  −





 f 2

2

2

l (C16)

where

 

P

c

= =+ − RR S l , . 12 1 T

The location of the freeze interface is given by

t

K

K

X K

K

KX K

K = −  +



 





















1

2

3

2

2 3 ln

3

(C17)

31
32

K R

1 1 R

2

=+ + ( ) 1 3/

KU S 2 = + ( ) 1 T

K RR 3 = + α ( ) 1 .

Note that if U is zero, the phase change interface is

X R R

R R

t 2 1

1 21 3 = 2 ( ) +

+ + / / ( ) α . (C18)

This is identical to the well known Stefan solution given in Lunardini (1991). We may compare the closed

form solution (for which G = 0) with the numerical quadrature of eq C12 by letting G be very small. Table

C1 shows that the results are quite good even for very long freeze times.

Method 2

We can examine the same problem with a different approximation method by referring to Figure 13. For

region 3, a quasi-steady approach will be used, leading to a linear temperature profile. The basic equations

for heterogenetic growth are valid except that the surface temperature will be replaced by a transient function T t ′( ) s . Equations 1–5 are valid but the temperature profiles are changed as follows. Quadratic temperature profiles in regions 1 and 2 and a linear temperature in region 3 that satisfy the boundary conditions are

T T aX x X

X

aX T x X

X 11 1

2

= +  −





 + − ( ) −





 f 3 ∆ (C19)

TT G X T x X GX T x X

2

2

2 =+ + [ ] ( ) 2 2 + − − + ( ) ( ) − f δ

δ δ

∆ ∆ (C20)

 

T T T MRX x

X 3 1 =+ + 1 



 

 s 

d

∆ (C21)

where

a X T M

g

g

T GX X

GXT Rk g

M

R T T T t R TM k g d

1

1 21

13

3 1 13

2 2

1 21

1

1

2 1

= = ( ) +

[ ] ( ) + +

+=− ( )

= + = − ′( ) = − [ ] ( )

∆ ∆

∆

∆ ∆

, , / /,

, / /.

α

δ δ

σ σ

σ σ

d

f s

Equation 5 can be used to find a relation between X and δ. In nondimensional form this is

β σβ φ M ρ β

g

k g

S

− + [ ] ( ) + = ( ) − 21

21 2 2

2 1

T

. (C22)

Equation 3, the energy integral equation, can now be written nondimensionally as

τ σ

σ

= ∫ K Hd /

0

(C23)

K b b A M g A Mk g g g g

b A Mg g A

=+ − ( ) − ( ) −+ + d ( ) ( ) ′ − − ( ) [ ]

+ − [ ] ( ) ′ + ′ ( ) − + ( )

1 2

2

13

2

2

2

3 6 2 05 1 2 1 3

3 6 3 05

σβ σ σ σ

σ σβ β

/ / / ./ / / / /

/ / / .. (C24)

Table C1. Comparison of closed

solution (G = 0) and numerical

quadrature (G = 0.0001).

ST = 0.144, α = 58.89 m2/yr, ∆T1 =

10°C, U = 1 mm/yr.

Freeze Time Time

depth (yr) (yr) Percent

(m) eq C12 eq C17 difference

1000 55,867 54,778 1.99

2000 206,996 203,653 1.64

3000 437,935 428,448 2.20
33

H M= − ( ) g kM g k [ ] + + ( ) +

σ

1 2 1 05 1 3 ψ σ 13 12 / ./ / (C25)

where

AC b M

S

C bC

M

R

G

G

X

= =− + + 



 



 = −

= + =

21 1 21 2 21

1

3

1 1

3

1

1

, , ,

, .

t

d d

φ φ

σ

∆

The derivatives of β and g can be found from the following equations

d

d

a aa

a aa

β

σ = β′ = ++

5 1 4

3 24

(C26)

dg

d gaaσ = ′ = − β′ 1 2 (C27)

where

a m m 1

21 1

2 = −  ( ) + ( ) +



 



 α σ φββ

a m 2

21

2 = 2 12 ( ) + [ ] ( ) + + α σφ σβ φ

a

M

g

g

S

k Ak g M

gH 3

21

21

13

2 2 1 21

3 = − ( ) − − + ρ ( ) − σ ψ β

T

/

a S

M

g

M

g

Mk g M

gH g Mk k 4

21

2

2

3

13

2

13 13

2 2 1

3

= + 05 1 2 



 



 + ( ) − [ ] ( ) + + +

















ρ β β

σ

ψ β σ σσ

T

d d

/ ./ /

a k M

g

k g

g

5 21 P

2

13

2 2

1 2

= + ( ) + ( ) − +





 



β  β

σ

σ

σ

d

P k g

H

A S Mg

Mk g

g = ( ) − − + ( ) + + 



 



 −+ + ( ) ( ) − − 



 

















ψ

σ

β σφ σφ σ

σ

13 2 1 13

3

2 1

3

05 1

1 2

1 / / ./ T

d

m = + βσβ φ [ ] ( ) 2 2. +

The problem has now been reduced to a simple numerical quadrature of eq C23 using the auxiliary relations of eq C24–C26. A FORTRAN program to carry out the integration is listed in Appendix E as

PFTSYN.FOR.

This approximation is inferior to method 1 since the variables in eq C23 are not strictly separable.

Nevertheless, predictions for modest times compare quite well to those of method 1.
35

Ts Tf

X (t)

δ (t)

x = 0

x = X (t),

Bottom of

Permafrost at t

Thawed, 2

Thawed, 3

x = X + δ,

Depth for

qu = qg

Frozen, 1

qf (t)

qu (t)

qg = ku G

Figure D1. Energy flows at base of permafrost.

APPENDIX D: ENERGY FLOWS AT THE PERMAFROST BASE

The heat flows at the base of the permafrost layer determine the rate of movement of the permafrost bottom. Examine Figure D1 to clarify the concepts involved. At any instant of time, an amount of energy,

qf(t), is conducted away from the phase-change interface through the frozen layer and to the interface from

the thawed zone 2, qu. It is important to realize that the constant geothermal heat flow qg = kuG will only

equal qu(t) at equilibrium. During movement of the permafrost base, qu(t) can be greater or less than qg. We

can examine the transient behavior of these terms. Let

Q q t

k G

k

k G

T Xt

x

f

f

u

f

u

= = ( ) 1 ( ) , 1 ∂

∂ (D1)

Q q t

kG G

T Xt

x u u

u

= = ( ) ( ,) . 1 2 ∂

∂ (D2)

In terms of dimensionless parameters

Q

k

f

21

=

+

[ ] ( ) + +

+





















1

2 2

1 21 σ α σφ

βσβ φ

( )

(D3)

Qu = σβ φ ( ) + +

σβ

2 2

(D4)

where we assume that ku = k2, i.e., the thermal conductivity of the thawed zone is constant. The results of

example 2 for Prudhoe Bay are tabulated in Table D1. Clearly, the heat flow from the deep layer greatly exceeds the geothermal heat flow for much of the permafrost formation period. Also, note how rapidly the

heat flow out of the frozen zone Qf drops to slightly more than the flow from the thawed layer Qu.

The temperature in the frozen zone adjusts very quickly toward the equilibrium result of a linear temperature distribution; however, the thermal zone approaches equilibrium very slowly.

Surface temperature increase

Suppose the permafrost has thickness Xo,

is growing and the surface temperature initially increases by a certain amount and is

held constant for several thousand years.

What is the effect on the permafrost? The

bottom of the permafrost will continue to

grow for several hundred years before starting to thaw, but Table D1 assures us that

this will be negligible. At t = 225,000 years,

Xo = 626.5 m and the permafrost growth

will only be 27.5 cm during the next 1000

years. Thus, we can assume that the permafrost thickness remains essentially constant

while the temperature in the frozen zone

Table D1. Freeze at Prudhoe Bay, Alaska.

∆T1 = 12.69°C φ = 0

k21 = 0.5812 ST = 0.1827

α21 = 0.3915 Xe = 763.5 m

Permafrost

Time depth

(yr) (m) β Qu Qf

1 4.44 1.9060 2.423 159.057 0.9251

350 79.9 1.6717 2.196 8.986 0.9405

3500 219.3 2.3380 1.855 3.352 0.9629

35,000 461.4 4.8154 1.4153 1.6351 0.9826

100,000 567.8 7.7133 1.2593 1.3375 —

225,000 626.5 11.1633 1.1792 1.2153 —

775,000 687.7 20.3168 1.0984 1.1090 —

q

q

f

fe

=

+ +

1

2

1 21 α

β β( )
adjusts to its new equilibrium value. Figure D2 sketches the problem. Obviously, after infinite time the new temperature profile is

as shown. However, the temperature will adjust to near the new

equilibrium in a relatively short time. This is a linear problem in

non-phase-change conduction and has been solved by Lachenbruch et al. (1982). The transient temperature is

TT T T x

X

T T

n

e n

x

X

o

n

M = ′ + − ( )′ + − ( )′ 



 



 =

∞

s fs ∑ o s s –

o

2 1

1 π

sin π (D5)

where

M n t

t

t

X = =

2 2 2

4

π

c α c of 4 , .

tc is a characteristic time for sensible temperature changes. The equilibrium temperatures are simply

T T TT x

X

o o o s ss o

=+− ( ) (D6)

T T TT x

X ∞ = ′ + − ( )′ s fs o

. (D7)

The change in sensible heat, going from the state at t = 0 to the state at t = ∞, is

Q CT To X

s fs s o ∞ = ( ) ′ − 2 . (D8)

The change in sensible heat at any time t is

Q C T T dx Q C T TX e

n

X

o

n

M n sf o s

f

2 s so

o

o

= − ( ) =+ − ( )′ [ ] − −( ) ∫ ∞

=

∞

∑ 2 1 1

1

2 π

–

. (D9)

Thus, the relative change in sensible heat is

Q Q

Q

e

n n

M s s n

s

∞

∞ =

∞ − = −− ∑ [ ] ( ) 4 1 1 2

1

2 π

–

. (D10)

Note that this quantity does not depend upon the surface temperatures. The relative change is shown in

Table D2. The sensible heat change attains 93% of its ultimate value at t/λ = 1.0 (t = 1666 years) and 99% at

t/λ = 1.78 (t = 2966 years). The sensible heat changes would be essentially completed after about 1670

years. From this time on, the bottom of the permafrost would slowly melt.

Change in frozen zone temperature gradient at bottom of permafrost

The time required for the temperature gradient, in the frozen zone, to

change is important since this quantity will determine the rate of change of the

permafrost bottom depth. The gradient at the depth Xo can be found from eq

D5 and is

δ ∂

∂ = = − ′

+ − ( )′ ( ) − = =

∞

∑ T

x

T T

X X

TT e

x X

o

n

n M

o

f s

o o

s s

2 1

1

– . (D11)

36

Table D2. Relative change

in sensible heat.

t/λ

0.25 0.438

0.50 0.236

1.0 0.069

1.50 0.020

1.78 0.010

2.0 0.006

Q Q

Q

s s

s

∞ −

∞

t = 0 X0

Tf Ts

 ' Ts

 0

t = ∞

Figure D2. Permafrost equilibrium temperature profiles.
37

The initial temperature gradient, at x = Xo, is

δ ∂

∂

o o T

x

X T T

X = ( ) = − o f s

o

, . 0 (D12)

Note that this simplified form of δo gives virtually the same

heat flux as the value from eq D1, when x = Xo. Now

δ δ

δ

o

o

o

o

n M

n

T T

T T

e − = ′−

− + −( ) 



 



 −

=

∞

∑ s s

f s

12 1

1

. (D13)

The time for a given change in the gradient can be closely approximated as

t

T T

T T

o o

o

λ π δ δ

δ

=

−  −









 −

′−









































4 2

1

2 ln .

f s

s s

(D14)

Table D3 shows the results for Prudhoe Bay.

It would take about 490 years for the bottom growth to cease and 900 years for the bottom gradient to

change significantly. Thus, the approximations used in the derivation of eq D5 are acceptable.

Bottom melt

The temperature in the frozen zone requires 1666 years for sensible heat adjustment, leaving ∆t = 15000 –

1666 = 13,334 years for bottom melt (the interglacial is 15,000 years long). The energy balance at the bottom

of the permafrost is

L dX

dt

k T T

X = Aq − ′ f − f s

g. (D15)

where A is the fraction of the geothermal energy that goes into melting; it can exceed 1.0.

The heat flow from the thawed zone is greater than the geothermal heat flow, as we

have seen. For the example discussed here, A = 1.179 at the beginning of thaw and will

decline towards 1 as thaw proceeds. The permafrost thickness Xf after ∆t years is given by

X b X b

X b f a X f

o

− o

+

+





 

 ln  = + (D16)

where

a A q t

L

b

kT T

Aq =− =− q kG ( ) − ′ = g ff s

g

g u

∆ , ,.

The final permafrost thickness is strongly dependent upon the value of A. Table D4 shows values for the

Prudhoe Bay example.

The heat flow from the thawed zone varies continuously during the thaw, denoted in eq D15 as Α(t) qg.

The heat flow at equilibrium is such that A = Ae = 1.0. Thus, let A(t) be a linear function of X, given by

AAA y y

y = − ( ) A −

−





 

 o e  + e

e

e 1 (D16a)

Table D3. Time for change in gradient.

T o

s = –13.69°C; T ′

s = –11.0; Tf = –1.0

t δ δ

δ

o −

t/λ (years)

δ δ

δ

o  −



 





*

0.01 0.301 502 0.0314

0.0297 0.3421 570 0.0440

0.1 0.540 899 0.1022

0.20 1.446 2408 0.20

0.212 ∞ ∞ —

* Calculated without eq D14 approximation.

Table D4. Permafrost thickness.*

Xf

A (m)

1.179 591.0

1.0895 605.5

1.0 620.0

*Final thickness

after surface temperature increase,

Prudhoe Bay.
38

X T T

k G e f s

uf

= − ′

Table D5. Permafrost

thickness after thaw.

Ao =1.1792.

Ae yf Xf (m)

1.0 0.9678 606.4

1.170 0.945 592.0

1.179 0.9434 591.0

where y = X/Xo

Ao = value at thaw commencement

Ae = equilibrium value

= new equilibrium permafrost thickness.

Equation D15 then has the following solution

ln ln ay a y

a a

a

a

y a

y a

a

a a 1 f

2

2 f

1 2

2

3

f 4

f 5

5

4

1 f

+ −

+ −











 − +

+





 



 +

+





 





















+ =

1

1

1

1

2 0 τ (D17)

where

a

q X

kT T

A A

y

1

g o

ff s

o e

e 1 = ( ) − ′

( ) −

( ) −

a Ay

A A

y

2 ee o e

1 e

= − ( ) −

( ) −

a aa 3 1 = + 2

2 4

a

a a

a

4

2 3

2 1

= −

a

a a

a

5

2 3

2 1

= +

τf

ff s

o

2 = kT T t ( ) − ′

LX

∆

and

y

X

X f

f

o

=

is the permafrost thickness after Dt years (Table D5).

Note that these results agree quite well with the values with a constant thawed zone heat flow, i.e., constant value of A(t). For this case, 15,000 years is nearly enough time to thaw back to the new equilibrium

thickness of 601.5 m.
39

APPENDIX E: FORTRAN PROGRAMS FOR NUMERICAL QUADRATURE

OF ENERGY EQUATION

Program PERM
40
41

Program PFTSYN
42
43
44
1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

6. AUTHORS

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION

 REPORT NUMBER

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITORING

 AGENCY REPORT NUMBER

11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

14. SUBJECT TERMS 15. NUMBER OF PAGES

16. PRICE CODE

17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION OF ABSTRACT

 OF REPORT OF THIS PAGE OF ABSTRACT

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and

maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information,

including suggestion for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington,

VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

April 1995

Permafrost Formation Time PR: 4A161102AT24

TA: BS/0045

Virgil J. Lunardini

U.S. Army Cold Regions Research and Engineering Laboratory

72 Lyme Road CRREL Report 95-8

Hanover, New Hampshire 03755-1290

Office of the Chief of Engineers

Washington, DC 20314-1000

Funding also provided by U.S. Army Strategic Environmental Research and Development Program 2.6.00 78 203

Approved for public release; distribution is unlimited.

Available from NTIS, Springfield, Virginia 22161.

48

UNCLASSIFIED UNCLASSIFIED UNCLASSIFIED

The age of permafrost is closely linked to the time required for soil systems to freeze, since the permafrost must be at least

as old as the formation time. Cycles of freeze–thaw will complicate the relation between the freeze rate and the age. A

model based on pure conduction heat transfer with freeze–thaw is used to predict the time required for a given thickness of

permafrost to develop, either heterogenetically or syngenetically. The formation time is a function of the long-term geothermal gradient (initial temperature of the thawed soil), the ratios of the frozen to thawed thermal properties, and the temperature history of the upper surface of the permafrost (higher than the air temperature). The simple theory allows universal

graphs to be produced that predict the formation time for a given thickness of permafrost. Realistic soil property ratios and

paleotemperature scenarios will then lead to estimates of the formation time of permafrost for a specific site. The model indicates that deep permafrost (more than 1500 m) requires formation times on the order of the complete Quaternary Period.

Age of permafrost Permafrost

Paleotemperatures Permafrost growth and melt