Introduction

During the past decade several authors have used the finite-element technique to model numerically the flow of an ice mass; Reference NguyenNguyen (unpublished), Reference SchmidtSchmidt (1977), Reference IkenIken (1977, Reference Iken1981), Reference Hooke, Hooke, Raymond, Hotchkiss and GustafsonHooke and others (1979), Reference Emery, Emery, Hanafy, Holdsworth and MirzaEmery and others (1979), Reference MacAyeal and ThomasMacAyeal and Thomas (1982, Reference MacAyeal and Thomas1984), Reference SikoniaSikonia (1982), Reference Fastook and SchmidtFastook and Schmidt (1982), Reference EchelmeyerEchelmeyer (1983), and Reference RaymondRaymond (1983). Applications have ranged from the overall motion of the Ross Ice Shelf (Reference MacAyeal and ThomasMacAyeal and Thomas, 1982, Reference MacAyeal and Thomas1984) to the development of a crack near an ice cliff (Reference IkenIken, 1977), and from a fast-moving, calving temperate glacier (Reference SikoniaSikonia, 1982) to a slow-moving, cold ice cap (Reference Hooke, Hooke, Raymond, Hotchkiss and GustafsonHooke and others, 1979).

The method is thus well established as a powerful tool to analyze many different topics in ice dynamics. Much of the reason for this popularity lies in its ability to handle boundary conditions easily, in particular, the complex geometric shapes of typical ice masses in the real world. It also lends itself to a modular programming approach, which makes it readily adaptable to different situations.

None of the published works, however, addresses specifically the time-dependent modeling of an arbitrarily shaped ice mass. Such generality, at least for plane (“two-dimensional”) flow, was one of two primary objectives of this paper. The other was to develop procedures to allow the model to be used by someone who had only a cursory knowledge of the details of the method, and to require, as much as possible, only glaciological data as input. These objectives were achieved and the resulting software will be described extensively in a forthcoming open-file report of the U.S. Geological Survey. This will make a reasonably general finite-element ice-flow model available to the glaciologist who has neither the time nor the inclination to commit himself to a long and often frustrating development period.

This model is based on an earlier one used by Reference SikoniaSikonia (1982) on the lowermost section of Columbia Glacier, Alaska. Numerous minor improvements were made to achieve the above objectives but only the two major ones will be described here: modeling of a “typical” glacier end, where the ice thickness tends to zero, and an automatic procedure for generating an appropriate “mesh” of “finite elements” for approximating the specified geometry. Neither of these features appears in previous works but their inclusion enables analysis of different problems to be done without any program modification or extensive data preparations.

The Finite-Element Method

A disadvantage of the finite-element method is that it is conceptually more difficult to understand than its finite-difference counterpart for solving differential equations. A description detailed enough to give a working knowledge of the method covers many chapters in a book and is beyond the scope of this paper. Many text-books describe the method in ample detail, for example, Reference ZienkiewiczZienkiewicz (1977) or Reference Lapidus and PinderLapidus and Pinder (1982). Nevertheless, a brief description of the method and its terminology is given to serve both as a basis for subsequent discussions and as a partial review of the subject for the glaciological literature, which has not been done previously.

The region of interest is subdivided into a set of discrete sub-regions, called “finite elements”, and the unknown functions, u(x_i ), are approximated within each element by some linear combination of interpolating polynomials, f_j (x_i ), referred to as “shape” or “basis” functions:

where x_i represents the spatial coordinates and the U_j are coefficients to be determined. The degree of these basis functions is specified, based on the desired physical size of the elements and upon some a priori knowledge of how the functions are likely to vary within each element.

A set of N points (or “nodes”) is then chosen within each element and the basis functions selected so that function f _j has a value of one at node j, and a value of zero at all other nodes. When this is done, U _j becomes the value of the unknown function at that node. The number and distribution of the nodes are dictated by the maximum degree of the basis functions chosen. Since polynomials are used, all derivatives in the original differential equations are now easily evaluated when Equation (1) is substituted for the unknown functions.

By utilizing a method of “weighted residuals”, the differential equations can be converted to a set of linear simultaneous equations, one for each unknown coefficient U _j , which can then be solved using standard methods. In the process, the differential equations are changed to an integral form, with the integrations performed over the entire domain of the element (or along its surface for boundary conditions). The solution obtained thus approximates the true one in an integral sense over each element, rather than only at discrete points as in a finite-difference formulation. This implies that the elements can be made bigger, that is, more “finite”, and thus hopefully fewer nodes are required for an adequate solution to the problem.

