## Introduction

Because of the structural complexity of glaciers and ice caps, direct measurements of them are necessary in order to reach a quantitatively pertinent flow law. As Lliboutry (1970) has pointed out, the measurements “are by no means simple since three distinct problems must be solved simultaneously (ice flow law, boundary conditions at the bottom, dynamic problem), and only the whole solution may be checked by field data”. In an effort to construct a model which could be used to relate solutions of three-dimensional, time-dependent glacier flow to observations, Campbell and Rasmussen (1970) assumed a viscous flow law for ice, a simple (Bodvarsson, 1955) basal flow law for the bottom boundary condition, and generated a system in which solutions could be obtained for any given glacier whose bedrock topography and mass balance are known.

Both the ice and basal flow laws were simple compared with existing flow laws (Glen, 1955, 1958; Lliboutry, 1968; Nye, 1960, 1963[a], [b], [c], 1963[a], [b], [c]; Weertman, 1957, 1964, 1969) and were chosen so as to render the equations mathematically tractable for numerical solution. Although the flow laws used were simple, the model was more complex than many existing models in that it was three-dimensional and time-dependent, it included the effect of longitudinal and lateral stress gradients, and it partitioned the ice viscosity from the basal friction allowing them to be parameterized independently. Although the solutions obtained were for an ideally smooth bedrock topography and mass balance versus altitude curve, the equilibrium profiles for selected sets of the ice and basal friction parameters resembled those of real glaciers. Campbell and Rasmussen (1969) further applied the model to glacier waves, surges, and surge recoveries, and again found solutions that resembled the behavior of real glaciers. At the time this model was developed the authors were aware that although the solutions resembled the behavior of real glaciers, the model would not be quantitatively accurate until a more complex flow law, one in keeping with the findings of the above named authors, was incorporated into it. Their *modus operandi* was to get the simpler model working and then go on to more complex ones.

The problem facing the authors was in what way the first model ought to be improved—should they adopt a more complex ice viscosity, or a more complex basal friction? Various measurements of temperate glacier ice, for example those of Colbeck (unpublished), suggest that at the stresses normally found in temperate glaciers, the Newtonian viscous approximation for ice flow is acceptable, whereas considerable controversy exists over the mechanisms of basal flow. Therefore, it was decided to compare solutions, using the same idealized valley glacier and mass balance to which the earlier model was applied, for three models incorporating different basal flow laws: the Bodvarsson law used in the earlier model, a vertically integrated form of Glen’s (1955) flow law used extensively by Nye (1960, 1963[a], [b], [c]) and others, which for the sake of simplicity we will refer to throughout the text as the Glen law, and another form of Glen’s law, involving certain assumptions concerning the stress distribution, used by Budd (1969), Nye, and others, which again for the sake of simplicity we will refer to as the Budd law, because in the literature he has most frequently used it.

These flow laws, which involve three parameters, are in turn coupled with the Campbell and Rasmussen momentum equation, obtained by vertically integrating the Navier–Stokes equation, which adds one more parameter to the problem, the ice-to-ice viscosity *v*. Each of these momentum equations is then solved simultaneously with the continuity equation. The greater part of this study involved variation of parameters associated with basal friction. The value of *v* we used (10^{15} cm^{2}/s) proved to be the optimum value in our earlier studies and is close to the value obtained by Raymond (1971 [a], [b]) in his study of the Athabasca Glacier.

The influence of each of the three basal flow parameters was first studied by modelling a steady-state valley glacier and noting how the glacier shape and flow varied for the equilibrium solution for each set of parameters. Then the effect on non-steady behavior was examined by depressing the basal friction parameter to five per cent of its initial steady-state value for a period of one year, thereby inducing a glacier surge. Finally the basal friction parameter was restored to its initial value and the recovery of each glacier was studied for one hundred years.

## Mathematical Model

As shown by Campbell and Rasmussen (1970), the Navier–Stokes equation, from which inertial and acceleration terms have been omitted, may be written

