Dynamical processes involved in the retreat of marine ice sheets

Abstract

Marine ice sheets with mechanics described by the shallow-ice approximation by definition do not couple mechanically with the shelf. Such ice sheets are known to have neutral equilibria. We consider the implications of this for their dynamics and in particular for mechanisms which promote marine ice-sheet retreat. The removal of ice-shelf buttressing leading to enhanced flow in grounded ice is discounted as a significant influence on mechanical grounds. Sea-level rise leading to reduced effective pressures under ice streams is shown to be a feasible mechanism for producing postglacial West Antarctic ice-sheet retreat but is inconsistent with borehole evidence. Warming thins the ice sheet by reducing the average viscosity but does not lead to grounding-line retreat. Internal oscillations either specified or generated via a MacAyeal–Payne thermal mechanism promote migration. This is a noise-induced drift phenomenon stemming from the neutral equilibrium property of marine ice sheets. This migration occurs at quite slow rates, but these are sufficiently large to have possibly played a role in the dynamics of the West Antarctic ice sheet after the glacial maximum. Numerical experiments suggest that it is generally true that while significant changes in thickness can be caused by spatially uniform changes, spatial variability coupled with dynamical variability is needed to cause margin movement.

1. Introduction

The main proposition of this paper is that our understanding of how marine ice sheets retreat is rather less complete than has previously been recognized. Recent results (e.g. Licht and others, 1996; Conway and others, 1999; Anderson and Shipp, 2001) strongly suggest that the postglacial retreat of the Siple Coast, West Antarctica, occurred very late and was therefore not directly related to sea-level rise. Evidence from moraines (Ackert and others, 1999) shows that the West Antarctic ice sheet (WAIS) thickened after the glacial maximum, presumably due to increased accumulation, but subsequently has thinned slightly. This thinning is consistent with the later retreat of the WAIS but remains to be shown to be coeval.

It has been thought that the combined action of sea-level rise coupled with decay of ice shelves and the consequent loss of buttressing was sufficient to explain any retreat or thinning of the WAIS. However, it has become increasingly clear that buttressing is not as important as has been supposed (Hindmarsh, 1993; Whillans and Van der Veen, 1993; Echelmeyer and others, 1994; Paterson, 1994; Bentley, 1997), and indeed mechanically self-consistent marine ice-sheet models which do not need to compute the evolution of ice shelves have been constructed (Hindmarsh, 1996; Le Meur and Hindmarsh, 2001). The latter paper shows that the observed retreat of the Siple Coast (500 km) is “large”; if the thickness gradient of the WAIS were 10⁻³ the response to a sea-level rise of 100 m should have been 100 km and should have been synchronous with postglacial sea-level rise.

In this paper, marine ice-sheet models will refer to a class of models uncoupled across the grounding line. There is some theoretical background to this class of models in papers by Hindmarsh (1993, 1996) and Le Meur and Hindmarsh (2001). These models of marine ice sheets are atypical in that they do not include mechanical coupling across the grounding line. That this coupling is important has been an axiom of glaciology. However, if the flow of the grounded ice is well described by the shallow-ice approximation (Hutter, 1983), longitudinal stresses are not significant and coupling with the shelf is not important. The models also suggest that the dynamics of marine ice sheets are likely to be rather peculiar, and are organized around the neutral equilibrium property. We can contrast land-based ice sheets, where one or a few discrete margin positions corresponding to equilibrium profiles exist, to marine ice sheets, where the equilibrium position is continuous, meaning that an infinity of equilibrium configurations exist. In a land-based ice sheet, ablation is a function of ice-sheet geometry, and only a few margin positions permit total accumulation and ablation to balance and the profile to remain unchanged. In contrast, in a marine ice sheet, the grounded sheet is drained by a shelf. Any flux (more properly, the mean velocity) is a valid (Dirichlet) boundary condition for the elliptical shelf momentum-balance equations, which can, in consequence, remove any amount of ice flowing across the grounding line. Thus, any margin position can correspond to zero total mass balance for the grounded sheet. It may transpire that the shelf cannot exist in steady state without further grounding or shoaling, but this does not obviate the basic principle of neutral equilibrium.

Hindmarsh (1996) has shown that the neutral equilibrium of marine ice sheets is an attracting state. Consider a temporary perturbation in the forcing of an ice sheet in equilibrium. The ice sheet will move away from equilibrium, then return to equilibrium once the forcing perturbation terminates. For a land-based ice sheet, where there are distinct equilibrium points, this means that the ice sheet will return to the same span and divide thickness (ignoring, for the moment, difficulties associated with multi-stability; see Oerlemans and Van der Veen, 1984). For a marine ice sheet, a consequence of neutral equilibrium is that this need not necessarily be so. Since any span is possible, a departure from equilibrium can imply, once the same climate has been re-established, a decay onto a different equilibrium position, even though the final climate is identical to the initial climate.

A further consideration is that systems in neutral equilibrium forced by noise exhibit “noise-induced drift”. A good example is an electric razor placed on a table. When the motor is switched off, the device is in neutral equilibrium; it could be in any position on the table. When the motor (noise source) is switched on, the razor migrates across the table; this is noise-induced drift. It could be that fluctuations in ice-sheet average viscosity (corresponding roughly to stream switching) cause the marine ice sheet to exhibit noise-induced drift (i.e. grow or shrink). Noise-induced drift is an inherently non-linear phenomenon, and arises from an asymmetry in the response of the system. In this paper we consider whether noise-induced drift may be a significant mechanism in controlling marine ice-sheet retreat.

