A. Background and Governing Equations
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
balance laws of mass, momentum, and energy,
and the
constitutive relations for heat flux and stress deviator.
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
a mass-balance condition (kinematic surface equation),
prescription of shear and normal traction (stress-free conditions),
prescription of surface temperature,
and at the fixed base:
prescription of the geothermal heat flow,
application of a basal sliding law.
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 (σ0, ρgd 0), velocities (q m), a typical depth (d 0), a characteristic length (d 0/∊), a typical strain-rate (D 0), a temperature range ΔT, etc. (Table II). For instance, σ0 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 0, 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 × 105 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[acc (·) 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).
B. A Brief Description of the Code, its Versatility and Limitation
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 0, W 0, T 0, 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 HD is at X = XD = 0 and the right margin is at XR. So ice-sheet semi-length is L = X – XR = XR. 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) = d2H/dX 2(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 acc 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 acc 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.
C. Results – A Parameter Study
As mentioned previously, a parameter study must incorporate variation of the scales (s, δ, α, β), the constitutive properties of the ice and the external forcings through the climate and the geothermal 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
1, and p
2 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
1 = p
2 = 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
1 and/or p
2 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 acc
= 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 (1) < 0.1, as well as 0 < T M (2)< 0.1. When T M (2) = 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−5 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 1 > 0 and/or p 2 > 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 = 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 = 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 ∊2, 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 1 which is the slope of the snow-equilibrium height, parameterized for a range of values of equilibrium height H 0 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 1. Nevertheless, the dependence on p 1 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 HD (left) and ice-sheet semi-length XR (right) plotted against p1 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 p1 where iterative computations failed or were not pursued, or required some interpolation. Ice-divide heights and extents depend critically on both 
and p1.
table v. ice-sheet semi-lengths and divide heights, and their absolute and relative changes with the accumulation parameter p1, 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 1 = 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 1 → 0. Extrapolated curves in Figure 4 must therefore be taken with care. The case p 1 = 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 
, p1 = 0.1, TM = –0.1, TM = (1)0.8, Te
M
(2) = 10−2. 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 – UF), 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(= ∊ 2 ) and θ assume values as shown in Table VI. Evidently, the ∊2 matrix is symmetric with respect to the main diagonal, whereas the θ matrix enjoys this symmetry with respect to the second diagonal. Large ∊2 values lie at the right lower corner of the ∊2 matrix, large θ values at the left lower corner of the θ matrix.
table vi. values of ∊2 (left matrix) and θ (right matrix) for ranges of values for s and δ. the ∊2 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 ∊2 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
1 = 0.1, p
2 = 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 ∊2 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/∊2 > 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 HD, the ice-sheet semi-length XR, and the aspect ratio 
plotted against μRef/∊2 in semi-logarithmic representation for Case III computations. The scale for HD is shown on the left-hand side; those for XR 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/∊2. 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/∊2. The shallowness of an ice sheet now depends on e and μ
Ref/∊2 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/∊2 was kept constant (= 0.05) and s < 1, δ < 1 were arbitrarily varied such that δs–1 < 103. 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/∊2 ration and higher surface temperatures.
It thus appears that in the range of the validity of this model, μ
Ref/∊2 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), 
p1 = 0.1, TM = 0.5, TM
(1) = 0.8, and TM
(2) = 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). (102,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 (1) = 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/∊2 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.
(a) Accumulation.
Patterns according to Equations (20a, b, c) were analysed under various thermal conditions (keeping everywhere the temperature below or at melting) and varying μ
Ref/∊2 within the interval (5 × 10–2, 2 × 10–1). The snow-equilibrium height was chosen according to Equation (21) with 
, p
1 = 0.1, and p
2 = 0.
Consider Equation (20c) first, because it is the simplest. Varying the parameters a and b (Case V) permits us to analyse changes in ice-sheet flow due to global climatic influences. The value of b is a measure of the solid precipitation at high altitudes. Alternatively, a measures how fast melting decreases with altitude variation of the surface-energy budget.
Figure 9 displays H D (left) and X R (right) as a function of μ Ref/∊2 for the case that a = 0.5 is held fixed while b is varied. Symbols represent points that were calculated; curves are eye-fitted. Evidently, H D grows with increasing a; but this same statement is not necessarily true for X R. Qualitatively, the larger a is, the more melting will occur per unit distance in the ablation zone. Thus, for a fixed extent of the ablation zone, the total ablation will increase with growing a, resulting in an imbalance of mass and thus needing an increase in ice-divide height until mass balance is re-established. Because each change of H D also induces a change of X R and, consequently of the extent of the ablation zone, such a plausibility argument does not apply to X R.
Fig. 9. Graphs of the ice-divide height HD (left, (a)) and the ice-sheet semi-length XR (right, (b)) plotted against μRef/∊2 as obtained from Case V computations using the accumulation function. Equation (20c), as shown in the inset of figure (b). Symbols represent computed values, curves (solid, dashed, and dotted) are interpolated by eye. Question marks indicate ranges of μRef/∊2 values for which no numerical solution was found (loss of convergence).
Important to observe in Figure 9 is the consistency of the computed data, as they permit a clear identification of a trend. (The interpolated curves are well defined by the data and the curves are clearly separated.) This is an indirect a posteriori proof that the computer code is well behaved and is reliable. Notice, however, that at fixed μ Ref/∊2 the difference of two values of X R for two values of a is larger than for H D, indicating a relatively poor sensitivity of the secant root-finder to the different H D values.
Variations of a are also manifest in changes of surface geometry. Besides the aspect ratio, the “bluntness” of the ice sheet also suffers some changes. Figure 10 gives some information on these parameters. For fixed μ Ref/∊2, the aspect ratio increases with increasing a (Fig. 10a) and so does the bluntness of the ice sheet (Fig. 10b). By this, we mean that, for two ice sheets with the same ice-divide height and ice-sheet semi-length, one ice sheet may appear more or less tapered while the other may be more or less blunt. Intuitively, decreasing a should make ice sheets more and more tapered. Figure 10b corroborates this but the effect of changing a by a factor of 4 is hardly visible. The case with a = 6.25 was also run but the profile was hardly distinguishable from the solid line in Figure 10b.
Fig. 10. (a) Graphs of the aspect ratio 
plotted against μRef/∊2 evaluated for the Case V computations on the basis of Figure 9a and b. Symbols indicate computed values; curves are interpolated, (b) Graphs of the surface geometry normalized as described in the legend to Figure 5 for Case V computations and μRef/∊2 = 0.05. The symbols refer to the values of a, and the accumulation function is the same as that shown in figure (a).
We have also varied the accumulation parameter b in Equation (20c), leaving the ablation parameter a = 6.25 fixed. Results are summarized in Figure 11. Evidently, changes of b by a factor of 2 (= 100%) result in very small absolute changes of H D but rather large changes of X R (cf. Table VIII). All these findings have obvious practical bearings and reflect the unfortunate sensitivity of the model. An analogous behaviour in a different context has already been observed by Reference NyeNye (1963) in his kinematic wave theory. There, it was deduced that relatively short temporal records of glacier-snout movement permitted reliable evaluation of the local climate function (the ablation at the snout) but that long temporal records of accumulation were needed to predict reliably short glacier-snout movement.
table viii. ice-divide heights H D, their differences ∆H D, ice-sheet semi-lengths X R, and their differences ∆X R for the three runs of figure 9 with μRef/∊2 = 0.08
Fig. 11. Graphs of the ¡ce-divide height HD (left, (a)) and the ice-divide semi-length XR (right (b)) plotted against μRef/∊2 as obtained from Case V computations using the accumulation function Equation (20c). as shown in the inset to figure (b). Symbols represent computed values curves (solid, dashed, and dotted) are interpolated by eye. Question marks indicate ranges of μRef/∊2 value for which no numerical solution was found (loss of convergence) or the accuiacy was questionable.
The considerable variation of ice-sheet semi-length with accumulation conditions also has its influence in the temperature distribution and the flow. In Figure 12, we display the isotherm plots for six different cases of the accumulation function, Equation (20c), given by the variations of a and b, as indicated in the figure. The graphs indicate that monotonic changes of b or a do not necessarily yield monotonic changes in the temperature distribution. This is so because ice-sheet semi-lengths may equally suffer non-monotonic changes under the same conditions (see Figs 9b and 11b).
Fig. 12. An explanation of the considerable temperature change with accumulation. Shown are isotherm plots determined by using the accumulation function, Equation (20c). Thermal and basal conditions are the same as for Case I (see Table IV). Figures (a), (b). and (c) (left) show the isotherm distribution when the accumulation parameter b is varied; figures (d), (e), and (f) show those when a is varied.
As a last test of the model to variations in the accumulation condition, results obtained with the accumulation functions, Equations (20a–c), were compared, when a = 12.5 and b = 0.5. The three accumulation functions are nearly equal in this case. Their function values agree at the ice divide, provided that H D > 0.75 and in the ablation area below the snow-equilibrium line but differ otherwise. For a prescribed H(X), a cc(Equation (20c)) is largest in these instances, followed by a cc(Equation (20b)) and a cc (Equation (20a)). So, one would expect H(a cc (Equation (20c))) > Ha cc(Equation (20b))) > H(acc(Equation (20a))). However, as we can infer from Figure 13, this string of inequalities is not borne out by the model, the obvious reason being that each ice sheet selects its own profile H(X). From a practical point of view, these results are rather disturbing, because they imply that a reliable prediction of the ice-sheet geometry requires a prescription of the climate more accurate than can be provided by the best climatological studies. That the ice-sheet extent is more sensitive to small changes in the accumulation pattern is yet more disturbing.
Fig. 13. The same as Figures 9 and 11 but now comparing ice-divide heights and ice-sheet semi-lengths for the three different accumulation functions as indicated in the inset. Symbols represent computed values; curves (solid, dashed, and dotted) are interpolated by eye. Question marks indicate ranges of μRef/∊2 values for which no or uncertain numerical values were obtained.
(b) Surface temperature and geothermal heat. Thermal conditions of the climate are manifest in the temperature function T s; the magnitude of the activity of the Earth’s core and the thickness of the mantle are expressible as function values of the geothermal heat. Here, we are not so much interested in spatial variations but in the influence of these quantities by changing their order of magnitude. We thus vary in Equations (13) and (7a) the magnitudes of T M (1) and θz. For Case II computations (see Table IV, in which T M (1) now varies), results are shown in Figure 14a–c, For very cold surface temperature and “normal” basal conditions, the temperature distribution has the usual inversion pattern. Warming up the surface temperature weakens this inversion property but does not destroy it. For this, β values had to be increased. This is a behaviour one might have expected.
Fig. 14. Graphs of isotherm distributions for Case II computations (see Table IV) keeping all external parameters fixed except the surface temperature, (a) is for “cold”; (b) for “intermediate”; and (c) for “warm” atmospheric conditions.
A qualitatively comparable effect is also exhibited by changes in the geothermal heat. While ice-sheet geometries and velocities are hardly affected, the temperature distribution is strongly influenced. For calculations with α = 0.5, β = 0.02, T M (1) = 0.6, μ Ref/∊2 = 0.2, and the accumulation function, Equation (20c), a = 6.25, b = 0.5, isotherms for two different values of the geothermal heat are shown in Figure 15. For θz = 1 (geothermal temperature gradient of 1°C/100m; Fig. 15a and b), temperatures are everywhere negative and profiles show a pronounced inversion pattern. When the geothermal heat is 2.5 times higher and all other conditions are kept fixed, then the entire ice sheet is substantially warmer, temperature inversion is less pronounced, and close to the snout basal temperatures are even positive (Fig, 15c and d). This case is not realistic and would have to be described by a polythermal model.
Fig. 15. Isotherms and vertical temperature profiles for computations using Equation (20c) as accumulation function and α = 0.5, β = 0.02. TM (1) = 0.6, a = 6.25, and b = 0.5 for two different values of the geothermal heat. The top two figures are valid for θ2 = 1 (geothermal temperature gradient 1°C/100 m); bottom figures hold for a geothermal heat flow 2.5 limes larger. Notice in the profiles (figures (b) and (d)) that the temperatures are everywhere negative for the small value of geothermal heat but locally positive (above melting) in the second case.
5. Rheological non-linearities: creep-response function and rate factor
Rheological non-linearities enter the formulation through Equations (4a–c) and are manifest in the stress-stretching relationship, Equation (2e), and in the viscous dissipation term of Equation (2h).
(a) Creep response. Parameters which govern the creep response are the shear stress τ and the scaling θ; the temperature T determines the values of the rate function a(T) and will be analysed below. Figure 16a shows a perspective representation of g(τ, θ), the coordinate scale for τ being linear, those for θ and g being logarithmic. Function values grow rapidly for growing τ and θ spanning several log cycles in the displayed τ and θ range. Figure 16b delimits the domains of the (τ, θ) plane in which the nonlinear terms of g(τ, θ), Equation (4c), are significant. Large θ values make the influence of non-linearities more likely. This is the reason why we chose s and δ values often in ranges where θ = 1, 10, 100, and seldom θ = 1000, and cannot be surprised that for 100 < θ < 1000, approximately, we were often facing overflow problems. The graph in Figure 16c displays the function values of ω(θτ2). Here, too, large θ values enhance the non-linearities.
Fig. 16. (a) Axonometric representation of the creep-response function g(τ, θ). The scales for θ and g are logarithmic, that for τ is linear. Function values vary substantially, in particular when θ is large. Dashed and dotted lines delimit the ranges when linear terms in τ alone, or linear plus cubic terms, or all terms of the g function are significant (cf. Equations (4b, c) and figure (b)). (b) Domains in the (θ.τ) plane within which non-linearities in the creep-response function g(θ, τ) are significant. To the left of the line “only linear terms significant”, the non-linear terms contribute less than 10%. To the left of the curve “linear and cubic terms significant”, the term in g(θ, τ) involving τ5 has less than 10% influence, (c) Graph of ω(θτ–2) in doubly logarithmic representation. For linear rheology, ω has the value 0.3336. Non-linearities are almost always active.
Computations were thus performed for 1 < θ < 200. To reduce the role of basal friction as much as possible, μ Ref/∊2 = 0.2 was chosen which is right at the upper bound of the applicability of the model. Four runs were performed, in which all essential parameters were kept constant with the exception of θ (cf. Table IX). Dimensionless ice-divide heights (and ice-sheet semi-lengths) were essentially unaffected by these variations. However, the velocity and also to a lesser extent the temperature distributions were affected by changes in θ values. With
we define the maximum of the ratio of the gliding velocity to the basal velocity. It represents a measure of the importance of the viscous deformation in comparison to the sliding mechanism. For the computations of Table IX, the results revealed
corroborating the intuition that large θ values will enhance viscous deformations. Parallel with an increase of θ goes a growth of the total velocity. Figures 17 and 18 summarize the results for θ = 1 and θ = 200 (all other parameters being as in Table IX). The temperature distribution (graphs (a) and (b) in Figures 17 and 18) is only slightly affected by changing θ from 1 to 200, but the total velocity (graphs (c)) and the gliding velocity (graphs (d)) show considerable changes. For θ = 1, gliding is much less conspicuous than for θ = 200. In addition, the distribution of the gliding part of the velocity is quite different in the two cases. For large θ values, it is concentrated to the outer region of the ice sheet. The profiles of the transverse velocities show no surprising features (graphs (e)) but flow lines (graphs (f)) and velocity scales indicate that the magnitudes of the velocity vectors are enhanced when θ is increased.
Fig. 18. The same as Figure 17 but now for a large value of θ(= 200). Note, in particular, that gliding is here enhanced in comparison with sliding.
table ix. characteristic parameters s, δ, θ, (μ Ref/∊2, α, β 0 for the scalings and θZ, T M (1) , a, b for the external forcings. note that systematic changes in θ values leave ice-divide heights nearly unchanged
Having found the conditions for which gliding becomes appreciable, additional computations were conducted, in which θ was kept fixed and large (θ = 100), but the thermal parameters α and β we have plotted values of ;were varied. It was found that the effects of α and β may compete. For instance, large α (corresponding to large strain heating) and small β (corresponding to low thermal convection) may cause similar temperature and velocity profiles. For this to occur, θ values must be large. Proof that enhancing α at large θ values raises the temperature everywhere in the ice sheet and thus also makes gliding more important is given in Figure 19, where results are shown for computations of Table IX (with δ = 10–1) but α = 0.1 and α = 0.05, respectively. This figure shows isotherms (top), vertical profiles of longitudinal velocities (middle), and difference (gliding) velocities (bottom), on the left for a dissipation number α = 0.1, which is twice as large as on the right, α = 0.05. On the left, temperatures and velocities are larger than on the right. A similar effect can also be observed when α is kept fixed but β is varied. In this case, temperatures seem to change more conspicuously than velocities.
Fig. 19. Distributions of isotherms ((a) and (b)), total longitudinal velocities ((c) and (d)), and gliding velocities ((e) and (f)) for a set of computations, using the accumulation function, Equation (20c), H°e = 0.5, p1 = 0.1. TM (2) = 0, and scaling parameters as indicated in Table IX for θ = 100. but α = 0.1 (left-hand figures) and α = 0.05 (right-hand figures).
(b) Rate factor. Having so far found almost no influence of the temperature variation on the ice-sheet geometry and to some extent also on the flow field, it seems compelling to analyse the role played by the rate factor at some greater depth. In particular, more should be known about the notions “weak” and “strong” thermal coupling. To this end, Reference Smith and MorlandSmith and Morland’s (1981) rate factor, which is based on Reference GlenGlen’s (1955) and Reference Mellor and TestaMellor and Testa’s (1969) uniaxial compression data, is replaced by the exponential relationship
with various different values of β*. Notice that T in Equation (38) is the dimensionless temperature (a negative number for cold ice) which vanishes at 273.15 K. Moreover, at T = 0 the value of the rate factor evaluated by Equation (38) agrees with that of Smith and Morland’s equation (4a).
In Figure 20 we have plotted values of a(T) for various values of β* between 0°C and –40°C. Also shown is the graph of a(T) using the Smith and Morland function. Since T is negative, large values of β* will lead to small values of a(T). Thus, in this temperature range β* > 3.51 corresponds roughly to a thermo-mechanical coupling which is stronger than that of ice, while β* ˃ 3.51 makes the coupling even weaker. Directions into strong and weak coupling are also indicated in Figure 20. Note, however, that in a smaller temperature range, larger β* values separate domains with weaker and stronger coupling than ice.
Fig. 20. Plots of the rate factor a(T) against temperature in semi-logarithmic representation. Solid lines show the graphs for a(T) according to Equation (37), the dashed–hatched curve gives the graph for the Smith and Morland relationship. The values of β* used in Equation (37) are also indicated. The large arrows indicate the direction into which graphs of rate factors would move if thermal coupling were to be made stronger and weaker, respectively.
In an attempt to push conditions as far as possible to the limit where thermo-mechanical coupling effects would be maximal, several runs were conducted under the conditions stated in Table IX (but with θ = 100 only) and using Equation (38) as the expression for the rate factor and varying β* in the range stipulated by Figure 20. Results for β* = 3.5177, which is roughly representative for ice, and β* = 1.215, which corresponds to very strong thermo-mechanical coupling, are summarized in Figures 21 and 22, The graphs of the isotherms and vertical temperature profiles (plots (a) and (b) in Figures 21 and 22) indicate that enhancing the coupling raises the temperature slightly throughout the ice sheet. Longitudinal velocities are also enhanced when β* values are lowered (compare graphs (c) in the figures) but, most conspicuously, strong thermo-mechanical coupling enhances the part of gliding in comparison to the total velocity (compare graphs (d)). In fact, while M, defined in Equation (37) is less than 10% for Figure 21, it amounts to more than 30% for Figure 22. Differences in thermo-mechanical coupling are hardly seen in the profiles of the vertical velocities except close to the margin (graphs (e)). Finally, the flow lines also suggest an enhancement of the outward flux in the near-margin zone when β* is small, However, ice-divide heights and semi-lengths remain unchanged.
Fig. 21. Distributions of temperature (isotherms (a), vertical temperature profiles (b)) and scaled velocities (graphs (c)–(f)) for computations using the accumulation-rate function, Equation (20c), and rate factor, Equations (38) with β* = 3.5177. as well as H°e = 0.5, p1 = 0.1, TM (2) = 0, and scaling parameters as indicated in Table IX for θ = 10. For a detailed explanation of the figure, see the legend to Figure 5.
Fig. 22. The same as Figure 21 but for θ = 100, β* = 1.215.
D. Discussion and Conclusions
In the preceding sections we have outlined the reduced continuum model proposed by Morland and Hutter; we have described a computer program with the aid of which symmetric plane, steady ice-sheet flow was analysed by using interactive methods and rather sophisticated graphics hard- and software in order to find deductions implied by the model. Two types of parameter were varied, first, characteristic numbers describing the mechanical and thermal behaviour of isotropic ice and, secondly, the external forcings responsible for the climatic input (surface temperature, accumulation rate, and the geothermal heat). Attention was focused on gross behaviour but not on the finer details. We could, for instance, also have investigated spatial variations; however, in view of the results, it is felt that these are of lesser significance.
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.
-
(i) The most important result of this paper is probably the recognition that the Morland–Hutter model with our integration procedure is only applicable if sufficient sliding is present. This is expressed in inequality (30) (or (31)), which states that the basal drag must be small: μ Ref /∊2 0.2. Because decreasing μ Ref is tantamount to increasing the slipperiness of the basal bed, the model requires a minimum amount of sliding. There is some competing factor against this slipperiness, namely in the aspect ratio, so that for thinner ice sheets more frictional resistance can be permitted. The permitted range of μ Ref /∊2 is, however, somewhat restrictive as no-slip (a realistic limit) for cold ice sheets can never be attained. Nevertheless, the soft ice in a small basal layer, which is often of sub-grid size, will lead to μ Ref values that do not violate inequality (30) for a particular ice sheet. So the range of applicability of this model is larger than if sliding would solely have to represent true slippage along the base. Reference Morland, Morland, Smith and BoultonMorland and others (1984), in their isothermal Greenland ice-sheet model, used : μ Ref /∊2 = 9.9 in violation of inequality (30). All subsequent inferences are drawn, subject to the satisfaction of inequality (30).
-
(ii) Of all characteristic parameters that describe the thermo-mechanical behaviour of the ice sheet (namely ∊2, θ, α, β, μ Ref ), the ice-divide height and the ice-sheet semi-length only depend on μ Ref/∊2 but not on the other possible independent combinations, provided the external forces remain unchanged. Neither does the shape nor the velocity distribution depend on characteristic parameters other than μ Ref /∊2. This statement requires qualification insofar as an increasing contribution from sliding to the total velocity can be observed when θ is large. The influence is, however, seldom larger than a few per cent.
Apart from surface temperature and geothermal heat, the temperature distribution depends on both α and β, but dependence on βis much stronger than that on α. For small β, diffusion is small and vertical advection of heat dominates except in a basal boundary layer. Increasing β makes diffusion more and more important. Also, for small β, temperature profiles may show inversions, which are typical and (probably) caused by longitudinal advection. The role of the parameter a is marginal unless large α values are paired with large θ values and relatively uniform temperature profiles. In those circumstances, we have observed a certain influence of α and θ on the flow and temperature distribution, mostly insofar as the enhanced strain heating gave rise to higher temperatures close to the margin (where the τs are large) which, alternatively, caused enhanced gliding there.
-
(iii) Besides μ Ref /∊2, the accumulation-rate function is the dominant factor in determining the ice-divide height, the ice-sheet semi-length, and the taperedness of the sheet. Spatial variations of the accumulation-rate function may change the geometry drastically. It was, however, also often observed that small changes in ice-divide heights were accompanied by fairly large changes in ice-sheet semi-lengths. This implies that high numerical accuracy is required if relatively subtle changes of ice-sheet extents are to be predicted. This is unfortunate as double precision has already been used in this program and better accuracy necessarily requires higher-order finite differencing.
-
(iv) Rheological non-linearities are hardly ever seen to be important. This is so, because a θ-dependence coupled with a thermo-mechanical coupling was seen in the velocity field only to have an effect of a few per cent. At stronger thermo-mechanical coupling (much beyond that of ice), the influence of the temperature on the flow was substantial.
2. Implications of the restriction μRef/∊2 ˂ 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 /∊2 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.
E. Acknowledgements
The plotting program for our graphics terminal was written by B. Rutherford. We are grateful for his help in performing the plotting routines. The art work was done by C. Bucher and I. Wiederkehr, and the manuscript was typed by M. Staub.
While performing this work we were supported in part by the National Science Foundation through contract No. DPP 8219439.
We thank L.W. Morland and an anonymous referee for their critical reviews of this manuscript.
