Three-Dimensional 

Heat Conduction in 

Permafrost Beneath 

Heated Buildings

GEOLOGICAL SURVEY BULLETIN 1052-B

Three-Dimensional 

Heat Conduction in 

Permafrost Beneath 

Heated Buildings

By ARTHUR H. LACHENBRUCH

EXPERIMENTAL AND THEORETICAL GEOPHYSICS

GEOLOGICAL SURVEY BULLETIN 1052-B

A discussion of three-dimensional heat 

conduction in permafrost with special 

reference to thermal disturbance caused 

by heated structures

UNITED STATES GOVERNMENT PRINTING OFFICE, WASHINGTON : 1957
UNITED STATES DEPARTMENT OF THE INTERIOR 

FRED A. SEATON, Secretary

GEOLOGICAL SURVEY 

Thomas B. Nolan, Director

For sale by tbe Superintendent of Documents, U. S. Government Printing Office , Washington 25, D. C. - Price 75 cents (paper cover)
CONTENTS

Page 

Symbols_______________________________________________ Facing 51

Abstract____________________________ __ 51

Introduction. ________________________________ 51

Acknowledgments. ___._________________________________ 52

Outline of the heated-building problem_________________________ 52

The principal disturbance, Oi_____________________________ 56

The seasonal disturbance, QZ- _______________________:__ 62 

Effect of idealizing the problem__________________________ 66

Summary_________________________________________ 68 

Literature cited___________________________________ 69

ILLUSTRATIONS

[Plates in pocket] 

PLATE 1. Solid-angle grid.

2. Solid-angle grid showing correct positioning for calculation of 

the thermal effect of a 40- by 100-foot region at two positions.

3. Weighting factor, V, for calculation of seasonal disturbance.

Page 

FIGURE 22. Weighting factor, <i>, as a function of r2/at for calculation of

principal disturbance..__________________ ______ 59

23. Principal disturbance 20 feet beneath a 40- by 100-foot region

with temperature B° above surrounding surface._______ 62

24. Calculated temperature deviation from the undisturbed mean

for large values of time at depth of 20 feet__________ 67

TABLES

Page 

TABLE 1. Calculation of principal disturbance._____________________ 61

2. Summary of principal-disturbance values______________ 61

m
SYMBOLS

6 or T temperature (degrees C).

t time since initiation of temperature disturbance in S (sec)

t0 time of initiation of temperature disturbance in S (sec)

t' =t+t0

a thermal diffusivity (cm2 per sec)

w angular frequency (w= 1.99X 10~7 radians per sec for a period

equal to 1 year)

Q solid angle (steradians) 

S finite region of plane z=Q 

p, \ polar coordinates in plane z=Q 

