General theory

Consider a laterally unconfined ice shelf or glacier tongue with axes as shown in Figure 1. The x and y axes are horizontal at sea-level and z is measured vertically. The ice shelf flows out to sea in the positive x direction and its thickness H and velocity u vary slowly over its length. The height of the upper surface above sea-level is denoted by h.

Fig. 1. Longitudinal section of an ice shelf.

To describe the deformation processes taking place within the ice mass we shall use the strain-rates ⋵ _ij and the stresses σ _ij (i, j = x, y, z), and we shall assume these to be related by the flow law for polycrystalline ice (Reference PatersonPaterson, 1969, p. 82):

where ⋵ and τ are defined through the second invariants of the strain-rate and stress-deviator tensors,

and

with stress deviators σ _ij ’ defined by

and

The parameter B is dependent on temperature and may therefore be expected to vary with depth in the ice shelf. Reference ThomasThomas (1973[a]) found that n is about three for Antarctic ice shelves; this is also the value expected if creep processes are controlled by dislocation glide.

For an unconfined ice shelf there is complete symmetry for all stresses and strain-rates in the x and y directions, at least along the centre line, so that horizontal shear quantities are zero. Where the thickness changes slowly, vertical shear quantities are also negligible (Reference Sanderson and DoakeSanderson and Doake, 1979).

In order to find the equilibrium profile we propose to repeat for ice shelves the exercise Reference NyeNye (1959) carried out for ice caps, that is, to find expressions for velocity and strain-rate as functions of distance x and thickness H, apply a steady-state requirement, and then solve the equations to determine the resulting profile.

Velocities and strain-rates may be found by considering the analysis given by Reference Weertman.Weertman (1957). We consider the static equilibrium of stresses acting on any vertical element far from the ice front, and require that the total force within the ice shelf acting on any vertical column must balance the total force exerted by sea-water pressure. If the density of sea-water is ρ _w this means that

where s and b refer to the upper and lower surfaces of the ice shelf. For a variable ice-shelf density ρ(z) this leads to

and therefore

where

the average of B over depth.

If we know, or can model, the functions B(z) and ρ(z), then Equation (1) gives ⋵ _xx as a function only of thickness H. For clarity, let

Now, the velocity of the ice shelf may be written as

where u(o) is the velocity at some chosen point of origin, for instance the hinge zone.This may be written as

so that we now have velocity and strain-rate in terms only of x and H. We now consider the condition for steady state. If the ice shelf undergoes accumulation rate a and melt rate m (positive for melting), then the condition for steady state may be found by considering a moving, deforming column of ice. Thickness gradients in the y direction are negligible so that the continuity equation is

where a and m are expressed in metres of pure ice per unit time, ρ _i is the density of pure ice, and

the average density through the ice shelf.

Now, by symmetry ⋵ _xx = ⋵ _yy , and substituting Equations (2) and (3) into Equation (4) and re-arranging we have a differential–integral equation solely in terms of x and H,

This equation cannot be solved analytically, but in principle it gives thickness as a function of distance; it enables an equilibrium profile to be calculated. It is most readily solved numerically by tracing the course of the curve using small steps, for, if we take some values of H(o) and u(o) at a starting point, we can calculate the slope ∂H/∂x by Equation (5) and hence can calculate the increment δH for some appropriately small interval δx. We can then carry out the same process for our new point H(δx), u(δx), and so on to the values at the rth step, H _r and u _r . In such a process of finite steps we have to calculate the integral term in the denominator as a finite sum; for the rth step

and as we make each further step we add the increment

The full expression for the increment in passing from one step to another is, in the limit of small intervals,

Before we can use this, we shall need to consider the form of the density integrals within the equation.

Density–depth relationship

Reference ThomasThomas (1973[a]) noted that we should not be tempted to assume that density is constant through an ice shelf. Such an assumption leads to an overestimate of τ by a factor of about two and therefore an overestimate of ⋵ by a factor of about eight.

Reference SchyttSchytt (1958) has considered different ways of modelling the variation of snow density with depth. Using data gathered at Maudheim he first fits a simple decaying exponential

where d and b are empirical constants, and then goes on to provide two other possible functions which both follow the Robin suggestion for a more realistic model of the processes controlling densification (Reference RobinRobin, 1958, p. 92) where “the proportional change in air space in firn is in linear proportion to the change in stress due to the weight of material above a given level”. This assumption leads to a model of the form

where β and γ are empirical constants.