Furthermore, since an integral over a region is equal to the sum of integrals over a set of sub-regions, the equations for all elements can be combined into a “global” set of linear simultaneous equations. Although this results in a large number of unknowns and equations (hundreds or thousands), the solution for the entire ice mass can still be obtained with standard methods for such equations.

Basis functions are also used to handle curvilinear element shapes rather than being restricted to purely rectangular ones. This model utilizes “isoparametric” elements, in which basis functions of the same degree are used for both the unknown function interpolation and the coordinate transformation (Reference ZienkiewiczZienkiewiez, 1977).

The heart of the finite-element method is the basis functions. They specify how the functional variation across the element is to be approximated. The reader is cautioned to keep Equation (1) in mind and not to interpret a basis function directly as a “conventional” interpolating polynomial. This is particularly true for two-dimensional basis functions (which are used here); these are formed by combining basis functions defined separately for the two coordinate directions. This allows “mixed-order” basis functions to be used, for example, linear in one direction but quadratic in the other. In fact, the order can be different along opposite edges of an element, which allows a transition to be made, if desired, between different degrees of functional variation.

The basis functions are also defined separately for each unknown function. Thus, by varying the number and distribution of nodes in an element and, whether or not a particular function is assigned to a particular node, many complex physical situations can be simulated. In addition, the elements and nodes can be distributed non-uniformly over the entire region of interest so that the calculations can be concentrated where the unknowns are varying more rapidly.

Another advantage of the particular manner in which basis functions are developed is that it is possible to write computer code which can automatically select different orders of basis functions, depending upon the input (Reference Taylor and ZienkiewiczTaylor, 1977). The model does not have to be tied rigidly to any particular type; the one described here, in fact, can alternate between various basis functions by changing a single input variable.

Application to Ice Flow

The particular formulation of the method applicable to the flow of viscous materials has been given by Reference ZienkiewiczZienkiewiez (1977), and the adaptation of this for ice flow that is used here has been given by Reference SikoniaSikonia (1982). Velocity and pressure are treated as the dependent variables; pressure must be added because in a fluid the deformation rates are related only to the deviatoric stresses and so the pressure term is necessary to complete the stress solution.

Mesh Generation

A set of “mesh-generation” routines was developed to create automatically from glaciological data not only suitable coordinates for the mesh nodes but also the boundary conditions and all the other input required by the finite-element routines. The input data consist of the bed elevation, the surface elevation, and the balance rate. This is a significant improvement over earlier models in which such information must be generated manually. It simplifies enormously the preparation of the data for the program. It eliminates any human bias which might be introduced in deciding where to place the mesh nodes. Also, if the technique is suitably chosen, it can optimize, at least to some extent, the node distribution specifically for the analysis of ice flow, and can then maintain this optimization regardless of how the geometry of the ice mass changes with time.

To achieve this optimization, several methods were tried. Most were abandoned because they were not robust enough to work on a wide variety of surface and bed geometries. The one finally adopted is outlined in detail in (Figure 1) the Appendix. The X -coordinates of the nodes are determined by adjusting iteratively the widths of each element column so that the elements are narrower where the geometry changes are greater, using as criteria a combination of the mean thickness change and mean bed curvature across the column. The Z-coordinates of the nodes are determined by requiring each element to have an equal fraction of the total deformational ice velocity occurring across it, assuming that plane laminar flow were taking place. The algorithm allows specification of the overall amount of adjustment that the element widths undergo by using a parameter, δ (see Appendix). The maximum value of δ that can be obtained will vary from one glacier to another since the amount of adjustment must remain less than the typical wavelengths in the surface and bed topography.

Figure 1 shows the results of such an adjustment on South Cascade Glacier for a value of δ = 100%, which is the value used throughout the remainder of this paper. Values up to 200% were easily obtained but higher values could only be approximated; this is thus the limit for this particular topography. The slope of the element boundaries is in general discontinuous between elements since the transformation from curvilinear coordinates is done separately for each element. Since quadratic basis functions are used for the transformation, each element boundary is effectively approximated by some quadratic shape.

Zero-Thickness Ends

If loss of material at the terminus of an ice mass is dominated by ablation, rather than by calving, the ice thickness tapers gradually to zero. This condition also occurs at the upper end provided the ice is not thick enough to “cap” the high point of the bed topography and flow away in either direction. These situations are quite common, particularly for temperate glaciers, and are referred to here as “zero-thickness ends”.