r = [(x-x'Y+(y-y'Y+Z^=[p^z^ ^ 

A amplitude of temperature variation in undisturbed portion

of plane z=Q (degrees C) 

B mean value of surface temperature in disturbed region, S, of

plane z Q (degrees C) 

C amplitude of temperature variation in disturbed region, S, of

plane z Q (degrees C) 

D =C A, amplitude of seasonal temperature disturbance in

disturbed region of plane z=0 (degrees C)

erf x =-f=- I e~u* du

V"" Jo

2 C*°

erfc x =-7= I e~"2 rfw

weighting factor for principal-disturbance calculation

2 r ~i^ , r

weighting factor for seasonal-disturbance calculation

sm

IV
EXPERIMENTAL AND THEORETICAL GEOPHYSICS

THREE-DIMENSIONAL HEAT CONDUCTION IN 

PERMAFROST BENEATH HEATED BUILDINGS

By ARTHUR H. LACHENBRUCH

ABSTRACT

The general Green's function solution has been integrated for the case of heat 

conduction in a homogeneous semi-infinite medium in which the temperature at 

the surface varies sinusoidally with time but the mean temperature and amplitude of the variation are different within and outside an arbitrarily shaped region 

at the surface. The amplitude and mean temperature can be treated as functions 

of position within the arbitrary surface region. For certain simple surface regions 

the results can be expressed in terms of tabulated functions. Numerical results 

for the general case can be obtained by simple graphical procedures.

The results can be applied to the study of disturbances in ground temperature 

caused by the presence of bodies of water or by engineering surface modifications 

such as those produced by erecting a heated building. The primary application 

of such studies is in high-latitude regions where much of the undisturbed ground 

is perennially frozen. In such areas, a method of predicting the extent of thawing 

induced by various modifications of the temperature of the ground surface is 

important in problems of engineering design and logistics.

INTRODUCTION

The temperature distribution in permanently frozen ground is of 

considerable importance in Arctic and Antarctic regions. Inasmuch 

as engineering operations in these regions are carried out on normally 

permanently frozen mediums (snow, ice, and permafrost) a knowledge 

of the temperature distribution in these mediums under various disturbing influences is important to the engineer so that he can, for 

example, design in such a way as to keep critical parts of structural 

foundations frozen, or water-supply sources and conduits unfrozen. 

Several problems of more scientific than practical nature, such as the 

relation between the age of a lake and the transient disturbance of 

ground temperature it produces, can also be approached through the 

study of idealized models such as those considered in this paper.

One of the important general problems is to determine the disturbance of subsurface temperature that results when the temperature

51
52 EXPERIMENTAL AND THEORETICAL GEOPHYSICS

at the ground surface within a finite region, S, is different from the 

surface temperature characteristic of the area outside S. The conditions correspond to the presence of natural features such as a lake or 

river, or a region in which the thermal properties of the surface cover 

are appreciably different from those characteristic of the area in general, and to modifications of the surface such as those resulting from 

erecting a building, stripping the vegetation, or emplacing a gravel 

fill.

Inasmuch as probably the most important application is to study 

of the temperature distribution beneath a heated building, this general 

problem will be referred to for convenience as the heated-building 

problem,

ACKNOWLEDGMENTS

The work reported here was done as part of an extended program 

of investigation of the thermal properties of permafrost by the U. S. 

Geological Survey. The work has been supported jointly by the Survey and the Office of Naval Research and the Bureau of Yards and 

Docks of the U. S. Navy, and during part of the time by the Snow 

Ice and Permafrost Research Establishment of the Corps of Engineers, 

U. S. Army.

Francis Birch and D. C. Hilton offered valuable comments on the 

manuscript from the physical and engineering points of view respectively. S. H. Lachenbruch provided a thorough check of the paper's 

mathematical content.

OUTLINE OF THE HEATED-BUILDING PROBLEM

For study of the general aspects of ground temperatures the simple 

two-parameter sinusoidal approximation for the undisturbed surface 

temperature, A sin at, usually serves quite well. The temperature 

origin is taken as the mean annual temperature just beneath the 

ground surface in the undisturbed area. The angular frequency, w, 

corresponds to a period of 1 year. In the finite region, S, of the surface in which the surface temperature is disturbed, the significant 

parts of the problem can usually be treated by representing the temperature by B(x,y)-\-C(x,y) sin ut. Here B(x,y) is the difference in 

mean annual temperature between the undisturbed area and the surface point (x,y) in the disturbed region S. C(x,y) is the amplitude of 

temperature variation at (x,y) in S.

It is generally convenient to assume that the surface temperature 

presently characteristic of the undisturbed region (outside S) formerly 

obtained over the entire surface, and that any initial transient effects 

had died out by the time, #0 , that the disturbing condition was introduced in the region, S.
THREE-DIMENSIONAL HEAT CONDUCTION IN PERMAFROST 53

Thus the heated-building problem can be stated formally as follows: 

Find the temperature, Q(x,y,z,t'), in the semi-infinite solid z^>Q, subject 

to the following conditions:

d2 e . d2 e . d2 e i de . .. m

(1) x Cty2 d02 a dt

'=0 (2)

all

, ^x., , x ,. 0 ,0v b) £o<£', (x,y) not ui £ (3)

Q(x,y,0,t')=B (x,y) + C(x,y) sin «*', 2=0, t*<V , (x,y) in 5. (4)

Here tQ is taken sufficiently large that transient initial effects have 

died out by the time (t'=t0) of initiation of the disturbance in S. 

f (z) represents the normal geothermal profile assumed independent 

of time.

It has been implicitly assumed above that the ground surface has 

no appreciable relief, and thus can be represented by the plane 2=0. 

In addition, we shall assume that the ground is homogeneous and 

neglect the effects of latent heat resulting from the change of state 

of interstitial water or ice. The physical implications of these limiting 

assumptions will be discussed briefly later.

With latent heat neglected the relations become linear and it is 

generally possible to consider B and C to be constant over finite subregions of S and solve the problem for each of these subregions. 

These solutions are later superposed to give the collective effect of 

the entire region, S. Thus the necessity of considering B and C as 

functions of position is eliminated. A further simplification resulting 

from neglecting latent heat arises from the fact that it is then possible 

to disregard the normal static geothermal profile, f(z), and consider 

the surface effects separately. To get the complete solution, f(z) 

can be added to the final result. It can be neglected, however, for 

most problems involving small heated structures.

With the assumption of homogeneity and no latent heat, the 

problem is reduced to

d2 e . d2 e . 52 e i de , _ ,. /CN (5) v ' dx2 ch/2 dz2 aftt 

Q(x,y,Q,t)=A sin ut, z=Q, £>0, (x,y) not in S (6) 

Q(x,y,Q,t)=B+C sin at, 2=0, *>0, fotf) in S (7)
54 EXPERIMENTAL AND THEORETICAL GEOPHYSICS

and

0(z,2/,2,0)=00 (2,o), 3>0, *=0 (8)

where t=t' to.

This solution can be expressed as the sum of three independent

solutions as follows:

(9)

where 60 is the limit, as t' becomes infinite, of the solution to the 

problem

WT 1 ftT ^£=1 ^£, s>o, r>o (10) d22 adt '

T(Q,t)=A sin orf', 2=0, t' >0 (11) 

T(z, 0)=0, 2>0, «'=0. (12) 

GI is the solution to

5*0! d2 0t d2©! 1 50! , . .

^ g + ^ 2 + ^. 9 ==~ ~5TT' 2 x>0, fP>U (16) do;2 dy2 ds2 a d£

0i(a;, y, 0, 0=5, <>0, s=0, (x, y) in S (14) 

QI (*, y, 0, 0=0, C>0, «=0, (x, y) not hi ^ (15) 

0X (x, y, 2,0) =0,^=0, 2>0. (16) 

and 02 is the solution to

d2 02 l d0

©2 (*, y, 0, f)=D sin coi, s=0, C>0, (^, y) in ^ (18)

02 (x, y, 0, 0=0, s=0, i>0, (x, y) not hi S (19)

e2 («,2/, 2, 0)=0, 2>0, i=0 (20)

where D= C- A. (21)

The solution for 00 is well-known (see, for example, Carslaw and 

Jaeger, 1947, p. 47, eq 3) and is given by

00 (2, t)=Ae~* sin ut-z' (22)
THREE-DIMENSIONAL HEAT CONDUCTION IN PERMAFROST 55

The solutions for Qt and 62 can be approached through the use of 

Green's functions in heat conduction. The following method has 

been used by Birch (1950) for a similar problem.

The temperature T (x, y, z, t) in a semi-infinite medium initially at 

zero, due to a surface temperature <j> (x, y, f) in S and zero outside S, 

is given by Carslaw and Jaeger (1947, p. 308, eq 3) as

**"'

By making the substitution

and reversing the order of integration, equation 23 is reduced to

(24)

where r2=(x-xr') 2 -\-(y-y') 2 -\-zi (that is, r is the distance between the 

field point (x, y, z) and the element.of surface dx'dy') . 

Equation 24 can be rewritten

,y,z)=^ JJ / J^*^V- (25)

^- 2-^/of

where dQ is the element of solid angle subtended by dx'dy' at the 

field point (x, y, z). By setting #(Vy,£ ^32)=^ and integrating 

with respect to /3, an expression for 61 is obtained:

e1 (*)^,e)=!JJfe-S+erfc^}rfa. (26)

By setting

B',J/',£ j =3 ) D sin cof t - 

an expression for 62 is obtained:

428961 57 2
56 EXPERIMENTAL AND THEORETICAL GEOPHYSICS

'. The sum of these two temperature fields gives the subsurface 

temperature disturbance due to a heated building (or other surface 

disturbance) in the region S. 6t (equation 26) will be referred to as 

the principal disturbance inasmuch as it gives the temperature disturbance due to a building (or other feature) except for the effects 

of seasonal fluctuations. It can be thought of as giving a mean 

annual value of temperature disturbance at any subsurface position 

at any time after the surface region, S, was disturbed (for example, 

after the building was erected) . The sum of 60 and 92 (equations 22 

and 27) will be referred to as the seasonal effect as it gives the seasonal oscillation about the principal disturbance. 92 will be called 

the seasonal disturbance. These three equations (22, 26, and 27) 

provide a complete solution to the heated-building problem for a 

homogeneous medium in which latent heat is neglected. The next 

two sections are discussions of equations 26 and 27 respectively.

THE PRINCIPAL DISTURBANCE, 0!

The principal disturbance as defined above is the temperature 

resulting in the ground (z^>0) from a constant temperature difference, 

B, between a point inside the disturbed region, S, and the mean 

annual temperature of the undisturbed region outside S.

In equation 26 the disturbed region may have an arbitrary configuration (for example, the shape of the foundation of a heated 

building). It is evident, however, that this equation cannot be 

integrated directly for arbitrary S. The integration can be performed 

for certain simple situations, for example under the vertex of a circular 

sector of radius R and central angle X in the plane z=Q. In this 

case

(28)

(29)

For a circular region this reduces to

^0, 0, z, t)=B ferfc -JL- 1= erfc ^±^1. (30) 

L 2-y/^ T/z2+R* 2^ J

The first term in equation 30 corresponds to the effect that would 

result if S covered the entire plane z=0, and the second term may be 

thought of as a correction for the finite radius, R, of S. By setting 

X=TT and 7T/2 in equation 29, results are obtained for a semicircle and 

a quarter circle respectively. These results may be used to estimate
THREE-DIMENSIONAL HEAT CONDUCTION IN PERMAFROST 57

the principal disturbance under the edges and corners of a building 

for depths that are small with respect to building dimensions. Under 

these conditions the correction term is negligible and the disturbance 

is respectively one-half and one-quarter of that produced by the 

"infinite building" (that is, S= en tire plane z=Q). 

From equation 26 we obtain the relation

T>

limit Gifo y, z, *)= Q(x, y, z ; S) (31)

£-»oo 61T

where fi is the solid angle subtended by S at the field point (x, y, 2). 

This special case is in agreement with results from potential theory 1 

A useful result for computing the steady-state principal disturbance 

at a point (x, y, z) beneath or to one side of the rectangular region 

b<y<-}-b is given by

