Stability of ice-sheet grounding lines

Richard F. Katz, M. Grae Worster


Recent observations of the West Antarctic Ice Sheet document rapid changes in the mass balance of its component glaciers. These observations raise the question of whether changing climatic conditions have triggered a dynamical instability in the ice-sheet–ice-shelf system. The dynamics of marine ice sheets are sensitive to grounding-line position and variation, characteristics that are poorly captured by most current models. We present a theory for grounding-line dynamics in three spatial dimensions and time. Our theory is based on a balance of forces across the grounding line; it is expressed as a differential equation that is analogous to the canonical Stefan condition. We apply this theory to the question of grounding-line stability under conditions of retrograde bed slope in a suite of calculations with different basal topography. A subset of these have basal topography inspired by the Pine Island glacier, where basal depth varies in both the along-flow and across-flow directions. Our results indicate that unstable retreat of the grounding line over retrograde beds is a robust feature of models that evolve based on force balance at the grounding line. We conclude, based on our simplified model, that unstable grounding-line recession may already be occurring at the Pine Island glacier.

1. Introduction

The glaciers of the West Antarctic Ice Sheet comprise a volume of ice equivalent to a rise in sea level of approximately 3.3 m (Bamber et al. 2009). This ice sheet is, in large part, grounded on rock that is below sea level. Un-grounding of a relatively small portion of the ice sheet, and the attendant change in mass balance, could result in a significant change in global sea level. Indeed, recent geological observations suggest that the sea level rose by 2–3 m on a human time scale (approx. 100 years) during the terminal phase of the last interglacial period (Blanchon et al. 2009). Understanding the dynamics of Antarctic grounding lines and how they may respond to climatic perturbations is a critical component in predictions of future environmental change. A particularly important question regards the stability of grounding lines: can gradual modulation in the environmental conditions of Antarctica lead to massive, irreversible grounding-line recession and disintegration of the West Antarctic Ice Sheet? Novel theories of grounding-line dynamics may help us to answer this question and to interpret observations of rapid change in West Antarctica.

Recent observations of glaciers that flow into the Amundsen Sea suggest that they are undergoing accelerating change. Thomas et al. (2004) described altimetry measurements that document glacial thinning, and estimated that glaciers in this region are discharging almost 60 per cent more ice per year than they are accumulating. More recently, Wingham et al. (2009) documented rapid changes in ice thickness of the Pine Island glacier. Rignot (2008) used synthetic-aperture interferometry on radar observations collected over a period from 1974 to 2007 to determine changes in flow rates of the Pine Island and Smith glaciers. According to this work, the two glaciers sped up by 42 and 83 per cent, respectively, between 1996 and 2007; both also experienced a significant grounding-line recession during this period. The measurements show that acceleration of the Smith glacier was substantially larger in 2006–2007 than in previous years. Rignot et al. (2008) estimates that the combined flux of ice across the grounding line of the Pine Island and Thwaites glaciers in 2006 was 85±26 gigatonnes per year more than their combined ice accumulation rate at that time; this imbalance is a factor of two larger than that measured 10 years earlier. In Greenland, geochemical evidence reported by Briner et al. (2009) indicates a rapid deglaciation of the Sam Ford Fjord in the Early Holocene, with inferred topographic control of the rate of glacial retreat.

These, and other observations, again raise a crucial question (e.g. Weertman 1976; Bentley 1997): do recent changes represent a mounting, catastrophic collapse of the West Antarctic Ice Sheet associated with unstable retreat of grounding lines (e.g. Weertman 1974) or, rather, a steady and controlled adjustment to changing climatic forcing (e.g. Van der Veen 1985)? Since the dominant resistance to flow in ice sheets is provided by basal drag (e.g. Van der Veen 1999), a reduction in the grounded area of an ice sheet would plausibly lead to more rapid flow of ice. Ice-sheet buttressing by embayed ice shelves could also contribute a stabilizing resistance to flow; disintegration of an ice shelf would relieve this buttressing and allow for flow acceleration (Dupont & Alley 2005). While buttressing may be important to the overall dynamics of ice sheets, we neglect it in the current analysis and focus instead on the dynamical changes associated with grounding-line migration.

Although other authors had proposed the idea earlier, Weertman (1974) published the first mathematical theory for grounding-line instability. The theory is based on the analysis of a two-dimensional ice sheet and its junction with a two-dimensional ice shelf; the conditions applied at the junction are the flotation condition and continuity of ice thickness. Weertman (1974) concluded that glaciers grounded below sea level on beds that slope upward in the flow direction are inherently unstable to grounding-line recession (such beds are referred to here as being retrograde). Other authors have developed ice-sheet–ice-shelf models in which the grounding-line position is governed by the flotation condition or its total differential (e.g. Thomas & Bentley 1978; Hindmarsh 1996; Hindmarsh & LeMeur 2001; Vieli & Payne 2005; Pattyn et al. 2006). Comparative analysis of these models by Vieli & Payne (2005) suggested that differences in the predicted grounding-line motion could be attributed to issues with the numerical methods, as well as to differences in the formulation of grounding-line motion. This inaccuracy should not come as a surprise: as noted by Schoof (2007a) and Robison et al. (in press), solving for the position of the grounding line is a free-boundary problem (like the canonical Stefan problem) and, as such, it requires two independent conditions at the moving boundary. The flotation condition can only be used as one of these.

