Reference RobinRobin (1955) was the first to propose the possibility that thermal instability could provide a trigger for the occurrence of surges. He argued on intuitive grounds that the strong stress-temperature coupling in the equations of motion could initiate a sufficient temperature rise in cold ice masses to melt the basal ice, and thus produce basal sliding. More recently, Reference ClarkeClarke and others (1977) and Reference Yuen and SchubertYuen and Schubert (1979) have proposed that ice sheets and glaciers may exhibit multiple steady states, and thus a runaway instability could occur if the viscous heating became sufficiently “large”.

Our purpose here is to draw attention to certain limitations in these authors’ models, which may render their conclusions inaccurate. If this is the case, then the question of the existence of thermal instability in large ice masses will require further study with more realistic models.

The model of Reference ClarkeClarke and others (1977) considers a temperature equation for a glacier of the form

where θ, τ, q,, ζ are dimensionless temperature, time, velocity, and vertical coordinate perpendicular to the slope; f(ζ, θ) is a non-linear function of the rheological properties of ice, and β is a dimensionless parameter which measures the importance of the viscous heating. The boundary conditions are

(prescribed temperature on the surface) and

(prescribed geothermal heat flux on the base).

Essential constituents of this model are that the depth is constant and prescribed (= 1, dimensionlessly), and the velocity field q is prescribed, and is independent of the lengthwise coordinate. In this case the authors obtain qualitatively similar results for zero and non-zero q, so it suffices to discuss here the case q = 0 (or alternatively q = (u(ζ), 0), θ = θ(ζ, τ)). The steady-state equation

with the boundary conditions (2) and (3), and an appropriate form for f, then admits multiple (up to three) solutions for various ranges of β. From this analysis, the possible occurrence of a creep instability is deduced if β becomes sufficiently large.

Now, this is only physically meaningful if the “control” variable β is determined purely from the input data to the system: alternatively, if we consider the glacier as a control system, then the control variable should depend only on the inputs, whereas the output variable (e.g. θ(1)) should be completely determined from the effect of the system on the input. Now the parameter β is in fact dependent on h, the glacier depth, which in a more complete model would itself be determined from the ice dynamics; thus there is a certain logical inconsistency in the model, and it is fair to question the applicability of the results.

A more general model, which incorporates the dependence of velocity and temperature on each other, was analysed by Reference Yuen and SchubertYuen and Schubert (1979), who used the depth h as the control variable and the surface velocity u ₀ as the output. They, too, obtained multiple solutions, but the physical meaning of these is similarly open to question. This is not to doubt the validity of the results, but simply to ask, given a set of physical inputs, how many solutions can there be? Clearly neither the depth nor the surface velocity is prescribed for a real glacier: in fact the only inputs are the bedrock slope and form, the surface temperature, the geothermal heat flux, and the accumulation and ablation rates. All other parameters occurring in any model are constant properties of the dynamics of ice. All the above are accounted for in the model of Yuen and Schubert, except the accumulation/ablation rate, and it is this which should therefore be used in determining a proper control variable.

In a steady state, conservation of mass requires that the flux through a vertical section at a point x measured from the head of the glacier be equal to the integral of the accumulation rate from 0 to x. Denoting this “flux function” by s(x), we see that s is a positive function, and is zero at the head and snout. If we consider the slab model as approximating locally the almost parallel-sided flow of a glacier, then it is clear that what should be prescribed in determining both h and u ₀ is the dimensionless flux

. Thus we should use s^*as a control variable, and not h: then, if multiple solutions were obtained, it would be reasonable to consider thermal instability as a possible surging mechanism.

The effect of doing this can be seen by examining Figure 3 of Yuen and Schubert’s paper. From their Figure 2, a reasonable approximation to s (∞s ^*) is obtained by putting

at least for high values of E^* , the activation energy. It is clear from Figure 3 that the lines log u₀ + log h = constant intersect the solution curves only once. This suggests that a control map of u ₀ versus s would be single-valued, and thus that real multiplicity of solutions for given external input functions would be precluded. In this case, the runaway instabilities envisaged by Yuen and Schubert will not occur.

We therefore suggest that a more realistic model of the shear flow occurring in glaciers and ice sheets would prescribe the flux rather than the depth, and that this may make the solutions single-valued. This is an effect of considering the top surface as a free boundary rather than a fixed one, and shows the stabilizing effect of consideration of such a free boundary in the problem. Since this appears to have such a dramatic effect, it seems reasonable to treat it in future models with as much priority as other “essential” features, for example advective terms, viscous heating, temperature-flow coupling, and so on.