These considerations suggest that the dynamics of marine ice sheets require further investigation by modelling. However, capturing neutral equilibrium with numerical models is a delicate issue (Hindmarsh, 1993, 1996), and a model which behaves reliably has only recently been developed (Le Meur and Hindmarsh, 2001). Hindmarsh (1996) showed dynamical consistency between a flowline finite-difference (FD) model and the governing equations by demonstrating the existence of a neutrally stable mode and showing that a corresponding mode also existed in the FD algorithm. Demonstrating that a numerical model possesses a neutrally stable mode requires a moving margin, and the algorithmic difficulties inherent in including this approach in a three-dimensional model make this next step non-trivial.

Nevertheless, certain facts about the dynamics of marine ice sheets seem reasonably certain. The models of Hindmarsh (1993, 1996) and Le Meur and Hindmarsh (2001) allow us to state that a marine ice sheet will not change its geometry in areas that remain grounded while sea level rises, provided (i) its mechanics are well described by the shallow-ice approximation (Hutter, 1983); (ii) there are no evolving interior fields such as temperature or basal water pressure; (iii) there is no sheet–shelf coupling; and (iv) it experiences an unchanging accumulation rate. In brief, we do not expect an automatic thinning on sea-level rise. Furthermore, we note that the dominant processes of heat transport in an ice sheet are vertical conduction/advection and horizontal, downstream advection. Therefore, any part of ice that is floated off will not affect the flow in the upstream area, which, under the shallow-ice approximation, is at a higher elevation and in general remains grounded. In consequence, we do not expect direct coupling between sea level and ice-sheet temperature, except possibly from the fact that 100 m sea-level rise will cause temperatures at any point to rise by the product of the lapse rate and the sea-level rise, which under usual circumstances is a bit less than 1 K.

Ensuring that the interior ice-sheet geometry remains unchanged under these conditions is actually quite difficult to capture numerically; numerical schemes propagate numerical error from grounding line to divide, and very high-accuracy FD schemes are needed to avoid this problem (Le Meur and Hindmarsh, 2001). The contention of this paper, which we demonstrate by modelling, is that this decoupling between ice thickness and grounding-line position is in a more generalized sense the “expected” behaviour of marine ice sheets, even when one includes interior fields which affect the rheology. Correlations of grounding-line retreat and thinning are not necessarily expected. Furthermore, it turns out that causing the grounding line of a marine ice sheet to retreat presents rather more of a theoretical challenge than has previously been suggested. This paper shows, for example, that large warming, when propagated through to the base of the ice sheet, does not cause significant grounding-line retreat.

Current Antarctic models do show thinning and retreat during the most recent deglaciation (Huybrechts, 1992), but the two mechanisms which have been incorporated in these models and which are believed to be crucial in causing the modelled WAIS to undergo substantial retreat (ice-shelf buttressing and ice streams with static hydrologies) are now widely believed to be wrong, and we still have the problem of how a marine ice sheet retreats. The most likely candidate for mechanisms responsible for promoting ice-sheet retreat is ice-stream activity. Basal hydraulics are of importance because it is widely believed that the effective pressure (difference between ice pressure and water pressure) affects the sliding viscosity. Two earlier models (Budd and Smith, 1981; Huybrechts, 1992) used a hydrostatic model of basal water pressure, which implies that any well bored into an ice stream will only be filled with water up to sea level (this is the so-called “height above sea level” model). Numerous observations of water level in drillholes (Engelhardt and Kamb, 1997) contradict this assumption, nor is it any theoretical surprise that water pressures are found not to be statically controlled (Lingle and Clark 1985; Boulton and Hindmarsh, 1987; Hindmarsh, 1998). The significant aspect of such static hydrologies (i.e. where the interfacial effective pressure is determined by hydrostatic principles) is that they link ice-sheet dynamics directly to sea level; a rise in sea level will by construction lower the effective pressure, increase ice-stream velocities and cause the ice sheet to thin. However, since real ice-stream hydraulics are essentially unaffected by sea level, raising sea level should have no direct effect on ice-stream dynamics. In short, static hydrologies would seem to automatically cause ice-sheet thinning when sea level rises. We examine this in this paper by modelling and find that the amount and timing of retreat is very sensitive to the sliding-model specification. Maximum retreat can occur several thousands of years later than the initiation of sea-level rise using this model.

In brief, this paper is concerned with (a) extending a theme of the companion paper (Le Meur and Hindmarsh, 2001), seeking to demonstrate that mechanisms previously used to explain retreat are inadequate; and (b) understanding how the neutral equilibrium property of marine ice sheets organizes their dynamics, and whether this can lead to a different hypothesis about marine ice-sheet retreat. In this paper, we show that (i) use of a static hydraulics does lead to substantial retreat upon sea-level rise, which implies that marine ice-sheet models which use this hydraulics may be computing erroneous retreats; (ii) a warming on its own does not lead (at least automatically) to a retreat of the grounding line; in fact, while it leads to a considerable thinning, grounding-line position is scarcely affected at all; (iii) fluctuations in the ice-sheet viscosity can lead to a noise-induced drift phenomenon; and (iv) internal thermo-viscous oscillations of the type proposed by MacAyeal (1992) and Payne (1995) can contribute to ice-sheet retreat via a noise-induced drift mechanism.

