r j

« m? j T\ *

 »'' *- 

'iA

! OC

O i^

O Q

. X.
Continuity Equation Model of the Predicted 

Drastic Retreat of Columbia Glacier, Alaska

By L A. RASMUSSEN and M. F. MEIER

STUDIES OF COLUMBIA GLACIER, ALASKA

GEOLOGICAL SURVEY PROFESSIONAL PAPER

A discussion of a one-dimensional numerical 

model developed to predict the retreat and the 

iceberg discharge from a large iceberg-calving 

glacier

UNITED STATES GOVERNMENT PRINTING OFFICE, WASH INGTON: 1982
UNITED STATES DEPARTMENT OF THE INTERIOR 

JAMES G. WATT, Secretary

GEOLOGICAL SURVEY 

Dallas L. Peck, Director

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

604 South Pickett Street, Alexandria, VA 22304
CONTENTS

Symbols and abbreviations. 

Abstract _________. 

Introduction _______. 

Glacier geometry ____. 

Flux calculations _____. 

Calving mechanism ___. 

Columbia Glacier data __.

AIV 

1 

1 

3

3

4

6

Longitudinal profiles ______.

Calculated terminus retreat ___.

Influence of error in data or model.

Effect of individual errors _.

.____________ A7

.____________ 9

.____________ 9

.____________ 9

Effect of superimposed errors _______________ 15

Effect of numerical error __________________ 17

Discussion _______________________ _ 17

References cited _________________________ 23

ILLUSTRATIONS

	Page

FIGURE 1. Index map of Columbia Glacier, Alaska _____________________________________________ A2

2. Surface and bedrock maps of lower reach of Columbia Glacier, Alaska ___________________________ _ 5

3. Diagram showing geometry of a typical glacier cross section (at Xj), showing half of the symmetric structure _____ 6 

4-22. Graphs showing:

4. Longitudinal geometry of lower reach of Columbia Glacier, Alaska _____________________ _ 8

5. Position of the terminus as a function of time ____________________________________ _ 10

6. Prime set (case 1) of postulated longitudinal profiles ________________________________ __ 10

7. Sensitivity of the model results to the value of the calving coefficient c _________________________ 11

8. Features of the calculations for c=28.3 ___________________________________________ 12

9. Longitudinal profiles assumed by the glacier in the calculations using the case 1 data and c =28.3 ___________ 14

10. Volume balance functions, V(z), for cases 1,2a, 2b, 2c, and 2d _______________________________ 15

11. Centerline bed altitude profiles, U(x), for cases 1,4a, and 4b _______________________________ 16

12. Shoal bathymetry, U(x), for cases 1,5a, and 5b ___________________________________ __ 17

13. Longitudinal profiles for case 6a that delay the transition from continuing vigorous thinning to beginning to approach

fm(x) _______________________________________________________________________________________ 18

14. Longitudinal profiles for case 6b that hasten the transition from continuing vigorous thinning to beginning to approach

fm(x) _______________________________________________________________________________________ 18

15. Longitudinal profiles for case 7a that are compatible with an initial retreat rate ofXt = - 36 m/a ____________ 19

16. Longitudinal profiles for case 7b that are compatible with an initial retreat rate ofXt = - 54 m/a ____________ 19

17. Final longitudinal profile,/m(x), for cases 1,8a, 8b, and 8c _________________________________ 20

18. Effect on the solution of varying the bed topography, U(x) _________________________________________________ 20

19. Effect on the solution of varying the shoal bathymetry, U(x) ________________________________________________ 21

20. Effect on the solution of varying the initial retreat rate,Zi ________________________________ 21

21. Effect on the solution of varying the final longitudinal profile,/m(x) ___________________________ 22

22. Resultant summation of the individual errors that may be present in the data used to describe Columbia Glacier ___ 22

TABLES

TABLE 1. Effect on the calving calculation of varying one or another of the components of the data while holding all others at the values

used in the prime case ___________________________________________________ A13

2. Effect on the calving calculation of varying one component and simultaneously adjusting c ___________________ 14

3. Effect on the solution due to truncation error ___________________________________________ 17

III
SYMBOLS AND ABBREVIATIONS

Symbol Name

Units 

(where applicable)

A 

a

Of

b

6'

c

c

D

F

f 

G

9 

y 

h

Glen's flow law parameter.

Year.

Glacier surface slope

__ a"'bar""

Local ice balance rate, in ice equivalent _______________________________

Ice balance rate function, in ice equivalent _______________________________

Calving coefficient related to average water depth _______________ 

Calving coefficient related to maximum water depth _______________________

Valley cross-section power-law parameter.

Intermediate quantity (eqs. 21,22) _____________________________________

Longitudinal surface profile __________________________________________

Standard deviation of data error _______________________________________

Acceleration due to gravity ___________________________________________

Shape factor for scaling surface velocity to average velocity through thickness __ 

Glacier thickness, in the z direction _____________________________________

Water depth at centerline at terminus _________________________________

Average water depth across width at terminus ____________________________

Subscript denoting terminus position _________________________________

Subscript denoting x coordinate ________________________________________

Subscript denoting upglacier end of lower reach _.

Kilometer.

Measure of error in predicted retreat __________

Meter.

dimensionless

m/a

m/a a'1

a'1

km/a

km

a

km/a2

dimensionless

km

km

km

dimensionless

dimensionless

dimensionless

3

k

km

M Measure of error in predicted retreat ________________________________________ dimensionless

m

ra Subscript denoting retreated terminus position ___________________________________ dimensionless

n Glen's flow law exponent _______________________________________________ dimensionless

P Probability level of data error ____________________________________________ dimensionless

Q Glacier flux ______________________________________________________ kmVa

Q» Balance flux _____________________________________________________ km3/a

Qe Calving flux _____________________________________________________ km3/a

Q» Thinning flux ____________________________________________________ km3/a

R Volume loss above lower reach of glacier ______________________________________ km3

r Valley cross-section power-law exponent.

p Ice density.

S Area of glacier cross section _____________________________________________ km2

s Shape factor _____________________________________________________ dimensionless

t Time _________________________________________________________ a

U Centerline bed altitude _________________________________________________ km

V Volume of lower reach of glacier ___________________________________________ km3

v Glacier surface speed, horizontal component ____________________________________ km/a

vc Calving speed _____________________________________________________ km/a

W Glacier width ____________________________________________________ km

X Terminus position __________________________________________________ km

X Rate of advance of terminus _____________________________________________ km/a

XE Early estimated retreat _______________________________________________ km

XL Late estimated retreat ________________________________________________ km

x Horizontal coordinate, positive in direction of flow _________________________________ km

Z Glacier surface altitude _______________________________________________ km

z Vertical coordinate, positive upward _________________________________________ km

IV
STUDIES OF COLUMBIA GLACIER, ALASKA

CONTINUITY EQUATION MODEL OF THE PREDICTED DRASTIC 

RETREAT OF COLUMBIA GLACIER, ALASKA

By L. A. RASMUSSEN and M. F. MEIER

ABSTRACT

A one-dimensional numerical model based on the continuity equation 

is developed to predict the retreat rate of the terminus of Columbia 

Glacier, a large (1,100 km2) grounded iceberg-calving glacier, and to 

estimate the time distribution of the attendant iceberg discharge. An 

extensive field program, using September 1,1977, through August 31, 

1978, as the principal data year, obtained ice surface altitudes and 

velocities, mass-balance values and thickness changes, and bed 

topography estimated from radio-echo sounding.

For the terminus itself, the continuity equation is used to express the 

retreat rate as the difference between the glacier flux to the terminus 

and the iceberg flux from it. The continuity equation also is applied to 

the lower 14 km of the glacier, which constitutes about 5 percent of its 

1978 area, to determine the ice flux to the terminus as the sum of the 

ice flux into the lower reach of the glacier and the flux increment 

resulting from the thinning of the lower reach not due to ablation. The 

iceberg flux is assumed to be proportional to the average water depth 

at the terminus. The seasonal variations of calving and mass balance 

are averaged through the year so that all calculations are performed on 

annualized data.

Instead of independently calculating the glacier flow dynamics, the 

model requires that a sequence of longitudinal profiles be supplied. On 

the basis of the calving relation and the other applications of the continuity equation, the model then determines the times at which the 

glacier assumes each of the supplied profiles. The effect on the 

predicted retreat rate caused by the arbitrariness of the longitudinal 

profiles used to describe the lower reach of Columbia Glacier appears 