Chugunov & Wilchinsky (1996) and Schoof (2007b) applied asymptotic analyses to solve for the stress conditions in a boundary layer that connects the ice sheet to the ice shelf. In this narrow region, shear stresses (that dominate in the sheet) and longitudinal stresses (that dominate in the shelf) are of similar magnitude. Both Chugunov & Wilchinsky (1996) and Schoof (2007b) showed that by collapsing the width of this narrow transition zone to zero, the boundary layer can be approximated with a boundary condition on the ice sheet. There are several key difference between these two studies, however: Chugunov & Wilchinsky (1996) applied a no-slip boundary condition at the base of a Newtonian ice sheet and retained the vertical shear term in the horizontal momentum equation, while Schoof (2007b) used Glen’s flow law, applied a friction condition at the base and neglected vertical shear.

New work on a Newtonian ice-sheet model by Robison et al. (in press) began with the simplifying assumption that the transition zone between the ice sheet and the shelf can be collapsed to a boundary of zero thickness. The authors derived a grounding-line equation based on the balance of forces between the ice sheet and the ice shelf at this boundary. Steady-state solutions for the position of the grounding line were determined for downward sloping (prograde) beds. However, Robison et al. (in press) did not consider the existence or stability of steady solutions for upward sloping beds. Furthermore, none of the authors who have developed the free-boundary approach to deriving a governing equation for the grounding-line position have considered a three-dimensional ice sheet. It is thus the purpose of the present work to extend the theory of Robison et al. (in press) to three dimensions and, furthermore, to model ice-sheet flow and grounding-line position over a bed with topographic variation in both the along-flow and cross-flow directions.

The stability of a grounding line governed by the balance of forces was studied by Schoof (2007a). His work showed that, under the assumptions particular to his model (i.e. rapid sliding, no internal shear), there are no steady-state solutions for a two-dimensional ice sheet with the grounding line on an upward sloping bed. Work by Nowicki & Wingham (2008) using numerical solutions of the two-dimensional Stokes equation supports this result, but is inconclusive when a no-slip condition is applied at the base of the ice sheet. Like the case investigated by Schoof (2007a), non-slipping ice represents an end-member case that is worth investigating. While our analysis neglects basal sliding and non-Newtonian ice viscosity, it incorporates internal shear and three-dimensional effects that are absent from the work of Schoof (2007a).

In the next section of this paper, we provide details of the model derivation, scaling and numerical solution procedure. The theoretical development in this section is an extension of the ice-sheet and grounding-line model of Robison et al. (in press) from two to three dimensions. In the two-dimensional model, the dynamics of the ice shelf are described by the leading-order terms in lubrication theory: the stress is hydrostatic and the strain is purely extensional. For the general case of a three-dimensional ice shelf bounded by a curved grounding line, the dynamics are more complex. In the theory developed below, we ignore this additional ice-shelf complexity and adopt the two-dimensional theory without modification. This is a deficiency in our model and a target for future work.

Having described the theory, we present numerical solutions for various physical scenarios in §3. In two-dimensional model runs, we concentrate our attention on the stability of an ice sheet resting on the basal topography proposed by Schoof (2007a) that contains a continental shelf with an upward bulge (and retrograde slope) near its margin. We present results on three-dimensional model runs in §3b,c; the former considers the problem of steady-state grounding-line profiles for flow down shallow valleys, while the latter considers the dynamics and stability of grounding-line position over more complex basal topography. Section 4 contains a summary of results, plus further calculations relevant to the stability of the Pine Island glacier in West Antarctica. Conclusions are drawn in §5.

2. Model derivation

Following Robison et al. (in press), we consider an ice sheet of density ρ and viscosity μ flowing into sea water of density ρw and negligible viscosity. The ice sheet is grounded on a bed given by the surface z=−b(x,y,t), where x and y are horizontal coordinates, z is a vertical (upwards) coordinate and t is time. The time dependence of b will represent changes in sea level; we do not consider erosion or sedimentation processes, although these may be important (Alley et al. 2007). Figure 1 is a schematic diagram of the domain and coordinate system.

Figure 1.

Schematic of the model domain. Periodic boundary conditions are applied at the faces of the domain with the normal in the y-direction. The surface shown in grey scale is an example of the basal topography applied to the model. Solid black lines are depth contours; the dashed black line shows the intersection of basal surface with the edges of the domain.

The surface of the ice sheet is given by z=h(x,y,t). In the present work, we will be interested only in the dynamics of the ice sheet and hence will neglect the ice shelf, except in considering the force it exerts on the sheet at the grounding line. The grounding line has a position x=xG(y,t) that we require to be a single-valued function of y.

(a) The ice sheet

Force balance in the ice sheet is governed by leading-order terms of the shallow-ice equation (e.g. Fowler & Larson 1978; Hutter 1983; Mangeney & Califano 1998), giving a pressure within the ice sheet that is hydrostatic, p=ρg(hz). Horizontal pressure gradients are balanced by vertical shear according toEmbedded Image 2.1 where ν=μ/ρ is the kinematic viscosity, v=(u,v) is the horizontal vector flow velocity, g is the acceleration due to gravity and x is the horizontal gradient (the subscript is dropped below, but it is implied unless otherwise noted). In writing equation (2.1), we have neglected the downslope component of gravity, which restricts the validity of this equation to gently sloping beds with |b|≪1.

