A Plasticity Theory Approach to the Steady-State Shape of a Three-Dimensional Ice Sheet

Abstract

The differential equation determining the elevations of a perfectly plastic three-dimensional steady-state ice sheet is set up. Analytical solutions of the equation are obtained in two simple case, (1) an ice sheet on a horizontal base with an arbitrary edge curve, and (2) an ice sheet on a plane sloping base with a rectilinear ice margin. The solutions are discussed, particularly with reference to the development of ice divides and ice streams.

For arbitrary base and ice-margin geometries, solutions are obtained by means of the method of characteristics, which reduces the problem to solving simultaneously three ordinary first-order differential equations. The integration, which is performed by numerical methods, is generally commenced at the ice margin, where the necessary boundary conditions are known.

The method has been applied to model the elevation contours and the flow pattern of the central Greenland ice sheet, using the bottom topography revealed by radio echo soundings and the present ice margin geometry. The result is in surprisingly good agreement with our knowledge of the ice-sheet topography and flow pattern, all significant ice divides and ice streams being reproduced. This suggests, that the method can be applied to model the shape and flow pattern of ice sheets under glacial conditions, using information about former ice-margin positions.

Une approche par la théorie de la plasticité de la forme d’équilibre d’une calotte glaciaire tridimensionnelle. L’équation différentielle déterminant les altitudes d’une calotte glaciaire en équilibre sur trois dimensions, parfaitement plastique a été établie. Des solutions analytiques de l’équation sont obtenues dans deux cas simples: (1) une calotte glaciaire sur une base horizontale avec des bords de forme arbitraire, et (2) une calotte glaciaire sur une base plane inclinée avec des bords reclillgnes. Les solutions sont discutées particulièrement en fonction au développement des diffluences glaciaires et des courants de glace.

Pour des bases et des géométries de bordure arbitraires, les solutions sont obtenues par la méthode des caractéristiques qui réduit le problème à la résolution simultanée de trois équations différentielles ordinaires du premier ordre. L’intégration qui est obtenue par des méthodes numériques est généralement commencée sur la bordure du glacier où les conditions aux limites sont connues.

La méthode a été appliquée pour modéliser les lignes de niveau et le mode d’écoulement de la calotte glaciaire du Groënland central, en utilisant la topographie du fond révélée par les sondages à écho-radio et la géométrie actuelle des bordures. Les résultats sont èlonnemment cohérents avec ce que nous savons de la topographie de la calotte glaciaire et de son écoulement, toutes les diffluences glaciaires et les courants de glace significatifs ont été reproduits. Ceci suggère que la méthode peut être appliquée pour modéliser la forme et l’écoulement des calottes glaciaires au moment des glaciations en utilisant les informations disponibles sur ¡es dernières positions de leurs bordures.

Zusammenfassung

Eine Herleitung der stationären Gestalt einer dreidimensionalen Eisdecke mil Hilfe der Plastizilälslheorie. Es wird die Differentialgleichung zur Bestimmung der Höhe einer vollkommen plastischen, dreidimensionalen Eisdecke in stationärem Zustand aufgestellt. Analytische Lösungen der Gleichung erhält man für zwei einfache Fälle, nämlich (1) für eine Eisdecke auf horizontalem Untergrund mit beliebiger Umrandung, und (2) für eine Eisdecke auf ebenem, geneigtem Untergrund mit geradlinigem Eisrand. Die Lösungen werden besonders hinsichtlich der Entwicklung von Eisscheiden und -strömen diskutiert.

Haben Untergrund und Eisrand unregelmässige Form, so ergeben sich Lösungen mit der Methode der Charakteristiken, die das Problem auf die simultane Lösung von drei gewöhnlichen Differentialgleichungen 1. Ordnung zurückführt. Die Integration, die mit numerischen Methoden ausgeführt wird, beginnt im allgemeinen am Eisrand, wo die erforderlichen Randbedingungen bekannt sind.

Die Methode worde zur Berechnung der Höhenlinien und des Fliessverhaltens des zentralen grönländischen Inlandeises angewandt, wobei die durch Radar-Echolotungen ermittelte Topographie des Untergrundes und die gegenwärtige Lage des Eisrandes eingeführt wurden. Das Ergebnis stimmt erstaunlieh gut mit unserer Kenntnis der Form der Eisdecke und ihres Fliessverhaltens überein; alle wichtigen Eisscheiden und -ströme werden wiedergegeben. Die legt es nahe, die Methode zur Berechnung der Form und des Fliessverhaltens von Eisdecken unter eiszeitlichen Bedingungen, fussend auf den Kenntnissen über frühere Eisrandlagen, heranzuziehen.

1. Introduction

Recent years have shown an increasing interest in modelling the shape and flow of ice sheets under glacial conditions: Much information obtained from analysing deep cores from ice sheets can be properly interpreted only if past ice sheet elevations and ice flow can be reconstructed (Dansgaard and others, 1973, p. 36; Paterson. 1977). This is also a condition for understanding the widespread landscapes formed by glacial erosion (Sugden, 1978) and finally attempts on modelling palaeo-climates require information about the topography of former ice sheets (Hughes. 1979).

Two different methods have been applied in three-dimensional (3-D) ice sheet modelling. One method is based on prescribing the distribution of accumulation rate over the potential ice covered area, and letting the ice sheet build up governed by the differential equations for non-stationary glacier flow (Campbell and Rasmussen, 1970; Mahaffy, 1976; Jenssen, 1977). This method has mainly been developed to tackle problems of glacier response to changes in external parameters such as accumulation rate and temperature and will yield steady state shapes of ice sheets in an indirect fashion only. The other method is based on extrapolating estimated steady-state ice-sheet profiles—obtained either from inspection of existing ice sheets or from two-dimensional ice-sheet theory—up-stream either at right-angles to the ice margin (Sugden, 1977), or along curved flow lines of assumed course (Boulton and others, 1977: Hughes and others. 1977; Andrews, 1979). A major shortcoming of the Sugden method is that the important influence of bottom topography on ice flow is not taken into account, and although basal topography has been considered in the reconstructions of the latter authors, this is done on a rather subjective basis.