Reference SchyttSchytt (1958, p. 124) compares these models with field data. All three models must be considered to give adequate agreement in view of the wide scatter of the data. In evaluating the double integrals in Equation (6) it is found that Equation (7) is integrable while Equation (8) is not. We shall therefore follow Reference ThomasThomas (1973[a]) in adopting Equation (7) with constants fitted appropriately.

How should the constants b and d be fitted? The constant d is readily related to the average density of surface snow ρ(h), which we may assume to remain constant along the ice shelf. We have

b is less easy to fit. Schytt found a value for b of 0.025 8 m⁻¹ at Maudheim, for an ice-shelf thickness of about 200 m. This value implies that the average density through the ice shelf is 0.827 Mg m⁻³. In extrapolating Schytt’s results to other ice shelves we therefore have two alternatives: either we assume that the attenuation factor b remains constant as thickness varies, in which case the average density of the ice shelf must vary, or we assume that the average density of the ice shelf remains constant as thickness varies, in which case b must vary. Studies on Erebus Glacier tongue (Reference HoldsworthHoldsworth, 1974) showed no systematic variation of average density while the thickness varied by a factor of two, and we shall therefore assume average density to be constant. I have found that trial calculations using the alternative assumption do not lead to radically different profiles.

Let us require then that the ice-shelf density be a constant,

; thus,

and on performing the integration in Equation (7) we find that

which we wish to solve for b as a function of H. An exact solution cannot be given analytically, but a simplification can be made if we consider the magnitude of exp (—bH). Reference SchyttSchytt (1958) quotes a value of b = 0.025 8 m⁻¹ for Maudheim, with H ≈ 200 m, so exp (—bH) ≈ 0.006 and is negligible. We can therefore write

Now consider the double integral in Equation (6). Using Equation (7) it becomes

using our approximation. We find then that Equation (6) becomes, after substitution and simplification,

with C defined through

We shall now examine how well this formula agrees with experimental data.

Agreement with reality

The Erebus Glacier tongue, McMurdo Sound, Antarctica, has been examined in some detail by Reference HoldsworthHoldsworth (1974). It extends 12 km into the Ross Sea, is about 1 500 m wide, and is between 120 and 360 m in thickness. It is believed to be floating freely. Using a variety of data sources Holdsworth determined thicknesses, velocities, and strain-rates along the ice tongue and used these to calculate melt-rates after making the assumption of steady state and taking a value for the flow-law parameter B . In using these values to calculate an equilibrium profile we are therefore using a circular argument; however, the calculation of the profile does at least demonstrate whether the values found are internally consistent and give the expected behaviour between data points. We shall also see that it provides some information on behaviour in the vicinity of the hinge zone.

Holdsworth gives the following values:

We take

and

We have adjusted d so that Equation (7) gives a density of about 0.60 Mg m⁻³ at a depth of 10 m. This accords with density measurements made by Reference Stuart and BullStuart and Bull (1963) on the nearby Ross Ice Shelf. We find that d has the value 0.41 Mg m⁻³.

We must choose a starting point in order to calculate an equilibrium profile. First, we try fitting to data at the hinge zone. Extrapolating Holdsworth’s data we take the values

The resulting profile, found by taking incremental steps 25 m long in the x direction, is shown in Figure 2. It shows poor agreement in that it thins too rapidly. We therefore try starting at alternative points which are 1, 2, and 3 km from the hinge zone. For starting points near the hinge zone, agreement is still poor but it is found that starting points further than about 3 km from the hinge zone all give acceptable agreement. At 3 km we take and we compute the profile by working both forward and backward from the 3 km mark. Figure 2 shows that the agreement here is very good from the 3 km mark onwards, both in thickness profile and in velocity. Any discrepancies fall well within the error bars of the measured values. This shows that the model is essentially correct and that the measured data and calculated parameters are internally consistent over the greater part of the glacier tongue.

Fig. 2. Comparison of models with (a) ice thicknesses and (h) ice velocities measured by Reference HoldsworthHoldsworth (1974) for Erebus Glacier tongue. The solid line shows Holdsworth’s assumed fit to his data points. The dashed lines show model A, in which the model is fitted to data at the hinge gone, and model B, in which the model is fitted to data at the 3 km mark.

Behaviour near the hinge zone