We apply no-slip and no-shear-stress boundary conditions at the bottom and top of the ice sheet, respectively,Embedded Image 2.2 andEmbedded Image 2.3 Integrating equation (2.1) in z twice, subject to these boundary conditions givesEmbedded Image 2.4 where H=h+b is the thickness of the ice sheet. A third integration in z gives the ice flux,Embedded Image 2.5 Using equation (2.5), the mass-balance equation can be written asEmbedded Image 2.6 where x⋅ is the horizontal divergence operator (again, the subscript is dropped below). This equation states that a change in the thickness of the ice sheet at a point in (x,y) is due to the divergence of the horizontal volume flux at that point. In writing equation (2.6), we have neglected the accumulation and ablation of ice, which are clearly important in large-scale models of ice sheets (e.g. Pollard & Deconto 2009).

We consider domains that are periodic in the y-direction, so equation (2.6) requires two boundary conditions. The first of these specifies a fixed flux at the inflow boundary,Embedded Image 2.7 while the second is the standard flotation condition at the grounding line,Embedded Image 2.8

As discussed in the introduction, the equations above do not provide enough information to determine the position of the grounding line xG. To do so, we require another condition at xG; this will take the form of a force balance between the ice sheet and the shelf. To obtain the force exerted by the ice sheet at the grounding line, we begin by considering stresses within the ice. Using the same arguments as Robison et al. (in press), we obtain the horizontal stresses asEmbedded Image 2.9 Embedded Image 2.10and Embedded Image 2.11 Integrating these stresses over the thickness of the ice sheet and using equation (2.4), we obtain the components of the vertically integrated horizontal stress tensor Nsheet as Embedded Image 2.12 Embedded Image 2.13and Embedded Image 2.14

(b) The ice shelf

Although the ice sheet described above has variations in both the x- and y-directions, for simplicity we adopt the same two-dimensional model of the ice shelf that was proposed by Robison et al. (in press). In this model, the ice shelf is floating in local hydrostatic equilibrium; the air above and water below impose no tangential stresses. The horizontal velocity of the ice is taken to be independent of depth. The depth-integrated longitudinal stress (viscous plus cryostatic) in the shelf is balanced by the hydrostatic pressure from the ocean underneath, acting on the sloped lower surface of the shelf. This depth-integrated stress is given by combining eqns (2.7) and (2.10) of Robison et al. (in press) asEmbedded Image 2.15 and is equal to the depth-integrated hydrostatic pressure of the ocean at the grounding line. We adapt the model by assuming that the calculated stresses are aligned in the direction normal to the edge of the ice shelf, at the grounding line; i.e. we take Embedded Image. The direction normal to the grounding line is defined below in equation (2.16).

It is important to note here that the use of a two-dimensional model for the ice shelf is strictly valid only when the grounding-line position is independent of y. This will not be the case in some of the calculations to follow. The assumption of a two-dimensional shelf is thus a deficiency of our model; addressing this deficiency is a target for future work.

(c) The grounding line

The equation governing the position of the grounding line will be determined by balancing stresses between the ice shelf and the ice sheet. The ice-sheet normal stress at the grounding line is Embedded Image, where n is a unit vector normal to the line xG(y,t) given byEmbedded Image 2.16 Here, ∂y is a partial derivative in the y-direction, and i and j are unit vectors in the x- and y-directions, respectively.

Expanding the ice-sheet normal stress at xG givesEmbedded Image 2.17 whereEmbedded Image 2.18

Assuming a balance between the forces normal to the grounding line from the ice sheet and the ice shelf (see above, §2b) givesEmbedded Image 2.19 In the limit Embedded Image, this equation reduces to the one-dimensional version,Embedded Image 2.20 which appears in Robison et al. (in press) as eqn (2.30).

Combining the mass-balance equation (2.6) with the time derivative of the flotation condition (2.8) gives an expression for the flux divergence,Embedded Image 2.21 where Embedded Image, and we have defined ρ′=ρw/ρ. Finally, this can be substituted into equation (2.19) to give the grounding-line equation,Embedded Image 2.22 at x=xG(y,t). I represents the identity matrix.

(d) Non-dimensionalization and change of coordinates

We define the following scales Embedded Image 2.23 Embedded Image 2.24and Embedded Image 2.24 where q0 is a constant flux with dimensions of area per time. These scales are used to non-dimensionalize variablesEmbedded Image 2.26 and Embedded Image. In terms of dimensionless variables, the grounding-line equation becomesEmbedded Image 2.27 where Embedded Image is the ice thickness at the grounding line. Below, we drop hats from dimensionless variables.

To facilitate discretization of the problem on a regular grid, we map the domain to a rectangular one using the transformationEmbedded Image 2.28 where ξ goes from zero at x=0 to one at x=xG. With this change,Embedded Image 2.29 Embedded Image 2.30and Embedded Image 2.31 where xG′ is the y-derivative of the x-position of the grounding line.

With dimensionless variables and transformed coordinates, the mass-balance equation (2.6) becomesEmbedded Image 2.32 whereEmbedded Image 2.33 The boundary conditions (2.7) and (2.8) becomeEmbedded Image 2.34 andEmbedded Image 2.35

The grounding-line equation (2.26) was transformed similarly in the numerical code.

(e) Numerical solution

We discretize the governing equations with implicit time stepping and a flux-conservative, finite-volume scheme. The mass-balance equation (2.32) is discretized onto a Cartesian two-dimensional grid with ni points in the ξ-direction and nj points in the y-direction; the grounding-line equation is discretized onto a one-dimensional grid with nj points. This yields a system of (ni+1)×nj nonlinear algebraic equations in (ni+1)×nj variables. Experience has shown us that a solution strategy based on iterative, separate updates of the grounding-line position and the ice-sheet height at each time step is numerically unstable. To avoid this problem, we invert the entire, coupled system for {hij,xGj} at each time step. This inversion is accomplished using a Newton–Krylov method (Generalised Minimum Residual), preconditioned with an incomplete Lower–Upper matrix factorization; both are provided by the Portable Extensible Toolkit for Scientific Computation (PETSc; Balay et al.2001, 2004).

