Introduction
The loss of continental ice, especially the large ice sheets of Greenland and Antarctica, makes a substantial contribution to current observed sea-level rise, and one that is accelerating more rapidly than was predicted in the Fourth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC) (Reference SolomonSolomon and others, 2007). Warming-induced ice-shelf loss has caused the flow of major glaciers and ice streams of Antarctica to speed up (Reference Hellmer, Kauker, Timmermann, Determann and RaeHellmer and others, 2012; Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van den Broeke and PadmanPritchard and others, 2012). Although a series of mechanisms for ice loss and associated grounding-line retreat have been identified, making predictive models remains an arduous task and projections of 21st-century sea-level rise require reliable models (Reference Alley and JoughinAlley and Joughin, 2012).
Over the past few years, a new generation of full-stress ice-sheet models has emerged, incorporating the physics needed to reproduce such processes (Reference Nowicki and WinghamNowicki and Wingham, 2008; Reference PattynPattyn and others, 2008; Reference Morlighem, Rignot, Seroussi, Larour, Ben Dhia and AubryMorlighem and others, 2010; Reference Favier, Gagliardini, Durand and ZwingerFavier and others, 2012; Reference Gillet-ChauletGillet-Chaulet and others, 2012; Reference Larour, Seroussi, Morlighem and RignotLarour and others, 2012), as well as being coupled to ocean circulation (Reference GladstoneGladstone and others, 2012; Reference Goldberg, Little, Sergienko, Gnanadesikan, Hallberg and OppenheimerGoldberg and others, 2012a,Reference Goldberg, Little, Sergienko, Gnanadesikan, Hallberg and Oppenheimerb). However, as models become more and more complex, validating them against observational evidence and verifying them against analytical solutions become more difficult to achieve. In particular, the complexity of the interactions at the ice/ocean boundary hampers a correct interpretation of observed thinning rates. Not only is there a lack of direct grounding-line change observations, but also feedback effects related to sub-shelf melting, decreased buttressing and increased basal lubrication all lead to ice-flow acceleration. In turn, flow acceleration may result in grounding-line retreat, leading to further acceleration of inland ice flow, and increased mass flux into the ocean and sea-level rise. Owing to this complexity, direct validation of grounding-line migration computations in numerical ice-sheet models is hampered.
Where direct validation is difficult, model verification is another means that can be accomplished by (1) testing models under simplified conditions against available analytical solutions or (2) comparing models in an intercomparison exercise that may eventually lead to a benchmark (Reference Huybrechts and PayneHuybrechts and others, 1996; Reference PattynPattyn and others, 2008, Reference Pattyn2012).
Recently, Reference PattynPattyn and others (2012) presented an evaluation of grounding-line migration in marine ice-sheet models along a flowline, where results could be directly compared with analytical solutions based on boundary layer theory (Reference SchoofSchoof, 2007a,Reference Schoofb, Reference Schoof2011). The major drawback of such flowline experiments, however, is the difficulty of including buttressing effects, which are considered a major factor in understanding current grounding-line retreat of both Greenland and Antarctic ice sheets (Reference Dupont and AlleyDupont and Alley, 2005; Reference Gagliardini, Durand, Zwinger, Hindmarsh and Le MeurGagliardini and others, 2010; Reference Docquier, Perichon and PattynDocquier and others, 2011).
The ice2sea MISMIP3d intercomparison is based on MISMIP, using a similar experimental set-up, for which analytical solutions under simplified conditions are available. A perturbation experiment is performed whereby a curved grounding line evolves. Resetting the perturbation is a test for the reversibility criterion (as shown by Reference SchoofSchoof, 2007a), predicting the existence of discrete steady states dependent only on accumulation rates, bedrock geometry and parameters describing ice deformation and basal sliding. Model differences and uncertainties are further analysed in a diagnostic experiment based on a geometry produced by a full-Stokes model.
In total, 33 realizations of 17 distinct plan-view ice-sheet models participated in the project. (Realizations of models are distinguished through their main contributor, physical approximation, numerical approach and spatial resolution.) This level of participation is significantly higher than for the ISMIP-HOM intercomparison project for higher-order ice-sheet models (Reference PattynPattyn and others, 2008) and the MISMIP intercomparison for marine ice sheets (Reference PattynPattyn and others, 2012). Moreover, there is a sufficient spread in numerical approaches to allow for a broadly based intercomparison. The next section of this paper describes the basis of marine ice-sheet models and their boundary conditions, and is followed by a description of the experiments. Then results are presented and analysed, and the final section is the discussion of these results.
Model Descriptions
The underlying ice-flow model
The basic problem is to solve the gravity-driven flow of an isothermal, incompressible and nonlinear viscous ice mass. Glen’s law is used as a constitutive equation relating stress to strain rates, i.e.
where τ is the deviatoric stress tensor and Dij are the components of the strain-rate tensor, defined by
where 
is the velocity vector. The effective viscosity, η, is expressed as
where the strain-rate invariant, De
, is defined as 
(using the usual summation convention). Most models add a small factor to De
to keep the term nonzero, and experiments have shown that this does not affect the overall result (Reference PattynPattyn, 2003; Reference CornfordCornford and others, 2013). We use a spatially uniform coefficient A, and set n = 3 (Table 1). The velocity (and pressure field) of an ice body is computed by solving the Stokes problem,
where p is the isotropic pressure, ρ
i the ice density and 
the gravitational acceleration.
List of parameters and values prescribed for the experiments, and some other symbols
The boundary conditions we use are essentially the same as those of Reference Favier, Gagliardini, Durand and ZwingerFavier and others (2012). Although the geometry of the marine ice sheet is three-dimensional (3-D), we consider ice flow essentially in the x-direction, without lateral variations (at least for the standard experimental setup). The ice is, therefore, delimited in the vertical by two free surfaces, i.e. the top interface z = z s(x, y, t) between ice and air, and the bottom interface z = z b(x, y, t) between ice and bedrock or sea. We denote bedrock positive above sea level b(x, y), which is assumed to be fixed. The transverse y-axis is perpendicular to the x-z plane (Fig. 1). The domain is bounded transversally by two lateral boundaries, both parallel to the x-z plane. The length of the ice sheet remains the same over time, which means that the calving front has a fixed position throughout the simulation.
Geometry of the standard experiment (modified, after Reference Favier, Gagliardini, Durand and ZwingerFavier and others, 2012).
At x = 0 the ice divide is a symmetry axis, so the horizontal velocity component is ux (0, y, 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, but that the thickness of the cliff evolves over time (free surface). 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 
be the outward-pointing unit normal to the ice surface and 
the two associated unit tangent vectors. Denoting total normal and shear stress by
we can write boundary conditions at the shelf front as σnn
= p
w and 
. The upper surface z = z
s(x, t) is stress-free, implying that
The lower interface, z = z b(x, y, t), is either in contact with the ocean or the bedrock, resulting in two different boundary conditions. Where in contact with the sea, z b(x, y, t) > b(x, y), the same conditions as at the calving front are applied:
Where in contact with the bedrock, z b(x, y, t) = b(x, y), 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 (Reference SchoofSchoof, 2005; Reference Gagliardini, Cohen, Råback and ZwingerGagliardini and others, 2007; Reference Durand, Gagliardini, de Fleurian, Zwinger and Le MeurDurand and others, 2009a):
and the free boundary between regions that are about to lift off the bed and those that remain in contact must be determined as part of the solution. In the case of Eqn (10), we specify a nonlinear friction law as
where C and m are parameters of the friction law (Table 1). A kinematic boundary condition determines the evolution of upper and lower surfaces:
where 
is the accumulation/ablation (melting) term, with 
and j = (b, s).
The lateral boundaries of the domain are parallel planes (Fig. 1). The first plane (y = W) is an actual border of the domain, and the second (y = 0) is a plane of symmetry, which makes it possible to model half the geometry in the y-direction. In both cases, no flux is considered through the surfaces and the boundary condition prescribed is uy (x, 0, z) = uy (x, W, z) = 0.
Approximations to the Stokes flow model
The model above represents the most complete 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 (Reference Favier, Gagliardini, Durand and ZwingerFavier and others, 2012). 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. Higher-order Blatter–Pattyn type models consider the hydrostatic approximation in the vertical direction by neglecting vertical resistive stresses (Reference BlatterBlatter, 1995; Reference PattynPattyn, 2003; Reference Pattyn, Huyghe, de Brabander and De SmedtPattyn and others, 2006). Reference HindmarshHindmarsh (2004) labels these models as LMLA (multilayer longitudinal stresses). A particular case of this type of model is a depth-integrated hybrid model, combining both membrane and vertical shear stress and is of comparable accuracy to the Blatter–Pattyn model (Reference Schoof and HindmarshSchoof and Hindmarsh, 2010). Vertical shearing terms are included in the calculation of the effective viscosity, but the force balance is simplified compared to LMLA models. Such models are labeled L1 L2, or one-layer longitudinal stresses, using membrane strain rates at the surface computed by solving elliptic equations (Reference HindmarshHindmarsh, 2004). A further approximation, known as the shallow-shelf approximation (SSA), is obtained by neglecting vertical shear (Reference Morland, Van der Veen and OerlemansMorland, 1987; Reference MacAyealMacAyeal, 1989), and called SSA or L1 L1 models, i.e. one-layer longitudinal stresses using membrane stresses at the surface computed by solving elliptic equations (Reference HindmarshHind-marsh, 2004). This is valid for ice shelves and ice streams characterized by low basal drag.
The most common approximation in large-scale ice dynamic 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) in combination with a significant traction at the bedrock. Its main advantage is that all stress and velocity components are locally determined. The approximation is not valid for key areas, such as ice divides and grounding lines (Reference HutterHutter, 1983; Reference Baral, Hutter and GreveBaral and others, 2001), since it excludes membrane stress transfer across the grounding line (Reference PattynPattyn and others, 2012). None of these models participated in this intercomparison. However, a more elaborate type of SIA model is a hybrid SIA–SSA model, in which the SSA model is used as a basal sliding law for the SIA model (with zero basal friction for the ice shelf), so this is essentially an SSA model with some vertical shearing terms involved. Coupling of the two models was first proposed through a heuristic rule as a function of sliding velocity (Reference Bueler and BrownBueler and Brown, 2009) and is now done by simply adding both contributions (Reference WinkelmannWinkelmann and others, 2011). We call this type of model HySSA. Although they include both membrane and shearing terms across the grounding line, they are less accurate than the L1 L2 type of model.
The fact that SIA is not valid at grounding lines is remedied by some HySSA models through use of grounding-line flux or grounding-line migration parameterizations based on solutions obtained using matched asymptotics (Reference SchoofSchoof, 2007b, Reference Schoof2011). They are classified as asymptotic models, i.e. A–HySSA, and are described in more detail below.
Discretization
The fundamental numerical issue with 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, stretched or ‘moving’ grid and adaptive techniques (Reference Docquier, Perichon and PattynDocquier and others, 2011). They essentially differ in the way grounding lines are represented. In fixed-grid models, the grounding-line position is not defined explicitly, but must fall between gridpoints where ice is grounded and floating. Large-scale ice-sheet models (Reference HuybrechtsHuybrechts, 1990; Reference Ritz, Rommelaere and DumasRitz and others, 2001) used this strategy to simulate grounding-line migration. We can divide them into two categories, i.e. models with a regular grid, or fixed regular grid (FRG) models, and those on an irregularly spaced grid with a smaller grid size near the grounding line (either based on finite-difference or finite-element methods), or fixed irregular grid (FIG) models. Some models, while not using a refined mesh around the grounding line, can be adapted in such a way that subgrid grounding-line position and migration can be achieved through local interpolations, thereby altering basal friction near the grounding line at subgrid resolution (Reference Pattyn, Huyghe, de Brabander and De SmedtPattyn and others, 2006; Reference Gladstone, Lee, Vieli and PayneGladstone and others, 2010a; Reference WinkelmannWinkelmann and others, 2011). These models are marked with a ‘+’ sign in Table 2.
List of participating models. Numerical approximations: FRG = fixed regular grid model, FIG = fixed irregular grid model, AG = adapted grid model, PSMG = pseudo-spectral moving grid model, and + denotes that grounding-line interpolation within gridcells is applied. Physical approximations: HySSA = shallow-ice approximation with SSA model as sliding condition, SSA = shallow-shelf approximation, L1 L2 = higher-order model, FS = full-Stokes model. A–HySSA explicitly simulates grounding-line migration using a flux condition at the grounding line, usually numerically approximated with a heuristic due to Reference Pollard and DeContoPollard and DeConto (2009) at the grounding line. The minimum grid size for each model is given in column min(Δx) and the time-step in column Δt
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 gridpoint (Reference Hindmarsh, Morland, Boulton and HutterHindmarsh and others, 1987; Reference Hindmarsh and Le MeurHindmarsh and Le Meur, 2001). Only one MG model participated using a pseudo-spectral method (PSMG). They are, however, harder to implement in plan-view models (Reference Hindmarsh, Morland, Boulton and HutterHindmarsh and others, 1987).
Adaptive grids (AG) 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 (Reference Durand, Gagliardini, de Fleurian, Zwinger and Le MeurDurand and others, 2009a). Adaptive refinement, for instance, permits grid elements to be refined or coarsened, depending on the level of resolution required to keep numerical error locally below an imposed tolerance (Reference Goldberg, Holland and SchoofGoldberg and others, 2009). The main difference between AG and FIG models is that grid re-meshing occurs iteratively with the moving grounding line.
Asymptotic flux conditions
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 (Reference Durand, Gagliardini, Zwinger, Le Meur and HindmarshDurand and others, 2009b; Reference PattynPattyn and others, 2012). In large-scale models, this can lead to unacceptably small time-steps and costly integrations. Reference Pollard and DeContoPollard and DeConto (2009, Reference Pollard and DeConto2012) incorporated the boundary layer solution of Reference SchoofSchoof (2007a) directly in a numerical ice-sheet model at coarse grid resolution, so the flux, q g, across model grounding lines is given by
and a similar equation can be written for q gy . This yields the vertically averaged velocity, u g = q g/h g, where h g is the ice thickness at the grounding line. The middle term in Eqn (13) accounts for back-stress at the grounding line due to buttressing by downstream pinning points or side-shear effects, where τxx is the longitudinal stress just downstream of the grounding line (or τyy in the y-direction), calculated from the viscosity and strains in a preliminary SSA or SIA solution with no Schoof constraints. The ‘unbuttressed’ stress, τ f, is the same quantity in the absence of any buttressing, given by τ f = ρ i gh g (1 − p i/p w)/4 (Reference WeertmanWeertman, 1957; Reference HindmarshHindmarsh, 2006; Reference Goldberg, Holland and SchoofGoldberg and others, 2009), so that
Implementation of Eqn (13) in large-scale models is generally based on a heuristic rule (Reference Pollard and DeContoPollard and DeConto, 2009, Reference Pollard and DeConto2012), but the exact details of the implementation may be model-dependent. In general, if the analytically computed boundary layer flux across the actual grounding line, q g from Eqn (13), is greater than the modelled flux through the last grounded gridpoint, qi , then q g is imposed at that gridpoint. Otherwise, q g is imposed one gridpoint further downstream (i.e. the first floating gridpoint). The former is usually associated with grounding-line retreat and the latter usually with grounding-line advance.
Participating models
The complete list of participating models is given in Table 2. Since the experimental set-up is kept rather simple and the computational cost relatively low, most participants submitted their results at the highest spatial resolution possible. However, those resolutions are generally not applied for prognostic simulations of (parts of) the Antarctic ice sheet, either due to the lack of appropriate data (e.g. BEDMAP has a nominal spatial resolution of 5 km; Reference Lythe and VaughanLythe and others, 2001), or the high computational cost. Only a limited number of models (4 out of 33) were actually applied to simulate (parts of) the Antarctic ice sheet and submitted results for that resolution as well. Of these four models three are published: DPO2 (Reference Pollard and DeContoPollard and DeConto, 2009, Reference Pollard and DeConto2012), TAL5 (Reference WinkelmannWinkelmann and others, 2011) and SCO4 (Reference CornfordCornford and others, 2013). An evaluation of their performance is discussed at the end of this paper. Finally, we note that not all models performed the different tests, as reflected in Figures 4 and 5.
Experimental Set-Up
All information and documentation concerning the ice2sea MISMIP3d experiments can be found on the MISMIP3d website (http://homepages.ulb.ac.be/fpattyn/mismip3d).
Standard experiment, Stnd
The initial set-up is comparable with the MISMIP experiment, i.e. a simple bed shape with a constant downward slope, but with a different slope and origin. There is no lateral variation in y, to allow comparison with the flowline case. The aim of this experiment is to verify the model with analytical solutions from Reference SchoofSchoof (2007a,Reference Schoofb). The results of the run should be similar to those of the flowline case, since this is a laterally extruded version of it. The bedrock elevation is defined by a sloping plane in the flow (x) direction:
where b(m) is the bedrock elevation (positive above sea level) and x is given in kilometres. The standard parameters are summarized in Table 1.
The rectangular domain stretches from 0 to 800 km in x and from 0 to +50 km in y (Fig. 1), in which the ice sheet is grown by applying the boundary conditions described above. Starting from the initial set-up (either an initial slab of 10 m of ice or an extruded version of the converged result of the flowline case), the model is run until a steady state is reached (experiment Stnd). Numerical parameters (e.g. grid size, time-step and integration period) were chosen by each participant (Table 2).
Basal sliding perturbation, P75S and P75R
The perturbation experiment starts from the geometry obtained from experiment Stnd. At t = 0, a perturbation in the basal sliding parameters is introduced, precisely at the grounding line, centred on the axis of symmetry (at y = 0 km), defined by a Gaussian bump, i.e.
where 
is the perturbed basal sliding coefficient, a = 0.75 is the perturbation amplitude (a maximum perturbation of 75%), x
b is the precise position of the grounding line at y = 0 km obtained from Stnd, y
b = 0 km, x
c = 150 km and y
c = 10 km. These parameters determine the spatial extent of the Gaussian perturbation. The model is then run forward in time for 100 years (experiment P75S). The grounding-line displacement according to this experiment is shown in Figure 2 for the Elmer/Ice model. It results in a curved grounding line, with maximum forward displacement along the axis of symmetry and a slight backward motion near the free-slip wall (y = 50 km). While this perturbation can be seen as a sudden large-amplitude event, it did not result in any shock response among the models, and the produced velocities and fluxes at the onset of the perturbation were comparable with those of large ice streams currently losing ice at a fast rate.
Description of the perturbation experiment on the numerical domain: position of the initial (black) and perturbed (red) grounding line produced with the Elmer/Ice model.
Based on Reference SchoofSchoof (2007b), we know that for an ice sheet resting on a downward-sloping bedrock, a perturbation in flow parameters (e.g. basal sliding) is reversible. The P75R experiment aims to prove this for the numerical model solution. Starting from the configuration obtained in the prognostic experiment, P75S (perturbation after 100 years), the perturbation in 
is removed and the original value of C is applied. The model is run to steady state, or at least long enough that the grounding line is stationary. According to the theory, the resulting configuration should be equal to the original set-up obtained with Stnd, or x
g(Stnd) ≡ x
g(P75R). The steady-state criteria differ among the participants. Most models apply a finite time criterion ranging from 10 000 years (DPO) to 30 000 years (FPA, HSE, HGU1, TAL, VUB, MTH, RHI2). For LFA, GJO and TKL, the normalized annual difference of volume between two states is <10−5. For DGO and SCO, a limit on the change in ice thickness is set, depending on the experiment and ranging between 0.043 m a−1 for SCO0 and 10−5 m a−1 for DGO. While some experiments did not reach a steady state in terms of ice thickness change, for all models the grounding line was static. Only RHI1 may not have reached steady state in terms of the definitions above, as the integration time was 400 years for the retreat experiment.
Diagnostic experiment, P75D
The aim of the diagnostic experiment is to directly test the numerical models approximating the full-Stokes equations with a high-resolution full-Stokes model (Elmer/Ice). The experiment uses the P75S solution produced with the full-Stokes model Elmer/Ice, for which each model calculates the flow field corresponding to this fixed geometry. This experiment will be discussed first.
Results
Diagnostic experiment
A test for the performance of the different approximations to the Stokes equations is given by a diagnostic experiment, based on the geometry of the perturbation experiment, P75S. This geometry is provided by the Elmer/Ice model, one of the full-Stokes models, and can be downloaded from the MISMIP3d website. Such experiments enable detection of possible inconsistencies between models and model set-ups. Analysis of normal stresses perpendicular to the grounding line shows extension across most of the grounding line, with the exception of the area near the symmetry axis (Fig. 3): according to the Elmer/Ice model, compression occurs where the grounding line is perturbed furthest downstream across the slippery spot. This result contradicts a premise on which Eqn (13) is built, which is only valid for extensional flow (Reference SchoofSchoof, 2007b) as this is the only case that permits steady state.
Normal stress (MPa) distribution perpendicular to the grounding line with the Elmer/Ice model for the P75S experiment. Ice flow is in the direction of the reader, parallel to the x-axis. Compression is observed near the symmetry axis (x = 0), where the grounding line moves out furthest.
Participating models used this geometry in a diagnostic way to reproduce the 3-D ice-flow field, for which the surface component perpendicular to the grounding line is shown in Figure 4. Generally, the flow velocity is correctly represented by a curve that reaches a maximum flow speed between 700 and 1000 m a−1 on the symmetry axis and between 250 and 500 m a−1 near the free-slip wall. Exceptions are either significantly higher (GJO1, HSE1) or lower flow speeds (MTH1, TAL3, TAL4) at the symmetry axis. These differences do not seem to be related to the physics incorporated, but could be due to the interpolation of the geometry produced by the Elmer/Ice model. The latter is clearly shown by VUB2 and TAL1, with velocities in agreement at the symmetry axis, but a wider spread elsewhere: here the model resolution is too coarse to capture the exact curvature of the Elmer/Ice grounding line. These discrepancies do not, however, occur in the prognostic runs.
Horizontal surface velocity perpendicular to the grounding line according to experiment P75D. This diagnostic experiment is based on a geometry of experiment P75S produced by the LFA1 model. Model types: A–HySSA (red); HySSA (orange); SSA (black); L1L2 (green); FS (dark blue). Models are ranked according to grid resolution within each type. The horizontal scale is the same for each of the participants, but may be shifted in some cases to facilitate comparison. ‘No data’ means that participants did not submit results for this experiment.
Reversibility test
This experiment is performed in order to find out whether participating models produce reversible grounding-line positions under simplified conditions. According to theory (Reference SchoofSchoof, 2007a,Reference Schoofb), a marine ice sheet without lateral variations (absence of buttressing) and resting on a monotone downward sloping bed is characterized by a unique steady-state grounding line. Consequently, perturbing the system and subsequently resetting the parameters to their initial values should result in exactly the same initial grounding-line position in steady state (comparison of results from experiment Stnd with P75R). This is shown in Figure 5, where the black line is the initial position of the grounding line, Stnd, and the light blue line the final position (after resetting the perturbation, P75R). The straightness of the Stnd and P75R lines arises from the lack of lateral variations in the model set-up. Although this is a very simple test that does not challenge the 3-D nature of the models, it is an experiment that may be compared with analytical solutions (Reference SchoofSchoof, 2007b; Reference PattynPattyn and others, 2012). Failure to produce the reversibility is related to a too coarse grid resolution (Reference PattynPattyn and others, 2012), which is the case for TAL3 and TAL4 (Δx=1 km without subgrid interpolation), TAL5–TAL8 (Δx = 16.6 km), FPA1 (Δx = 2.0 km), GJO1 (Δx = 1.5 km) and SCO0 (Δx = 6:25 km). They are therefore excluded from the analysis below. Although the SCO4 model (Δx = 0:39 km) returns better than the previous models, generally higher resolutions (Δx ≤ 0:30 km) are needed to produce convincing results.
Steady-state grounding-line position (Stnd, black line), P75S grounding-line position (red) and P75R position (light blue) in the x-y plane. Perturbation occurs on symmetry axis (x = x g, y = 0). Model types: A–HySSA (red); HySSA (orange); SSA (black); L1 L2 (green); FS (dark blue). Models are ranked according to grid resolution within each type. Models marked with ‘R’ show reversibility of the grounding line after perturbation. The x- and y-axes have the same scale for all models, but to facilitate comparison, shifts in x may occur. ‘No data’ means that participants did not submit results for this experiment.
Most significantly, models that do exhibit reversibility are those that implement asymptotic flux conditions at the grounding line (Reference PattynPattyn and others, 2012; Reference Pollard and DeContoPollard and DeConto, 2012). Their P75R grounding line is generally within hundreds of metres of the Stnd position, irrespective of the horizontal spatial model resolution. The second category of models that show reversibility are models that have a sufficiently small grid size (either nominal or through local refinements near the grounding line). In this case, spatial resolution near the grounding line should be well below 1 km (generally <500 m and preferably <300 m) in order to guarantee reversibility. Finally, a third group consists of models that apply an interpolation of the grounding line at subgrid scale. However, this last category remains hampered by the nominal resolution, as it influences the span of the marine ice sheet; too coarse resolution leads to a too small steady-state ice cap, irrespective of the grounding-line interpolation applied. This is illustrated by comparing FPA1 with FPA2, or SCO0 with SCO4 and SCO6 (same numerical/ physical model, different spatial resolution).
Some models have their P75R grounding line situated upstream of the initial Stnd position (DMA6, RHI1, DGO2, LFA1). This is most likely related to the fact that the model has not yet reached steady state. Both a higher resolution and a longer physical calculation time may lead to a closer match between P75R and Stnd, as indicated by RHI2.
Steady-state grounding-line positions
As shown by Reference PattynPattyn and others (2012), steady-state grounding-line positions for the Stnd experiment are different for all model approximations. Only models that passed the reversibility test were taken into account (flagged ‘R’ in Fig. 5), as results from others may deviate for different reasons. Models that use a prescribed flux condition at the grounding line according to asymptotic theory have their steady-state grounding-line position the furthest downstream, i.e. x g ≈ 600–640 km. Such models are largely based on the shallow-shelf approximation, with no vertical shearing. Therefore, the SSA model is the next category having a grounding-line position further downstream, albeit less than for the previous category, i.e. x g ≈ 600–620 km (Fig. 5).
Introducing vertical shearing reduces the effective viscosity at the grounding line, resulting in faster flow, a smaller ice cap and a grounding-line position further upstream. All models including some effect of vertical shearing are characterized by a smaller span and have grounding-line position x g Pz:i 540 km (Fig. 5). This accounts for the FS, L1 L2 and HySSA models. This result is quite robust and independent of the discretization scheme used, as DMA6 (SSA) and SCO6 (L1 L2) are the same numerical model at the same spatial resolution, and both clearly show this difference in steady-state grounding-line position. Nevertheless, there are two exceptions of SSA models having x g in a position corresponding to a model with lowered viscosity at the grounding line, although vertical shearing is not included (HSE1 and TAL2). We did not analyse this in more detail and the difference may be a coincidence. In fact, four types of models based on three sets of equations are compared here, i.e. SSA, HySSA, L1 L2 and FS, and it is to be expected that they produce different results. However, it is beyond the scope of the intercomparison to analyse why, for instance, a L1 L2 model behaves similarly to a FS model. It is worth noting that the analogous situation of a transversely integrated narrow ice shelf also exhibits lower velocities compared with an equivalent full two-dimensional model (Reference HindmarshHindmarsh, 2012).
Perturbation analysis
All models that passed the reversibility test reproduce a similar response to the basal sliding perturbation, of comparable magnitude and in the shape of the curved grounding line. After 100 years, the advance of the grounding line is most pronounced along the symmetry axis, and is ∼20 km. The direction of ice flow into the more slippery area results in a retreat of ∼5 km observed on the free-slip boundary (y = 50 km). This results in a curved grounding line, produced by all models. While this response is typical for the whole range of models, HySSA, SSA, L1 L2 and FS, the response of A–HySSA models is somewhat different: firstly, their response amplitude is definitely larger (∼25 km). This may be because their grounding line is further downstream than in the other models, resulting in a larger ice sheet, characterized by a higher flux at the grounding line. Secondly, none of these models produces the retreating grounding line on the free-slip wall. Both factors may be related to the way buttressing is introduced in a parameterized way, and the fact that Eqn (13) is not valid for compressive flow.
Transient response
The transient response of the grounding line is analysed by plotting the advance and retreat during the 100 year period of the experiment for the basal sliding perturbation, P75S, and the return to the initial state, P75R, for the axis of symmetry (y = 0) and the free-slip wall (y = 50 km). Results are shown in Figure 6.
Time-dependent plot of grounding-line position during (P75S; red) and after (P75R; light blue) the basal sliding perturbation, on the symmetry axis (y = 0: top curves) and on the free-slip boundary (y = 50 km: bottom curves). Model types: A–HySSA (red); HySSA (orange); SSA (black); L1 L2 (green); FS (dark blue). Models are ranked according to grid resolution within each type. The x- and y-axes are on the same scale for all models, but to facilitate comparison, shifts in y may occur. ‘No data’ means that participants did not submit results for this experiment. ‘No R’ means that the model did not pass the reversibility test and results are therefore discarded.
Computed grounding lines in most models take ∼30 years to reach their maximum position, after which they either remain stable or exhibit a slightly retrograde translation. Exceptions are the FPA2, SCO0, SCO4, SCO6 and LFA1 models, which require a longer time to reach the maximum position (50–80 years). It is not clear what delays their response, but in general SSA models are faster in their response than models including both membrane stresses and vertical shearing. This is illustrated by DMA6 (SSA) and SCO6 (L1 L2); both are the same numerical model, but with different physics, and their response time to the perturbation is different. With the exceptions of MTH1 and MTH2, A–HySSA models produce similar results that are consistent with the boundary layer theory. However, none of them shows a grounding-line retreat near the free-slip wall. Grounding-line migration is fastest with the A–HySSA models, since they cover a longer distance in about the same time compared with the other models. This makes them somewhat more sensitive to grounding-line migration and they might therefore overestimate ice discharge due to, for instance, grounding-line retreat (Reference DrouetDrouet and others, 2013). Moreover, single-cell dithering (the flipping back and forth between upstream and downstream points; Reference Pollard and DeContoPollard and DeConto, 2012) occurs when reaching steady state at higher spatial resolutions, as shown by DPO1, DPO2 and VUB1 (Reference Pollard and DeContoPollard and DeConto, 2012). Besides these differences and those reported above on the initial grounding-line position and the overall migration, all remaining models show a similar transient behaviour, but with a slower response time when more physics is included.
Discussion and Conclusions
For the first time, a comprehensive intercomparison of plan-view ice-sheet/ice-shelf models (both two and three dimensions) was carried out, focusing on the migration of curved grounding lines. Verification of the results was done using a semi-analytical solution based on boundary layer theory (Reference SchoofSchoof, 2007a,Reference Schoofb) under simplified conditions. The ice2sea MISMIP3d experiments are clearly in line with results obtained from a previous marine ice-sheet intercomparison project (MISMIP) in which grounding-line migration in flowline models (absence of buttressing) was investigated (Reference PattynPattyn and others, 2012). The main differences are a larger participation and a larger spread in physical model approximations to the Stokes flow. This allowed us to discern major differences between approximations and their effect on grounding-line position and migration, while the discussion by Reference PattynPattyn and others (2012) focused on the impact of different numerical approximations on grounding-line migration and marine ice-sheet instability.
Reference PattynPattyn and others (2012) report that investigated unstable situations on reverse bedrock slopes could be inherently stable when 3-D buttressing effects were taken into account. Reference Gudmundsson, Krug, Durand, Favier and GagliardiniGudmundsson and others (2012) give examples of stable grounding lines on such reverse slopes, illustrating that marine ice sheets are not unconditionally unstable in two horizontal dimensions. In this experiment we explicitly refrained from using complex 3-D bedrock topography (Reference Gudmundsson, Krug, Durand, Favier and GagliardiniGudmundsson and others, 2012) in order to verify models with semi-analytical solutions, while still incorporating elements of buttressing.
Although clear conclusions (listed below) can be drawn with respect to requirements for plan-view models to compute grounding-line migration accurately, some models produce results that are not completely in line with the majority of the participating models. These differences could be due to coding errors, errors in the set-up of the experiment (parameters), or just plausible results that we cannot address based on the analysis of the submitted data alone. The diagnostic experiment is a way to eliminate obvious errors, as the magnitude of the velocity field should converge to the same order of magnitude.
In any case, we expect that models based on different physics will yield different results, and that results from an SSA model will differ from those obtained with a L1 L2 or a FS model. As stated above, the fact that results of L1 L2 models are in close agreement with results from FS models has not been analysed in detail. It is suspected that the lower viscosity at the grounding line due to vertical shearing increases ice flux and makes steady-state grounding lines appear further upstream than for models that only include membrane stresses in their force budget. Furthermore, it is not expected that L1 L2 and FS models will produce similar results under all circumstances. If a transition from no slip to free slip at the base occurs over a horizontal length scale comparable with the ice thickness, then normal stress in the vertical direction can no longer be expected to be hydrostatic, as required by L1 L2 models (Reference Schoof and HindmarshSchoof and Hindmarsh, 2010). Higher-order models therefore cannot necessarily be applied if such transition zones are present, and the full-Stokes equations may need to be solved. Similarly, the full-Stokes equations may have to be solved close to any contact lines to resolve the velocity field there (Reference Schoof and HindmarshSchoof and Hindmarsh, 2010).
The minimum requirements for a numerical ice-sheet model to cope with grounding-line migration are (1) the inclusion of membrane stresses across the grounding line and (2) a freely moving grounding line that is sufficiently resolved. The steady-state grounding-line position is, among other factors, a function of ice viscosity (Reference PattynPattyn and others, 2012), and lowering the viscosity at the grounding line leads to steady-state ice sheets of slightly smaller spans. Therefore, pure membrane models (SSA, A–HySSA) produce larger ice sheets than models that also incorporate an amount of vertical shearing (FS, L1 L2, HySSA). It is noteworthy to add that A–HySSA models are intrinsically pure membrane models at the grounding line, since they are explicitly forced by a heuristic rule based on the boundary layer theory valid for SSA.
High spatial accuracy, usually achieved by high resolution, is a fundamental requirement to resolve a freely moving grounding line. Low-resolution models (>1 km) permit grounding-line advance, but fail to return to the initial position. Resolving grounding lines requires a horizontal spatial resolution of <500 m (preferably <300 m). However, there are several ways of improving model performance without descending to such small grid sizes. The first technique is to resolve grounding lines at subgrid scale and to scale the basal friction field accordingly (Reference Pattyn, Huyghe, de Brabander and De SmedtPattyn and others, 2006; Reference Gladstone, Lee, Vieli and PayneGladstone and others, 2010a,Reference Gladstone, Payne and Cornfordb; Reference WinkelmannWinkelmann and others, 2011). Nevertheless, the ice-sheet span remains a function of the nominal grid size, and too coarse grid sizes (>10 km) lead to grounding-line positions that are far upstream compared with the higher-resolution solutions. A second technique is to adjust the ice flux at the grounding line using a theoretical value obtained from boundary layer theory, as given by Eqn (13), based on a heuristic rule. This ensures unique grounding-line positions to coincide within hundreds of metres.
Heuristic rules in A–HySSA models, however, force them to react immediately to a condition given by boundary layer theory, valid for steady-state conditions. The boundary layer theory is not valid for a very short transient (Reference SchoofSchoof, 2007a,Reference Schoofb; Reference PattynPattyn and others, 2012). While the theory was intrinsically developed for flowline models, Reference SchoofSchoof (2007b) introduced a buttressing factor, θ, to cope with 3-D effects of grounding-line migration. The applied correction directly compares plain strain membrane stresses across the grounding line with those that are influenced by buttressing effects and calculated by the basic model (HySSA in our case). The parameterization is valid for extensional flow, but as we have shown here, compressional forces may apply in buttressed cases (Fig. 3), invalidating Eqn (13). This may be the origin of the faster response times to grounding-line perturbations, as well as the longer distance over which the grounding line migrates, compared with the other models. While models extended with heuristic treatment of the grounding line guarantee reversibility, they are not valid for short transients, as shown by Reference PattynPattyn and others (2012). However, they may be applied for longer-time-scale evolution of ice sheets (>100 years).
Over the past few years considerable effort has been put in place to develop models that are capable of treating grounding-line dynamics in both two and three dimensions. None of the models that participated in the ice2sea MISMIP3d intercomparison were at hand at the time of the publication of the Fourth Assessment Report of the IPCC (Reference SolomonSolomon and others, 2007). All participating models are capable of producing grounding-line migration to a perturbation in basal sliding and produce coherent results as long as the freely moving grounding line is sufficiently resolved.
The MISMIP3d intercomparison cannot circumscribe the grounding-line positions and migration rates that Antarctic models should produce. We can only verify positions under simplified conditions and against a boundary layer theory, which is in itself an approximation of the Stokes flow. However, those models that are applied in a prognostic way to (parts of) the Antarctic ice sheet should at least produce grounding-line positions and migration valid for simplified conditions in accord with what semi-analytical models suggest (Reference SchoofSchoof, 2007b). Both the grounding-line reversibility criterion and the steady-state position of the grounding line are therefore basic elements in such evaluation, and grid-size sensitivity should be tested in any case, since it influences both. None of these models (DPO2, TAL5, SCO4 or VUB2) completely satisfies the criteria listed above. While DPO2 and VUB2 clearly satisfy grounding-line position and reversibility criteria, their transient response is hampered by the invalidity of the Schoof boundary condition for short transients (Reference PattynPattyn and others, 2012). However, their response may be valid over longer time-spans. TAL5 and SCO4, though, still suffer to a certain extent from under-resolution issues at the grounding line, hampering a correct representation of the steady-state position (TAL5) and the reversibility property (SCO4), although the latter model comes rather close.
Ultimately, the accuracy required of a grounding-line representation in a prognostic Antarctic model depends upon the accuracy required in the results. Our results suggest that the numerical error associated with predicting grounding-line motion can, in practice, be reduced significantly below the errors associated with parameter ignorance (e.g. bedrock properties) and uncertainties in future scenarios (e.g. basal melting underneath ice shelves).
Acknowledgements
This work was supported by funding from the ice2sea project from the European Union 7th Framework Programme, grant No. 226375. This is ice2sea contribution No. 112. LGGE was granted access to the high-performance computing resources of CINES (Centre Informatique National de l’Enseignement Supérieur, France) under allocations 2011- 016066 and 2012-016066 made by GENCI (Grand Equipement National de Calcul Intensif) to perform the Elmer/Ice simulations. E. Larour and M. Morlighem are supported by the NASA Cryospheric Sciences and Modeling Analysis and Prediction Programs. H. Seroussi was supported by an appointment to the NASA Postdoctoral Program at the Jet Propulsion Laboratory, administered by Oak Ridge Associated Universities through a contract with NASA. We wish to acknowledge the helpful comments of J. Johnson and F. Saito, as well as the Scientific Editor R. Greve.
Author Contributions
F.P. wrote the paper. L.P. analysed the different submissions and prepared the figures. G.D., L.F., O.G., R.H., F.P. and T.Z. devised and tested the experimental set-up. All authors contributed to the submission of experimental results.
