Results of the third Marine Ice Sheet Model Intercomparison Project (MISMIP+)

. We present the result of the third Marine Ice Sheet Model Intercomparison Project, MISMIP+. MISMIP+ is intended to be a benchmark for ice-ﬂow models which include fast sliding marine ice streams and ﬂoating ice shelves and in particular a treatment of viscous stress that is sufﬁcient to model buttressing, where upstream ice ﬂow is restrained by a downstream ice shelf. A set of idealized experiments ﬁrst tests that models are able to maintain a steady state with the grounding line located on a retrograde slope due to buttress-ing and then explore scenarios where a reduction in that but-tressing causes ice stream acceleration, thinning, and grounding line retreat. The majority of participating models passed the ﬁrst test and then produced similar responses to the loss of buttressing. We ﬁnd that the most important distinction between models in this particular type of simulation is in the treatment of sliding at the bed, with other distinctions – notably the difference between the simpler and more complete treatments of englacial stress but also the differences between numerical methods – taking a secondary role.

Abstract.We present the result of the third Marine Ice Sheet Model Intercomparison Project, MISMIP+.MISMIP+ is intended to be a benchmark for ice-flow models which include fast sliding marine ice streams and floating ice shelves and in particular a treatment of viscous stress that is sufficient to model buttressing, where upstream ice flow is restrained by a downstream ice shelf.A set of idealized experiments first tests that models are able to maintain a steady state with the grounding line located on a retrograde slope due to buttressing and then explore scenarios where a reduction in that buttressing causes ice stream acceleration, thinning, and grounding line retreat.The majority of participating models passed the first test and then produced similar responses to the loss of buttressing.We find that the most important distinction between models in this particular type of simulation is in the treatment of sliding at the bed, with other distinctions -notably the difference between the simpler and more complete treatments of englacial stress but also the differences between numerical methods -taking a secondary role.

Introduction
A number of ice-flow models have been developed in the last decade that simulate fast-flowing ice streams and ice shelves as well as larger, slower-moving ice masses.The key difference between this generation of models and the previous generation is their choice of viscous stress balance equations.All ice sheet models are based upon Stokes flow or, more commonly, one of several approximations to Stokes flow (Hind-Published by Copernicus Publications on behalf of the European Geosciences Union.marsh, 2004; Greve and Blatter, 2009;Pattyn and Durand, 2013).Models designed to simulate the creeping flow of continental ice sheets over glacial cycles are typically based on the shallow ice approximation (SIA), which considers only vertical shear stresses, while more complete approximations are needed for ice shelves and ice streams.The simplest model that can be applied is the shallow-shelf/shelfy-stream approximation (MacAyeal, 1989), which includes horizontal normal and shear stresses and requires the solution of vertically integrated, two-dimensional stress balance equations.More complete models include the L1Lx class of vertically integrated models (Hindmarsh, 2004;Schoof and Hindmarsh, 2010), which resemble the SSA in many respects; the higher-order (HO) models (Pattyn, 2003), which require the solution of simplified three-dimensional stress equations; and the complete Stokes models that include all viscous stresses (Le Meur et al., 2004).
There have been several community exercises comparing ice sheet models where ice stream and shelf dynamics are important.These can be divided into two types: exercises involving real-world ice flows, perhaps forced with climate inputs from sophisticated atmosphere and ocean models (Seroussi et al., 2019;Goelzer et al., 2018), and exercises involving idealized settings with simple forcings.This paper describes the results of an idealized exercise, which can be regarded as sequential to three previous exercises.The Ice Sheet Model Intercomparison Project for Higher-Order Model (ISMIP-HOM, Pattyn et al., 2008) quantified the differences between SIA, SSA, higher-order, and full-Stokes models in time independent settings with periodic bedrock topography and slipperiness.The first Marine Ice Sheet Model Intercomparison Project (MISMIP, Pattyn et al. (2012)) considered a time-dependent but laterally unvarying problem and highlighted the technical challenges faced by numerical models of an ice stream with both grounded and floating portions, that is, with a grounding line.Many models failed to reproduce theoretically well understood properties of such systems.Notably, that a laterally unvarying ice stream on a bedrock that slopes monotonically down in the direction of flow has a single equilibrium state where ice flux across the grounding line matches the total accumulation upstream.The second Marine Ice Sheet Model Intercomparison Project (MISMIP3D, Pattyn et al., 2013) extended a subset of the MISMIP experiments to include perturbations with some lateral variation and once again demonstrated that many models of the time did not exhibit the expected unique equilibrium state.It did not, however, feature strong buttressing of the kind that is important in Antarctic ice shelf and ice stream systems.Recent real-world cases include applications to Pine Island Glacier (Joughin et al., 2019) and Thwaites Glacier (Yu et al., 2018) in western Antarctica and in Antarctica as a whole (Gudmundsson et al., 2019;Martin et al., 2019).
MISMIP+ explores the ability of ice sheet models to simulate coupled ice sheet and ice shelf systems where the ice shelf strongly buttresses the flow upstream and to respond to a loss of buttressing caused by ice shelf ablation.All of the experiments are based around an idealized ice stream, adapted from Gudmundsson (2013), with ice sliding into an ice shelf along a bedrock trough with steep walls (Sect.2.1).Part of the bedrock trough is retrograde -it slopes upward in the direction of flow -and key model parameters are chosen so that, in the absence of sub-ice-shelf melting, models should form a stable equilibrium state with the grounding line crossing the center of the channel on this retrograde slope.The equilibrium state is only stable because of lateral variation in the flow field: it is well known that without such stresses stable steady states form only when the grounding line lies on a prograde slope (Schoof, 2007).The experimental design was stated in Asay-Davis et al. (2016); in this paper we recap the design for convenience and report the results from 15 distinct participants, several of whom carried out the experiments with multiple model configurations.

Experimental design
Three groups of experiments (Fig. 1) were carried out.The Ice0 experiments (Sect.2.5) were designed to show that models were close to steady state at the start of the experiments (time t = 0).The Ice1 experiments (Sect.2.6) saw the ice shelves subjected to ablation at the base of the ice, with a simple formula intended to resemble the gross features of melt rates computed by ocean circulation models, with maximum melt rates close to (but not at) the grounding line.Ice shelf ablation also drives the Ice2 experiments (Sect.2.7), but in this case the imposed ablation is concentrated at the calving front and does not evolve over time.

