## Abstract

Membrane stresses act along thin bodies which are relatively well lubricated on both surfaces. They operate in ice sheets because the bottom is either sliding, or is much less viscous than the top owing to stress and heat softening of the basal ice. Ice streams flow over very well lubricated beds, and are restrained at their sides. The ideal of the perfectly slippery bed is considered in this paper, and the propagation of mechanical effects along an ice stream considered by applying spatially varying horizontal body forces. Propagation distances depend sensitively on the rheological index, and can be very large for ice-type rheologies.

A new analytical solution for ice-shelf profiles and grounded tractionless stream profiles is presented, which show blow up of the profile in a finite distance upstream at locations where the flux is non-zero. This is a feature of an earlier analytical solution for a floating shelf.

The length scale of decay of membrane stresses from the grounding line is investigated through scale analysis. In ice sheets, such effects decay over distances of several tens of kilometres, creating a vertical boundary layer between sheet flow and shelf flow, where membrane stresses adjust. Bounded, physically reasonable steady surface profiles only exist conditionally in this boundary layer. Where bounded steady profiles exist, adjacent profile equilibria for the whole ice sheet corresponding to different grounded areas occur (neutral equilibrium). If no solution in the boundary layer can exist, the ice-sheet profile must change.

The conditions for existence can be written in terms of whether the basal rate factor (sliding or internal deformation) is too large to permit a steady solution. The critical value depends extremely sensitively on ice velocity and the back stress applied at the grounding line. High ice velocity and high stress both favour the existence of solutions and stability. Changes in these parameters can cause the steady solution existence criterion to be traversed, and the ice-sheet dynamics to change.

A finite difference model which represents both neutral equilibrium and the dynamical transition is presented, and preliminary investigations into its numerical sensitivity are carried out. Evidence for the existence of a long wavelength instability is presented through the solution of a numerical eigenproblem, which will hamper predictability.

## 1. What are membrane stresses and why are they important?

When a thin body experiences low traction conditions on its two larger surfaces, stresses tend to direct along the body. These are known as membrane stresses. Such stresses can be of significance in ice shelves and glaciers, where they act essentially in the two horizontal planes (HP). This is due to two special factors. A thin Newtonian body attached to its base, but with a free surface, does not have such strong membrane effects as a glacier. This is because the glacier ice has a nonlinear rheology, and may be slipping over its bed. Both of these factors increase the significance of membrane stresses in thin body mechanics. Glaciologists have considered the role of membrane stresses mostly in plane flow, for which reason they are usually referred to as ‘longitudinal stresses’. This is ambiguous in three-dimensions, and for this reason the term ‘membrane stresses’ is used in this paper. They can be viewed as the three-dimensional equivalent of longitudinal stresses.

Strictly, membrane stresses are those which occur in membranes with tractionless surfaces, with minimal shear stresses acting in planes parallel to the surface, with the implication of low shear-strain rates. However, when the role of the horizontal stresses in glaciers becomes significant their role is very similar to membrane stresses, and it is therefore convenient and appropriate to call these horizontal stresses ‘membrane stresses’. Where we need to be strict, where there is no tangential traction on the bed, we will call the horizontal deviator stress fields ‘pure membrane stresses’ and other situations ‘quasi-membrane stresses’.

The role of membrane stresses in deglaciation of marine ice sheets has been emphasized by glaciologists for several decades (e.g. Weertman 1974). Membrane stresses play an important role in glaciological fluid dynamics as a consequence of the nonlinear rheology of ice and the presence of sliding at the glacier bed. Scenarios have been presented whereby ice shelves transmit membrane stresses long distances up sliding ice streams, ‘buttressing’ the ice sheet and preventing rapid outflow of grounded ice to the sea (Hughes 1973). The implication is that removal of ice shelves will result in large increases of ice flow into the ocean, with significant effects on global sea level.

Despite the widespread currency of these conceptual models, there is little theoretical backing for them in the form of quantitative models, and indeed the best-known Antarctic model (Huybrechts 2002) apparently relies on an entirely different mechanism (Hindmarsh & Le Meur 2001) to simulate the Holocene deglaciation of West Antarctica. However, there is a lot of evidence for changes in ice-sheet dynamics associated with shelf removal.

Even the controlling qualitative dynamics of marine ice sheets remain unknown, with different authors suggesting conditionally no equilibrium positions (Weertman 1974), unique equilibria depending on climate (Oerlemans 2002), or neutral equilibrium (Hindmarsh 1993). Such knowledge is fundamental in designing appropriate numerical models.

This paper reviews the role of membrane stresses in various issues of ice-sheet stability. The transmission of membrane stresses up perfectly slippery laterally resisted streams is considered. As has been speculated, the distance of transmission of forces is governed by the ice-stream width (Gudmundsson 2003), and scaling analyses agree with this (Hindmarsh 2006). In agreement with some recent models (Schmeltz *et al*. 2002; Sugiyama *et al*. 2003; Payne *et al*. 2004; Thomas 2004; Dupont & Alley 2005; Schoof, 2006*a*,*b*), it is shown that the initial response to changes in forces is fast, but that the total response is not large. In other words, the rapid changes associated with the transmission of membrane stresses are not necessarily sustained for long periods. There are many observations of rapid change (Rignot 1998, 2002; Wingham *et al*. 1998; Shepherd *et al*. 2001–,2003; Joughin & Tulaczyk 2002; Rignot *et al*. 2002, 2004; De Angelis & Skvarca 2003; Joughin *et al*. 2003; Scambos *et al*. 2004) and these have been argued to be caused by downstream changes in abutting ice shelves. It is not yet clear that we have adequate theoretical knowledge to make such claims, and one of the aims of this paper is to establish some basic results about the propagation of effects up an ice stream. In particular, near the boundary grounding line a membrane stress boundary layer a few tens of kilometres in length can be shown to exist, where changes in membrane stresses transmitted across the grounding line affect the dynamics. This idea has been developed independently by C. Schoof (2005, personal communication). These theoretical effects are consistent with recent observations of change and are not necessarily indicative of instability.