where **V** = *u*
**i** + *v*
**j** + *w*
**k** is the velocity vector, *v* is the kinematic viscosity, **g** = *g*_{x}
**i** + *g*_{y}
**j** + *g*_{z}
**k** is the gravity vector (with magnitude *g*), *p* is the pressure, *ρ* the density. If Equation (1) is applied to a glacier and ∇^{2} is the two-dimensional Laplacian operator, and the coordinate system is chosen so that the *z*-axis is normal to the bed of the glacier, then Equation (1) may be written in component form as

where *τ*
_{
x
} and *τ*
_{
y
} are the shearing components parallel to the bed, and the velocity component *w* normal to the bed has been neglected.

The hydrostatic equation (3) may be integrated vertically, yielding the hypsometric equation

where *h* is the glacier thickness, measured normal to the bed, and *p*
_{
a
} is the atmospheric pressure.

A volume transport vector **Q** = *Q*_{x}
**i** + *Q*_{y}
**j** may be defined as

where

Term by term, the *x*-component equation (2) may be integrated vertically to produce a transport equation. By using Equation (4), the first term may be integrated as follows:

The gravity term can be integrated directly:

Assuming *v* to be independent of *z*, the horizontal viscous term is integrated,

Since the shearing stress vanishes at *z* = *h*, the vertical viscous term is integrated,

where *τ*
_{
x
}′ is the shearing stress at the bed.

The vertical integration of the *y*-component equation (2) leads, then, to the component form of the transport equation:

If the transport equation is to be applied to an arbitrarily shaped glacier bed, it may conveniently be transformed so that *h* is measured along the Earth-vertical. This facilitates the plotting of computation results, and accumulation or ablation is specified in this direction. The transformation is accomplished by invoking the relation

The symbol *h* in Equations (4)–(6) indicates *h*
_{(normal)} ; in all appearances hereafter it indicates *h*
_{(vertical)}, and the relevant coordinate system is that in which *z* is measured along the Earth-vertical, positive upward, with *x* and *y* in the Earth-horizontal plane.

The transport equation, so transformed, is

If tan *α*
_{
x
} and tan *α*
_{
y
} are the components of the bed’s inclination to the horizontal, the gravity components are given by

where

If a functional relationship can be established between **Q** and **
***τ*
′, then Equation (7) can be solved for **Q**, presuming *h* to be known. Suggested forms of this relationship, in which the bed friction coefficients *A* are specified constants, include

estimates of the exponent *n* ranging from 1 to 3. Setting **Q** = *h*
**V**, these may be generalized to

by which is meant

where sgn (*x*) is + 1, 0, −1 as *x* is positive, zero, negative.

The Bodvarsson law has *k* = 0, the Glen law has *k* = 1, and the Budd law has *k* = 2. Also *A* is seen to have the units g^{
n
} cm^{
k−n−2} s^{1−2n
} and, therefore, should be denoted *A*
_{
k
}(*n*). Substituting for *τ*′ from Equation (6) yields

If *a*(*x, y, t*) is the specified mass-balance forcing function, the continuity equation for incompressible flow can also be integrated vertically to yield

Equations (11) and (12) now constitute a system of three algebraic equations in three unknowns *Q*
_{
x
}, *Q*
_{
y
}, and *h*. Assumption of a subglacial geometry determines the *x, y* distribution of tan *α*
_{
x
} and tan *α*
_{
y
} and therefore, through Equations (8), of the components of the gravity vector.

Mathematically, then, the glacier is treated as an array of vertical columns of unit area and varying height. The ice-to-ice friction on the sides of the columns is controlled by the viscosity coefficient *v*, whereas the ice-to-bed friction is controlled separately by the bed friction coefficient *A*.

## Numerical Solution

Equations (11) and (12) are solved by a forward-difference method on a square finite-difference grid with spacing *S*, and whose grid points are denoted by standard matrix notation. If (*x, y*) are the coordinates of grid point (*i, j*), then (*x*+*S, y*) are the coordinates of point (*i, j*+1), and (*x*, *y*–*S*) of point (*i*+1, *j*).

and the continuity equation (12) can be approximated by

where

The equations for *Q*
_{
x
} and *Q*
_{
y
} are solved by over-relaxation, in about three passes when *v* ≠ 0, in (exactly) one pass when *v* = 0. The boundary conditions were conveniently managed by choosing the solution region to be sufficiently large that the glacier lay wholly within it. Then, since *h* = 0 all around the glacier, and by requiring **Q** = 0 whenever *h* = 0, it is possible to impose a no-slip boundary condition precisely at the glacier margin without specifying the location of the margin in advance.