The calculations are initialized with an asymptotic solution to the two-dimensional governing equations at t≪1 (Robison et al. in press). For two-dimensional calculations (e.g. figure 2a), we use 2000 grid cells, evenly distributed in the ξ-direction. For three-dimensional calculations, unless noted in the figure caption, we use a grid of 120×120 cells, evenly distributed in ξ,y; the grounding-line position is discretized over the same number of cells as in the y-direction of the ice sheet. The time-step size is adjusted at each step to give a maximum proportional change in the grounding-line position of approximately 0.05 per cent. We accept the solution at any time step when the L2 norm of the nonlinear system residual is less than 10−8. The code has been verified in the simple case of a two-dimensional ice sheet with constant bed slope by comparison with an analytical solution from Robison et al. (in press). This work demonstrated that the theory admits a single solution for the grounding-line position on any contiguous section of prograde bed. As a qualitative check for the case of more complex bed topography, we compare the results with the model of Schoof (2007a) in the next section.

Figure 2.

Hysteresis in a two-dimensional ice-sheet model with the changing sea level (SL). (a) Steady-state curves h(x) representing the height of the ice sheet in metres for different values of sea level. The dashed curves represent a series in which the sea level is sinking; the solid curves represent a series with the sea level rising. The bedrock profile is as used by Schoof (2007a). (b) The hysteresis loop between the position of the grounding line xG and sea level. Calculations shown in (a,b) were performed with parameters μ=1×1013 Pa s, ρ=900 kg m−3, ρw=1000 kg m−3 and q0= 0.02 m2 s−1. (c) The changes in the hysteresis loop for different values of q0. Right-pointing open triangles mark the value of sea level at which the grounding line advances; left-pointing filled triangles mark the value at which the grounding line recedes.

3. Results

(a) Stability in two-dimensional models

It is instructive to consider grounding-line instability in comparison to the results of Schoof (2007a). To do so, we adopt the profile of basal topography proposed in eqn (10) of that paper (shown in our figure 2a). This basal topography profile mimics the form of many of the ice-sheet beds around western Antarctica: a long, shallow, prograde continental slope curves upward into a shallower outer rise prior to dropping sharply into the deep ocean. The presence of retrograde slopes on the upstream side of this outer rise creates a region where grounding lines are expected to be unstable. Our results affirm this expectation.

 Figure 2a shows equilibrium profiles of ice-sheet surface height at various sea-level heights. Each curve was calculated by stepping sea level in increments (approx. 1 m) from the previous value and allowing the height and grounding-line position to evolve according to the governing equations until a new steady state was reached. The critical step over which unstable grounding-line recession occurred was refined to 0.4 mm change in sea level.

Like Schoof (2007a), our model produces a hysteresis loop, shown in figure 2b, for grounding-line position as a function of sea-level height. With the grounding line on the outer prograde slope and a rising sea level, the grounding line recedes in small steps, proportional to the change in sea level. When it reaches the crest of the outer topographic rise, an incremental increase in sea level leads to a large leap backwards in grounding-line position. The grounding line settles on the inner prograde slope at a position where the basal depth is approximately equal to that of the crest of the outer rise. Further small increases in sea level lead to proportional recession of the grounding line up the inner slope. As sea level is dropped in small steps, a reversal of this process occurs, except that the unstable leap in grounding-line position occurs when the grounding line has reached the trough between the inner slope and the outer rise. There, the grounding line leaps forward to the prograde outer slope.

 Figure 2c demonstrates that, in calculations with different values of the input flux q0, the unstable grounding-line transitions occur at different values of sea level. A larger flux of ice results in a thicker ice sheet and therefore requires a higher sea level for both prograde and retrograde transitions of the grounding-line position. More importantly, for a fixed sea level, a transition in grounding-line position can result from a change in the ice flux. Although our model does not explicitly include basal melting of ice at the grounding line as does Walker et al. (2008), this melting may be similar in effect to a reduction in flux q0. As figure 2c shows, such a reduction can lead to unstable grounding-line retreat.

To investigate the dynamics of a retrograde transition event, we consider a fixed flux q0 and begin the calculation with a steady-state ice-sheet profile that is on the verge of a transition. A step increase in sea level is applied and the system is allowed to evolve according to the governing equations. Figure 3 shows the results of an example calculation with q0=0.02 (other parameters as in figure 2a). In this case, the step in sea level triggered a retrograde transition that shifted the grounding line xG by almost 500 km. Figure 3a shows the grounding-line position and the ratio of the flux across the grounding line to that entering the domain qG/q0, as a function of time, from the time of the change in sea level. The first phase of the transition is slow owing to the shallow retrograde slopes at the crest of the outer rise. As the grounding line reaches the steeper retrograde slope of the outer rise, the transition enters the second phase. In this phase, the basal depth at the grounding line increases, leading to an accelerating change in ice flux across the grounding line and a thinning of the ice sheet, which, in turn, causes further grounding-line recession. Figure 3b shows that qG/q0 reaches a maximum when the grounding line is at the deepest point between the inner prograde slope and the outer rise. In phase 3 of the transition, the grounding line climbs the inner slope. As the basal depth at the grounding line decreases, so too does the flux ratio; the transition ceases as the flux ratio approaches unity.

Figure 3.