2. Model Formulation

All the modelling in this paper is carried out with a two-dimensional plane-flow model, with dimensions chosen to represent a typical marine ice sheet. Clearly, better realism could be attained with a map-plane model, but no such model exists where dynamical consistency has been demonstrated. We expect some of the dynamical features of the plane-flow model to extend to map-plane models. Figure 1 is a sketch illustrating the configuration and the notation.

Fig. 1. The ice-sheet configuration modelled in this paper. The distance s(x, t) − Z(x, t) represents the freeboard; Z =.(ρ _u /ρ _i )(f − b).

2.1. Isothermal model

The mechanical model is as described by Hindmarsh (1996), which is a variant of the shallow-ice approximation (Hutter, 1983; Morland, 1984; Fowler, 1992). Dimensional quantities are indicated with a tilde. The governing equation we are considering is, in dimensionless form,

(1)

where H(x, t) is the thickness of the ice sheet, s(x, t) is the upper surface and a is the surface mass-balance exchange. A summation is implied over j, representing N different modes of ice-sheet flow by deformation and sliding. The horizontal, vertical and time coordinates are (x, z, t). This equation describes the evolution of ice-sheet thickness where the flow mechanism is by internal deformation according to some non-linearly viscous flow law and by sliding according to some Weertman-type law. We can include both effects by having a summation over j in Equation (1).

The quantity C_j is either directly related to a rate factor A _d used in the viscous relationship

(2)

where E is a second invariant of the deformation rate and σ is a second invariant of the deviator stress (Glen, 1955), or comes from a sliding relation of the form

(3)

where A _s is a sliding rate factor, σ _b is the basal shear stress and p _e is the effective pressure. The exponents n, ℓ and λ are phenomenological variables. We construct the following quantities for use in the general evolution equation:

(4a)

(4b)

(4c)

(4d)

The Glen exponent n = 3, while the other exponents are varied. In this paper we set the effective pressure P _e = ρ _i gH − ρ _w g(f − b), where f(t) is sea level (we ignore distortion of the geoid in this paper), b is the elevation of the ice-sheet base, g is the acceleration due to gravity and ρ _i , ρ _w are the densities of ice and water, respectively. This effective-pressure model contradicts observations, but we seek to show that this sliding law, a feature of many previous models, when used in computations with a rising sea level, produces retreat. In some of the experiments we set λ = 0, i.e. we exclude the basal hydrology from consideration. Hindmarsh (1998) has discussed the role of scale in determining whether static hydrologies are important or not, showing that at shorter length scales, determined by the hydrological properties of any subglacial aquifer, static hydrological effects may be important. Typically, such length scales are of the order of a few kilometres. There exists the possibility that static hydrological effects may be important in the area immediately adjacent to the grounding line. Such effects are not explicitly considered in this paper; we assume either that hydrology has no effect on sliding, or that static effects determine the effective pressure.

The derivation of the evolution equation (1) follows Hindmarsh (1996), who uses the same scaling as Hutter (1983), Morland (1984) and Fowler (1992). In this paper, vertical distances have been scaled by a thickness magnitude , horizontal distances by a magnitude , accumulation rates by [ã], time by and rate factor by

(5)

The magnitudes used in this study were , [ã] = 1 m a⁻¹, , ε =1/750, and . These imply that with n = 3 for internal deformation . Two sliding configurations were used: (i) ℓ = 3, λ = 0, with , and (ii) ℓ = 3, λ = 1, with . We also set g = 1. Note that the results (section 3) are reported in physical units. We used dimensionless ice and water densities equivalent to a density of 917 and 1030 kg m⁻³.

The flotation condition is

(6)

where quantities are taken considered at the grounding line, and again following Hindmarsh (1996) we can obtain an expression for the grounding-line migration rate G _r

(7)

where all the quantities on the righthand side are evaluated at the grounding line.

The algorithms used for steady and time-dependent calculations were as described by Hindmarsh (1996), with the exception that the margin slope was described by a fifth-order formula (see also Le Meur and Hindmarsh, 2001). This more accurate formula does not affect the property possessed by the FD algorithm of respecting the zero eigenvalue. An important point to note is that the equations were transformed onto a moving domain, which permits accurate computation of grounding-line motion. Rather than having a fixed timestep as in Hindmarsh (1996), time-steps were selected automatically, such that

(8)

where Δ_ξ = 1/(M−1) and M is the number of gridpoints. The diffusivity is evaluated at grid centres, as are the thickness, surface slope and effective pressure. The flux discretization is the one-dimensional version of method I or II described by Hindmarsh and Payne (1996). This stability criterion is considerably tighter than needed for ice-sheet models discretized on fixed grids, and is believed to be necessary because of the destabilizing influence of the central difference first-order operator (Hindmarsh, 1996) which appears as a result of grid motion. We have avoided the obvious stabilizing remedy of using upwinded differences, as they introduce an inbuilt asymmetry which could influence our investigations into noise-induced drift.

Finally, we wish to reintroduce the useful concept of “freeboard”. Consider ice with elevation Z(x, t) = (ρ _w/ρ _i) (f(t) − b(x, t)). This is the elevation of ice which is just about to float. The freeboard is given by s(x, t) − Z(x, t), and represents the thickness of ice above flotation. By construction, ice at the grounding line has zero freeboard.

2.2. Thermal coupling