The calculations were performed on a CDC-6400 computer by a modular program written in FORTRAN IV. An empirical stability examination showed

where ∆*t* is given in days for *A*
_{
k
} (*n*), and *h* in c.g.s. units but *S* in m. The calculation rate is given approximately as

where *M* is the number of (*h* > 0) points in the solution region, and *N* is the average number of relaxation passes per time step. For the cases described in this paper, the time step varied from 10 to 40 d for steady-state solutions and surge recoveries, and from to 2 d for surges. The calculation rate was about one second per time step.

## Steady-State Solutions

Shown in Figure 1 are the glacier bed and annual mass balance used in all calculations. The glacier bed has a steep (30°) upper slope joined by a short, smooth transition to a shallow (8°) lower slope, with a constant, U-shaped transverse cross-section that is symmetric about the glacier center-line. The annual mass balance, unchanging in time, varies from +2 m/year at the top of the accumulation zone to about –8 m/year at the terminus, depending on its precise position. A 200 m grid spacing was used for all calculations.

## Steady-state solution for Glen flow law (n = 3).

For any rheology (*k, n* pair) considered, the steady-state solution is strongly dependent on the numerical value of the coefficient *A*. Therefore, to examine the differences between rheologies (separately from differences due to variations in the value of *A*), the *A*-values were chosen so that the glacier volumes of the steady-state solutions were the same for all seven rheologies. The value for *A*
_{0}(1) (Bodvarsson, *n* = 1) was taken to be 10^{6}, as by Campbell and Rasmussen (1970), which produced a steady-state volume of 3.375 km^{3}. Shown in Table I are the values for the other six rheologies that produce this same steady-state volume. Also shown in Table I are the sensitivities of each steady-state volume to unit variations both in *A*
_{
k
}(*n*) and in *v*.

## Table. I. Dependence of steady-state volume on *A*_{k}
(*n*) and *v*

Figure 1 gives the steady-state distributions of glacier thickness and volume transport for the rheology defined by the use of the Glen flow law (*k* = 1) with *n* = 3. It exhibits the familiar pattern of concave transverse cross-sections in the accumulation zone, with converging streamlines, and convex transverse cross-sections in the ablation zone with diverging streamlines.

When constrained to have the same volume, the steady-state solutions for the other six rheologies have glacier thickness distributions so similar that their differences are not readily perceived when plotted at the scale of Figure 1. Hence these differences are shown in Figure 2, which gives the center-line glacier thickness profiles for each of the seven rheologies.

## Steady-state glacier thickness center-line profiles comparing seven different rheologies.

At steady-state, Equation (12) becomes

which, when integrated over the areal extent of the glacier, gives **Q** solely as a function of *a*(*x,y*), independent of any direct influence from Equation (11). The slight differences of the steady-state **Q** distributions from rheology to rheology are due only to the differences in the areal extent of the seven steady-state solutions which, as shown in Figures 5–11, are quite small. Figure 1 also contains the distribution of **Q** for the rheology *k* = 1, *n* = 3.

## Surge and Recovery Solutions

Each steady-state solution was caused to surge by carrying out the calculation for one year with *A*
_{
k
}(*n*) depressed to 0.05 of the value that produced the steady-state solution. Then, *A*
_{
k
}(*n*) was restored to its former value, and the first 100 years of its reversion to steady-state was examined.

The reduction of *A*
_{
k
}(*n*) increases the transport for a given basal stress, permitting the glacier to flow more rapidly downhill. Hence, the surge is marked by height falls in the upper part of the glacier, height rises in the lower part, and a slight advance of the terminus. Also, since

over the extent of the steady-state glacier, its advance appends a region in which the mass-balance is negative, resulting in a net loss of glacier mass.

When the steady-state *A*
_{
k
}(*n*) value is restored, the glacier begins to reassume its steady-state distribution, but in a more complicated manner than the simple pattern of height changes that occurred during the surge. In fact, the glacier continues to lose mass for several more tens of years. The resumption of the steady-state distribution begins at the top and proceeds down the glacier.