However, the problem of how the West Antarctic deglaciated after the Last Glacial Maximum remains (e.g. Conway *et al*. 1999). The influence of longitudinal stresses on the marine ice-sheet equilibria is examined, in particular in the case of neutral equilibria, where it is shown that the presence of membrane stresses imposes a condition on the neutral equilibrium property and the existence of physically reasonable steady solutions. This has fundamental significance for marine ice-sheet dynamics, as if the boundary conditions do not permit a steady solution, the ice sheet must change. Numerical evidence is presented that the dynamics where no steady solution exists lead to large-scale change rather than, for example, small oscillations.

## 2. The Stokes equations

In this section, the basic mechanical field equations and boundary conditions are presented. Three flow configurations are considered: (i) cases where an ice sheet in plane flow is symmetric around the divide and terminates at the grounding line, where a resisting force is provided either by water or by an ice shelf; or (ii) flow down the infinite plane, with the ice fixed to the bed or sliding over it with non-zero traction; and (iii) flow down an infinitely long ice stream with a perfectly slippery bed, with the ice fixed to the sides. In the latter two cases, flow is driven by a horizontal body force component.

The coordinates are , where *z* is perpendicular to the base plane and *x* is in the zeroth-order flow direction. The *z*-direction is called ‘vertical’ and the plane ‘horizontal’. The upper and lower surfaces are given by , , respectively, the thickness of the ice is given by and *t* represents time. Sea level (constant in this paper) is taken to be at *z*=0. Subscripts (b) and (s) indicate evaluation at the surface or base, respectively. The operators represent the horizontal gradient and divergence, respectively.

The three-dimensional velocity field is conveniently represented by the vertical velocity *w* and the horizontal velocity vector , and we also use . The governing equations, which apply to all * r*, are(2.1)(2.2)(2.3)(2.4)Here, equation (2.1) expresses conservation of mass in the ice; equations (2.2)–(2.4) describe conservation of momentum in the ice, where σ is the stress tensor,

*ρ*is the density of ice, is the gravitational acceleration vector and

*is the outward normal vector at the indicated surface. The horizontal component of gravity is functionally equivalent to a slope of . The constitutive relations comprise of (i) a nonlinear viscous relationship within the ice(2.5)where e is the strain-rate tensor,*

**n***τ*is a second invariant of the deviator stress tensor τ,

*n*is the Glen index and

*A*

_{c}is a rate factor, and (ii) an isotropic sliding relation of the form(2.6)where is the sliding velocity, is the basal tangential traction, is the sliding index, is the sliding rate factor, is the effective pressure, is the normal traction, is the sub-glacial water pressure and

*ν*is a further index. The upper and lower surface kinematical conditions are(2.7)

In some cases (§5), the flow of an infinitely long section is considered. It is necessary in general to consider the effect of the non-zero longitudinal strain rate has on the flow, especially on the viscosity. The assumption of quasi-uniform flow states that such a flow , where superscript QU indicates a quasi-uniform flow, is given bywhere is treated as being independent of ** r**, but that may be non-zero. Essentially, it is being assumed that over one wavelength is sufficiently small that the

*-independence of is a valid assumption. This assumption is applied to the base flow about which linearization is performed. In the flows along infinite sections, relevant field variables are linearized about a base case solution (steady uniform flow), for example, (2.8)and so on. These are used to derive a set of zeroth- and first-order equations expressing conservation of mass and momentum and a Fourier transform in the HP is then applied to the first-order equations.*

**r**In other cases, an ice sheet with finite span is considered. At the divide, one applies an even or odd symmetry condition on the velocity and stress fields as appropriate, while at the grounding line, we consider that arbitrary forces can be applied by either the water or an abutting ice shelf onto the grounded ice. This boundary condition will only be used in vertically integrated (force balance models) in this paper, so attention is restricted to force balance. The force exerted on the grounded ice is given by(2.9)where is a normal stress. Subscript G indicates evaluation at the grounding line. Letting the vertical normal stress be approximated by , which is true for both types of mechanical approximation needed in this paper, the usual formula for the mean horizontal deviator stress is found to be(2.10)in the unconfined case.

## 3. Shallow scaling showing membrane stresses

In this section, a scaling of the Stokes equations is briefly presented, which shows how the influence of membrane stresses depends upon the basic traction. These results are known and based on one presented by Hindmarsh (1993). This scaling is used in this section only. The scaling is fairly similar to most shallow ice approximation (SIA) scalings (Hutter 1983; Morland 1984; Fowler 1992). Variables without a tilde are taken to be dimensionless counterparts of the physical quantities defined above. Letting subscript * represent a scale magnitude, we scale , vertical coordinates by and horizontal coordinates by . A small Stokes number *δ* may be defined by(3.1)Following Hindmarsh (1993), one may scale the shear stress by and the (basal) tangential traction . The quantity *λ* is termed as the *traction number*, and it represents the ratio of the basal shear stress to the surface longitudinal stress. The choice is made, which implies that(3.2)The horizontal velocity scale is .

In two-dimensions, the momentum balance equation (2.2) can be manipulated to eliminate the hydrostatic pressure to produce one integro-differential equation. In scaled form the momentum balance equation can be written(3.3)and Hindmarsh (1993) shows that provided , the Stokes equations may be rewritten as(3.4)where the ‘bridging terms’ of (3.3) have been replaced by their gauge function. Hindmarsh also shows that the boundary conditions provide bounds to the traction number *λ*, such that , the lower bound corresponding to perfect slip, the upper bound to a fixed bed. These bounds are known from shallow theory of ice shelf and ice-sheet flow, respectively (e.g. Hutter 1983; Weiss *et al.* 1999). It can also be shown (e.g. Morland 1984, where the same notation is used) that(3.5)and the momentum balance equations (3.4) may therefore be written as(3.6)This is the plane flow (quasi-)membrane stress approximation.

The main point to be derived from this scaling is how the influence of the membrane stress term (second term on the LHS) depends upon the basal traction. When the traction is high (), the membrane stress term is small, while when the traction is low (), the membrane stress is large. Similarly, *ϵ* can be regarded as a wavenumber, showing that membrane stress terms become larger at shorter wavelengths.

