1. Introductory motivation

There have been several attempts in the past few years to deduce an acceptable rationale for the description of the effects of long-time climatic variation of large ice masses. The need for a full thermo-mechanical treatment has already been emphasized by Reference Budd and RadokBudd and Radok (1971) but a full quantitative analysis for wholly cold, wholly temperate, or polythermal ice masses, which are partly grounded and partly afloat, is still beyond reach. Restrictions to simplified physical situations are needed or mathematical approximations must be developed that make the complex problem tractable.

Early developments were fraught with various ad hoc simplifications, often physically inferred, but equally often hiding the essential content that is being dismissed (see e.g. Reference Paterson and ColbeckPaterson (1980, Reference Paterson1981), Reference ThomasThomas (1979), and introductions to chapters 4, 5, and 6 in Reference HutterHutter ([^c1983])). What is needed is a presentation of consistent treatment of ice-sheet flow through a proper scale analysis of the fluid equations and boundary conditions, which permits estimation of the importance of the arising parameters and thus leads to systematic simplifications. Table I lists works that are based on such “rational” deductions and describes the physical situations to which they apply. Only non-surging flow conditions are presumed, and ice masses are either cold or isothermal throughout. The constitutive model in all these papers is that of a non-linear viscous heat-conducting fluid that obeys a creep law which is akin to Glen’s flow law. Further, the temperature-dependence of the stress tensor (when present) is included through a temperature-dependent rate factor in the deviatoric stress–strain-rate relationship. While it is probably difficult to judge priorities, it is fair to say that, of the authors listed in Table I, Fowler and Larson pioneered the scaling analysis for glaciers, whereas Morland was the first to present the coupled thermo-mechanical scale analysis.

table 1a.isothermal conditions

table 1b.non-isothermal, cold conditions.

Some of the reduced models of Table I apply to time-dependent and spatially three-dimensional situations, but computations have been almost exclusively restricted to steady, plane, or axisymmetric flows and homothermal conditions. The corresponding numerical analysis (Morland and associates) showed that the geometry of an ice sheet and the velocity distribution within it depend to comparable magnitudes on both the sliding motion and the internal viscous deformation. Moreover, the thermo-mechanical analysis (Reference MorlandMorland, 1984) corroborated the conjecture that temperature variations are equally significant. However, the coupled thermo-mechanical steady-plane problem of simultaneous determination of ice-sheet geometry, velocity, and temperature distributions offers considerable numerical difficulties which must be overcome before a full three-dimensional steady or time-varying analysis, or a physically more complex situation (for instance, polythermal conditions) can be attacked.

In this paper we present the computational procedure that we used to determine the geometry, the velocity, and temperature distributions in plane grounded ice sheets under steady conditions. We briefly state the governing equations, explain the computer code “THEMPIF”, but refrain from presenting numerical peculiarities and a convergence analysis (see Reference Yakowitz, Yakowitz, Hutter and SzidarovszkyYakowitz and others, in press). Instead, we aim at a deeper physical understanding and thus explore the roles of the scaling parameters, the external forces (accumulation, surface temperature, geothermal heat), and the creep law.

We shall not dwell on applications of concrete physical situations because experience has shown that we and the model are not yet ready to predict reliably actual glacier situations. We would rather like to demonstrate here how interactive computer techniques with some sophisticated graphics hard- and software can be constructively used, first to demonstrate to the glaciologist how useful a proposed flow model is and, secondly, to explore it in the range of its suitability. In this particular case, we uncover a rather restrictive limitation of what might be called the Morland–Hutter model, and find in this way hints as to its improvement. Clearly, the earlier and less rigorous models are fraught with similar, or even greater, restrictions than the Morland–Hutter model. In fact, it was the consistency of the scaling and perturbation procedure (that led to the leading order model) which laid the fundamentals for this discovery. We also demonstrate that in the range of its applicability the strong thermo-mechanical coupling suggested by the Morland–Hutter model is rather weak and probably atypical for many glaciological applications. Suggestions for improvements will be given in the last section.

2. The continuum model and its simplifications

Governing equations of the rheologically non-linear incompressible viscous heat-conducting fluid, upon which our cold ice-sheet flow model is based are the

and the

These relations have been amply explained in Reference HutterHutter ([^c1983], chapter 1) and Reference MorlandMorland (1984).

Boundary conditions that must be invoked are: at the free surface

and at the fixed base:

We ignore effects of basal drainage, stress-induced recrystallization, and changes in constitutive relations due to changes in contents of impurities, gravel, etc. (Reference Hutter and VullietHutter and Vulliet, 1985).

The union of the above-mentioned equations forms an initial boundary-value problem in a space–time domain of which the determination of the extent is part of the solution of the problem.

Non-dimensionalizing the equations that describe a particular physical problem is advantageous as it discloses the physically important parameters through the introduction of typical scales. Because the latter are to a certain extent arbitrary, the form of the resulting equations may differ slightly from one author to another. Often, differences are only semantic and restricted to notation and procedure. For instance, Reference HutterHutter’s ([^c1983], chapters 3 and 5) and Reference MorlandMorland’s (1984) glacier and ice-sheet analyses are the same as a one-to-one correspondence between the two sets of equations can be established (see Reference HutterHutter, [^c1983], chapter 5). Here we use Reference MorlandMorland’s (1984) set of equations.

Scales are introduced for stresses (σ₀, ρgd ₀), velocities (q _m), a typical depth (d ₀), a characteristic length (d ₀/∊), a typical strain-rate (D ₀), a temperature range ΔT, etc. (Table II). For instance, σ₀ is a typical (shear) stress range to which the ice under natural conditions may be subjected. In terms of these quantities and the physical parameters of the model (Table III), the variables in physical space (denoted by asterisks) and the corresponding dimensionless variables are related by