The finite-element method allows a transition to be made between different element shapes. Hence this condition could be realistically modeled by using triangular elements at the end, with quadrilaterals elsewhere. An easier way, however, which appears to work equally well in practice, is to truncate the ice artificially so that it has a small, but non-zero, specified thickness at the end of the mesh. This short vertical face is assigned the same number of nodes as elsewhere, which allows the same element types to be used throughout the ice mass. These highly disproportionate quadrilateral elements do not cause any problems, which is ample testimony to the enormous flexibility of the basis-function concept.

Hydrostatic ice pressure is applied along this face as a boundary condition to approximate the “missing” triangular wedge of ice. Otherwise, the nodes here are free to move just as they can elsewhere.

For geometrically static modeling, this would be sufficient. When time-dependency is added, however, additional considerations are necessary. First, the computed velocity at the end of the ice cannot be used to calculate the amount of advance or retreat of the end in the same way that it was used to project the surface at each time step, because this fails to account for the effective advance or retreat caused solely by the balance rate. This effect, in fact, is the only way a zero-thickness no-sliding end can advance or retreat, since the velocity is zero by definition in this case. The correct procedure is a “bed/surface-intersection” method in which the new surface is formed by using the solution for

as outlined previously, and then the intersection of this with the bed is located and used as the new position. This requires extrapolation for an advancing end and interpolation for a retreating one.

Such considerations apply equally to the head or terminus of the glacier, provided they are both free to advance or retreat. Often, however, it is desired to keep the head of the ice mass fixed. This can be easily done by requiring

to be zero here, or by maintaining a fixed X-coordinate for the start of the mesh, but, in either case, the balance rate must also be zero at this point. If not, an unrealistic bump forms near the end. This point is easily overlooked because measured balance data typically do not extend right to the very head of an ice mass and positive non-zero values measured in the accumulation area are used to extrapolate to the end. If the head is to remain stationary, the balance curve must be tapered off to zero at the end. A plausible justification for this is that wind effects probably dominate at this point and work in such a way as to maintain the net balance close to zero.

Even this condition must be treated cautiously. Maintaining a fixed head is only justified if the adjacent ice surface is already close to equilibrium with the balance data. If not, the surface will attempt to adjust to a significantly different configuration as time progresses. Fixing the head effectively holds back the ice artificially and eventually an oscillation in the surface is induced. The correct way of modeling such a situation is, of course, to free the head completely and allow it to move like the terminus, which is, in effect, exactly what a bergschrund represents.

Time-Dependency Considerations

The discussion in the previous section illustrates an important point. Time-dependent modeling subjects the numerical calculations to very rigorous testing. In many cases, the results are very sensitive to small errors or slight inconsistencies in the data, and a solution to the equations is immediately impossible. There are also, however, cases where the opposite is true, and only after the effect has accumulated over many time steps does the error become apparent.

An example of this occurred while developing the surface/bed-intersection procedure mentioned above. The surface and bed data sets consist of a series of discrete points in which interpolation is performed to obtain a value at an arbitrary X-coordinate. Originally, linear interpolation was used. It was some time before it was realized this was the cause of a perplexing, erratic advance/retreat of the end of the ice, including abrupt “jumps” off the end of the defined data. Cubic spline interpolation, in which there are no discontinuities in slope, corrects this problem and produces smoothly varying advance or retreat. This method should be used for time-dependent modeling whenever interpolation is required in a set of tabular data points.

A far more serious example is shown in Figure 2, a time simulation of South Cascade Glacier when the balance rate is uniformly decreased (a Case II run described in the next section). This undesirable “washboard” effect, which is clearly a numerical artifact since it is tied to the node locations, is similar to oscillations which are known to occur in other finite-element flow problems (Reference Walters and CareyWalters and Carey, 1983), and which have their counterpart in finite-difference methods. Large internode oscillations in the values of the dependent variables, such as velocity and pressure, can develop when certain combinations of basis functions are used. In particular, solutions of the Navier–Stokes equation in two-dimensions result in spurious oscillations when both velocity and pressure have quadratic basis functions.

Fig. 2. The numerical oscillation induced by using the same order of basis functions for the

and velocity unknowns (quadratic). It can be detected at the end of the Case 1 simulation (t = 25 a) but it is really only obvious after allowing the model to run for many time steps. The insets show the calculated unknowns at t = 25 a. U is the X-component of velocity, W the Z-component, P the pressure, and