Geometry
The MISMIP+ ice stream is set in a rectangular domain, spanning 640 km in the x direction and 80 km in the y direction.Ice flows in a direction roughly parallel to the x axis from x = 0 toward x = 640 km, so we will refer to x as the along-flow direction and y as the lateral-or across-flow direction.A no-slip boundary applies at x = 0 and free-slip boundaries apply at both lateral boundaries, while calving front boundary conditions apply at x = 640 km.This choice of boundary conditions, together with the bedrock geometry and mass sources and sinks, results in solutions that have mirror symmetry about the lateral center of the ice stream.For that reason, some participants chose to run their experiments in only one half of the domain, a perfectly acceptable practice for the MISMIP+ and MISMP3d (Pattyn et al., 2013) experiments, though not the related ISOMIP+ and MI-SOMIP experiments (Asay-Davis et al., 2016), where ocean circulation results in nonaxisymmetric melt rates.For convenience in this paper, we define y such that −40 ≤ y ≤ 40 km,  in contrast with Asay-Davis et al. (2016), so that the axis of symmetry is the x axis.
The bedrock elevation z b , measured in meters, is given by where , with the parameters value given in Table 1. Figure 2 illustrates the main features of Eq. ( 1): a steep-walled channel around 48 km wide running parallel to the x axis, and a ridge around x = 505 km.Ice flows from the divide at x = 0 toward the calving front at x = 640 km, so that portion of the ridge between x = 390 km and the summit at x = 505 km is retrograde: it slopes upward in the direction of ice flow.

Initial state
Participants were asked to compute a steady state given the bedrock geometry and boundary conditions above, a constant rate of accumulation a = 0.3 m a −1 , and no melting at the ice shelf base.No particular method was prescribed.Participants simply needed to show that their initial state was sufficiently close to equilibrium that model drift in experiment Ice0 was small compared to the response to the perturbations of the Ice1 and Ice2 experiments.An obvious if time-consuming method is to carry out a spin-up, evolving the ice sheet from https://doi.org/10.5194/tc-14-2283-2020 The Cryosphere, 14, 2283-2301, 2020 One point of departure from the previous MISMIP designs is the specification of the initial-steady-state grounding line position, rather than a full set of model parameters.Participants were asked to produce a steady state where the grounding line crossed the center of the channel (y = 0) at a point x gl (y = 0) = 450 ± 10 km, that is, on the retrograde slope.In order to meet this requirement, participants were free to set any value at all for the rate factor A, which affects stresses within the ice (Sect.2.3), and for the basal friction coefficient β 2 , which affects stresses at the ice base (Sect.2.4).This is a key part of the design for two reasons.First, we wanted all models to begin with their grounding line in steady state on the partly retrograde bed to test that all models could achieve this basic result.Second, we wanted to more closely emulate real-world applications of ice sheet models, where the present day geometry of the ice sheet might be well known, but parameters such as A and β 2 are unknown and are found by some kind of calibration.

Englacial stresses
Participants were permitted to choose any approximation at all for the englacial stresses.That said, the shallowshelf/shelfy-stream approximation (SSA, MacAyeal, 1989) is the simplest approximation that includes the horizontal normal and shear stresses that describe the coupling between floating and grounded ice, so that in particular we did not expect submissions based solely on the shallow ice approximation (SIA).See, for example, Hindmarsh (2004); Greve and Blatter (2009) for a discussion of stress approximations.Whichever approximation was chosen, it was expected to involve Glen's flow law.In the general case, strain rate components D ij and deviatoric stress components τ ij satisfy where n = 3 and 2D 2 e = D ij D j i .Participants were free to chose any constant A in order to realize a steady state within the specified tolerance, though a suggested value was given: A = 2.0 × 10 −17 Pa −3 a −1 (6.34 × 10 −25 Pa −3 s −1 ).

Basal friction
Participants submitted results from simulations carried out with one or more of three basal friction laws.All three ensure that the basal friction τ b = 0 for floating ice, but differ in their approach to ice close to flotation, where the effective pressure at the base, N , is low.The simplest of the three is based on the Weertman (1957) rule for sliding over hard beds, with where u is the horizontal ice velocity at the bed.
Although well known, the Weertman (1957) sliding law is certainly not the final word on glacier sliding (Fowler, 2010;Joughin et al., 2019).Taking a pragmatic (and modelcentered) view, it may not be applicable close to the grounding line, where τ b will be discontinuous.The two alternative rules considered here ensure that τ b is continuous (but not differentiable) 1 .Schoof (2005) and, later, Leguy et al. (2014) considered the case of sliding with cavitation, over hard beds, leading to a rule, which, among other features, ensures that |τ b | does not exceed some fraction α 2 of the effective pressure but can be approximated by the Weertman (1957) rule far from the grounding line where N is large.Tsai et al. (2015) considered a case where friction is due to either sliding over hard beds at high N and deforming beds when N is low.The resulting rule is We refer to both modified rules as Coulomb-limited, because the Schoof (2005) and Tsai et al. (2015) rules both imply |τ b | ≤ α 2 N : in other words, that the magnitude of the basal friction cannot exceed that given by a Coulomb law, where α 2 < 1 is a coefficient of friction.Joughin et al. (2019) uses the term regularized Coulomb, reflecting another benefit of this class of rules: that they permit Coulomb sliding but do not insist upon it over the whole domain, which would be problematic for two reasons.The first problem is physical: across much of the ice sheet N 1 MPa, leading to a magnitude of friction larger than observed.The second problem is mathematical, or at least pragmatic: if |τ b | does not increase with |u| anywhere in the domain, then a Dirichlet condition for u must be imposed on at least some part of the horizontal boundary.The modified rules avoid both of these problems by reverting to the familiar Weertman (1957) rule of Eq. ( 3) where N is large (or ice flows slowly in the case of Joughin et al., 2019).In all three cases participants were free to modify the parameter β 2 in order to achieve the desired steady state: the suggested value was β 2 = 10 4 Pa m −1/3 a 1/3 (3.16 × 10 6 Pa m −1/3 s 1/3 ).

The Ice0 experiment
The Ice0 experiment is simply a test of the initial steady state.Simulations ran from t = 0 a to t = 100 a, with zero sub-ice-shelf melt and an upper surface accumulation rate a = 0.3 m a −1 .Since these are the melt and accumulation rates defining the the initial steady state, models should exhibit little or no variability during this process, or, if they do exhibit some variation, it should result in little long-term drift.That is, fluctuations in ice thickness were acceptable provided that the grounding line did not advance or retreat and the ice volume did not grow or shrink to any great extent compared to the other two experiments.Participants were free to produce the initial state by any method; the majority of participants chose to evolve their models with the stated parameters for ∼ 10 ka in order to approach steady state.