61(3, y, 3)=s- -< tan -

tan

This result is obtained by substituting the appropriate expression 

for J2 in equation 31.

From equation 29 we obtain

erf c . _ 2-yW

(33)

which represents the principal disturbance at a depth z beneath the 

vertex of the sector of a circular annulus with central angle X, inner 

radius Ri} and outer radius R2 . The principal disturbance for a region 

of arbitrary shape may thus be evaluated by dividing the region into 

concentric annular sectors and applying equation 33. This method 

can be used to give any desired degree of accuracy by choosing X and 

R2 Ri sufficiently small near the region's periphery. It can also be 

used to give rapid approximations in many cases. Generally, however, 

the calculation is simpler if the integration with respect to equation 

26 is performed graphically. The graphical methods described below 

have been used previously by Birch (1950).

To carry out the graphical integration, a grid is constructed by 

gnomonic projection of meridians and parallels of a hemisphere onto
58 EXPERIMENTAL AND THEORETICAL GEOPHYSICS

its north polar tangent plane. The meridians and parallels are so 

selected that the plane of projection is subdivided into concentric 

annular sectors each subtending equal elements of solid angle at the 

sphere's center. The field point (x, y, 2) is identified with the point 

at the center of the sphere, the projection plane is identified with 

the ground surface. Thus the radius of the sphere corresponds to 