Ratio of grounding-line flux qG to influx q0 over an episode of unstable grounding-line recession caused by a step change in sea level. When qG/q0=1, the ice sheet is in a steady state. (a) The grounding-line position (dashed line) and the flux ratio (solid line) as a function of time t since the step change in sea level. The three phases of the transition are delimited above the panel (see main text for explanation). (b) The basal topography (dashed line) and the flux ratio (solid line) as a function of horizontal distance x. The grounding-line position is initially at the right end of this plot; with time it moves to the left. The parameters for this calculation are as stated in the caption to figure 2. The step change in sea level was 5 m.

For the basal topographic profile adopted from Schoof (2007a), the rapid phases (2 and 3) of the transition in grounding-line position occur over about 100 kyr, a period that is nearly independent of the influx q0 and the step size in sea level, for values considered. The maximum flux ratio during the transition decreases by about 10 per cent with changing background flux q0, from 1.57 at q0=0.02 to 1.42 at q0=0.045. These results are specific to the basal profile proposed by Schoof (2007a) and used here (as well as the chosen model and parameter values). In §4b, we demonstrate that a different profile of basal topography yields a grounding-line transition that is significantly more rapid. Beyond this dependence, the dynamical time scale associated with unstable grounding-line recession will be affected by the rheology of the ice and the basal conditions (Joughin et al. 2009). We expect that our use of Newtonian rheology and a no-slip condition at the base of the grounded ice yields an overestimate of the dynamical time scale for recession.

(b) Steady states in three dimensions

The governing equations derived above permit the calculation of grounding-line evolution in three dimensions. With the addition of the across-flow (y) direction, we have the opportunity to specify variations in basal topography in both x and y. The calculations presented in this section use a particularly simple basal topography: a constant, prograde slope that is modulated by a sinusoidal function in the y-direction,Embedded Image 3.1 Here, α is the mean slope in the x-direction, β is the amplitude of modulation in slope and k is the wavenumber of that modulation in the y-direction; variables in equation (3.1) are dimensional. The ratio B=β/α gives the proportional variation of basal depth at fixed x relative to the mean depth at that distance. Our assumption that the downslope component of gravity is negligible (see above) requires that we consider only α≪1, which is an acceptable restriction because this is uniformly true for natural ice sheets.

 Figure 4 shows an example of a calculation with mean slope α=3.3×10−3, giving about 3 m of basal-depth change per kilometre in the x-direction, and β/α=0.8, giving up to 80 per cent deviation in basal depth at any value of x. Setting Embedded Image gives a single valley with a width equal to the width of the domain. Although this bedform does not correspond with a specific natural ice sheet, it is reminiscent of the Pine Island glacier, where the basal topography confines much of the flow to a broad ice stream underlain by a basal valley (e.g. Vaughan et al. 2006).

Figure 4.

Results from an example calculation with variations in basal depth in the y-direction as described by equation (3.1), with α=3.3×10−3 and β= 2.7×10−3. (a) The top and bottom surfaces of the ice sheet in its steady-state configuration. Sea level is at zero on the z-axis. The ice shelf is not shown. (b) The grounding-line position (black curves) at 50 year intervals, as it approaches a steady state. Coloured curves are contours of basal depth. (c) The magnitude of the vertically integrated, steady-state ice flux. The black curves are streamlines. Parameters for this calculation are as follows: μ=1×1013 Pa s, ρ=900 kg m−3, ρw=1000 kg m−3 and q0= 0.04 m2 s−1.

 Figure 5a shows the characteristics of the ice flux at the grounding line for the calculation depicted in figure 4. The grounding line is curved, but it does not follow a depth contour (figure 4b). Hence, the basal depth at the grounding line varies substantially, although less than it would if xG were a constant, independent of y. Given this variation in basal depth, it is not surprising to find that the ice flux across the grounding line is concentrated into the centre of the valley. Figure 5a shows that, for this value of β/α, there is almost zero flux at the edges of the domain, along the ridges that flank the valley (recall that the domain is periodic in the y-direction).

Figure 5.

Characteristics of the ice flux at the grounding line, x=xG, for the simulation shown in figure 4 (β/α=0.8). (a) The flux normal to the grounding line, and the flux in the x- and y-directions (narrow lines, y-axis at left), as well as the depth of the ice-sheet base below sea level (wide line, y-axis at right). (b) The flux divided by the cube of the grounding-line ice thickness Embedded Image and normalized by the maximum value of the qn curve. Using equation (2.5), we can deduce that the upstream dynamics of ice flow influence the gradient in surface height at the grounding line, and thus influence the flux at the grounding line. The asymmetry in the qn curve about y=25 km is a numerical artefact.

The components of the ice flux at the grounding line are plotted in figure 5b. Because we expect a strong variation in the flux with HG (equation (2.5)), the curves are scaled by the cube of the ice thickness there and then normalized. The heavy curve representing Embedded Image has little variation across the centre of the grounding line, between y≈5 and y≈45 km, whereas outside of this region, it decays sharply towards zero. This indicates that focusing of flow redistributes ice towards the centre of the valley, and that the flow along the ridge crests is anomalously sluggish relative to the thickness of ice there. Unsurprisingly, the flow has a three-dimensional aspect that is not negligible.

