# Full Stokes modeling of marine ice sheets: influence of the grid size

Gaël Durand, Olivier Gagliardini, Thomas Zwinger, Emmanuel Le Meur, Richard C.A. Hindmarsh

► To cite this version:

Gaël Durand, Olivier Gagliardini, Thomas Zwinger, Emmanuel Le Meur, Richard C.A. Hindmarsh. Full Stokes modeling of marine ice sheets: influence of the grid size. Annals of Glaciology, 2009, 50 (52), pp.109-1144. 10.3189/172756409789624283. insu-00420761

HAL Id: insu-00420761

https://insu.hal.science/insu-00420761

Submitted on 7 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Full Stokes modeling of marine ice sheets: influence of the grid size

Gaël DURAND,<fnref n="1" /> Olivier GAGLIARDINI,<fnref n="1" /> Thomas ZWINGER,<fnref n="2" /> Emmanuel Le MEUR,<fnref n="1" /> Richard C.A. HINDMARSH<fnref n="3" />

**ABSTRACT.** Using the finite-element code Elmer, we show that the full Stokes modeling of the ice-sheet/ice-shelf transition we propose can give consistent predictions of grounding-line migration. Like other marine ice-sheet models our approach is highly sensitive to the chosen mesh resolution. However, with a grid size down to <5 km in the vicinity of the grounding line, predictions start to be robust because: (1) whatever the grid size (<5 km) the steady-state grounding-line position is sensibly the same (6 km standard deviation), and (2) with a grid-size refinement in the vicinity of the grounding line (200 m), the steady-state solution is independent of the applied perturbation in fluidity, provided this perturbation remains monotonic.

## INTRODUCTION

Marine ice sheets rest on beds that lie below sea level and their drainage takes place through surrounding ice shelves. Grounded ice-sheet flow is dominated by horizontal shearing while the ice-shelf flow is dominated by longitudinal stretching and lateral shearing. The two types of flow couple together across a transition zone near the grounding line, where longitudinal and shear stresses are of the same order of magnitude. A long debate on the dynamics of such ice sheets was initiated in the 1970s, when Weertman (1974) proposed that a marine ice sheet which lies on an upward-sloping bed is unstable. Recently, the instability hypothesis has been strongly reinforced, based on a boundary-layer theory due to Schoof (2007b). Moreover, Vieli and Payne (2005) showed the poor ability of marine ice-sheet models to give consistent prognostic results and, more particularly, they highlighted the influence of the grid size on model results. One of their main conclusions was that no reliable model was able to predict grounding-line dynamics at the time of their study.

As a consequence, there is an urgent need to improve marine ice-sheet models in order to (1) corroborate recent theoretical predictions, and (2) obtain confident simulations of the grounding-line dynamics. We recently proposed a full Stokes resolution of the ice-sheet/ice-shelf transition (Durand and others, in press). This approach has been built on abundant literature dealing with the coupling between a grounded ice sheet and a floating ice shelf and identifying this transition zone as a crucial control of the marine ice-sheet dynamics (e.g. Weertman, 1974; Van der Veen, 1985; Chugunov and Wilchinsky, 1996; Hindmarsh, 1996; Vieli and Payne, 2005; Schoof, 2007a,b). The main equations developed by Durand and others (in press) are recalled here. We then emphasize the impact of the mesh resolution on the calculated steady state with the proposed method. Finally, we show that consistent results can only be obtained if a small grid mesh is prescribed in the vicinity of the grounding line.

## EQUATIONS

### Notations and main hypothesis

The problem to be solved consists of the gravity-driven flow of isothermal, incompressible and non-linear viscous ice. The geometry is restricted to a two-dimensional plane flow perpendicular to the y direction. Ice flows along the x direction, and the z axis is the vertical upward-pointing axis. Notation is detailed in Figure 1. The ice sheet flows over a rigid bedrock, $z = b(x)$, assuming a non-linear friction law for the grounded part, and extends further as an ice shelf over the ocean. In the case of the ice shelf, the ice slides perfectly over the sea. The last grounded point defines the grounding line and is denoted $x_G$ hereafter. The left-hand side of the domain is assumed to be symmetric and the shelf ends at the right-hand side of the domain.

