Finite-Difference Model Predictions

of the Drastic Retreat

of Columbia Glacier, Alaska

y-vz Ij**

^sa 

^~-~«»£5

GEOLOGICAL PROFESSIONAL'PAPER 1258-D
Finite-Difference Model Predictions of the Drastic 

Retreat of Columbia Glacier, Alaska

By R. A. Bindschadler and L. A. Rasmussen

STUDIES OF COLUMBIA GLACIER, ALASKA

GEOLOGICAL SURVEY PROFESSIONAL PAPER 1258-D

Discussion of a one-dimensional dynamic model 

developed to predict the retreat and the iceberg 

discharge from a large iceberg-calving glacier

UNITED STATES GOVERNMENT PRINTING OFFICE, WASH I N GTON : 1 983
UNITED STATES DEPARTMENT OF THE INTERIOR 

JAMES G. WATT, Secretary

GEOLOGICAL SURVEY 

Dallas L. Peck, Director

Library of Congress Cataloging in Publication Data

Bindschadler, R. A. (Robert A.)

Finite-difference model predictions of the drastic retreat of Columbia Glacier, Alaska.

(Studies of Columbia Glacier, Alaska) (Geological Survey professional paper; 1258-D)

"Discussion of a one-dimensional dynamic model developed to predict the retreat and the iceberg discharge from a large icebergcalving glacier."

Bibliography: p.

Supt. of Doc. no.: I 19.16:1258-0

1. Columbia Glacier (Alaska) Mathematical models. 2. Finite differences. I. Rasmussen, L. A. II. Title. 

III. Series. IV. Series: Geological Survey professional paper; 1258-D. 

GB2425.A4B561983 551.3'12'097983 83-600208

For sale by the Distribution Branch, U.S. Geological Survey, 

604 South Pickett Street, Alexandria, VA 22304
CONTENTS

Page

Symbols and abbreviations 

Abstract 

Introduction 

Acknowledgments 

Previous modeling of Columbia Glacier - 

Finite-difference model 

Parameterization of glacier geometry 

Continuity equation 

Glacier dynamics 

Base shear stress 

Deformation velocity 

Sliding velocity 

Volume flux 

 iv

 1

 1

 3

 3

 3

 3

 4

 4

_ 4

 4

 5

 5

Finite-difference equationBoundary conditions

Temporal behavior of terminus 

Advance or retreat 

Calving relationship 

Prevention of floating terminus 

Summary of variables and parametersColumbia Glacier data 

Prediction of Columbia Glacier behavior - 

Water law versus ice law 

Effect of floating ice -

Comparison of predictions with observed retreatReferences 

Page

- 5

- 6

- 7

- 7

- 7

- 8

- 9

- 9

- 10

- 10

- 12 

_ 12 

_ 16

ILLUSTRATIONS

Figure 1. Map of Columbia Glacier ---- -------- 

2. Schematic of finite-difference approximation of a longitudinal section of a glacier 

3. Schematic depicting modeled representation of calving terminus 

4. Longitudinal profiles of important variables for Columbia Glacier 

5. Predictions of terminus retreat for various calving parameterizations 

6. Predicted terminus retreat using modified water law of calving for various values of c and 6 

7. Predicted terminus retreat with modified water law of calving for various treatments of floating ice and values of c -

8. Profiles of predicted thickness change since 1978.2 for c=16.9 a" 1 and 6=0 

2

6

8

11

12

13

14

15

TABLE

Table 1. Initial retreat rates for various values of coefficient in water and ice calving laws 10

CONTENTS III
SYMBOLS AND ABBREVIATIONS

Symbol Name

A Coefficient in ice-flow law

a Year

b Mass balance (ice equivalent)

c Coefficient in calving law (water law)

c' Coefficient in calving law (ice law)

f Velocity shape factor

f* Flux shape factor

g Gravitational acceleration

hw Water depth at terminus

H Ice thickness (vertical)

H* Ice thickness at terminus

i Gridpoint index

m Time step index

n Exponent in ice-flow law

P Coefficient of transverse channel shape

Q Volume flux

Q Calving flux

R Coefficient of transverse channel shape

S Channel cross-section area

t Time

Af Time step

Unit

bar"" a" 1

m a" 1

aa'1

dimensionless

dimensionless

m s"2

m

m

m

dimensionless

dimensionless

dimensionless

m1'2

m3 a"

m3 a'1

dimensionless

m2

a

a

Symbol

u

V

V

w

X

X*

Y

a

ft

8

f.

<t>

X

Q

Qw

T

T'

T*

r*

e

Name Unit

Deformation velocity m a"1

Total velocity m a' 1

Calving velocity m a"1

Channel width m

Longitudinal coordiante m

Terminus position m

Bed elevation m

Surface slope dimensionless

Thickness gradient dimensionless

Advance of terminus m

Strain rate a"1 

Weighting factor for effective stress dimensionless 

Ratio of sliding to total velocity dimensionless

Density of ice kg m"3

Density of water kg m"3

Effective base shear stress N m"2

Effective deviatoric stress N m"2

Base shear stress N m"2

Average base shear stress N m"2 

Weighting factor in calving law dimensionless

IV SYMBOLS AND ABBREVIATIONS
STUDIES OF COLUMBIA GLACIER, ALASKA

FINITE-DIFFERENCE MODEL PREDICTIONS OF THE DRASTIC 

RETREAT OF COLUMBIA GLACIER, ALASKA

By R. A. BiNDSCHADLER1 and L. A. RASMUSSEN

ABSTRACT

A finite-difference model of Columbia Glacier is used to study the 

possible retreat of the terminus. The model calculates the glacier motion from the glacier geometry via relationships of glacier dynamics 

and updates the glacier geometry in time increments via the continuity equation and a prescribed mass balance function. The sliding 

velocity is assumed to have the same functional relationship to 

glacier geometry as does the velocity due to deformation. Values of 

model parameters are fit to Columbia Glacier data gathered during 

the 1977-78 mass balance year. Approach to flotation of the terminus is implied by the model as the cause of the predicted rapid 

retreat and is a direct consequence of the rapid thinning observed at 

the terminus. Alternative schemes are examined to treat a floating 

terminus. In all cases, a rapid retreat eventually occurs.

INTRODUCTION

Columbia Glacier is one of the larger Alaska 

glaciers. Located a few kilometers west of Valdez, the 

glacier covers a 1,100-km2 area with the longest 

tributary stretching more than 66 km (fig. 1). Many 

large tributaries coalesce to form the lower trunk of 

the glacier, which fills a deep fiord and ends in a calving terminus towering 90 m above sea level at the head 

of Columbia Bay. Although the altitude of the bed 

beneath the glacier reaches depths of more than 400 m 