To explore the sensitivity of three-dimensional models to basal topography, we have performed a suite of calculations for different values of the mean bed slope α and the sinusoidal perturbation to the bed slope β; these are summarized in figure 6. Figure 6a shows how the mean grounding-line position (averaged over the y-direction) depends on the mean bed slope. There is a linear relationship between the inverse slope 1/α and the mean grounding-line position. This relation was predicted for the two-dimensional case (β=0) by Robison et al. (in press). Our results demonstrate that the mean grounding-line position is somewhat sensitive to the perturbation in bed slope; larger perturbations yield smaller mean values of xG. Figure 6b shows the difference between the largest and the smallest value of xG as a function of the relative perturbation in bed slope. As expected, this difference grows with increasing β/α. It is largest for the most steeply sloping bed, where a given amount of relative perturbation translates to the largest amount of actual bed-depth variation at a fixed distance x. Finally, figure 6c shows how the distribution of flux across the grounding line varies with increasing bed-slope perturbation. As shown in figures 4c and 5a, ice flux concentrates above the valley. With increasing β/α, the ratio of the largest flux to the smallest flux increases, reflecting increased focusing of ice flux towards the centreline in y where the bed is deepest. For the largest value of α, the grounding line is too close to the constant-influx boundary for efficient focusing to occur.

Figure 6.

Variations in ice-sheet characteristics with α and β. (a) The mean distance to the grounding line as a function of the inverse of the mean slope of the bed 1/α. Steeper slopes are at the left of the diagram. Different lines represent different values of the relative perturbation in bed slope β/α. The grey line is the analytical solution of Robison et al. (in press). (b) The difference in distance between the furthest point on the grounding line and the nearest as a function of the perturbation in bed slope. A larger slope perturbation typically results in a more strongly curved grounding line. Different lines represent different values of α. (c) The ratio of maximum (max) to minimum (min) flux in the x-direction across the grounding line as a function of the perturbation in bed slope. The different lines are as in (b). In (b,c), α represents kilometre height per kilometre x-distance. (a) Shaded thick line, β/α=0; circle with continuous line, β/α=0.1; circle with dashed line, β/α=0.2; circle with dashed dotted line, β/α=0.4; circle with dotted line, β/α=0.6; diamond with continuous line, β/α=0.8. (b,c) Circle with continuous line, α=8.3×10−4; circle with dashed line, α=1.7×10−3; circle with dashed dotted line, α=3.3×10−3; circle with dotted line, α=1.3×10−2; diamond with continuous line, α=5.3×10−2.

(c) Stability in three-dimensional models

In the previous sections, we examined the unstable behaviour of two-dimensional ice sheets with oscillatory topography (i.e. an outer rise) and the stable behaviour of three-dimensional ice sheets flowing over topography that is uniformly prograde in the x-direction. Next, we consider a three-dimensional scenario with basal topography that is a combination of these two: it is an internal valley, bounded in the y-direction and in the x-direction, that is superposed on a gentle, monotonic slope in x (colour contours in figure 7). In this case, there is a competition between the monotonic prograde slope, which would yield stable grounding-line recession, and the internal valley bounded by a retrograde slope, which would yield unstable grounding-line recession. Such a configuration is a highly idealized version of the basal topography beneath the Pine Island glacier (Vaughan et al. 2006) and is hence of special relevance here. Its behaviour cannot be deduced from the calculations described above.

Figure 7.

Grounding-line positions at steady state for two different basal topography maps that differ only in their value of σ from equation (3.2). The background slope α is equal to 0.52 m km−1 in both cases. Each grounding line, shown in black, is a steady-state solution for a different value of sea level. (a) A narrow central valley (σ=10 km). Sea level rises in increments of 3.2 m. (b) A wider central valley (σ=20 km). Sea level rises in increments of 6.3 m. The grounding line recedes unstably in this case, while it recedes stably in (a). An examination of the time-dependent transition between grounding lines in (b) indicates that the difference in behaviour between (a,b) is not due to the difference in sea-level increments. The grid size was 102×102 in these calculations.

We parametrize the above configuration with the dimensional functionEmbedded Image 3.2 where x, y and b are measured in metres, G(y,σ) is a Gaussian centred at the domain’s midpoint in y with unit amplitude and standard deviation σ, and TS is the topographic profile proposed by Schoof (2007a) and used above in §3a.

 Figure 7 shows results from two calculations of steady-state grounding-line position for a range of sea-level values. With rising sea level, the grounding line retreats towards the left in each panel of the figure. The basal topography in these two calculations differs only in the width of the central valley, given by the standard deviation σ of the Gaussian in equation (3.2). The difference in behaviour of the grounding line is striking, however. In the case of the narrower central valley in figure 7a, there are steady, stable grounding-line positions throughout the range in x between the valley’s trough and outer rise. In contrast, figure 7b shows that, for a wider valley, the grounding line undergoes unstable recession through this same range in x. This is the same instability as described in §3a, but now mediated by three-dimensional flow.

A set of calculations similar to (and including) those shown in figure 7 are summarized in figure 8. Here, the grounding-line position has been averaged in the y-direction and plotted on the x-axis. The independent variable, sea level, is plotted on the y-axis. With increasing time, the curves trace a trajectory that goes from the right-hand side of the figure to the left. Only the grey curves (α=0.52 m km−1) with narrow central valleys (σ=5 and σ=10) show stable grounding-line recession. A narrower valley drains the grounded ice less efficiently as sea level rises, reducing its destabilizing effect. Similarly, for larger α, grounded ice flows more rapidly on the monotonic slope outside the central valley. This diminishes the relative contribution of drainage through the valley and hence diminishes its destabilizing effect on the grounding line.

Figure 8.

Summary of grounding-line position (averaged in the y-direction) as a function of sea level for two sets of calculations with different values of α. Mean grounding-line position Embedded Image is on the x-axis and the independent variable, sea level, is on the y-axis. Each curve is labelled with the corresponding value of σ from equation (3.2). The grey curves labelled 10 and 20 are derived from the calculations shown in figure 7. Black line, α=0.17 m km−1; grey line, α=0.52 m km−1.