The significance of thermal coupling lies in the fact that the rate factor in the viscous relationship (Equation (2)) and the occurrence of sliding are strongly temperature-dependent. There have been numerical calculations which indicate that this can lead to oscillatory behaviour (Payne, 1995; Greve and MacAyeal, 1996; Pattyn, 1996; Payne and Dongelmans, 1997). The model for thermal coupling is based on the scalings of Morland (1984), who showed that the heat-transport equation can be written

(9)

(10)

(11)

(12)

where

(13)

(14)

and where θ is the temperature, θ ^s is the temperature at the upper surface, θ _m is the melting-point temperature of ice at the ambient pressure, u is the ice velocity, u ^b is the ice velocity at the base, Q _G is the geothermal heat gradient, is the tangential traction at the base, c is the specific heat capacity of ice, k is the thermal conductivity of ice, κ is the thermal diffusivity of ice and S _y is the number of seconds per year. We set , , and . The pressure-melting point deviation is computed by

(15)

(Paterson, 1994). We set . The dependence of the rate factor on temperature was modelled by

(16)

Here R(θ, p) incorporates a dual-Arrhenius relationship

(17)

where (c ₁, c ₂)= (3.7 × 10⁶, 5.4 × 10²⁶), (E ₁, E ₂) = (60, 140) kJ mol⁻¹, and is the universal gas constant. Here, θ _H is the homologous temperature θ/θ _m(p), and this relationship thus introduces an implicit dependence on the pressure p. The quantity is the melting point at zero pressure. This dual-Arrhenius relationship is analogous with and is quantitatively very similar to that proposed by Hooke (1981) and is preferred here simply for its dual-Arrhenius form. In practical terms it produces different activation energies above and below 263 K. Surface temperature θ ^s was prescribed as a linear function of elevation, i.e. , where is the temperature at z = 0 and Γ is the lapse rate. Since the upstream velocity is zero and the grounding line is a hyperbolic outlet, there is no need to specify thermal boundary conditions on the vertical faces at x = 0 and x = S(t).

A serious deficiency of the present model is that it does not deal with the problem of temperate layers within the ice (e.g. Greve, 1997). These layers generally arise from internal dissipative heating, although in principle pressure release can cause ice below the melting point to reach the melting point. However, once the pressure-melting point is reached at the base, sliding can occur, and in this case the dissipation necessarily occurs at the base through frictional heating. If a situation is achieved where the dissipation is largely due to sliding rather than internal dissipation, we can further conclude that the temperature maximum must be at the base, since there is no significant internal heating. This further implies that the ice above the base will be below the pressure-melting point providing the temperature gradient is greater than the pressure-melting gradient. In the absence of a proper model of temperate layers, computations which predicted ice at the pressure-melting point at the base only were retained, while computations which yielded predictions of temperate ice (due ultimately to internal heating) within the ice sheet were rejected. For reasons described above, this happened when there was low sliding and large dissipative fluxes leading to large dissipation within the shear layer.

A heuristic which we used to increase the numerical stability of the scheme, which also has some physical motivation, is to increase the range of temperature over which sliding starts. Thus, the sliding rate factor, rather than being a Heaviside function at the pressure-melting point, is changed to a function defined by

(18)

where is the rate factor at the pressure-melting point. The physical motivation stems from work by Kleman and others (1999) and Wilch and Hughes (2000), who argue that melting occurs in sub-grid patches at a glacier bed as the temperature approaches the melting point. As the proportion of warm-based area increases, the basal resistance decreases as slip can occur over more and more warm-based areas, until at the pressure-melting point, full sliding can occur. The constant γ was set to 1 K⁻¹ in the calculations reported here; this does not have a theoretical basis, but was found to avoid numerical instabilities. At 2 K below the pressure-melting point, sliding is 0.14 of its value at the pressure-melting point.

2.3. Numerical solution procedure

We use the transformation of Hindmarsh and Hutter (1988) to render the heat-transport equation into normalized domains in the vertical and horizontal (this is Jenssen’s (1977) σ coordinate with a horizontally extending mesh). The vertical advection is computed using the formula due to Hindmarsh (1999). While the ice-sheet profile evolution equation is solved using an explicit matching scheme, the heat equation is solved sequentially (i.e. subsequently to the profile equation) with a semi-implicit method (i.e. with velocities, dissipation, etc., computed using “old” values of the temperature). It is not solved as often as the profile equation (i.e. the time-step is longer), as stability problems are not so severe. The temperature equation was solved on a grid staggered with respect to the profile grid.

The ice-sheet profile is second-order accurate when in steady state, and first-order (with accuracy determined by the grid spacing at the margin) when running in time-dependent mode. The heat equation was solved in the vertical with pseudo-spectral methods (Fornberg, 1996; Hindmarsh, 1999) and in the horizontal with first-order FD methods. Use of higher-order schemes in the horizontal frequently led to the appearance of “modes” and deleterious numerical effects on the solution. Long-period oscillations analogous to those reported by Payne (1995) were found under certain parameter combinations. These were robust under changes of the time-step, although the period of the oscillation did change with changes in time-step.