At distances less than 3 km from the hinge zone our model departs from the measured data. We find that the ice shelf has to thicken very rapidly towards the hinge zone and its velocity falls dramatically. At about 0.75 km from the hinge zone the model breaks down completely, giving an infinite thickness and zero velocity. This discrepancy must in some way be connected with the process of transition between a land glacier and a floating ice shelf, a process which is not at present adequately understood. There are two principal changes taking place during this transition. First, the velocity distribution through the ice must change from the form characteristic of a land-based glacier (in which velocities decrease from the surface downwards becoming zero at the ice-bedrock interface) to that characteristic of an ice shelf (in which negligible shear occurs and the velocity is uniform with depth). We expect a complicated stress and velocity field in the area where this change is taking place (Reference HoldsworthHoldsworth, 1977). Secondly, the temperature at the bottom generally undergoes a large change when it comes into contact with sea-water rather than rock. This change may well be accompanied by bottom freezing.

If the change in velocity distribution is the cause of the anomaly, then our model gives us information about its extent and general character. It extends for about 3 km which corresponds to a period of 20 to 25 years. We expect these figures to be dependent on the thickness and velocity of the ice shelf in question. The nature of the effect seems to be that as ice flows from the land into the sea it initially undergoes less rapid thinning than would be expected. This is surprising, since a qualitative consideration (Reference HoldsworthHoldsworth, 1977) of the velocity field involved leads us to expect thinning which is more rapid as the slow-moving basal ice accelerates to the same velocity as the main body of ice. Data from Maudheim Is-shelf (Reference RobinRobin, 1958) and Ross Ice Shelf (Reference Robin, Robin, Swithinbank and SmithRobin and others, 1970; Reference HughesHughes, 1975) appear to confirm Holdsworth’s picture.

Let us consider whether bottom freezing could be responsible for the thickness anomaly. It is found that in order to make the model satisfy the condition that H = 370 m and u = 102.5 m a⁻¹ at the hinge zone while H = 272 m and u = 137.5 m a⁻¹ at the 3 km point we have to assume bottom freezing at the rate of 3.0 m a⁻¹ of ice, which amounts to a total accretion of some 75 m over the distance in question. Considering that oceanographic conditions near the hinge zone are unlikely to be very different from conditions over the rest of the glacier tongue, we must account for the freezing by looking to the source of cold provided by the cold land glacier as it becomes afloat.

Reference RobinRobin (1955) calculates the temperature T(z) through a land-based ice sheet as a function of accumulation a, surface temperature T _s, and constant temperature gradient (∂T/∂z)_b at the base. His formula is

where ξ = z−H + h, the height above the bottom of the ice shelf, and k is the thermal diffusivity. We have used

Let us take the average inland accumulation rate to be 0.2 m a⁻¹ (Reference Stuart and BullStuart and Bull, 1963). We take

and

The resulting temperature profile is shown in Figure 3. In calculating temperature changes we shall be interested in the temperature profile near the bottom of the ice shelf. This portion is so close to linear that, for the purposes of the calculations here, we shall treat it as a straight line, as shown by the dashed line. The temperature at the base is −13.4° C and the gradient is −0.023 deg m⁻¹. Now, in order to calculate the maximum amount of water freezing which could be caused by placing this cold base in sea-water we need to calculate the way in which temperature change is transmitted through the ice shelf. Consider a simple, constant-thickness conduction model in which the temperature at the base is suddenly altered from −13.4° C to −1.8° C, the temperature of the ice equilibrium with typical Antarctic sea-water (Reference WexlerWexler, 1960). Reference Carslaw and JaegerCarslaw and Jaeger (1959, p. 99) treat this general problem, and for our particular case their formula reduces to:

where T ₀ = −13.4° C, T _b= −1.8° C, and

the extrapolated surface temperature for our linear model. Figure 3 shows the temperature profiles for the first 100 years. Looking at the temperature profile after 25 years, which is the lime taken to reach the 3 km point, we see that the bottom 75 m of the ice shelf has undergone an average temperature rise of about 6 deg. The heat capacity of ice is 2.09 kJ kg⁻¹ deg⁻¹ (Reference HobbsHobbs, 1974, p. 362), so that this warming requires 8.15 × 10⁵ kJ of heat per unit area of ice shelf. Making the extreme assumption that all this heat is supplied through the freezing of sea-water onto the base we find, taking the latent heat of fusion of ice to be 333.6 kJ kg⁻¹ (Reference HobbsHobbs, 1974, p. 361), that it would result in the freezing of 2.7 m of ice over the total period of 25 years. This amounts to an average of only 0.11 m a⁻¹, though we expect the freezing to be more intense at the hinge line. This amount, which may be an overestimate, is not nearly large enough to account for the anomaly in the thickness profile. Unless some other effect, is causing intense freezing near the hinge zone we must therefore conclude that the thickness anomaly is a result of changes in the velocity field.