below sea level, it average less than 100 m below sea 

level at the terminus. This is because the terminus 

rests on a high ridge of material deposited by the 

glacier itself as it advanced into what was once a larger 

Columbia Bay (Post, 1975). Similar shoal features are 

found at the outlets of other glacially scoured fiords in 

Alaska (Post, 1980a-f), but Columbia Glacier is unique 

as the only glacier which still extends out to its 

neoglacial shoal.

1Code 912.1, Ice Station, Goddard Laboratory for Atmospheric 

Sciences, NASA/Goddard Space Flight Center, Greenbelt, MD 

20771.

Considerable evidence now suggests that the terminus of Columbia Glacier, in near equilibrium since 

the first recorded observation in 1794, is on the verge 

of a rapid and drastic retreat from its shoal (Post, 1975; 

Sikonia and Post, 1979). An intensive study of Columbia Glacier, begun in 1977 by the U.S. Geological 

Survey, has documented an accelerated rate of thinning in the terminus region and a trend of increasingly 

larger embayments being formed as an ever greater 

fraction of the terminus retreats from the crest of the 

shoal each summer. These partial retreats are believed 

to be precursory to the drastic, wholesale retreat soon 

to commence when the small fraction of terminus resting in shallow water can no longer prevent irreversible 

embayment growth. Observations of calving glaciers 

(grounded glaciers which terminate in tidal waters and 

discharge icebergs) have established that stable or 

nearly stable terminii occur at water depths of less 

than 80 m and that, when terminii have had to retreat 

through deeper waters, the retreats have been very 

fast (Meier and others, 1980).

A direct consequence of a rapid retreat is an increased 

flux of icebergs. In the case of Columbia Glacier, the 

presence of icebergs in Columbia Bay is a concern owing to the proximity of the shipping lanes in Valdez 

Arm, which is less than 20 km from the terminus (fig. 

1), used by tankers transporting Alaska pipeline oil 

stored in nearby Valdez. Even at current ice discharge 

rates, small icebergs drift into the shipping lanes 

(Kollmeyer and others, 1977). At increased calving 

rates, more numerous, larger icebergs might create a 

hazard to shipping.

Because of this potential hazard, the U.S. Geological 

Survey initiated its investigation of Columbia Glacier 

(Meier and others, 1978). Details of each phase of this 

study are being published as a series of scientific 

papers, of which this is one. Those already published 

discuss the relation between iceberg calving and water
2, STUDIES OF COLUMBIA GLACIER, ALASKA 

147°00'

 61°00'

10 KILOMETERS

FIGURE 1. Map of Columbia Glacier. Arrows indicate direction of ice flow. Largest tributary is marked every 2 km from ice divide to 

terminus. Modeled region, beginning at the 53-km mark, is hachured.
PARAMETERIZATION OF GLACIER GEOMETRY

depth (Brown and others, 1982); the adjustment and 

interpolation of point data (Sikonia, 1983); a finiteelement model of the lower glacier (Sikonia, 1982); and 

a continuity-equation model of the lower glacier 

(Rasmussen and Meier, 1982).

This paper discusses the results of applying a finitedifference numerical model of glacier flow to the lower 

reach of Columbia Glacier (fig. 1) to predict its future 

behavior. Modeling a tidewater glacier such as Columbia Glacier presents a new and challenging problem in 

glacier dynamics. A major aspect of this program is the 

time-dependent behavior of the terminus. This modeling is also important because of the potential hazard of 

discharged icebergs to shipping. The parameterization 

of calving was an important element of each model 

(this finite-difference model as well as the continuity 

and finite-element models referenced above) but was 

limited by knowledge of the calving process and its controls. Thus, any model is restricted by gaps in our 

understanding and can only predict a range of possible 

future behaviors. Available data help reduce the range 

of the predictions, but further refinements can only 

result from future observations. Nevertheless, it is encouraging that, although each of the three numerical 

models took different approaches and made different 

assumptions, the general results were similar.

ACKNOWLEDGMENTS

The authors are grateful to Dr. R. H. Thomas and 

Dr. E. D. Waddington for carefully reviewing the original version of the manuscript and for making numerous 

constructive suggestions.

PREVIOUS MODELING OF COLUMBIA GLACIER

The continuity model of Rasmussen and Meier (1982) 

incorporated a minimum of glacier dynamics into the 

predictions of terminus retreat. By requiring mass conservation, the time intervals between a set of sequential elevation profiles were calculated. The series of 

profiles was interpolated between the measured elevation profile existing in 1978 and an hypothesized elevation profile corresponding to a 14-km retreat of the terminus. The model was one-dimensional, simulating 

the retreat along a central flowline, and was independent of any transverse variations in any of the parameters. The calving flux in this model was specified to 

depend linearly on the water depth at the terminus. 

The details of this relationship are discussed in a later 

section. The results of this continuity model predicted 

that the terminus will begin a drastic retreat in the 

1980's. The speed of the predicted retreat depended on

the value of a parameter in the calving term but was 

on the order of 4 km a" 1 , with a maximum calving flux 

of about 10 km3 a'1 (compared with about 1 km3 a"1 in 

1978).

The finite-element model of Sikonia (1982) analyzed 

the flow of a vertical section along a central flowline 

over the lower 14 km of the glacier. A thin layer of elements with a very low viscosity at the base of the 

glacier was used to simulate basal sliding. The calving 

flux in this model was related to the ice thickness at 

the terminus unsupported by buoyancy and included a 

superimposed seasonal variation due to the apparent 

dependence of calving on freshwater runoff. The 

results of this modeling were extremely sensitive to 

the precise values of the calving law parameters; seemingly minor variations in these parameters caused 

variations ranging from no retreat to drastic retreat 

before 1980.

FINITE-DIFFERENCE MODEL

The finite-difference model, which is applied here to 

Columbia Glacier, was originally developed as a general 

glacier model capable of simulating the dynamics of 

valley glaciers whose motion is predominantly due to 

internal deformation of the ice. It was first applied to 

predict the behavior of a surge-type glacier during the 

quiescent phase of motion (Bindschadler, 1982). The 

model is described in Bindschadler (1982), and additional detail can be found in Bindschadler (1978); what 

will be presented here is a summary of the approach 

used and the assumptions made in applying this model 

to Columbia Glacier.

PARAMETERIZATION OF GLACIER GEOMETRY

The model is one-dimensional; the profile of glacier 