## 4. Stokes equation approximants and their accuracy

The Stokes equations are rarely solved in the form above, and usually some kind of approximation procedure is used to simplify the equations, enhancing analytical and computational tractability. The purpose of this section is to list the most widely used approximations, so that their relationships are clear and the accuracy of the approximations used later on in this paper are apparent.

The approximants all in some way start from the SIA (Hutter 1983), and, like this approximation, involve dropping terms from the momentum balance equations and simplifying the strain-rate definitions. However, membrane stress approximants contain more terms than the SIA under the supposition that this decreases the error of the approximation. This has been discussed in §3.

The SIA shows that for thin flows where sliding is small, the vertical gradient of the HP shear stresses balances the horizontal pressure gradients. The accuracy of the SIA decreases as the amount of basal slip increases (Gudmundsson 2003). The shallow theory for ice-shelf spreading, where the proportion of slip is very high, balances the vertical gradient of the HP shear stresses with the horizontal gradients of the longitudinal stresses as well as the pressure gradient. These gradients must play a role in well-lubricated grounded stream flows (Muszynski & Birchfield 1987; MacAyeal 1989). This can be represented through one of two equivalent parameters, the traction number (the ratio of basal tangential traction to the surface stress deviator) or the slip ratio (the ratio of the sliding velocity to the difference between surface and basal velocity). Hindmarsh (1993) and Wilchinsky & Chugonov (2000) have shown that the approximation error depends on the degree of slip and also shows that a shallow long wavelength theory suffices for all possible slip ratios.

Membrane stresses can be included with a vertically integrated model Hindmarsh (2004) has produced a single-layer scheme which is formally of the same order of error as the scheme of Blatter (1995), who constructed a scaling which permits free vertical variation of shear and membrane stresses. Such ‘longitudinal stress schemes’ are used widely in ice-sheet modelling (Hubbard *et al*. 1998; Ritz *et al*. 2001; Huybrechts 2002; Pattyn 2003; Saito *et al*. 2003; Payne *et al*. 2004; Dupont & Alley 2005). Quantitative comparisons show that some vertically integrated schemes are of comparable accuracy to the three-dimensional schemes (Hindmarsh 2004).

Reminding readers that in this section we have now returned to a dimensional system of equations, we now define several Stokes equations approximants. The continuity equation is satisfied in all the approximations. The following notation is used(4.1)

(4.2)

(4.3)

### (a) Lubrication theory approximations

Vertical HP shear stress gradients are balanced by the pressure gradients and the downstream body force. Thus, the field equations are approximated by(4.4)where *p* is the pressure, with boundary tractions given by(4.5)

### (b) Membrane stress approximations

Horizontal gradients of membrane stress components are incorporated into momentum balance equations in this class. The main differences between the class members are the way in which these stresses are approximated. The field equations are approximated by(4.6)with boundary tractions given by(4.7)

#### (i) Vertically integrated membrane stress approximation (1) using at surface

Here, the surface velocities used in computing the membrane stresses are computed using the shear stresses in the shear strain relationship and in the sliding relationship. If one excludes the within-ice deformational component of velocity this is identical to the approximation used by MacAyeal (1989). The momentum balance equations are thus approximated by (4.6) and (4.7), and the stretching relationships use the surface strain rates to compute the membrane stresses(4.8)This approximation is used in §§5 and 6.

#### (ii) Vertically integrated membrane stress approximation (2) using at surface, vertical correction of to account for stress softening

This approximant is not used directly in calculations in this paper, but is referred to in §4*a* to assess the accuracy of vertically integrated membrane stress models. Surface velocities used in computing the non-HP stresses are computed using the shear stresses in the shear strain relationship and in the sliding relationship. The momentum balance equations are thus approximated by (4.6) and (4.7), and the stretching relationships are(4.9)where is given by its SIA approximant, and for the HP shearing(4.10)where is now given by its membrane stress approximant.

#### (iii) Free vertical variation of membrane stresses

This is the classic ‘longitudinal stress’ scheme as used by Blatter (1995), Pattyn (2003) and Payne *et al*. (2004). The momentum balance equations are thus approximated by (4.6) and (4.7), and the stretching relationships are(4.11)Compared with the vertically integrated schemes, the membrane stresses use the strain rate at the corresponding elevations rather than at the surface, and that the stress invariant calculations use local strain-rate values. This increase in accuracy comes at the expense of turning a two-dimensional computational problem into a three-dimensional one.

Hindmarsh (2004) carries out a quantitative comparison of the accuracy, showing that formally the vertically integrated membrane stress approximant (2) has the same accuracy as the permitting free vertical variation, and in detailed quantitative comparisons performs very well. Hindmarsh (2004) also shows that the vertically integrated membrane stress approximant (1) is somewhat less accurate, but has the right qualitative features. In particular, it is a legitimate asymptotic expansions in *ϵ* for the range of traction numbers (bed slipperiness) encountered in glaciology. This approximant is used in a number of models (e.g. MacAyeal 1989; Schmeltz *et al*. 2002; Joughin *et al*. 2004; Dupont & Alley 2005) as well as some calculations in this paper where the plane and radial flow of marine ice sheets are considered.

The thrust of the argument is that since the vertically integrated membrane stress approximations are qualitatively and quantitatively similar to three-dimensional membrane stress approximations, the analyses which follow are not artefacts of the vertical integration procedure.

### (c) The stream equations

In the case of perfect slip, all the membrane stress approximants become identical, and vertically integrated models can be used without loss of accuracy.

The stream equations are(4.12)

(4.13)

## 5. The perfectly slippery stream