Oceanic forcing and the Ice1 experiments
The Ice1 experiments were intended to examine the response of models to intense ablation at the base of the ice shelf, with a spatial distribution reflecting the results of typical cavity circulation models (Asay-Davis et al., 2016).The melt rate, measured in m a −1 , was where z d is the ice shelf draft and z d − z b is the cavity thickness, both measured in meters.The resulting melt rate will generally increase with draft, reaching values of around 100 m a −1 when z d ≈ 700 m, but will vanish at the grounding line where the cavity thickness is zero.In the Ice1r experiment, the melt rate was applied to floating ice over the course of 100 a, starting from the steady state at t = 0.This results in the loss of much of the original ice shelf thickness over the course of the experiment and was expected to result in grounding line retreat, assuming that the ice shelf had a role in buttressing the ice upstream.The two follow-on experiments, Ice1rr and Ice1ra, were designed to start from the end of the Ice1 experiment at t = 100 a and terminate at t = 200 a.Ice1rr is simply a continuation of Ice1r, with the same melt rate applied, while Ice1ra imposes zero melt rate, allowing the ice shelf to thicken, so that the grounding line should readvance.

Calving and the Ice2 experiments
The Ice2 experiments follow the same basic structure as the Ice1 experiments but impose a different melt rate, applied to ice shelf regions only.This choice of melt rate in Ice2 results in something like a large calving event, where a downstream portion of the ice shelf, amounting to around half the total area, is removed over a short period of time.
The majority of the ice removed lies in a zone that provides little buttressing.By allowing a thick ice shelf to form in the wake of the retreating grounding line, the Ice2 experiments test a model's ability to form a new stable steady state with the grounding line on a retrograde slope.They also test the numerical implementation of the model, because there is often an abrupt increase in melt rate immediately downstream from some portions of the grounding line, in contrast to the smooth increase in melt rates seen in the Ice1 experiments.

Participating models
The participating models cover the same variety of englacial stress approximations as the earlier MISMIP (Pattyn et al., 2012) and MISMIP3d (Pattyn et al., 2013) exercises, with each model including some approximation of the horizontal normal and shear stress (membrane stress).The most complete models make use of the full-Stokes equations, but the computational expense entailed by solving the full 3D stress balance equation limits both the number of participants able to run such a model and the number of submissions from those participants, so that there are only two full-Stokes submissions.The most common class of model is based upon a 2D, vertically integrated hydrostatic stress balance equation, either through the shallow-shelf/shelfy-stream approximation (SSA), which neglects shear strains D xz and D yz above the bed, or an approximation that assumes a simple form for D xz and D yz , for example, the models of Schoof and Hindmarsh (2010) and Goldberg (2011).We will follow Pattyn and Durand (2013) in labeling this second class of vertically integrated model "L1Lx".Intermediate in complexity are the higher-order (HO) or first-order models, which include the hydrostatic Blatter-Pattyn models.HO models exploit the low aspect ratio of ice sheets to reduce the four 3D Stokes equations to a pair of 3D equations in the horizontal velocity components.Finally, the HySSA models solve the 2D SSA stress balance equation but adapt that analytic expression (Schoof, 2007) derived for flow-line models with no buttressing to the general case by means of a heuristic buttressing factor: such models were able to produce results in https://doi.org/10.5194/tc-14-2283-2020 The Cryosphere, 14, 2283-2301, 2020 earlier exercises that, in contrast to all other types, did not depend strongly on mesh resolution.All of the participating models construct and solve their stress balance and mass transport equations though a limited choice of methods.There are several finite-volume methods and finite-difference methods based on rectangular meshes, extruded vertically, and several finite-element methods based on unstructured triangular meshes, also extruded vertically.One model (MALI) takes a mixed approach, combining a finite-element discretization of the stress balance equation with a finite-volume discretization of the mass transport equation.The majority of finite-volume and finite-difference methods employ spatially uniform meshes, with two exceptions: WAVI employs a wavelet-based adaptive grid to reduce the computational expense in solving the stress balance equation, while BISICLES makes use of a time-evolving adaptive block structured mesh in both the stress balance and mass transport equations.The finite-element methods tend to employ spatially nonuniform meshes that do not change over the course of the simulation, with the exception of two of the several Ice Sheet System Model (ISSM) submissions, which update their meshes over time.
Models differ in their discretization of the stress balance equations in the region close to the grounding line.One common class of techniques is the modification of the discretized basal friction term around the grounding line.These techniques have been given at least two different names in the literature -grounding line parameterization (Gladstone et al., 2010;Leguy et al., 2014) and sub-element parameterization (Seroussi et al., 2014b) -but they all represent a similar approach.We will use the terms preferred by the individual model authors when referring to their submissions.The simplest type constructs a piecewise linear approximation to the thickness above flotation (h − h f ) and uses that to evaluate a weight, w ∈ [0, 1], associated with each grid cell or element that reduces the discrete approximation to the friction term accordingly.None of these schemes introduce additional degrees of freedom, and they also have the same order of accuracy with respect to mesh resolution, O( x), as the unmodified scheme.Nonetheless, they have been seen to improve accuracy for a given mesh resolution in several cases (Feldmann et al., 2014;Seroussi et al., 2014a).A related -but more controversial (Seroussi and Morlighem, 2018) -modification to standard methods applies a similar weighting to the basal melt rate term in the mass balance equation.Several participants employed such a modification in these exercises, with clear consequences in the Ice2 experiments (Sect.4.3).
The participating models are described briefly below, ordered alphabetically by model name in Sect.3.1 to Sect.3.11 and listed in Table 2.The supplement also contains a data sheet for each model.