the depth, 2, of the field point beneath the surface.

A grid of this type is shown in plate 1. It is divided into twelve 

concentric zones designated by letters A through L. The region in 

zone A subtends 2 percent of the total solid angle (2ir steradians) 

subtended by the plane 2=0. Zone B subtends 4 percent; zone C, 

6 percent; zone D, 8 percent; and zones E through J each subtend 

10 percent of that subtended by the entire plane. Zones K and L 

each subtend 5 percent. Thus, according to equation 31, the steadystate effect of a circular building whose periphery coincided with the 

outer limit of zone L would be 90 percent of the effect of a building of 

infinite dimensions covering the entire plane 2=0. The zones are 

cut by radial lines in such a way that each element of area of the grid 

subtends 0.1 percent of the solid angle subtended by the entire plane 

2=0. The depth, z, of the field point corresponds to 1 iuph. on the 

grid (pi. 1). Columns 2 and 3 of the table in plate 1 give the inner 

and outer radius of each zone in units of 2. This grid can be used to 

calculate the effects of surface disturbances at horizontal distances 

from the field point as much as about 10 tunes the depth. Grids 

with other specifications are easily constructed.

The integrand in equation 26 depends only upon the ratio 

{where a is the thermal diffusivity of the medium). In figure 22, the 

value, $, of the integrand in equation 26 is plotted against r^/at. For 

Any assumed at, the value of the integrand corresponding to any r is 

read directly from this graph. For purposes of the graphical integration it will be assumed that all points in a given zone are at the same 