the rate of change of thickness. The vertical scale for each inset is: U, from 0 to 16 m a ⁻¹; W, from –4 to +1 m a⁻¹ ; P, from –0.5 to +2.0 bar: and

from –1 to +1 m a⁻¹. The inter-element node columns are indicated by the vertical lines in the main figure and by the tick marks in the insets.

In this model, pressure has a different basis function (linear) in one direction, so the situation is not exactly the same. Instead, the oscillation develops in the

unknown (see inset), which, in the original formulation (Reference SikoniaSikonia, 1982), utilizes quadratic basis functions. If linear basis functions are used for the oscillations vanish (Fig. 3). Any variations in surface slope now correlate well with the bed topography. It is, of course, only natural that this change should eliminate the problem since the intra-element nodes are not used for when linear interpolation is used. Unfortunately, the change means the resolution of surface features is diminished; however, this can always be made up for, if necessary, by increasing the number of elements longitudinally.

Fig. 3. Case I and Case II simulations. Case I shows the adjustment that the original surface (dashed line at t = 0) must undergo to be in equilibrium with the 1975 balance. The result is the heavy solid line at t = 25 a. Case II shows the changes that this shape undergoes after the balance rate is then uniformly decreased by –0.25 m a⁻¹. The lines are 50 years apart; the new equilibrium shape takes about 300 years.

Care must be taken, therefore, to keep the order of the basis functions sufficiently “mixed” that such oscillations do not appear. Time-dependent modeling towards steady-state conditions is the only way of testing for this condition. Reference SikoniaSikonia (1982) did not detect this problem because conditions were not close enough to steady-state on the rapidly changing Columbia Glacier.

The Flow of South Cascade Glacier

South Cascade Glacier, in Washington State, is a small, steep valley glacier which has been well studied over the past few decades. Data for the model are known over its entire length and it is a good example of an “arbitrary” ice mass. These data are shown in Figure 1; the bottom topography is from radar sounding done in 1975 (Reference HodgeHodge, 1979), the surface topography is from a high-resolution map made from aerial photography flown in 1977, and the balance rate is from the 1975 balance year (personal communication from R. Krimmel, 1984).

All simulations assume a constant temperature (0° C) and density (917 kg m⁻¹), which are valid conditions for this glacier. They also assume no sliding and a constant shape factor of 0.5; the latter is introduced merely as a velocity-scaling factor to estimate the drag exerted by the valley walls (Reference PatersonPaterson, [^c1981]). Both of these conditions were required so that calculated velocities agreed approximately with observed ones; the implications of this are discussed later. The flow-law parameters used are those recommended by Reference PatersonPaterson ([^c1981], p. 39) for ice at the melting point: n = 3 and A = 0.167 a⁻¹ bar^-3. These values are likewise assumed to be constant throughout the ice mass.

Two simulations, referred to as Case I and Case II, are shown in Figure 3. In Case I, the model is allowed to converge to the steady-state shape that is required for the given bed topography and balance rate (Fig. 1). This took about 25 years and produced some thinning in the upper half and some thickening in the lower half, except near the terminus where a total recession of about 200 m took place. (The term “steady-state” is used here in a practical sense: the glacier is considered to have reached an “equilibrium” shape when the calculated annual-thickness changes become of the same order as the error in measuring these data in the field, generally a few tens of centimeters in this case. Alternatively, a similar criterion based on surface speed could be used; since the speed is very sensitive to thickness, this takes a longer time, typically by a factor of two or so.)

In Case II, this steady-state shape is subjected to a uniform decrease in the balance rate of 0.25 m a⁻¹, and the model run to determine the new equilibrium shape. This took approximately 300 years.

Figure 4 shows a series of simulations, all similar to a Case II run, in which the glacier is subjected to different decreases in the balance rate. The approximate time taken to achieve these shapes ranges from about 300 to 1400 years. As expected, the overall change is greater as

approaches the peak of the balance curve (Fig. 1), at which point there would be no positive balance anywhere and the glacier would eventually vanish.

Fig. 4. The response of South Cascade Glacier to various uniform step decreases in the balance rate,

. Only the final steady-state shape for each case is shown. The values of

, in m a⁻¹, are shown on each line, and the approximate times it took each case to reach equilibrium are given in the table in the upper right. The dashed line shows the last time step obtained (t = 700 a) when the balance was decreased by –0.75 m a⁻¹