BISICLES
Six submissions are based on BISICLES (Cornford et al., 2013(Cornford et al., , 2015)), a finite-volume model that employs time-evolving adaptive mesh refinement to maintain fine resolution close to the grounding line.Three vertically integrated stress approximations are included: the shallow-shelf approximation (SCO_BISICLES_SSA_Schoof_250m, SCO_BISICLES_SSA_Tsai_250m), the Schoof and Hindmarsh (2010) L1L2 approximation (SCO_BISICLES_L1L2a_Tsai_500m), and a modified L1L2 approximation that includes vertical shear in the effective viscosity but neglects it in the mass flux (the remainder).All three basal friction rules are represented.Mesh spacing at the grounding line is set to 250 m in most cases, but two coarser-resolution cases are included with x ≥ 500 m and x ≥ 1 km.All of the submissions apply a one-sided difference when evaluating the gravitational driving stress at the grounding line, with no other parameterization, although a subgrid friction scheme has been useful in other cases (Cornford et al., 2016).

CISM
Two submissions are based on the CISM (Lipscomb et al., 2013) model.Both employ the shallow-shelf approximation to describe englacial stresses, discretized according to the finite-element method, on a uniform mesh with 1 km horizontal resolution.Thickness transport is effected with an incremental remappping scheme (Dukowicz and Baumgardner, 2000): ice thickness and velocity data are stored at staggered locations.The difference between submissions is the choice of basal friction rule: one submission uses the Schoof (2005) scheme and the other uses the Weertman (1957) scheme.Both make use of a grounding line parameterization (Leguy et al., 2014).
IME_ElmerIce_FS_Schoof_250m is a full-Stokes model, corresponding to the majority of Elmer/Ice publications, for example Durand et al. (2009); Seddik et al. (2012).IME_ElmerIce_L1L2b_Schoof_250m is a vertically integrated model.Both are finite-element models and apply a horizontal resolution of 250 m (indicated by convergence studies to be adequate) over the region swept out by the grounding line during the Ice1 and Ice2 experiments, and a vertical discretization of 7 layers, with finer resolution toward the base.Both models employ the modified basal friction law of Eq. ( 4).

ISSM
Several contributions are based on the ISSM (Ice Sheet System Model), which can treat englacial stresses with a variety of approximations, from the SSA to the full-Stokes equations, using the finite-element method (Larour et al., 2011).The submissions included include one set of full-Stokes simulations (HYU_ISSM_FS_Weertman_500m), HSE_ISSM_SSA_Tsai_1km, HSE_ISSM_SSA_Tsai_500m, and HSE_ISSM_SSA_Weertman_1km all make use of a fixed-in-time, nonuniform mesh of triangular elements that is refined to either 500 m or 1 km across the region swept out by the grounding line and employ the SEP1 subelement parameterization.CBO_ISSM_SSA_Tsai_500m also relies on a fixed, nonuniform mesh but chooses a different parameterization, SEP2 (Seroussi et al., 2014b), when evaluating friction in partly grounded elements.Two further submissions, TDI_ISSM_SSA_Tsai_500m and TDI_ISSM_SSA_Weertman_500m, differ from the others in their use of an evolving adaptive mesh (dos Santos et al., 2019), which is updated throughout the simulations to maintain 500 m resolutions close to the grounding line and coarser resolution elsewhere.

MALI
The MALI (MPAS-Albany Land Ice) (Hoffman et al., 2018) submission treats englacial stresses with the Blatter-Pattyn approximation and basal friction with the Weertman rule.The stress balance equation is discretized horizontally in space on an unstructured mesh of triangular finite elements, while the mass conservation equation is discretized, using the finite-volume method, on the corresponding hexagonal Voronoi tessellation.Time discretization is accomplished with the forward (explicit) Euler scheme.Basal friction around the grounding line is evaluated by computing it and the flotation criterion at the quadrature points of a fifthorder scheme, leading to a treatment comparable to SEP3 of Seroussi et al. (2014b).The experiments were carried out with 500 m spatial resolution, while the Ice1r experiment used resolutions from 4 km to 250 m.

PISM
PISM (Bueler and Brown, 2009;PISM authors, 2019;Bueler et al., 2007) employs a uniform-mesh finite-volume method to discretize the mass conservation equation and a finitedifference method to discretize the stress balance equation.It uses either the shallow-shelf approximation (SSA), or the SSA + SIA approximation, which complements SSA sliding with SIA internal deformation.All eight submissions employ a mesh spacing x = 1 km and differ in subgrid friction interpolation versus none (Feldmann et al., 2014), SSA versus SSA + SIA, Tsai versus Weertman friction rules, and the use or otherwise of an ice thickness transformation, η = H (2n+2)/n , to achieve better numerical results at ice sheet margins (Bueler et al., 2005, Sect. 5.2).All submissions apply a one-sided difference when evaluating the gravitational driving stress at the grounding line.Time integration is explicit, with the time step satisfying both an advection Courant-Friedrichs-Lewy (CFL) criterion determined from the SSA sliding speed and an additional constraint found by expressing the SIA as a diffusion equation (and so proportional to x 2 ).

PSU3D
The two PSU3D (Pollard and DeConto, 2012) submissions represent the only HySSA model that took part in MISMIP+.This class of models is distinct from the more common types of vertically integrated models in their direct imposition of a flux across the grounding line, derived from the analytic expression of Schoof (2007), which provides consistent performance between ∼ 10 and ∼ 1 km.As a result, PSU3D has been able to simulate simulations of Antarctica lasting millions of years (Pollard and DeConto, 2009).It is a finitedifference model, based on a staggered horizontal grid, and Runge-Kutta time integration.Note that the ice-cliff failure mechanisms included in current version of this model (Pollard et al., 2015) do not arise in MISMIP+.

STREAMICE
STREAMICE (Goldberg and Heimbach, 2013) is a physical package of the MITgcm climate model (Marshall et al., 1997).It solves velocities via a finite-element method, using bilinear basis functions on quadrilateral elements on a regular grid; while its thickness is evolved via an explicit finite-volume method.Its vertically integrated stress model is based on Goldberg (2011).Near the grounding line, basal drag is regularized using a sinusoid profile where thickness is within 5 m of flotation (either above or below) -as this treatment is found to yield reversible grounding line movement in the MISMIP3d experiments.However, melt is applied to cells only where cell-averaged thickness is below flotation, i.e., no melting is applied to any cells in which basal drag is nonzero.