distance, r, from the field point and that this distance is given by

r=VsMT2J that is, P=s2 [l+(p/*) 2]

where ~p is the mean radius of that zone. The error introduced by 

this simplification should be taken into account in designing grids 

for special purposes.

The method of integration is as follows. First a depth, z, is selected 

at which temperature information is desired. Next a plan drawing of 

the outline of the building foundation is made on tracing paper. The 

drawing is in units of z-, that is, if the depth under study is 10 feet, 

then, as 2 corresponds to 1 inch on the grid (pi. 1), the building will 

be drawn on the scale of 1 foot=0.1 inch. Next, the tracing paper
THREE-DIMENSIONAL HEAT CONDUCTION IN PERMAFROST 59

 i s

£

I

 It

II <e>
60 EXPERIMENTAL AND THEORETICAL GEOPHYSICS

is positioned on the grid in such a way that the horizontal coordinates 

of the field point coincide with the center of the grid. If we want to 

calculate the principal disturbance under a corner of this building, the 

drawing is placed with this corner over the center of the grid. Next 

a value of od, is selected that depends upon the time for which temperature information is desired and the assumed diffusivity of the 

medium. The number of solid-angle elements subtended by the 

building in zone A is then counted. This number multiplied by 

0.001 gives the contribution that this zone will make after equilibrium 

is established. Next multiply z2/at by the value of l-j-(p/s) 2 for zone 

A given in column 4 of the table in plate 1. This product gives r*/od 

for zone A. The corresponding value of the ordinate in figure 22 

gives the amount by which this zone contribution must be weighted. 

This weighting factor multiplied by the steady-state contribution of 

the zone gives the contribution of zone A to the principal disturbance 

at (x, y, z) for the at assumed. The process is repeated for all zones 

intersected by the building's outline. The contributions of all zones 

are then added to get the principal disturbance.

The accuracy of the graphical integration with respect to £2 can be 

checked quickly by comparison of the effect under the center of a 

circular region calculated by equation 30 with that obtained by the 

use of the grid. This was done for several cases and for a variety of 

assumed values for R, z, and at. The effects calculated by the two 

methods agreed within 0.005 B. The error introduced by counting 

solid-angle increments on the grid can be checked by comparing the 

total with the computed solid angle subtended by the region at the 

field point if the geometry: of the region permits (for example, for a 

rectangular region equation 32 may be used).

As an example of the application of the graphical method, we shall 

calculate the principal disturbance at 2 points 20 feet beneath a surface assumed to be disturbed oyer a [rectangular region 40 feet by 

100 feet. To fix ideas we; may jthink of the effect of a building 40 

feet by 100 feet in wMch the m0an annual floor temperature is constant and exceeds the mean annual teniperatiire of the ground surface 

outside the building by 5°C. Let position 1 correspond to a point 

20 feet beneath the center of such a building, and position 2 to a 

point 20 feet beneath a point on the transverse center line 10 feet 

from the building. The appropriate scale and positioning of the 

region's outline for these two calculations are shown in plate 2.

Table 1 illustrates the steps in the calculation of the principal disturbance for an assumed value «Z=5X105 cm2 . Table 2 gives a 

summary of the results of table 1 and a comparison of the equilibrium principal disturbance as computed with the grid and with the 

solid-angle formula (eq 32). The discrepancy corresponds to a net
THREE-DIMENSIONAL HEAT CONDUCTION IN PERMAFROST 61

error in counting of about one solid-angle unit on the grid (0.002r 

steradians).

The third column in table 1 shows the degree of advance toward 

equilibrium for each zone for the at assumed. The portion of the disturbed region in zone J contributes on the order of only 10 percent 

of its full equilibrium effect while that in zones A through D contributes more than 90 percent. Because in position 2 the region subtends 

more of the outer zones, it is approximately 15 percent farther from 

equilibrium for the assumed at (table 2). To contour the temperature 

at the 20-foot depth for a given at, simply shift the tracing over the 

grid, changing only the values corresponding to columns 4 and 5 in 

table 1. Although the first calculation might take about 20 minutes, 

subsequent ones may each be made in about 5 minutes. To study 

the temperatures at the 10-foot depth, the scale of the region's outline is doubled and the process repeated. To study the change of 

temperature with time at a given point, values corresponding to those 

in columns 2 and 3 are changed but those corresponding to column 4 

are not, and the calculation is again accomplished speedily. The 

change in principal disturbance with time, as calculated with the 

grid method, is shown graphically for positions 1 and 2 in figure 23.

TABLE 1. Calculation of principal disturbance 

[See pi. 2; fig. 22]

Zone

A....... ____ ..

B... .............

O... ..____ .......

D...... .......

E....... .........

P... ..........

G................

H...... .......

I..... . .....

J..... ...

Principal disturbance, 

units 10-3B...

<rf=5X105 cm2

S-IK!)']

0.752 

.804 

.896 

1.053 

1.326 

1.773 

2.491 

3.752 

6.309 

12. 852

Weighting 

factor, *

0.945 

.939 

.930 

.914 

.880 

.830 

.740 

.600 

.370 

.093

Position 1 

2=20ft,z=0, y=0

Number of 

Q units 

subtended

20 

40 

60 

80 

99 

67 

46 

34 

9 

0

455 

(equilibrium)

Zone contribution 

(units of 

KH.B)

18.9 

37.6 

55.8 

73.1 

S7.1 

55.6 

34.0 

'20.4 

3.3 

.0

385.8 

(transient)

Position 2 

2=20 ft, z=30ft, y=0

Number of 

£J units 

subtended

0 

0 1*A 

MH 

31 

36 

39J4 

42 

27 

2

193M 

(equilibrium)

Zone contribution 

(units of 

WHS)

0 

0 

1.6 

13.3 

27.3 

29.9 

29.2 

25.2 

10.0 

.2

136.7 

(transient)

TABLE 2. Summary of principal-disturbance values

Position

l_._.-_. _______________

2

at*

Table 1

0. 4555

.1945

..

Equation 32

0. 4559B

. 19455

at=5X105 cm»

0. 3865

.1375

Percent of 

equilibrium at at-5X10' cms

84.8

70.6
62 EXPERIMENTAL AND THEORETICAL GEOPHYSICS

3 YEARS 

(«=0.01)

FIGTTKE 23. Principal disturbance 20 feet beneath a 40- by 100-foot region with temperature B° 0 

above surrounding surface. Abscissa is atom*.

In the example given, the mean annual temperature disturbance, 

B, was assumed to be independent of position throughout the disturbed surface region. If it were different over some subregion (for 

example, in the vicinity of a furnace or an unheated wing), the number of solid angle units subtended by such a subregion would be 

weighted proportionally and the calculation carried out as before.

THE SEASONAL DISTURBANCE, 62

The seasonal disturbance, 62 , as previously defined, is the temperature resulting in the ground (2>0) from a temperature in S of 

D sin ut while the temperature of the surface outside S is maintained 

at zero.

Equation 27 can be rewritten

e,G», y, 2, «)=~if/{ Jf /"in

The first integral in equation 34 represents the steady seasonal disturbance that will obtain after the surface disturbance in S has persisted for some time (£>rYa). The second integral describes the 

manner hi which the seasonal disturbance builds up for small values 

of time. The expression in braces hi this integral will have to be 

evaluated numerically before the grid of the previous section can be 

used to perform the integration with respect to fl.
THREE-DIMENSIONAL HEAT CONDUCTION IN PERMAFROST 63

For most problems the first integral in equation 34 is the more 

important quantity, as the second may eventually be neglected. The 

expression in braces in the first integral can be evaluated as follows:

f* ^ dyr2 C m 6>r2

=sin«* cos -P-& e-*pdfi- cos <at\ smf^e~^d^. (35) J o lap* J o 4af}*

By integrating each of the expressions on the right side of equation 35 

by parts and using the relations

2g+sin 2g] 