Figure 3 shows the variation of glacier mass during both the one-year surge and the 100-year recovery for each of the seven rheologies. Equation (11) may be roughly regarded as a polynomial in *h* of degree *k* + *n*, which serves as a crude index characterizing the dynamic response of the glacier. The speed of recovery generally increases with *k* + *n*.

## Variation of glacier mass during surge and recovery comparing seven different rheologies.

Another indication of differences in dynamic response is given by Figure 4. It shows, for each of the seven rheologies, the centerline profile of the down-slope volume transport at the end of the one-year surge, while the *A*
_{
k
}(*n*) value is still depressed. The volume transport during the surge decreases with *k*+*n*.

## Down-slope volume transport center-line profiles at the end of one year surge comparing seven different rheologies.

Figures 5–11 show, each for one or another of the rheologies, the height changes from steady-state during surge and recovery. The transient behavior of the height changes along the glacier center-line is exhibited both during the surge and during the recovery. Also shown is the distribution, over the entire glacier, of the net height change at the end of the surge, as well as the glacier margin both at steady-state and at the end of the surge.

## Height changes during surge and recovery for Bodvarsson flow law (n = 1).

## Height changes during surge and recovery for Glen flow law (n = 1).

## Height changes during surge and recovery for Glen flow law (n = 2).

## Height changes during surge and recovery for Glen flow law (n = 3).

## Height changes during surge and recovery for Budd flow law (n = 1).

## Height changes during surge and recovery for Budd flow law (n = 2).

## Height changes during surge and recovery for Budd flow law (n = 3).

The surges follow a relatively simple pattern of continuously falling heights in the upper part of the glacier, continuously rising heights in the lower part, and a slight continuous advance of the terminus. The glacier bed at the top of the glacier is exposed during the surge for rheologies Bodvarsson, *n* = 1, and Glen, *n* = 1. Presumably it would also be exposed for the other rheologies, were the surge continued for a sufficient time. The magnitude of the terminus advance, as well as of the maximum height falls and rises, generally decreases with *k*+*n*. Also, as *k*+*n* increases, the height rises tend more to a double maximum, one occurring at the base of the steep part of the bed, the other at the terminus.

The recovery is more complicated than a simple reversal of the pattern of height changes taking place during the surge. At steady-state the down-slope transport values in the lower part of the ablation zone are in equilibrium with the mass-balance. The excess of mass here at the end of the surge produces higher transport values, which continue to supply mass to the region in advance of the terminus, where the mass balance is strongly negative. This results in losses of glacier mass until the transport falls below the steady-state values, in 45–65 years depending on the rheology.

At the end of the surge, the wave formed by the excess of mass at the down-stream end of the steep part of the bed and the deficiency of mass in the accumulation zone move down the glacier with a wave speed of about 100 m/year, rapidly diffusing. The mass deficiency part of the wave is followed by a minor, short-lived positive anomaly at the base of the steep part of the bed. When the mass deficiency passes the terminus, after 50–70 years depending on the rheology, the entire glacier surface is at or below its steady-state position.

The final stage of recovery consists of a simple pattern of continuous height rises until steady-state is attained, which occurs first at the top of the accumulation zone and then progresses down-glacier. The down-glacier speed of the point of attainment of steady-state increases with *k*+*n*.

The erratic behavior at the terminus that occurs during recovery does not occur during the surge, when the glacier is moving an order-of-magnitude faster. Neither does it ever produce a glacier surface at the terminus that is irregular when plotted as glacier thickness, rather than as a departure from steady-state. This suggests that it is a real, crack-of-the-whip phenomenon at the terminus, not simply a symptom of numerical volatility.

*
*A*
_{
k
}(*n*) value producing 3.375 km^{3} steady-state and *v* = 10^{15} cm^{2}/s.

## References

Bodvarsson, G.
1955. On the flow of ice-sheets and glaciers. Jökull, Ár 5, p. 1–8.

Budd, W. F.
1969. The dynamics of ice masses, ANARE Scientific Reports. Ser. A(IV). Glaciology. Publication No. 108.