If the largest balance change shown in Figure 4 (−0.7 m a⁻¹) is decreased slightly, to −0.75 m a⁻¹, the model runs successfully for 700 years and reaches the configuration shown by the dashed line. After this, it breaks down as a section of ice near the terminus becomes stranded, a consequence of the local maximum in balance at this point (see Fig. 1). Note that this break-down occurred in roughly half the time that it took the –0.7 m a⁻¹ run to achieve about the same shape, demonstrating the rapidly accelerating response as

approaches the peak balance.

For the case where

= –0.1 m a⁻¹, the overall change in thickness at the terminus is of the order of 90 m, approximately three times larger than that obtained by Reference NyeNye (1963) for a zero-frequency perturbation. This is not surprising, considering that Nye’s theory only dealt with small perturbations, and, as can be seen from Figure 4, the perturbation near the terminus, even for the –0.1 m a⁻¹ case, is anything but “small”.

Figure 5 shows a different type of application of the model: using it to determine topographic effects on the ice flow. The bedrock sill in the middle of the ablation zone is the dominant feature of the bed topography of South Cascade Glacier. It is natural to ask what effect this has on the flow. Figure 5 shows the equilibrium shapes with and without the sill, for both the Case I and the Case II simulations. Qualitatively, the results are as expected. The greatest difference in thickness (between the “with” and “without” cases) occurs on the up-glacier side of the sill, where the longitudinal stress gradients would be expected to be affected most. Furthermore, the relative thickness change in this region is much greater for the Case II shape (about 35%) than for the Case 1 shape about (15%), since the sill blocks proportionately more of the overall ice thickness in the former. In both cases, the final length of the glacier is the same, independent of the topography, since the overall length of the glacier is determined by the condition that

when the glacier is in equilibrium with the balance rate.

Fig. 5. The effect of the bedrock sill on the flow of South Cascade Glacier. The end results for a Case I and a Case II simulation are shown.

Sensitivity of the Results to Calculation Parameters

One fundamental question which must always be asked of any numerical technique is how sensitive are the results to changes in the calculation parameters, such as the number of elements or the magnitude of the time step. When this is known, one can then proceed to study the sensitivity of the results to changes in the physical parameters, which, in this case, could be the flow-law constants, the balance rate, the density, and so on.

Figure 6 shows two such tests. In the first the number and distribution of elements is varied. If the number of elements is increased in the vertical direction, then, on the scale of the plot, any differences are completely insignificant. Increasing the number horizontally, however, does produce a noticeable change, but this is to be expected since smaller-scale variations in bed topography can now manifest themselves in the solution for the unknowns. Most of the differences, in fact, occur around the bedrock sill and the overall shape is still much the same over the remainder of the ice mass. This demonstrates an important result of finite-element ice-flow simulations using quadratic basis functions: two rows of elements are quite sufficient to model this complex non-linear flow.

Fig. 6. The sensitivity of the results to calculation parameters. The uppermost lines show the end result of a Case I simulation for various mesh sizes. No significant difference results when the number of elements is increased vertically (solid line), but slight changes occur when the number of elements is increased horizontally (short dashed line). N_ec is the number of element columns and N_er is the number of element rows. The lower set of lines shows the end result of a Case II simulation for various time steps, Δt. Values of 1, 5, and 10 years give the same results, but instability clearly results when a value of 25 years is used.

The second test illustrated in Figure 6 shows the effect of various time steps. The use of too large a step is readily apparent in the results; the solution starts to oscillate almost immediately.

A third test, not illustrated, was made to determine the effect of various values of the mesh-adjustment factor, δ. This will be tied to the number of element columns used, as well as to the complexity of the topography. For ten element columns, variations of up to 5–10 m in the final thickness occurred when δ was varied over a range from 0 (no adjustment) to 500%. For 20 element columns, the (figure 6) effect is much less and is only significant immediately down-glacier from the sill. Such test runs also reveal an optimum range for δ which appears to give more consistent results; for most of the simulations in this paper, this was between 100 and 200%.

Another test of the model was made by comparing with the analytic solution for laminar flow of a parallel-sided slab on a bed of uniform slope. The surface speed agreed to less than 1% difference, as found by Reference Hooke, Hooke, Raymond, Hotchkiss and GustafsonHooke and others (1979). However, since quadratic basis functions are used here, this agreement was, in contrast to their results, independent of the number of nodes distributed over the depth, and two element rows (five nodes) gave just as good agreement as seven element rows (15 nodes).

Use of a Small Computer