to be less than the effect of oftier uncertainties in the data used.

The predicted retreat slowly accelerates from the initial observed 

rate of 45 m/a until 1983.0 ±0.9, when the terminus reaches deep 

water, and the retreat rate rapidly increases to about 4 km/a. By some 

time between 1983.8 and 1986.5, the terminus is expected to have 

retreated 9.5 km. A maximum iceberg calving flux of about 10 kmVa, 

which is 6 to 8 times the initial amount, is expected to occur in about 

1983. It is assumed, in making this prediction, that the annual mass 

balance and other factors do not substantially differ from the 1977-78 

observation.

INTRODUCTION

Nearly all grounded, iceberg-calving glaciers experience large-scale asynchronous advances and retreats. 

This behavior apparently is not directly related to 

climatic variations. A critical factor is the water depth 

at the terminus, as instability results when the glacier 

retreats even a short distance into a deep fiord or basin.

The glacier then may retreat rapidly and irreversibly as 

unusual volumes of ice break away (Post, 1975). Observations of such glaciers reveal that they can advance 

only by moving a moraine shoal forward, so that the terminus rests in shallow water.

Columbia Glacier, near Valdez, Alaska (fig. 1), is a 

large calving glacier; it is 67 km long and 1,100 km2 in 

area. It ends on a moraine shoal in shallow water, but, 

upglacier from the terminus, the bed is about 400 m 

below sea level. Except for some small areas associated 

with ice-dammed lakes, the glacier is grounded throughout; none of it is floating at the terminus. Although the 

position of the terminus has been in a state of near 

equilibrium since the first recorded observation in 1794, 

evidence now suggests that rapid, drastic retreat may 

be imminent (Post, 1975). Small icebergs drift from Columbia Glacier toward and occasionally into Valdez Arm 

(Kollmeyer and others, 1977). Drastic retreat would 

substantially increase the discharge of ice and, thus, 

would increase the iceberg hazard to shipping. To determine when this might happen and how much the iceberg 

discharge might increase, the U.S. Geological Survey 

began an intensive study in 1977 (Meier and others, 

1978). Much of the field data has been reported (Mayo 

and others, 1979), and a preliminary prediction was 

issued in 1980 (Meier and others, 1980).

This report describes a newly developed time-dependent, one-dimensional numerical model, based mainly on 

an equation of continuity that can be used to calculate 

the terminus retreat and the rate of iceberg discharge 

for a calving glacier. This simple model was developed 

and used for the Columbia Glacier prediction because 

more complex models, which fully utilized glacier 

dynamics theory, could not be developed by the time the 

prediction had to be given. The simple model described 

here has been tested on one other calving glacier and 

should have general application to the study of the asynchronous variations and the rapid retreats of other calving glaciers. Because the model does not include indeAl
A2 STUDIES OF COLUMBIA GLACIER, ALASKA

Boo

Area shown 

in figure 2

 61°00'

10 KILOMETERS

FIGURE 1.-Index map of Columbia Glacier, Alaska. Arrows show direction of flow. Main ice stream is indicated by longer

arrows and dots at 2-km intervals along the longitudinal coordinate system.
CONTINUITY EQUATION MODEL, PREDICTED DRASTIC RETREAT, COLUMBIA GLACIER, ALASKA A3

pendent calculation of the glacier flow dynamics, it can 

be applied only to a glacier that, on an annual basis, is 

already retreating.

The report begins with a description of how the threedimensional shape of the lower Columbia Glacier is 

described in terms of one independent space dimension, 