TIMFD3
TIMFD3 (Kleiner and Humbert, 2014) solves the stress balance equations with the LTSML (Hindmarsh, 2004) higherorder approximation together with the modified basal friction rule, Eq. ( 5), discretized according to the finite-difference method.In comparison to the Blatter-Pattyn approximation, the LTSML approximation considers the vertical resistive stress R zz (van der Veen and Whillans, 1989) in the momentum balance, and thus vertical longitudinal stresses are not hydrostatic.The MISMIP+ experiments were carried out on a 1 km horizontal grid, with the vertical extent of the ice sheet treated as nine terrain-following layers, more closely spaced at the base.Coarser resolutions do not exhibit a stable grounding line.The initial state is found by taking the output from an SSA model (one of the BISICLES submissions with the same selection of A and β 2 ) and performing a 500-year relaxation from that point.The model does not employ any kind of subgrid interpolation.

Úa
Úa (Gudmundsson, 2019;Gudmundsson et al., 2019) is an open-source finite-element ice-flow model based on a vertically integrated formulation of the momentum equations.The ice-flow equations are solved on an unstructured mesh consisting of linear, quadratic, or cubic triangular elements.Various meshing options are available, including automated mesh refinement and coarsening.When simulating the flow of marine ice sheets these meshing options allow, for example, the areas around grounding lines to be automatically highly resolved as the grounding lines migrate through the computational domain.Elements can be activated and deactivated.This enables the computational domain itself to change in the course of a run, for example when simulating the growth and decay of a large group of mountain glaciers.Ice thickness positivity is enforced using the active set method.Forward time integration can be done in a fully coupled manner, and the resulting nonlinear system is solved using the Newton-Raphson method.

WAVI
WAVI (Arthern et al., 2015) is a finite-volume model that makes use of a wavelet-based adaptive grid to accelerate the solution of the stress balance equation.Its vertically integrated stress model is derived from Goldberg (2011) and treats both membrane and simplified vertical shear stresses.Subgrid interpolation of the basal drag, gravitational driving stream, and sub-ice-shelf melt rates are deployed in the finite volumes immediately adjacent to the grounding line.Two complete submissions are included in this paper, differing only in their uniform mesh resolution.A further two partial submission are included, covering only the Ice2r experiments.These restrict nonzero sub-ice-shelf melt rates to finite volumes whose cell center is floating.

Results
The majority of participants completed each experiment, leading to a large volume of results, many of which are in close agreement.The results of each experiment (Ice0, Ice1, and Ice2) are summarized below, concentrating on the spread of results rather than any individual model.In most cases, a substantial portion of variability in the results can be explained by straightforward groupings of the participating models: for example models that make use of the Weertman (1957) friction rule see slower grounding line migration than models employing either the Schoof (2005) or Tsai et al. (2015) rules.At the same time, distinctions that have been seen to be important in other experiments, such as the use of grounding line parameterization in MISMIP3d (Leguy et al., 2014;Seroussi et al., 2014a;Feldmann et al., 2014), appear unimportant here.

Ice0 experiments
All participating models produced a similar initial grounding line, crossing the channel in a region where the bedrock slopes upward.Figure 3 plots the grounding line for all models with x ≤ 1 km at the start (t = 0 a) and the grounding line position x gl (y = 0, t) in the center of the channel, mean averaged over all models.At the start of the experiments, each grounding line intersects with the channel centerline close to x c = 450 km, and the domain edge close to x e = 500 km, covering most of the distance between x c and x e around the channel walls where the bedrock slopes sharply in the y direction.Most of the models also exhibit a narrow prominence at the top of this lateral bedrock slope, extending downstream from x e .There is little apparent change over time in any model, as required.Not every submission included this test, although the majority did.
Variation between the models is fairly minor.SSA, L1Lx, and HO models that employ the Coulomb-limited basal friction laws are grouped most closely, largely because there was no need to modify the default A = 2.0 × 10 −17 Pa −3 a −1 given in Asay-Davis et al. (2016).SSA, L1Lx, and HO models that include the Weertman basal friction law are more widely spread, with some models increasing A by as much as 25 % to achieve the specified x gl = 450±10 km and others submitting results with the grounding line further upstream.Neither of the two full-Stokes models elected to modify A, but their grounding lines are close to others.Finally, the one HySSA model does produce x c ≈ 450 km, setting A = 3.5×10 −17 Pa −3 a −1 to do so, but differs from the remaining models over the shallower bedrock with x e ≈ 540 km.