4. Discussion

(a) Review of the new model

In the foregoing sections, we describe a model of ice-sheet flow and grounding-line dynamics extended from that of Robison et al. (in press). Our model differs from that of Schoof (2007a) in that we neglect basal sliding of the glacier, assume a constant viscosity for the ice and account for vertical shear stresses. We collapse the transition zone between grounded and ungrounded ice proposed by Chugunov & Wilchinsky (1996) and Schoof (2007b) to a single interface on which normal forces between the ice sheet and the ice shelf are assumed to be in balance.

The results presented above support the theory of unstable recession of ice sheets grounded on retrograde beds (Weertman 1974). In this paper, our focus is on the effect of changes in sea level, but in our model, changes in the influx of ice can also trigger instability. Such flux changes may be qualitatively interpreted as changes in the accumulation rate of ice on the glacier surface, or as changes in the ablation rate of ice at or upstream from the grounding line; both of these have been documented on the West Antarctic Ice Sheet (Rignot et al. 2008).

Our predictions of unstable grounding-line retreat are similar to the findings of Schoof (2007a). These results reinforce the conclusion that, for ice-sheet/shelf models which enforce boundary conditions at the grounding line, the position of this line should be modelled as a free boundary, and it should be determined with a dynamic boundary condition (e.g. equation (2.22)) similar to that used in Stefan-type solidification problems. In contrast, models that solve the full Stokes equation without vertical integration, throughout a unified domain comprising of the ice sheet and the shelf, do not require boundary conditions at the grounding line; in that case, the conclusion above is not applicable. These latter type of models could be used as an independent test of our results and conclusions. Below, we discuss the implications of our model for understanding the current behaviour of the Pine Island glacier.

(b) Stability of Pine Island glacier

The theory and calculations described above are idealized in order to highlight the fundamental behaviour of the grounding line; we assumed Newtonian viscosity, one-dimensional ice-shelf flow and smoothed basal topography. With these assumptions, the model is not well suited to a detailed study of the dynamics of the Antarctic ice cap as a whole, or that of its component glaciers. However, it is interesting to consider calculations that incorporate a basal topography map inspired by observations of the Pine Island glacier. These observations are available thanks to the Airborne Geophysical Survey of the Amundsen Embayment (AGASEA) project (; Vaughan et al. 2006).

 Figure 9a shows a map of the basal topography beneath and around the Pine Island glacier. A transect along the approximate flow axis of the glacier is given in figure 9b. A cubic polynomial was fit to the transect using least squares and became the basis for two-dimensional calculations of the grounding-line motion described below (details are given in the caption of figure 9). To constrain the influx parameter q0, we take the estimate of Rignot (2008) for the average rate of ice input into the Pine Island glacier over the years 1980–2004 of 61±9 gigatonnes per year; assuming a grounding-line length of 50 km and an ice density of 900 kg m−3, this gives an average flux of ice across the grounding line of q0≈0.04 m2 s−1. As above, we assume a Newtonian ice viscosity of μ=1×1013 Pa s. To generate an initial condition for instability calculations, we adjusted sea level until the calculated steady-state grounding-line position was located just downstream of the crest of the outer rise (approx. 360 km in figure 9b). A different initial condition was determined for each relevant value of q0. The instability calculations were then initialized using these steady-state solutions, but with sea level raised by several centimetres from the value used to generate the initial condition.

Figure 9.

The Pine Island glacier as observed by the AGASEA project (Vaughan et al. 2006). (a) A map of basal topography with respect to the current sea level. The dashed line is the current grounding line from data by the National Snow and Ice Data Center (Haran et al. 2005). The black line is a transect through the map that intersects the grounding line at the point marked with a red star. (b) The basal topography −b, ice thickness H and calculated surface height h along the transect from A to A′ in (a). Where the ice is grounded, h=Hb; where the ice is floating h= (ρwρ)/ρwH with ρ=900 and ρw=1000 kg m−3. The blue curve is a polynomial fit given by y(x)=−586−15.1x+(7.89×10−2)x2−(1.01×10−4)x3, where x is distance in kilometres and y is the height above the current sea level in metres.

Output from instability calculations is shown in figure 10. In that figure, we compare the predicted flux ratio qG/q0 as a function of time with data for the Pine Island glacier (the data is derived from table 1 of Rignot (2008)). The data points form an array that bends sharply upward, mirroring the behaviour of the curves. The curves, however, are a poor fit to the slope of the data array. This may indicate that the Pine Island grounding line is presently undergoing unstable recession, but at a rate that is much more rapid than our simple models predict. This discrepancy may be due to our modelling assumptions (e.g. Newtonian viscosity, no basal slip), or may be due to a difference in the forcing that is causing recession. For the model curves, the forcing is a small step in sea level that triggers a change from one stable state to another. Over the time range of the data series from Rignot (2008), approximately 30 years, the global mean change in sea level was approximately 6 cm (Solomon et al. 2007). It is possible that the system crossed a threshold that triggered instability at some sea-level point within this range. It is more likely, however, that the natural system is being forced by a larger single perturbation such as the change in sea level since the last glacial maximum, or by a set of perturbations including reduced ice accumulation, increased ablation at the grounding line, more basal melt water, warmer ice and sea-level change. The comparison between data and the model indicates that the current rate of recession of the Pine Island grounding line is well above what might be considered a minimum estimate from our models.

Figure 10.