If the scales are appropriately chosen, the dimensionless variables X, Y, V, W, and T on the right-hand side, when determined by the field equations, are of order unity. A sketch of an ice sheet on a gently sloping base is shown in Figure 1.

Fig. 1.

Geometry of a plane ice sheet on a base. Cartesian coordinates (X, Z) are in the direction of a mean basal plane sloping at small angle γ. Basal surface is defined as Z = F(X) and steady-state free surface as Z = H(X). Upper and lower margins (or left and right margins when γ = 0) are defined as the intersections H(X) = F(X) and the ice divide is where dH/dX = γ.

table II.scales used to non-dimensionalize the field equations and boundary conditions. numerical values are typical for one situation. a range of values is discussed in the text

table III.physical parameters

In dimensionless form, the field equations and boundary conditions can be explicitly expressed in terms of the aspect ratio ∊ which is small (see below). On length scales which are larger than the ice-sheet thickness d ₀, one may therefore restrict attention to the limit equations as ∊ becomes vanishingly small. This limit is called the shallow-ice approximation; its steady-state field equations are (Reference MorlandMorland, 1984)

in which p is the pressure, σ _xz the shear stresses, J the second-stress deviator invariant, and Z = H(X) the equation of the surface profile (see Fig. 1). Moreover, s, δ, α, β, v, and θ are the dimensionless quantities

Four of these are independent. We call s and δ mechanical parameters, α the thermal dissipation number, and β the diffusion number. Finally, a(T) is a temperature-dependent rate factor, and ω(J) and g(τ, θ) are creep-response functions of which the explicit form depends on thes creep law that is used. Reference Smith and MorlandSmith and Morland (1981) based their expressions on Reference GlenGlen’s (1955) and Reference Mellor and TestaMellor and Testa’s (1969) uniaxial compression data and obtained

with excellent data fitting in the stress range 0 → 5 × 10⁵ N m^–2 and the temperature range 213 → 273 K. If the constitutive relation in physical space has the form D* = A(T*)Ω(σ_II*)σ*, then

Often A(T*) obeys an Arrhenius relationship and is a power law.

Before we proceed, a few clarifying remarks are in order. For one, the rheological properties stated in Equations (4) are one of many possible forms. Our code allows whatever relationship the user prefers. Moreover, the second stress-deviator invariant in Equation (2e) is given by the shear stress squared, which is inaccurate close to the divide. Difficulties are therefore expected to arise, but a model that handles ice-divide behaviour more appropriately is far more complicated than this one. We regard it therefore as advantageous to collect experiences with this simplified model.

In the shallow-ice approximation, the boundary conditions become, at the free surface:

and at the base:

in which Z = F(X) is the equation of the rigid basal surface and θ_z (= 1 for values of Table II) is the dimensionless geothermal temperature gradient. The functions on the right-hand side of Equations (6) and (7) are, in general, order-unity functions which describe the effect of the climate[a_cc (·) and T _s(·)], and the conditions of the substratum [Q (·), U _F(·)] on the ice sheet.

These functions may depend on a number of additional variables and will be specified in detail below. For instance, the sliding velocity U _F depends on the shear traction, the local pressure, and perhaps the basal temperature; the geothermal temperature gradient may vary with position. Alternatively, the accumulation rate is determined by the snow-equilibrium height H _e and the position on the surface, and the free-surface temperature may depend on a mean temperature at the grounding line, T _M, and on position, in short

Reference Morland, Morland, Smith and BoultonMorland and others (1984) used a sliding law in which the tangential traction t _s and the sliding velocity u _s are linearly related, with a coefficient which depends on the normal traction t _n : t _s = t _n μ(–t _n)u _s. Explicit dependence on t _n is required for finite-margin slopes, if μ(·) is a non-vanishing and finite function of the pressure. In terms of the dimensionless variables and in the shallow-ice approximation, the sliding law proposed by Reference Morland, Morland, Smith and BoultonMorland and others (1984) becomes

where is a scale for μ and an order-unity function. Details on the specific choices of the functions listed in Equations (8) and (9) will be given later.

As a matter of completeness and for later use, we mention that the continuity equation and the boundary conditions. Equations (6a) and (7b) imply the integrated mass-balance statement

which will be used in the ice-divide analysis.

Evidently, the general ice-sheet flow problem is described by three classes of dependences:

• (i) the scales s, δ, α, β (see Equations (3)),

• (ii) the constitutive properties of our ice model through the rate factor and the creep-response functions (see Equations (4)),

• (iii) the external forcing through the accumulation-rate function a _cc (·), surface temperature T _s(·), the sliding law U _F(·), and the geothermal heat θ_z Q(·) (see Equations (8)).

It is our intention to investigate the role of these in greater detail. Ensuing developments will therefore mostly be restricted to situations for which X is horizontal and Z vertical (y = 0).

1. Description of the solution procedure

The idea is the recognition that Equations (2) can be solved by marching forward in the X-direction, solving for each X a set of (ordinary) differential equations in Z. Marching must be in the direction of U, for otherwise the problem is ill-posed. In other words, integration must start at the ice divide and always proceed “down-glacier”. This is a limitation as the location of the divide X = X _D and its depth H _D must be known or estimated and iteratively improved. So far we have not handled geometries with multiple domes; but symmetric situations are easy as the location of the divide is usually at the symmetry line.