In this paper the problem is approached by developing a theory of a perfectly plastic, three-dimensional, steady-state ice sheet. We shall discuss a few simple cases with idealized bed and ice-margin topographies that can be treated analytically. In spite of their simplicity, they throw light on some fundamental aspects of three-dimensional ice sheet flow. Moreover, a method of solution applicable in the case of an arbitrary bed and ice-margin topography will be discussed. Next, the efficiency of the model will be demonstrated by comparing the reconstructed surface topography and flow pattern of the actual ice sheet of central Greenland with observations, and finally the limitations and potentialities of the model will be discussed.

2. Perfect Plasticity Theory

Application of plasticity theory in ice-sheet dynamics is not of recent date. In fact, in the very first attempt to predict the surface profile of a steady-state ice sheet theoretically (Orowan, 1949), the ice was considered a perfectly plastic material, i.e. the shear stress at the ice-bedrock interface was assumed not to vary along a flow line.

Neglecting the influence of longitudinal stress gradients and of transverse shear stress, this assumption leads to the following differential equation for the surface elevation E of the ice sheet (Nye, 1952):

(1)

where s is a coordinate measured horizontally along the flow line, positive in the opposite direction to ice flow, τ ₀ the yield stress of glacier ice (τ ₀ ≈ 1 bar = 10⁵ N m⁻²), ρ the density of glacier ice (ρ = 920 kg m⁻³), g the gravitational acceleration (g = 9.81 ms⁻²), and H the ice thickness. In the case of a two-dimensional ice sheet on a horizontal base. Equation (1) predicts the maximum ice-sheet thickness H ₀ to be

(2)

where L is the distance between the ice divide and the ice margin.

In central West Greenland the bed underlying the ice is almost horizontal, ice flow is to a good approximation two-dimensional and L = 500km. Applying this value of L, the formula above gives H ₀ ≈ 3200 m. to be compared with an observed ice thickness of 3150 (Gudmandsen, 1973). The perfect agreement is to some extent caused by choosing the yield stress of ice to be equal to 1 bar, which is somewhat arbitrary, since observations as well as more advanced theory show, that the mean value along a flow line of the basal shear stress varies in the range 0.5–1.5 bars, depending on mean accumulation rate and mean basal ice temperature. Nevertheless, the result encourages us to apply perfect plasticity theory to obtain first order approximations to ice-sheet profiles, in particular because observed or assumed regional variations of accumulation rate and basal ice temperature can be accounted for by changing the value of the yield stress for the various regions of the ice sheet. Guidance as to how to make this change, is provided by comparing ice-sheet profiles calculated by perfect plasticity theory with profiles obtained by applying a more realistic ice flow law.

If, for instance, the stress–strain-rate relationship for glacier ice is supposed to follow Glen’s law

where

is a strain-rate component, σ′ _ij is a stress deviator component, τ _e the effective shear stress defined by τ ² _e = ½σ′ _ij σ′ _ij , and T the ice temperature. B ₀ ≈ 0.5 bar³ year⁻¹, k ≈ 0.2 deg ¹ and n ≈ 3 are constants in the flow law of ice (Paterson, 1980, p. 13–16.) Furthermore, if the accumulation rate a (metres of ice per year) and the basal ice temperature Tare assumed to be constant along the flow line, then the following expression is obtained for the maximum ice thickness of a two-dimensional ice sheet on a horizontal base (Haefeli, 1961; Paterson. 1980, p. 20–21):

Comparing this expression with Equation (2) shows that the same value of the ratio of H ² ₀/L will be obtained from the two formulae, if we put

(3)

Now consider two regions of an ice sheet, where mean accumulation rates are a ₁ and a ₂ . and mean basal ice temperatures are T ₁ and T ₂ , respectively. In order approximately to account for these differences, the yield stresses to be applied for the two regions should have a ratio of

(4)

(see Equation (3)), assuming that the ice flow is strictly two-dimensional. Otherwise, the value of a applied should be decreased or increased in magnitude, depending on whether the flow is diverging or converging.

For example. Equation (4) shows that a surplus by a factor of two of the accumulation rate in one region of an ice sheet as compared to another can be accounted for by applying a 20% higher yield stress. A similar increase of τ ₀ is required to account for a decrease of the basal ice temperature of 4 deg.

Changing the yield stress in accordance with Equation (4) should ensure almost correct positions and almost correct relative elevations of ice divides determined by perfect plasticity theory. In order to make such changes, however, the distributions of accumulation rate and basal ice temperature must be known. As regards the present large ice sheets, the distribution of accumulation rate is fairly well known, whereas basal ice temperatures must in general be estimated because they have been measured at only very few locations. In the case of a glacial age ice sheet, both parameters obviously have to be estimated.

Even though positions and relative elevations of ice divides can thus be correctly determined, at least in principle, the elevation profiles obtained from perfect plasticity theory will show certain deficiencies. This is illustrated by comparing the profile of a two-dimensional ice sheet obtained by perfect plasticity theory, with the corresponding profile obtained by applying Glen’s law and assuming a horizontal base and uniform accumulation rate and basal ice temperature along the profile (Figure 1), If the two profiles are adjusted to the same elevation at the ice divide, it appears that the former profile on the average underestimates ice sheet elevations by 12% (Paterson, 1972, p. 891).

Fig. 1. Theoretical two-dimensional ice-sheet profiles. Assumptions: horizontal base and uniform distribution of accumulation rate and basal ice temperature.