We now turn out attention to how far forces at the grounding line can be transmitted along glaciers, motivated by the intuition that the more membrane-like the configuration (e.g. the more slippery the bed), the further the forces will be transmitted. In this section and the following, flow along an infinite stream with semi-width *W* is considered using the perturbation (2.8), and apply boundary conditions . We consider perfectly slippery streams, where . Here, the membrane stress approximation is accurate (MacAyeal 1989). Although the rate of surface lowering in, for example, Pine Island Glacier is large compared with the accumulation rate, the longitudinal strain rate associated with this lowering is around 0.001, compared with transverse vertical plane strain rates of around 0.1 (Payne *et al*. 2004). As will be seen, this means that the perturbation is still ‘small’ in a mathematical sense, and offers encouragement that even rapidly changing ice streams can be treated as linear systems, which is of immense use in data assimilation algorithms (e.g. Arthern & Hindmarsh 2003).

We scale the equation set (2.5), (4.12) and (4.13) as follows; thicknesses *H* and elevations are scaled by , a representative thickness; horizontal distances by the ice stream semi-width ; and other quantities are scaled as follows:(5.1a)(5.1b)(5.1c)(5.1d)to get(5.2a)(5.2b)and the velocity scale is(5.3)Hindmarsh (2006) considers the case where a series expansion in the width to length ratio is obtained. As he shows, and as has been shown many times before (e.g. Echelmeyer *et al*. 1993; Raymond 2000), when downstream stress gradients are small(5.4)which serves here as a base case for a linearization. For this, we expand(5.5a)(5.5b)(5.5c)(5.5d)(5.5e)(5.5f)(5.5g)to get at zeroth-order(5.6a)(5.6b)and at first-order(5.7a)(5.7b)Use of the fact that at zeroth-order, downstream changes are zero, which gives(5.8a)(5.8b)and defining(5.9a)(5.9b)(5.9c)gives the first-order equations momentum balance equations(5.10a)(5.10b)Applying a Fourier transform(5.11a)(5.11b)(5.11c)(5.11d)gives(5.12a)(5.12b)The linear rheology statement is(5.13)while a nonlinear rheology relates stresses to strain rates as(5.14a)(5.14b)The first-order continuity equation is(5.15)There is some dispute over the correct boundary condition to apply to the continuity equation. We shall not apply one on the basis that there is no flow into this stream from the flanks. In other words, is free to vary on the boundary, in contrast to Schoof (2006*a*), who fixes *H*.

In Fourier space, the continuity equation is(5.16)At the boundary, we have(5.17)and the continuity equation becomes(5.18)In steady state, we have(5.19a)(5.19b)(5.19c)At the boundary, we have(5.20)which implies(5.21)

In these cases where the flow is down the infinite plane, the force at the grounding line is not directly representable, and it is simulated by including a non-uniform horizontal body force component. Forces are applied by constructing a Gaussian bump in the perturbed body force. The force is also Gaussian in the transverse direction, as applying a transversely uniform body force has, for reasons yet to be explored, an unduly large influence on the flow near the lateral margins of the stream. The body force is given a distribution(5.22)where , are bump length parameters. Sincethe total force is(5.23)

Some results are shown in figures 1–3, for different ice rheological indices. The depth of the ice is 2 km, the slope *ϵ* is 0.001 and the flow speed is around 1.8 km yr^{−1}, which are roughly consistent with Pine Island Glacier (Payne *et al*. 2004). The rate factor *A* has been tuned to obtain these values, and should not be expected to be the same is in Payne *et al*. (2004) as that study included basal friction. It can be seen that rheological index has a dramatic influence on the distribution of stress/velocity effects along the trunk of the stream (right-hand panels), which in turn affects the distribution of final elevation along the stream. Note that, in common with Payne *et al*. (2004), that the strain-rate/thickness feedback causes the final effects to be more spatially extensive than the initial response.

With the nonlinear rheologies, the viscosity of the ice varies markedly across the width of the ice stream. This means that there is a central spine of viscous ice in the low shear zone, analogous in many ways to the more viscous ice near the upper surface of a glacier. Stress/velocity effects are propagated much further along this spine, meaning that perturbed strain rates are lower and that therefore the final response is less. This effect has been mitigated against by the choice of lateral variation of the forcing, but is very evident in the *n*=4 case.

## 6. Analyical solutions for steady ice shelf profiles in plane flow

Much insight into the behaviour of grounded ice-shelf profiles immediately upstream of the grounding line can be gained from solutions to ice-shelf profiles. The first solutions for steady ice-shelf profiles experiencing spatially uniform surface mass transfer were presented by Oerlemans & van der Veen (1984, henceforth OVdV), for the case of the unconfined ice shelf in plane flow. These solutions can be extended to include the case of where a ‘back pressure’ (i.e. a force different from that provided the pressure of seawater) exists downstream of the shelf. In this case, analytical solutions only appear to exist where the index in the flow-law is integral. We present integrations for *n*=1 and 3. In this paper, we present the OVdV result again but in a slightly different form, and then present the extended solutions for where back pressure exists.

### (a) Steady solution for an unconfined shelf

We consider the plane flow of an ice shelf whose mechanics are given by the usual reduced set of equations (Weertman 1957; Sanderson 1979; Morland 1987; Morland & Zainuddin 1987). The momentum balance is(6.1)where is the longitudinal deviator stress, , where , are the densities of ice and water, respectively, *g* is the acceleration due to gravity, is the thickness of the ice shelf which is a function of the horizontal coordinate *x*. For an unconfined shelf this integrates to(6.2)and if we use the form of Glen flow law appropriate to the case where one stress/deformation rate dominates, *viz.*(6.3)assuming that HP shear is negligible and remembering that by definition , we find(6.4)

The next step is to consider mass conservation. Continuity for the case of a steady-state requires(6.5)and *a* is the accumulation taken to be independent of position. This is not a particularly good approximation for ice shelves, where melting underneath the ice mass can vary rapidly with space; but the purpose is to illustrate concepts. We take the origin of the ice shelf to be at *x*=0 and suppose that a flux is entering the ice shelf, which at this point has thickness .

The reduced mechanical model for ice shelves also gives us the result that horizontal velocity varies negligibly over the depth of the ice shelf. Using and expanding (6.5), we haveNoting that and using the formula (6.3) for the longitudinal deformation rate gives us a first-order differential equation for ,(6.6)