Assume for the moment that sufficient conditions at the ice divide can be prescribed which permit integration away from the ice divide. We may then assume that the functions listed in Equations (2), (4), and (8), and the fields V, W and T, etc. are known at a particular column, say X = L (Fig. 2). Let these be V ₀, W ₀, T ₀, etc. At the neighbouring column X + ∆X = R, these same fields can be evaluated from FD-representations in X of the last three Equations (2f–h), viz.

These are linear differential equations for U, W, and T in the variable Z at column R, which are subject to the boundary conditions, Equations (6b) and (7b, c). They can be solved by the standard two-point-boundary value problem (TPBVP) routines (Reference Szidarovszky and YakowitzSzidarovszky and Yakowitz, 1978). This procedure requires that H at column R is known.

Fig. 2.

Right portion of a symmetric plane ice sheet with γ = 0. The X-coordinate is horizontal, the Z-coordinate vertical, and the flow is from left to right. The ice divide with heights H_D is at X = X_D = 0 and the right margin is at X_R. So ice-sheet semi-length is L = X – X_R = X_R. Two columns at X (denoted by L) and X + ∆X (denoted by R) indicate the finite-difference approximation used in the X- and Z-directions.

Once these equations are solved for U, W, T, etc., the boundary condition, Equation (6a) may be used to evaluate dH/dX; the new height at the column to the right of column R may thus be computed. This column will now become the new column R, whereas the old column R will become the new column L. In this fashion, one can proceed until H = 0 is obtained, which marks the outer margin X = X _R. Finally, the total accumulation can be computed,

which, because of steady mass balance, must vanish. However, a _cc ^int ‡ 0, in general; so ice-divide depths H _D must be varied until a _cc ^int(H) = 0, to a sufficient degree of accuracy.

2. Construction of the solution at the ice divide

The ice divide is defined as the location X = X _D at which H’ = dH/dX = γ. Assume that the height at the divide is prescribed. Then, from Equations (12), (4), and (9), it follows that U = τ = σ _xz = g = 0. Thus, from Equations (2f) and (9), we may deduce by integration with respect to Z and subsequent differentiation with respect to X that

in which H”(X _D) = d²H/dX ²(X) is unknown and ω(0) = ω(J = 0) (cf. Equations (4)). Similarly, by using a Taylor series expansion in X of the mass-balance statement in Equations (10), we obtain

With the aid of Equation (13b), this yields an expression for H”(X _D), namely

Evidently, the curvature of the free surface at the ice divide, H”(X _D), is determined by the ratio of the value of the accumulation-rate function and a “normalized flux” which consists of a gliding and a sliding contribution. In principle, Equation (15) applies to the no-slip situation but a careful ice-divide analysis shows that longitudinal stretching is significant under such circumstances; this is a case not treated here.

The same inferences can also be drawn from the numerical computations which follow. Thus, we assume viscous sliding. With the aid of Equation (13), Equations (2g, h) may be written as

Thus, introducing the auxiliary variables V = ∂W/∂z, S = ∂T/∂z, the following TPBVP for W and T, valid at X = X _D, may be deduced

subject to the boundary conditions

which can be inferred from Equations (6), (7), and (2g). The TPBVP in Equations (17) and (18) can be solved by using the shooting method or quasi-linearization. Return now to Equations (11), in which quantities bearing the index zero (at the left column = ice divide in this case) are presumed to be known. Once the ice-divide analysis is performed, this is indeed so; thus Equations (11) can be solved tor U, W, and T, etc. in the next column R, and H’, H can be computed in the column to the right of it, and so on, as explained in the text following Equations (11).

3. A brief description of the code

THEMPIF is a FORTRAN code which permits numerical computation of THErmo-Mechanical balances of Plane Ice-sheet Flows. It consists of a secant root-finder SECA which searches for H _D such that the equation a_cc ^int (H _D) = 0 is satisfied. The main body of the program consists of two sets of sub-routines (see Fig. 3). The first set provides the computational procedures for the integration of the differential equations and the evaluation of the integrated accumulation rate a_cc ^int(H _D); the other set serves to define the physical conditions and is kept flexible to adjust the physics to the situation at hand.

Fig. 3.

Chart of the structure of the computational routine THEMPIF.

Sub-routine INPUT is designed for the user to prescribe the scales and physical parameters listed in Tables II and III. While the physical parameters are well defined, the glaciologist must choose the scales best suited to his example. Sub-routines RATE, OMEG, and CREEP define and compute the rate factor and the creep-response functions of Equations (4). Here is the place to vary stress-deviator–strain-rate relationships. Sub-routine TOP calculates the dimensionless accumulation rate and the surface temperature; their functional relationships must be prescribed by the user. The geometry of the base, the geothermal heat flow, and the drainage are defined in sub-routine BOTTOM; however, the sliding conditions Equations (7b, c) with the specifications, Equation (9), are accordingly accounted for by sub-routine SLIDE. In its execution it makes use of calls to BOTTOM and SHEAR for the evaluation of basal geometry and basal stresses.

The computational routines constitute a sub-routine FF, which simply computes the value of a _cc ^int(H _D); however, this requires integration of Equations (2), subject to the boundary conditions, Equations (6) and (7), and is achieved by calls to the sub-routines DIV, SHEAR, COEF, and possibly COMPARE. Sub-routine DIV constructs the solution for U, W, T, and H” at the ice divide (given its location and height) and thus provides the initial conditions for marching away from it. SHEAR permits evaluation of the shear stress, the second stress-deviator invariant, etc. listed in Equations (2a–d), and COEF computes the coefficients of the differential equations arising on the right-hand side of Equations (11). The main routine FF integrates these equations column by column, as described earlier, stores the velocity and temperature profiles, and the values of H and the integrated accumulation between the divide and the actual position of the column until H vanishes. In the case that the non-linear variant of Equations (11) is used and quasi-linearization is employed, sub-routine COMPARE decides whether quasi-linearization iterates have converged or not.Footnote *