In spite of such imperfections, it seems worthwhile to design a three-dimensional perfectly plastic ice-sheet model as an objective method of reconstructing the three-dimensional topography and ice-flow pattern of former ice sheets, using bed and ice-margin topographies as an input. Moreover, considered as a second-generation model, it may serve as guidance on how to develop more advanced third-generation three-dimensional models able to account for regional variations of accumulation rate and basal ice temperature.

3. The differential equation

The basic assumption of the theory is that flow lines may be defined as trajectories to the elevation contours of the ice-sheet surface, and that Equation (1) holds along such flow lines. Along a trajectory the surface gradient dE/ds is given by the equation

(5)

where x and y are orthogonal coordinates in a horizontal plane.

Introducing a quantity of dimension length by the equation H _f = τ ₀/ρg and substituting H for E − B, where B = B(x, y) is the elevation of the base of the ice sheet. Equations (1) and (5) may be combined to give the following differential equation for the elevations of a perfectly plastic three-dimensional ice sheet:

(6)

Equation (6) may be solved by means of the method of characteristics. Applying the notation and p = ∂E/∂x and q = ∂E/∂y, Equation (6) may be rewritten

(7)

assuming that the x-axis is oriented in such a way, that p ≥ 0.

According to Kamke (1965, p. 66–67), the characteristic equations of the partial differential Equation (7) are

(8)

(9)

and

(10)

where p may be expressed in terms of x, y, E, and q through Equation (7).

The problem of solving the non-linear partial differential Equation (6) is thereby reduced to solving simultaneously three first order differential equations, of which Equation (8) defines the course of the flow-line projections on the xy-plane (from now on these projections will be simply referred to as flow lines), and Equation (9) and (10) define the variations along the flow lines of the surface elevation and its gradient in the direction of the y-axis. It should be pointed out. that the flow lines are the only set of curves in the x–y plane along which the elevations of the ice sheet can be determined by integration (in mathematical terms: the basic characteristics along which the integral surfaces of the partial differential equation can be build up). Consequently, in order to determine a specific ice sheet surface, a set of initial values of y, E, and q must be known at one point of each flow line. Since the origin and the orientation of the coordinate system may be chosen arbitrarily, this means that the surface elevation and one surface gradient must be specified at one point of each flow line.

It can be shown that the ice-sheet elevations are uniquely determined within the area confined by any specified elevation contour line. On the other hand, this is not necessarily the case as regards the area outside a contour line, where divergence regions may develop down-stream of points of intersection between the contour line and possible ice divides, see Figure 2. In such divergence regions, the determination of the surface elevations will require further specification. However, no such troubles will arise, if we choose a point of the ice margin as the starting-point for the integration. At the ice margin the initial values are E = B and q = ∂E/∂n = ∂B/∂n, if the differentiation is performed in the direction of the marginal curve.

Fig. 2. Ground plan of ice sheet. Definition of local coordinate system.

In general, the ice margin will not coincide with an elevation contour line. Nevertheless, the flow lines will intersect the ice margin at right angles, because the ice-sheet surface has a vertical tangent plane at any point of the ice margin, see Equation (1).

4. Ice Sheet on a Horizontal Base

In the case of an ice sheet on a horizontal base we may choose B = 0 and consequently E = H, without losing generality. Then combining Equations (9) and (10), we get

which condition must hold along any flow line. By integration we get

(11)

where C is a constant. Now consider an arbitrary flow line, which intersects the marginal curve of the ice sheet at a point O and let us introduce a local cartesian xy -coordinate system with origin at this point and x- and y-axes along the tangents to the flow line and the marginal curve respectively (see Fig. 2). Then at point O we have H = 0 and q = ∂H/∂y = 0 and consequently, by Equation (11), C = 0. Since at other points of the flow line, H is not equal to zero, renewed application of Equation (11) will show that q is everywhere equal to zero. Then from Equation (8) it appears, that the direction of the flow line is determined by the equation dy/dx = 0. which shows (when combined with the above definition of the coordinate system) that the flow lines of a perfectly plastic three-dimensional ice sheet on a horizontal base are straight lines running normal to the ice margin. Moreover, it can be deduced from Equation (9). that the elevation profile along any flow line is the well-known second order parabola, valid for the two-dimensional case,

(12)

These conclusions are in agreement with the results deduced by Nye (1952) for an ice sheet on a horizontal base. It follows that the surface elevation contours are parallel to the ice margin. Hence, convex parts of the marginal curve create divergence regions of the ice sheet up-stream, which may develop ice divides running inland from the margin, whereas concave parts of the marginal curve create convergence regions, which may feed ice streams, see Figure 3a. Except for a few minor corrections, this ice-sheet topography is identical to the one suggested by Sugden (1977) in his reconstruction of the morphology of the Laurentide ice sheet. But. the present investigation shows that the method applied by Sugden neglects the influence of bottom topography on ice-sheet elevations and ice flow. This may be justified where the relief of the base is small compared to the ice thickness. However, this is seldom the case close to the ice margins where ice thicknesses are small, and where bottom hollows or ridges may create convergence or divergence regions, respectively, which may greatly influence flow directions even in the interior parts of the ice sheet (see Section 6). Hence, Sugden’s approach is at best a rough approximation even in the case of smooth bottom relief.

Fig. 3. Central Greenland ice sheet. Heavy full lines are ice divides, medium full lines are elevation contour lines with elevations in 100 m (upright figures), thin full lines are flow lines. Dotted lines are elevation contour lines of ice sheet base with elevations in 100 m (sloping figures). (a) shows surface contours and flow pattern obtained by perfectly plastic ice-sheet modelling with no consideration of basal topography, (b) is observed surface and bed contours according to radio echo soundings (Gudmandsen, 1978).

Figure 3a shows the result of a three-dimensional modelling of the ice sheet in central Greenland when the influence of bottom topography is neglected. As shown in Appendix B, this does not necessarily require the base of the ice sheet to be a horizontal plane. A base given by the equation

(13)