Numerical modeling has, up to now, been done exclusively with large “mini” or “mainframe” computers. The trend in computing, on the other hand, is towards small, desk-top machines. The computer code described by Reference Taylor and ZienkiewiczTaylor (1977) and Reference SikoniaSikonia (1982) was extensively modified to see whether finite-element ice-flow modeling could be successfully accomplished on such a computer.

The particular machine used here was a Digital Equipment Corporation PDP 11/23, with 256KB of memory. This is probably roughly equivalent to the newer and smaller IBM PC/AT. For the simulations involving a 10 × 2 element mesh, each solution iteration took approximately 1 min and the entire run about 1 h. Since no direct cost is involved, such times are perfectly acceptable.

Execution time is the limiting factor. Memory size is now sufficiently great on most computers that this becomes prohibitively long well before the memory of the computer is exceeded. The above runs, for example, used less than half the available memory.

Doubling the number of element rows results in approximately the same number of unknowns as doubling the number of element columns; however, the total execution time when the number of element columns is doubled is only about 60% of that when the number of rows is doubled. The reason for this is that increasing the number of unknowns vertically increases the band width of the coefficient matrix, whereas increasing the number horizontally does not. This is a favorable result, since, as shown earlier, it is the number of columns that usually needs to be increased.

Conclusions

A time-dependent finite-element model has been developed which successfully describes the flow of an arbitrarily shaped ice mass, where the ice thickness can approach zero at both ends. It is two-dimensional, partially temperature-dependent, and can accommodate sliding in several ways. It employs an automatic mesh-generating algorithm which optimizes the mesh for ice flow and greatly reduces the data-preparation time. It has been successfully implemented on a small computer and detailed analyses of ice flow can be made effectively within very reasonable time frames.

The use of quadratic basis functions, instead of the linear ones usually used, allows modeling to be done with only two element rows (five nodes vertically). More detailed resolution of the flow only requires increasing the number of nodes horizontally. This results in a considerable decrease in execution time, because not only are there fewer unknowns but also the equations can be solved faster when the nodes are distributed horizontally rather than vertically. Such basis functions permit a relatively coarse grid to be used and allow the full power of the finite-element method to be utilized. However, this is not without its pitfalls and care must be taken to ensure that the order of the basis functions is sufficiently “mixed” that numerical oscillations are not induced. Finally, interpolation in data sets must be done with a method such as cubic splines which does not produce discontinuous first derivatives; otherwise, unrealistic behaviour can result.

Extension of the model to three spatial dimensions would be straightforward, largely involving only input/output, since the basic structure of the program involves an arbitrary number of degrees of freedom. Execution times, however, would increase considerably and a fast mini- or mainframe computer would be required.

South Cascade Glacier is currently close to equilibrium. When the current balance is applied to the current geometry, the terminus recession and maximum change in thickness are only about 5%. In fact, even this 5% figure could be a result of using only a single arbitrary year for the balance data, rather than a long-term average over a number of years, which would be more appropriate. The relatively large changes right at the terminus might also be because the bed topography is not well known down-glacier from the bedrock sill; it has been approximated by a uniform slope but it is obvious from the initial surface topography (Fig. 1) that this is probably far from reality.

The bedrock sill, which has a height of about 20% of the ice thickness, produces a maximum increase in the overall ice thickness of only about 6 or 7%, compared to what the thickness would have been if the sill were not present. This maximum increase occurs along the up-glacier side of the sill. The effect of the sill, however, extends almost all the way to the head of the glacier.

South Cascade Glacier does not appear to be sliding along its bed very much. It was only by freezing the ice to the bed, and by lowering the shape factor to a minimally acceptable value of 0.5, that the calculated surface velocities could be reduced enough to bring them more or less into agreement with observed ones. (The region between the sill and the terminus is an exception to this and some sliding is possible here; on the other hand, the bedrock, as noted above, is so poorly known here that little should be inferred from this.) Why this should be the case is puzzling, since South Cascade Glacier is in a warm, wet maritime climate and flows down a reasonably steep slope. Possibly, drag exerted by the valley walls has a much greater effect than conventional shape-factor concepts imply; only full three-dimensional modeling can resolve this question. Alternatively, it could be that this glacier has a particularly efficient conduit system which channels water away rapidly and keeps most of the bed without much lubricating material, an inference which does agreed with the author’s frustrating attempts to measure sliding and basal water pressures on the same glacier (Reference HodgeHodge, 1979).