A detailed description of analysis of the computational algorithm and its stability and convergence conditions has been given by Reference Yakowitz, Yakowitz, Hutter and SzidarovszkyYakowitz and others (in press). A user’s manual has also been produced (Reference Hutter, Hutter, Yakowitz and SzidarovszkyHutter and others, unpublished). The program is designed to use interactive computing to advantage. The user can experimentally prescribe values of ice-sheet parameters of interest and solve the associated problem without re-compiling the program. We have found that our graphics terminal enhances these experimental studies by actually plotting the profiles and the velocity and temperature fields, as opposed to presenting tables of numbers. Through readily understood pictorial displays of solutions, we are developing our intuition regarding influence and sensitivity of ice-sheet response to variations of model parameters and boundary conditions.

1. Selection of forcing functions

Standard calculations are based on the following explicit expressions for the forcing functions: the accumulation pattern is based on an extension and an alteration of Reference Morland and SmithMorland and Smith (1984).

Footnote ^† and

in which a and b are any combination of the sets a = (12.5, 6.25, 3.125) and b = (0.5, 0.75, 1). Moreover,

in which a _Ref, , p ₁, and p ₂ are constant parameters. Evidently, a _cc varies with height through ã_cc(·) and with horizontal position through the function of snow-equilibrium height H _e. The constant a _Ref allows for absolute changes in accumulation while keeping the scale fixed. Typical values are, perhaps

With p ₁ = p ₂ = 0 and the choice of Equation (20a) the accumulation-rate function has a continuous derivative throughout, is constant and positive for H > H _e, decreases according to a cubic law when H _e + 0.25 003E; H > H _e, and increases linearly with height when H < H _e. When p ₁ and/or p ₂ are non-zero, ã _ccis still continuous with a continuous derivative everywhere except at H = H _e + 0.25. Values of ã_Cc computed according to Equation (20a) are between –6.25 and 0.5, so a_cc = 0(a _Rrf). Equation (20c) is a simpler, continuous variant of Equation (20a) which allows us to estimate the sensitivity of the ice sheet to variations in accumulation distributions. Also, the variation (21) has been introduced, because computations indicated that a constant accumulation rate close to the ice divide resulted in very long (perhaps infinitely long) ice sheets. In fact, margin positions were never found when . Furthermore, Equation (20) for H < H _e implies relatively rapid changes of ã _?? with H(X), which implies caution in the numerical evaluation of the integral a _cc ^int(H _D) in Equation (12). A simple Riemann sum may require a smaller-mesh ∆X whenever H < H _e.

Computations were also performed with the underlined terms in Equation (20a) being omitted; ã _cc was then discontinuous at H = H _e + 0.25, with a sudden increase as one moves down-slope, permitting (at least in a gross fashion) a description of sudden, spatial climatic changes. Equations referring to this case will be numbered as Equation (20b).

The question why horizontal equilibrium heights caused difficulties will be addressed later.

For the surface temperature, we take the linear representation

with –2.0 ≼ T _M ≼ 0, corresponding to –40° C < T* _M < 0°C, and 0 < T _M ⁽¹⁾ < 0.1, as well as 0 < T _M ⁽²⁾< 0.1. When T _M ⁽²⁾ = 0, then T _M is the value of the surface temperature at H = 0, agreeing with the margin temperature when F(X) = 0.

Standard computations are also performed with a flat base and for a constant geothermal heat flow corresponding to

with θ_Z = 1 for the values of Table II and a range of realistic dimensionless geothermal temperature gradients 0 < θ_Z < 10. Moreover, the sliding law deduced by Reference Morland, Morland, Smith and BoultonMorland and others (1984) for their isothermal computations of the Greenland ice sheet is adopted, with

in which x = cos γ. (H – F), μ _Ref = 2.5 × 10⁻⁵ and

In the domain 0 ≼ x ≼ 1.6, the range of is approximately . Computations show that the value of μ _Ref is critical.

It is clear that there is no universal sliding law, and the choice of a basal sliding condition will not only depend on the particular condition but also on the discretization and the amount of basal sub-grid shearing. We use Equations (25)–(27) as reasonable examples found in the literature that have proved useful for the Greenland data but keep the program flexible enough to account for any other relationship.

2. Preliminary computations. The limitation of the model

Initial computations with the above forcing functions proved that a finite-length ice sheet could numerically only be obtained with p ₁ > 0 and/or p ₂ > 0 in Equation (21). A separate analysis of the full equations close to the ice divide has further shown that an accumulation function without an explicit X-dependence, a _cc ≠ a _cc(X, H(X)), implies H”(X _D) = 0, corresponding to a completely flat divide region. Moreover, it transpired that satisfaction of the surface mass-balance condition, Equation (6a), immediately off the divide was possible only with a sufficient amount of sliding. We regard it as a significant “computer-aided finding” that under the present Morland-Hutter reduced ice-sheet model, the magnitudes of the accumulation functions and sliding functions must be related. Having seen it computationally, we are able to justify this assertion analytically. The model is scaled in s and δ so that the ice-sheet height is approximately 1 and its semi-length is approximately 1 or more. As H”(X) monotonically decreases for most accumulation functions, one would presume its magnitude at the divide is far less than 1. A value of about 0.1 or less would be more appealing to experience and intuition. This being the case, one is forced to conclude, in view of Equation (15), that sliding has to be far larger than gliding. For if a _cc(X _D , H _D) is about 0.5 and H”(X _D) is about 0.1, then the denominator must be about 10. But the gliding or shearing term obeys the inequality