Where E ₀ is a constant level and F is a constant factor, will result in the same simple ice sheet morphology as the one connected to a horizontal base. Equation (13) was introduced by Weertman (1964), who considered an ice sheet with its base, originally horizontal, in local isostatic equilibrium, so that F = ρ/(ρ − ρ _r), where ρ and ρ _r are densities of ice and rock respectively. If ρ = 900 kg m⁻³ and ρ _r = 2700 kg m⁻³. we get F = −0.5.

However, Equation (13) may be used with any F < 1 to approximate the shape of the bedrock. In Figure 3a we have chosen E ₀ = 500 m and F = −0.20, which results in a bottom elevation close to sea-level in the central part of the ice sheet, in accordance with observations. Furthermore, H _f in Equation (12) has been assigned a value of 10 m corresponding to a yield stress of τ ₀ = 0.9 bars.

The deficiencies of the model applied appear from a comparison between Figure 3a and b, of which the latter shows observed elevation contours in central Greenland (Gudmandsen, 1978). In Figure 3a the central ice divide is located midway between the western and eastern ice margins, whereas the ice divide is actually offset more than 100 km to the east (Figure 3b), mainly due to the influence of the bottom topography. Also the secondary ice divides separating the drainage basins feeding the major ice streams are malpositioned, these drainage basins being in general underestimated whereas the drainage basins of the intervening divergence sectors are overestimated.

As illustrated by Figure 12, which shows the result of applying the model discussed in Section 6, a much better approximation to the ice-sheet morphology is obtained when the influence of bottom topography is taken into account. Before dealing with this model, however, the influence of bottom topography will be discussed in more general terms, by considering an ice sheet with a rectilinear ice margin resting on a plane sloping base.

Fig. 12. Reconstruction of the central Greenland ice sheet by perfectly plastic ice-sheet modelling with due consideration of the topography of the base as given in Figure 3b. Heavy full lines are ice divides, medium full lines are elevation contour lines with elevations in 100 m, thin full lines are flow lines. Dotted lines are ice divides, deduced from ice sheet elevations, as obtained during radio echo-sounding flights (Gudmandsen, 1978).

5. Ice Sheet on a Plane Sloping Base with a Rectilinear Ice Margin

A coordinate system is introduced with origin at an arbitrary point of the ice margin, the y-axis along the projection of the ice margin onto a horizontal plane, and the positive direction of the x-axis pointing into the ice sheet. In this coordinate system the plane sloping base of the ice sheet is supposed to be given by the equation

(14)

where the maximum upward slope of the plane is β > 0 and makes an angle with respect to the x-axis of

(see Fig. 4).

Fig. 4. Ice sheet on a plane sloping base.

Since neither the conditions at the ice margin proper nor the topography of the ice-sheet base relative to the margin varies from one flow line to another, all flow lines and all ice-sheet profiles along the How lines must have the same shape. This means, however, that the ice-sheet surface must be insensitive to a translation in the direction of the ice margin, i.e. H must be independent of y. The equation of the surface elevations may therefore be written

(15)

where the elevations of the base are given by Equation (14).

Introducing dimensionless quantities by

,

, the differential Equation (7) may be rewritten

(16)

whereas Equation (15) should be replaced by the equation

(17)

The solution for an arbitrary value of α in the range

will be given in Appendix A. For the discussion of the fundamental aspects of the flow of a three-dimensional ice sheet over a non-horizontal base, it is sufficient to consider the case of α = 0, i.e. the case where the ice margin follows the direction of maximum slope of the base. Putting α = 0 in Equation (17) and combining this equation with Equation (16), we get

(18)

which defines the variation of the ice thickness with

. Equation (18) may be integrated to give

(19)

where

must fall within the range 0 ≤

≤ 2, to ensure a non-imaginary square root. In fact, the range is reduced to 0 ≤

≤ 1, as shown below.

Consequently, the ice-sheet profile of any section normal to the ice margin follows that of a quarter ellipse. Now from Equations (17) and (19) we get

= 1 and

= (1–

)/(2

−

) respectively. Substituting these expressions into the dimensionless equation analogous with Equation (8), the following differential equation of the flow line is obtained:

from which we may deduce that

and

are related by the equation

(20)

Furthermore, for the increment of the dimensionless distance

along a flow line, we have

from which we may deduce that

and

are related by the equation

(21)

Here

is measured from the ice margin, positive in the opposite direction of ice flow. Substituting in Equation (19) for

by means of this equation, we find that the variation of the ice thickness with distance along the flow line is

(22)

For

→ ∞ Equations (20), (21), and (22) show that ỹ → ∞,

→ 1, and

→ 1.

Hence all flow lines have a common asymptote parallel to the y-axis at a dimensionless distance from the ice margin of

= 1. The corresponding non-dimensionless distance is x ₁ = H _f /β ². Furthermore, the ice thickness approaches asymptotically the value H = H _f /β (which is the ice thickness in uniform ice flow down a constant slope β), as the direction of the flow lines approaches that of the y-axis, i.e. as the slope of the bottom in the direction of the flow line approaches the value β.

The equation of the elevation contour lines are obtained by combining Equation (17) for α = 0 and Equation (19):

Taking Ẽ = Ẽ ₀ in this equation we get after a little manipulation

(23)

which shows that the elevation contour lines are arcs of circles with dimensionless radius 1 and centre at point (1, Ẽ ₀). In Figure 5a, flow lines and elevation contour lines are shown for the ice sheet considered. Also shown by Figure 5b are the surface, bottom, and ice thickness profiles along a flow line and by Figure 5c the ice-sheet profile of a section normal to the ice margin.

Fig. 5. Ice sheet on a plane sloping base with an ice margin following the direction of the maximum slope. Dimensionless representation. (a) shows ground plan with surface-elevation contours (heavy full lines with upright figures), flow lines (thin full lines), and bed elevation contours (dotted lines with sloping figures), (b) shows surface, bottom, and ice thickness profiles along a flow line, and (c) shows ice-sheet profile of section normal to the ice margin.