The constitutive law for the ice behavior is a Norton-Hoff type law (Glen's flow law in glaciology):

$$S_{ij} = 2\eta D_{ij}, \quad (1)$$

where $S$ is the deviatoric stress tensor, $D_{ij} = (u_{i,j} + u_{j,i})/2$ are the components of the strain-rate tensor and $u$ is the velocity. The effective viscosity, $\eta$, can be expressed as

$$\eta = B^{-1/n} \gamma_e^{(1-n)/n}, \quad (2)$$

where the strain-rate invariant, $\gamma_e$, is defined as

$$\gamma_e^2 = 2D_{ij}D_{ij}. \quad (3)$$

In the following applications, $n$ is set to 3 and ice is assumed isothermal, such that the fluidity parameter, $B$, remains constant (see values in Table 1).

### State equations and boundary conditions

The ice flow is computed by solving the Stokes problem, expressed by the mass-conservation equation in the case of incompressibility

$$\mathrm{tr} \, D = \mathrm{div} \, u = 0, \quad (4)$$

<footnote n="1">¹Laboratoire de Glaciologie et Géophysique de l'Environnement du CNRS (associé à l'Université Joseph Fournier-Grenoble I), 54 rue Molière, BP 96, 38402 Saint-Martin-d'Hères Cedex, France</footnote>

<footnote>E-mail: durand@lgge OBS.UJF-Grenoble.fr</footnote>

<footnote n="2">²CSC - Scientific Computing Ltd, PO Box 405, FIN-02101 Espoo, Finland</footnote>

<footnote n="3">³British Antarctic Survey, Natural Environment Research Council, Madingley Road, Cambridge CB3 0ET, UK</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Notation for the problem to be solved.</figcaption></figure>

and the momentum conservation equations

$$ \mathrm{div} \, \boldsymbol{\sigma} = \rho_i \boldsymbol{g}, \qquad (5) $$

where $\boldsymbol{\sigma} = \boldsymbol{S} - p\boldsymbol{I}$ is the Cauchy stress with $p$ the isotropic pressure; $\rho_i$ is the ice density and $\boldsymbol{g}$ the gravity vector.

In the present problem, the lateral boundaries of the domain are a dome and a calving front. The dome is assumed to be axisymmetric for the flow problem, which implies that $u_x(0, z) = 0$. The calving front is an artificial cutting of the shelf, which can be seen as the point where icebergs are calved. The exact position of the calving front is not relevant to our problem because it does not influence the upstream flow. This surface undergoes the sea-ice pressure, $p_w(z, t)$, which evolves vertically as

$$ \begin{cases} p_w(z, t) = \rho_w g (l_w(t) - z), & z < l_w(t) \\ p_w(z, t) = 0, & z \ge l_w(t) \end{cases} \qquad (6) $$

where $p_w$ is the sea-water density and $l_w$ the sea level.

In addition to the lateral boundaries, the ice body is bounded by two free surfaces, namely the stress-free upper surface, $z = z_s(x, t)$, and the bottom surface, $z = z_b(x, t)$, at the interface between the bed or sea and the ice. The evolution of the two free surfaces is determined by solving a local transport equation. Note also that the length of the sea/ice interface, starting from the grounding line, $x_G(t)$, is not known in advance and is part of the solution.

The upper surface is a stress-free surface, which implies that $\boldsymbol{n} \cdot (\boldsymbol{\sigma} \cdot \boldsymbol{n})|_s = p_{\mathrm{atm}} \approx 0$, where $\boldsymbol{n}$ is the unit normal vector of the surface pointing outward and the subscript 's' denotes the value taken at the ice/air interface. The equation governing the upper stress-free surface evolution reads:

$$ \frac{\partial z_s}{\partial t} + u_x(x, z_s) \frac{\partial z_s}{\partial x} - u_z(x, z_s) = a(x, t), \qquad (7) $$