as a(T) is less than 1.0 for temperatures below freezing. (In fact, we would anticipate a(T) to be far smaller than 1.0 down most of the divide.) In view of these numbers, from Equation (15) it follows that

or

Thus, in order that surface mass balance can be satisfied at the divide, we necessarily need

In a situation for which Equation (31) is not satisfied, our reduced equations are inappropriate, at least close to the divide. In such an event, longitudinal stretching must play a non-negligible role close to the divide. A divide analysis of the full theory is necessary and results obtained with it must be matched with those of the reduced model and valid in a region distant from the divide. More on this will be said in the conclusions section.

Reference Morland and SmithMorland and Smith (1984) used μ _Ref = 2.5 × 10^–5, ∊² = 2.75 × 10^–6, explaining why computations of results obtained with their parameters led to violations of surface mass balance. On the other hand, μ _Ref = 2.5 × 10^–7, ∊² = 2.75 × 10^–6 led to acceptable results.

Physically, Equation (31) sets a minimum amount to the “lubrication” of the bed, which according to Equations (3) is basically controlled by the scales for accumulation, glacier thickness, and stretching. The larger the typical accumulation, and the smaller the ice-sheet thickness and the typical scale for the stretching are, the larger will be the value of ∊², and consequently the less can be the lubrication of the bed. The no-slip condition can however never be accepted; this is obviously a disadvantage.

Initial computations were based on the physical parameters of Tables II and III and the forcing functions of the last paragraph with the accumulation function, Equations (2a, b) and the specifications as shown in Table IV.

table iv.values of the forcing function parameters for two sets of case studies, in which surface conditions are varied

Figure 4 shows results for the Case I study, in which accumulation conditions are varied. In Figure 4a we display the ice-divide height H _D as a function of p ₁ which is the slope of the snow-equilibrium height, parameterized for a range of values of equilibrium height H ⁰ _e. Figure 4b shows the corresponding graph for the half-width length X _R of the ice sheet. Evidently, H _D and X _R depend conspicuously on the height of the snow-equilibrium line and to a lesser extent also on p ₁. Nevertheless, the dependence on p ₁ is surprisingly strong, if physical quantities are looked at (see Table V). A change of the snow-equilibrium height within a distance of 1200 km by only 20 m (= 2%) changes the ice-divide height by 260 m (corresponding to 10.7%). Consequently, ice-divide height and ice-sheet extent depend critically on the accumulation function and small changes of it.

Fig. 4.

Ice-divide height H_D (left) and ice-sheet semi-length X_R (right) plotted against p₁ for various values of the snow-equilibrium height (see Equation (21)) for the Case I study (cf. Table IV). Dots are computed values: solid lines connect them smoothly. Dashed lines are extrapolations into ranges of the parameter p₁ where iterative computations failed or were not pursued, or required some interpolation. Ice-divide heights and extents depend critically on both and p₁.

table v.

ice-sheet semi-lengths and divide heights, and their absolute and relative changes with the accumulation parameter p₁, corresponding to values and scales and physical parameters as listed in tables ii and III

Keeping the accumulation conditions fixed by varying the surface-temperature distribution (margin temperatures are varied from –20 °C to –0.25 °C in Case II) led to the result that ice-divide heights and ice-sheet semi-lengths were unaltered to two significant figures (the only reliable ones anyhow). This remained so when the geothermal heat flow was varied by several orders of magnitude (0 < θ_Z < 5). The reason for this is obviously the fact that with μ _Ref = 4.15 × 10^–8 sliding over-rides gliding by two orders of magnitude, so that temperature can only marginally affect the geometry and flow pattern.

An anonymous referee was kind enough to point out that a horizontal equilibrium line H _e (no X-dependence) corresponds to an unstable ice sheet. He states that “[for p ₁ = 0] … any perturbation from a steady-state solution should have led to a solution of an ice sheet that either blew up because it became so large or shrank to nothing. What factors are in the computer code that take this into account ….” At least, in parts, this statement is misleading. To answer the question, we remark that with an accumulation function that is completely flat close to the divide, it can rigorously be proven that a marching procedure with O(∆X)-accuracy can never be initiated at the divide. In other words, our numerical scheme fails for p ₁ → 0. Extrapolated curves in Figure 4 must therefore be taken with care. The case p ₁ = 0 should still be computable but a totally different numerical integration technique is required.

With regard to stability, we can only state that this analysis cannot answer it. A small perturbation analysis about a steady state that is already a known solution could yield a statement in favour of or against stability. Inferences from other, simpler, models such as those of Weertman or Oerlemans are dangerous and most likely wrong, because the models are different in more than trivial details.

Figure 5 summarizes the results of a typical run, for conditions as described in Table IV and the figure legend. In this figure the top two graphs (a, b) display the temperature distribution in the form of isotherms and vertical profiles, respectively. They show the pattern one would expect, given the available data from observations and earlier approximate models (see, for instance, Reference JonesJones (1978)). The aspect ratio of the plotted ice sheet differs from that obtained from computations because we have chosen to scale the horizontal and vertical coordinates such that the screen of our graphics terminal would show the results optimally. Thus, the aspect ratio of the graphs is constant (= 0.54) but the true aspect ratio can easily be inferred from the coordinate scales shown on the graphs (compare legend to the figure).

Fig. 5.

Distributions of temperature and (scaled) velocities for the Case I computations (cf. Table IV) using , p₁ = 0.1, T_M = –0.1, T_M = ⁽¹⁾0.8, T^e _M ⁽²⁾ = 10⁻². The temperature distribution, indicated in °C, is shown for the top two graphs, (a) displaying a few selected isotherms, and (b) showing vertical profiles, with the temperature scale drawn as an inset.