Certain procedures were adopted since steady states of thermally coupled ice sheets are notoriously difficult to obtain reliably (Hutter and others, 1986). This difficulty is partially due to the multi-stable properties of these flows (Clarke and others, 1977; Fowler and Larson, 1980). For example, iteration sequences often fail to converge, and running a model computation to steady state can be enormously computer-intensive since thermal time-scales are much longer than ice-sheet time-scales (Hindmarsh, 1990). To obtain steady states, we adopted a procedure of running the ice-sheet profile evolution equation prognostically, while solving the heat equation diagnostically (i.e. solving for the steady state). This automatically ensures that the ice sheet is in thermal equilibrium, but does not guarantee convergence to steady state of the profile equation. If a steady-state solution is obtained, the procedure yields a physically meaningful steady solution. In cases where this procedure did not lead to a steady profile, we then ran both temperature and profile equations prognostically, while keeping the span fixed, until the system had settled into a pattern of repeating oscillations. This was used as an initial condition for the full time-dependent problem, with the convenience of having a specified span. Since the ice-sheet and temperature equations are initial value problems, this is a mathematically valid procedure. We started the full time-dependent equations during periods when the state was changing slowly, which meant that the “initialization shock” was rather small in the thermally coupled cases.

The procedure for dealing with the switch in the basal boundary conditions was as follows. If the ice at the base went above pressure-melting point, the basal boundary condition was changed in the next time-step so that basal temperature was prescribed (a Dirichlet condition). These cases were checked after the time-step, and if it was found that the ice was cooling, the boundary condition was set to the Neumann condition for the ensuing thermal time-step.

3. Model Experiments

The models are designed to be representative of the Siple Coast/Ross Sea area of the WAIS. Currently, the divide to grounding-line distance is about 500 km, while at the glacial maximum it may have reached 1000 km. These models are plane-flow idealizations of marine ice sheets, with simplified basal geometry (a flat bed) and spatially uniform accumulation rate. In all cases, the ice sheet was set to be symmetric about the divide (which is therefore a zero-flux boundary condition under the shallow-ice approximation), the bed was at −500 m and the sea level was initially at 0 m. In one set of experiments, sea level was raised. Initial steady surface topography for the isothermal model was computed using the integration technique described in Hindmarsh (1996). In other cases the computation of the initial conditions is described in the text. The initial accumulation rate was set at 0.1 m a⁻¹. It was changed in some of the experiments with the isothermal model. The ice-sheet rate factor at pressure-melting point was set to be 10⁻¹⁶ Pa³ a⁻¹ _, which corresponds to , apart from two thermally coupled experiments, where . In the isothermal experiments the rate factor of the ice corresponded to that at melting point. In all cases where sliding occurs (i.e. ) we set the basal tangential traction index for sliding ℓ = 3. The results are reported in physical units.

3.1. Isothermal model

We first report some experiments which do not include thermal coupling. These are designed to look at cases where the issue of thermomechanical coupling is subordinate. Firstly, we investigate the response of the ice sheet to step changes in accumulation rate and viscosity. A second set of computations investigate the influence of hydrology on ice-sheet and ice-stream response to sea-level rise. Specifically, we are interested in whether “static” hydrologies automatically produce ice-sheet retreat once sea level rises, as argued above. The final issue is a stochastic stability issue: what influence does noise-induced drift have on grounding-line motion, and what characteristics of the noise affect the rate of migration of the grounding line? We represent noise by periodic variations in the sliding viscosity in the downstream part of the ice sheet.

3.1.1. Step changes

Here we computed steady-state profiles, applying a step change in either accumulation or rate factor at t = 0₊. The evolving profile was computed until steady state was reobtained. The initial span was 750 km, and the base case corresponds to ice at the pressure-melting point with unit accumulation rate. Forty-one gridpoints were used. Table 1 shows computed changes in the re-equilibrated divide thickness and span at large t. The uniform changes in rate factor do not have a physical interpretation; we are simply demonstrating the consequences of such changes for the ice-sheet thickness and span. It can be seen that these spatially uniform changes in accumulation rate and rate factor have a far greater influence on thickness than on span. This should be compared with the land-based case, where H ∼ S ^1/2.

Table 1. Change in steady thickness H and span S in response to step changes in accumulation rate a and rate factor C Parameter value at time t = 0₊ is indicated by subscript +

3.1.2. Hydrostatic hydrologies

These experiments comprised the computation of an initial steady profile, and subjecting it to a sea-level rise. This should make the ice flow faster, as the effective pressure at a point decreases with the increase in sea level. We used a sliding law with ℓ = 3, λ = 1, and with four cases of dimensionless sliding rate factors with a scale factor [A _s] = 3.91 × 10⁻³ Pa⁻² m a⁻¹. The steady span was set to be 500 km. A stable steady profile was computed, and this was subject to a sea-level rise of 0.01 m a⁻¹ from . The computation was carried out until the span became < 150 km. Forty-one gridpoints were used. Results are shown in Figure 2. Figure 2a shows the initial profile, and also shows Z(x, t ₁₅₀), where t ₁₅₀ represents the time the span reached 150 km for each of the four cases. Where in x-space this line intersects the initial profile shows what the span would have been had there been no coupling between sea level and sliding. Figure 2b shows the span plotted against time. In all cases, there is a period of slow retreat followed by a very rapid retreat. The time to collapse of the ice sheet depends very sensitively on the dimensionless rate factor . This number can therefore be tuned to produce a significantly delayed grounding-line retreat. Note also that very small differences in initial profile can be associated with very different retreat histories.