Fig. 3. Simple model of temperature changes as a land-based glacier becomes afloat. The solid lines show the temperature profiles after 0, 5, 25, and 100 years, plotted against height ξ above the bottom of the ice shelf. The dashed line shows the linear initial profile used to simplify the conduction calculation.

(i) General theory

Consider an ice shelf confined on two sides by parallel, vertical rock faces and flowing out to sea in the positive x direction (Fig. 4a). This is the model treated by Reference ThomasThomas (1973[b]). To derive an expression for the strain-rate at a point on the ice shelf we consider the combined forces due to sea-water pressure and to shear at the sides of the bay. The sea-water force is the same as in the simple unconfined model of Reference Weertman.Weertman (1957) used in Section II. A reasonable form for the force due to shear is derived by Thomas in the following way: The force per unit volume due to shear stress acting on an element at (x, y, z) is given by −∂σ _xy /∂y at the point in question, and therefore the total up-stream force due to shear along the section between x = x and x = X acting on a vertical element of unit width in the y direction is

The line x = X represents the limit of the influence of the bay sides on the ice shelf. Beyond this line the ice shelf can be considered to be unconfined.

Fig. 4. Plan views ice shelves in (a) a parallel-sided bay, (b) a bay with diverging sides, and (c) a bay with converging sides.

Thomas makes two further assumptions:

• (i) that ∂σ _xy /∂y is a constant across the ice shelf in the y direction, so that shear stress varies linearly. Shear stress is therefore given by

  (11)

  
  since it is zero along the centre line by symmetry considerations;

• (ii) that shear stress at the sides, averaged over z, reaches some limiting value

  
  independent of x. We expect
  
  to be of the order of 1 bar, although recent work (Reference HughesHughes, 1977) has suggested that it may fall to 0.4 bar if preferred crystal orientations are established by fast-moving ice. If λ is the half-width of the ice shelf then by Equation (11) we have at the sides

so that

We can therefore write

and Equation (10) becomes:

Now, Thomas’s general expression for the strain-rate in an ice shelf may be written

where

and is the total force due to sea-water pressure and shear which opposes the movement of a unit vertical section of ice shelf. θ is a number which takes account of the relative contributions of shear and longitudinal stresses to the effective stress term in the flow law; if α and β are defined through ⋵ _yy = α⋵ _xx and ⋵ _xy = β⋵ _xx then θ may be expressed as

For a bay with parallel sides there is no component of velocity in the y direction, and hence ⋵ _yy = 0 and α = 0. Further, along the centre line, ⋵ _xy = 0 so that β = 0; therefore θ = 2⁻ⁿ. On the centre line of a bay ice shelf the strain-rate is then given by

Our earlier formula, Equation (1), for strain-rates in a glacier tongue gave ⋵ _xx as a function only of thickness H; Equation (14) for a bay ice shelf gives ⋵ _xx as a function of thickness H and horizontal distance (X—x) from the ice front. In calculating the steady-state profile of an ice shelf we therefore need to set up boundary conditions at the ice front and integrate backwards towards the hinge zone. We cannot start at the hinge zone, since we should then need to know in advance the form of thickness over the whole distance between the hinge zone and the ice front—this is what we are trying to calculate.

For the case of an ice shelf restricted by parallel sides we have the condition that ⋵ _yy = 0, and hence the condition for steady state, Equation (4), becomes

for the centre tine of the ice shelf.

If the velocity of the ice shelf at x = X is u(X) then we can calculate increments of height in the negative x direction through the formula

where

with the assumption, as in Section II, that the mean density of the ice shelf is constant along its length.

In constructing the best possible model for a particular ice shelf we should consider all data available for such parameters as accumulation and melt, mean density, and mean temperature, and we should examine the way in which these vary along the length of the ice shelf. If the ice shelf does not have parallel boundaries we should take account of any transverse strain-rates due to the divergence of flow lines (Reference RobinRobin, 1958, p. 121). Initially we look at an ice shelf for which all these parameters are constant.

The ice shelf forms in a rectangular bay of width 2λ and length L. It is fed by a glacier at its landward margin, is subject to melt and accumulation, and discharges freely into the sea once it has escaped the influence of the bay. We have to choose the boundary conditions H(X) and u(X) at the seaward margin of the bay if we are to use Equation (15). This is inconvenient, as it is preferable to work in terms of conditions which arise at the hinge zone, such as the total mass input of the feeding glacier and the depth of bedrock where the glacier becomes afloat. We can however use these quantities if we consider the overall mass balance of the ice shelf. Let the volume input of the feeding glacier be M (in units of m³ s⁻¹ of ice at average density

), and let the net accumulation (a—m) be uniform over the entire area 2λL of the ice shelf. Then, to balance input and output we must have