The plots (c), (d), (e), and (f) depict the non-dimensional and scaled velocities. The first figure (c) shows vertical profiles of the (total) horizontal velocity; figure (d) shows the same profiles for the difference (U – U_F), called gliding velocity, while figure (e) gives vertical profiles of the vertical velocity W. Dimensionless scales for all three are given as insets. Finally, figure (f) gives the vector plot for the velocities, indicating the flow pattern within the ice sheet. Scales are not indicated for the reasons explained below.

In these graphs (and alt similar ones which follow) the aspect ratio ( = 0.54) of the plotted ice sheet is always the same and usually differs from the computed aspect ratio . The value of can, however, be inferred from the graph as horizontal and vertical scales for X and Z are indicated by the scales along the coordinate axes.

Also, the lengths of the vectors in the plot of figure (f) are not proportional to but proportional to , where . In this particular case, so λ = 0.54/0.216 = 2.5. Thus figure (f) gives the correct stream-line pattern on the basis of the aspect ratio . This procedure was adopted because it provided the optimal use of the screen of our graphics terminal and led to the standardized form of figures like this one. Notice, finally, that positive values of temperature and velocity profiles are drawn to the right, while negative values are drawn to the left.

Figure 5c, d, e, and f summarize the results obtained for the dimensionless velocity distribution. Graph (c) shows vertical profiles for the total longitudinal velocity U, graph (d) the difference between U and the sliding velocity U _F, characterizing the flow component due to viscous deformation. This difference velocity will be called gliding velocity. In view of the scales shown as inserts on these graphs, we see that the gliding velocity is approximately 0.5% or less of the sliding velocity, which corroborates the statement made about the role of sliding. Note the continuous growth of the gliding velocity as one moves upwards away from the bed. This behaviour is, of course, expected but is nevertheless surprising, because it means that at least as far as horizontal velocities are concerned, our program guarantees more than two places of accuracy, for otherwise inconsistent results for the gliding velocity would be expected. Ensuing computations will have to explore situations in which gliding plays a more important role.

Figure 5e displays vertical profiles of the dimensionless vertical velocity. The linearity of the profiles has often been conjectured in glaciology, and was first used by Reference RobinRobin (1955) to explain the contribution of vertical convection to the temperate distribution. Here, it is a proven result of the computation. Its accuracy is nevertheless somewhat surprising. Notice also that W is downward everywhere including the ablation zone, contrary to what one might expect. The reason is that in the ablation area |dH/dX| and U are comparatively large, so that in the surface boundary condition, Equation (6a), the product UdH/dX < 0 will outweigh (–a _cc) > 0.

The flow pattern along the free surface is still as expected, namely into the ice within the accumulation area and out of the ice in the ablation zone. This is demonstrated in Figure 5f, which is a vector plot of the velocity distribution (see legend for interpretation) that gives a fairly reliable view of the stream-line pattern.

In this description we have restricted attention to the scaled dimensionless variables, the only exception being temperature. Physical variables can easily be deduced with the aid of Equations (1) and the characteristic parameters shown in Equations (3) and Table II. In this particular case multipliers are as follows:

This yields an ice-sheet extent of the order of 4000 km and gigantic horizontal velocities of 200–300 km a^–1. The next section will indicate why this somewhat unrealistic case is still meaningful.

3. Variation of scaling parameters

Consider Equations (3) and substitute appropriate values for the scales; then it is found that realistic ranges for s, δ α, and β are, approximately

so that v(= ∊ ² ) and θ assume values as shown in Table VI. Evidently, the ∊² matrix is symmetric with respect to the main diagonal, whereas the θ matrix enjoys this symmetry with respect to the second diagonal. Large ∊² values lie at the right lower corner of the ∊² matrix, large θ values at the left lower corner of the θ matrix.

table vi.values of ∊² (left matrix) and θ (right matrix) for ranges of values for s and δ. the ∊² matrix is symmetric with respect to the main diagonal; the θ matrix has symmetry with respect to the second diagonal. also shown are the formulae for ∊² and θ in terms of the basic scale

To explore the significance of the parameters in the ranges mentioned in relations (32), two sets of computations were performed:

For Case III, 76 runs were performed in the intervals

and using the accumulation function, Equations (20a) and (21), with p ₁ = 0.1, p ₂ = 0, , with the remaining boundary parameters as listed in the Case I study (Table IV). In so doing, it was strictly observed that μ _Ref , and ∊² satisfied inequality (31).

It was found that dimensionless ice-divide heights H _D and dimensionless ice-sheet semi-lengths X _R did not depend in any recognizable way on all three parameters s, δ, and μ _Ref but only on one single quantity, namely the π-product

Only when the second parameter, θ = δs^{– 1}, was large (θ > 100 or θ > 1000, in general) could a slight trend of dependence of H _D and X _R on θ be observed but it was inconclusive, and often computations resulted in overflow due to strong rheological non-linearities and/or high temperatures (above melting). Table VII documents this behaviour for a few runs; Figure 6 illustrates the dependence graphically. “Error bars” mark the range of values for X _R and H _D, respectively, which were obtained with different values of the parameters s, δ, and μ _Ref. The width of the error bars is more the reflection of the accuracy of our numerical scheme than an indication of a dependence on other variables (e.g. θ). We also observed that with μ _RAF/∊² > 0.1 or with the secant root-finder for a _cc ^int(H _D) sometimes failed to converge unless initial guesses for H _D were close to the solution value.

Fig. 6.

