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 Reference LliboutryLliboutry (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, Reference Campbell and RasmussenCampbell and Rasmussen (1970) assumed a viscous flow law for ice, a simple (Reference BodvarssonBodvarsson, 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 (Reference GlenGlen, 1955, Reference Glen1958; Reference LliboutryLliboutry, 1968; Reference NyeNye, 1960, Reference Nye1963[a], Reference Nye[b], Reference Nye and Kingery[c], Reference Nye1963[a], Reference Nye[b], Reference Nye and Kingery[c]; Reference WeertmanWeertman, 1957, Reference Weertman1964, Reference Weertman1969) 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. Reference Campbell and RasmussenCampbell 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 Reference GlenGlen’s (1955) flow law used extensively by Reference NyeNye (1960, Reference Nye1963[a], Reference Nye[b], Reference Nye and Kingery[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 Reference BuddBudd (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 (1015 cm2/s) proved to be the optimum value in our earlier studies and is close to the value obtained by Reference RaymondRaymond (1971 [a], Reference Raymond[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 Reference Campbell and RasmussenCampbell 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 = gx i + gy j + gz 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 = Qx i + Qy 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 s1−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.
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 106, as by Reference Campbell and RasmussenCampbell and Rasmussen (1970), which produced a steady-state volume of 3.375 km3. 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.
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.
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.
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.
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.
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.
