Edinburgh Research Explorer Results of the Marine Ice Sheet Model Intercomparison Project, MISMIP

. Predictions of marine ice-sheet behaviour require models that are able to robustly simulate grounding line migration. We present results of an intercomparison exercise for marine ice-sheet models. Veriﬁcation is effected by comparison with approximate analytical solutions for ﬂux across the grounding line using simpliﬁed geometrical conﬁgura-tions (no lateral variations, no effects of lateral buttressing). Unique steady state grounding line positions exist for ice sheets on a downward sloping bed, while hysteresis occurs across an overdeepened bed, and stable steady state grounding line positions only occur on the downward-sloping sections. Models based on the shallow ice approximation, which does not resolve extensional stresses, do not reproduce the approximate analytical results unless appropriate parameterizations for ice ﬂux are imposed at the grounding line. For extensional-stress resolving “shelfy stream” models, differences between model results were mainly due to the choice of spatial discretization. Moving grid methods were found to be the most accurate at capturing grounding line evolution, since they track the grounding line explicitly. Adaptive mesh reﬁnement can further improve accuracy, including ﬁxed grid models that generally perform poorly at coarse resolution. Fixed grid models, with nested grid representations of the grounding line, are able to generate accurate steady state positions, but can be inaccurate over transients. Only one full-Stokes model was included in the intercomparison, and consequently the accuracy of shelfy stream models as approximations of full-Stokes models remains to be determined in detail,