We will now scale this, essentially for reasons of clarification. It will emerge that there are two length scales, one given by and the other given by . We could use either of these to scale the problem, but we shall choose the former because it seems slightly more natural; why this is will be seen in due course.

Thus, representing the normalized horizontal coordinate by *ξ* and the normalized profile of the ice shelf by , defined according to the following scalings(6.7)and we note that the normalized thickness of the ice shelf at the grounding line is 1.

The differential equation for the profile (6.6) becomes(6.8)where(6.9)is a parameter defining the ratio of the surface mass-exchange rate to the creep thinning rate at the grounding line. The solution to the ODE (6.8) is monotone, so for comparison with the confined case, we treat *η* as the independent variable, and write the solution obtained by OVdV as(6.10)Rearranging,and therefore whenthe solution blows up. Clearly, this happens a finite distance upstream of the grounding line . There is an issue of whether this solution is somehow untypical, and that perhaps more generally blow up does not occur, so we need to consider the more realistic case where, for example, of an ice stream resting on a flat bed, the solution has to be generalized to include a back pressure.

In all the cases, considered here, provided and thinning can compensate accumulation at the grounding line, the steady solution is well behaved downstream of the grounding line and the far-field solution is a shelf of constant thickness which can easily be deduced by setting to zero in the governing ordinary differential equation.

### (b) Solutions for floating shelves and grounded streams experiencing downstream resistance

Because ice shelves experience no tangential traction upon their upper and lower surfaces, any forces which oppose their spreading must come from either their downstream surfaces or their lateral surfaces. In the previous case, the confining pressure was that provided by seawater; this is now extended to include a case where the resistance is different to that provided by the pressure of seawater acting as a horizontal force. Such a case could arise when there is a downstream obstruction to ice-shelf flow. In addition, the case where a perfectly slippery ice stream is resting on a flat bed but opposed by seawater pressure at the grounding line creates an expression for the membrane stress which appears to have a back stress.

In plane flow, with a perfectly slippery bed and no downstream body force, the ice-shelf reduced model (4.12) givesNotice that we now multiply the driving stress (the RHS) by the factor *r*, which is the ratio of surface elevation to ice thickness, i.e. *s*=*rH*. For an ice shelf, , while for an ice stream lying on a flat bed, *r*=1. We treat *r* as a constant; it is a simple device for introducing sloping beds. For a marine ice stream, we must have , and one may have *r*>1 for a bed dipping in the direction of flow. We apply one force boundary condition at the grounding line *x*=0. This force can either come from seawater, or, if the abutting shelf is confined, can also include a force from the shelf. The differential equation integrates to(6.11)where(6.12)The important term is the second one on the RHS, which has units of force. If we take the force to be equal to that provided by the pressure of seawater, then .

Using the appropriate (reduced) form of the flow law (6.3) in (6.11) gives(6.13)where *α* is nowwhere is a back-stress parameter defined by(6.14)Substitution for yields an ordinary differential equation for ice-shelf thickness for the confined case, corresponding to (6.6)(6.15)which, upon using the scalings (6.7) and writing gives an ODE for in dimensionless form,(6.16)with boundary condition .

This can be rewritten as(6.17)A formula for arbitrary *n* does not appear to exist. Instead, for each *n*, the term to the power *n* in the denominator of the LHS must be expanded. The resulting polynomial is factorized and the whole of the LHS may be expressed as partial fractions; this may now be integrated. This implies that *n* be integral to get a finite expansion, and that the structure of the solution is different for different *n*. We present solutions for *n*=1 (because this is easy) and for *n*=3.

For *n*=1, equation (6.17) may be written as(6.18)which has a solution(6.19)which is of very similar form to the solution obtained for the unconfined case, except that *β* is replaced by *β*−*δ*.

For *n*=3, expansion of (6.17) yields(6.20)which is best rewritten as(6.21)where the represent the roots of the polynomial forming the denominator of the LHS of (6.20). Because this sixth-order polynomial has only even terms, there is an analytical solution, but this is somewhat unwieldy.

Let us adopt the notation . Differential equation (6.21) may be rewritten in terms of partial fractions(6.22)which has the solutionwhere the boundary condition has been used. For the case *δ*=0 (unconfined and floating), we have , and this solution is equivalent to the solution for the unconfined shelf (6.10). I have not been able to find conditions describing the blow up of this solution in the region , but straightforward calculations indicate that it does. The reason seems to be that generally , which means that a wide range of values of the RHS are explored for small distances from the grounding line.

To demonstrate this blow up further, for sufficiently large (actually, for ) from (6.17) one may writewhich has the solutionwhich gives finite for . In practice, the ice-stream profile blows up a few kilometres upstream of the grounding line for tractionless streams.

Intuitively, one can see that this will hold even when there is some basal friction. Adding basal or lateral friction increases the back pressure in a way that is proportional to distance from the grounding line (Hindmarsh 1993). However, if the blow up happens within a short distance of the grounding line, then this extra back stress cannot be large, so the conditions which permitted the blow-up are maintained. However, we can speculate that if the back stress is large, blow-up will not occur. To proceed with this analysis, we need to understand the conditions under which the zone near the grounding line is shelf-like in character. Some plots of shelf profiles are shown in figure 4. If the origin is regarded as the grounding line, upstream of the grounding line these plots could be regarded as stream profiles with reverse-sloped beds. Qualitatively, the profiles are similar to cases where there is no back pressure. Further investigations show that as one might expect, streams with flat beds have even steeper profiles.

## 7. The membrane stress boundary layer at the grounding line

In general, the stress fields in the shelf immediately downstream of the grounding line are expected to be dissimilar from those upstream of the grounding line, and in between there is a transition zone or boundary layer. By scaling using the shallowness of ice §3, it can be shown that for a nonlinear rheology, or for a sliding bed, that the membrane terms in the momentum balance equation are larger than the ‘bridging’ terms (i.e. the horizontal gradients of , ). This implies that the boundary layer associated with the membrane stresses will be greater in extent than that associated with the bridging terms, which is known to be of the order of the ice thickness in horizontal extent as this term is of (see §3).