Fig. 2. Response of an ice sheet in plane flow with an effective-pressure-dependent sliding law and a statically determined effective pressure. Initial profiles are show in (a) for different sliding coefficients (indicated in legend). (b) shows evolution of span in an experiment where at t = 0₊ sea level starts to rise at 0.01 m a⁻¹ The computations were terminated when the span reached 150 km. The ice elevation with zero freeboard at final sea level is also indicated in (a). Note how sensitive the behaviour is to sliding coefficient, and that total retreat is much more than can be attributed to sea-level rise. The static model contradicts borehole observations.

This example shows that having effective-pressure-dependent sliding, with a statically determined effective-pressure-dependent distribution, can easily produce significant retreat in model computations. If a model relies on this mechanism (e.g. Budd and Smith, 1981; Huybrechts, 1992), then a significant model component is in conflict with borehole observations of water pressure (Engelhardt and Kamb, 1997).

3.1.3. Noise-induced drift

We now consider the influence of fluctuations as a means of stimulating migration of the grounding line through a noise-induced drift mechanism. We considered only internal deformation, and we specified the dimensionless rate factor A _d to vary periodically in time over the downstream half of the ice sheet such that

The function F(t) represents a square wave with an initial taper; this prevented unreasonably large initial fluxes during the fast period. We considered two periods representing possible ice-stream fluctuation periods, , 500 years.

The computations, which used 121 gridpoints, were carried out as follows. The rate factors A _f were chosen to produce reasonable balance velocities (around 500 m a⁻¹). The systems were allowed to settle down into a steady oscillation without the span (750 km) being allowed to change. Once a steady oscillation was reached, the span was allowed to evolve at t _{mod T} /T = 0. The integrations were carried out for 20 kyr or until the span decreased to 150 km. Computations with 21 gridpoints were also initiated by using the computed profiles from the 121 gridpoint calculations at 1000 years, in order to assess the influence of truncation error on computed drift rates. (The initial transient response was strongly dependent on grid interval, which is why we adopted this procedure.)

The coefficients chosen were A _f = 2.35 × 10⁵ and 2.65 × 10⁵ for and 500 years, respectively. These values produced mean velocities comparable with velocities measured in modern ice streams (around 500 m a⁻¹). The results are shown in Figure 3. Figure 3a shows the evolving span, demonstrating that there is an initial response, after which the oscillations settle into a pattern. In the two cases we report in detail, this leads to a mean negative drift in the span after an initial response lasting around 1000 years. Figure 3b and c show the maximum mean velocity during a time interval of 2–3 years for the two different forcing periods. It fluctuates significantly, but, apart from the initial transients, remains at values comparable with those observed in modern ice streams. Figure 3d shows the evolving profile.

Fig. 3. Effect of noise due to internal oscillations of an ice sheet on grounding-line position. (a) Evolving span for internal oscillations with period 100 and 500 years. These were computed with 121 gridpoints; the dotted line shows the same computations with 21 gridpoints. The low-resolution computation for the higher-frequency forcing is coincident with the higher-resolution forcing (b, c) The maximum mean velocities in the ice sheet over short time intervals (2.5–3 years). (d) The evolving ice-sheet profile over two 500 year oscillations.

This drift effect could be large enough to be a significant factor in the dynamics of the WAIS. Changing the coefficient A _f caused different sorts of behaviours: sometimes the span increased; sometimes it oscillated about a steady value.

3.2. Thermal coupling

In this subsection we compute coupled evolutions of profile, span and temperature. Interest in coupling temperature and ice-sheet evolution principally relates to the dependence of ice viscosity on temperature and the fact that sliding starts when the base of the ice approaches melting point. We consider whether this coupling can increase the rate of retreat of the ice sheet. For example, one might suggest that the warming at the end of the glacial maximum would have propagated through in a few thousand years (Whillans, 1981), increasing flow rates, resulting in a flatter glacier and possibly more retreat. Another possibility is that thermal oscillations (MacAyeal, 1992; Payne 1995) might provoke a noise-induced drift. We investigate these hypotheses by modelling.

3.2.1. Specification

Twenty-one gridpoints were used in the horizontal for both elevation and temperature, and eleven points were used in the vertical. Since we are using pseudo-spectral methods for the vertical discretization, the dominant spatial truncation error emanates from the horizontal discretization. With these grid spacings, it has a relative error of 0.25% in steady-state calculations and 5% in time-dependent calculations. Computations were done with a 10 year thermal time-step and a time-step in the ice-sheet equation varied to avoid blow-up using Equation (8).

We carried out a number of computations, trying to reach a steady state as described above, varying parameters such as sea-level temperature and the amount of sliding. In some cases, a steady solution could be computed, and this was used as an initial condition, representing a physically valid steady state, in subsequent moving-margin computations. In other cases, no such solution was found when using the procedure described above of running the ice-sheet profile evolution equation prognostically, while solving the heat equation diagnostically. In this case, subsequent calculations ran both equations prognostically with the span fixed. This continued until the ice sheet settled into a steady oscillation, which was used as an initial condition for true moving-margin experiments. As described above, this is a physically valid initial condition. This procedure allowed us to start experiments from a given span. We report on the results of a few of these steady and oscillatory calculations below.

3.2.2. Non-oscillatory computations