Once reasonable values for volume input and net accumulation have been selected we can obtain a value for the product H(X)u(X). There is then just one degree of freedom in our choice of boundary conditions; choosing H(X) determines u(X).

We can use this last degree of freedom to make a choice of H(X) that will lead to an appropriate thickness of ice at the hinge zone, i.e. that thickness which will be just beginning to float when resting on the bedrock which is present at the hinge zone. The depth of bedrock at the hinge zone is thus the last remaining parameter we need to make the model determinate. For simplicity we shall consider the hinge zone to be in the form of an abrupt vertical edge at distance L from the seaward margin of the bay, submerged at a depth D below sea-level. We shall not consider here the possibility of instabilities in the position of the hinge zone (Reference WeertmanWeertman, 1974; Thomas and Bentley, 1978); such instabilities depend on the dynamics of the feeding glacier and on the slope of the bedrock, and are beyond the scope of this paper.

(ii) Behaviour of the model

It is found, when computing models of ice shelves, that the thickness calculated at the hinge zone is extremely sensitive to the values chosen for H(X) and u(X) at the seaward margin, but that, by making fine adjustments, the hinge-zone conditions can be exactly satisfied. It may seem unreasonable that an adjustment of a fraction of a metre in the value taken for H(X) can lead to differences of 100 m and more in the value the model gives for thickness at the hinge zone. However, this peculiarity is only a result of the process of working backwards from the seaward margin, and does not reflect a physical instability. Quite the contrary, it merely means that large changes in the thickness and velocity at the hinge zone lead to relatively small changes at the seaward margin, provided that the total volume input is not altered. This is a physically reasonable result.

We shall look at the profile of an ice shelf in a bay 100 km wide and 150 km in length. We shall take the volume input to be 1.2 × 10¹⁰ m³ a⁻¹, which corresponds, for instance, to an ice shelf 600 m thick moving at 200 m a⁻¹ at the hinge zone. For flow and shear stress parameters we shall adopt the values found by Thomas (1973[a]) in his analysis of data from Amery Ice Shelf (Reference BuddBudd, 1966). He found that

and we shall use the values

Three cases will be considered: an ice shelf undergoing zero net accumulation, one undergoing net accumulation of 0.5 m a⁻¹, and one undergoing net melting of 0.5 ma⁻¹. The model need make no distinction between the mass-balance contributions of snow-fall or ablation at the surface and freezing or melting at the bottom.

Figure 5 is intended to illustrate the general behaviour of the models. The essential feature of a bay ice shelf is the delicate balance between the driving force of weight of ice above sea-level and the restraining force due to shear at the bay sides. These two forces are represented by the two opposing terms in Equation (16):

which we shall call the “Weertman term”, and

which we shall call the “shear term” (

is negative). The choice we make in the initial conditions H(X) and u(X) is critical in any determination of which term dominates. If we start the model with too large a thickness at the seaward margin we find that the Weertman term dominates and increases rapidly; the ice shelf has large increasing extending strain-rates which result in rapid thickening as the model moves towards the hinge zone. If the chosen initial thickness is much too large then the ice shelf diverges to infinite thickness before the hinge zone is reached. If, on the other hand, we choose too small an initial thickness then the shear term dominates; strain-rates decrease inwards towards the hinge zone and can eventually become compressive. The thickness of the ice shelf decreases and may vanish before the hinge zone is reached. There are solutions of this kind which do not vanish before the hinge zone, but they give a very thin, fast-moving ice shelf at the hinge zone, and compressive strain over much of the ice-shelf length.We shall ignore these solutions as unrealistic, since most known ice shelves are thick and slow-moving at the hinge zone. In order to model conditions which arise in practice, we vary conditions at the seaward margin within the conditions prescribed by the general mass-balance formulation of Equation (17) until the model gives an appropriate thickness at the hinge zone. We shall take this thickness to be 600 m, which corresponds to hinge-zone bedrock at a depth 496 m.

Fig. 5. Behaviour of ice-shelf models for a parallel-sided bay with zero net accumulation: (a) ice thicknesses, (b) velocities, (c) strain-rates, and (d) the opposing terms, the Weertman term, and the shear term. The shear term is plotted as a positive quantity. The shorter dashes show the behaviour of a model with too great an initial thickness, the longer dashes one with too small an initial thickness.

