## Abstract

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 (2007*a*) 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 (2007*b*) 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 (2007*b*) 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 (2007*b*) 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 (2007*a*). 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 (2007*a*), 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 (2007*a*).

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 (2007*a*) that contains a continental shelf with an upward bulge (and retrograde slope) near its margin. We present results on three-dimensional model runs in §3*b*,*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.

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*=*x*_{ G}(*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*(*h*−*z*). Horizontal pressure gradients are balanced by vertical shear according to
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,
2.2
and
2.3
Integrating equation (2.1) in *z* twice, subject to these boundary conditions gives
2.4
where *H*=*h*+*b* is the thickness of the ice sheet. A third integration in *z* gives the ice flux,
2.5
Using equation (2.5), the mass-balance equation can be written as
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,
2.7
while the second is the standard flotation condition at the grounding line,
2.8

As discussed in the introduction, the equations above do not provide enough information to determine the position of the grounding line *x*_{ G}. To do so, we require another condition at *x*_{ G}; 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 as
2.9
2.10and
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 **N**^{sheet} as
2.12
2.13and
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) as
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 . 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 , where **n** is a unit vector normal to the line *x*_{ G}(*y*,*t*) given by
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 *x*_{ G} gives
2.17
where
2.18

Assuming a balance between the forces normal to the grounding line from the ice sheet and the ice shelf (see above, §2*b*) gives
2.19
In the limit , this equation reduces to the one-dimensional version,
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,
2.21
where , and we have defined *ρ*′=*ρ*_{w}/*ρ*. Finally, this can be substituted into equation (2.19) to give the grounding-line equation,
2.22
at *x*=*x*_{ G}(*y*,*t*). **I** represents the identity matrix.

### (d) Non-dimensionalization and change of coordinates

We define the following scales
2.23
2.24and
2.24
where *q*_{0} is a constant flux with dimensions of area per time. These scales are used to non-dimensionalize variables
2.26
and . In terms of dimensionless variables, the grounding-line equation becomes
2.27
where 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 transformation
2.28
where *ξ* goes from zero at *x*=0 to one at *x*=*x*_{ G}. With this change,
2.29
2.30and
2.31
where *x*_{ G}′ is the *y*-derivative of the *x*-position of the grounding line.

With dimensionless variables and transformed coordinates, the mass-balance equation (2.6) becomes 2.32 where 2.33 The boundary conditions (2.7) and (2.8) become 2.34 and 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 *n*_{i} points in the *ξ*-direction and *n*_{j} points in the *y*-direction; the grounding-line equation is discretized onto a one-dimensional grid with *n*_{j} points. This yields a system of (*n*_{i}+1)×*n*_{j} nonlinear algebraic equations in (*n*_{i}+1)×*n*_{j} 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 {*h*_{ij},*x*_{ G}_{j}} 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 2*a*), 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 (2007*a*) in the next section.

## 3. Results

### (a) Stability in two-dimensional models

It is instructive to consider grounding-line instability in comparison to the results of Schoof (2007*a*). To do so, we adopt the profile of basal topography proposed in eqn (10) of that paper (shown in our figure 2*a*). 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 2*a* 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 (2007*a*), our model produces a hysteresis loop, shown in figure 2*b*, 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 2*c* demonstrates that, in calculations with different values of the input flux *q*_{0}, 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 *q*_{0}. As figure 2*c* 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 *q*_{0} 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 *q*_{0}=0.02 (other parameters as in figure 2*a*). In this case, the step in sea level triggered a retrograde transition that shifted the grounding line *x*_{ G} by almost 500 km. Figure 3*a* shows the grounding-line position and the ratio of the flux across the grounding line to that entering the domain *q*_{G}/*q*_{0}, 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 3*b* shows that *q*_{G}/*q*_{0} 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.

For the basal topographic profile adopted from Schoof (2007*a*), 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 *q*_{0} 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 *q*_{0}, from 1.57 at *q*_{0}=0.02 to 1.42 at *q*_{0}=0.045. These results are specific to the basal profile proposed by Schoof (2007*a*) and used here (as well as the chosen model and parameter values). In §4*b*, 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,
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 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 5*a* 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 4*b*). Hence, the basal depth at the grounding line varies substantially, although less than it would if *x*_{ G} 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 5*a* 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).

The components of the ice flux at the grounding line are plotted in figure 5*b*. Because we expect a strong variation in the flux with *H*_{G} (equation (2.5)), the curves are scaled by the cube of the ice thickness there and then normalized. The heavy curve representing 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 6*a* 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 *x*_{ G}. Figure 6*b* shows the difference between the largest and the smallest value of *x*_{ G} 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 6*c* shows how the distribution of flux across the grounding line varies with increasing bed-slope perturbation. As shown in figures 4*c* and 5*a*, 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.

### (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.

We parametrize the above configuration with the dimensional function
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 *T*_{S} is the topographic profile proposed by Schoof (2007*a*) and used above in §3*a*.

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 7*a*, there are steady, stable grounding-line positions throughout the range in *x* between the valley’s trough and outer rise. In contrast, figure 7*b* 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 §3*a*, 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.

## 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 (2007*a*) 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 (2007*b*) 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 (2007*a*). 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 (http://www.ig.utexas.edu/research/projects/agasea/; Vaughan *et al.* 2006).

Figure 9*a* 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 9*b*. 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 *q*_{0}, 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 *q*_{0}≈0.04 m^{2} s^{−1}. As above, we assume a Newtonian ice viscosity of *μ*=1×10^{13} 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 9*b*). A different initial condition was determined for each relevant value of *q*_{0}. 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.

Output from instability calculations is shown in figure 10. In that figure, we compare the predicted flux ratio *q*_{G}/*q*_{0} 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.

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 §3*c*, 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 (2007*a*). Incorporation of this profile into a three-dimensional model would mean modification of the basal profile of the central valley, *T*_{S}(*x*) in equation (3.2). Based on the basal topography shown in figure 9*a*, 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 7*b*: 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 (2007*a*) and calculated the response of an idealized ice sheet to changes in sea level. Our results are similar to those of Schoof (2007*a*), 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.

## Acknowledgements

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.

## Footnotes

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

- © 2010 The Royal Society