Ice1 experiments
Figure 4 plots the grounding line for all models with x ≤ 1 km over the course of the Ice1 experiments.All models see retreat of their grounding line in the center of the channel while the melt rate of Eq. ( 7) is imposed (the Ice1r and Ice1rr experiments) and little change outside the channel walls beyond the erosion of the prominences.The vast majority of the models also see their grounding lines readvance once the melt rate is reduced to zero (the Ice1ra experiment), albeit at a much lower rate so that the initial grounding line is not regained by the end of the experiment at t = 200 a.Although all of the models are closely grouped at the start of the Ice1r experiments, a considerable spread of grounding line contours is evident after 100 years.The median retreat of midchannel grounding line position x gl (y = 0) is in the region of 40 km, but the range is rather larger -close to 100 km.Within groups of models, there is a much smaller spread: around 20 km between the Coulomb-limited models, a rather larger spread between the Weertman models, and an obvious outlier in the one HySSA model with x ≤ 1 km , whose midchannel grounding line position retreats by around 100 km. https://doi.org/10.5194/tc-14-2283-2020 The Cryosphere, 14, 2283-2301, 2020 Grounding lines migrate upstream while the melt rate of Eq. ( 7) is applied until t = 100 a (b).There are two branches for t > 100 a: the Ice1rr branch which continues with the same melt rate (c) and the Ice1ra branch with zero melt rate (d).Yellow contours correspond to SSA, L1Lx, and HO models with the Coulomb-limited basal friction law; orange contours to SSA, L1Lx, and HO models with the Weertman basal friction law; red (ISSM) and cyan (Elmer/Ice) contours to full-Stokes models; and magenta contours to HySSA models.The color map and black contours depict bedrock elevation.
We will present the rest of the Ice1 results in terms of the midchannel grounding line position x gl (y = 0), since it represents the other bulk quantities of interest well enough.The rate of change in x gl (y = 0) over time represents both the rate of change in volume above flotation V f and the rate of change in grounded area A g (Fig. 5).A single proportional relationship links x gl (y = 0) and A g in the Ice1r, Ice1ra, and Ice1rr experiments.Ice1r and Ice1rr also see proportional relationships between x gl (y = 0) and V f , albeit with distinct constants.The exception is the Ice1ra experiment, where the readvance of the grounding line following the retreat in Ice1r is associated with a continued, mild loss of volume above flotation.
The majority of models behave similarly in the Ice1 experiments, so we define a number of subsets to represent the spread of results.The first of these is simply the set of all models, but since there are some clear outliers we also define a main subset comprising models which (1) completed all experiments, (2) place their Ice0 steady-state grounding line at x = 450 ± 15 km2 , and (3) see their grounding lines retreat along the center of the channel to x = 385 ± 35 km. Figure 6 shows the spread in midchannel grounding line position against time for all models, the main subset, and the remainder.Of the remainder, two are the HySSA models, and the remaining two are L1Lx models.We also separate  the two full-Stokes models from the lower-order approximations.The main subset then comprises HO, SSA, and L1Lx models.
One major division within the main subset is the distinction between models that employ the Weertman basal friction law (Eq.3) and models that utilize either of the Coulomblimited rules (Eq. 4 or Eq. 5) (Fig. 7).The mean rate of retreat seen in the Weertman models is around 0.5 km a −1 versus 0.7 km a −1 in the Coulomb-limited models, and indeed the ranges of these two subsets barely overlap.The Weertman subset also exhibits a much larger range, but we attribute this to the larger range of initial states.That larger range is due to only some participants altering the rate factor A to place the initial grounding line within the specified range, rather than any inherent difficulty with implementing the Weertman friction law.
Although the two HySSA simulations achieve an initial state close to the other models, at least in terms of grounding line positions, their transient behavior differs.Figure 8 shows both sets of HySSA results with the main subset and the Weertman subset.Their rate of retreat in Ice1r is around 1 km a −1 , compared to 0.6 km a −1 for the main subset but it is more appropriate to compare to the 0.5 km a −1 of the Weertman subset given that both HySSA calculations employ that rule.Something along the same lines was seen in the MIS-MIP3d experiments (Pattyn et al., 2013).Note though that the HySSA simulations were computed with same model (PSU), albeit at different resolutions (1 and 10 km), so this may not be a typical result.Note also that the 1 and 10 km HySSA simulations are essentially the same, as has been the case for this class of models in other cases (Pattyn et al., 2013).
The two full-Stokes models see their grounding lines retreat at the same average rate as other models while the ice shelf is ablated, but they see a lower rate of advance when the ice shelf regrows.Figure 9 compares the two full-Stokes models with the main subset, and the Ice1r retreat rate of one (IME_ElmerIce_FS_Schoof_250m) lies at the mean of the main subset.The other (HYU_ISSM_FS_Weertman_500m) retreats rather more rapidly to begin with but slows to obhttps://doi.org/10.5194/tc-14-2283-2020 The Cryosphere, 14, 2283-2301, 2020  tain the same position after 100 years -this may be related to this model's position outside the main subset with regard to its initial state.Neither model sees rapid readvance when ice shelf ablation ceases in the Ice1ra experiment.There is little variation between the higher-order and vertically integrated hydrostatic models.Figure 9 shows that the mean rate of retreat in the Ice1r experiment is close to 0.6 km a −1 over 100 years for both higher-order (HO) and vertically integrated (non-HO) models.Likewise, there is essentially no difference between HO and non-HO models in either the mean rate of further retreat in Ice1rr or the mean rate of readvance in Ice1ra.There is rather more variability between the non-HO models than between the HO-models, but that can be attributed to the smaller number of HO models (or the larger number of non-HO models).
At least for the x ≤ 1 km models included in the main subset, the choice between finer resolution or a subgrid friction scheme is unimportant.Dividing the main subset into subgroups with and without a subgrid friction scheme results in the same mean and range of grounding line migration rates, and the same is true if the subset is divided into subgroups with x = 1 km and x < 1 km (Fig. 11).That is not to say that neither a subgrid scheme nor a fine mesh is consequential in general.Apart from the fact that all models in the main subset have x ≤ 1 km, of the 15 models that employ a subgrid scheme, only 5 have x < 1 km, while of the 11 models that do not employ a subgrid scheme, 6 have x < 1 km.The major difference between models in this case is related to numerical methods.The main subset can be split into two subsets: those that exhibit no grounding line retreat along the domain wall and those that see some retreat.
The first group, Ice2A, consists solely of models that do not apply a subgrid interpolation scheme when computing melt rates.The second group, Ice2B, includes models that do apply such a scheme and some that do not but have some treatment special to calving fronts (PISM variants with "eta" in their name).Figure 14 shows the variation in grounding line position at the domain wall between these two groups.Some of the Ice2B submissions have their grounding lines retreat to exactly x = 480 km -the limit of nonzero melt rate in these experiments -along the wall.All of these do employ a subgrid interpolation scheme when computing melt rates and, notably, see a retreat rate that grows with their nominal mesh spacing x.For example, the complete WAVI submissions (with subgrid interpolation of melt rates) exhibit this retreat, whereas the supplementary results (without subgrid interpolation) do not.
The impact of this grounding line retreat along the domain wall, when it does occur, is limited in this case (but is not necessarily limited in others).Figure 15 shows both the grounding line position in the center of the channel and the grounded https://doi.org/10.5194/tc-14-2283-2020 The Cryosphere, 14, 2283-2301, 2020 area as they evolve over the course of the experiment.While the loss of grounded area is significantly greater in the Ice2B models, with the mean rate of loss over the Ice2B group comparable to the maximum rate in the Ice2A group, the vast majority of the additional loss in the grounded area is restricted to the thin ice outside of the channel, and there is little difference between the two subgroups in the center of the channel.