Table I for a parallel-sided bay shows the values calculated for H(x)u(X) using Equation (17), and the values of H(X) and u(X) which were found to satisfy the hinge-line condition. They are shown for three cases: net accumulation of 0.5 m a⁻¹, which might correspond to 0.2 m a⁻¹ surface accumulation and 0.3 m a⁻¹ bottom freezing; zero net accumulation, which corresponds to surface accumulation exactly balanced by bottom melting; and net ablation of 0.5 m a⁻¹, which might correspond to 0.2 m a⁻¹ surface accumulation and 0.7 m a⁻¹ bottom melting. The thickness profiles are shown in Figure 5 (zero accumulation) and Figure 6 (net melt and freeze), together with the velocity, strain-rate, and the Weertman and shear terms. As expected, an ice shelf subject to melting is thinner and slower at the ice front than an ice shelf subject to freezing. However, it is an interesting feature of the solutions that the thickness in all three cases appears to change approximately linearly as we move back from the ice front, at least for some 100 km. Furthermore, it appears that the thickness gradient is very nearly the same for the three cases, even though the conditions of thickness, velocity, and melt are radically different. Changes in thickness gradient occur almost entirely within the last 50 km before the hinge zone: in the model with melting we find large slopes at the hinge zone, while with freezing the ice shelf is almost flat. Near the ice front the thickness gradients supplied by the model are: freeze, 2.08 × 10⁻³; zero accumulation, 2.06 × 10⁻³; melt, 2.23 × 10⁻³. These show a spread of only eight, per cent, whereas thicknesses in the three cases differ by a factor of two, velocities by a factor of three and strain-rates by an order of magnitude.

Fig. 6. Behaviour of ice-shelf models for a parallel-sided bay with 0.5 m a⁻¹ melt (dashed line) and with 0.5 m a⁻¹ freezing (solid line).

This happens because, if we refer to the graphs of the shear term (Figs 5(d) and 6(d)), we see that in the first 100 km from the ice front this term has very nearly the same form for the three different cases. That is, the restraining force due to shear at the sides seems to be practically independent of the thickness of the ice shelf. This can be shown mathematically. Consider an ice shelf of thickness H(X) at the ice front with a small thickness gradient ζ, so that thickness as a function of distance is

Then the shear term is given by

Now

can be treated as a small quantity for distances within 50 km from the ice front and we can expand binomially to get

This means that the restraining force due to shear increases roughly linearly with distance in from the ice front, and has a rate of increase of the order of, or slightly less than,

. This quantity is clearly independent of thickness, velocity, or melt-rate. Our parameters give it the value of 6.47 × 10⁻⁹ s^−⅓ m⁻¹ which compares with values from our models of: freeze, 6.21 × 10⁻⁹ s^−⅓ m ⁻¹; zero accumulation, 6.18 × 10⁻⁹s^−⅓m⁻¹; and melt, 5.98 × 10⁻⁹ s^−⅓ m⁻¹. This agreement is satisfactory.

Let us see what this means for the dynamics of a stable ice shelf. Shear force increases steadily, and almost linearly, as we move inwards from the ice front and has to be counterbalanced by the driving force due to the weight of ice above sea-level. It may be seen from the graphs of these two opposing forces that in order for a large stable ice shelf to exist they must just keep pace with each other. If the Weertman term begins to dominate then the ice shelf diverges to infinite thickness before the hinge zone; if the shear force dominates the ice shelf converges and vanishes. The difference between the two forces must remain approximately constant. Thus, in the three models we find that the Weertman term has quite different absolute values and yet approximately the same rate of increase with distance. In our models the Weertman term increases at the following rates: freeze, 3.70 × 10⁻⁹ s^−1/3 m⁻¹; zero accumulation, 3.66 × 10⁻⁹ s^−1/3 m⁻¹ and melt, 3.84 × 10⁻⁹ s^−1/3m⁻¹. These are some 40% less than the rates of increase of the shear parameter given above, a feature which is reflected in the fact that the strain-rate (which is given by the cube of the difference of the two quantities) decreases as we travel inland. Nonetheless, the quantities are nearly enough equal that we can use this as the basis for a calculation of an expected value for the ice-shelf gradient.

Let us assert that the rates of increase of the Weertman term (Equation (20)) and the shear term (Equation (21)) must be approximately equal for a large, stable ice shelf to exist. The rate of increase of the Weertman term is

and the rate of increase of the shear term is

. Equating these leads to