Budd, W. F., *and*
Jenssen, D. In press. Numerical modelling of glacier systems.

Campbell, W. J., *and*
Rasmussen, L. A.
1969. Three-dimensional surges and recoveries in a numerical glacier model. Canadian Journal of Earth Sciences, Vol. 6, No. 4, Pt. 2, p. 979–86.

Campbell, W. J., *and*
Rasmussen, L. A.
1970. A heuristic numerical model for three-dimensional time-dependent glacier flow. [Union Géodésique et Géophysique Internationale. Association Internationale d’Hydrologie Scientifique.] [International Council of Scientific Unions. Scientific Committee on Antarctic Research. International Association of Scientific Hydrology. Commission of Snow and Ice.] International Symposium on Antarctic Glaciological Exploration (ISAGE) Hanover, New Hampshire, U.S.A., 3–7 September 1968, p. 177–90.

Colbeck, S. C. Unpublished. The flow law for temperate glacier ice. [Ph.D. thesis. University of Washington Seattle, Washington, 1970.]

Glen, J. W.
1955. The creep of polycrystalline ice. Proceedings of the Royal Society, Ser. A, Vol. 228, No 1175, p. 519–38.

Glen, J. W.
1958. The flow law of ice: a discussion of the assumptions made in glacier theory, their experimental foundations and consequences. Union Géodésique et Géophysique Internationale. Association Internationale d’Hydrologie Scientifique. Symposium de Chamonix, 16–24 sept. 1958, p. 171–83.

Lliboutry, L. A.
1968. Théorie complète du glissement des glaciers, compte tenu du fluage transitoire. Union de Gédésie et Géophysique Internationale. Association Internationale d’Hydrologie Scientifique. Assemblée générale de Berne, 25 sept.–7 oct. 1967. [Commission de Neiges et Glaces.] Rapports et discussions, p. 33–48.

Lliboutry, L. A.
1970. Current trends in glaciology. Earth Science Reviews, Vol. 6, No. 3, p. 141–67.

Nye, J. F.
1960. The response of glaciers and ice-sheets to seasonal and climatic changes. Proceedings of the Royal Society, Ser. A, Vol. 256, No. 1287, p, 559–84.

Nye, J. F.
1963[a]. On the theory of the advance and retreat of glaciers. Geophysical Journal of the Royal Astronomical Society, Vol. 7, No. 4, p. 431–56.

Nye, J. F.
1963[b]. The response of a glacier to changes in the rate of nourishment and wastage. Proceedings of the Royal Society, Ser. A, Vol. 275, No. 1360, p. 87–112.

Nye, J. F.
1963[c]. Theory of glacier variations. (In Kingery, W. D., ed. Ice and snow; properties, processes, and applications: proceedings of a conference held at the Massachusetts Institute of Technology, February 12–16, 1962. Cambridge, Mass., M.I.T. Press, p. 151–61.)

Nye, J. F.
1965[a]. The flow of a glacier in a channel of rectangular, elliptic or parabolic cross-section. Journal of Glaciology, Vol. 5, No. 41, p. 661–90.

Nye, J. F.
1965[b]. The frequency response of glaciers. Journal of Glaciology, Vol. 5, No. 41, p. 567–87.

Nye, J. F.
1965[c]. A numerical method of inferring the budget history of a glacier from its advance and retreat. Journal of Glaciology, Vol. 5, No. 41, p. 589–607.

Raymond, C. F.
1971 [a]. Determination of the three-dimensional velocity field in a glacier. Journal of Glaciology, Vol. 10, No. 58, p. 39–53.

Raymond, C. F.
1971 [b]. Flow in a transverse section of Athabasca Glacier, Alberta, Canada. Journal of Glaciology, Vol. 10, No. 58, p. 55–84.

Weertman, J.
1957. On the sliding of glaciers. Journal of Glaciology, Vol. 3, No. 21, p. 33–38.

Weertman, J.
1964. The theory of glacier sliding. Journal of Glaciology, Vol. 5, No. 39, p. 287–303.

Weertman, J.
1969. Water lubrication mechanisms of glacier surges. Canadian Journal of Earth Sciences, Vol. 6, No. 4, Pt. 2, p. 929–42.