(Bierens de Haan, 1939, table 369, formula 1)

fa e-z2 cos^ ^=^e-2'[cos 2<z-sin 2q] 

Jo x2 x2 4g

(idem, table 369, formula 2)

fV^sm^<fo=^<r*sin 22 Jo « 2

(idem, table 263, formula 12)

rV*2 

Jo x2 2

(idem, table 263, formula 13) 

we obtain

(36)

By substituting this result in equation 34,

dti transient term. (37)
64 EXPERIMENTAL AND THEORETICAL GEOPHYSICS

As was true for the integral in equation 26, the integration in 

equation 37 is easily performed for a point (x,y,z) under the vertex of 

a circular sector, S, of radius R and central angle X.

r^ / 

2«sinf ut lz2+R2

transient term. (38)

For a circular region this reduces to

transient term. (39)

In direct analogy to equation 30 the first term in equation 39 corresponds to the effect that would result if S covered the entire plane 

s=0, and the second term represents edge effects arising as a result of 

the finite radius, R, of S.

As explained above, the seasonal effect is obtained by adding the 

seasonal disturbance, 62 , and the normal undisturbed seasonal variation, 60 (eq 22). Thus, after the transient effects have disappeared, 

the seasonal effect under the vertex of a circular sector of central 

angle X and radius R is