This gives the expected thickness gradient of an ice shelf in terms of its density, its half-width, and the shear stress at the sides. In this approximation the thickness gradient is independent of melt or accumulation, flow-law parameters, thickness, and velocity. We only expect the formula to be approximate as we have equated two quantities which, in our modelled cases, differ by 40%; indeed, Equation (22) predicts a thickness gradient of 3.7 × 10⁻³ for an ice shelf with the parameters used in our models, which is greater than the modelled values by 75%. We should therefore trust the formula only to within a factor of about two.

(iii) Agreement with reality

The important feature of Equation (22) is its inverse dependence on half-width. The factor

is approximately constant for different ice shelves and has a value of 1.85 × 10² m when the parameters of equalities (18) and (19) are used. Let us see how the rule agrees with available data. Figure 7 shows thickness gradient plotted against half-width for various bay ice shelves. The error bars shown for each set of data do not represent the experimental errors made in determining the quantities shown, they represent the range of values typically found for the quantities on each ice shelf. Thickness gradients have been taken for the middle region of the ice shelves, well away from the bay walls, the hinge zone, and the ice front, and the Ross Ice Shelf has been considered to be a single ice shelf, even though it is dominated by major ice streams (Reference RobinRobin, 1975). These streams give rise to a large scatter in the thickness gradients. The solid line in the graph shows the theoretical curve based on Equation (22), the agreement is satisfactory, both in form and in magnitude. The curve has been drawn using the data of equalities (18) and (19), yet it coincides almost exactly with a hyperbola fitted by regression.

Fig. 7. Relationship between thickness gradient and half-width for ice shelves. The solid line shows the theoretical curve based on Equation (22). Data are plottedfor eight ice shelves: Bach Ice Shelf (J . F. Bishop, personal communication); Jelbartisen (Reference Autenboer and DecleirAutenboer and Decleir, 1975) ; George VI Ice Shelf, southern end (J. F. Bishop, personal communication ) ; Maudheim Is-shelf (Swithinbank, 1957); Amery Ice Shelf (Thomas, 1973[a) ); Filchner Ice Shelf (Scott Polar Research Institute radio-echo flights ) ; Ronne Ice Shelf (Reference SwithinbankSwithinbank, 1977) ; and Ross Ice Shelf (Reference RobinRobin, 1975). The bars show the range of values typically found on each ice shelf.

The slightly anomalous results for the Maudheim, Jelbartisen, and Ronne ice shelves can be explained in a number of ways: a different mean temperature through the ice shelf would very likely lead to a different value for the limiting shear

at the sides; different density parameters would alter the density term in the solution; and, more importantly, not all the ice shelves lie in bays which have almost parallel sides. Maudheim Is-shelf clearly converges while the Jelbartisen, and indeed the Amery Ice Shelf, diverge. Nevertheless, all the ice shelves studied show behaviour which lies within a factor of two of that predicted.

(i) Theory

Consider now an ice shelf which has formed in a bay with diverging, though still vertical, bay walls. Each wall diverges at an angle ψ (Fig. 4b). The ice shelf now widens as it flows out to sea since it is no longer restrained in the y direction. If the divergence of the bay is sufficiently gentle, then the ice shelf will spread sufficiently rapidly in the transverse direction to fill the entire width of the bay. This means that, for an ice shelf moving forward with an average velocity u , the average strain-rate over width in the y direction must be

Reference BuddBudd (1966, p. 352) shows theoretically that, with the flow parameter n = 3, shear strain is concentrated towards the sides of an ice shelf; the effect is observed to be even more pronounced (Reference SwithinbankSwithinbank, 1963; Reference HughesHughes, 1977). Velocity changes little over the width of the shelf, so the average velocity is close to that of the centre line. Here I shall consider the average velocity to be equal to that calculated for the centre line. The mass-balance equation (4) then becomes:

A more precise definition of “sufficiently gentle divergence” can be given by looking at the maximum possible rate at which an ice shelf can creep in the transverse direction. For a given thickness of ice shelf this maximum rate occurs for an unconfined glacier tongue, and is, by Equation (9),

The maximum permissible divergence of the ice shelf sides ψ _max is then given by combining Equations (23) and (25) to give

The substitution of values from equalities (18) and (19) indicates that

This equation shows that as an ice shelf flows in a diverging bay, becoming both faster and thinner towards the front, there will come a point where the ice shelf is unable to sustain sufficient transverse creep for it to remain intact. It should then either come adrift from the bay sides or develop crevasses and areas of weakness. It is then likely to break up. The most important factor in Equation (26) is the thickness term. The thickness of an ice shelf typically varies by a factor of about three between the hinge zone and the ice front so that H ³ is expected to vary by a factor of about 30. The velocity and the half-width vary by factors of only two or three, and both increase towards the ice front, so that their ratio remains relatively stable. For an ice shelf of half-width 50 km and velocity 300 m a⁻¹ we can therefore write