We could not find steady solutions with sub-zero internal temperatures for but were able to do so for . This is because the higher rate factor results in a lower slope and less dissipative heating. A steady solution, with a sea-level temperature of −30°C, a lapse rate of 8 K km⁻¹ and no sliding, is shown in Figure 4a. The ice thickness at the divide is around 2.5 km, which means that the Peclet number is 6.9. This leads us to expect the observed basal boundary layer in the temperature around two-fifths of the ice thickness (Morland, 1984), with evidence of advection dominating in the top part of the ice sheet. The basal warming towards the margin is due to the increased dissipation in this area. This state was used as the initial condition for a time-dependent computation where the sea-level temperature was raised by 10 K at t = 0₊. This represents the temperature change at the termination of the Last Glacial Maximum, allowing us to investigate whether warming can, after propagating to the base and making the ice more fluid, produce a grounding-line retreat. Figure 4b shows the internal temperature field after integration for 60 000 years, while Figure 4c shows the evolving basal temperature. The ice sheet thins by around 250 m (Fig. 4b), but grounding-line retreat is negligible. Substantial warming of the base and lowering (not shown) occur about 10 000 years after the warming at the upper surface has propagated to the base. (More realism would have included an increase in accumulation rate, but we wish to isolate the effects of different processes in this study.) This shows that the results of isothermal step-change calculations reported above, which predict very small grounding-line movement, represent a more general type of behaviour. Calculations with the rate factor set at the lower level used in the other computations reported in this paper showed the same pattern, but with some ice above the pressure-melting point.

Fig. 4. Temperature fields within a marine ice sheet and the evolution of basal temperature, non-oscillating case. (a, b) Ice-sheet geometry and temperature contours; (a) shows the initial steady condition, while (b) shows the computed geometry and internal temperature 60 000 years after the configuration in (a) had a surface warming of l0 K applied. (c) Evolving basal temperatures.

3.2.3. Oscillatory computations

We now ask whether free thermal oscillations can cause grounding-line retreat through a noise-induced drift phenomenon. We attempted to compute a steady state as in the previous experiment, but also included sliding with ℓ = 3, λ = 0 and A _s = 2.25. We computed a steady state for a sea-level temperature of −20°C (not shown) but were unable to find a steady state when the sea-level temperature was −30°C. As described above, to initialize the run the system was allowed to decay onto a steady oscillation with both equations in prognostic mode, but with the span fixed. (The features of this oscillation were similar to those computations we performed where grounding-line motion was permitted; see below.) Once the oscillation was steady, the span was permitted to evolve at time t = 0₊. Two computations were carried out: (i) a case where there was no warming, and (ii) a case where sea-level temperature was raised by 10 K at t = 0₊. Integrations were carried out for 20 000 years. In both cases, the oscillations were maintained as the span moved, and there was retreat of the grounding line. Such thermally driven oscillations have previously been modelled by MacAyeal (1992) and Payne (1995), neither of whom considered the effect on grounding-line motion.

The oscillation appeared to be of relaxation type, with a rapid warming, localized at a point on the flowline, which was driven backwards by the coupling between the dissipation and the ice-sheet profile. This occurred because the thinning associated with faster flow in the downstream area necessarily caused higher slopes further upstream. The faster flow and thinning advect cold ice nearer to the base, which is followed by basal cooling and freezing and the build-up of the ice sheet. We show the details of the calculations for the case with warming. Figure 5a–e show the evolution of the internal structure of the temperature field. The sequence starts immediately prior to the rapid phase. The slope is high at the margin, and the consequent increased heating leads to the development of fast flow shown in Figure 5b and c and the upstream propagation of the basal warm zone driven by increased dissipation, caused by the increased slope. Figure 5d shows the state immediately after the end of the fast-flow phase, where the ice sheet is low and cold because of the downward advection of cold ice, while Figure 5e shows the situation after the build-up of the ice sheet immediately prior to the next surge. The temperature inversions arise from cooling of surface ice as the elevation increased. Figure 6a and b show the evolution of the basal temperature profile for both cooling and non-cooling cases, which show (a) warm conditions associated with fast flow, and (b) the complex cooling patterns during the slow-flow phase. Figure 6c shows the evolution of the surface profile. Note how quickly the ice-sheet profile evolves during the fast-flow phase.

Fig. 5. Evolution of temperature in an ice sheet with internal, thermally generated oscillations. Ice-sheet geometry and temperature contours are shown. (a) is immediately prior to start of fast-flow phase; (b) and (c) show how the profile and temperature distribution change during the fast-flow phase; (d) shows state immediately after the end of the fast-flow phase; (e) shows the build-up phase. The sequence follows features noted by Payne (1995); heating propagates backwards in an area of greaterslope, and a cold plume of ice is “pulled” down by the increased flow. The evolving basal temperatures are shown in Figure 4.

Fig. 6. Evolution of basal temperature and surface profile for cases where the ice sheet is experiencing a thermo-viscous oscillation. (a) is a case where the surface was warmed by 10 K at t = 0₊; (b) is experiencing a constant climate; (c) shows the evolving surface profile for the warming case. The warming case is also illustrated in Figure 5.

A point of interest is that while a warming has virtually no effect on the span of a non-oscillating ice sheet (Fig. 4c), it has a significant effect on an oscillating ice sheet, causing faster retreat, possibly through its effect on the period of the oscillation. A further point of interest is that even though a steady state could be found for a sea-level temperature of −20°C (not shown), there seems to be no sign of the oscillations diminishing in the warming case, implying that there may be more than one possible state (steady and oscillating) for a given climate. We conclude that thermo-viscous oscillations can cause grounding-line retreat through a noise-induced drift phenomenon.

4. Discussion: Relationship to Previous Work