The flow pattern of Figure 5 is clearly divergent. We shall now analyse this flow a little closer:

In the case of three-dimensional flow, the following differential equation is valid along a flow line (personal communication from Sigfus Johnsen):

(24)

In this equation q is the ice flow function q(s) = HV _m where V _m is the average over the ice thickness H of the horizontal velocity component in the direction of ice flow, s is the distance along the flow line, positive in the opposite direction of ice flow, a is the rate of accumulation, and R = R(s) is the radius of curvature of the elevation contour lines at their points of intersection with the flow line, positive or negative depending on whether the centre of curvature is located up-stream or down-stream relative to the flow-line point considered. Assuming the rate of accumulation to be constant, we introduce a dimensionless ice flow function by the equation

Noting then, that

(see Equation (23)), Equation (24) may be written in dimensionless form as

With the boundary condition

= 0 for

=

, this equation may be integrated to give

from which the value of the dimensionless ice flux at the ice margin may be determined by putting

= 0:

This relationship is plotted as the heavy curve in Figure 6. For comparison the corresponding relationship for the case of parallel flow

is shown as a thin line in Figure 6.

Fig. 6. Dimensionless representation of the ice flux q_e at the ice margin as a function of flow line length s_e, assuming uniform accumulation rate along flow line. Heavy curve corresponds to the flow pattern of Figure 5. Thin line represents parallel flow.

The slower increase of the heavy curve reflects the divergence of the ice flow pattern shown in Figure 5.

As previously mentioned (see p. 441), there is an upper limit to the extent of the portion of an ice sheet that can possibly drain towards a margin parallel to the maximum slope of its plane base. Taking H _f = 10 m. Table 1 presents the limiting distance x ₁ = H _f/β ² as a function of β.

The regions of an ice sheet at distances greater than x ₁ from the margin considered, must necessarily drain towards other margins. In the case of an ice sheet unlimited in both x- and y-directions, for example, ice flow in the zone x > x ₁ will be parallel flow in the negative y- direction, with an ice thickness corresponding to uniform flow down a slope of inclination β.

In the case of a band-shaped ice mass limited by two parallel margins, there are three possible flow patterns, depending on the ratio of the width W to the limiting distance x ₁.

• (a) W/x ₁ < 2, see Figure 7a.