### (a) Velocity fields in the membrane stress boundary layer

We now consider the likely extent of the of this boundary layer. One motivation is the result from the last section, which shows that shelf profiles have a tendency to blow-up—in other words a physically reasonable steady profile cannot exist. If a boundary layer with membrane stresses exists upstream of the grounding line, the question arises as to whether it might show the same blow-up.

Horizontal force balance for a membrane in plane flow, with a traction-free upper surface resting on a flat bed, can be written(7.1)where the basal traction is related to the velocity by(7.2)and(7.3)and depending on whether we are considering plane flow or axisymmetric radial flow. At the grounding line, the force balance is given by(7.4)which serves to define *κ*, the confinement ratio (this construction for the stress in terms of the confinement ratio is quite general and does not apply only at the grounding line).

We follow Fowler (1992) in using plug flow asymptotics, meaning that vertical variations in the horizontal velocity *u* are ignored. This is the vertically integrated membrane stress approximation, and is also similar to the one used by Dupont & Alley (2005) and in some of the time-dependent calculations used by Payne *et al*. (2004) and in the static calculations by Schmeltz *et al*. (2002).

The relationship (7.2) can either be regarded as a sliding law or representing shear in the basal layer of ice, in which case it would be appropriate to set . The rate factor *A* in (7.3) should be appropriate for temperatures in the upper part of the ice.

First of all, we consider a flat slab governed by equations (7.1–7.3), with a far-field stress given byand the normal grounding line horizontal force boundary condition (2.6) applied at . For simplicity, we treat the cases where *u* and are positive. The force balance equation can be writtenwhere and writing(7.5)we obtain the equation,(7.6)The quantity is defined(7.7)which parametrized the effect of lateral resistance in ice streams through rewriting equation (7.2),This can be derived by noting thatwhencewhere *W* is the semi-width of the ice stream. Then, choosing allows us to write (7.7). At the grounding linewhich can be written

It can be seen that a length scale *L* is given byFor cold ice, an upper surface value for *B* is 10^{6} Pa yr^{−1/3} and for an ice stream . With a velocity of 100 m yr^{−1} the boundary layer extent is between 20 and 30 km. For a temperate tidewater glacier (*B* is 10^{5} Pa yr^{−1/3}, *ϵ*=0.1, same velocity), the boundary layer is 200 m in extent, while for a cold glacier the membrane stress boundary layer is 1 km in extent.

For a linear rheology and a linear sliding relationship alternative expressions arewhere *η* is the Newtonian viscosity. Substitution of typical values indicates that the boundary layer size depends strongly on the ratio between top and bottom viscosities.

A solution for the force balance equation (7.6) with , is given bywhile for nonlinear rheologies use of a numerical procedure seems unavoidable.

We now introduce the confinement ratio, which is the proportion to which the longitudinal membrane stress is reduced from the value it would have if the resistive force were entirely from water at flotation depth. A confinement ratio of unity implies is zero, and a grounding line can in principle be overconfined.

Figure 5 shows some solutions for velocity fields for different rheologies, different flow geometries and different confinement ratios. The nonlinear rheologies decay over somewhat longer dimensionless distances than the linear rheology. Flow geometry only has a minor effect on velocity fields in the boundary layer, and will not be considered further.

### (b) Steady solutions in the boundary layer

Setting *ϵ*=0 and writing to rewrite (7.1) and using the steady-state continuity equation allows us to write(7.8a)(7.8b)where the flux *q*=*ax* and . For a traction-free bed, we obtainwhere is the back force exerted by an abutting ice shelf. The extra term compared with a normal ice shelf arises from the fact this is a grounded ice stream buttressed by water at the grounding line. If the ice ended in a calving cliff, this term would be zero. For convenience, we write

### (c) Conditional existence of steady profile solutions in the membrane stress boundary layer

As explained above, ice-stream profiles with very low resisting forces have a tendency to blow up a very short distance upstream of the grounding line. The problem stems from the fact that the strain rate increases with the force *F*, which itself is proportional to the thickness for case of zero basal traction. From the continuity equationand differentiation yieldswhere we have writtenWhen , it is straightforward to see that if increases monotonically with ice thickness or , the slope must simply become more and more negative as one moves upstream of the grounding line (), and will blow up. More precisely,(7.9)

We now seek quantitative expressions for . Nowwhich implies(7.10)

We now consider in more detail conditions for . The following argument assumes that is positive at the grounding line and is therefore negative—otherwise the grounded zone would be to the right of the grounding line. If we can then find a condition for , we can guarantee that decreases upstream from the grounding line since elevation increases by construction. From (7.8*a*), clearly can only happen if since by construction. Then from (7.8*a*), and from (7.10) and (7.9)(7.11a)(7.11b)(7.11c)whence, after using (7.2), i.e. , we findWriting , one can show that equivalent conditions are(7.12a)(7.12b)

At the grounding line using (7.4), one can rewrite (7.12*a*) as(7.13)where for completeness, we state that(7.14a)(7.14b)

If inequality (7.13) is respected, the slope becomes increasingly negative as one moves upstream from the grounding line. Note that for a perfectly slippery bed , the inequality will always be respected, and a steady solution can never exist. It must be emphasized that the inequality is a *local* condition, and in principle, as one moves upstream from the grounding line, the inequality could change from not being satisfied to being satisfied, and the profile therefore not blow up.

To investigate this further, let us assume for the moment that is uniform in space. Let us suppose that at the grounding line the inequality is respected, seeming therefore that the solution will blow up with the slope increasing (i.e. becoming more negative) as one moves upstream. If the RHS of (7.13) were to increase, then might no longer exceed it. As one moves upstream from the grounding line, the numerator of the RHS decreases, as flux decreases and thickness increases. Turning now to the denominator, the term increases. Since the slope is negative, the term is positive. Notice that implies that , so the term must therefore increase as well. Finally, we consider the term , where we recall the is the ratio of surface slope to thickness. The ratios *κ* and *r* are both functions of position, but let us consider the case where *r* is a constant. Since at the grounding line and we have only considered cases where between the point of interest and that the strain rate must be positive throughout this section, which implies that . The behaviour of the *κ* is complicated, and it is better to note that in general , and that since *κ* is effectively bounded by 0 and 1, the variation of this term is unlikely to be large, even though, intuitively, one expects confinement to increase with distance from the grounding line and the term to diminish. The two terms already in the RHS of the inequality already discussed are likely to diminish as one moves from the grounding line, and if is constant, and the existence condition will not be traversed. If is not constant, and decreases from the grounding line, one may have situations where the flow appears supercritical at the grounding line but is in fact subcritical.