where $u_i(x, z_s)$ denotes the upper surface velocity in the horizontal ($i = x$) and vertical ($i = z$) directions and $a(x, t)$ is the accumulation/ablation function given as a vertical flux at the upper surface. In what follows, the accumulation is supposed constant both in space and in time (see Table 1).

The bottom sea stress-free surface obeys the following equation:

$$ \frac{\partial z_b}{\partial t} + u_x(x, z_b) \frac{\partial z_b}{\partial x} - u_z(x, z_b) = 0, \qquad (8) $$

<table><caption>**Table 1.** Values of the parameters used in this study, which correspond to steps 5 and 6 of the Marine Ice-Sheet Intercomparison Project (MISMIP) benchmark, but are expressed differently, as here the fluidity parameter is $B = 2A$. However, numerically the constitutive relations are rigorously the same</caption><thead><tr><th>Parameter</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>a</td><td>0.3</td><td>m a<sup>-1</sup></td></tr><tr><td>B<sub>1</sub></td><td>3.16 &times; 10<sup>-18</sup></td><td>Pa<sup>-3</sup> a<sup>-1</sup></td></tr><tr><td>B<sub>2</sub></td><td>6.31 &times; 10<sup>-18</sup></td><td>Pa<sup>-3</sup> a<sup>-1</sup></td></tr><tr><td>n</td><td>3</td><td></td></tr><tr><td>C</td><td>2.41 &times; 10<sup>4</sup></td><td>Pa m<sup>-1/3</sup> a<sup>1/3</sup></td></tr><tr><td>m</td><td>1/3</td><td></td></tr><tr><td>&rho;<sub>w</sub></td><td>1000</td><td>kg m<sup>-3</sup></td></tr><tr><td>&rho;<sub>i</sub></td><td>900</td><td>kg m<sup>-3</sup></td></tr><tr><td>g</td><td>9.8</td><td>m s<sup>-2</sup></td></tr></tbody></table>

where accretion of sea water by refreezing or melt of bottom ice is neglected. Note that this equation is still valid for the points on the bottom surface which are in contact with the bedrock. Assuming a rigid, impenetrable bedrock, $z = b(x)$, the following topological conditions must be fulfilled by $z_s$ and $z_b$:

$$ z_s(x, t) > z_b(x, t) \geq b(x) \quad \forall x, t. \qquad (9) $$

As a consequence, the unilateral link between the ice and the bedrock can be treated as a contact problem: the ice cannot penetrate the bedrock but is allowed to move away from it (Lestringant, 1994). Resolutions of the contact problem have been inspired by previous studies on subglacial cavities, namely the works of Schoof (2005) and Gagliardini and others (2007). At a given point, $x$, ice is assumed to be in 'true' contact with the bedrock (and corresponding boundary conditions are applied; see below) if the ice touches the bed and the stress exerted by the ice is larger than the sea-water pressure. Conversely, the ice is assumed to be in contact with the sea if the bottom surface is above the bed or if the ice touches the bed and the sea-water pressure remains larger than the normal stress, $\sigma_{nn}$. In other words:

(1) the ice/bedrock boundary condition applies if

$$ z_b(x, t) = b(x) \quad \text{and} \quad -\sigma_{nn}|_b > p_w(z_b, t), \qquad (10) $$

(2) the ice/sea boundary condition applies if

$$ \begin{cases} z_b(x, t) > b(x), \\ \text{or} \\ z_b(x, t) = b(x) \quad \text{and} \quad -\sigma_{nn}|_b \le p_w(z_b, t), \end{cases} \qquad (11) $$

where the subscript $|_b$ denotes the value taken at the bottom surface. For the flow problem, two different boundary conditions have to be applied, depending on whether the ice is in contact with the bedrock or with the sea. For the ice/bedrock boundary condition (i.e. condition (1) above), a non-linear friction law is applied:

$$ \left\{ \begin{array}{l} \tau_b = t \cdot (\boldsymbol{\sigma} \cdot \boldsymbol{n})|_b = Cu_b^m, \\ \boldsymbol{u} \cdot \boldsymbol{n} = 0. \end{array} \right. \qquad (12) $$

where $\tau_b$ is the basal shear stress, $u_b = \boldsymbol{u} \cdot \boldsymbol{t}$ is the sliding velocity at the base and $\boldsymbol{t}$ is the tangent vector to the surface $z_b$. The parameters $C$ and $m$ entering the friction law are given in Table 1. For the ice/sea boundary condition (i.e. condition (2) above), the ice is sliding perfectly over the sea,
<figure><img src="image_3.png" /><figcaption>Fig. 2. Steady-state surface profiles for different grid resolutions. Filled areas correspond to simulations with a constant horizontal grid size over the whole domain: 15, 10, 7.5 and 2.5 km (from black to light gray). Black thin curve depicts results obtained using the adaptive mesh refinement method with Δx₀ = 2.5 km (light gray). The inset shows an enlargement near the grounding line of the upper surface. Black curve with symbols corresponds to the local steady surface profile obtained using the adaptive mesh refinement method with Δx₀ = 2.5 km; the gray curve with symbols to that with Δx₀ = 200 m.</figcaption></figure>

i.e. $t \cdot (\sigma \cdot n)|_b = 0$, and the normal stress is equal to the buoyancy sea-water pressure (Equation (6)).

The equations presented above have been implemented in the finite-element code Elmer (http://www.csc.fi/elmer/). For gravity-driven ice flow we solve the set of the Stokes equations neglecting inertia terms. More details on the numerics are given by Durand and others (in press). Note that a similar approach, i.e. full Stokes finite-element modeling of marine ice sheets, was initiated by Lestringant (1994) and used more recently by Nowicki and Wingham (2008).

## SET-UP OF THE SIMULATIONS

Numerical experiments were carried out using the over-deepening bed presented by Schoof (2007a) which, expressed in our reference frame, is given by:

$$b(x) = 720 - 2184.8 \left( \frac{x}{750 \times 10^3} \right)^2 + 1031.72 \left( \frac{x}{750 \times 10^3} \right)^4 - 151.72 \left( \frac{x}{750 \times 10^3} \right)^6, \quad (13)$$

where x and b are in meters. The values of the different parameters used in this study are detailed in Table 1. Note that the set-up of the experiments corresponds to steps 5 and 6 of the Marine Ice-Sheet Intercomparison Project (MISMIP, http://homepages.ulb.ac.be/fpattyn/mismip/). We used a constant accumulation, $a = 0.3 \text{ m a}^{-1}$, over the whole domain, whereas accretion–melting is neglected below the ice shelf. Initial simulations start with an initial 10 m layer of ice resting on land up to the position where this layer begins to float ($x_G = 482.4 \text{ km}$ for the prescribed $\rho_i$ and $\rho_w$; see Table 1). Downstream of $x_G$, a 10 m thick ice shelf extends the grounded part. The vertical extension of the domain is discretized using $N_{by} = 30$ equal thickness layers, whereas the number of elements in the horizontal direction $N_{bx}$ is adjusted for the different simulations and given below, for

<figure><img src="image_4.png" /><figcaption>Fig. 3. Evolution of $x_G$ through time for runs with various constant horizontal grid size over the whole domain: 20 km (thick black curve), 15 km (thick gray curve), 10 km (thick gray dashed curve), 7.5 km (thin gray dashed curve) and 2.5 km (thick light gray curve). The fine black line indicates results using the adaptive mesh refinement method ($\Delta x_0 = 2.5 \text{ km}$). Horizontal dotted line corresponds to the result obtained with the boundary layer theory developed by Schoof (2007a) for the same prescribed fluidity, $B_1$. This latter is insensitive to grid resolution.</figcaption></figure>

each case. Transient simulations are run until steady state is reached (once the relative variation of volume is $<1 \times 10^{-6}$ between two successive time-steps).

## SENSITIVITY TO THE HORIZONTAL MESH RESOLUTION

Using the settings described in the previous section, different simulations for fluidity, $B_1$ (see Table 1), were driven for a regular horizontal mesh with horizontal grid sizes of 20, 15, 10, 7.5, 5 and 2.5 km. Steady state was reached (for durations of 26.9, 30.0, 24.2, 21.2, 20.1 and 21.2 ka, respectively), and the corresponding ice-sheet and ice-shelf surface profiles are plotted in Figure 2 (neglecting grid sizes 20 and 5 km for the sake of clarity). Evolution of $x_G$ through time is shown in Figure 3.

As shown by Vieli and Payne (2005), marine ice-sheet models are highly sensitive to their horizontal grid size. The results of the present experiments clearly demonstrate that the method described above does not derogate from the Vieli and Payne conclusions, in the sense that important sensitivity to the grid resolution is clearly highlighted as $x_G$ decreases significantly with increasing horizontal resolution (see Figs 2
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Steady-state grounding-line position as a function of the horizontal extension of the grid, $\Delta x_0$, used for the simulations. Filled circles indicate results from the constant-grid mesh simulations, and empty circles those from the adaptive mesh refinement method described in the text. Crosses depict results obtained with a vertically integrated model (Hindmarsh, unpublished data), and the horizontal dotted line represents the result obtained with the semi-analytical formula due to Schoof (2007a).</figcaption></figure>

and 3). Predictions of the model can be dramatically affected, particularly in the presence of an overdeepening bed, as the grounding line is not stable on reverse bed slope (Schoof, 2007a; Durand and others, in press). Indeed, in a reverse-slope region, surrounded by two areas of normal slopes that are two branches of stable solutions, for large grid size, the branch corresponding to the larger solution is sometimes selected (probably erroneously). Such a situation is illustrated in Figure 2, as the larger grid sizes (>15 km) lead to steady positions of the grounding line beyond the unstable area. In practice, typical grid resolutions used for ice-sheet modeling (>10 km) are unable to predict a steady position for the grounding line. It should be noted that computing times increase from 2 to 20 days for 20 and 2.5 km grid resolutions, respectively, using a 2.6 GHz processor. This severely constrains further investigations of the model behavior with high-resolution grid sizes.

## ADAPTIVE MESH REFINEMENT AROUND $x_G$

A crucial validation of the present model would be to find a maximal horizontal grid size below which grounding-line dynamics are no longer mesh-dependent. As mentioned above, accessing fine grid resolutions (<1 km) would require large computational facilities and rapidly becomes unmanageable. Therefore, we developed a method that allows an adaptive refinement of the mesh around the grounding line. A constant grid size, $\Delta x_0$, is set over a distance, $L_f$, centered around $x_G$. Upstream of $x_G - L_f/2$ and downstream of $x_G + L_f/2$, a geometric progression of the horizontal extension of the mesh is prescribed. The total number of elements in the horizontal direction is adjusted to obtain reasonable increases (i.e. <10% increase from one element to the next). Readjustment of the mesh is repeated after each displacement of the grounding line. Note that with this method, the total number of nodes is constant; only the horizontal distribution of the nodes is modified. Sensitivity tests have shown very similar behavior of the grounding-line dynamics

whatever the value of $L_f$. Despite the limited impact of $L_f$, the horizontal mesh extension around $x_G$ was kept constant, as this allows regular displacement of the grounding line during the simulation (i.e. a multiple of $\Delta x_0$ at each time-step).

In Figure 2, the steady surface profile for a simulation with $\Delta x_0 = 2.5$ km over a $L_f = 10$ km range around $x_G$ is shown. A good match between this last surface profile and that obtained with a 2.5 km regular mesh is observed. Relative differences between the two simulations in terms of steady-state grounding-line position and volume are <2%. Note also that the two simulations show very similar transient behavior, as illustrated by the evolutions of $x_G$ versus time plotted in Figure 3. Some conclusions can already be drawn from the comparison of these two simulations. First, this fully validates the adaptive mesh refinement method used here, which then allows for exploration of the model behavior for further refinements of the grid around the grounding line. Second, this shows that the state of a marine ice sheet is strongly influenced by the behavior of a narrow column centered over $x_G$, confirming the theoretical assertion of Schoof (2007b).

Simulations were performed for $\Delta x_0$ ranging from 100 to 2500 m. For each run, the steady-state position of the grounding line is plotted versus $\Delta x_0$ in Figure 4. Results from simulations with a regular mesh size are also shown. As observed in Figure 2, horizontal extension of the ice sheet significantly decreases with increasing grid resolution. However, for $\Delta x_0 \le 5$ km there is not such a clear trend. Unfortunately, because the computing resources required for simulations with $\Delta x_0 \le 100$ m are huge (>2 weeks), we were unable to establish whether the model converges toward a unique value of $x_G$ when $\Delta x_0$ tends to 0. If such convergence exists, it is obviously not monotonic. Despite this limitation, the model seems to give consistent results for resolution finer than 5 km as in this case, the standard deviation of $x_G$ at steady state being 5.5 km. A similar exercise with grid resolutions ranging from 2.5 to 20 km was performed using a vertically integrated stress approximation model with a differentiated margin condition to calculate the margin velocity (Hindmarsh, unpublished information). These results are also plotted in Figure 4. Sensitivity to grid resolution appears to be less important than with the full Stokes finite-element model; for example, although $x_G$ does not converge to a particular value with a smaller grid, it never crosses the upward-sloping area for large grid sizes in that case. Variability from one model to the other is large; in particular, the full Stokes finite-element approach presents a difference of ~100 km in the grounding-line steady position with the boundary layer results. Note that this discrepancy decreases to ~40 km when $x_G$ at steady state moves away from the unstable region (Durand and others, in press).

Further validation of the model requires us to test whether, for a given grid size, the steady state is influenced by the way the ice sheet has been built. Following Schoof (2007a), we chose to vary the fluidity rather than other parameters (e.g. sea level). In what follows, we chose to use $\Delta x_0 = 200$ m. This is a compromise between reasonable calculation times and fine enough resolution to properly describe surface oscillations induced by changes in basal conditions in the vicinity of the grounding line (see inset of Fig. 2 and Durand and others, in press, for more details). A first steady state is obtained starting from a 10 m thick ice layer with a prescribed fluidity of $B_2$ (see Table 1). Starting from the obtained steady state with fluidity $B_2$ at $t = 0$, three simulations were run
with different variations of the fluidity from $B_2$ to $B_1$ (1) a steady decrease from $t=0$ to $t=30$ ka followed by 5 ka of relaxation with $B_1$; (2) ten regular step changes every 3 ka followed by 5 ka of relaxation; and (3) fluidity set to $B_1$ at $t=0$ and kept constant during the 35 ka of the simulation. Evolution of the fluidity with time is plotted in Figure 5a for each run. The corresponding evolutions of $x_G$ with time for the three simulations are plotted in Figure 5b. A notable feature is the ability of the model to react after even a small perturbation. Indeed, following values given by Paterson (1994), each small step of the fluidity in run (2) corresponds to a temperature shift of ~0.5°C. Such sensitivity to small perturbations has also been observed for small modifications in sea level, basal sliding or accumulation. For example, a perturbation of the steady state obtained with $B_2$ with a 1 m sea-level drop leads to a 1.8 km advance of $x_G$ before a new steady state is found (data not shown). This is a very important result as some marine ice-sheet models have shown very poor sensitivity to much more dramatic perturbations (e.g. a few kilometers of grounding-line retreat after a sea-level rise of 125 m; see Vieli and Payne, 2005). Also note that all the simulations have reached their steady state after the 35 ka of the run. Notably, the positions of the grounding line after 35 ka are exactly the same: $x_G = 822.8$ km. This is also the steady position reached if the simulation starts from a 10 m thick ice layer with $B_1$.

The sensitivity of the model to perturbations as well as the reproducibility of the grounding-line positions, whatever the rate of decrease of fluidity, gives confidence in the model results.

## CONCLUSIONS

We recently proposed a new approach to solve the contact problem at the interface between the grounded part of an ice mass and its extension over the sea (Durand and others, in press). In the proposed approach, the ice flow is determined by solving the full Stokes problem implemented in the finite-element code Elmer. Here we have shown that the proposed method, despite its high sensitivity to the horizontal mesh resolution, can give consistent results as long as the horizontal extension of the element is small enough (<5 km) in the vicinity of the grounding line. Indeed, standard deviation of the grounding-line positions at the steady state for simulations with $\Delta x_0 < 5$ km is ~6 km. This value is small in comparison with the previous similar experiments of Vieli and Payne (2005) (40 km is a classical grid size for models using the shallow-ice approximation; Ritz and others, 2001). However, this shows that full Stokes finite-element modeling of marine ice sheets should be handled with care, and that intensive grid sensitivity tests should be conducted. Furthermore, steady-state solutions are independent of the applied perturbation in fluidity, provided this perturbation remains monotonic. It is also worth noting that the model is sensitive to small variations of fluidity which correspond to a perturbation of <1°C.

## ACKNOWLEDGEMENTS

This research is a part of the DACOTA project (ANR-06-VULN-016-01) funded by the Agence National de la Recherche, France. Most of the computations presented in this paper were performed at the Service Commun de Calcul Intensif de l'Observatoire de Grenoble.

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** (a) Evolution of the fluidity versus time for run (1) in light gray, run (2) in dark gray and run (3) in black. All simulations start from the steady state obtained with fluidity $B_2$ (see text and Table 1 for details). (b) Evolution of the grounding line versus time for the different runs. For a given run, gray tones are the same in (a) and (b).</figcaption></figure>

## REFERENCES

* Chugunov, V.A. and A.V. Wilchinsky. 1996. Modelling of a marine glacier and ice-sheet-ice-shelf transition zone based on asymptotic analysis. *Ann. Glaciol.*, **23**, 59–67.

* Durand, G., O. Gagliardini, B. de Fleurian, T. Zwinger and E. Le Meur. In press. Marine ice-sheet dynamics: hysteresis and neutral equilibrium. *J. Geophys. Res.* (10.1029/2008JF001170.)

* Gagliardini, O., D. Cohen, P. Råback and T. Zwinger. 2007. Finite-element modeling of subglacial cavities and related friction law. *J. Geophys. Res.*, **112**(F2), F02027. (10.1029/2006JF000576.)

* Hindmarsh, R.C.A. 1996. Stability of ice rises and uncoupled marine ice sheets. *Ann. Glaciol.*, **23**, 105–115.

* Lestringant, R. 1994. A two-dimensional finite-element study of flow in the transition zone between an ice sheet and an ice shelf. *Ann. Glaciol.*, **20**, 67–72.

* Nowicki, S.M.J. and D.J. Wingham. 2008. Conditions for a steady ice sheet-ice shelf junction. *Earth Planet. Sci. Lett.*, **265**(1–2), 246–255.

* Paterson, W.S.B. 1994. *The physics of glaciers*. Third edition. Oxford, etc., Elsevier.
Ritz, C., V. Rommelaere and C. Dumas. 2001. Modeling the evolution of Antarctic ice sheet over the last 420 000 years: implications for altitude changes in the Vostok region. *J. Geophys. Res.*, **106**(D23), 31,943–31,964.

Schoof, C. 2005. The effect of cavitation on glacier sliding. *Proc. R. Soc. London, Ser. A*, **461**(2055), 609–627.

Schoof, C. 2007a. Ice sheet grounding line dynamics: steady states, stability, and hysteresis. *J. Geophys. Res.*, **112**(F3), F03S28. (10.1029/2006JF000664.)

Schoof, C. 2007b. Marine ice-sheet dynamics. Part 1. The case of rapid sliding. *J. Fluid Mech.*, **573**, 27–55.

Van der Veen, C.J. 1985. Response of a marine ice sheet to changes at the grounding line. *Quat. Res.*, **24**(3), 257–267.

Vieli, A. and A.J. Payne. 2005. Assessing the ability of numerical ice sheet models to simulate grounding line migration. *J. Geophys. Res.*, **110**(F1), F01003. (10.1029/2004JF000202.)

Weertman, J. 1974. Stability of the junction of an ice sheet and an ice shelf. *J. Glaciol.*, **13**(67), 3–11.