Discussion
One clear result of the MISMIP+ exercise is the degree of similarity between models.The majority of models demonstrated a similar and apparently stable equilibrium with a grounding line crossing a partly retrograde bed, an outcome which should not occur in laterally unvarying ice streams (Schoof, 2007) but does occur in buttressed systems (Goldberg et al., 2009;Gudmundsson et al., 2012).Those that did not had simply omitted the Ice0 experiment.Every model then exhibited ∼ 1 km a −1 grounding line retreat in response to widespread ablation of the shelf, and nearly every model exhibited a less pronounced response to ablation localized at the calving front.Likewise, nearly every model exhibited grounding line readvance once ablation ceased and ice shelves thickened.One obvious reason for such agreement is the similarity of the models, with most models featuring similar physics and mesh spacing x ∼ 1 km -which most demonstrated to be adequate.This stands in contrast to the previous MISMIP (Pattyn et al., 2012) and MISMIP3d (Pat-  Inasmuch as the MISMIP+ experiments are representative, the choice of basal friction law has more impact on ice stream models than the choice of englacial stress model.Higherorder models deliver essentially the same results as vertically integrated treatments, and there is also little variation between the shallow-shelf/shelfy-stream approximation (SSA) models and vertically integrated (L1Lx) models that include simplified vertical shear streams.The MISMIP+ experiments are, however, rather biased towards this conclusion, with fast sliding evident across much of the domain.This result appears to extend to the full-Stokes models, though we note that only two models of this type were included.On the other hand, models that employ the Coulomb-limited basal friction laws produce faster grounding line retreat and advance compared to models based on the Weertman rule.A similar conclusion in the context of Thwaites Glacier is reached by Yu et al. (2018) The level of agreement between full-Stokes models and other types is less clear.Both participating full-Stokes models exhibited time-averaged grounding line retreat in line with other models in response to ice shelf ablation over 100 a, but one model (ISSM in full-Stokes mode) saw more variation over that period than other models.The other model (Elmer/Ice in full-Stokes mode) saw its grounding line continue to retreat, albeit slowly, for 100 a after ablation ceased, while ISSM saw readvance that was slower than the average across all models but within their range.In other words, where the full-Stokes models disagree with the approximate models, they also disagree with one another.
There is one more notable exception to the rule that basal friction physics matter more than englacial physics.The HySSA submissions showed the same qualitative behavior as but greater rates of grounding line retreat and advance than all of the other model types.This was also observed in the earlier MISMIP3d exercise (Pattyn et al., 2013), but the reasons are not necessarily the same.The HySSA models produced the same steady-state grounding line position as properly resolved SSA models in MISMIP3d with the same value for A, because the boundary layer approximation they rely upon was derived for that exact case, with no buttressing.In this exercise, the HySSA submissions have A around twice as large, presumably because the buttressing formulation is https://doi.org/10.5194/tc-14-2283-2020 The Cryosphere, 14, 2283-2301, 2020 an informal extension to the boundary layer approximation.It may well be the case that the faster retreat is at least in part due to the slacker grounded ice (with the HySSA models setting the rate factor A = 3.5 × 10 −17 Pa −3 a −1 as opposed to the value A = 2.0×10 −17 Pa −3 a −1 typical of other models), but we note that the result might still carry over to realistic problems because it will be equally necessary to tune A or its equivalent when matching observations.However, recent modifications to the HySSA model bring its MISMIP+ grounding line migration rates into the envelope of the results presented here (Pollard and DeConto, 2020).Model initial state -in this case the initial ice thicknessis also a key determinant of the response to ice shelf thinning.Models using the Weertman law had a larger spread of initial states than models based on the Coulomb-limited laws, due to model tuning, and they also showed a larger range of retreat rates when subject to strong melt rates.This is within a group of models that all began in similar states, with similar grounding line shapes crossing a retrograde slope, separated by less than 20 km: a greater difference still would be anticipated if any of the models had started from a more obviously distinct state, for example with the grounding line positioned on an entirely pro-grade slope.
Subgrid treatment of basal friction appears to be of minor importance in these experiments but does offer considerable benefit in other circumstances.Around half of the participating models employ some sort of modification to the discretization of the basal friction term in the stress balance equation, and around half do not, but there is essentially no difference between the results of these two groups.This stands in contrast to the MSIMIP3d experiment, where numerical error associated with the abrupt change in basal friction at the grounding line is a major source of differences between models at mesh resolutions comparable to and finer than the x ∼ 1 km chosen by most MISMIP+ participants (Pattyn et al., 2013).Convergence studies submitted by participants to MISMIP+ (see the model data sheets in the supplement) typically indicate that the choice of x ∼ 1 km is adequate for MISMIP+, subgrid friction scheme notwithstanding, while similar studies based around MISMIP3d reach the opposite conclusion: that a subgrid friction scheme permits x ∼ 1 km and even coarser resolutions while the lack of such a scheme required far finer resolutions (Seroussi et al., 2014a;Feldmann et al., 2014;Gladstone et al., 2010;Leguy et al., 2014).The most obvious difference between the MISMIP3d and MISMIP+ experiments is in the shape of the grounding line and ice shelf, with MISMIP+ having a larger zone of stress transfer between floating and grounded ice, but another possible cause is the colder, stiffer ice and greater basal friction of MISMIP3d.MISMIP3d imposed A ≈ 3.0 × 10 −18 Pa −3 a −1 versus A ≈ 2.0 × 10 −17 Pa −3 a −1 in MISMIP+ and β 2 ≈ 3 × 10 4 Pa m −1/3 a 1/3 as opposed to β 2 ≈ 1 × 10 4 Pa m −1/3 a 1/3 , so that shorter length scales in strain rates are to be expected within the grounding zone.
Subgrid treatment of melt rates -as opposed to subgrid treatment of basal friction -can result in major numerical errors, even at the moderately fine resolutions employed by the majority of participating models.The Ice2 experiments show that imposing melt on even a single cell or element that is partly grounded leads to an erroneous melt rate proportional to the mesh spacing x, which is not small compared to the errors caused by under-resolution of the ice dynamics.In other words, if the mesh is fine enough to obviate the error caused by subgrid melt schemes, it is also fine enough to resolve the dynamics.On the basis of these experiments, then, such treatments should be avoided.They are not to be confused with attempts to treat, for example, tidal variation in the grounding line, which may well result in strong melt rates upstream from the mean grounding line but should not produce a retreat rate that vanishes as x → 0.

Summary
We have presented the results of the third Marine Ice Sheet Model Intercomparison Project, MISMIP+.It is distinct from its predecessors in that its initial state includes a floating ice shelf that buttresses upstream ice, so that when the ice shelf is thinned the grounded ice upstream accelerates and the grounding line retreats, and the converse when the ice shelf is allowed to regrow.More than thirty distinct submissions using 11 model programs produced similar results: an apparently stable equilibrium with a grounding line located on a retrograde slope, rapid grounding line retreat on the ablation of the ice shelf, and slower readvance on ice shelf regrowth.We found that the most important distinctions between models were the basal friction model and the initial state.In contrast, the distinction between models employing several small aspect ratio approximations to the Stokes equations was minor, at least with the limited range of model resolutions tested.The two full-Stokes models themselves do exhibit some distinct behavior, but where they agree with one another they also agree with other models.Another exception was the HySSA model, which makes use of an explicit expression for the flux across the grounding line to permit coarse (∼ 10 km) resolution and thus simulation over longer time intervals.It exhibited substantially faster grounding line migration than any other model, but recent modifications bring it within the range of other models (Pollard and DeConto, 2020).

Figure 1 .
Figure 1.Relationship between the MISMIP+ Ice0, Ice1, and Ice2 experiments.An upward sloping curve indicates an advancing grounding line, a downward sloping curve a retreating grounding line.All three experiments start from a common steady state.The Ice0 experiment provides a control for the Ice1 and Ice2 experiments, which induce retreat with ice shelf melting or calving respectively for 100 a.Later stages of the experiments remove or continue the melting or calving.

Figure 2 .
Figure 2. MISMIP+ domain showing bedrock elevation z b (x, y) and boundary types.The spot height (−644 m) indicates the beginning of a retrograde slope in the center of the channel.
treating basal friction with the Weertman power law model, and two sets of Blatter-Pattyn approximation simulations (HSE_ISSM_HO_Weertman_1km, JCH_ISSM_HO_Tsai_200m), treating basal friction with the Weertman power law model and Coulomblimited basal friction rules respectively.The remaining submissions are all SSA configurations, which use either the Weertman or the Coulomb-limited basal https://doi.org/10.5194/tc-14-2283-2020The Cryosphere, 14, 2283-2301, 2020 S. L. Cornford et al.: MISMIP+ results friction rules and differ in their treatment around the grounding line.

Figure 3 .Figure 4 .
Figure 3. Grounding line contours for all models with x ≤ 1 km at the start of the Ice0 experiments (a) and variation over 100 years (b).In panel (a), yellow contours correspond to SSA, L1Lx, and HO models with the Coulomb-limited basal friction laws; orange contours to SSA, L1Lx, and HO models with the Weertman basal friction law; red (ISSM) and cyan (Elmer/Ice) contours to full-Stokes models; and magenta contours to HySSA models.The color map and black contours depict bedrock elevation.In panel (b), the bold line shows the grounding line position in the center of the channel, x gl (y = 0), mean averaged over all models, and the shaded region depicts the maximum and minimum values.The dashed lines show the range specified in the experimental design, x gl (y = 0) = 450 ± 10 km

Figure 5 .
Figure 5.Rates of change in midchannel grounding line position x gl (y = 0) plotted against rates of change in volume above flotation V f (a) and grounded area A g (b).For all three of the Ice1r, Ice1ra, and Ice1rr experiments each point x gl (y=0) t , V f t lies close to a straight line, as do the points x gl (y=0) t , A g t , despite the considerable variation between models.

Figure 6 .
Figure 6.Midchannel grounding line position plotted against time for the Ice1 experiments.Panel (a) shows the mean (solid curves) and the range (shaded regions) over all models.Panel (b) shows the mean (solid curves) and range (shaded regions) for the main subset, with individual curves (lines and symbols) for the remainder.

Figure 7 .
Figure 7. Ice1 midchannel grounding line position plotted against time for the main subset (a) and the Weertman and Coulomb-limited models (b).Weertman (W) models retreat at two-thirds the rate of Coulomb-limited (S/T) models.The larger spread in the Weertman models is attributed in large part to the spread in initial states.

Figure 8 .
Figure 8. Ice1 midchannel grounding line position plotted against time for the PSU HySSA models compared to other models.The HySSA models (solid lines with symbols) both exhibit a retreat rate around twice as fast as the main subset mean (shaded regions, panel a), and 3 times as fast as the Weertman subset mean (shaded regions, panel b).
Figure 10 compares the Stokes models with their nearest vertically integrated counterpart: the Elmer/Ice full-Stokes and L1Lx models, which both use the Schoof (2005) sliding law, and the ISSM full-Stokes and SSA models, which both use the Weertman (1957) sliding law.The ISSM full-Stokes model exhibits a similar total decrease in grounded area and both a similar regrowth and similar rate of regrowth to its SSA counterpart in Ice1ra.In contrast, the Elmer/Ice full-Stokes model sees its grounded area abate in Ice1r in a similar fashion to its L1Lx counterpart but is the only model that continues to show retreat in the Ice1ra experiment.All four models lose volume above flotation in Ice1r and Ice1ra, but the Elmer/Ice full-Stokes model sees more rapid loss in Ice1ra than the other three.

Figure 9 .
Figure 9. Ice1 midchannel grounding line position plotted for Stokes, higher-order, and other models.The two full-Stokes models (symbols and lines, left panel) see a comparable rate of retreat to the main subset (shaded regions, panel a) in Ice1r but a much lower rate of advance in Ice1ra.Within the main subset, higher-order models (HO, panel b) behave in essentially the same way as the other models (non-HO, panel b).

Figure 10 .
Figure 10.Comparison of full-Stokes and vertically integrated models in the Ice1r (t ≤ 100 a) and Ice1ra (t ≥ 100 a) experiments.Panel (a) shows the change in grounded area A g and panel (b) the change in volume above flotation F f ."FS" indicates a full-Stokes model, "(S)" the Schoof (2005) sliding law, and "(W)" the Weertman (1957) sliding law.

Figure 11 .
Figure 11.Ice1 midchannel grounding line position plotted for models in the Weertman subgroup with varying numerical treatments.(a) Subgroups with (SG) and without (non-SG) subgrid friction schemes show similar rates of retreat and advance.(b) Likewise, the subgroup with x = 1 km is little different from the subgroup with x < 1 km

Figure 13 .
Figure 13.Midchannel grounding line position plotted against time for the Ice2 experiments.Panel (a) shows the mean (solid curves) and the range (shaded regions) over all models.Panel (b) shows the mean (solid curves) and range (shaded regions) for the main subset.

Figure 14 .
Figure 14.Domain edge grounding line position plotted against time for the Ice2 experiments.Panel (a) shows the mean (solid curves) and the range (shaded regions) for the main subset.Panel (b) shows the mean (solid curves) and range (shaded regions) for two subgroups, Ice2A and Ice2B.Models in subgroup Ice2A see little or no grounding line retreat along the domain boundary.
Figure 15.Midchannel grounding line position (a) and grounded area (b) plotted against time for the Ice2A and Ice2B models.The Ice2B models exhibit greater change in grounding line area than Ice2a, but that change is restricted to thin ice outside the channel.

Table 1 .
Parameter values for the MISMIP+ bedrock geometry

Table 2 .
Details of the participating models