Graphs of the ice-divide height H_D, the ice-sheet semi-length X_R, and the aspect ratio plotted against μ_Ref/∊² in semi-logarithmic representation for Case III computations. The scale for H_D is shown on the left-hand side; those for X_R and are given on the right-hand side. “Error bars” indicate the range of values of the respective variable that were obtained by using different values of s, δ, and μ_Ref as explained in the text. Solid lines are eye-fitted interpolations. The aspect ratio has been computed by using the interpolated curves.

table vii.some values of the characteristic parameters s, δ, μ _Ref, θ, and corresponding computed non-dimensional ice-divide heights H _D and ice-sheet semi-lengths for the case iii study. computed values are reliable only to two significant figures

Figure 6 and Table VII further suggest that we may set

in which the asterisk refers to physical variables and the dot in the argument indicates functional dependences that were omitted in this run. Thus, the aspect ratio is a function of the form

The function is also plotted in Figure 6. Accordingly, basically increases for increasing values of μ _Ret/∊². Since for ∊ as fixed, decreasing μ _Ref.means enhanced sliding in comparison to gliding, and Figure 6 implies the following: leaving all other variables fixed, the more slippery a bed is, the shallower the corresponding ice sheet will be. Of course, this is in agreement with intuition. Note, finally, that the aspect ratio in physical space, , is a function of two parameters, e and μ _Rrf/∊². The shallowness of an ice sheet now depends on e and μ _Ref/∊² but it is still true that increasing basal friction reduces the aspect ratio .

In the computations for Case IV, the thermal parameters α and β were varied such that 10^–2 < (α, β) < 1, while μ _Ref/∊² was kept constant (= 0.05) and s < 1, δ < 1 were arbitrarily varied such that δs^–1 < 10³. The thermal-boundary conditions were those of Case I (see Table IV) and the accumulation function was that of Equations (20a, b). The result: H _D and X _R were essentially unaltered by the variations of α and β (20 runs). This feature was also corroborated in a few additional runs using different μ _Ref/∊² ration and higher surface temperatures.

It thus appears that in the range of the validity of this model, μ _Ref/∊² is the only dimensionless parameter of the field equations which affects the dimensionless ice-divide height H _D, the ice-sheet semi-length X _R, and the aspect ratio . Because comparison of surface profiles computed for various parameter sets but with the same aspect ratio has not led to noticeable differences of profile geometries, these results imply that temperature does not affect the geometry. The latter could have been computed by ignoring temperature variations. This result is surprising and a consequence of Equation (31).

That the geometry remains essentially unaltered, when α and β are varied, does not imply that the flow and temperature distributions would not depend on these parameters. In fact, intuitively, we might anticipate that the dissipation number will hardly affect the flow and the temperature distribution unless the temperature in the entire ice sheet is close to melting. The temperature distribution is then primarily governed by convection and/or diffusion. Furthermore, for small diffusion numbers β, the temperature distribution must be dominated by convection except in a basal boundary layer. Larger diffusion numbers will (most likely) not give rise to such boundary layers. For Case IV runs with all combinations of α = (10^–2, 10^–1, 1), β = (10^–2, 10^–1, 1), except α = 1, β = 1, it was found that the dimensionless velocity fields hardly differed from each other. Plots for these are therefore only shown for α = 10^–1, β = 10^–2 (compare the lower part of Figure 7). The temperature distribution is, however, mainly governed by the relative weights of diffusion versus vertical convection (compare graphs (a) and (b) in Figure 7 with graphs (a) and (b) and (c) and (d), respectively, in Figure 8). Evidently, when β is small, vertical convection dominates; vertical temperature profiles change slowly in the upper part of the ice sheet but relatively quickly close to the base. The boundary layer can clearly be seen in the isotherm plot (Fig. 7a). Isotherms have the typical shape known from other studies. One detail in these temperature distributions should be emphasized: over most of the ice sheet the temperature profile for fixed X shows an inversion; in other words, along a vertical line, the temperature is coldest not at the surface but at a certain depth. The prediction of this feature has a relatively long history and in steady state has been demonstrated to be due to longitudinal advection (Reference RobinRobin, 1955; see also Reference HutterHutter, [^c1983], p. 171–73). The location of the inversion point relative to the surface varies with position (it is close to the surface towards the snout). Its existence is to a large extent the result of the fact that thermal diffusivities are small. Figure 8 corroborates this statement. When β = 0.1 (top of Figure 8), vertical temperature profiles are still curved but more tapered than before. There are no inversion points, as can be seen from the isotherm plot. The basal boundary layer has disappeared, advection no longer dominates over diffusion, but both compete with comparable amounts. Finally, when β = 1 (bottom of Figure 8) diffusion over-rides advection. This is why isotherms are essentially horizontal and temperature profiles linear in this case.

Fig. 7.

Distributions of temperature and (scaled) velocities for Case IV computations, using the accumulation function, Equation (20b), p₁ = 0.1, T_M = 0.5, T_M ⁽¹⁾ = 0.8, and T_M ⁽²⁾ = 0, and for α = 0.1, β = 0.01. For a detailed explanation of the figure see the legend to Figure 5.

We mention that the graphs for velocity as shown in figures (c), (d), (e), and (f) are identical to those when (α, β) = (10^–2, 10^–2). (10^–2, 10^–1). (10²,1), (10^–1,10^–1), (10^–1,1), (1.10^–2), (1.10^–1), and (1.1), reflecting the importance of sliding.

Are situations like those displayed in Figure 8 realistic? In view of Equation (3d), probably yes! One simply needs to lower the accumulation scale q _m by a factor of 10 and the representative thickness by a factor of 2 (q _m = 10 cm/a, d _g = 1000 m, both realistic values) to shift β into a range typical of Figure 8. This simply indicates that no single ice sheet is typical of all and scales must be carefully selected.