Introduction
Recent years have seen significant effort invested in developing numerical marine ice sheet models that are intended to have the capability of accurately representing the motion of an ice-sheet grounding line. This effort has seen a variety of models emerge. However, it is unclear to what extent these models agree with one another and how well they are able to model real marine ice sheets. The purpose of the MISMIP Marine Ice Sheet Model Intercomparison Project is to address the first issue, examining how well marine ice sheet models agree with one another. A second purpose of the experiments is to assess how well the numerical schemes used solve the partial differential equations on which they are based. Not all the models used in marine ice sheet simulations employ the same governing equations, but many of them share basic characteristics.
Theoretical arguments, due to Weertman (1974) and Thomas and Bentley (1978), indicate marine ice sheets have one or more distinct equilibrium shapes and that only those with grounding lines located on downward slopes could be stable to small perturbations. That is, only if the grounding line is located on a downward slope should a slightly perturbed ice sheet return to its original steady state over time. By contrast, later work suggested that they could exhibit "neutral equilibrium" (Hindmarsh, 1993), that is, a perturbation in grounding line position away from steady state could result in a distinct, new steady state close to the original one. This issue has continued to dominate studies of marine ice sheets and numerical simulations in particular. The impetus for the intercomparison has been providedamongst others -by the papers of Vieli and Payne (2005), Pattyn et al. (2006), Hindmarsh (2006) and Schoof (2007aSchoof ( ,b, 2011. These have demonstrated not only the strong possibility of numerical artifacts in marine ice sheet simulations, but also the importance of grid resolution and of accurate representation of the sheet-shelf transition zone. The papers above have also shed some light on the question of steady states and their stability. No stable steady states were found on upward-sloping ice sheet beds, in line with the earlier theory of Weertman (1974) and Thomas and Bentley (1978), but some of the numerical results in Vieli and Payne (2005) and Pattyn et al. (2006) have left open the possibility of "neutral equilibrium". Using an analytical approach based on asymptotic expansions, Schoof (2007a) concluded that these results were likely to be numerical artifacts for parameter regimes in which there was significant sliding, predicting instead the existence of discrete steady states dependent only on accumulation rates, bedrock geometry and parameters describing ice deformation and basal sliding. However, in the absence of basal sliding, Nowicki and Wingham (2008) raise the possibility of non-unique steady states.
In addition to the question of whether steady states are distinct and whether their stability is controlled by bed slope, other important issues to address are (i) how do equilibria de-pend on bed geometry and the physics of sliding at the bed, ice viscosity and gravitational forces as well as accumulation rates, (ii) is hysteresis under changes in forcing and internal physical properties possible when the bed is overdeepened, (iii) to what extent is high grid resolution, especially near the grounding line, necessary to obtain reliable results, and (iv) is high resolution particularly important when modelling transients?
Here we use the development of a boundary layer theory for sheet-shelf interactions with rapid sliding (Schoof, 2007a) to facilitate an intercomparison exercise. This theory takes a complementary approach to numerical marine ice sheet models: it uses a systematic set of approximations to "parameterize" the sheet-shelf transition zone in a simpler "shallow ice" model, yielding a relationship between flux across the grounding line and ice thickness at the same location. The difference with other models is that the theory is not specific to a particular numerical method. Importantly, it allows steady states to be computed semi-analytically and at low computational cost, and this has been used to guide the experimental design for the present intercomparison.
The boundary layer refers to a zone a few kilometres to tens of kilometres upstream of the grounding zone (Hindmarsh, 2006;Schoof, 2007a) where the flow changes from shelf-like at the grounding line (with significant extensional stresses) to sheet-like further upstream (with flow dominated by shear stresses). It is expected that producing an accurate solution within the boundary layer is crucial to computing grounding line motion accurately in numerical simulations (Schoof, 2007a;Goldberg et al., 2009), and this idea lies behind much of the experimental design of this intercomparison.
This paper describes the results of the first Marine Ice Sheet Model Intercomparison Project (MISMIP), based on an idealised two-dimensional ice sheet geometry. In total, 27 flow-line models participated in the project. This level of participation is comparable to the ISMIP-HOM intercomparison project for higher-order ice sheet models (Pattyn et al., 2008). Moreover, there is a sufficient spread in numerical approaches that allow for a broadly-based intercomparison. The first section describes the basis of marine ice sheet models and their boundary conditions, followed by the description of the experiments. In the next section, results are presented and analysed for both steady state conditions and transient states, and the final section is the discussion of these results. to strain rates, i.e., where τ is the deviatoric stress tensor and D ij are the components of the strain rate tensor, defined by and where u is the velocity vector. The effective viscosity η is expressed as where the strain-rate invariant, D e is defined as 2D 2 e = D ij D ij (using the usual summation convention). We use a spatially uniform coefficient A, whose value is varied between experiments, and set n = 3. The flow of an ice body is computed by solving the Stokes problem (Stokes, 1845), where p is the isotropic pressure and g the gravitational acceleration.

Boundary conditions for the underlying Stokes flow model
The boundary conditions we use are essentially the same as in Durand et al. (2009b). We consider an ice sheet along a flow-line (x-direction) in plane flow, without lateral variations. The ice is, therefore, delimited in the vertical by two free surfaces, i.e., the top interface z = z s (x, t) between ice and air, and the bottom interface z = z b (x, t) between ice and bedrock or sea. We denote bedrock above sea-level b(x), which is assumed to be fixed. At x = 0, the ice divide is a symmetry axis, so that the horizontal velocity component is u x (0, z) = 0. The other end of the domain is a calving front kept at a fixed position x = x c . We assume that the ice ends in a vertical cliff there. The front boundary is subject to a normal stress due to the sea pressure p w (z, t) that depends on elevation as: where ρ w is the sea water density and z w is sea level. Let n be the outward-pointing unit normal to the ice surface and t a unit tangent vector. Denoting total normal and shear stress by σ nn = τ ij n i n j − p, σ nt = τ ij n i t j we can write boundary conditions at the shelf front as σ nn = p w and σ nt = 0. The upper surface z = z s (x, t) is stress-free, implying that The lower interface z = z b (x, t) is either in contact with the ocean or the bedrock, resulting in two different boundary conditions. When in contact with the sea, z b (x, t) > b(x), the same conditions as at the calving front are applied: Where the ice is in contact with the bedrock, z b (x, t) = b(x), different conditions are applied depending on whether the ice is about to lift off from the bed or not: either the compressive normal stress −σ nn is larger than the water pressure that would otherwise exist at that location and the ice has zero normal velocity, or alternatively the normal velocity is only constrained not to point into the bed. This leads to a pair of complementary equation/inequality pairs as in Schoof (2005); Gagliardini et al. (2007); Durand et al. (2009a) and the free boundary between regions that are about to lift off the beds and ones that are not must be determined as part of the solution. In either of the cases in Eq. (10), we specify a nonlinear friction law as where C and m are parameters of the friction law (Table 1). Lastly, a kinematic boundary condition determines the evolution of upper and lower surfaces: where a j is the accumulation/ablation (melting) term, with a b = 0 and j = (b, s).

Approximations to the Stokes flow model
The model above is valid in the absence of any further approximations and, therefore, represents the most accurate mathematical description of marine ice sheet dynamics within the intercomparison exercise. Numerical models labelled as full Stokes (FS) below solve this full system of equations (Durand et al., 2009a,b). Owing to the considerable computational effort, approximations to these equations are often used, such as higher-order, shallow-shelf and shallow-ice approximations. They involve dropping terms from the momentum balance equations as well as simplifying the strain rate definitions and boundary conditions. Higherorder Blatter-Pattyn type models (HOM) consider the hydrostatic approximation in the vertical direction by neglecting vertical resistive stresses (Blatter, 1995;Pattyn, 2003;Pattyn et al., 2006;Schoof and Hindmarsh, 2010). A further approximation, known as the shallow-shelf approximation (SSA), is bed elevation (m a.s.l.) q ice flux q g grounding line ice flux u = (u x , u z ) ice velocity obtained by neglecting vertical shear (MacAyeal, 1989). This is valid for ice shelves and ice streams characterised by a low basal drag and is the one used in the prototype experiments for this intercomparison (Schoof, 2007a). The most common approximation in large scale ice dynamics simulations is the shallow-ice approximation (SIA). This approximation incorporates only vertical shear stress gradients opposing the gravitation drive, which is valid for an ice mass with a small aspect ratio (i.e., thickness scale much smaller than length scale). Its main advantage is that all stress and velocity components are locally determined. However, the approximation is not valid for key areas such as ice divides and grounding lines (Hutter, 1983;Baral et al., 2001); models that remedy this through the use of grounding line flux or grounding line migration parameterizations based on matched asymptotics (Schoof, 2007b(Schoof, , 2011 are classified as asymptotic models (ASY) here. An overview of different Stokes approximations and their validity is given in Hindmarsh (2004). SSA, SIA and asymptotic models based on FS/SSA are described in more detail in Appendix A, and our terminology is further clarified in Table 2. Note that the only forces applied to the shelf downstream of the grounding line are due to hydrostatic pressure from the seawater. For vertically integrated models such as those based on the shallow-shelf approximation, this allows direct application of the sea-water pressure at the grounding line, and in consequence it is not necessary from a mechanical point of view to model the shelf. In consequence, many of the simplified models in the intercomparison do not have a shelf included.

Discretization
The most important numerical issue in marine ice sheet models is that the grounding line is a free boundary whose evolution must somehow be tracked. There are several numerical approaches in ice sheet models to simulate grounding line migration: fixed grid (FG), stretched or "moving" grid (MG) and adaptive techniques (Docquier et al., 2011). They essentially differ in the way grounding lines are represented. In FG models, the grounding line position is not defined explicitly, but must fall between grid points where ice is grounded and floating. Large-scale ice-sheet models (Huybrechts, 1990;Ritz et al., 2001) used this strategy to simulate grounding line migration. Moving grid (MG) models allow the grounding line position to be tracked explicitly and continuously by transforming to a stretched coordinate system in which the grounding line coincides exactly with a grid point (Hindmarsh, 1993;Hindmarsh and Le Meur, 2001). Adaptive grids apply a mesh refinement around the grounding line without necessarily transforming to a coordinate system in which grounded ice occupies a fixed domain as in an MG model (Durand et al., 2009a). Adaptive refinement, for instance, allows grid elements to be refined or coarsened depending on the level of resolution required to keep numerical error locally below an imposed tolerance (Goldberg et al., 2009). Some FG models, while not using a refined mesh around the grounding line, can be adapted in such a way that sub-grid grounding line position and migration can be achieved through local interpolations. For instance, the FG model of Pattyn et al. (2006) determines grounding line position at sub-grid resolution using a flotation criterion, while the adapted grid (AG) models in Gladstone et al. (2010a) use several interpolation schemes in combination with locally increased resolution for the same purpose.
Previous studies have indicated that it is necessary to resolve the transition zone/boundary layer at sufficiently fine resolution in order to capture grounding line migration accurately. In large scale models, this can lead to unacceptably small time steps and costly integrations. Pollard and DeConto (2009) incorporated the boundary layer solution in Schoof (2007a) directly in a numerical ice sheet model at coarse grid resolution by applying a heuristic rule: if the analytically computed boundary layer flux across the actual grounding line q g from (Eq. A12, see Appendix) is greater than the modelled flux through the last grounded grid point q i , then q g is imposed at that grid point (Docquier et al., 2011). Otherwise, q g is imposed one grid point further downstream (i.e., the first floating grid point). The former is usually associated with grounding line retreat and the latter usually with grounding line advance: The Cryosphere, 6, 573-588, 2012 www.the-cryosphere.net/6/573/2012/  The asymptotic models ASY are described in Appendix A3. These are SIA models that include a systematic parameterization of grounding line migration based on near-grounding line boundary layer in which the shallow ice approximation fails. An SIA/SSA model in which the SSA is used to compute basal sliding rates, and a SIA correction is applied to compute flux due to vertical shearing (Bueler and Brown, 2009 (2008) a FPA1, FPA2, FPA3: non-staggered grid, b FPA4, FPA5: staggered grid, c FSA1: uses a 4th order difference scheme for computing velocity gradients, d FSA2: uses a 2nd order difference scheme for computing velocity gradients, e RGL: while many more model results were submitted, we limit the analysis to only two models.
The above-described Pollard and DeConto heuristic is implemented in a number of models that we describe as "fixed grid with heuristic" models (FGH). When applied to SSA models, we identify those as being "SSA-H" as introducing the heuristic (Eq. 14) changes the mathematical structure of the model. When applied to an SIA model, the heuristic Eq. (14) in fact results in a discretization of the asymptotic (ASY) model B * described in Appendix A3, and we consequently classify these as being of type ASY. A complete list of participating models and their characteristics is given in Table 2. Note, however, that most models are SSA and SSA-H, which makes it difficult to compare different approximations to the Stokes equations (as opposed to different discretizations of the same approximation). How-ever, as is shown below, discrepancies between models appear to be dominated by the discretization method chosen, indicating that some of these may not perform well at solving the differential equations on which they are based.

Experimental setup
All information and documentation concerning the MISMIP experiments can be found on the MISMIP website (http: //homepages.ulb.ac.be/ ∼ fpattyn/mismip). The experimental specification document is also included as supporting online material.

Relaxation to steady state on a downward-sloping bed (EXP 1)
The theory in Weertman (1974) suggests that, for a fixed accumulation rate, there should be a single stable equilibrium profile for a marine ice sheet along a flowline on a downward sloping bed (Fig. 1). To compare model output against this result, we used a simple bed shape with a constant downward slope given by Schoof (2007b) b where b is the bed elevation (m a.s.l.). In order to attain different equilibria, Glen's flow law rate factor A is varied as shown in Table 3: larger A corresponds to smaller ice sheets. Table 1 lists the symbols used for various parameters that are prescribed, and the numerical values to be used for these parameters are given in Tables 3 and 4. We use fixed parameter values for accumulation rate a, ice and water densities ρ i and ρ w , acceleration due to gravity g and Glen's law exponent n in Table 4. The sliding law exponent (m = 1/3) and the coefficient C in Table 4 is chosen to give sliding velocities around 35 m a −1 for a driving stress of 80 kPa. Each model run (with the values of A and C, m as described above) starts from a ten metre layer of ice on land, extended up to the position where this ten metre layer becomes afloat (x = 694.5 km). Models are run to steady state for each setting of A (Fig. 1).

Reversal of parameter changes (EXP 2)
The discussion about neutral equilibria in marine ice sheets motivates EXP 2, which is a follow-up to EXP 1. We wish to establish whether the ice sheet relaxes back to the same equilibrium profile during a retreat that follows the advance. EXP 2 starts from the final ice sheet profile obtained from EXP 1. With this initial condition, A is set to the value in run no. 8 in Table 3, and the ice sheet is allowed to relax to a new equilibrium. Subsequently, A is reduced to the value in run no. 7 in Table 3, and again the ice sheet is allowed to relax to equilibrium. This proceeds until the value of A in run no. 1 in Table 3 is attained. Sample results for both EXP 1 and 2 from one of the participating models are displayed in Fig. 1.

Hysteresis (EXP 3)
The results in Schoof (2007a) indicate that ice sheets can undergo hysteresis under parameter variations when the shape of the bed has an overdeepening. This experiment aims at assessing whether other ice sheet models exhibit the same behaviour, and to assess further how transients differ between different models and discretizations. We deliberately use step changes in parameters to isolate the threshold behaviour expected (and to allow comparison with semi-analytical steady states by allowing relaxation to steady state following each step change). We use the same overdeepened polynomial bed shape as in Schoof (2007a) (Fig. 2): The values for A to be used in EXP 3 are given in Table 5. The experiment is conducted as follows: as in EXP 1, the model domain is set up with a 10 m thick grounded ice layer up to the location where this becomes afloat (x = 479.1 km). The model is then run with the value A given for "step 1" in Table 5, as appropriate, for the time interval given for "step 1". Subsequently, the value of A is changed to that given for "step 2" thereby leaving the ice sheet geometry unchanged. Step no. A (s −1 Pa −3 ) time interval (years) 1 3 × 10 −25 3 × 10 4 2 2.5 × 10 −25 1.5 × 10 4 3 2 × 10 −25 1.5 × 10 4 4 1.5 × 10 −25 1.5 × 10 4 5 1 × 10 −25 1.5 × 10 4 6 5 × 10 −26 3 × 10 4 7 2.5 × 10 −26 3 × 10 4 8 5 × 10 −26 1.5 × 10 4 9 1 × 10 −25 1.5 × 10 4 10 1.5 × 10 −25 3 × 10 4 11 2 × 10 −25 3 × 10 4 12 2.5 × 10 −25 3 × 10 4 13 3 × 10 −25 1.5 × 10 4 The model is run for the time interval given for "step 2". This is continued until the end of the table is reached. The different time intervals have been chosen to account for the fact that some changes in A result in irreversible transitions across the overdeepening, and these transitions can be anticipated to take longer to complete than the smaller ice sheet adjustments that occur for other steps in A (see also the sample results in Fig. 2).

Grid spacing and time steps (Mode M1, M2 and M3)
Each of the experiments EXP 1-3 are conducted at three grid resolutions. Mode M1 consists of a computational domain of 100 equally spaced horizontal grid points ( x = 12 km). Nested grids or local grid refinements are not allowed in this mode. However, a stretched coordinate system of the form σ = x/x g is permitted for MG models, in which case a uniform grid in σ with 100 grid points is applied. Mode M2 is similar to Mode M1, but with 1000 equally spaced horizontal grid points in the domain ( x = 1.2 km). Finally, Mode M3 has a grid definition at will. This mode allows for local grid refinements near the grounding line, which may be crucial to obtain accurate results. The aim of this mode is to produce the most accurate results for every user. Time steps are specified by the user as needed to satisfy stability criteria. They are chosen such that results are insensitive to changes in the time step. (EXP 1 and 2)  for EXP 1 and 2. We deliberately left out the few other models that have other approximations to the Stokes flow problem (SIA, SIA/SSA, HOM, FS, Asymptotic) in order to facilitate the comparison. An analysis of the other physical approximations is reported in the sections on the overdeepened bed experiment (EXP 3) and the discussion. A reassuring observation is that most models are capable of producing the advance of the grounding line when a stepchange perturbation in flow parameter A is invoked. The only models that fail to advance the grounding line (or to cause it to retreat in EXP 2) are fixed grid models at too low resolution (red circles in Fig. 3, Mode M1). At higher spatial resolutions, FG models simulate an advance/retreat of the grounding line in closer accordance with the boundary layer theory (Fig. 3, Mode M2 and M3), albeit that retreat simulations show less agreement. Models that incorporate adaptive grid refinement (AG models) generally agree better with the theoretical steady state solutions than FG models, which may be attributed to the higher spatial resolution in the boundary layer.

Steady state analysis for SSA and SSA-H models
More robust results are obtained with MG models. This is not surprising as they follow the grounding line precisely and do not rely on any interpolation of the grounding line be-tween grid points. However, differences between SSA MG models and the theoretical curve are of the order of tens of kilometres. This difference tends to get smaller when a higher resolution is applied.
Results that are the most in agreement with theory are produced with SSA-H FGH models and the HPOMG model. SSA-H FGH models are fixed grid (FG) models that explicitly incorporate the boundary layer theory by setting the grounding line fluxes to the theoretical flux according to Eq. (A12). They are, therefore, not strictly speaking selfconsistent discretizations of the SSA equations, as the latter do not explicitly include the flux relation Eq. (A12). However, as expected, the steady state grounding line position is within metres to tens of metres from the theoretical value required by the boundary layer theory, and this result is independent of the grid resolution applied.
The HPOMG model is a consistent numerical discretization of the SSA equations (without the heuristic "fix" of the SSA-H FGH models) that has very high accuracy within the boundary layer. The agreement between this model and the boundary layer theory is interesting because the latter is based on an asymptotic limit that is not attained by any finite parameter values, and the boundary layer results, therefore, remain approximate, albeit with an error that goes to zero in the asymptotic limit. Consequently, the numerical and theoretical curves need not overlap precisely, but the HPOMG model with its high accuracy in the boundary layer produces results in close agreement with the boundary layer theory, giving added confidence in the latter.
In general, models converge better to the steady state Schoof (2007a) solution with increasing grid resolution (e.g., Mode M2 and M3), but the difference between each model steady state grounding line position and the Schoof (2007a) solution largely depends on the numerical approach (order of hundreds of meters for SSA-H FGH and HPOMG models; order of 50-100 km for MG and AG models; more than 200 km for FG models).
According to the theory in Schoof (2007a,b), steady grounding-line positions during the retreat experiments should be exactly the same as during advance, as steady state grounding line positions are unique on a downwardsloping bed. Any difference should then be a numerical artifact, most likely due to under-resolution. It appears that FG and AG models are most susceptible to this; for coarse grid sizes, grounding line advance and retreat is not always guaranteed. Agreement between advance and retreat in these models improves when smaller grid sizes are applied. SSA-H FGH models, with their explicit incorporation of the flux formula Eq. (A12) very accurately reproduce the same advance and retreat steady state grounding line positions whenever a staggered grid is used. For each model we can calculate the difference between the starting steady state of EXP 1 for the first value of A in Table 3 and the final steady state of EXP 2 for the same value of A. Both grounding line positions should be the same. The difference is of the order of 10 m for SSA-H FGH and SSA HPOMG models; 100-1000 m for MG models; tens of kilometres for AG models and hundreds of kilometres for FG models. It is noteworthy that SSA-H FGH models are accurate irrespective of the grid size used, which is not the case for SSA FG or AG models.
While most fixed grid (FG) models perform well in advancing the grounding line, retreat only occurs when the resolution is sufficiently high (EXP 1 vs. EXP 2). This was already observed in earlier intercomparisons of grounding line motion with FG models (Huybrechts et al., 1998). Therefore, large-scale FG models of the Greenland and Antarctic ice sheet can be questioned to accurately reproduce grounding line retreat, unless refinements at the grounding line are applied (such as AG refinements).

Transient behaviour for SSA and SSA-H models
The Model B boundary layer theory developed by Schoof (2007aSchoof ( ,b, 2011 is valid over long time periods. However, the theory omits brief transients in which the ice surface inside the boundary layer adjusts after changes in parameters such as A and we expect models that do not use the Schoof formula (i.e., all models apart from SSA-H FGH) to differ from the theoretical curve in that respect. To compare the different models, we plotted the transient change in ice flux and grounding line thickness for each SSA and SSA-H model as well as the same curve for the boundary layer theory in a logarithmic plot (Fig. 4). According to the boundary layer theory Eq. (A12), the relation between grounding line ice thickness and flux takes the form (17) Since the experiments all predict a continuing advance into deeper water, we replace the time coordinate by the thickness and plot flux against thickness as a representation of the evolution. Using the values for m and n according to Table 4, the exponent on the right-hand side equals (m + n + 3)/(m + 1) = 4.75. In a logarithmic plot, this relationship results in a straight line with a slope corresponding to this exponent of 4.75 (black line in Fig. 4). The different model results are plotted every 50 yr. Therefore, dots that lie close together represent a solution that converges toward a steady state. By contrast, dots lying far apart along the line represent an abrupt change, generally corresponding to the initial sudden change in A at the start of each step in the experiment. Since the experiment describes a grounding line advance, the time evolution in Fig. 4 can be established by following the dots from left to right along each curve. To make a sensible evaluation possible, only SSA models capable of producing an advance of the grounding line following the forcing of EXP 1 are considered. Furthermore, we only selected the model resolution that produces the "best" results (either Mode M2 or M3) 1 .
Again, reassuringly, all models produce a similar relationship between grounding line flux and ice thickness with an overall power-law dependence of flux upon thickness (Fig. 4). However, upon closer inspection, the slope of the curves -equivalent to the exponent in Eq. (17) -differs slightly between the different numerical approaches. As expected, SSA-H FGH models follow the theoretical line with slope 4.75, since they are prescribed to do so by Eq. (17). Apart from an initial transient after a step change in A, the SSA HPOMG model closely follows the predictions of the boundary layer theory. This initial transient is expected and we can be relatively confident that the SSA HPOMG model produces correct transient behaviour. Other MG models exhibit similar behaviour, but differ more noticeable from the predictions of boundary layer theory. Presumably, these initial transients are associated with the readjustment of ice geometry in the boundary layer described above, during which the grounding line and, hence, ice thickness h g have not changed significantly yet. During that phase, ice flux differs considerably from the theoretical prediction. Similar behaviour was previously found in Schoof (2007a), Fig. 10. Other SSA models also deviate from a straight line in Fig. 4 by displaying spikes and these seem to depend on the numerical method employed. SSA-H FGH models fit much closer to the theoretical line than other SSA models as expected. Deviations from the straight line for SSA-H FGH models are not necessarily associated with sudden changes in A (only the larger spikes are), but related to grid jumps as the grounding line advances of the finite difference grid. Similarly, SSA FG models are largely biased by grid jumps, a bias that is reduced with increasing resolution (Fig. 4) and, therefore, indicative of numerical error. The same goes for AG models, which are also largely biased by the grid resolution (a higher resolution gives better results).

Hysteresis (EXP 3)
The main purpose of EXP 3 is to find out whether numerical marine ice sheet models fail to reproduce steady state grounding line positions on an upward-sloping bedrock, as shown by theoretical considerations (Weertman, 1974;Schoof, 2007b). Forcing the model by decreasing A (in-creasing viscosity) leads to an advance of the grounding line across a downward sloping bedrock until the upward-sloping bed portion is reached. When the ice flux at the grounding line becomes sufficiently high, the grounding line will traverse this "unstable" sector until a new steady state position is found on the next portion of downward sloping bedrock. The retreat part of the experiment should then inevitably lead to hysteresis, as a sufficient reduction in grounding line ice flux (increase in A) is needed to make the jump over the "unstable" part again. This is illustrated in Fig. 2. In view of this hysteresis, the advance and return path of grounding line positions is no longer along a single line, but follows an Sshaped curve (Fig. 5), characterised by two stable branches (the upper and lower branch of the S-curve) and an unstable one (the central portion of the S-curve). Steady state grounding line positions are expected to be found only on the upper and lower branch of that curve.
Contrary to the previous section (EXP 1 and 2), we do not limit our analysis to SSA models only, but consider the entire range of both physical and numerical approximations. Furthermore, only the "best" (Mode M2 or M3) result of each of the models was retained. However, not all participating models carried out this experiment. Results of those that did so are displayed in Fig. 5. Obviously, the same conclusions apply as to the previous section, but here we will focus on the ability of the different models to produce the jump across the "unstable" part, leading to hysteresis. We also anticipate that different models will have different thresholds. This would be an important observation as it indicates that observations of thresholding behaviour in numerical simulations of real ice sheets may be subject to significant uncertainty.
Most models are capable of producing the advance and retreat of the grounding line and in doing so generate a hysteresis loop (Fig. 5). Only FG models (SIA and SSA) at too low horizontal resolution (> 1 km) are not capable of even advancing the grounding line across the unstable part. None of the participating SIA models is capable of reproducing the hysteresis. For those models that do advance across the overdeepening, the exact position at which the "jump" across the unstable branch happens is again model dependent. This dependency can be largely attributed to the discretization, as the bedrock slope changes in space (and was constant in EXP 1 and 2). This is one of the reasons why low resolution results are discarded in this analysis (for M1 the spread is simply too large). For the high-resolution models, the jump is generally at the position where theory predicts it.
It remains very difficult at this point to distinguish between the performances of the different physical approximations. The participating AG full-Stokes model produces the hysteresis loop as predicted, but the exact location of the steady state grounding line position as well as the location of the "jump" is different from the position predicted by boundary layer theory. This could be either due to the difference in physical model, but also due to numerical under-sampling in the vicinity of the grounding line.

Discussion and conclusions
For the first time, a large-scale intercomparison exercise has been conducted to test numerical marine ice-sheet models. Verification of the results was done with a semi-analytical solution based on boundary layer theory (Schoof, 2007a,b). The intercomparison comprised a large variety of models, both on the level of physical approximations to the Stokes flow, as on the level of numerical approaches (fixed grid, moving grid, adaptive grid). However, despite the large number of submissions, the majority of models were of the SSA type (Shallow-Shelf or "shelfy stream" Approximation), which hampers a detailed comparison of different levels of approximation to the Stokes equation (only one full-Stokes models based on the AG method participated). This shifted the scope of the present paper toward discussing the impact of different numerical approximations on grounding line migration and marine ice sheet instability.
Moving grid (MG) models are probably the most reliable from a numerical point of view, since the grounding line is part of the solution and no interpolations are involved; all other models introduce interpolations of the grounding line between grid points, and results will be biased by the imposed grid resolution. However, MG models are not so easily translated into two-dimensional planform models of grounding line migration, hampering their use for large scale simulations of Antarctic and Greenland ice sheets.
While MG models give the most accurate steady state results on a downward-sloping bedrock, FG models give comparable results as long as the grid size is sufficiently small. Adapted grid models perform in a similar way, but are less time consuming since they use smaller grid sizes for only a part of the domain, around the grounding line. For both model types the horizontal grid resolution should be of the order of hundreds to several hundreds of metres. Models that produce steady state grounding line positions in close agreement with the boundary layer theory are not surprisingly the FGH models, in which the boundary layer theory is used to control fluxes at the grounding line.
As for time-dependent responses, the boundary layer theory developed by Schoof (2007b) is expected to be valid except for a short transient. In the MISMIP experiments, sudden jumps in A drive the migration of the grounding line. This leads to transient discrepancies between predicted and computed fluxes for approximately 50 to 100 yr, and these are expected to result from the adjustment of ice geometry in the boundary layer to the new forcing. For slower and, therefore, possibly more realistic, perturbations, these transients discrepancies may well be smaller. While the boundary layer theory predicts a simple power-law relationship between grounding line flux and ice thickness, transient effects produce a deviation from this function. The MG models produce the most robust results, and all participating SSA models produce similar transient results. SSA-H (FGH) models fail to reproduce the short transients following step changes in forcing parameters, as they use boundary layer theory results to control fluxes across the grounding lines. Instead, deviations from the boundary layer theory are then numerical artifacts due to jumps in the grounding line from one grid point to another. Similar effects due to jumps between grid points are observable in FG and AG models that (unlike SSA-H FGH models) do not explicitly use the boundary layer theory. The effect of these grid jumps is considerably reduced in amplitude when a sufficiently small grid size is employed near the grounding line, typically with spacings < 200 m.
The only full-Stokes model that participated in the intercomparison is the Elmer/Ice model (Durand et al., 2009a). This is an AG model, which like other AG models has the grounding line lying between grid points. As shown in (Durand et al., 2009b,a), results from the model are grid size dependent, with node spacings of less than 100 m necessary to obtain self-consistent results. However, with the smallest grid spacings used, the steady state grounding line positions in this model still lie several tens of kilometres from the position predicted by boundary layer theory, which is considerably further than the discrepancy between highly resolved SSA models and the boundary layer theory. This was already observed in Durand et al. (2009a). The obvious conclusion is that full-Stokes models include all the higher order corrections neglected by SSA models, and that this leads to a noticeable correction in the grounding line position of tens of kilometres. Importantly, the boundary layer theory can be derived either directly from the Stokes equations (Schoof, 2011) or from the SSA equations (Schoof, 2007a), and may be a better approximation to the latter for the parameter choices used in the MISMIP experiments. Schoof (2011) discusses, in further detail, the additional restrictions on parameter choices required to make the boundary layer theory a good approximation to a full-Stokes model, as opposed to the looser restrictions required to obtain a good approximation to an SSA model. These results may seem to indicate a need to solve the full-Stokes equations near the grounding line to obtain fully accurate results. Whether this is a worthwhile undertaking in view of the additional cost is a question of available computing resources, and of other uncertainties in the model: for instance, no matter whether a full-Stokes model or a simpler one such as the SSA or the boundary layer theory are employed, assumptions have to be made about material properties, and these may introduce bigger errors than the choice of full Stokes versus simpler ice flow models. This is particularly relevant with regard to sliding behaviour, whose physics is likely to be poorly constrained by observations for many ice flow simulations.
Nonetheless, an appropriate representation of grounding line migration in numerical models is essential. Although these large outlet glaciers and ice streams near marine margins often exhibit fast flow dominated by longitudinal stress effects, this does not exclude a significant contribution of vertical shearing. The boundary layer theory, on which we have built the present intercomparison, is based on the assumption of rapid sliding near the grounding line (though this need not imply rapid sliding further upstream see Schoof, 2011), and cannot accurately represent grounding lines in whose vicinity sliding is slow compared with shearing, or even absent altogether. These are also regions in which the SSA can be expected to fail, and for them an alternative theory along the lines of Nowicki and Wingham (2008) (with appropriate scalings developed by Chugunov and Wilchinsky, 1996) is more appropriate. Even in the presence of rapid sliding, we note that the asymptotic Model A (Appendix A3) may be a better representation of transient grounding line dynamics than the flux prescription Eq. (A12) used in the FGH models considered here (Schoof, 2007b).
Next-generation marine ice sheet models should be capable of producing adequate response of grounding line migration to external forcing, and the present-day series of largescale FG models is not ideally suited to this task (IPCC, 2010). Implementation of grounding line migration requires resolving the transition zone at sufficient high resolution, ei-ther by using a moving grid approach (following the grounding line directly) or by sufficient fine sub-sampling around the grounding line. Importantly, SIA models that do not include an explicit condition on grounding line flux or the migration of the grounding line as a moving boundary fail to pass all tests based on the boundary layer theory of Schoof (2007aSchoof ( ,b, 2011.
These experiments are deliberately omitting lateral variations to enable comparison with boundary layer theory. However, a new intercomparison exercise for three-dimensional marine ice sheet models is on its way.

A1 SSA model
The simplest possible model for a marine ice sheet that accounts for the longitudinal stresses that couple ice shelf to ice sheet is the Shallow-Shelf Approximation (SSA). It describes a rapidly sliding ice sheet in which there is no vertical variation in horizontal ice velocity over the thickness of the ice (Muszynski and Birchfield, 1987). Furthermore, normal stress at the base and upper surface of the ice is cryostatic, which allows the normal stress conditions in Eqs. (8) and (10) into a simple flotation condition, see also Schoof (2007bSchoof ( , 2011. For a grounded ice sheet that occupies the region 0 < x < x g , where the grounding line position x g can change over time t, the mass (Eq. 13) and momentum (Eq. 5) conservation are described by Upstream of the grounding line x < x g , ice thickness must exceed or at least attain its flotation value (this results from the constraints on normal stress in Eq. 10 coupled with the approximation of cryostatic stress), while downstream x > x g , the ice is freely floating and, hence, its thickness is less than the maximum flotation value. Hence, At the grounding line itself, x = x g , the ice, therefore, just becomes afloat, while coupling with the ice shelf (assumed unbuttressed) imposes a longitudinal stress: Alternatively, the model can easily be adapted by adding an ice shelf. In principle, the boundary conditions (Eqs. A3 and A4) arise from integrating equations for an attached ice shelf. An alternative formulation is to resolve the attached shelf.
Let the shelf occupy a domain x g < x < x c , where x c is the calving front (or the downstream edge of the domain). Momentum conservation for the shelf then becomes while Eq. (A2) is retained as descriptions for the grounded ice sheet. At the calving front, there is an imbalance between hydrostatic pressures in ice and water due to the buoyancy of ice. This imbalance must be compensated by a longitudinal stress, and Eq. (A4) is applied at x = x c .

A2 SIA model
Even simpler than the marine ice sheet model above Eq. (A1) is the Shallow-Ice Approximation (SIA) that takes into account vertical shearing, but omits longitudinal stresses that couple ice shelf to ice sheet. For a grounded ice sheet that occupies the region 0 < x < x g , where the grounding line position x g can change over time t, the mass (Eq. 13) and momentum (Eq. 5) conservation are described by whereū x is a vertically averaged horizontal velocity given bȳ The simplicity of the model lies in the fact that u x is defined locally and only depends on local geometric characteristics, such as surface slope and ice thickness. However, as formulated, this model lacks boundary conditions at the grounding line (or the ice divide, where we require u x = 0). In models classified as "SIA" here, only a flotation condition is applied at the grounding line x = x g , but no other condition. As Eqs. (A6) and (A7) combine to give a parabolic problem for h, a single boundary condition is not sufficient to describe a moving grounding line and the intercomparison results bear this out. Models of this type labelled "SIA" in the main text and show no migration of the grounding line. It is possible to add an additional condition that controls grounding line migration, and this is done in the asymptotic models described next. Note that adding the Pollard and DeConto heuristic rule (Eq. 14) to a shallow ice model as described above (in order to take into account longitudinal stress coupling across the grounding line) leads to model B * below.

A3 Asymptotic models (ASY)
In this section, we describe four shallow ice-type models combine the simplicity of an SIA model (which does not require an elliptic equation to be solved for velocity) with the ability of SSA models to account for longitudinal stress coupling across the grounding line. These shallow ice models are based on the boundary layer theory due to Schoof (2007a). As described in Schoof (2007b), they are in good agreement with numerical results for the SSA model, at least for the particular discretizations employed in that paper. All the shallow ice models consider only the grounded part of the ice sheet 0 < x < x g , and apply a moving boundary condition at x = x g to evolve the grounding line position.

A3.1 Model A
This is model A of Schoof (2007b), and is intended to provide a good approximation to the SSA model. The ice sheet interior is described by a shallow ice-type equation in which ice flux q = u x h arises purely through sliding: Grounding line migration in terms of local bed and ice geometry at x = x g then becomes This is combined with a flotation condition at x = x g :

A3.2 Model B
This is model B of Schoof (2007b), and is similar in its scope and formulation to model A. The only difference from model A is that the grounding line migration prescription (Eq. A10) is replaced by a flux condition q g = A (ρ i g) n+1 (1 − ρ i /ρ w ) n 4 n C www.the-cryosphere.net/6/573/2012/ This flux condition forms the basis of a heuristic rule that is implemented in some ice sheet models (Pollard and De-Conto, 2009;Docquier et al., 2011).

A3.3 Models A * and B *
As shown in Schoof (2011), the boundary layer theory in Schoof (2007a) can be extended to deal with ice sheets that also have some shearing across the ice. This version of the boundary layer theory allows for representations of vertical shear. The equivalent of models A and B in this case are the same, except that the flux prescription in Eq. (A9) is altered to Note that this is the same as the flux hū x computed by the shallow ice model (Eq. A7), and model B * is what the "heuristic rule" of DeConto (2009) andDocquier et al. (2011) implements.