A particular feature of the inequalities (7.12) is the enormous sensitivity they have to the various parameters—for unconfined flow the RHS is proportional to . When comparing values of with those of the basal rate factor of ice (typically 10^{−16} yr^{−1} Pa^{−3} it needs to be multiplied by the ice thickness, so typical values for when the ice is not sliding might be around 10^{−13} m yr^{−1} Pa^{−3}.

Failing to satisfy inequalities (7.13) at the grounding line does not guarantee the existence of a bounded solution for *H* and *F*. If such a solution exists, one supposes the far-field solution is expected to converge to the SIA solution as in the stress fields in the flat slab example. In this case, we expect a steady solution for the whole ice sheet to exist as for the SIA (Hindmarsh 1996), and do not anticipate any change in the neutral equilibrium property. This could change if the narrower boundary layer of increased bridging stress terms ( terms) added some extra properties.

However, when a bounded steady profile solution does not exist, it appears likely that the neutral equilibrium property for marine ice sheets no longer exists. In principle, the bridging stress boundary layer could introduce other properties, or less likely, re-establish the existence of the steady solution. This is a matter for further research, but further discussion will be based on the assumption that failing to satisfy (7.13) is a necessary condition for the existence of neutral equilibrium. In contrast, if condition (7.13) is respected, the ice sheet cannot achieve a steady state and it will change. It may, by changing its geometry, move into a parameter regime where neutral equilibrium is achievable, and stop changing, it might simply disappear, or it might oscillate.

One glaciological mantra questioned by this analysis is the belief that reverse slopes destabilize marine ice sheets. Since then, if , . Since, *r* appears in the denominator of (7.13), *r*<1 allows greater values of before no solution can exist—in other words, the bed can become warmer and/or more slippery. However, the situation is more complicated. If we find ourselves in a parameter range where no solution exists, and the ice retreats into deeper water, this might increase thinning rates at the grounding line and move the ice sheet even further away from stability. Even this statement has caveats, because such changes will surely change shelf behaviour and the back stress applied across the grounding line.

### (d) Computation of steady solutions in the membrane stress boundary layer

Numerical integration of (7.8) from the grounding line is not a numerically stable procedure, so the equations for *F* and *H* need to be solved as a two-point nonlinear boundary value problem. However, if the solution blows up, the nonlinear solution can never be found numerically. Fortunately, with the condition (7.12) we can at least avoid seeking steady solutions where they do not exist. Steady solutions are found by integrating the time-dependent equations forward in time.

A finite difference algorithm which solves the partial differential form of (7.1) and the time-dependent plane flow continuity equation(7.15a)(7.15b)(7.15c)in one-dimension has been coded. The details are very similar to the algorithm described in detail by Hindmarsh (1996). The thicknesses *H* are solved on a grid that has points aligned with the grounding line and divide, and the other variable solved for is the velocity, on a staggered grid. Symmetry conditions are applied at the divide, and a virtual point beyond the grounding line ensures that the grounding line force boundary condition is respected. No boundary condition is applied to the evolution equation, which is the correct way to treat a hyperbolic equation when there is no boundary domain where the velocity is directed inwards.

Margin motion is computed in the same way as in Le Meur & Hindmarsh (2001). The algorithm can also be specified that it seeks a steady solution for a specified margin position. The margin position evolution equation ensures that the height of the margin is always at flotation provided the initial conditions also respects this condition. This is not really a boundary condition on the thickness, merely a means of marking where the thickness reaches the flotation point. We could model a notional shelf, but since in plane flow this transmits whatever downstream back force we apply to it upstream to the grounding line *whatever its geometry* this is simply an unnecessary decoration for the present purposes. Of course, in reality, shelf dynamics and therefore backpressure are a function of grounding line position, amongst other things.

Writing , , where *S* is the span, equations (7.15) may be written as(7.16)(7.17)(7.18)and linearizing , , we obtain at zeroth-order(7.19)(7.20)(7.21)At first-order, we have the continuity equation(7.22a)the momentum balance equation(7.22b)the first-order downstream boundary condition is(7.22c)and the margin motion equation(7.22d)

This system is discretized and numerical eigenvalues for this system are computed in a similar way to that described for the SIA by Hindmarsh (1996). Programming was checked by comparison with a Jacobian operatorcomputed from numerical differentiation of the time-dependent equations (7.15).

We concern ourselves with some issues of convergence of the solution as the grid size is decreased. A steady state is spun up, with a particular back pressure known to permit a steady profile, i.e. the flow is subcritical. The back pressure is then reduced, so that the ice sheet is closer to the transition in parameter space or has traversed it to a supercritical state. Figure 6 shows the calculation. The results show that the numerical model permits large changes of the ice sheet in some cases where the confinement ratio is subcritical, and further experiments show that this behaviour gets *worse* as the grid size is reduced. Other numerical experiments confirm the existence of neutral equilibrium where steady solutions exist.

This is not entirely unexpected; after all, the condition is framed so that a sufficiently large basal rate factor guarantees non-existence of a steady solution, and we cannot be sure that subcritical basal rate factors guarantee the existence of a solution. However, a second issue is the stability of the solution. An investigation of the linear stability using equation set (7.22) shows that an unstable mode appears near to where the ice-sheet configuration is traversing from subcriticality to supercriticality. When the steady profile can exist, this mode is the neutral mode (see Hindmarsh 1996), and is a *long* wavelength mode. This means that any growth of this mode does not add high frequency information to the ice-sheet profile. Near or at the point in parameter space where the flow becomes supercritical, the numerical evidence suggests that this mode bifurcates to an unstable mode, but remains a long-wavelength mode. Growth rates are of the order of centuries to millenia.

This instability has so far only been demonstrated numerically as it is computed using a numerical discretization of (7.22). Its existence is the reason for the poor convergence and predictability of supercritical marine ice sheet configurations. It is no particular surprise that accurate computation of grounding line motion exhibits difficulties (Hindmarsh 1996; Vieli & Payne 2005).

### (e) Implications for the past and future retreat of the Antarctic ice sheets

A first question is whether the change from subcritical to supercritical occurs at physically attainable values. The answer is almost certainly yes, owing to the extreme sensitivity to the parameters. An example is the sensitivity to the back force, which is approximately with an exponent of .

One test is whether an apparently steady body like Berkner Island is sub- or super-critical. The north end of Berkner Island is essentially unconfined, so one can set the confinement ratio to 0. Taking an upper surface rate factor appropriate to ice at −30 °C, an ice-thickness of 700 m, an accumulation rate of 0.13 m yr^{−1}, a surface temperature of −30 °C and a flow-path length of 100 km (R. Mulvaney 2005, personal communication) gives a critical basal rate factor of around 1.3×10^{−17} Pa^{−3} yr^{−1}. If the basal rate factor is greater than this, no steady solution exists. This rate factor corresponds to a basal temperature of around −10 °C, which is a plausible value for Berkner Island. In other words, Berkner Island is apparently near the transitional value. Altimetric evidence from Berkner Island suggests that it is not changing fast at the moment (Arthern & Hindmarsh 2003). An interesting conjecture is that Berkner Island has retreated until the existence criterion became satisfied again, perhaps through the bedrock rising.

The larger the RHS of (7.12) the more easily the condition for existence of a steady profile is respected. The maintenance of subcriticality is maintained if the grounding line flux increases, either because the accumulation rate increases or the ice sheet gets bigger. The latter is important, as it suggests bigger ice sheets are stable ice sheets. Increased accumulation after the glacial maximum should have helped maintain the ice sheets. The condition is more sensitive to increased strain rates at the grounding line, which can either happen because of a change in the back pressure—implying a change in the ice shelf—or because the rate factor in the upper layer has increased due to warming. This could be responsible for retreat in areas where there were not large shelves, for example the Antarctic Peninsula or East Antarctica.

## 8. Concluding comments

This paper presents some new results as well as being a review of membrane stress effects in ice sheets. The three main features of interest are: (i) the very long distances over which effects can be transmitted ice streams for nonlinear rheologies; (ii) the existence of a boundary layer in the membrane stresses extending some tens of kilometres in from the grounding line, even for ice adhering to the bed; and (iii) the conditional existence of steady profiles in the boundary layer and the conditional existence of a long wavelength instability.

The approach to examining the very long distance transmission of stresses in ice streams is essentially another way of formulating the theories of Thomas (2004), Payne *et al*. (2004) and Dupont & Alley (2005). As it considers infinite sections, the overall length-scale of transmission is more clearly emphasized. The fact that it is formulated as an eigenvalue problem means that the stability can be examined. As expected (cf. Gudmundsson 2003; Hindmarsh 2004), perturbations are stable. Even though changes in the stress can create rapid changes, the eventual response is to re-equilibrate.

Ice streams are not quite the same things as the membrane stress boundary layer, a concept developed independently by C. Schoof (personal communication). This is a vertical boundary layer over which membrane-like stresses at the grounding line decay onto those consistent with the SIA. Typically, this boundary layer is a few tens of kilometres in extent, and can be shorter than an ice stream.

The existence of this boundary layer has the following important consequence. If the ice can, in the region near the grounding line, move over too easily to be replaced by interior flow, a bounded steady profile cannot exist. This can happen if the bed is too slippery, the downstream resistance is too small or the influx of ice too small. The conditional existence of steady solution suggests many ways of causing ice-sheet change, as it can be expressed as a requirement for the basal rate factor or sliding viscosity to be smaller than a quantity which depends upon the velocity, upper surface viscosity and back pressure at the grounding line. If no steady solution exists, and the ice sheet were to retreat, enhanced flow would be seen in the membrane stress boundary layer. The retreat could be likened to a process of spallation.

This condition does not have to be explicitly inserted into ice-sheet simulators—any model which includes membrane stresses should demonstrate its consequences. We have used a very simple vertically integrated model, and more investigations using better representations of membrane stresses will doubtless change the quantitative features of the results.

A further complication is the apparent presence of an unstable mode when steady solutions do not exist. This mode is of long wavelength, comparable with the span of the ice sheet, and has slow growth rate. The instability has been demonstrated by solving a numerical eigenproblem. One can be certain that the instability, which is displayed in time-dependent finite difference calculations, is not introduced by the time-stepping procedure, but it remains to be determined whether it has been introduced by the spatial discretization scheme. It should not be surprising that an instability exists when no steady solution is possible.

Finally, the fact that a ‘tipping point’ for marine ice sheets has been identified should be stressed. Perturb the system too far, and changes start which may be irrevocable. This tipping point can be passed by a number of different perturbations, principally the removal of shelf buttressing, but also warming of the upper surface and decrease of accumulation. Since retreat leads to loss of catchment area, it can lead to the ice sheet moving away from the transition point in parameter space as well. The numerical evidence presented, both through time-dependent finite difference calculations and the solution of a numerical eigenproblem, suggest that the loss of existence of a steady solution does not lead to an oscillatory solution, but to large-scale change. Schoof (2006, personal communication) argues that the inequalities (7.12, 7.13) can infact be used to derive similar equalities, which form boundary conditions relating the horizontal deviator stress at the grounding line to the mean velocity. Steady solutions only exist when the equality holds. This would have entirely different consequences for the ice-sheet dynamics.

## Footnotes

One contribution of 14 to a Discussion Meeting Issue ‘Evolution of the Antarctic Ice Sheet: new understanding and challenges’.

- © 2006 The Royal Society