Fig. 7. Band-shaped ice sheet on plane sloping base with ice margins in the direction of the maximum basal slope. (a) is for the case W/(H/β²) < 2, (b) is for the case W/(H//β²) = 2, and (c) is for the case W/(H/β²) > 2. Heavy full curves are surface-elevation contours, thin full curves are flow lines, and dotted curves are bed elevation contours.

Two regions of the type just mentioned will develop with flow towards either margin. The two regions are separated by a central sharp ice divide, where elevation contour lines from either region intersect at an angle different from 180°. The cross-sectional profile of the ice sheet consists of two elliptically shaped arcs, intersecting at a point of the sharp ice divide.

Table I. Limiting distance x ₁ as a function of maximum base slope β for H _f = 10 m

(b) W/x ₁ > 2, see Figure 7b.

Flow pattern as described above, except that the ice divide is smoothed, elevation contours intersect the ice divide without a bend, and the cross-sectional profile is a full semi-ellipse.

(c) W/x ₁ > 2, see Figure 7c.

In this case three different flow regions will develop, the two regions mentioned at (b) being separated by a region of parallel flow in the direction of the negative y-axis W −2x ₁ wide and with an ice thickness corresponding to uniform flow down a slope of inclination β.

It is questionable if the ideal-plastic band shaped ice mass model has more than illustrative significance. When considering the shape of its ground plane, one might think of applying it to a glacier tongue. However, one of the basic assumptions of the model, viz. that the flow lines coincide with the trajectories of the elevation contours, will in general not be fulfilled in the case of a glacier tongue, where the two sets of curves may deviate due to bottom sliding in the longitudinal direction of the glacier. In the case of a cold glacier tongue frozen to its base, the model may be of some relevance.

The solution obtained for a three dimensional ice sheet on a plane sloping base—in particular the general solution for an arbitrary value of the angle α given by Equations (A–4), (A–7), and (A–8) of Appendix A, opens up the possibility of numerous exercises of modelling ice sheets with various shapes of the marginal curves. For example. Figure 8a shows an ice sheet with a rectangular ground plane. The ice sheet is divided into four flow regions, two with parallel flow towards north and south respectively and two with curvilinear divergent flow towards east and west respectively. The regions are separated by sharp ice divides bisecting the angles between surface contour tines as well as flow lines. Since the ice-sheet surfaces of the different regions are determined by simple analytic expressions, simple equations may also be deduced in this case for the ice divides. Even in the case of a curvilinear ice margin, Equations (A–4), (A–7), and B(A–8) can be used to construct ice-sheet surfaces, since a specific pair of values of Ẽ ₀ and α can be attached to each point of the ice margin, determining the course of and the ice thickness variation along the flow line terminating at the point considered.

Fig. 8. Ice sheet on a plane sloping base. Dimensionless representation. Heavy full curves are ice divides, medium full curves are surface-elevation contours with elevations indicated by upright figures, thin full curves are flow lines, and dotted curves are bed elevation contours with elevations indicated by sloping figures. (a) Ice sheet with rectangular ground plan. (b) Ice sheet with circular ground plan.

The flow lines and elevation contour lines of the ice sheet with a circular ground plan shown in Figure 8b have been constructed this way.

A more interesting use of the theory than constructing ice sheets on plane sloping bases with hypothetical ground plans, however, is to study the influence on a three-dimensional ice sheet flow of basal ridges and troughs running transversely to the ice margin.

Consider an ice sheet resting on a base composed of two sloping planes forming a ridge and a trough respectively (see longitudinal sections of Figures 9 and 10). The rectilinear ice margin is supposed to run in the direction of the maximum slope of the planes, so that Equations (20). (22), and (23) will apply for each of the planes composing the respective ice-sheet bases.

Fig. 9. Ice-sheet flow around a basal transverse ridge. Dimensionless representation. In ground plans the thick full line is an ice divide, medium full lines are surface elevation contours with elevations indicated by upright figures, and thin full lines are flow lines. Dotted lines are base elevation contours with elevations indicated by sloping figures. (a) Basal ridge composed of planes of equal inclinations. (b) Base composed of one horizontal and one sloping plane.

Fig. 10. Ice-sheet flow in a basal transverse trough. Dimensionless representation. In ground plan thick full lines indicate flow lines separating the central convergence region from neighbouring divergence regions, medium full lines are surface-elevation contours with elevations indicated by upright figures, and thin full lines are flow lines. Dotted lines are base elevation contours with elevations indicated by sloping figures. Shaded area is catchment area used for determining the ice flux function shown in Figure 11.

(a) Basal ridge

If the two basal planes have equal inclination β, a sharp ice divide will form above the junction of the basal planes (Figure 9a). Ice flow on either side of the ice divide will be divergent flow of the type discussed on p. 442. If the inclinations of the basal planes are different, the ice divide will be gradually displaced towards the plane with the smaller slope. This is illustrated in Figure 9b for the case of β ₁ = 0 and β ₂ = β ≠ 0, The equation of the ice divide can in this case be determined by means of Equations (20), (21), (23), and (12) as

(b) Basal trough

In this case three flow regions will develop (Figure 10). Two of them, which are divergence regions of the type considered in Figure 5, are bounded by flow lines, terminating at the point of intersection O between the ice margin and the bottom line of the sub-ice valley. This point is a singularity, being the common end point of all the flow lines occupying a convergence sector of the ice sheet located between the above mentioned limiting flow lines. The bottom line of the trough is a flow line itself, along which the ice-sheet elevation profile is given by Equation (12). All other flow lines of the convergent sector are tangents to the valley bottom flow line, following its course for some shorter or longer distance before terminating at the singularity point. The course of the flow lines of the convergent sector and the ice-sheet elevation profiles along these flow lines cannot be expressed by simple analytic equations. The flow lines shown in Figure 10 have been determined by numerical integration of Equations (8), (9), and (10), as described in Section 6. Figure 10 illustrates the significance of basal transverse troughs as generators of sectors of converging ice flow, which may feed ice streams.

This is further illustrated by Figure 11, which in dimensionless form shows the total ice flux Q at point O as a function of the “inland” extension x _e of the convergence sector (heavy curve). If the annual accumulation rate a is assumed to be uniformly distributed over the convergent region, the dimensionless total ice flow Q is simply determined as the shaded area shown in Figure 10. i.e.

For comparison, the corresponding

relationship for the case of parallel flow is shown as a thin line in the figure,

being the total ice flow in a band of dimensionless width

_0.02 = 0.0054, which is the width of the convergence sector shown in Figure 10 at a dimensionless distance from the ice margin of

= 0.02.

Fig. 11. Dimensionless representation of the total ice flux Q at the ice margin as a function of the “inland” extent x_e of the ice-sheet sector considered. Heavy curve corresponds to the sector of converging flow shown in Figure 10. Thin curve represents parallel flow in a band of dimensionless width

= 0.0054 corresponding to the convergence sector width at a dimensionless distance from the ice margin of

= 0.02.

Even transverse valleys with slopes of very moderate inclinations will have a considerable effect on the flow pattern of an ice sheet.

Taking for example H _f = 10 m, β = 0.002 (200 m/100 km), a = 0.3m of ice/year and x _e = 500 km corresponding to a dimensionless

of 0.2, the total ice flux at point O of Figure 10 may be calculated to be 27 km³ of ice/year, a figure of the same order of magnitude as the total ice flux of the Jakobshavn ice stream in central West Greenland (Carbonell and Bauer, 1968), which has a catchment area with a and x_e values of the same order of magnitude as the values applied in the above calculation.

For comparison, in the case of parallel flow, the total ice flow through a section of width W_0.02 = 13.5 km would be 2 km³ of ice/year only, if the same values of a and x _e were applied. (A sector width of 13.5 km occurs in the convergent flow case at a distance from the ice margin of 50 km where the ice thickness according to the model is 1000 m.) In addition to demonstrating the important role of sub ice ridges and troughs as generators of sectors of divergent and convergent ice flow respectively, Figures 9 and 10 also illustrate a condition which is quite general for three-dimensional ice-sheet flow, viz. that a flow line, when followed up-stream from the ice margin, will always be diverted towards the direction of maximum up slope of the bed. The degree of diversion for a given magnitude of the slope is the larger the less the ice thickness.

Therefore, in general, the curvature of the flow lines will be greatest in the marginal zones of an ice sheet, and the flow lines consequently become more and more rectilinear towards the centre. In an actual ice sheet, this statement, which is based on perfect plasticity theory, has to be modified as regards flow lines approaching secondary ice divides. Such flow lines will often show large curvatures due to strong divergence of the flow round the divide, see the discussion in Section 7.

Also deviations from the above-mentioned general flow pattern may occur, for example in the case of large elevation changes of the sub-ice surface in the interior of an ice sheet.

6. Arbitrary Sub-Ice and Ice-Margin Topography

The analytic solutions obtained in Sections 4 and 5 for simple cases of three-dimensional ice-sheet flow have mainly illustrative significance, but are rather useless for modelling surface elevations and flow patterns of real ice sheets with their far more complicated sub-ice and ice-margin topographies. In order to tackle this more interesting problem, a method of solution based on numerical integration of Equations (8), (9), and (10) has been developed.

As discussed in Section 3, the ice-sheet elevations and flow patterns are constructed as a synthesis of calculated variations along individual flow lines. Also it was mentioned, that the calculation was most conveniently commenced at the ice margin, thus ensuring a unique solution for the entire ice-covered area.

The choice of the ice margin as the starting-point of the integration is the more convenient, since a potential application of the model is to reconstruct surface elevations and flow patterns of ice-age ice sheets, for which the most reliable geometric feature is probably the very position of the ice margin.

However, when initiating the numerical integration at the ice margin, one is faced with a problem due to the ice thickness H _m there being equal to zero and as a consequence the surface gradient p _m = [∂E/∂x]_m being infinitely high.

This is a rather unpleasant boundary condition for numerical calculations. But the difficulty is overcome by taking the first step of integration by means of the equations for an ice sheet on a plane sloping base, i.e. Equations (A–4), (A–7), and (A–8) of Appendix A. The maximum upward slope β of the basal plane and the angle α between the direction of this slope and the ice margin, which are needed for applying these expressions, are determined as the corresponding values of the tangent plane to the bed surface at the ice margin point considered. As soon as this first step has been performed, the integration can be continued without further troubles by standard numerical integration methods.

For each flow line the calculations are performed in the local xy-coordinate system shown in Figure 2. Along each flow line the positions of points with specified elevations (e.g. in 200 m steps) are marked, and eventually the topography of the ice sheet is constructed by drawing surface contour lines through the points of equal elevation of the various flow lines. Locations of ice divides are in principle determined by connecting intersection points of flow tines where elevations match.

A more practicable method, however, is to connect the sharp bend points of the surface contour lines. In the case in which the same H _f-value has been applied for the two ice-sheet sectors meeting at an ice divide, the divide will bisect the angle between the respective elevation contour-line segments. Otherwise the ice divide will be displaced towards the sector with the greater value of H _f and the consequent greater ice-sheet surface gradients.

In connection with the numerical integration procedure, it remains to be mentioned that the function B = B(x, y) describing the topography of the ice-supporting surface, is evaluated, by two-dimensional spline interpolation in x and y directions respectively between discrete bed surface elevations obtained by digitalizing the bed surface in a rectangular grid. The mesh-width of the grid can be varied, allowing the point density to be adjusted in accordance with the roughness of the bed surface relief.

The model has been applied to construct the surface topography and flow pattern of the ice sheet of central Greenland, see Figure 12, applying a constant value of H _f = 10 m over the entire region. The sub-ice topography applied in the model is shown in Figure 3b. It is based on radio echo-soundings carried out during the Greenland Ice Sheet Program (GISP) (Gudmandsen, 1978). The bed topography is extrapolated through the regions close to the ice margins, where soundings are scarce, if not missing, so as to give a smooth transition to the relief of the ice-free surface ahead of the ice margin. The mesh-width of the grid applied, varies from 40 to 100 km. with the denser point distributions in the marginal areas. Also shown in Figure 12 by dashed lines are the positions of the ice divides from Figure 3b, deduced from the ice-sheet elevations obtained during the radio echo-sounding flights (Gudmandsen, 1978).

As appears on comparing Figures 3b and 12, there is a general agreement between the constructed and the observed flow pattern of the ice sheet, all significant ice divides and ice streams being reproduced, though not always in exactly correct positions. The largest descrepancies between constructed and observed positions occur in the case of secondary ice divides in East Greenland. This is probably due to the fact that the 40 km mesh width of the grid, applied for digitalizing the bed topography, is too coarse to provide a sufficiently detailed representation of the very rough relief of the bed supporting this part of the ice sheet. However, the relatively modest density of radio echo-sounding flight lines in the area hardly justifies application of a grid of finer mesh-width.

As regards elevations, predicted values at the central ice divide are in close agreement with the observed ones (around 3200 m). This of course is due to choosing H _f = 10 m, corresponding to an average bottom shear stress of 0.9 bar. What is the more interesting, is the close agreement between modelled and observed elevation changes along the central ice divide, the tendency to a local dome formation around Crête, and the very modest increase of elevation from Crête to Summit (Mock, 1976) being correctly reproduced.

As to the elevations on the slopes of the ice sheet, these are generally underestimated by up to 200 m, as was to be expected (see discussion in connection with Figure 1).

One aspect of the modelling should be specially emphasized, viz. the predominant role of the convergence sectors for determining the over-all surface elevations of the ice sheet (see e.g. the drainage basins of Jakobshavns Isbræ, Daugaard Jensen Gletscher, and Kangerdluqssuaq in Figure 3c) as opposed to the far more subordinate role of divergence sectors. Since convergence sectors are generally associated with the lower points of the ice margin, we may consequently conclude, that surface elevations in the interior parts of an ice sheet will in general be determined in relation to the parts of the ice margin with the lowest elevations.

7. Discussion and Conclusions

7.1. Flow pattern around ice divides

As already discussed in Section 2, the ice-dynamics problem is entirely separated from mass balance considerations in perfect plasticity theory. As a consequence, only the flow pattern (i.e. the course of the flow lines and the ice thickness distribution along them) will be determined by the theory, whereas the ice flux will be left undetermined in the first place. However, once the flow pattern is constructed, the ice flux can be determined by one of the methods described on p. 441 and p. 446, assuming that the distribution of mass balance over the ice sheet surface is known.

The differential equation method (p. 441) and the area method (p. 446) will respectively, provide spot values and average values over ice-sheet depth sections of the ice flux, from which ice flow velocities averaged over depth can be determined by dividing by the ice thickness.

Rigorous application of the differential-equation method, however, will lead to absurdities. This is due to the fact that in constructions of ice sheet surfaces using three-dimensional perfect plasticity, not only central ice divides, but also secondary divides separating individual drainage basins (the counterparts of watersheds in unglacierized areas) will be reproduced as sharp ridges, from which flow lines will originate (see Figure 12). Since the surface slope in the direction of an ice divide will always be less than the slope along a flow line at its starting-point at the divide, the basal shear stress in the direction of an ice divide will always be less than the yield stress τ ₀ and, consequently, the perfect plasticity theory allows no flow along a divide.

This is obviously unrealistic when considering “watershed” ice divides, such as those separating the drainage basins of the central West Greenland ice sheet (see Figure 12). Along these divides the basal shear stress is only diminished by an insignificant amount, as compared to the shear stress along the almost parallel flow lines, and actual flow velocities along these weakly developed ridges will not differ significantly from the flow velocities of the adjacent ice.

Therefore, ice fluxes (and thereby velocities) obtained from a three-dimensional perfectly plastic ice-sheet surface by the differential-equation method, have to be smoothed around ice divides, the more so for less distinctly developed ridges. The essential problem can be referred to the basic assumption of perfect plasticity theory, viz. that the basal shear stress τ is constant in the direction of a flow line, and that ice will not flow if τ < τ ₀. As discussed above, this assumption will cause not only primary ice divides but also secondary drainage divides to be sharp ridges with zero flow in the direction of their crests. In actual ice-sheet flow all divides are smoothed, ice will flow along the divides, and all flow lines will originate at ice domes or horizontal ice ridges. However, the smaller the slope along a divide, the larger will be the divergence of the flow, and the lower will be the flow velocity. One might say that ideal plasticity theory provides a caricature of the surface topography and the flow pattern of an ice sheet, all features being exaggerated.

In order to reproduce the actual flow pattern more closely a more realistic model should be applied, e.g. a model analogous to that mentioned on p. 433 (i.e. a steady state Mahaffy model (Mahaffy, 1976)).

This however, is not trouble-free either, first because the differential equation to be integrated is far more complicated, and secondly because flow dynamics will then be directly coupled to the flux of ice. which in turn depends on the mass-balance distribution and the flow pattern along the entire up-stream part of a flow line. This means, that the (in many respects) very convenient technique of integrating from the ice margin inland along a flow line, will not apply. The integration will have to be initiated at ice domes or at horizontal ice divides, the positions of which are not always accurately known, not even as regards existing ice sheets, but certainly not in the case of an ice-age ice sheet.

Furthermore, one may question, whether the method of integration by means of characteristics can be applied at all to a steady-state version of the Mahaffy equation.

7.2. The ice margin

Although the construction of the surface topography and the flow pattern of a perfectly plastic three-dimensional ice sheet is most conveniently initiated at the ice margin, this does not mean, of course, that the ice margin is the determinant of the flow in the interior of the ice sheet. On the contrary, the position of the ice margin is controlled by the mass-balance distribution at the ice-sheet surface and the relief of the bed underlying the ice.

Therefore, when the flow pattern of an ice sheet has been determined on the basis of estimated ice-margin positions, the total mass balance of the constructed drainage basins should be calculated, which implies the distribution of the mass balance over the ice-sheet surface is known. If imbalance is found for any drainage basin, the ice marginal position should be changed and the flow pattern recalculated with the new position of the ice margin. Only when all drainage basins are in balance with the surface mass-balance distribution, can the ice-sheet construction be claimed to be self-consistent. Unfortunately, in the case of ice sheets under glacial conditions this checking of the model can only be tentative, due to the lack of a reliable model for predicting the distribution of precipitation over past ice-sheet surfaces. A further difficulty in performing mass-balance checks is that large masses of ice may be ablated from the ice sheet through ice streams which may either terminate in iceberg-calving fronts or may be feeding ice shelves. Consequently, the mass balance of the entire ice-sheet/stream or ice-sheet/shelf system should be considered. This requires several additional ice-dynamics problems to be solved, among others how to predict reliably the positions of the seaward margin and the grounding line of former iceberg-calving ice streams or ice shelves respectively.

7.3 Depression of ice-sheet base

A potential application of the theory is to reconstruct the topography and the flow pattern of ice sheets under glacial conditions from information about former ice margin positions, and a paper on this is in preparation for the cases of the Late Wisconsinan Greenland and Innuitian ice sheets. In this connection the deflection of the Earth’s crust due to the increased load from the ice sheet must be considered.

Several models describing the deflection history of the crust in response to changing load conditions caused by growing or decaying ice sheets have been developed (see e.g. Brotchie and Silvester, 1969; Farrell and Clark, 1976; Peltier and Andrews, 1976).

From these models it may be concluded that in the equilibrium stage, i.e. so long a time after loading that rates of deflection are insignificant, the magnitude of the depression can be reasonably well approximated by the value corresponding to local isostatic equilibrium, even though small discrepancies will occur in particular close to the boundaries of the loaded area.

On this basis, it is concluded, that at the present moderate stage of sophistication of steady state three-dimensional ice sheet modelling, crustal deflections are sufficiently accurately accounted for by considering local isostatic equilibrium, and this the more so, since such a crustal deflection model can be immediately built into the ideal-plastic three-dimensional ice-sheet model (see Appendix B). Based on the present relief of the formerly ice-covered landscape, the topography of the ice-sheet surface and the topography of the depressed bed may be calculated simultaneously.

7.4 Conclusion

The above discussion has aimed at clarifying the potential and the limitations of the three-dimensional perfect-plasticity ice-sheet model. To summarize: The model does not directly account for the influence of the ice flux on the topography and on the flow pattern of the ice sheet, although it is known to have a significant, though not predominant influence. If, on the other hand, the distribution of surface mass balance and bottom ice temperature is known or can be estimated, regional variations of these parameters can be accounted for by changing the assumed basal shear stress from one region to another, see p. 433. Since some uncertainty in such estimates will probably always remain, the three-dimensional perfect-plasticity ice-sheet model will provide the gross elevation and drainage pattern of the modelled ice sheet but not the exact ice-surface elevations, nor the exact positions of ice and drainage divides.