The stability of grounding lines on retrograde slopes

. The stability of marine ice sheets grounded on beds that slope upwards in the overall direction of ﬂow is investigated numerically in two horizontal dimensions. We give examples of stable grounding lines on such retrograde slopes illustrating that marine ice sheets are not unconditionally unstable in two horizontal dimensions. Retrograde bed slopes at the grounding lines of marine ice sheets, such as the West Antarctic Ice Sheet (WAIS), do not per se imply an instability, nor do they imply that these regions are close to a threshold of instability. We therefore question those estimates of the potential near-future contribution of WAIS to global sea level change based solely on the notion that WAIS, resting on a retrograde slope, must be inherently unstable.


Introduction
Large parts of the West Antarctic Ice Sheet (WAIS) rest below sea level on a bedrock sloping downwards towards the ice sheet's centre. It has been argued that such marine-type ice sheets on retrograde slopes are inherently unstable and susceptible to rapid disintegration (e.g. Weertman, 1974). The question of WAIS being potentially subject to a marine ice-sheet instability (MISI) is of utmost importance for any quantitative estimates of future sea level change, because it raises the possibility of WAIS rapidly contributing several meters of sea level change.
The possibility of a MISI and the conditions under which such an instability could manifest itself have been investigated by a number of authors. To date most analytical and numerical work has focused on analysing the stability regime of MISI in situations where the problem geometry and the flow field varies primarily in along-flow direction. The problem can then be treated as a problem in one horizontal dimension (1HD).
In what can be considered to be the first quantitative study of this problem, Weertman (1974) concluded that a marinebased ice sheet resting on a retrograde slope is inherently unstable. Further work by Thomas and Bentley (1978) reached similar conclusions and gave increased weight to the idea of a MISI and its relevance to WAIS. Hindmarsh (1993Hindmarsh ( , 1996 revisited this problem and found, in contrast to previous work, a marine-type ice sheet resting on a retrograde slope to be neutrally stable and therefore not susceptible to an unstable retreat or advance once perturbed away from a steady state. Schoof (2007a), in a further study of this same problem, however concluded that rapidly sliding marine-based ice sheet do not exhibit neutral equilibrium and that steady grounding lines cannot be stable on retrograde bed slopes.
The difficulty in settling the question of whether marinetype ice sheets are subject to MISI, has been further confounded by the fact that the numerical solution of this problem has turned out to be more difficult than originally envisioned. Vieli and Payne (2005) showed that calculated grounding-line positions, using some commonly used numerical procedures, exhibited a strong dependency on grid size and this question has recently been addressed further in Gladstone et al. (2012). A recent model intercomparison experiment (Pattyn et al., 2012) focusing on a 1HD problem furthermore demonstrated that models based on the shallowice approximation -a class of models frequently used to simulate the dynamics of large ice masses -do not allow for correct description of grounding-line motion. Encouragingly though, the same intercomparison experiment also showed a good overall agreement between numerical models that both (a) account for horizontal transmission of stresses between the grounded and the floating parts, and (b) employ very high spatial resolution (less than one ice thickness) in the vicinity of the grounding line. These models produced in turn steady-state solutions that were in a close agreement with a semi-analytical solution (Schoof, 2007a) of this problem. This combination of theoretical work (Schoof, 2007a) and communal effort by the numerical community (Pattyn et al., 2012) has therefore now provided a general consensus that rapidly sliding 1HD marine-based ice sheets are indeed subject to the MISI. However, the more general case in two horizontal dimensions (2HD) is still unanswered.
The possible stabilising effect of ice-shelf buttressing on a grounding-line position has been investigated in several recent studies. Dupont and Alley (2005) presented an example of a stable grounding line on a retrograde slope. Their model was a 1HD flow-line model using a parameterised side drag and a varying degree of imposed buttressing at the calving front. The resulting buttressing at the grounding line was, thus, in effect prescribed and it is unclear what three-dimensional geometrical configuration, if any, gives rise to buttressing used by Dupont and Alley (2005) in their example. A similar flow-line modelling approach was used by Jamieson et al. (2012) to study transverse flow effects. Jamieson et al. (2012) state that lateral drag can give rise to "transient stabilizations" on reverse bed slopes. Stability is a property of steady states and the meaning of the term "transient stabilizations" as used by Jamieson et al. (2012) is unclear.
To our knowledge the only study showing stable grounding lines on retrograde slopes, while resolving stresses in both horizontal dimensions, is that of Goldberg et al. (2009). Goldberg et al. (2009) performed a series of numerical experiments using a vertically integrated model that accounted for all the stress terms in 2HD otherwise missing in vertically integrated 1HD models. One of their numerical experiments suggested that ice-shelf buttressing might restore stability and they concluded that a further investigation was warranted.
Recently, Katz and Worster (2010) investigated the stability of ice-sheet grounding lines in 2HD numerically using a vertically integrated formulation of the force balance at the grounding line. They concluded that unstable retreat of grounding lines over retrograde beds is a "robust" feature of such models. Katz and Worster (2010) assumed the stress vector at the grounding line to be aligned in the direction normal to the grounding line. Their model study did therefore not fully account for the variation of stresses in 2HD. The fact that Katz and Worster (2010) did not find an example of a stable grounding line on a retrograde slope in 2HD clearly does not exclude the possibility that such an example may exist.
It is, hence, currently not clear if the MISI applies unconditionally in 2HD. On the other hand, to date neither an an-alytical or a numerical example of a stable grounding-line configuration of a marine ice sheet resting on a retrograde slope has been found. The purpose of this work is to provide such an example, and to show that the MISI does not hold unconditionally in 2HD.

Problem definition
The model setup is motivated by the "Marine Ice Sheet Model Intercomparison Project" (MISMIP) experiment 3 (Pattyn et al., 2012), and can be considered to represent a three-dimensional (i.e. two horizontal dimensional) extension of that problem. The MISMIP experiment 3 was in turn inspired by Schoof (2007a) who, in his analysis of grounding-line dynamics in one horizontal dimension, used a particular type of bed profile to demonstrate the marine instability mechanism. This MISMIP experiment 3 has become a benchmark case for assessing and testing the capability of numerical models to simulate the advance and retreat of the grounding line. The experiment also provides a striking example of the marine instability mechanism in 1HD.
The model domain used here stretches from 0 to 1800 km in x direction, and from −120 to 120 km in y direction. For x = 0, both horizontal velocity components are set to zero, i.e. u(0, y) = v(0, y) = 0, and along the left and right hand sides where y = ±120 km, the y component is set to zero, i.e. v(x, ±120 km) = 0. No shear stresses are applied, and the ice is allowed to flow freely, along the left and right side margins, i.e. for y = ±120 km. Along the ice shelf edge, at x = 1800 km, the ocean pressure is applied. The form of this shelf-edge boundary condition is specific to the two models used and is described below.
The bed geometry used here can be written as where and where the units are meters, and B(x, y) stands for the topography of the ocean floor. The bed is in the form of a channel incised into a slowly undulating plane with an overall downward slope in x direction (see Fig. 1). The parameter f c controls the slope of the flanks, d c controls the depth of the channel, and w c is the half-width of the channel. The only parameter that is varied in the model runs described below is the half-width (w c ).
As can be seen from Eqs. (1)-(3), the slope of the bed in x direction is independent of y, i.e. ∂ 2 yx B(x, y) = 0. The slope in x direction is zero at x = 0, x = x a and x = x b where x a = 973.7 km and x b = 1265.7 km, and for x a < x < x b the bed slopes upwards (∂ x B(x, y) > 0) and is retrograde across the entire domain, i.e. ∂ x B(x, y) > 0 for x a < x < x b independently of y. The values for the geometrical parameters B 0 , f c , d c used in all runs are listed in Table 1.

Numerical models
We use two different numerical models: (1) a threedimensional "full Stokes" model, and (2) a vertically integrated model of a type commonly used in glaciology to model the flow of ice shelves and ice streams (in the glaciological literature, this second model type is often referred to as either "SSA", "SSTREAM", or as a "shelfy stream" model). Both numerical models are based on the finite elements method and both are capable of using unstructured grids with spatially variable element sizes. The full Stokes solver uses the open source finite element software package Elmer and the solver will be referred to here as Elmer/Ice. The equations representing the vertically integrated model are solved using a finite element code developed by one of the authors (G. Hilmar Gudmundsson), and the code will here be referred to asÚa. Most of the calculations presented here were obtained withÚa. The Elmer/Ice code was primarily used to demonstrate that the some key aspects of the results obtained withÚa are not specific to that model. The vertically integrated model used here is one of the most commonly used models in glaciology for describing the flow of ice shelves and ice streams. Derivation of the model can be found, for example, in Hutter (1983), Morland (1984), Muszynski and Birchfield (1987), MacAyeal (1989), and Baral and Hutter (2001) and properties of the linearised equations are discussed in some detail in Gudmundsson (2008). The full Stokes model includes all the terms of the Stokes equations, i.e.
where σ pq are the components of the Cauchy stress tensor, and f p the components of the body force, with f = ρg(0, 0, −1) T . The vertically integrated model solves a simplified stress balance that can be written in the form where T = 2τ xx + τ yy τ xy lτ xy 2τ yy + τ xx , with where h, when used as a subscript, is a mnemonic for "horizontal". In the above equation τ ij are the components of the deviatoric stress tensor, s is the surface topography, h is the ice thickness, ρ is the ice density, g is the gravitational acceleration, and t bh is the horizontal part of the bed-tangential basal traction t b where withn being a unit normal vector to the bed pointing into the ice. Both models employ the Weertman sliding law where is the basal sliding velocity, C is the basal slipperiness, and m a stress exponent. Both models also employ Glen's flow laẇ where τ is the second invariant of the deviatoric stress tensor ij are the strain rates, A is the rate factor, and n a stress exponent. In addition to the equations above (and not listed here) both models use appropriate forms of the massconservation equation for an incompressible medium, and employ the kinematic boundary conditions along free upper and lower surfaces.
Along the calving front a pressure boundary condition is applied that corresponds to ocean surface level at z = 0, i.e.
wheren is a unit normal vector pointing horizontally outward from the ice front, and p w is the hydrostatic pressure in the ocean. The surface of the ocean (S) is at S = 0, and the oceanic pressure is p w = ρ w g(S − z) for z < S and zero otherwise, where ρ w is the ocean density. For vertically integrated model the boundary condition at the ice-shelf edge at x = 1800 km takes the form provided the ice is fully floating at the calving front, and where In the vertically integrated model the floating condition takes the form h < h f where where H is the ocean depth. The quantity h f represents the maximum possible ice-shelf thickness at any location. A detailed description of the contact condition that arises in the treatment of the grounding line in the full Stokes model, and its implementation, is given in Favier et al. (2012). The two numerical codes, i.e. Elmer/Ice andÚa, have been used extensively in the past to calculate grounding-line migrations for various different setups (Durand et al., 2009a;Gagliardini et al., 2010;Pattyn et al., 2012). Both models participated in the recently finalised "Marine Ice Sheet Model Intercomparison Project" (MISMIP) (see Pattyn et al., 2012) and results from these two models were submitted to the ongoing "Marine Ice Sheet Model Intercomparison Project for Planview Models" (manuscript submitted to Journal of Glaciology).
Automated remeshing was used to ensure a highly dense discretisation of the area around the grounding line. Calculations of grounding-line migration are known to require a finely resolved mesh around the grounding-line area (e.g. Durand et al., 2009b). How accurate the mesh needs to be is a question that is difficult to answer a priori, but can be answered easily by performing a number of runs using increasingly fine meshes until no resolution dependent behaviour is seen.
Runs performed with solverÚa were done using linear, quadratic, and cubic triangular elements, however the set of results fromÚa shown here were obtained using quadratic elements. The Elmer/Ice solver used linear prismatic elements stabilised using the bubble-stabilising method. The exact sizes and number of elements varied between runs, but, as an example, results shown below usingÚa for a halfwidth of w c = 50 km used a mesh of 228 537 elements with 457 340 nodes, with median, maximum and minimum element sizes of 569 m, 26 862 m, and 86.7 m, respectively. The region around the grounding line, and areas where ice thickness changed most markedly was resolved much more finely than other areas.Úa uses an automated remeshing algorithm where the mesh is refined based on several possible userdefined criteria. The criteria used in the runs shown here were (1) distance from flotation, and (2) second derivative of ice thickness. The first criterion gives small elements in the vicinity of the grounding line, and is similar to the one used by Goldberg et al. (2009).
A fully implicit forward integration procedure was used bý Ua, with the prognostic and the diagnostic equations being solved simultaneously using the Newton-Raphson method. Further details of the solution procedure of the full Stokes model Elmer/Ice can be found in Durand et al. (2009a).

Results
A number of runs using the parameters listed in Table 1 and the bed geometry defined by Eqs. (1)-(3) were performed and the runs continued until steady state was reached. Steady state was considered to have been reached once both the mean rate of surface elevation change over the surface nodes was less than 0.001 m a −1 .
The values used for surface accumulation a, the slipperiness parameter C and the stress exponent m in Weertman's sliding law (Eq. 9), and the stress exponent n in Glen's flow law (Eq. 11) are identical to those used in MISMIP experiment 3a (see Pattyn et al., 2012). The surface accumulation, a, is set at a = 0.3 m a −1 across the whole surface. Note that  Table 1. The two magenta lines at x = 973.7 km and x = 1265.7 km indicate where bedrock slope in x direction is equal to zero. The strip 973.7 km < x < 1265.7 km is an area of retrograde slope where the bed slopes upwards in x direction. The grounding line marking the division between grounded and floating ice is shown as a thick black line. The arrows represent surface velocities, with the corresponding colorbar shown to the right. To make the direction of the shortest velocity vectors visible in the figure, the speed has been set to 100 m a −1 wherever calculated velocities are less than that value. A and n are the rate factor and the stress exponent of Glen's flow law, respectively, C and m are the basal slipperiness and the stress exponent of Weertman's sliding law, and ρ and ρ w are the specific densities of ice and ocean. The variable a is the surface mass balance in the units of ice equivalent. The number of days in a year is 365.25. An example of calculated steady-state grounding-line position is given in Fig. 2. This particular run was obtained with Ua for w c = 50 km. Note that in Fig. 2 only a part of the model domain is shown and that the whole domain stretches in x direction from 0 to 1800 km. A perspective plot of the same model output, showing the whole model domain, is given in Fig. 3.

Parameter Value Units
As Fig. 2 shows the grounding line (shown in black) curves considerably within the xy plane. At the left and righthand side margins of the numerical domain (y = ±120 km), the grounding line is located at approximately x = 1400 km. Where the grounding line crosses the medial line (y = 0), it is located at about x = 1100 km. As a consequence, a sizeable confined ice shelf is formed. The upper and lower limits of retrograde bed slope, x a = 973.7 km and x b = 1265.7 km, are marked by two magenta lines in the Fig. 2. As the figure illustrates, sections of stable grounding line are located on retrograde slopes. A perspective plot showing the steadystate model geometry is given in Fig. 3. Inspection of Figs. 2 and 3 reveals that the highest velocities along the grounding line, and the largest ice thicknesses, are found in the area of retrograde slope.  Table 1. The perturbations were applied for 100 yr, and the model then run towards steady state. As the figure shows both models converged back to the original steady state, showing that the initial steady state is stable with respect to such perturbations in C.
The steady-state geometry shown in Figs. 2 and 3 was arrived at by starting with an initial constant ice thickness of 100 m. In the span of several thousand of years, the numerical solution converged slowly to the steady-state solution shown in the figures, and the possibility of the steady-state being unstable with respect to perturbations in ice thickness can be excluded. To illustrate that the steady-state solution shown in Figs. 2 and 3 is stable with respect to perturbations in other key model parameters as well, the steady state solution was subjected to perturbations in both basal slipperiness (C). Starting from steady state, the slipperiness was increased/decreased by a factor of 2/0.5, respectively. The perturbations were applied for 100 yr and the models then run towards steady state. The resulting grounding-line positions at the medial line are shown in Fig. 4 as functions of time. As the figure shows, the grounding lines reverted with time back to the original position once the perturbations were no longer applied. The grounding line is, hence, stable with respect to perturbations in slipperiness.
Comparison between results obtained by Elmer/Ice andÚa is given in Fig. 5 for w c = 50 km. As can be seen the surface profiles obtained with the two models closely agree, and in both cases the grounding line is located at a section of the bed with a retrograde slope. Another model-intercomparison experiment using w c = 40 km resulted in a similarly close agreement. The results obtained with Elmer/Ice were arrived at by starting with the final geometry as given byÚa and run until a new steady-state was reached.  Runs were performed with half-widths (w c ) ranging from 20 to 70 km. Each run was started by taking an ice-thickness distribution from a previous run for a similarly shaped bedrock as a starting point. Again the possibility of the models runs converging with time to an unstable steady-state can be excluded.
The steady-state grounding-line position (x gl ) at the medial line (y = 0) for several different models with half-widths (w c ) ranging from 20 to 70 km is shown in Fig. 6. These calculations were performed for the vertically integrated model (see Eq. 5) with the numerical codeÚa. The results indicate a gradual change in grounding-line position as a function of channel width. Several cases of stable steady-state grounding lines located on retrograde slopes were found. As expected, the grounding line is no longer located on retrograde slopes for both large and small w c values. Figure 6 suggest that, for the particular parameter set used (see Table 1), grounding lines are located on retrograde slopes for w c within the range of 35 to 55 km.
A close inspection of Fig. 2 reveals that at around (x, y) = (1400, ±50) km, the grounding line (shown in black) extends some distance downstream and appears to break up in a few isolated regions of grounded ice. This is not an artefact of the figure. In this area the ice is either grounded or close to being grounded, and in the course of the transient run towards steady state various smaller patches of ice become either grounded or ungrounded.
A further, in our view an entirely realistic, feature found in a number of our model runs are small (on the order of one ice thickness in horizontal dimension) regions of thin ice that form downstream of the grounding line within the unconfined section of the ice shelf. For numerical purposes, the thickness was forced to be slightly positive. In the runs shown here, the minimum thickness was set at 10 m. This thickness constraint was implemented in a similar manner in both Elmer/Ice andÚa. The method used in Elmer/Ice is described in Durand et al. (2009a). TheÚa solver uses for this purpose an active-set method where any violated thickness constraints are activated using Lagrange multiplier approach. The modified system is then solved again, and the sign of the Lagrange multipliers (or "slack variables" as these variables are sometimes referred to in this context) used to determine if an active constraint should subsequently be inactivated, or if further constraints need to be included in the active set. This procedure is repeated until the active set no longer changes, before progressing to the next time step. The number of active constraints differed between runs and changed in the course of the transient calculations, but the number of active nodal constraints was usually less than 10. For the final steady-state solution shown in Fig. 2, for example, the thickness constraint was activated at 7 nodes in total (the total number of nodes in the finite-element mesh at the end of this run was 457 340). These 7 nodes were located at approximately (x, y) = (1460 ± 88) km. This is a region of the ice shelf that is subjected to strong longitudinal extension and comparatively small concomitant transverse compression. Runs with different minimum ice thickness constraints of 1 m, and 0.1 m, showed that the exact value of the thickness constraint used was, for the purpose of the model exercise, irrelevant, and had no effect of whether the grounding line was located on a retrograde slope or not.

Discussion
The basis for the MISI is the idea of ice flux at the grounding line being a monotonically increasing function of ice thickness (h gl ) at the grounding line. If this holds, then a simple heuristic argument (Weertman, 1974) shows steady-state grounding-line positions of marine-type ice sheets on retrograde slopes to be unstable. As shown by Schoof (2007b), ice flux can be expressed as a function of both h gl and longitudinal stress, with ice flux being an increasing function of both of these variables. It turns out that in 1HD, longitudinal stress is also an increasing function of thickness and, hence, ice flux always an increasing function of h gl in 1HD. It is important to realise that this argument does not hold in 2HD. In 2HD it is conceivable, in principle, to have a configuration where longitudinal stress decreases with thickness with the overall effect of flux being a decreasing function of thickness. In that case, the MISI is not expected to operate and a stable grounding-line position on a retrograde slope becomes a possibility. However, hitherto, no such example of a stable grounding line in 2HD has been provided.
The absence of such a specific counter-example may have been taken to suggest that the MISI applies unconditionally to the 2HD situation. Observations of grounding lines located on reverse bed slopes have been used repeatedly in the past as the sole basis on which to argue that large sections of WAIS are either unstable, or close to a threshold of instability (e.g. Ross et al., 2012). Some recent estimates (e.g. Bamber et al., 2009) of potential sea-level rise from WAIS implicitly assume that MISI operates in 2HD. By providing a numerical example of a marine-type ice sheet in a stable steady-state configuration with the grounding line located on a retrograde slope, we have now shown marine ice sheets to be conditionally stable. Our finding that marine-type ice sheets on retrograde slopes are not unconditionally stable do not contradict, and are fully compatible with, theoretical and numerical work showing that, in 1HD, rapidly sliding marine ice sheets on retrograde slopes are unconditionally unstable.
For our particular parameter set (see Table 1), we find stable grounding-line positions for half-widths between the approximate range of 35 to 55 km. We do not expect this range of half-widths values (w c ) to remain the same for other parameter sets. Neither do we expect the existence of stable grounding-line positions on retrograde slopes to be a simple function of channel widths. The stability regime can be expected to be primarily determined by the level of ice-shelf buttressing at the grounding line, with the buttressing in turn being a (potentially complicated) function of bed geometry and ice rheology. The stability of a given marine ice-sheet configuration can most likely not be assessed from geometrical considerations (bedrock slope, ice-stream width, etc.) only, and each situation may require targeted modelling efforts.
In many situations the sensitivity of ice sheets to external forcing is of key interest. We stress that here we have only addressed the very different question of the stability of marinetype ice sheets. Our results do not, for example, exclude the possibility that the sensitivity of grounding-line position to changes in basal melt rates depends on bed slopes. Note in this respect that in 1HD (and provided the horizontal variation of shear stresses can be ignored), the grounding-line position is independent of any forcing to which an unconfined ice shelf is subjected to. In other words, in 1HD the sensitivity of (stable) grounding-line positions to oceanic forcing is not only independent of bed slope, but identically equal to zero. Although we have not investigated this possibility here, we find it likely that the 2HD situation is, in this respect, also sharply different from the 1HD case.
Numerical modelling has shown ice sheets to respond to changes in ice-shelf melting in complicated and at times nonintuitive ways (Walker et al., 2008;Gagliardini et al., 2010). We do not discard the possibility that further work, using both 2HD and 3-D models, may show grounding lines on retrograde slopes to be stable yet at the same time highly sensitive to changes in external forcing such as ice-shelf melt rates.

Conclusions
We have provided specific numerical examples of model geometries for which sections of stable grounding lines are located on retrograde slopes. We conclude that marine ice sheets on retrograde slopes are not unconditionally unstable.
Our confidence in our model results is strengthened by the fact that we have solved the flow problem using two different numerical codes (Úa and Elmer/Ice), based on two different set of model assumptions.
The conditional stability in 2HD does not contradict previous findings of unconditional instability in 1HD. These differences simply show that the 1HD and the 2HD problems differ fundamentally in this respect. Real ice sheets are of course three dimensional features and we caution against applying arguments to ice sheets that are inherently tied to 1HD. These stark differences between the 1HD and the 2HD cases should also be considered in any applications of 1HD flow-line models, and the use of 1HD flow-line models to investigate the stability regime of marine ice sheets is problematic.
Observations of retrograde bed slopes at the grounding lines of marine ice sheets, such as the West Antarctic Ice Sheet, do not as such imply an unstable mode of retreat or advance, nor do they imply that such regions are close a threshold of instability. Questions regarding stability of a given marine ice-sheet configuration are unlikely to be answered without targeted modelling efforts.