Flux ratio qG/q0 as a function of time. The lines represent the flux ratio calculated as in figure 3, except with basal topography as given by the blue curve in figure 9b, for four different values of q0. Data points are derived from the ratio of glacial outflow to the mean accumulation rate of the Pine Island glacier, as reported in table 1 of Rignot (2008). Error bars represent uncertainty propagation from the data in Rignot (2008). The curves and data are shifted in time such that the grounding-line position xG≈260 km (figure 9b) occurs at the present time t=0. Open circles, data from Rignot (2008); blacked dashed line, model, q0=0.02 m2 s−1; solid black line, model, q0=0.04 m2 s−1; grey dashed line, model, q0=0.06 m2 s−1; grey solid line, model, q0=0.08 m2 s−1.

Given the complex, three-dimensional nature of the real Pine Island glacier, with its convergent feeder streams and subglacial hydrology, it should be clear that the above model is a very crude representation of reality. Our current model implementation is incapable of capturing this complexity, although in §3c, we presented the results for a basal topography inspired by that of the Pine Island glacier. Figure 10 demonstrates that the Pine Island basal profile yields an instability in two-dimensional models that is qualitatively identical to, but more rapid than, that associated with the basal topographic profile introduced by Schoof (2007a). Incorporation of this profile into a three-dimensional model would mean modification of the basal profile of the central valley, TS(x) in equation (3.2). Based on the basal topography shown in figure 9a, this model would use σ≈20 km and α<0.17 m km−1. Judging from the behaviour of curves in figure 8, a model with these parameter values would produce a result similar to those in figure 7b: unstable grounding-line recession. However, owing to the complexity of the natural Pine Island glacier, these results would not provide additional insight over the two-dimensional models. In fact, similar patterns between data and models (figure 10) suggest that our simple two-dimensional models may capture an essential aspect of the Pine Island glacier.

5. Conclusions

The theory developed in §2 can be used to describe the dynamics of ice sheets and their grounding lines in three dimensions. We have applied this theory to develop models of grounding-line stability and instability in two and three dimensions. In particular, inspired by the questions posed by early investigators (e.g. Weertman 1974; Thomas & Bentley 1978) and by recent observations (Blanchon et al. 2009), we have considered the instability of marine ice sheets grounded on a bed with an outer rise separating the interior of the ice sheet’s bed from the steep continental slope. As a benchmark of our model, we have considered a basal profile adopted from Schoof (2007a) and calculated the response of an idealized ice sheet to changes in sea level. Our results are similar to those of Schoof (2007a), and again confirm the hypothesis of unstable grounding-line recession over retrograde beds (Weertman 1974).

We have also examined the stable shape of grounding lines for three-dimensional ice sheets with bed topography that is variable in two dimensions. The results of this work show that ice-sheet dynamics play a significant role in determining the position and shape of the grounding line, and the ice flux across it; these features are not simply related to the map of basal elevation.

Three-dimensional models were then extended to a case in which an internal valley with a retrograde outer rise is superposed on a prograde bed of constant slope. These calculations showed that the entire ice sheet can be destabilized by the presence of the internal valley, even if it is narrow with respect to the width of the domain. These results bear on sections of an ice sheet affected by strong variations in basal topography, such as the Pine Island glacier.

Calculations of grounding-line recession over a profile of basal topography derived from AGASEA data (Vaughan et al. 2006) were considered in order to compare the ratio of grounding-line ice flux with accumulation or influx of ice. Data from Rignot (2008) allowed us to compare model results with the observed behaviour of the Pine Island glacier. The data follow a trajectory with a rate of change that is larger than model curves. If the Pine Island glacier is indeed undergoing unstable grounding-line recession, the discrepancy between the model and data may be due to two key factors that will affect the dynamical adjustment time scale in our model: our use of Newtonian viscosity, and our application of zero basal sliding. The latter, especially, is known to be important for the dynamics at Pine Island (Joughin et al. 2009).

These results reinforce the conclusion that grounding-line dynamics are a critical component of dynamical models of the West Antarctic Ice Sheet. While most continental-scale models of ice sheets still use kinematic (and potentially inconsistent; Vieli & Payne 2005) boundary conditions at the grounding line, the current results indicate that the grounding line should be treated as a free boundary, analogous to the solidification front in a Stefan problem.

The results also underscore the importance of continued research on ice-sheet dynamics. This research should include the development of models that incorporate the detailed bedform and rheology of the West Antarctic Ice Sheet, such as the work of Pollard & Deconto (2009). It should also include the development of idealized models, such as those presented here, that explore the consequences of model formulation and assumptions, and probe the basic physics of ice sheets. In particular, more work is needed to predict the time scale of unstable grounding-line retreat.

Finally, our results suggest that, in contrast to earlier assessments (e.g. Van der Veen 1985; Vaughan & Spouge 2002), the scenario of unstable grounding-line recession on retrograde beds in West Antarctica is likely. Indeed, in the case of the Pine Island glacier, it may be presently occurring.


R.F.K. is grateful to R. Hindmarsh and R. Peltier for discussions and encouragement, to O. Sergienko for comments on an early version of the manuscript, to B. Smith and the PETSc team at Argonne National Laboratory for PETSc support and enhancements, and to B. Raup and J. Bohlander from the NSIDC for help with Antarctic grounding-line data. The authors acknowledge helpful reviews by R. Gladstone and an anonymous reviewer. R.F.K. received support for this work under award 0602101 from the US National Science Foundation International Research Fellowship Program, and as an Academic Fellow of the Research Councils UK.


    • Received August 17, 2009.
    • Accepted December 2, 2009.


View Abstract