This work and the companion paper (Le Meur and Hindmarsh, 2001) address substantially the same issues as some previous modelling studies of marine ice sheets (Alley and Whillans, 1984; Van der Veen 1985, 1987). We term these models AW, VdV1 and VdV2, respectively. It is not straightforward to compare the results. Firstly, this paper considers that the whole issue of longitudinal stresses is of secondary importance where the freeboard is large; it cannot be overemphasized that we consider the exclusion of them a feature and not a failing of the model. The models AW, VdV1 and VdV2 consider longitudinal stresses at the grounding line, but with only a very coarse grid, where gradients at the grounding line in particular were unlikely to have been represented properly.

More pertinently, the experience of Hindmarsh (1993) shows that ice-sheet numerical models which correctly represent the zero-eigenvalue (neutral-equilibrium) property are not obtained without exercising some care. If this property is not obtained, the zero eigenvalue is approximated by a finite eigenvalue. This is equivalent to replacing an infinite time constant by a finite one. The corresponding mode is of long wavelength (Hindmarsh, 1996); in other words, we obtain a long-wavelength, stable or unstable mode with a long time constant. This is not a typical numerical instability and can appear deceptively similar to plausible physically based behaviour. It is quite likely that some or all of the models AW, VdV1 and VdV2 did not capture the neutral stability properties correctly. A diagnostic feature of this might be sensitivity of dynamics to the model formulation, as small errors in what should be the zero eigenvalue affect the spurious long time constant significantly. Thus, in Van der Veen (1985, figs 1 and 2), we find that the model which includes a longitudinal stress-gradient computation in the whole ice sheet, rather than in the transition zone only, is very much more sensitive, even though we know from the shallow-ice approximation that longitudinal stresses are irrelevant upstream of the grounding line. This effect is not shown in Van der Veen (1987). According to our analysis, the thinning associated with sea-level rise reported by Alley and Whillans (1984) and Van der Veen (1985) is a numerical artefact.

Van der Veen (1987) computed initial steady profiles by integrating from zero thickness. No mention is made in that paper of possible neutral equilibria. Van der Veen’s figures 5, 6 and 8 show ice-sheet profiles with steady span strongly dependent on the sliding model. The neutral-equilibrium hypothesis would argue that steady profiles are possible for all the spans, meaning that the differences between the Van der Veen (1987) steady spans are purely a function of history and not of physics. Of course, for a given span, the profile is affected by the basal physics.

Finally, Hindmarsh (1996) and Van der Veen (1999) both agree that the Thomas and Bentley (1978) model (which is not truly time-dependent) ignores the stabilizing influence of flux from the grounded area into the transition zone, which counteracts any destabilization from creep-thinning.

5. Conclusions

This paper has been concerned with understanding some general features of the dynamical processes which can cause retreat of marine ice sheets. A companion paper (Le Meur and Hindmarsh, 2001) shows that sea-level rise and the incorporation of bedrock dynamics cannot provide the mechanisms for observed marine-ice-sheet retreat. Ice-shelf disappearance as a cause of ice-sheet retreat is discounted on mechanical grounds. Step changes in accumulation rate or rate factor produce far larger changes in ice-sheet thickness (roughly compatible with scale-model predictions) than they do in span. Ice-stream models based on static hydrologies can produce retreat, and can produce substantial retreat much later than the initiation of sea-level rise, but are inconsistent with borehole observations. Internal noise sources (e.g. ice-stream fluctuations) can induce a noise-induced drift of the grounding line. The rate of drift could be high enough to be a significant factor in ice-sheet dynamics. Such drift arises essentially as a consequence of the asymmetric response of the ice sheet to forcing.

Warming an ice sheet in a way that would have occurred at the termination of the Last Glacial Maximum will thin the ice sheet, but does not in itself stimulate grounding-line retreat. Under certain parameter combinations internal thermally driven oscillations can promote retreat, and the combination of warming and internal oscillations promotes faster retreat. Internal oscillations do significantly affect grounding-line motion.

These results seem to suggest that while significant changes in thickness can be caused by spatially uniform changes, spatial variability coupled with dynamical variability is needed to cause significant margin movement. With drift mechanisms we seem to be able to cause a marine ice sheet to retreat as rapidly as Conway and others (1999) suggest it must, and we have no problem in principle in avoiding an early retreat associated with postglacial sea-level rise. The price for this success is that we probably need to invoke a change in ice-stream fluctuations sometime during the Holocene in order to stimulate the retreat, perhaps due to the postglacial warming reaching the base. We have not established with certainty that drift does play a significant role in marine ice-sheet dynamics, but if this is the case, processes with relatively short time-scales and spatial scales may ultimately control large-scale marine ice-sheet dynamics.

Many previous modelling efforts have been directed towards modelling the evolution of marine ice sheets and the WAIS in particular (e.g. Budd and Smith, 1981; Huybrechts, 1992), but none of them possesses an extremely important feature of the present model, which is that dynamical consistency has been shown to exist between the numerical model and the corresponding differential equations. Technically, the differential equations and difference schemes have been shown to be Lyapunov stable in the same way (Hindmarsh, 1996). It is the use of this model which has shown that the retreat of the WAIS is an issue in the sense that it is “large”, and it is not straightforward to explain. The retreat does not appear to have been caused directly by postglacial forcing and this makes predicting the future behaviour of the WAIS rather more difficult.

Acknowledgement

The authors would like to thank R. Greve for his thorough editing and his determination to make these papers readable.