by using a simple analytic approximation of the crosssectional shape (see section entitled "Glacier 

Geometry"). The ice flux (discharge), consisting of timeand distance-varying combinations of balance flux and 

thinning flux components, is examined in the section 

"Flux Calculations." A terminus boundary condition, a 

relation between calving speed and water depth, is 

discussed in "Calving Mechanism."

The one-dimensional numerical model requires that 

the field data, obtained over three-dimensional space, be 

averaged, smoothed, or otherwise transformed to form 

a one-dimensional data set; this is described in the section "Columbia Glacier Data." A crucial aspect of this 

model is the need for future longitudinal surface 

profiles. The section entitled "Longitudinal Profiles" explains how these profiles are constructed. The profiles 

are not initially identified with time; the model 

calculates when the glacier reaches each profile and, 

thus, predicts the time-variation of terminus retreat and 

the iceberg flux (see "Calculated Terminus Retreat").

The section entitled "Influence of Error in Data or 

Model" describes the influence of error in the data or in 

the assumptions of the model on the calculated terminus 

retreat and iceberg flux and includes the effect of 

numerical error. The results of this sensitivity analysis 

are summarized in a set of tables. Finally, the quality of 

the prediction and its sensitivity to data error and model 

assumption are examined and discussed in the last section (see "Discussion").

GLACIER GEOMETRY

The lower reach of the glacier is represented by a 

longitudinal coordinate system (x,z) in which x is the 

horizontal slightly curvilinear axis (fig. 2), positive in the 

direction of flow, and z is the vertical axis, positive upward. The terminus position is denoted X, and the 

altitude of the glacier surface, Z. Data are provided at a 

set of Xi that are equally spaced and numbered in the 

direction of negative x (during retreat this is the direction X changes with time), so that xt+i=xt -hx and

Xi is the position of the terminus X at time ti,

xm is the position of the terminus X at time tm, and

xk is the farthest point considered up the glacier.

When the terminus position X=xt , at time tt, the 

longitudinal surface profile is denoted f,(x), which is the 

piecewise linear function defined by the points (xhZt^ 

(xM,ZlM\ . .. , (xk,Zik\ That is, Xy =/,fo). In general,

where a variable has two subscripts, the second refers to 

the x location of the value, and the first refers to the x 

location of the terminus. Thus, where 6 is any variable, 

0y is the value of 0 at Xj when the terminus is at X,.

The valley cross section is assumed to follow a power 

law; that is, the width at altitude z is

W(z)=D(z-Uy (1)

where U is the bed altitude at the valley centerline and 

the power law parameters are D and r. The area of the 

cross section is obtained by integrating the width

= \W(z)dz=D(z- 1)= W(z)h(z)l(r+ 1) (2)

in which a horizontal transverse glacier surface profile is 

assumed and where

h(z)=z-U (3)

is the glacier thickness at the centerline (fig. 3). If, for 

some z, the width, area, and centerline depth are 

specified, then the power law parameters are obtained:

r=W(z)h(z)IS(z)-l

D = W(z)lh'(z). 

For a concave cross section, 0<r< 1, and thus

;i

l<W(z)h(z)IS(z)<2.

(4)

(5)

If the bed at the centerline is below sea level (C7<0), 

then the centerline water depth is ^,= -[7, and the 

average water depth across the width of the glacier is

T^ - S(0)/ W(Q) = hj(r+1) = - Ul(r + 1). (6)

Water depth is defined as the vertical distance below sea 

level to the bed; it is taken to be positive.

The geometry of the lower reach of the glacier is 

determined by specifying the initial longitudinal profile, 

fi(x), and, for 1 £>j<k, the valley centerline altitudes, UJt 

and the initial cross section widths, Wv, and areas, <SV. 

From these, the power law parameters Dj and r, are 

found by using equations 3 and 4. For any other longitudinal profile, f{x\ equations 1 and 2 are used for 

i<,j<,k to get the corresponding widths, WtJ, and areas, 

Sv.

FLUX CALCULATIONS

As the terminus recedes from X=Xi at time tt to 

X=xt+i at time tM , the cross-section area change at any 

is

 St+ -S (7)
A4 STUDIES OF COLUMBIA GLACIER, ALASKA

and the total volume change over the entire lower reach 

is calculated, by integrating trapezoidally,

 ^tf fo-a^

2 jmi

(8)

where (AS,), is defined to be zero. This volume change 

contributes a flux component at the terminus:

(9)

The volume -jk^- ,+i) is lost to calving (see

equation 18). When the terminus is at X=xh the balance 

over the lower reach also contributes a flux component 

at the terminus

^b'(Z^} (10)

where b'(Z) is a balance function, in ice equivalent.

The total flux at xk , the top of the lower reach, is 

assumed to obey Glen's flow law (Glen, 1954), and the 

flow is assumed to be laminar (longitudinal stresses are 

neglected). Also, both the deformational and the sliding 

components of velocity (as in Budd and others, 1979) are 

taken to follow the same form

2Ah 

n+1 (spgh sin a)" (ID

in which v is the horizontal component of the total velocity at the surface, A and n are the flow law parameters, s 

is a shape factor, p is the density of glacier ice, g is the 

acceleration due to gravity, and tan a = - dZ/dx is the 

surface slope. This equation is not strictly correct 

because h should be measured perpendicular to the 

glacier surface. However, the surface slope of the lower 

Columbia Glacier is rather small (0.002 < tan a < 0.005); 

as a result, whether h is measured vertically or perpendicular to the surface makes little difference (Budd and 

Jenssen, 1975). The total flux through the cross section 

of area <S is then

Q=yvS (12)

where 7 is a shape factor that scales down the surface 

velocity to give the average velocity through the section 

(Nye, 1952; Bindschadler, 1978). If, at the top of the 

lower reach, xk , equations 2 and 3 are used to express S(k 

and hik in terms of Zik , then equations 11 and 12 can be 

combined and divided to eliminate A, s, and y, yielding

= [(Zik - Uk)l(Zlk - Uk)Y"^ (sin a,.k/sin alk)» (13)

where sin «,k = tan <xj(l + tan2 a,k)"2 and tan a,fc *(Z«- 

Zlk_,)/A», the Zik being taken from the supplied longitudinal profiles/.(x). Dividing out A, s, and 7 is done on the 

assumption that the flow characteristics at xk do note

change during the time the terminus retreats from xl to 

xm .

The total flux to the terminus consists of three components: the flux into the lower reach, the flux contribution caused by thinning of the lower reach, and the flux 

contribution due to the ice balance at the surface of the 

lower reach. The average value, as the terminus recedes 

from xt to x,+1 , of the total flux to the terminus is obtained as the sum of the three components, using equations 9, 10, and 13:

(AQ.U 2 . (14)

The total flux, QIk , at the top of the lower reach, xk , can 

also be determined by applying the continuity equation 

to the upper reach of the glacier (0<»<ajk) in the form

Qik= j j(i_i). ji+ j -/* s(Q6),k +(Q*),k (15)

where b is the local ice balance and h is the local glacier 

thickness change, both averaged annually, and the domain of the integrals is the entire area of the glacier 

above xk . At time tt, the total flux at xk consists of a 

balance component (Qb)ik and a thinning component

(Q0«-

If the times, tt, when the terminus positions are at the 

Xi are known for l<i<m and if it is assumed that 

(Qb)ik=(Qb)ik is constant in time, then the total volume 

removed from the reservoir of ice above xk (R) may be

calculated

tm tm

R= \ (Qh)*dt= Qikdt-(tm -tl)(Qb)lk (16)

(17)

where the time integral is approximately

j

m-l

equation 13 being used to compute the QIk . The determination of R is not necessary to calculating the retreat 

rate but is useful in estimating the effect the retreat of 

the terminus has on the remaining glacier.

CALVING MECHANISM

As developed in Meier and others (1980), a calving 

equation may be obtained by applying the continuity 

equation at the glacier terminus in the form

SX=Q-QC , (18)

FIGURE 2. - Surface and bedrock maps of lower reach of Columbia 

Glacier, Alaska. Contour interval is 200 m on exposed 

bedrock, 100 m on glacier surface, and 100 m 

(dashed) on subglacial terrain. Subglacial terrain map 

was derived by 2-dimensional continuity adjustment 

(Sikonia and others, written commun., 1981) from radioecho sounding and other field measurements. For location, see figure 1.
CONTINUITY EQUATION MODEL, PREDICTED DRASTIC RETREAT, COLUMBIA GLACIER, ALASKA A5

147

EXPLANATION

Exposed land 

Water

1 KILOMETER

CONTOUR INTERVAL 10O METERS 

(-20 m contour is dotted )

NATIONAL GEODETIC VERTICAL DATUM OF 1929

Surface contours

 -. .__? Glacier bed and bathymetric contours 

 66 Longitudinal coordinates in km
A6 STUDIES OF COLUMBIA GLACIER, ALASKA

WIDTH, W

FIGURE 3.-Geometry of a typical glacier cross section perpendicular to the x axis (at Xj), showing half of the symmetric 

structure. When the terminus is at Xi>Xj, the surface 

altitude is Ztf, the width is WUt and the area is 

Sv=2(SI+SII+Sai). When the terminus is at x}, the surface altitude is ZJJt the width is WJJt and the area is 

SJj=2(SI+Sa). The sea-level width is W^O), and the area 

is SJ (0)=2SI. The altitude of the bed centerline is Uj. 

When Uj<0, the centerline water depth is hw = - Uj, and 

the average water depth is hw =SJ (G)IWJ (0).

where S is the area of the projection of the terminus onto a vertical plane normal to the x axis, X is the time rate 

of change dX/dt of the position of the terminus, Q is the 

volume flux of ice to the terminus, and Qc is the iceberg 

calving flux from the terminus, which is taken to have 

the form

Qc =Schw (19)

where c is a coefficient with dimensions cr1 and hw is the 

average water depth across the width of the terminus. If 

equation 19 is substituted into equation 18, and if both 

sides are divided by S, then

X=Q/S-chw, (20)

in which the ratio Q/S is defined to be the average 

glacier speed, v, at the terminus, and the ratio QC/S = chw 

is defined to be the calving speed, vc . In the case of 

retreat, the calving speed exceeds the average glacier 

speed, and X is negative.

When the terminus is at X= xi} equation 20 is approximated, in terms of the data spaced equally in x,

where (21)

2Q.- + V2 [

Ut UM

that, when solved for tM, gives the time when the terminus reaches the position X=XM:

i i A «-/Z7l /OO\ *i+i = tt AX/r ,>i/2 . (^^) 

Equation 22 may then be used to calculate the times, tit

The coefficient c gives only the annually averaged calving rate and cannot be used to obtain an instantaneous 

value. Although the geometry of the entire glacier 

changes only negligibly throughout the year because of 

relatively small seasonal variations of thickness and terminus position, the calving rate does have a pronounced 

seasonal variation with almost no calving in winter and 

high calving during summer. This variation, which 

causes a corresponding annual variation in terminus 

position, may be related to seasonal runoff and lake outbursts (Sikonia, in press; Sikonia and Post, 1980). As a 

consequence, equation 21 produces an annually averaged X(t) result.

An alternative form of the calving equation is obtained 

by relating the calving rate to the maximum water 

depth, hw, rather than to the average water depth, hw, as 

considered by C. S. Brown and others (written 

commun., 1980),

Qc =Schw. (23)

Through equation 6, using the hw relation modifies the 

calculation only by removing the r, terms from the 

denominators in equation 21.

COLUMBIA GLACIER DATA

The principal data year for the Columbia Glacier field 

program was September 1, 1977, through August 31, 

1978 (Mayo and others, 1979). Because the present 

model does not account for seasonal variations, annually 

averaged data are used. The middle of the data year is 

chosen at ^ = 1978.2, at which time, the transversely 

averaged position of the terminus was X=xl = 66.6 km, 

as measured from the head of the glacier. The surface 

topography also was averaged transversely to form the 

initial longitudinal surface profile, f^x), up to xk =52.6 

km, where tan «lk = 0.039. Above this point, the glacier 

consists of a broad basin receiving the flow of three 

large tributary branches, as shown in figures 1 and 2.

The x,y bed topography calculated by W. G. Sikonia 

and others (written commun., 1981) is used to obtain the 

valley centerline altitudes, C7, and, by using the 

measured glacier widths, W^, to determine the initial 

cross-section areas, <SV.
CONTINUITY EQUATION MODEL, PREDICTED DRASTIC RETREAT, COLUMBIA GLACIER, ALASKA A7

The retreated position of the terminus is taken to be 

#m = 57.1 km, with A;c=0.1 km; thus, m=96 and fc = 141. 

Without specifying the time, tm, the calculation of which 

is one of the chief results of the modeling, a final longitudinal surface profile, fm(x), is postulated. These basic 

variables are shown in figure 4. The graphs also show 

the sea-level width, W(G); and area, S(0); the average 

water depth, hw ; and the width, Wa, and area, Sy> , at x}

In the section entitled "Longitudinal Profiles," values 

are reported with up to four significant figures; these 

values are used because the model requires internally 

consistent data and because they enable the reader to 

check results. This is not to imply that the data are this 

accurate; the fluxes Qlfc and Qn for instance, are known 

only to an accuracy of about 7 to 10 percent (M. F. Meier 

and others, written commun., 1981).

The initial flux, Qlfc at xk, is estimated to be 1.354 km3/a 

by transversely integrating the product of glacier thickness and the downglacier component of the surface 

velocity. In equation 12, the value 7 = 0.86 is used, as 

determined by Sikonia and others, (written commun., 

1981). According to equation 15, the Qlfc is partitioned 

into two components, (Qb )lk = 0.873 km3 /a and 

(Qh)ik = 0.481 km3/a, by integrating the x,y ice balance to 

get (Q6)lfc and subtracting (Qb)lk from Qlk to obtain (QA)lfc .

The measured ice balance over the lower reach is well 

approximated by a linear b'(z) with a value of -10.68 m/a 

at sea level and -1.42 m/a at z= 600 m. When integrated 

over the initial surface, fi(x), in equation 10, it yields a 

flux component at the terminus of (AQ6)i = - 0.475 km3/a 

for the balance over the lower reach of the glacier. The 

total balance flux at the terminus, thus, is 0.398 km3/a at 

time t^ The transverse integration of thickness and 

velocity at the terminus, where 7 is taken to be unity to 

represent pure sliding motion, gives an initial flux at the 

terminus of Qn = 1.343 km3/a. Therefore, the total thinning flux at the terminus is 0.945 km3/a, and the thinning flux component for the lower reach is (AQ/,)i = 0.464 

km3/a.

Photogrammetric analysis of the lower reach of the 

glacier (Sikonia and Post, 1980) reveals the pronounced 

seasonal variation of calving. Figure 5 represents the 

transversely averaged terminus position on 26 dates 

from 1976.6 through 1981.7 (M. F. Meier and others, 

written commun., 1981). The terminus position exhibits 

a roughly sinusoidal seasonal variation of about 300-m 

amplitude, superimposed on a 1977-81 average retreat 

rate of about X = - 45 m/a.

LONGITUDINAL PROFILES

Because the present model does not independently 

calculate glacier flow, longitudinal surface profiles, fi(x), 

/2(x), . . . , fm(x), must be supplied. When the terminus is

at X;, the longitudinal profile isf,{x). Through the calving 

relation (equation 22), the model then determines the 

times, tt , that the glacier assumes each of these m 

profiles.

The first two profiles are obtained from the field data, 

ft(x\ directly, and f2(x), by transversely averaging the 

observed thickness change rate and then scaling the 

averaged thickness change rate by the initial terminus 

retreat rate, Xt = - 45 m/a,

(24)

X,

so that, through equations 7 to 9,/1 (;c) and/2(«) produce 

the observed (AQA)j value. The effect of this extrapolation when A#>45 m is examined bf varying Asc from 50 

m to 500 m.

A final longitudinal profile, fm(x), is postulated by setting the flow law exponent, n=3, and arbitrarily assuming that Qmk = 2Qlfc . The effect of these arbitrary values is 

examined in the section entitled "Influence of Error in 

Data or Model." At #fc = 52.6 km, fm(x) has altitude 

Zmk = 430 m, compared with the initial altitude, Zlk = 452 

m, and, at xm = 57.1 km, it has cliff top altitude Zmm = 100 

m, compared with the initial cliff top altitude, Zn = 86 m. 

When Columbia Glacier has retreated to xm = 57.1 km, its 

surface area will be about 95 percent of its present 

value.

This fm(x) profile (figs. 4 and 6) is speculative but is 

glaciologically reasonable because it displays a high ice 

cliff and slope steepening near the terminus that is 

characteristic of rapidly retreating calving glaciers. 

Furthermore, the profile obeys equations 11 to 13 with 

n= 3 and has nearly constant flux from xk to xm. This/m(aj) 

profile has a fairly high slope and produces an increased 

discharge with a modest thickness change, as compared 

with the present profile. During the period 1957-78, h at 

xk was less negative compared with values farther downglacier, and, above xk , h changed progressively to almost 

zero at #=36 km, 16 km above xk (Sikonia and Post, 

1980). Also, Q at xk increased during the period 1977-80 

(Meier and others, written commun.); this increase is 

consistent with the concurrent increase in slope.

Obviously, the postulated fm(x) would not occur upon 

extrapolation of /,(#) and/2(«). The combination of pronounced thinning and slow retreat is believed to be a 

symptom of the final phase of the connection of the terminus to the shoal. This connection results in a lessening 

of surface slope over much of the lower reach, where the 

surface velocity has been observed to be slightly increasing. After a small amount of additional retreat behind 

the shoal, the change in the longitudinal profile is expected to undergo a transition; the profile changes will 

cease to produce vigorous thinning and will begin to approach fm(x). Experience with a dynamic model used
A8 STUDIES OF COLUMBIA GLACIER, ALASKA

400

CO

£ 200

UJ

UJ

O

-200

-400

flC 

O

CO

do

2 2

54 56 58 60 62 64 

DISTANCE, x, IN KILOMETERS

66
CONTINUITY EQUATION MODEL, PREDICTED DRASTIC RETREAT, COLUMBIA GLACIER, ALASKA A9

with Columbia Glacier data also suggests the existence 

of this transition (R. A. Bindschadler, written commun., 

1981). Figure 6 gives the set of longitudinal profiles interpolated between f2(x) and fm(x); shown are the ft{x) 

profile for i=l, 2, 3, 4, 5, 6, 11, 16, 21, 26, . . . , 91, 

96=m. The interpolation was made so that the fJ(x) 

profiles changed smoothly from f2(x) to fm(x), embodied 

the transition, and were free of waves; that is, at any x, 

the surface altitude was always decreasing with time.

The influence on the solution of the arbitrariness of 

this interpolation is examined by considering alternative 

interpolations; alternative fm(x) profiles also are considered. In fact, considerable uncertainty is involved in 

all the values used to describe the Columbia Glacier; the 

effect of such errors is the subject of the section 

"Influence of Error in Data or Model" of this report.

CALCULATED TERMINUS RETREAT

The section entitled "Longitudinal Profiles" describes 

the data adopted as the prime case (case 1). Setting 

Xi = - 45 m/a not only leads, through the effect of equation 24 on/2(a;), to the formation of a set of longitudinal 

profiles,/X&), but also provides a condition for determining the value of the calving coefficient, c. If X= -45 m/a 

is taken.to. apply not only at tlt but also until X= x2, then

t, = tt - AZ/XX = 1978.2 +100/45 = 1980.4 (25)

is the required condition, which is met when c=28.3 per 

year. This value is close to the average value of 27 ±2 

per year derived from observations of many Alaska 

glaciers (C. S. Brown and others, written commun., 

1981).

The sensitivity of the model to the value of c is shown 

in figure 7, which displays X(t) for three different values 

of c used with the case 1 data and displays the corresponding three Qe(t) curves.

The results of the calculations for c= 28.3 a'1 using the 

case 1 data are further illustrated in figures 8 and 9. 

Shown in figure 8 are the flux components at the terminus from above xk = 52.6 km and from the lower reach, 

as well as the retreat rate and calving and average 

glacier speeds. A maximum retreat rate of about 

J£=-4.1 km/a occurs in early 1983, and the total 

volume, R, removed from the reservoir of ice above xk is

FIGURE 4.-Longitudinal geometry of lower reach of Columbia 

Glacier, Alaska. At time tt = 1978.2, when the terminus is 

at a;, = 66.6 km, the surface profile is/j(x), the width is 

WM, and the cross-section area is Sy. At time, tm, when 

the terminus is at xm =57.1 km, the surface profile is 

fm(x). At sea level, the width is W(Q), and the area is S(0). 

The centerline bed altitude is U, and the average water 

depth is hw =S(Q)IW(Q). The Wa curve and SM curve give, 

at each xjt the width and cross-section area when the terminus is at Xj.

6.0 km3. A maximum calving flux of Qc =9.9 knvVa occurs at £=1983.5. Figure 9 shows the longitudinal 

profiles, ffa), assumed by the glacier at selected times 

during the retreat from id = 66.6 km to jcm = 57.1 km.

INFLUENCE OF ERROR IN DATA OR MODEL

EFFECT OF INDIVIDUAL ERRORS

To investigate the effect of possible errors in the data 

or assumptions used in the model on the calculated X(t) 

and Qc(t) results, several cases were created by varying 

one or another of the components of the data. Each 

case, except for the adjustment stated for it, consists of 

the same data used for case 1. The size of each data adjustment is approximately equal to the expected observational error (W. C. Sikonia and others, written commun., 1981).

Case 2 includes four subcases for different variations 

of the volume balance function, b'(z). Case 2a uses a 

1-m/a more positive V(z), along with an initial flux, Qlk, 

decreased by 0.072 knvVa to preserve the initial flux, Qu, 

at the terminus. Case 2b is the opposite, with b'(z), made 

more negative by 1 m/a and.Qlk increased by 0.072 

km3/a. Case 2c uses a steeper V(z), with -11.41 m/a at 

sea level and -0.51 m/a at z =600 m, which preserves 

the initial flux component, (&Qb)1 = -0.475 knWa, contributed by the balance over the lower reach. Case 2d 

uses a shallower U(z), with - 9.95 m/a at sea level and 

-2.32 m/a at z=600 m, which also preserves (AQ6)t . 

Figure 10 shows the original b'(z) function along with 

those for cases 2a, 2b, 2c, and 2d. Steeper is defined as a 

rapid increase in balance with increasing z; by 

shallower, a slower increase.

Case 3a uses a Qlk increased by 8 percent, and case 3b 

uses a Qlk decreased by 8 percent. The variation is applied to the initial value Qlk =1.354 km3/a and, through 

equation 13, proportionately scales all the Qlk values. 

The variation also additively modifies all the Qv by 

±0.08Qlk.

Case 4a uses a shallower bed topography, with Uj increased by 40 m for x ^65.6 km. Case 4b uses a deeper 

bed topography, with Uj decreased by 40 m for x^ 65.6 

km. Each of these variations is accompanied by a linear 

function of x (x> 65.6 km) that reduces the variation to 

zero at Xi - 66.6 km. These two cases are shown as figure 

11.

Case 5a uses a shallower shoal, raised by adding a 

linear function of x that decreases from 10 m at xv = 66.6 

km to zero at #=65.6 km. Case 5b uses a deeper shoal, 

lowered by the opposite application of the same linear 

function. These two cases are shown as figure 12.

Cases 2 to 5 all use the same set of longitudinal 

profiles, fjx), as case 1.
A10 STUDIES OF COLUMBIA GLACIER, ALASKA

66.8

cc

UJ 

UJ 66.7

z 66.6

Z 66.5

O

P

CO

O 66.4

Q.

(0

? 66.3

tr

UJ

66.2

1976.0 1978.0 1980.0 1982.0

FIGURE 5. - Position of the terminus as a function of time. Position X(t) is averaged over a 3.5-km width of active ice. The heavy line corresponds 

to.X"= -45 m/a. The measurement error in X at each data point is approximately equal to the diameter of the dot (15 m).

500

400

CO

oc

UJ 

UJ

5 300

N

n

UJ

Q

P200

100 -

VERTICAL 

EXAGGERATION 

20 :1

52 54 56 58 60 62 64

DISTANCE, x, IN KILOMETERS

FIGURE 6.-Prime set (case 1) of postulated longitudinal profiles ffa) for i=l, 2, 3, 4, 5, 6,11, 16,..., 91, 96.

66
CONTINUITY EQUATION MODEL, PREDICTED DRASTIC RETREAT, COLUMBIA GLACIER, ALASKA All

Z

vt

X

68

**66

64

O 62

H 

(O 

O

GL 

(/> 60

OC 

UJ 58

56

12

10

CC 

UJ 

Q_

(/> 

CC 

UJ

O 

O 

O

X

3

1978.0 1980.0 1982.0 1984.0 1986.0

0 O 

1988.0

FIGURE 7. - Sensitivity of the model results to the value of the calving coefficient c using the case 1 data. The terminus retreat, X(t), and the calving flux, Qc(t), are shown for c=25.5, 28.3, and 31.2.

Case 6a uses a set of longitudinal profiles, fj(x), that 

delays the transition from continuing vigorous thinning 

to beginning to approach fm(x). Case 6b uses a set that 

hastens the transition. These two sets of longitudinal 

profiles are shown as figures 13 and 14, which show the 

ft(x) for the same i values as in figure 6.

Case 7a uses a 20 percent smaller initial retreat rate, 

Xi = - 36 m/a, which produces t2 - 1981.0. Case 7b uses a 

20 percent larger initial retreat rate, Xv = - 54 m/a, 

which produces £2 = 1980.0. These varied initial retreat 

rates, through the effect of equation 24 on/2(#), require 

the formation of modified sets of longitudinal profiles. 

The profiles for cases 7a and 7b are shown in figures 15 

and 16, respectively; the /,(#) shown are for the same i 

values as in figure 6.

Case 8 includes three subcases for variation of the 

longitudinal profile, fm(x\ for the retreated terminus 

position. They were constructed, as were thefm(x) used 

for case 1, by using equations 11 and 12 and the assumed 

bed topography to produce nearly constant flux between 

#k =52.6 km and #m = 57.1 km. Case 8a uses an/m(ff) that 

assumes Qmk = 3Qlk instead of 2QU. Case 8b uses an/m(ff) 

that assumes the flow law exponent n=4 instead of 3. 

Case 8c uses an fm(x) that assumes Zmk = 420 m rather 

than 430 m. These cases are shown, along with the basic

fm(x), as figure 17. They each lead to a different set of 

longitudinal profiles,/,(&), that, except for i»m, closely 

resembles the set for case 1.

The results for these cases are summarized in table 1, 

which gives several significant points on the X(t) and 

Qc(£) curves for case 1 and also gives the differences 

from those values for cases 2 to 8. Because, from equation 19, the calving flux depends only on the calving 

coefficient, c, and on the cross-section geometry, Sh*, 

the dependence of Qc on x will be the same for all cases in 

which only data components other than c and Sh* are 

varied. Such variations may change X(t) and, therefore, 

Qc(0> but Qc(x) is unchanged. Also included isR, the total 

volume of ice removed from above xk = 52.6 km, as given 

by equation 16 and 17.

As is shown in table 1, several of the variation cases 

significantly affect £2 when the same value of the calving 

coefficient is used as in case 1. Because the condition 

£2 = 1980.4 is considered to be more reliable than the 

value c=28.3, the coefficient should be adjusted in the 

variation cases so that the condition is still met. Table 2 

shows the revised results for those cases that, for 

c=28.3, give a £2 significantly different from 1980.4; 

they are denoted by a prime.
A12 STUDIES OF COLUMBIA GLACIER, ALASKA

LIJ

> 4

OC 

LiJ 

Q.

£ 3 

LiJ

|

* 2

O

ffi

D 

O

X

u_

LU

2 

D

O

-1

oc

< in

OC 4 LiJ **

OC 

LiJ

5 o

D LiJ LiJ

Q.

-4

1978.0 1980.0 1982.0 1984.0
CONTINUITY EQUATION MODEL, PREDICTED DRASTIC RETREAT, COLUMBIA GLACIER, ALASKA A13

TABLE 1. -Effect on the calving calculation of varying one or another of the components of the data while holding alt others at the values used in 

the prime case (1). Except for case 9, the same calvingformula, Qc =2S.5SEW is used. Shown for each case are the differences, in years, from the 

results for case 1. The values in parentheses for case 1 are the elapsed time, in years, after 1978.2

Case Description or component varied

1 ____Prime data set, Qc = 28.3 Sh* _________ 

Component varied 

2a __ + ! m/afrXz) ... .._ . ,,....

2b_ -Im/a&Xz) _ _

2c steeper #(z) _

2d shallower V(z)__

3a _ 8% larger Q* _ _____

3b 8% smaller Q* __ _ __

4a +40 m bed (shallower) __ _ _

4b_ _- 40m bed (deeper) _ _ _

5a + 10 m shoal (shallower) _

5b _ _ - 10 m shoal (deeper) _

6a delayed/X*) transition _ _

6b hastened/^) transition _

7a __A>-36m/a _ _ __

7b *,- 54m/R

8a /-(o^withQ^-SQu

8b /_(*)withw=4 _ _ _ _

8c _/_(o;)withZ_»=420m_

9 _ -Prime, but Qc = 19.7 Sh* _

Time at which terminus 

position x (km) reaches:

66.5

1980.4 

(2.2)

0 

0 

0 

0

+0.7 

-0.4

+0.2 

-0.2

+ 1.1 

-0.6

0 

0

+0.6 

-0.4

+0.2 

0 

-0.1

0

65.6

1982.4 

(4.2)

0 

0 

0 

0

+0.9 

0.6

+0.6 

0.5

+ 1.5 

-0.8

+0.3 

-0.3

+0.9 

-0.7

+0.5 

0 

0.2

fi 9

57.1

1984.9 

(6.7)

-0.1 

+0.1 

0 

0

+ 1.0 

0.7

+ 1.9 

-0.9

+ 1.5 

-0.8

+0.3 

-0.3

+0.8 

n R

+ 1.0 

0 

0.2

0.6

Time at which calving flux 

Q« (knWa) first reaches:

3

1982.0 

(3.8)

0 

0 

0 

0

+ 0.9 

0.6

+0.8 

0.5

+ 1.5 

-0.9

+0.2 

-0.3

+0.9 

-0.7

+0.4 

0 

0 2

n 9

6

1982.6 

(4.4)

0 

0 

0 

0

+0.9 

-0.6

+0.9 

-0.7

+ 1.5 

-0.8

+0.3 

-0.3

+0.9 

-0.7

+0.5 

0 

-0.2

0.4

9

1983.0 

(4.8)

0 

0 

0 

0

+0.9 

-0.6

-0.8

+ 1.5 

-0.8

+0.3 

-0.3

+0.8 

-0.6

+0.5 

0 

-0.2

0.5

Maximum

1983.5 

(5.3)

0 

0 

0 

0

+0.9 

-0.6

+0.9 

-0.6

+ 1.5 

-0.8

+0.3 

-0.3

+0.8 

-0.6

+0.6 

0 

-0.2

-0.4

Maximum

Q.

(km'/a)

9.9

0 

0 

0 

0

0 

0

-1.8

+2.0

0 

0

0 

0

0 

0

0 

0 

0

+ 1.6

Volume 

loss R 

km3)

6.0

-0.8 

+0.7 

0 

0

+ 1.7 

-1.4

+ 1.5 

-0.9

+0.8 

-0.5

+0.2 

-0.1

+0.5 

-0.3

+4.3 

0 

-0.5 

ft R

Case 2 has negligible effect except that cases 2a' and 

2b', through their adjustments of Qrt, make a change in 

R of about ± 12 percent.

Case 3, if c is not adjusted to maintain k = 1980.4, 

results in a large initial change to X(t) when X is weak 

and, thereafter, a slight progressive growth of that 

change. In case 3a, for example, when the Qlk are increased by 8 percent, more flux is produced at the terminus; through equation 18, this makes X less negative 

and slows retreat. If, for either case 3a' or case 3b', c is 

adjusted to make k = 1980.4, then, for t>t2, the variation made to the Q,fc is slightly over balanced, as shown 

in tables 1 and 2.

Case 4, if c is not adjusted to maintain £2 = 1980.4, 

results in a progressive change toX(t). In case 4a, for example, when the bed U(x) is made shallower (see fig. 11),

FIGURE 8.-Features of the calculations for c=28.3, using the case 1 

data. Shown are the balance flux, (Qt)», and thinning 

flux, (Qh)k, components of the total flux, Qt, into the 

lower reach from above a;* = 52.6 km; the flux components due to the balance, AQfc , and the thinning flux, 

AQ», over the lower reach; and the total flux at the terminus, Q,. Also shown are the average glacier speed at 

the terminus, v; the calving speed, vc ; and their difference (the retreat rate), A".

the Qc term is reduced through equation 19, which 

makes X less negative and slows retreat. If, for either 

case 4a' or case 4b', c is adjusted to make k = 1980.4, 

then, for t> t2 , the effect is diminished to about half. The 

adjustment to c that compensates for the smaller change 

in U(x) for x2 ^x^x1 is overmatched by the larger 

change to U(x) for x<x2 . The X(t) and Qc(£) curves for 

these cases and the curves for case 1 are shown in figure 

18.

Case 5, if c is not adjusted to maintain k=1980.4, 

results in a large change to X(t) over x»^x^xit where 

the variation of U(x) occurs and, thereafter, a slight progressive growth of that change. In case 5a, for example, 

when the shoal is raised (see fig. 12), the Qc term is 

reduced through equation 19, which makes X less 

negative and slows retreat. If, for either case 5a' or case 

5b', c is adjusted to make £2 = 1980.4, then, for x<x2 , the 

adjustment operates on the original U(x) and produces a 

result of the opposite sign. That is, if c is increased to 

compensate for a shallower shoal and make the terminus 

reach x2 at k=1980.4, then the larger c causes the terminus to traverse the remaining distance, where U(x) 

has not been varied, at a greater rate than in case 1. The 

X(t) and Qe(t) curves for these cases and the curves for 

case 1 are shown in figure 19.
A14 STUDIES OF COLUMBIA GLACIER, ALASKA

500 i i i i i i i i i r i r

400

cc

LU 

LU

N

n

LU 

Q

300

200

100

VERTICAL

EXAGGERATION 

20 M

l___i___i I____I I____I

52 54 56 58 60 62 

DISTANCE, x, IN KILOMETERS

64 66

FIGURE 9.-Longitudinal profiles assumed by the glacier in the calculations using the case 1 data and c=28.3. The /.(a;) are shown for £= 1978.2,

1980.4, 1981.3, 1982.2, 1983.2, 1984.2, and 1984.9.

TABLE 2.-Effect on the calving calculation of varying one component of the data and simultaneously adjusting c, as indicated, to constrain the 

solution to pass through x=66.5 km at ^=19804. Shown for each case are the differences, in years, from the result for case 1

Case

3a'

3b'

4a'

4b'

5a'

5b' _

6a'

6b' _

c

30.6

26.1

29.0

27.6

32.1

25.2

. _ 29.0

28.0

Time at which terminus 

position x (km) reaches:

66.5

0 

0

0 

0

0 

0

0 

0

65.6

0.1

+0.1

+ 0.3 

n 9

0.2

+0.2

+0.2 

0.1

57.1

0.3

+0.3

+0.9 

-0.6

-0.6

+ 0.7

+ 0.6 

0.1

Time at which calving flux 

Q« (km3/a) first reaches:

3

-0.1 

+ 0.1

+ 0.4 

0.3

0.2

+ 0.2

+ 0.1 

-0.1

6

-0.2 

+0.2

+ 0.5 

0.4

0.3

+ 0.3

+ 0.2 

0.1

9

-0.2 

+ 0.3

-0.5

-0.4

+ 0.2 

-0.1

Maximum

-0.2

+ 0.2

+ 0.5 

0.4

-0.3 

+ 0.4

+ 0.2 

-0.1

Maximum

Qc

(kmVa)

+ 0.8 

0.8

-1.6 

+ 1.7

+ 1.3 

-1.1

+ 0.2 

-0.1

Volume 

loss-R 

(km3)

+ 0.6 

-0.7

+ 1.1 

0.7

-0.7 

+0.8

+3.8 

-0.4
CONTINUITY EQUATION MODEL, PREDICTED DRASTIC RETREAT, COLUMBIA GLACIER, ALASKA A15

600

-10 -5 0 

ICE BALANCE, b', IN METERS PER YEAR

FIGURE 10.-Volume balance functions, b'(z), for case 1, 2a, 2b, 2c, and

2d.

Case 6 does not affect t2 because the first two longitudinal profiles, fi(x) and f2(x), are the same as those in 

the jX*) used with case 1. The main effect of the variation 

of the/,(#), as shown in figures 13 and 14, is on the relation of volume change to terminus position, dV/dX. 

Through equation 9, this controls the flux component, 

AQfc , at the terminus contributed by the thinning of the 

lower reach, but the total volume lost to thinning is unchanged by this variation. Therefore, the retreat rate is 

decreased and then increased (case 6a), or it is increased 

and then decreased (case 6b); in either case, only a small 

change to X(t) results.

Case 7 varies the initial retreat rate, jflf and, thus, by 

definition, presumes a proportionately different t2. As in 

case 6, the variation produces a retreat rate that slightly 

decreases and then slightly increases (case 7a) or that 

slightly increases and then slightly decreases (case 7b). 

The resulting X(t) and Qc(t) curves and the curves for 

case 1 are shown in figure 20.

Case 8 demonstrates that little effect is caused by the 

variations to/,(se) for i**m. Case 8a, though, does show a 

pronounced slowing of the retreat rate. Unlike case 3a, 

which assumes a Qik that ranges from 1.08 to 2.16 times 

the Qlfc used with case 1, case 8a assumes a range from

1.00 to 3.00 times that Qlfc value. Although case 8a attains a much higher Q,fc than does case 3a, the increase 

occurs when i is large (when the retreat rate is great) 

and, thus, has only a small effect on X and X(t)\ this is 

the condition F>0 in equations 21 and 22. Because case 

8a has Qik only slightly larger than in case 1 when i is 

small, it produces a progressively slowing retreat rate 

when compared with case 1. When c is adjusted to maintain t2 = 1980.4, the effect is reduced to about half. The 

X(t) and Qc(t) curves for the c-adjusted case 8a', for case 

1, and for the c-adjusted case 8c' are shown in figure 21.

Case 9 uses the alternative form of the calvh-t; relation 

given by equation 23. If the coefficient is chosen to have 

£2 = 1980.4 when the glacier data for case 1 are used, 

then c= 19.7; the X(t) curve produced is not identical to 

that of case 1 because Jijh^, is not constant with x. That 

ratio varies according to the variation with x of the 

valley cross-section shape, as given by equations 4 to 6.

The volume loss, R, above sefc = 52.6 km differs from 

case to case mainly because of the time, tm -tlt required 

for the terminus to retreat to sem = 57.1 km. A second 

order effect on the response of R depends on whether a 

particular case differs more from case 1 during the early 

phase of retreat, when Q,fc is small, or during the late 

phase of retreat, when Qik is large. This difference in 

response exists because Qik is a function of longitudinal 

profile, f{(x), and, thus, of terminus position, Xt, rather 

than of time.

EFFECT OF SUPERIMPOSED ERRORS

In the first phase of the sensitivity analysis, each 

variation was considered separately to examine the 

nature of its effect singly on the X(t) and Qc(t) curves. A 

second phase of the analysis takes into account the 

possibility that several of these errors may actually exist 

together in the data used here and that the curves 

should be accompanied by an estimate of the resultant 

summation of the effects of the individual errors. Each 

of the cases in this second phase of analyses is constructed on the presumption that the magnitude of the 

variation used is representative of the likely error in the 

quantity varied, as based on field program experience 

and glaciological judgment. The term "likely" is used in 

place of normal statistical measures of error, such as 

standard or probable, because many kinds of errors due 

to measurement or assumption cannot be defined quantitatively on the basis of known, or knowable, statistical 

distributions.

A resultant error is estimated for X(t) by assuming 

that there exists some particular, but as yet unspecified, 

probability, p, to which the term "likely error" pertains. 

When one case at a time is considered, the response of 

X(t) is caused entirely by the variation imposed on the
A16 STUDIES OF COLUMBIA GLACIER, ALASKA

-100

(0 

cc

LU 

LU

N

LU

0

-200

-300

-400

VERTICAL

EXAGGERATION 

20 :1

52 54 56 58 60

DISTANCE, x, IN METERS

62 64 66

FIGURE 11.-Center-line bed altitude profiles, U(x), for cases 1, 4a, and 4b.

defining data; thus, the same probability can be imputed 

to the resulting X(t) curve as is assumed for the data 

variation that produced it. Because the model constitutes a mildly nonlinear, but continuous, transformation of the data (that is, as shown in tables 1 and 2, 

X(t) does not respond exactly linearly to variations in the 

defining data), the two distributions are actually slightly 

different. If, however, a normal distribution is assumed 

to be reasonably valid, then p can be expressed in terms 

of M(p) standard deviations. Provided that the other 

cases are c adjusted to maintain t2 = 1980.4, so that they 

share none of the Xl variation that is the subject of case 

7, the cases are assumed to be statistically independent 

of each other.

Considering, first, only those subcases (7) that produce 

faster retreat than case 1, let

t^-tj^MGjt, (26)

where tjt denotes the time when X- Xj for some subcase / 

and where Gfl is the standard deviation of the tn values

yi -* )>]">.

about tn for that subcase. Since the variance of the sum 

of uncorrelated variables is the sum of the variances of 

those variables, the time, tJR , for the resultant summation of the effect of the individual errors is given by

(27)

Then, using the tJR for l^j^m, the curve X$) is constructed as an estimation of likely early retreat. Without 

having explicitly stated p, or the corresponding M, the 

curve Xe(t) is at the same probability level as denoted by 

the term "likely error." The estimated likely late retreat 

curve, XL(t), is constructed similarly by summing over 

only those subcases that produce slower retreat than 

case 1. The curves X^f) andXL(t) are shown, and the J£(£) 

for case 1, in figure 22. This approach to estimatingX$) 

and XL(t) is undertaken as a simpler alternative to mak-
-100

-200

-300

-400

CONTINUITY EQUATION MODEL, PREDICTED DRASTIC RETREAT, COLUMBIA GLACIER, ALASKA A17

TABLE 3. -Effect on the solution due to truncation error

VERTICAL 

EXAGGERATION 

20 :1

65 66

DISTANCE, x, IN KILOMETERS

67

FIGURE 12.-Shoal bathymetry, U(x), for cases 1, 5a, and 5b.

ing an extensive stochastic simulation with data not well 

enough known to justify such treatment.

EFFECT OF NUMERICAL ERROR

The numerical treatment of the data for Ao?=100 m introduces a truncation error that is smaller than the 

resolution at which the values are reported in tables 1 

and 2. Those values are given to tenths, and the truncation error is a fraction of a tenth. The truncation error 

value was determined by making the calculations for 

case 1 with defining data specified at equally spaced x{ 

for Ao?=500 m, 250 m, 100 m, and 50 m. Several significant features of the solution are shown in table 3 for 

each of the four Aa? values.

Data spacing As

50m 100m 250m 500m

Time terminus reaches: 

z-65.6 km

z-57.1 km

Time calving flux first reaches: 

Q«-3kmVa

Q,-6kmVa T ,. _ ,

Q«-9km3/a

Total discharge (km1) into lower reach 

through z= 52.6 km

Volume loss, R (km1), above z=52.6 km

1982.40

1984.96

1982.04

1982.65

1983.04

11.90

6.00

1982.36

1984.93

1982.01

1982.62

1983.00

11.87

5.99

1981.90

1984.46

1981.55

1982.15

1982.54

11.26

5.80

1981.22

1983.79

1980.85

1981.45

1981.85

10.86

5.48

DISCUSSION

This model uses the equation of continuity to predict 

the retreat rate of the Columbia Glacier and largely 

avoids uncertainties involved in glacier dynamics such 

as calculating the magnitude of sliding velocity at the 

glacier bed or the flow law of ice deforming in a highly 

transient mode. However, the model avoids these uncertainties at the cost of needing to specify a somewhat arbitrary set of longitudinal profiles. The sensitivity 

analysis shows that the different geometrical assumptions involved in constructing these profiles have some 

effect on the timing of retreat but little effect on the 

form of the retreat curve or the calving flux curve. 

These effects are, in fact, less than the effects of some 

uncertainties in the data.

The model is based on information collected during the 

1977-78 measurement year, except for the value of terminus retreat, X, which is the average over 4.5 years of 

record. The observed change in the glacier since 1978 

generally is consistent with the change expected. For instance, the ice discharge, Qk, at the upglacier limit of the 

model was observed to increase by 7.8 percent from the 

1977-78 measurement year to the 1979-80 measurement year (Meier and others, written commun., 1981); 

the expected increase (see fig. 8) was 7.4 percent.

Calving at the terminus, however, was much less in 

1980 than expected. Calving appears to be related to 

meltwater runoff (Sikonia and Post, 1980; Sikonia, in 

press). In the summer of 1980, the coldness and heavy 

snowfalls presumably decreased summer runoff. Furthermore, the outlet stream moved east, out of the main 

channel, and, as a result, the embayment was not extended upglacier. During the winter of 1980-81, the 

seasonal advance appeared to be ahead of expectation. 

However, rapid embayment formation during the summer of 1981 caused unusually rapid terminus retreat,
A18 STUDIES OF COLUMBIA GLACIER, ALASKA

500

400

CO

cc lil

300

z 

isT

PI

LU

§ 200

100

52

VERTICAL 

EXAGGERATION 

20:1

54 56 58 60 62 64 66

DISTANCE, x, IN KILOMETERS

FIGURE 13.-Longitudinal profiles for case 6a that delay the transition from continuing vigorous thinning to beginning to approach/^). The/fa) are shown for t-1, 2, 3, 4, 5, 6,11,16,..., 91, 96.

500

400

CO

cc

LU

tD

300

N

UJ 

Q

ID 200

100 -

VERTICAL

EXAGGERATION

20:1

52 54 56 58 60 62 64 66

DISTANCE, x, IN KILOMETERS

FIGURE 14.-Longitudinal profiles for case 6b that hasten the transition from continuing vigorous thinning to beginning to approach/^). The/fr) are shown for t-1, 2, 3, 4, 5, 6,11,16,..., 91, 96.
CONTINUITY EQUATION MODEL, PREDICTED DRASTIC RETREAT, COLUMBIA GLACIER, ALASKA A19

QC 

LU

LU

N

III

Q

500

400

300

200

100

i i T I i

VERTICAL

EXAGGERATION

20:1

I I I I I i l I I____l____I____I

52 54 56 58 60 62 64 66

DISTANCE, x, IN KILOMETERS

FIGURE 15. - Longitudinal profiles for case 7a that are compatible with an initial retreat rate of Xi = - 36 m/a. The/, (a;) are shown

for i-1, 2, 3, 4, 5, 6,11, 16,..., 91, 96.

500

400

0) 

QC 

LU

300

N 

UJ

P 200

100

\ i i i i i r

VERTICAL

EXAGGERATION

20:1

I____i____i____i____i____l____i____l____i____I____l____I

52 54 56 58 60 62 64 66

DISTANCE, x, IN KILOMETERS

FIGURE 16. - Longitudinal profiles for base 7b that are compatible with an initial retreat rate of Xi = - 54 m/a. The/(a;) are shown

for i-1, 2, 3, 4, 5, 6,11,16,..., 91, 96.
A20 STUDIES OF COLUMBIA GLACIER, ALASKA

54 56 

DISTANCE, x, IN KILOMETERS

58

FIGURE 17. - Final longitudinal profile, fm(x), for cases 1,8a, 8b, and 8c.

and so the terminus position in September 1981 was consistent with the 4.5-year average of X.

The retreat curve, X(t), has the same characteristic 

shape for all the variation cases and for all other alterations to the calculation. Approximately the first half of 

the period is consumed in very slow retreat (about 100 

m/a), then the rate accelerates rapidly to about 3,000 

m/a. The critical feature of the prediction is the timing 

of the transition from the present phase of disconnection 

from the shoal to the phase of drastic retreat. Describing the process after the onset of rapid retreat is transcended in practical importance by the interest in when 

the transition itself will occur. As shown in figure 22, the 

interval of likelihood of the transition is about 1982 ± 0.9. 

The mathematical crux of the difficulty is equation 20, 

which gives the crucial retreat rate, X, as the difference 

between two large, but poorly known, terms. The 

configuration of the bed is not well defined above the 

1978 terminus, and knowledge of the ice thickness is 

vital to scaling the fluxes for those two terms.

Another problem is that calving, ice flow at the terminus, and rate of retreat at the terminus (the fluxes in 

equation 20) show a strong seasonal fluctuation; as 

shown by figure 5, the annual retreat is about 45 m, and 

the seasonal variation of the terminus position is about 

300 m. The model described in this report uses annually

CO 

CC

LU 

LU

O

68

cc 66

CO

O

0. 

CO

62

60

58

/ U/ ©/ Yfc

//&' NB

56

1978.0

12

10

Q.

CO 

OC 

UJ

UJ

O 

6 00

3 

O

4 a

 * 

X

<

O

1980.0 1982.0 1984.0 1986.0 1988.0

FIGURE 18.-Effect on the solution of varying the bed topography, U(x), as shown in figure 11. Case 4a' uses a shallower bed, and case 4b*, a 

deeper bed, than case 1; in both cases, c is adjusted to maintain k=1980.4. Solid lines refer to terminus position (left axis); dashed lines refer 

to calving flux (right axis).
CONTINUITY EQUATION MODEL, PREDICTED DRASTIC RETREAT, COLUMBIA GLACIER, ALASKA

68

£66

ui

ui

O 

=! 64

Z

><

Z 62

O 

a. 60

58

56

A21

12

OC

g

10

a.

w

oc

O 

6 5

=)

O

Z

4 rf

 *

X

3

U.

1978.0 1980.0 1982.0 1984.0 1986.0 1988.0

FIGURE 19. - Effect on the solution of varying the shoal bathymetry, U(x), as shown in figure 12. Case 5a' uses a shallower shoal, and case Stf, a 

deeper shoal, than case 1; in both cases, c is adjusted to maintain k = 1980.4. Solid lines refer to terminus position (left axis); dashed lines refer 

to calving flux (right axis).

68

oc 66

UI t

§ 62

O

Q. 

0)60

OC 

UI 58

56

1978.0

12

10 Ui

Q.

V) 

OC 

UI

o

6 CD

4 0°

< 

O

1980.0 1982.0 1984.0 1986.0 1988.0

FIGURE 20. - Effect on the solution of varying the initial retreat rate^. Case 7a has-X", = - 36 m/a, and case 7b has-X", = - 54 m/a, compared with 

Xi= -45 m/a for case 1. Solid lines refer to terminus position (left axis); dashed lines refer to calving flux (right axis).
A22 

68

STUDIES OF COLUMBIA GLACIER, ALASKA

ui

ui 

2 

O 

d 64

X

2? 62 

O

0)

O

Q. 

CO

60

Z

cc UI58

O 

56

1978.0 1980.0 1982.0 1984.0 1986.0 1988.0

FIGURE 21. - Effect on the solution of varying the final longitudinal profile, fm(x\ as shown in figure 17. Case 8a' uses a steepened profile and an 

increased discharge, Q,^, and case Sd, a lower altitude, Zmk, than case 1; in both cases, c is adjusted to maintain k=1980.4. Solid lines refer to 

terminus position (left axis); dashed lines refer to calving flux (right axis).

12

10 AC iw m 

Q.

0)

oc

UI

 i

O

6 1

O

Z

4 c?

68

DC 

UI

2 M

Z

><

262

60

56

1978.0

12

oc 

10 £

0)

oc

UI 

UI

O 

« 

O

O

O

Z

< 

O

1980.0 1982.0 1984.0 1986.0 1988.0

FIGURE 22.-Resultant summation of the individual errors that may be present in the data used to describe Columbk Glacier. The likely early 

retreat, X^f), and the likely late retreat, XL(t), are shown, along with X(t) for case 1. The probability assigned to-X^O and XL(t) is the same as 

that ascribed to the variations of the data in cases 2-8. Observed values of X(t) are indicated by dots (see fig. 5).
CONTINUITY EQUATION MODEL, PREDICTED DRASTIC RETREAT, COLUMBIA GLACIER, ALASKA A23

averaged data to describe the glacier. The method for 

properly annualizing seasonal fluctuations is not obvious. At the very least, the seasonal cycles introduce an 

additional uncertainty in prediction.

REFERENCES CITED

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

glacer flow and its application to predict changes in the surge-type 

Variegated Glacier during its quiescent phase: Unpublished PhD 

thesis, University of Washington Geophysics Program, Seattle, 

Washington, 244 p.

Budd, W. F., and Jenssen, D., 1975, Numerical modelling of glacier 

systems: Proceedings of the Moscow Symposium, August 1971; 

Actes du Colloque de Moscou, aout 1971: LAHS-AISH Publ. No. 

104, 1975, p. 257-291.

Budd, W. F., Keage, P. L., and Blundy, N. A., 1979, Empirical studies 

of ice sliding: Journal of Glaciology, v. 23, no. 89, p. 157-170.

Glen, J. W., 1954, The creep of polycrystalline ice: Proceedings Royal 

Society, Ser. A, v. 228, p. 519-538.

Kollmeyer, R. C., Motherway, D. L., Robe, R. Q., Plate, 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.

Mayo, L. R., Trabant, D. C., March, Rod, and Haeberli, Wilfried, 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.

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 Report 

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.

Nye, J. F., 1952, The mechanics of glacier flow: Journal of Glaciology, 

v. 2, no. 12, p. 82-93.

Post, Austin, 1975, Preliminary hydrography and historic terminal 

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

Hydrographic Investigations Atlas 559, 3 sheets.

Sikonia, W. G., Finite element glacier dynamics model applied to 

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

Paper 1258-B [in press].

Sikonia, W. G. and Post, Austin, 1980, Columbia Glacier, Alaska: 

recent ice loss and its relationship to seasonal terminal embayment, thinning, and glacier flow: U.S. Geological Survey Hydrologic Investigations Atlas 619, 3 sheets.

ilrU.S. GOVERNMENT PRINTING OFFICE: 1982-361-614/4