which yields the critical thickness at which a typical ice shelf must come adrift from bay walls of a certain angle of divergence Ψ. Table II shows the relationship. If an ice shelf of these assumed dimensions exists in a bay with sides diverging at 10° or 20° we expect it to become free from the sides once it has thinned to a thickness of 200–250 m. This might explain why the Amery Ice Shelf, for instance, extends no further than it does in its diverging bay; at the front its thickness is 270 m (Thomas, 1973[a]), its velocity is about 1 250 m a⁻¹ (Reference BuddBudd, 1966) and the half-width is about 100 km. According to Equation (26) this means that the maximum permissible divergence before the ice shelf comes adrift is is ψ _max = 13.4°. The actual divergence at the front is about 16°, increasing to 26° (angles taken from the plan view in Thomas (1973[a], p. 63)).

In order to derive an expression for the equilibrium profile of an ice shelf in a diverging bay we now need to use Equation (24) with a suitable expression for ⋵ _xx ; Equation (12) gives the general form of this expression. In order to use it we need to know a value for α in Equation (13), which expresses the effect of transverse stress deviators on the flow law; α is defined by ⋵ _yy = α⋵ _xx , and is therefore important only when the transverse strain-rate becomes comparable with the strain-rate in the x direction. For n = 3, we find that α = 1 gives θ = 0. 111, α = ½ gives θ = 0.112, and α = 0 gives θ = 0.125. For most of the ice shelf we shall find that ⋵ _yy « ⋵ _xx , so that if we take θ = 0.125, as in Equation (14), we shall be risking errors of at most ten per cent.

An exact treatment is much more complicated and not worthwhile since other parameters, such as m and B , are not known to better than ten per cent. In what follows we shall therefore use Equation (14) as before with one minor modification: because flow lines are diverging and, at the edges of the ice shelf, parallel to the bay walls, the shear stress

no longer acts in the same direction as the centre line of the ice shelf. Along the centre line the component of shear drag is cos Ψ which, for small angles, leads to very little change. For the equilibrium profile we therefore have:

where

The quantity λ in the shear integral is now a variable with the form

where λ(o) is the half-width of the ice shelf at the hinge zone.

(iii) Behaviour of the model

We retain the input conditions of our previous model, a volume input of 1.2 × 10¹⁰ m³ a⁻¹ over a hinge zone which is 100 km wide with ice 600 m deep. We retain the values previously used for

, B , , and d, but we now have a divergence angle Ψ of 15°. At the front the bay walls are therefore 180.4 km apart. Our equation for general mass balance becomes

Table I shows values of H(X)u(X) calculated for a diverging bay using this equation and the values of H(X) and u(X) which satisfy the hinge-line conditions, for the three cases of net freezing, zero net accumulation, and net melting. We can, using Equation (26), calculate the maximum permissible bay divergence for these three models if the ice shelf is to remain intact, it is 25° for net freezing, 14° for zero net accumulation, and 0.01° for net melting. Since the actual bay divergence is 15° this means that with net freezing the ice shelf is stable, with zero net accumulation it just becomes unstable as it nears the seaward margin, and with net melting it is quite unstable and the ice shelf would come adrift long before reaching the seaward margin of the bay. The model treatment is therefore not valid in this case, since it assumes contact with the bay walls over the entire length of the bay.

We conclude that an ice shelf which is undergoing melt cannot exist in the entire length of the bay for the assumed input conditions. A smaller ice shelf can form, but it terminates before the seaward margin of the bay is reached. By a lengthy process of trial and error we find that an ice shelf extending 38 km from the hinge zone can form and remain in contact with the bay walls. At the 38 km point it has a thickness of 220.7 m and a velocity of 366.8 m a⁻¹, for which the critical bay divergence is 15°. With further thinning, the ice shelf can no longer cope with the divergence of the bay and so comes adrift. The form of the solution is shown in Figure 8.

Fig. 8. Behaviour of ice-shelf models in a bay with sides diverging at 15°, for 0.5 m a⁻¹ melt (longer dashes); net balance of melt and accumulation (solid line); and 0.5 m a⁻¹ freezing (shorter dashes ). The melting ice shelf can extend 38 km before it comes adrift from the bay walls.