Fig. 8.

Temperature distributions for Case IV computations using the same conditions as described in the legend to Figure 6 but now for values α = 0.1, β = 0.1 (a, b) and α. = 0.1, β = 1.0 (c. d). respectively.

We also constructed solutions for α = 1, β = 1 and substantially smaller surface temperatures (T _M ⁽¹⁾ = 0.4), These solutions show a 10–30% contribution of gliding to the total velocity. The reason, in this case, is the occurrence of positive temperatures in the lower half of the ice sheet with an accompanying substantial enhancement of the rate factor a(T). The case is at most interesting but, clearly, not realistic. This does not mean that flow situations with considerable gliding would not occur. In fact, we will demonstrate later that, when μ _Ref/∊² is large (but still ⩽0.2), the influence of gliding is visible.

Finally, a note is in order regarding the magnitudes of heights, extents, and velocities in physical variables. We have not indicated the respective scales in Figure 7. The transformations obey Equations (1) and involve, among other things, the parameter ∊. Because dimensionless quantities do not (seem to) depend on this parameter, the values of physical variables can be adjusted largely by choice of ∊. It is not difficult to envisage situations for which longitudinal velocity components are considerably smaller than in the example of Figure 5.

4. Variation of the external forcing

Alterations in the external forces express themselves as variations of the accumulation-rate function, the surface temperature, and the geothermal heat.

1. Summary of results>

To judge properly the usefulness and validity of the model, let us briefly summarize the results. Inferences apply to dimensionless variables.

2. Implications of the restriction μ_Ref/∊² ˂ 0.2

We anticipate that in a model without the limitation of inequality (30) the thermo-mechanical coupling will be more evident and a sole dependence on μ _Ref /∊² may no longer be valid. While we will say a few words on such extensions in the next section, our goal here is primarily to make a few cautious remarks on earlier works (listed in Table I). We ignore steep glaciers and confine our attention to plane isothermal or non-isothermal ice sheets (cf. Reference Morland and JohnsonMorland and Johnson, 1982; Reference HutterHutter, [^c1983], chapter 7; Reference Morland and SmithMorland and Smith, 1984; Reference Morland, Morland, Smith and BoultonMorland and others, 1984; Reference Hutter, Alts, Niordson and OlhoffHutter and Alts, 1985). We have failed to find any statements as to the limitation of the models implied by inequality (30). In a personal communication to one of us, Morland has, however, expressed his opinion that close to the ice divide longitudinal stretching ought to be important. This is most likely the deeper cause of the imposition of inequality (30). It then remains a question of how reliable the constructed numerical solutions of Morland and associates can be.

A partial answer to this question is obtained if one notices that in the computational work of Morland and associates temperature was never an issue. It was therefore irrelevant whether integrations were performed with or against the flow, and these authors chose to start integrating at the margin, where U, H’, and a _cc were all finite nonzero numbers. Thus, local surface mass balance could be satisfied to a sufficient accuracy independent of whether sliding was substantial or not. (They only needed sliding linear in overburden pressure to have a finite margin slope.) Integration into the ice sheet must have proceeded for some distance without difficulty, but it is anticipated that satisfaction of the local surface mass balance became difficult in the vicinity of the ice divide, where H’ and U are small. The integrated mass balance must have been used to define the location of the ice divide. This was obviously at positions where H’ was numerically very small. So Morland and associates were able to obtain reliable solutions away from divides and proper overall geometries with their integration procedures, but solutions ought to be in some doubt close to the divide.

With temperature, integration must start at the divide, and then large errors are introduced right at the start, unless sufficient sliding is available.

3. Into new avenues

When solving the thermo-mechanical coupled problem, as we did with the Morland–Hutter model, there is no way around inequality (30). Thus, we may ask the questions: (i) How is the present work best extended? (ii) Are there glaciological situations for which the Morland–Hutter model is applicable without a restriction like inequality (30)? (iii) What model will have to be used to eliminate the inequality?

Answer (i). A direct and easy extension is the application of the Morland–Hutter model to axisymmetric situations. Cold cirque glaciers may be modelled this way. Equations must simply be expressed in cylindrical coordinates and azimuthal dependences be dropped. Restrictions in the form of inequalities are likely to persist and computations must be very similar. The physical range of inferences is, however, limited to the satisfaction of the inequality and may be limited.

A similar remark applies to time-dependent plane ice-sheet problems. Time rates of change of H(X, t) will emerge from the surface mass balance which must be accurately determined.

Answer (ii). Inequality (30) evolves from the satisfaction of the mass balance at the divide where U and H’ are small. For steep glaciers, it can be shown (see Reference HutterHutter, [^c1983], chapter 5) that UH’ at the head is non-zero and finite, independent of whether there is sliding or not. This makes it likely that plane glacier-flow problems can be solved with the Morland–Hutter model for a larger range of basal conditions permitting adherence or sliding.

Answer (iii). The Morland–Hutter model emerges from an asymptotic analysis of the basic viscous flow equations of isotropic ice (see Reference HutterHutter, [^c1983], chapter 3) by stretching coordinates. The “stretched” equations are valid approximations except for boundary-, edge-, and transition-layers. Since longitudinal stretching is ignored in the model but important in the vicinity of the divide, it follows that the transition layer is governed by the full unstretched equations. Local solutions of this transition layer must be matched, in principle, with the solutions of the Morland–Hutter model far away from the divide, and, since the steady-state system is now elliptic, the latter must be matched with local margin solutions. Numerically, such a solution procedure must be difficult. Our preliminary investigation indicates that solving the full original problem is just as easy.

Work of this kind will guide our future avenues.