or, using equation 21,

e0+92= [^1-A (4-C)] e^ sin (crf-g VS)+

*V-0)^.-^&*(^J?+*<J£). (41)

As in equation 30, equation 41 can be used to estimate effects under 

straight edges and corners of buildings for depths that are small 

relative to dimensions of the building. At such depths under the
THREE-DIMENSIONAL HEAT CONDUCTION IN PERMAFROST 65

edge of a building, the seasonal temperature variation would be 

approximated by that produced by linear heat flow under an undisturbed surface whose seasonal temperature variation had an amplitude of % (A-\-O). The surface amplitude in this case is the average 

of that inside the building and that outside. Under a square corner 

the effect would be approximated by that produced in linear heat 

flow by a surface temperature of amplitude % (3A-\-C).

The important special case of no seasonal temperature variation in 

S is obtained by setting (7=0. Thus the seasonal temperature effect 

under the center of a circular region of constant temperature is seen 

from equation 41 to be

/2*-KR «-»"- -V- 2asin ( a>^~~ o,f-V?+B* ,/£ («>

Corresponding to equation 33, the seasonal effect under the vertex 

of the sector of an annulus with central angle X, inner radius RI, 

and outer radius R2 , is seen from equation 38 to be

(43)

Equation 43 can be used to evaluate graphically the equilibrium 

seasonal disturbance at any point, resulting from a temperature 

D sin wt in an arbitrarily shaped region of the surface 2=0.

As in the case discussed in the previous section, the calculation is 

somewhat simpler, though less precise, if the integration with respect 

to fi is performed graphically with the aid of the grid of plate 1. 

The procedure is the same as that used to compute the principal 

disturbance, except that in place of figure 22, we now use plate 3, 

which gives the value of the integrand, SF, in equation 37 as a function 

of r2/a for four values of u>t. Inasmuch as the seasonal disturbance is 

sinusoidal, with a period of 1 year at any point, its evaluation for 

only 2 values of ut would be enough to determine its behavior throughout the annual cycle. The four curves are given in plate 3, so that 

additional points on the seasonal-effect-versus-time curve can be 

evaluated as a check.

As an extension of the example in the previous section, the equilibrium seasonal effect has been computed for positions 1 and 2, on the 

assumption that the temperature in the disturbed region is independent of time and that the temperature of the surface outside the region 

varies as A sin a>t. This seasonal effect is obtained by adding to GO
66 EXPERIMENTAL AND- THEORETICAL GEOPHYSICS.

(eq 22) the seasonal disturbance resulting from a temperature of 

 A sin ut in the disturbed region and zero outside (eq 37 with 

(7=0). This seasonal effect superposed on the equilibrium principal 

disturbance is shown graphically in figure 24. These results may be 

pictured as an approximation to the temperature beneath a surface 

supporting a heated building whose floor temperature does not vary 

seasonally and is everywhere B° C above the mean annual temperature of the surface outside. The numerical results in figure 24 are 

based upon an assumed diffusivity of 0.01 cm2 per sec. The bottom 

curve in this illustration represents the temperature to be expected 20 

feet beneath an undisturbed surface. A comparison of the curves 

representing the undisturbed case and position 2 in figure 24 reveals 

that the "building" has the effect of raising the mean annual temperature about 0.2B, diminishing the amplitude of seasonal variation 

about 20 percent, and causing the seasonal temperature wave to lead 

the undisturbed wave by a phase angle corresponding to about 1 week. 

At position 1 the building has the effect of raising the mean annual 

temperature about QA5B, diminishing the amplitude by more than 85 

percent, and causing the seasonal wave to lag the undisturbed wave 

by almost 2 months. As an example of the orders of magnitude involved, both A and B are close to 15° C for most heated buildings in 

the vicinity of Barrow, Alaska.

Cases in which the amplitude of seasonal variation of temperature 

in the disturbed region varies with position may be treated by the 

method discussed for variable B (p. 62).

Cases in which the phase of sinusoidal variation of surface temperature is different inside and outside the disturbed region may be treated 

by superposing the seasonal disturbance caused by a temperature 

D'sm ut in S on that of a temperature D" cos «£, where D' and D" are 

adjusted to give the required phase difference and amplitude.

EFFECT OF IDEALIZING THE PROBLEM

In order to obtain the above results it was necessary to assume that 

the medium was homogeneous and to neglect the effects of latent heat 

absorbed (or liberated) by the change of state of ice (or water).

In many engineering problems connected with permafrost, a thermal 

disturbance of the surface results in a descent of the freezing isotherm 

and the consequent thawing of frozen ground. Because the heat required to melt the ice causes no change in temperature, the establishment of equilibrium conditions is retarded, and because the rate of 

descent of the freezing isotherm is, in general, a function of position, 

the amount by which the buildup of the temperature disturbance is 

retarded is also a function of position. When equilibrium is finally 

reached, however, the freezing isotherm is stationary and no latent
THREE-DIMENSIONAL HEAT CONDUCTION IN PERMAFROST 67

+0.48

+.44 -

+.40 -

Principal disturbance = 0.4565

Principal disturbance = 0.1945

FIQUBE 24. Calculated temperature deviation from the undisturbed mean for large values of time at depth 

of 20 feet for C=0, a=*0.01 cm" per see. JB assumed equal to A to simplify presentation.
68 EXPERIMENTAL AND THEORETICAL GEOPHYSICS

heat is absorbed or liberated. Thus, if the effect of latent heat alone is 

considered, the equilibrium principal disturbance would be unaltered 

but the establishment of equilibrium would be somewhat slower than 

that described above (for example, in figure 23). The shape of the 

isotherms in the transient stages would also be altered. Latent heat 

will alter the seasonal effects by diminishing the amplitude and introducing some asymmetry for points near the freezing isotherm. It will 

be observed that for most problems of foundation engineering, neglecting the effects of latent heat provides a margin of safety.

Effective homogeneity is rarely found in nature, and in general the 

ground is stratified with respect to thermal properties. For purposes 

of calculating the buildup of the principal disturbance in time or the 

equilibrium seasonal effects, an engineering approximation can probably be obtained in most cases by assuming an average value for the 

diffusivity and by using the results determined for the homogeneous 

case. Bracketing situations can commonly be obtained by assuming 

extreme values of diffusivity. It will be observed from equation 31 

that the steady-state principal disturbance is independent of thermal 

properties, provided the material is homogeneous. In a stratified 

medium this is not so, although the homogeneous case will still probably give results satisfactory for engineering purposes in most cases.

SUMMARY

With the above results it is possible to evaluate the transient 

principal temperature disturbance and the equilibrium seasonal effects 

in a semi-infinite homogeneous medium in which surface temperature 

varies sinusoidally, but with different mean temperatures and amplitudes inside and outside a surface region of arbitrary shape. Two 

graphical methods are available. The first involves summing integrated contributions to the temperature; the second, a graphical 

integration with respect to solid angle. The first method is more 

precise but takes a little more time to apply. Since these methods 

would probably be primarily of interest to engineers who would use 

them to obtain approximate results, the second was discussed in 

greater detail. In the application of these methods to problems of 

foundation engineering in permafrost areas, some uncertainty is introduced by neglecting the effects of stratification and latent heat. 

The latter, however, provides a margin of safety in cases of thawing 

permafrost.

Although the discussion in this paper is devoted primarily to the 

thermal disturbance caused by a heated building, the methods described apply directly to several other important thermal problems 

in permafrost. For example, the geothermal disturbance caused by 

bodies of water in high latitudes is similar to that caused by a heated
THREE-DIMENSIONAL HEAT CONDUCTION IN PERMAFROST 69

building in that the mean annual temperature beneath the water is 

anomalously high and the amplitude of seasonal temperature variation is low. The thermal disturbance caused by the presence of a 

roadway usually results in a slight shift in the mean annual temperature but a large increase in the amplitude of seasonal variation. The 

increased amplitude causes a thickening of the active layer beneath 

the roadway, which frequently results in destruction of the road. 

Another useful application of these methods is to predict the extent 

to which ground will freeze beneath cold-storage vaults in temperate 

regions.

LITERATURE CITED

Bierens de Haan, David, 1939, Nouvelles tables d'inte"grales de"finies: New York,

G. E. Stechert and Co., 716 p. 

Birch, Francis, 1950, Flow of heat in the Front Range, Colo.: Geol. Soc. America

Bull., v. 61, p. 567-630. 

Carslaw, H. S. and Jaeger, J. C., 1947, Conduction of heat in solids: Oxford

Clarendon Press, 386 p.