depth, H(x\ along a central flowline (x axis) is the principal dependent variable. The transverse channel cross 

section is specified by the equation,

W(x,z)=P(x)zl'2+R(x)z, (1)

relating the channel width, W, to the height, z, above 

the glacier bed, and describing, the shape of the channel, P and R. Equation 1 represents a linear combination of channels with either a parabolic or V-shaped 

cross section, the shapes most common for mountain 

glaciers. The cross section area is obtained by integrating equation 1 over the depth:

(2)
STUDIES OF COLUMBIA GLACIER, ALASKA

For a known set of widths, cross-section areas and 

depths, equations 1 and 2 can be used to determine the 

values of P and R at each x. Once the parameters are 

determined and the bed elevation, Y(JC), is specified, 

changes in the geometry of the glacier depend only on 

changes in H(x). Thus, in strict mathematical terms, 

the model is one-dimensional; however, equations 1 

and 2 parameterize additional dimensions which affect 

the dynamic behavior of the glacier. The effect of a 

channel with sloping sides is to slow the speed of the 

glacier's response to changing conditions (Bindschadler, 1978, p. 127).

CONTINUITY EQUATION

For each time step, the ice motion is calculated from 

the glacier geometry and is used to update the depth 

profile through the vertically integrated continuity 

equation,

dS(x) 

dt (3)

where Q(x) is the volume flux along the #-axis, #(*) is 

the mass balance, W(x) is the surface width, and t represents time. Each of the variables in equation 3 must 

either be specified or calculated from the geometry.

GLACIER DYNAMICS

BASE SHEAR STRESS

Gravitational forces create stresses within the ice 

which cause it to deform. This deformational motion of 

the ice is one component of the ice velocity; the other is 

basal sliding and is discussed later. If it is assumed 

that the top and bottom surfaces of the glacier are 

parallel, that there are no longitudinal stress gradients, 

and that the glacier is infinitely wide, then the flow is 

laminar and the shear stress at the glacier bed is

T*(X)=QgH(x)cosa(x)8m<x(x), (4)

where Q is the ice density, g is the gravitational acceleration, a is the surface slope, and H(x)coaa(x) is the surface 

normal depth.

On Columbia Glacier, the surface and bed are nearly 

parallel, but longitudinal stresses vary along the 

width, and the channel is not infinitely wide. However, 

equation 4 can be modified to partially account for 

these effects. Longitudinal stress variations occur 

whenever the longitudinal gradient of velocity is not 

constant. By averaging the surface slope in equation 4 

over distances of 10 to 20 times the glacier depth, the

effects of longitudinal stress gradients are greatly 

smoothed, permitting a more accurate calculation of 

the longitudinal velocity profile (Budd, 1970; Meier 

and others, 1974). In this model, the entire right-hand 

side of equation 4 will be averaged over these 

distances, but, because the variation of a(x) is much 

larger than the variation in H(x), the differences between the two averaging methods are sjnall. These 

averaged base stresses are denoted T*, and the 

modified base shear stress is

r*(*)= (5)

A second modification to the basal shear stress takes 

account of the friction of the moving ice against the 

valley walls, which causes additional shear within the 

ice. The magnitude of this shear is less than at the 

centerline where the maximum depth is assumed to occur. Thus, the mean shear stress in the ice is some fraction, f< 1, of the shear stress in equation 5. Nye (1965) 

has calculated this fraction (referred to here as the 

"velocity shape factor") for various channel shapes. 

His calculations show that fis insensitive to changes in 

ice depth for a given cross-section shape.

A final modification to equation 5 was required to 

preserve numerical stability in the finite-difference 

solution (a discussion of the problem of numerical stability in glacier modeling can be found in Waddington, 

1981, Section 2.2.4). Values of T* were combined with 

unaveraged T*', the linear combination of these two 

stresses formed an "effective" average base stress,

(6)

where 0 is a positive weighting factor, less than unity. 

For values of <j> close to unity, this final modification 

has little effect on the overall dynamics of the glacier. 

A detailed study of the numerical stability showed that 

0 depended slightly on the averaging distance in T* but 

that 0=0.8 always ensured numerical stability (Bindschadler, 1978, Section 5.6.4).

DEFORMATION VELOCITY

The rate of deformation in the ice that results from 

an applied stress is usually expressed as the following 

nonlinear flow law:

e=A(rO" (7)

(Nye, 1957), where e is the effective strain rate, T' is the 

effective deviatoric stress, and n and A are empirically 

determined parameters. Using equation 6 and integrating over the depth, the ice deformation velocity at
FINITE-DIFFERENCE EQUATIONS

the glacier surface due to the effective average shear 

stress, T, is

SLIDING VELOCITY

The second component of ice velocity is the sliding 

velocity. Unfortunately, the understanding of this 

process is very poor. Expansive theoretical developments on this subject fall far short of being able to 

predict sliding velocity from glacier geometry (Meier, 

1968). Field measurements of Columbia Glacier show 

that sliding is a major fraction of the velocity near the 

terminus (W. G. Sikonia, written commun., 1979); 

therefore, some parameterization of sliding was required for the model. In lieu of using complex schemes 

of sliding with no sound physical basis, it was felt most 

expedient to use a simple parameterization which includes the known important variables of sliding. The 

most important variables controlling sliding are basal 

shear stress, bed roughness, and either subglacial 

water pressure (Bindschadler, 1983) or subglacial 

cavity growth (Iken, 1981). For the model, bed roughness, water pressure, and the growth rate of subglacial 

cavities were taken as constant in time at each point 

in space, but changes in the basal shear stress were calculated from equations 4, 5, and 6. The simplest conversion of the base shear stress to sliding velocity is to 

assume the same form of relationship as equation 8, including the same value for the exponent. There is additional justification for assuming a nonlinear sliding 

law; various theoreticians (Kamb, 1970; Weertman, 

1964; Lliboutry, 1968; Morland, 1976a, b) have derived 

such relationships, but there is some disagreement on 

whether the value of the stress exponent is the same n 

as in equation 7. This assumption is equivalent to postulating that the ratio of the sliding velocity to total 

surface velocity is a constant at each point in space. 

Thus, a parameter X is defined as

V(x)-U(x) 

V(x) (9)

where V(x) is the total surface velocity. Values of \(x) 

were determined from the field data by considering the 

annual variation of velocity and assuming the winter 

minimum velocity to be due entirely to deformational 

motion.

VOLUME FLUX

Finally, the velocity must be converted to a volume 

flux. Because the surface velocity at the channel center 

is the maximum velocity over the entire cross section, 

a final parameter, the "flux shape factor," /*^1, is 

defined as the ratio of the deformation velocity averaged over the entire cross-section area, S, to the maxinum deformation velocity, U (eq 8). By assuming that 

the sliding component of motion is constant over the 

ice-rock interface in the section, the volume flux can be 

written

Q(x)=[f*(x)U(x)+ V(x)- U(x)]S(x). (10)

As with the velocity shape factor, the flux shape factor 

is not sensitive to thickness changes (Nye, 1965); once 

the values are determined for /"and /"*, they are held 

constant in time.

FINITE-DIFFERENCE EQUATIONS

In the finite-difference method, the terms of the continuity equation (eq 3) must be represented by finite 

differences on a discrete grid. The gridpoints occur at #, 

(l<i<p); depth, H,; bed elevation, Y,; cross- section 

area, S,; channel width, W,; and mass balance, 6, are 

either specified or calculated at each xh and the volume 

fluxes, Qm/2, are calculated at the grid midpoints, xi+l/2 

(fig. 2). Thus, the finite-difference form of the continuity equation is

_/")m-H/2 - -i TTT = b' W"

m+l _

where Q!^y22 =j(Q!til/i +Q#i/2), and superscripts (12)

refer to the time step index m or (m+1).

Equations 11 and 12 are referred to as the CrankNicholson finite-difference scheme, where the term 

Q£t/22 is the average volume flux over the time step 

(McCracken and Dorn, 1964). The finite difference 

forms of the ancillary equations are

Wt=Pt

_2A 'n+l cosam/2

(13)

(14)

(15)

(16)
STUDIES OF COLUMBIA GLACIER, ALASKA

LLJ 

> 

LLJ

LLJ 

CO

LLJ

>

o

CD

o

hLLJ

_J 

LLJ

X;

DISTANCE FROM ICE DIVIDE

1+1

FIGURE 2. Schematic of finite-difference approximation of a longitudinal section of a glacier. Glacier is represented by a bedrock elevation, 

Yi, and ice thickness, H» at each gridpoint, x,. Continuity equation at x, is applied over the hachured region with mass balance, 6,, 

specified at x, and volume fluxes, Q.-i/2 and Q.+I/I, calculated between xt.i, xlt and x,+1 .

 "« ! .

sinaI+1/2

. a,+ i/2 =tan

(17)

/.. o\ (18)

(19)

where xk+ i-xh is the large-scale averaging distance 

centered at x^1/2 . As xi+1/2 approached either end of the 

grid, the averaging distance necessarily decreased but 

remained centered on JcI+ i/2 .

BOUNDARY CONDITIONS

Equations 11 through 19 represent a second-order 

system; thus, their solution requires two boundary conditions. The scheme used in the Columbia Glacier 

calculations was to specify volume fluxes at the upstream and downstream ends of the modeled region. At 

the upstream boundary, a constant flux, Qi/2, was 

specified that was approximately equal to the observed 

flux at this boundary. The use of this boundary condition assumes that the glacier behavior upstream of the 

boundary is unaffected by the response within the 

modeled region. For the modeling of Columbia Glacier,
ADVANCE OR RETREAT

this is probably a reasonable assumption over the short 

term (5-10 years) but limits the accuracy of predictions 

made over longer time scales.

At the downstream boundary, the terminus, a condition of constant flux is unreasonable. This is the most 

active region of the glacier during a retreat, so a boundary condition must be established that links this 

downstream volume flux to the glacier geometry. The 

method used to satisfy this second condition is to calculate a volume flux at the last gridpoint along the 

flowline and to use a backward-difference approximation for the flux gradient term in equation 11:

-t-l/2 -/^m+1/2

, _i w" "'

where xp is the terminus gridpoint (fig. 3). The other 

modified equations are

-]UP sp

ffHcosa

Tp =fpQgHpsinap

(15*)

(16*)

(17*) 

(19*)

No large-scale average of T* is used in equation 17* 

because at the endpoints of the grid any centered largescale average reduces to an average over the two gridpoints nearest the boundary.

This system of equations can be solved for the glacier 

depths at any time, given the geometry at any earlier 

time. Because the Crank-Nicholson approximation is 

implicit (that is, the Hr+1 cannot be solved for directly), 

an iterative scheme was employed to converge to the 

unknown Hr+1 at each time step. A Newston-Raphson 

technique was used, McCracken and Dorn (1964) give a 

general description of the method, and Bindschadler 

(1978) details its implementation in this model.

TEMPORAL BEHAVIOR OF TERMINUS

ADVANCE OR RETREAT

As mentioned earlier, the central issue in predicting 

the behavior of Columbia Glacier is predicting the 

time-dependent position of the terminus and the rate of 

calving. To calculate the new position of the terminus, 

once again, continuity is invoked, this time within the 

final segment of the glacier represented by the volume 

between the terminus positions at the beginning and

end of any time step, A£(fig. 3). The volume of this segment is

(20)

where Q is the calving flux (discussed below). The surface slope from JCP_I to xp is assumed for any ice advancing beyond xp, and the advance of the terminus calculated as

(21)

where /3=:

 D *^D 1

Equation 21 assumes a parallel-sided section of advancing or retreating ice with length 6 and width Wp, 

within the nearest few hundred meters behind the 

shoal, the parallel-sided assumption is reasonable. The 

quantity 6 is positive for an advance and negative for a 

retreat. The last gridpoint is then moved to the new 

terminus position, X*, where the new ice depth is

(22) "** i

CALVING RELATIONSHIP

Equations 20 and 21 require that the flux of ice calving from the terminus, Q, be known. The development 

of this relationship is the subject of a separate paper in 

this series on Columbia Glacier research (Brown and 

others, 1983). In summary, based on data from 13 

temperate calving glaciers, the calving relationship 

that provides the best empirical fit to the data is

V=ch (23)

where V is the calving speed, hw is the water depth at 

the center of channel, and c is a constant. Equation 23 

will be referred to as the "water law." The best fit, 

providing a 0.84 coefficient of determination, is for 

c=16.94 a'1 . As established by Post (1975), an irreversible, rapid retreat ensues when the terminus of a calving glacier enters deep water behind its shoal. An 

alternative formulation of calving speed, but with a 

smaller coefficient of determination, is

V=c'tfp, (24)

where Hp is the ice thickness at the terminus instead of 

the water thickness. Equation 24 will be referred to as 

the "ice law." For either law of calving, the corresponding calving flux is
STUDIES OF COLUMBIA GLACIER, ALASKA

FIGURE 3. Schematic depicting modeled representation of calving terminus. Volume flux, Qp , is calculated as described in text, and continuity equation applied at terminus is determined from calving flux (Q), QP, and continuity as described in text.

Q=vsp. (25)

Clearly equation 25 does not account for the strong seasonal variations in calving that have been measured 

on Columbia Glacier (Sikonia and Post, 1979); it quantifies the mean calving flux for a grounded terminus.

PREVENTION OF FLOATING TERMINUS

Although much of the physics involved in the calving 

process remains a mystery, there are still some 

physical constraints which can be incorporated into the 

model. The absence of any temperate ice shelves in the 

world strongly suggests that floating temperate ice 

quickly disaggregates into icebergs. This phenomenon 

was not incorporated into the empirical determination 

of the calving relationships (eqs 23, 24) because none of 

the tidewater glaciers studied had a floating or nearfloating terminus. Nevertheless, the model's treatment of the motion of the terminus must be modified to

take account of an expected dramatic increase in calving when the terminus is near flotation. Two alternative treatments were included in the model. The first 

was to locate the new terminus just upstream of any 

floating ice and to remove the floating ice from the 

glacier. This appeared rather arbitrary, however, and 

ignored the possible dynamic effects as the terminus 

approaches flotation.

A second, more involved, method was to modify the 

water law (eq 23) as follows:

V= 9

1 (26)

where Hp and hw are the ice and water thicknesses, 

respectively, at the terminus; e» is the density of water 

at the terminus; and Q is a parameter that determines 

how sensitive calving is to a near-floating situation. 

For 0=0, the calving speed is independent of buoyant
COLUMBIA GLACIER DATA

forces, and equation 26 reduces to the water law. To 

avoid a sudden transition to a floating terminus, the 

parameter 6 had to be set to some small positive value, 

but no data were known to exist which might be used 

to estimate a reasonable value for 6.

Equation 26 does not ensure that at the end of a 

given time step there would not be some floating section of ice at the terminus. For larger values of 6, a 

floating situation was less likely to occur (a fact that 

possibly can be used as an ad hoc guide to reasonable 

values of 0), but the model still needed to be equipped to 

deal with this eventuality of floating ice within the 

final segment of the glacier. Rather than remove any 

floating volume of ice remaining at the end of a time 

step, it was redistributed over the surface of the final 

glacier segment (fig. 3). It was assumed that this redistribution was accomplished by some unprescribed 

dynamics during the most recent time step. By redistributing this extra mass in the shape of a wedge, the 

thickness at the gridpoint adjacent to the terminus remained unchanged, and a reduced surface slope compensated for the increased thickness at the terminus. 

This scheme rarely needed to be used and then resulted 

in only slight changes to the volume flux at the terminus. These special features of this particular mass 

redistribution maintained a relatively smooth 

behavior of the dynamic response of the glacier in the 

terminus area.

SUMMARY OF VARIABLES AND PARAMETERS

In the description of the model, numerous variables 

and parameters have been specified as necessary inputs to the model. They can be classified into the 

following groups: geometric, physical, and dynamic. In 

the geometric class are included the arrays of gridpoints, jc,; bed elevation, Y» ice thickness, H» and transverse channel shape, P, and Rt (eqs 1, 2). The physical 

parameters include the constants of gravitational acceleration, g=9.8 m s"2 ; ice density, e=910 kg nr3 ; and 

the ice-flow law parameters, n and A (eq 7). More generally, A is a function of ice temperature, ice fabric, 

and even impurity and water content. The final class 

includes the mass balance, 6, at each gridpoint; the 

velocity shape factors, /|+ i/2; flux shape factors, /m/2; 

and sliding ratios, Xm/2 specified at the midpoints of 

the grid. Also included in this class of dynamic parameters are the constants of upstream volume flux, Ql/2', 

calving flux coefficient, c; averaging length, xk -xh', and 

weighting factor, <j> (eq 6). With these parameters and 

variables chosen to represent the geometry and 

dynamics of a given glacier, the model predicts the 

future behavior of the glacier. The accuracy of the 

prediction rests on how precise the parameterization

can be made and on the assumption that all of the important physical processes have been incorporated.

COLUMBIA GLACIER DATA

The model is general and may be applied to any 

glacier. What must be done for prediction of the timedependent behavior of Columbia Glacier is to choose 

the appropriate parameter values that will yield the 

most accurate model representation of Columbia 

Glacier and that will serve as an initial condition. The 

field data used were collected primarily during the 

mass balance year from September 1, 1977, through 

August 31, 1978 (Mayo and others, 1979).

Figure 1 shows the complex, dendritic structure of 

the tributary glaciers, all of which contribute ice to the 

calving terminus. Modeling the entire glacier system 

would be a difficult task for the numerical models 

available and would require considerable data adjustment to provide a smooth initial condition. Thus, this 

model concentrated on the lower trunk of Columbia 

Glacier, the final 14 km (fig. 1). This permitted a finer 

grid resolution and a more accurate prediction of the 

sequence of events just prior to and during the initial 

rapid retreat phase. The assumption of only minor influence of the glacier upstream of the modeled region 

on the terminus behavior is justified by the expected 

imminence of rapid retreat. The surface topography 

and spatial distributions of surface velocity, mass 

balance, thickness change, and bed topography resulting from the field program are given by Sikonia (1983), 

Fountain (1982), and Mayo and others (1979). Because 

the present model does not account for seasonal variations, annually averaged data are used; the middle of 

the principal data year is 1978.2, at which time the 

transversely averaged position of the terminus was 

X=jCi9 =66.60 km, as measured from the head of the 

main trunk glacier.

The parameters must be chosen so that the initial 

conditions are consistent with the observed data as 

well as internally consistent with respect to each of the 

equations used to describe the glacier's dynamic 

behavior. Otherwise, the model would rapidly redistribute the mass of the glacier, not as a realistic projection of the glacier's future behavior, but artifically, to 

reconcile the inconsistencies in the initial conditions. 

Because there are numerous sources of uncertainty, 

both in the equations and in the observed data, the 

parameterization and construction of initial conditions 

is an underdetermined problem.

Uncertainties in the equations arise from such 

sources as the approximation of a three-dimensional 

flow regime by using a one-dimensional formulation 

with shape-factors parameterizing the effect of the
10 STUDIES OF COLUMBIA GLACIER, ALASKA

glacier channel on the longitudinal flow (eqs 13-15, 

17); the representation of the longitudinal stresses by 

the slope-averaging scheme (eq 17); the assumption 

that the sliding component of motion follows the same 

functional form as the deformation component (eq 9); 

the treatment at the terminus of the calving 

mechanism (eqs 20-26); and the selection of particular 

numerical values for the flow law parameters (eq 7).

Uncertainties in the data include the bed topography, the spatial and temporal distributions of the 

fraction of the total velocity that is due to sliding, and 

the transverse averaging of the mass balance and the 

thickness change data, as well as the original observation error in all of the field data.

The underdeterminedness could be exploited by taking some of the variables to be determined exactly by 

the field data and then using the equations to calculate 

values for the other variables, but this could lead to unsatisfactory results. It could lead to values for the geometric or dynamic variables that might differ markedly from the field data. It could lead to values for the 

shape factors that might fictitiously represent features 

that could not be independently substantiated.

Instead, an error band was associated with the 

longitudinal profile of each variable, narrower for 

better-known variables, such as the surface elevation 

profile, and wider for less well-known variables, such 

as the bed profile. Then the profiles, beginning with 

the best estimates from the field data, were simultaneously adjusted, each within its own error band, so 

that they were smooth and free from extraneous features and were consistent through each of the equations.

A grid was established within 19 gridpoints equally 

spaced at intervals of 762.5 m. These gridpoints remained fixed in space with the exception of the last 

gridpoint, which moved with the terminus, as described 

in the previous section. The flow law parameters (eq 7) 

were set at A=0.14 bar~" a" 1 and n=3. These values for 

A and n, although not the result of any individual 

rheological study of ice (for example, Glen, 1955, or 

Nye, 1953), are typical of the values for temperate ice 

that have appeared in the literature. Any variation of 

A due to inhomogeneities of ice fabric or impurities is 

probably negligible. Figure 4 shows the adjusted profiles of the variables and represents a "best-fit" case. 

The lack of major features in the profiles is a result of 

the underdeterminedness of the fitting procedure.

One of the most important parameters remained outside this fitting procedure, the coefficient in either 

calving relationship (eqs 23, 24). The models of 

Rasmussen and Meier (1982) and especially Sikonia 

(1982) showed a marked sensitivity of the timing of 

rapid retreat to the calving relationship used. Thus, in

the simulations that follow for this model, the sensitivity of the predictions to the value of this coefficient 

is keen interest.

PREDICTION OF COLUMBIA GLACIER 

BEHAVIOR

For any prediction computer run, a calving relationship (eqs 23, 24, 26) had to be specified with the appropriate parameters (c, c', or c and 0), and a method 

identified for treatment of floating ice (bulk removal or 

wedge redistribution). Thus, a large number of computer runs were possible. The approach taken here was 

first to examine the predicted response of the terminus 

for the water and ice laws of calving without any treatment for floating ice, then to use the modified water 

law with both treatments for floating ice, and finally to 

compare these predictions to available data of terminus position and suggest a "most-probable" scenario 

for the future behavior of Columbia Glacier.

WATER LAW VERSUS ICE LAW

In using either law, the value of the coefficient must 

be determined. This is accomplished by specifying an 

initial retreat rate and, from the known geometry and 

terminus volume flux (Qp\ by solving equation 20 for c 

(water law) or c' (ice law). Table 1 gives the coefficient 

values for various initial retreat rates.

The predicted response of the terminus for the five 

cases of table 1 are shown in figure 5. In each case, the 

rate of retreat slows from the initial value with the 

more pronounced decelerations for the ice law. 

Analysis of terminus position data completed after this 

modeling study indicated that retreat rates of 800 m a' 1 

were typical for the maximum retreat rate period during midsummer to late fall; however, Rasmussen and 

Meier (1982) calculate the average retreat rate during 

the 1977-78 measurement year to be 47 m a'1 . Figure 5 

shows that the case W!(c=16.934 a"1 , very close to the 

preferred value of 16.94 a" 1 in eq 23) predicts an average retreat of 47 m in the first year. This result lends 

credence to this particular choice of water law and c in

TABLE 1. Initial retreat rates for various values of coefficient in 

water and ice calving laws

Initial retreat rate

Law 200 m 400 m aM 800 m a'1

Water: c 

Ice: c' 

14.739

8.817

16.934

10.130

20.000 a-'
WATER LAW VERSUS ICE LAW 11

CO 

m OC 500

8

w < -500 QJ

CO

i i r i I i i i i i i

DC 

< 

IU

>

DC 

LU 

Q_

LU 

O

LL 

O

CO 

DC 

LLJ HUJ

0.5

SURFACE

-0.05

-0.04

25 

0.03 c

<o

U0.02 

0.01

54 56 58 60 62 64 66 68 

DISTANCE FROM ICE DIVIDE, IN KILOMETERS

FIGURE 4. Longitudinal profiles of important variables for Columbia Glacier. Distances are measured from ice divide of main trunk glacier 

(fig. 1).
12 STUDIES OF COLUMBIA GLACIER, ALASKA

LU 

Q

1 1 1 1 1 1 1 1 1 1 1 1 1

lo

7 8 9 10 11 

TIME, IN YEARS

12 13 14 15

FIGURE 5. Predictions of terminus retreat for various calving parameterizations. Curves labeled lo and I, correspond to ice law of calving 

(eq 24) with c'=8.817 |f' for I0 and 10.130 a"1 for I,. Curves labeled W0 , W,, and W2 correspond to water law of calving (eq 23) with 

c=14.739, 16.900, and 20.000 a'1 , respectively.

the model and to the effective parameterization of the 

glacier dynamics by showing that the predicted retreat 

rates are controlled more by the dynamics of the 

glacier than by the initial retreat rate used.

EFFECT OF FLOATING ICE

The predictions of any of the curves in figure 5 are 

not reliable once the ice thickness at the terminus is 

near to the flotation thickness. Figure 6 illustrates 

how the behavior of the terminus is changed when the 

modified water law (eq 26) is used for various values of 

c and 6 and an initial retreat rate of 400 m a'1 . The 

curve labeled "Flotation ignored" corresponds to curve 

Wi in figure 5. In every case where flotation was

parameterized, the glacier progresses into a phase of 

very rapid retreat. This retreat could only be followed 

for a few hundred meters because, as the retreat accelerated, the computation time step shrank to a few 

minutes which lengthened the required computing 

time prohibitively.

Figure 6 also shows that for 0=0, the wedge redistribution of floating ice delays the timing of the drastic 

retreat from the case when the floating ice is removed. 

The delay, however, only amounts to a few months. 

The rapid retreat is dramatic because, once the terminus retreats into deeper water, so much of the ice is 

near the flotation thickness; continued thinning near 

the terminus causes increasingly larger volumes of ice 

to float.
u, 66.7 

O

Q

Hi 66.6 

O

66.5

O 

DC

CO 

DC 

HI

HI

O 66.4

| 66.3

H 

CO 

O 

GL

co 66.2

D 

Z

DC 

HI

66.1

COMPARISON OF PREDICTIONS WITH OBSERVED RETREAT 

 I I I I 

13

6=0.02 

=15.96 a'1

. **

Flotation

ignored 

 

6=0.01 

e = 16.42 a"1 r-'

-47 m a"' -

6=0

c = 16.9 a'1

1978 1979 1980 1981 1982 1983 1984 1985

YEAR

FIGURE 6. Predicted terminus retreat using modified water law of calving (eq 26) for various values of c and 9. Each case refers to an initial 

terminus retreat of 400 m a"1 . Dotted curve results from ignoring flotation (see curve W1( fig. 5); dashed curve results from a wedge 

redistribution scheme for floating ice; and all others result from removing any floating ice from terminus.

The other curves in figure 6 complete a consistent 

picture of the effect of different values for 6 on the timing of the rapid retreat. The general results are that 

the onset of rapid retreat occurs sooner, and, as 0 increases, the transition to rapid retreat is smoother. In 

each run with a different 0 value, the value of c was adjusted to preserve the same initial retreat of 400 m a" 1

An even more significant result is that, despite the 

variation in when the rapid retreat begins, the position 

of the terminus at that time is roughly the same in all 

of the cases shown. This suggests that even if the most 

reasonable value of 0 is not well known, the onset of the 

rapid retreat of Columbia Glacier might still be 

predicted by monitoring the position of the terminus.

When the time axis is expressed as calendar dates 

(fig. 6), £=0 corresponds to 1978.2 as this is the approximate midpoint of the 1977-78 balance year during

which most of the field data were collected. By including this dateline in figure 6, it can be seen that, if 

the model is valid, 6 must be very small because by late 

1981 the rapid retreat phase had not yet been observed.

COMPARISON OF PREDICTION WITH OBSERVED RETREAT

Figure 7 compares the observational time series of 

width-averaged terminus positions with many runs of 

the model using the water law for calving (c= 14.74, 

16.90, and 20.00 a" 1 ) and each alternative for treating 

floating ice (no adjustment, bulk removal, and wedge 

redistribution). Based on figures 5 and 6, the case of 

c=16.9 a" 1 and a wedge redistribution appears to be the 

most realistic parameterization, and, indeed, figure 7 

shows that the predictions for this case generally fit 

quite well within the observed seasonal fluctations of
66.8

66.1

1976 1978 1980 1982 1984 

YEAR

FIGURE 7. Predicted terminus retreat with modified water law of calving (eq 26 wherein 6=0) for various treatments of floating ice 

as follows: a, no flotation condition; b, wedge redistribution of floating ice; and c, removal of floating ice. Also shown are data 

minus position: solid line preceded model study, dotted line postdated completion of model study, and dots indicate measured 
COMPARISON OF PREDICTIONS WITH OBSERVED RETREAT 15

terminus position. Data confirmed that the retreat in 

fall 1980 was truncated by an abnormally cold and 

snowy August and that the pronounced retreat in summer 1981 finally ended in November when a readvance 

began (M. F. Meier, personal commun., 1981). Comparing those data to figure 6, it can be seen that the 

retreat measured in 1981 fits the case of c=16.9 a" 1 and 

0=0.01 well, except that in the model no recovery was 

possible.

Also from figure 7, it is clear that c=20.00 a" 1 is too 

high a value and that c=14.74 a' 1 is too low for the coefficient of the water law of calving. Again, although the 

time of onset of rapid retreat is sensitive to the value of 

c, the terminus position at that time falls within a relatively narrow range. This result might have significance in understanding more about the dynamics of

rapid retreat but would be difficult to use in practice 

due to the large seasonal fluctuation of terminus position and its variability across the terminus width when 

large embayments are present.

The predicted rate of thinning over the lower glacier 

for the preferred case c=16.9 a" 1 , 0=0 is shown in 

figure 8 up to the moment of drastic retreat. The rate of 

thinning at the terminus decreases slightly until flotation is approached. Observations indicate that the recent thinning which has occurred over the lower section of Columbia Glacier has been more uniform spatially than figure 8 indicates and that this thinning has 

been proceeding at an average rate of about 3.5 m a'1 

(M. F. Meier, personal commun., 1982). This might be 

due to the role of longitudinal stresses in the dynamics 

in this region or to a temporal change in the propor-

-10

CO 

DC 

LLJ -20

LLJ

^ -30

uf 

o

O

CO 

CO

O

-40

-50

-60

-70

-80

54 56 58 60 62 64 

DISTANCE, IN KILOMETERS FROM ICE DIVIDE

66 68

FIGURE 8. Profiles of predicted thickness change since 1978.2 for c=16.9 a'1 and 0=0. At 4.7 years (or date 1982.9^, flotation occurs at the 

terminus and rapid retreat begins.
16 STUDIES OF COLUMBIA GLACIER, ALASKA

tional contribution of sliding velocity, neither of which 

are explicity included in the model. In the model, the 

longitudinal gradient of thinning represents the 

glacier's adjustment to accommodate an increasing 

calving flux. With time, the surface slope of the glacier 

steepens to provide more mass at the terminus, but the 

required thinning is so extreme that approach to flotation and rapid retreat are unavoidable. If the volume 

flux at the terminus could be increased without increasing the surface slope, then flotation might not 

precede rapid retreat. An example of this might be if 

the sliding velocity were to increase. If this were the 

case, thinning would be more widespread and uniform 

but would decrease the ability of the glacier to increase 

the volume flux later. This process appears to be important in the Sikonia (1982) finite-element model when 

rapid retreat is predicted. In either case, because the 

increase of calving flux is so rapid for even a slight 

retreat (owing to the steep slope of the shoal on the upglacier side), the glacier must be able to respond quickly with a much-increased flux of ice to the terminus. If 

it cannot, then the retreat will quickly become irreversible.

Additional evidence that flotation may play an important role in rapid retreat can be drawn from studies 

of the embayments which have formed annually at the 

terminus in recent years (Sikonia and Post, 1979). The 

location of these embayments usually coincides with 

the position of the major subglacial drainage channel 

(Sikonia and Post, 1979). The pressure of the water in 

this channel has been calculated to be large enough to 

float the ice locally (Bindschadler, 1982) causing a 

localized catastrophic retreat which is impeded only by 

the stability of the surrounding ice. These embayments 

have been observed to form very quickly. Their size has 

increased annually, presumably because the terminus 

is thinner and, as a whole, is less able to retard the 

formation of expanding embayments.

The finite-difference model, however, cannot model 

this embayment behavior because the growth is caused 

by seasonal increases in the calving flux and their extent is limited by altered transverse stresses within 

the more stable surrounding ice. Nevertheless, the 

model predictions presented here suggest that Columbia Glacier will most likely begin an irreversible phase 

of very rapid retreat in the 1980's. This irreversible 

retreat will probably first manifest itself by a continued retreat during the winter months when the terminus has, heretofore, readvanced (fig. 7).

REFERENCES CITED

Bindschadler, R. A., 1978, A time-dependent model of temperate 

glacier flow and its application to predict changes in the surgetype Variegated Glacier during its quiescent phase: Unpublished 

thesis, University of Washington Geophysics Program, Seattle, 

Washington, 244 p.

 1982, A numerical model of temperate glacier flow applied to 

the quiescent phase of a surge-type glacier: Journal of Glaciology, 

v. 28, no. 99, p. 239-265.

 1983, The importance of pressurized sub-glacial water in 

separation and sliding at the glacier bed: Journal of Glaciology. 

(In press.)

Brown, C. S., Meier, M. F., and Post, Austin, 1982, The calving relation of Alaskan tidewater glaciers, with application to Columbia 

Glacier: U.S. Geological Survey Professional Paper 1258-C, 13 p.

Budd., W. F., 1970, The longitudinal stress and strain-rate gradients 

in ice masses: Journal of Glaciology, v. 9, no. 55, p. 19-27.

Fountain, A. G., 1982, Columbia Glacier Photogrammetric Altitude 

and Velocity: Data Set (1957-1981), U.S. Geological Survey 

Open-File Report 82-756, 225 p.

Glen, J. W., 1955, The creep of polycrystalline ice: Proceedings of the 

Royal Society, Series A, v. 228, p. 519-538.

Iken, Almut, 1981, The effect of the sub-glacial water pressure on 

the sliding velocity of a glacier in an idealized numerical model: 

Journal of Glaciology, v. 27, no. 97, p. 407-421.

Kamb, Barclay, 1970, Sliding motion of glaciers Theory and observation: Reviews of Geophysics and Space Physics, v. 8, no. 4, p. 

673-728.

Kollmeyer, R. C., Motherway, D. L., Robe, R. Q., Platz, B. W., and 

Shah, A. M., 1977, A design feasibility study for the containment 

of icebergs within the waters of Columbia Bay, Alaska: U.S. 

Coast Guard Final Report No. 128, 135 p.

Lliboutry, Louis, 1968, General theory of sub-glacial cavitation and 

sliding of temperate glaciers: Journal of Glaciology, v. 7, no. 49, 

p. 21-58.

Mayo, L. R., Trabant, D. C., March, Rod., and Haeberli, W., 1979, 

Columbia Glacier stake location, mass balance, glacier surface 

altitude, and ice radar data, 1978 measurement year: U.S. Geological Survey Open-File Report 79-1168, 72 p.

McCracken, D. D. and Dorn, W. S., 1964, Numerical methods and 

FORTRAN programming: John Wiley & Sons, Inc., New York, 

455 p.

Meier, M. F., 1968, Calculations of slip of Nisqually Glacier on its 

bed: No simple relation of sliding velocity to shear stress: Union 

de Geodesie et Geophysique Internationale, Association Internationale d'Hydrologie Scienfitique, Assemblee generale de 

Berne, 1967, Commission de Neiges et Glaces, Rapports et discussions, p. 49-57.

Meier, M. F., Kamb, Barclay, Alien, C. R., Sharp, R. T., 1974, Flow of 

Blue Glacier, Olympic Mountains, Washington, U.S.A.: Journal 

of Glaciology, v. 13, no. 68, p. 187-212.

Meier, M. F., Post, Austin, Brown, C. S., Frank, David, Hodge, S. M., 

Mayo, L. R., Rasmussen, L. A., Senear, E. A., Sikonia, W. G., Trabant, D. C., and Watts, R. D., 1978, Columbia Glacier progress 

report December 1977: U.S. Geological Survey Open-File 

78-264. 56 p.

Meier, M. F., Rasmussen, L. A., Post, Austin, Brown, C. S., Sikonia, 

W. G., Bindschadler, R. A., Mayo, L. R., and Trabant, D. C., 1980, 

Predicted timing and disintegration of the lower reach of Columbia Glacier, Alaska: U.S. Geological Survey Open-File Report 

80-582, 47 p.

Morland, L. W., 1976a, Glacier flow down an inclined wavy bed: 

Journal of Glaciology, vo. 17, no. 77, p. 447-462.

 1976b,, Glacier flow down an inclined wavy bed with friction:

Journal of Glaciology, v. 17, no. 77, p. 463-477. 

Nye, J. F., 1953, The flow law of ice from measurements in glacier

tunnels, laboratory experiments and the Jungfraufirn borehole
REFERENCES CITED 17

experiment: Proceedings of the Royal Society, Series A, v. 219, 

p. 477-489.

 1957, The distribution of stress and velocity in glaciers and ice 

sheets: Proceedings of the Royal Society of London, Series A, 

v. 239, p. 113-133.

 1965, The flow of a glacier in a channel of rectangular, elliptic,

or parabolic cross-section: Journal of Glaciology, v. 5, no. 41,

p. 661-690. 

Post, Austin, 1975, Preliminary hydrography and historic terminal

changes of Columbia Glacier, Alaska: U.S. Geological Survey

Hydrologic Investigations Atlas 559, 3 sheets.

 1980a, Preliminary bathymetry of Northwestern Fiord and 

neoglacial changes of Northwestern Glacier, Alaska: U.S. Geological Survey Open-File Report 80-414, 2 sheets.

 1980b, Preliminary bathymetry of Blackstone Bay and neoglacial changes of Blackstone Glacier, Alaska: U.S. Geological Survey Open-File Report 80-418, 2 sheets.

 1980c, Preliminary bathymetry of Aialik Bay and neological 

changes of Aialik and Pederson Glaciers, Alaska: U.S. Geological 

Survey Open-File Report 80-423, 1 sheet.

 1980d, Preliminary bathymetry of McCarty Fiord and neoglacial changes of McCarty Glacier, Alaska: U.S. Geological Survey 

Open-File Report 80-424, 4 sheets.

 1980e, Preliminary bathymetry Esther Passage to Eaglek

Island, Alaska: U.S. Geological Survey Open-File Report 80-426, 

1 sheet.

 1980f, Preliminary bathymetry Approaches to Unakwik Inlet, 

Alaska: U.S. Geological Survey Open-File Report 80-428, 1 

sheet.

Rasmussen, L. A., and Meier, M. F., 1982, Continuity equation model 

of the drastic retreat of Columbia Glacier, Alaska: U.S. Geological Survey Professional Paper 1258-A, 23 p.

Sikonia, W. G., and Post, Austin, 1979, Colu;mbia glacier, AlaskaRecent ice loss and its relationship to seasonal terminal embayments, thinning, and glacier flow: U.S. Geological Survey OpenFile Report 79-1265, 2 sheets.

Sikonia, W. G., 1982, Finite element model of drastic retreat of 

Columbia Glacier, Alaska: U.S. Geological Survey Professional 

Paper 1258-B, 74 p.

 1983, Two-dimensional adjustment of point data, Columbia 

Glacier, Alaska: U.S. Geological Survey Open-File Report. (In 

press.)

Waddington, E. D., 1981, Accurate modelling of glacier flow: Unpublished thesis, The University of British Columbia, Vancouver, 

British Columbia, 460 p.

Weertman, Johannes, 1964, The theory of glacier sliding: Journal of 

Glaciology, vo. 5, no. 39, p. 287-303.