Sea-level response to melting of Antarctic ice shelves on multi-centennial time scales with the fast Elementary Thermomechanical Ice Sheet model ( f . ETISh v 1 . 0 )

The magnitude of the Antarctic ice sheet’s contribution to g lobal sea-level rise is dominated by the potential of its marine sectors to become unstable and collapse as a response to ocean (and atmospheric) forcing. This paper presents Ant arctic sea-level response to sudden atmospheric and oceanic forci ngs on multi-centennial time scales with the newly develope d fast Elementary Thermomechanical Ice Sheet (f.ETISh) model. The f.ETISh model is a vertically integrat ed hybrid ice sheet/ice shelf model with an approximate implementation of ice sheet thermomechanics, making the model two-dimensional. Its ma rine 5 boundary is represented by two different flux conditions, co herent with power-law basal sliding and Coulomb basal frict ion. The model has been compared to a series of existing benchmark s. Modelled Antarctic ice sheet response to forcing is dominat ed by sub-ice shelf melt and the sensitivity is highly depend ent on basal conditions at the grounding line. Coulomb friction in the grounding-line transition zone leads to significantl y higher mass loss in both West and East Antarctica on centennial time scal s, leading to 2 m sea level rise after 500 year for a moder ate 10 melt scenario of 20 m a −1 under freely-floating ice shelves, up to 6 m for a 50 m a −1 scenario. The higher sensitivity is attributed to higher driving stresses upstream from the gro unding line. Removing the ice shelves altogether results in a disintegra ion of the West Antarctic ice sheet and (partially) marine b asins in East Antarctica. After 500 years, this leads to a 4.5 m and a 12 .2 m sea level rise for the power-law basal sliding and Coulom b friction conditions at the grounding line, respectively. T he latter value agrees with simulations by DeConto and Polla rd (2016) 15 over a similar period (but with different forcing and includ ing processes of hydrofracturing and cliff failure). The chosen parametrizations make model results largely ind ependent of spatial resolution, so that f.ETISh can potenti ally be integrated in large-scale Earth system models.

Abstract. The magnitude of the Antarctic ice sheet's contribution to global sea-level rise is dominated by the potential of its marine sectors to become unstable and collapse as a response to ocean (and atmospheric) forcing. This paper presents Antarctic sea-level response to sudden atmospheric and oceanic forcings on multi-centennial timescales with the newly developed fast Elementary Thermomechanical Ice Sheet (f.ETISh) model. The f.ETISh model is a vertically integrated hybrid ice sheet-ice shelf model with vertically integrated thermomechanical coupling, making the model twodimensional. Its marine boundary is represented by two different flux conditions, coherent with power-law basal sliding and Coulomb basal friction. The model has been compared to existing benchmarks.
Modelled Antarctic ice sheet response to forcing is dominated by sub-ice shelf melt and the sensitivity is highly dependent on basal conditions at the grounding line. Coulomb friction in the grounding-line transition zone leads to significantly higher mass loss in both West and East Antarctica on centennial timescales, leading to 1.5 m sea-level rise after 500 years for a limited melt scenario of 10 m a −1 under freely floating ice shelves, up to 6 m for a 50 m a −1 scenario. The higher sensitivity is attributed to higher ice fluxes at the grounding line due to vanishing effective pressure.
Removing the ice shelves altogether results in a disintegration of the West Antarctic ice sheet and (partially) marine basins in East Antarctica. After 500 years, this leads to a 5 m and a 16 m sea-level rise for the power-law basal sliding and Coulomb friction conditions at the grounding line, respectively. The latter value agrees with simulations by DeConto and Pollard (2016) over a similar period (but with different forcing and including processes of hydrofracturing and cliff failure).
The chosen parametrizations make model results largely independent of spatial resolution so that f.ETISh can potentially be integrated in large-scale Earth system models.

Introduction
Projecting future sea-level rise (SLR) requires ice sheet models capable of exhibiting complex behaviour at the contact of the ice sheet with the atmosphere, subglacial environment and the ocean. Some of these interactions demonstrate nonlinear behaviour due to feedbacks, leading to self-amplifying ice mass change. For instance, surface mass balance interacts with ice sheets through a powerful melt-elevation feedback, invoking non-linear response as a function of equilibrium line altitude, such as a positive feedback on ablation that can be expected as the ice sheet surface becomes lower (Levermann and Winkelmann, 2016). This feedback is also the main reason for the threshold behaviour of the Greenland ice sheet on multi-millennial timescales (e.g. Ridley et al., 2010). Typical for these self-amplifying effects is that they work both ways: the melt-elevation feedback equally allows for ice sheets to grow rapidly once a given threshold in positive accumulation is reached, resulting in hysteresis (Weertman, 1976).
Another powerful feedback relates to the contact of ice sheets (especially marine ice sheets with substantial parts of the bedrock lying below sea level) with the ocean. F. Pattyn: Antarctic sensitivity to sub-shelf melting Mercer (1978) and Thomas (1979) identified marine ice sheet instability (MISI) for ice sheets where the bedrock dips deeper inland from the grounding line (retrograde bed slopes) so that increased (atmospheric-oceanic) melting leads to recession of the grounding line. This would result in the glacier becoming grounded in deeper water with greater ice thickness. Since ice thickness at the grounding line is a key factor in controlling ice flux across the grounding line, thicker ice grounded in deeper water would result in increased ice discharge and further retreat within a positive feedback loop. Early numerical ice sheet models failed to reproduce this feedback due to the lack of physical complexity (e.g. neutral equilibrium; Hindmarsh, 1993) and the poor spatial resolution to resolve the process of grounding-line migration (Vieli and Payne, 2005;Pattyn et al., 2006). A major breakthrough was provided by an analysis of grounding-line dynamics based on boundary layer theory (Schoof, 2007a(Schoof, , b, 2011, mathematically confirming the earlier findings by Weertman (1974) and Thomas (1979), i.e. that grounding-line positions are unstable on retrograde bedrock slopes in absence of (ice shelf) buttressing. Schoof (2007a) showed that numerical ice sheet models need to evaluate membrane stresses across the grounding line, hence resolving them on a sufficiently fine grid of less than a kilometre, which was further confirmed by two ice sheet model intercomparisons (Pattyn et al., 2012. Since then several marine ice sheet models of the Antarctic ice sheet have been developed, with varying ways of treating the grounding line, i.e. by increasing locally spatial resolution at the grounding line (Favier et al., 2014;Cornford et al., 2015), by making use of local interpolation strategies at the grounding line (Feldmann et al., 2014;Feldmann and Levermann, 2015;Golledge et al., 2015;Winkelmann et al., 2015) or by parametrizing grounding-line flux based on boundary layer theory (Pollard and DeConto, 2009;Pollard et al., 2015;DeConto and Pollard, 2016).
Other feedbacks relate ice sheet dynamics to basal sliding through thermoviscous instabilities, which may lead to limitcycle behaviour in ice sheets (Payne, 1995;Pattyn, 1996) as well as ice stream development in the absence of strong basal topographic control (Payne and Dongelmans, 1997;Payne et al., 2000;Hindmarsh et al., 2009). More elaborate subglacial water flow models have since been developed, exhibiting similar feedback mechanisms in ice discharge (Schoof, 2010). For marine portions of ice sheets, the major subglacial constraint is governed by till deformation and observations have led to new insights in subglacial till deformation based on Coulomb friction controlled by subglacial water pressure (Tulaczyk et al., 2000a, b). In contact with the ocean, subglacial water pressure may therefore stem from the depth of the bed below sea level, which led to new characterizations of grounding-line dynamics (Tsai et al., 2015).
In this paper, I present a new ice sheet model that reduces the three-dimensional nature of ice sheet flow to a two-dimensional problem, while keeping the essential (or elementary) characteristics of ice sheet thermomechanics and ice stream flow. Processes controlling grounding-line motion are adapted in such a way that they can be represented at coarser resolutions. This way, the model can more easily be integrated within computationally demanding Earth system models. A new grounding-line algorithm based on the zero effective pressure conditions reigning at the contact with the ocean has been implemented (Tsai et al., 2015), which leads to a more sensitive grounding-line response. I start by giving a detailed overview of the model and its components. The initialization procedure for the Antarctic ice sheet is then given and, finally, the sensitivity of the Antarctic ice sheet to sudden atmospheric and ocean warming is presented on centennial timescales. The appendices further describe results of known benchmarks for grounded ice flow (Huybrechts et al., 1996;Payne et al., 2000) and marine ice sheet dynamics (Pattyn et al., 2012).

Model description
The model consists of diagnostic equations for ice velocities and three prognostic equations for the temporal evolution of ice thickness, ice temperature and bedrock deformation beneath the ice. Prescribed boundary fields are equilibrium bedrock topography, basal sliding coefficients, geothermal heat flux and sea level. Present-day mean surface air temperatures and precipitation are derived from data assimilation within climate models. Ablation is determined from a positive degree-day (PDD) model. A list of model symbols is provided in Tables 1-3. A general overview of the Cartesian geometry used is given in Fig. 1.
For the coupled ice sheet-ice shelf system the surface elevation h s is defined as where h is the ice thickness, b is the bedrock elevation, z sl is the sea-level height with respect to the chosen datum, ρ i and The ice sheet-ice shelf model has several modes of operation, depending on the boundary conditions that are applied. The most elementary flow regime of the grounded ice sheet is according to the shallow-ice approximation (SIA; Hutter, 1983), extended with either a Weertman type (or power-law) function or a linear/plastic Coulomb friction law for basal sliding. Ice shelf flow is governed by the shallow-shelf approximation (SSA; Morland, 1987;MacAyeal, 1989), defined by zero basal drag and extended by a water-pressure condition at the seaward edge. The transition between both systems is given by a flux condition at the grounding line DeConto, 2009, 2012a), either derived from boundary layer theory based on SSA (SGL; Schoof, 2007a) or given by a flux condition based on Coulomb friction at the grounding line (TGL; Tsai et al., 2015). A second mode of operation is the hybrid mode, in which the flow regime of the grounded ice sheet is governed by a combination of SIA, responsible for ice-deformational flow, and SSA for basal sliding (Bueler and Brown, 2009;Martin et al., 2011;Winkelmann et al., 2011). The hybrid model is used in combination with power-law sliding or linear/plastic Coulomb friction underneath the ice sheet. All components of the flow model are detailed in the sections below.

Shallow-ice approximation
The SIA by Hutter (1983) is commonly used in ice sheet modelling. This approximation is valid for ice sheets of small 1854 F. Pattyn: Antarctic sensitivity to sub-shelf melting aspect ratios h L, where L is the horizontal length scale of the ice sheet domain, and further characterized by a low curvature and low sliding velocities. The approximation is, however, not valid near grounding lines or for ice shelf flow, for which other approximations are applied (see below). According to SIA, the vertical mean horizontal velocity in an ice sheet is given by where τ d = −ρ i gh∇h s is the driving stress, A is the flow parameter in Glen's flow law (with n = 3), v b = (u b , v b ) is the basal sliding velocity and v SIA = (u, v) is the vertical mean horizontal velocity according to SIA. The flow parameter A is a function of ice temperature (see Sect. 2.4).

Hybrid shallow-shelf/shallow-ice approximation (HySSA)
The flow velocity in an ice shelf or an ice stream characterized by low drag is derived from the Stokes equations (Stokes, 1845) by neglecting vertical shear terms and by integrating the force balance over the vertical. The resulting equations are (Morland, 1987;MacAyeal, 1989) 2 and where τ d x = ρ i gh(∂h s /∂x) (similar for τ d y ).ε 0 = 10 −20 is a small factor to keep η finite and hence to prevent singularities when velocity gradients are zero. For the ice shelf, τ b = 0, while for the grounded ice sheet the basal drag is a function of the friction at the base. The SSA stress-equilibrium Eqs. (3) and (4) where n x and n y are the outward-pointing normal vectors in the x and y direction, respectively. The ice shelf velocity field is needed for determining the effect of buttressing in the grounding-line flux conditions (see below) as well as for the thickness evolution of the ice shelf. For the purpose of buttressing, velocity gradients downstream from the grounding line are used to determine the longitudinal stretching rate, which is compared to the stretching rate of a freely floating ice shelf to determine a so-called buttressing factor.
Both SIA and SSA velocities are combined to obtain the velocity field of the grounded ice sheet according to the hybrid model (HySSA; Bueler and Brown, 2009). While Bueler and Brown (2009) use a weighing function to ensure a continuous solution of the velocity from the interior of the ice sheet across the grounding line to the ice shelf, Winkelmann et al. (2011) have demonstrated that a simple addition (for the grounded ice sheet velocities) still guarantees a smooth transition for ice stream flow. Thus basal velocities for the grounded ice sheet are SSA velocities for the velocity field in the grounded ice sheet. In the ice shelf, the SIA velocity is kept zero throughout.

Power-law basal sliding
Basal sliding is introduced as a Weertman sliding law, i.e.
where τ b is the basal shear stress (τ b ∼ τ d for SIA), A b is a basal sliding factor, and m is the basal sliding law exponent. The basal sliding factor A b is temperature dependent and allows for sliding within a basal temperature range between −3 and 0 • C. It further takes into account sub-grid sliding across mountainous terrain (Pollard et al., 2015): where r = max[0, min[1, (T −T r )/(−T r )]], A froz is the sliding coefficient in case of frozen bedrock (chosen to be very small but different from zero to avoid singularities in the basal friction calculation), T is the temperature corrected for the dependence on pressure (see Sect. 2.4.3) and where σ b is the standard deviation of bedrock elevation within the grid cell (Pollard et al., 2015). Basal sliding factors A b are either considered constant in space-time or are spatially varying and obtained through optimization methods (see Sect. 4.1). Basal velocities in the hybrid model are defined through a friction power law, where

Coulomb friction law
Basal friction within the HySSA equations can also be calculated based on a model for plastic till (Tulaczyk et al., 2000a). Several variations of a basal till model can be found in the literature (Schoof, 2006;Gagliardini et al., 2007;Bueler and Brown, 2009;Winkelmann et al., 2011). Deformation of saturated till is well modelled by a plastic (Coulomb friction) or nearly plastic rheology (Truffer et al., 2000;Tulaczyk et al., 2000a;Schoof, 2006). The yield stress τ c satisfies the Mohr-Coulomb relation: where the term between brackets is the effective pressure of the overlying ice on the saturated till (Cuffey and Paterson, 2010) or the ice overburden pressure minus the water pressure p w , c 0 is the till cohesion (c 0 = 0 is further considered) and φ is the till friction angle. The latter can be either taken as a constant value or vary as a function of bedrock elevation (Maris et al., 2014): limited by φ = φ min for b − z sl ≤ −10 3 m and φ = φ max for b − z sl ≥ 0. O b is a spatially varying parameter used to optimize the basal friction field, in a similar way as A b in Eq. (10). Without optimization, it takes the value of O b = 1.
The most comprehensive approach to solve for the subglacial water pressure in Eq. (12) is due to Bueler and van Pelt (2015) by considering a hydrological model of subglacial water drainage within the till. However, Martin et al. (2011) propose to relate major till characteristics to bedrock geometry and allow till friction angle and basal water pressure to be a function of the bed elevation compared to sea level. This leads to zones of weak till and saturation in subglacial basins that are well below sea level Maris et al., 2014). Following their analysis, the subglacial water pressure is defined by Here, λ p is a scaling factor such that the porewater pressure is maximal when the ice is resting on bedrock at or below sea level. Below sea level, the pores in the till are assumed to be saturated with water so λ p is then equal to 1. The factor λ p is scaled with the height above sea level up until 1000 m. At and above 1000 m, λ p is equal to 0 (Maris et al., 2014). While there is no direct physical evidence for such water-pressure distribution in the interior of ice sheets, near grounding lines in direct contact with the ocean the subglacial water pressure of saturated till may also be approximated by (Tsai et al., 2015): which is valid for b − z sl < 0, otherwise p w = 0. By definition, p w = ρ i gh at the grounding line and underneath floating ice shelves so that the effective pressure becomes zero. Bueler and Brown (2009)  To link Coulomb friction to basal drag, the formulation proposed by Bueler and van Pelt (2015) is opted for, where τ c and v b combine to determine τ b through a sliding law, i.e.
where 0 ≤ q ≤ 1, and u 0 is a threshold sliding speed . The sliding law, Eq. (16), includes the case q = 0, leading to the purely plastic (Coulomb) rela- At least in the q 1 cases, the magnitude of the basal shear stress becomes nearly independent of |v b |, when |v b | u 0 . Equation (16) could also be written in a generic power-law form in the linear case q = 1, β 2 = τ c /u 0 (Bueler and van Pelt, 2015).
Alternatively, both the power-law sliding law Eq. (9) and the Coulomb friction law Eq. (16) can be combined (Tsai et al., 2015;Asay-Davis et al., 2016) by taking the lowest friction value of both. Since at the grounding-line basal sliding velocities are considered highest, this equally implies high basal drag in a traditional power-law sliding law. As a consequence, power-law sliding/friction still leads to a relatively sharp transition in τ b at the grounding line (Tsai et al., 2015). Coulomb basal conditions imply that basal drag vanishes towards the grounding line, thus ensuring a smooth transition between the ice stream and ice shelf. Expressing the basal traction as 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 (Durand et al., 2009;Pattyn et al., 2012Durand and Pattyn, 2015). In large-scale models, this can lead to unacceptably small time steps and costly integrations. DeConto (2009, 2012a) incorporated the boundary layer solution of Schoof (2007a) directly in a numerical ice sheet model at coarse grid resolution, so the flux, q g , across model grounding lines is given by This yields the vertically averaged velocity u g = q g /h g , where h g is the ice thickness at the grounding line. in Eq. (18) accounts for back stress at the grounding line due to buttressing by pinning points or lateral shear and is defined as where τ xx is the longitudinal stress just downstream of the grounding line, calculated from the viscosity and strains in a preliminary SSA solution without constraints given by Eq. (18), and τ f the free-water tensile stress defined by b f is an additional buttressing factor to control the buttressing strength of ice shelves and may be varied between 0 (no buttressing) and 1 (full buttressing). All experiments in this paper use b f = 1, except the sensitivity experiments on ice shelf de-buttressing where b f = 0. As in Pollard and DeConto (2012a), C s is Schoof's basal sliding coefficient and m s the basal sliding exponent so that C s is related to the sliding coefficients A b by C s = A b −m s , where m s = 1/m. Groundingline ice thickness h g is linearly interpolated in space by estimating the sub-grid position of the grounding line between the two surrounding floating and grounded h grid points. Therefore, the height above floatation is linearly interpolated on the Arakawa C grid between those two points to where it is zero. Subsequently, the bedrock elevation is linearly interpolated to that location, and the floatation thickness of ice for that bedrock elevation and current sea level is obtained (Pattyn et al., 2006;Gladstone et al., 2010;Pollard and DeConto, 2012a). The velocity u g is then calculated at the groundingline points and imposed as an internal boundary condition for the flow equations, hence overriding the large-scale velocity solution at the grounding line. u g = q g /h g is imposed exactly at the u grid grounding-line point when the flux q g is greater than the large-scale sheet-shelf equation's flux at the grounding line.
Equation (18) applies equally to the y direction, with v g and τ yy instead of u g and τ xx . Note that spatial gradients of quantities parallel to the grounding line, which are not included in Schoof's flow-line derivation of Eq. (18), are neglected here (Katz and Worster, 2010;Gudmundsson et al., 2012;. This parametrization was also found to yield results comparable to SSA models solving transient grounding-line migration at high spatial resolution of the order of hundreds of metres Durand and Pattyn, 2015) despite the fact that Eq. (18) applies to steady-state conditions.

Grounding-line flux condition for Coulomb friction (TGL)
The grounding-line parametrization based on the boundary layer theory by Schoof (2007a) is invalid when Coulomb friction near the grounding line is considered and the effective stress tends to zero. However, Tsai et al. (2015) offers such a solution for vanishing Coulomb friction at the grounding line and therefore independent of basal sliding coefficients: where Q o ≈ 0.61 is a numerical coefficient determined from the boundary layer analysis. The flux in the y direction is obtained in a similar fashion. As in Eq. (18), buttressing scales to the same power as (1−ρ i /ρ w ), which is n−1. The performance of both flux conditions is tested in Appendix C.
The TGL flux condition can be used in conjunction with power-law basal sliding. Indeed, Tsai et al. (2015) have shown that the crossover from Coulomb to power-law roughly occurs at stresses 100 kPa; hence the Coulomb regime occurs within 17 m above the floatation height. This is a very small height difference, which implies that in F. Pattyn: Antarctic sensitivity to sub-shelf melting most cases -with exception of ice plains -a narrow Coulomb regime exists, within a grid cell of a continental-scale model.

Ice thickness evolution
Ice sheet thickness evolution is based on mass conservation, leading to the continuity equation. For the general ice sheetice shelf system, this is written as whereȧ is the surface mass balance (accumulation minus surface ablation) and M is the basal melt rate (solely underneath ice shelves, as basal melt rates underneath the ice sheet are not accounted for). The treatments of the various local ice gains or losses (e.g. surface mass balance) are described in later sections. For the SIA model in the grounded ice sheet, Eq. (22) is written as a diffusion equation for ice thickness (Huybrechts, 1992): where h b is the bottom of the ice sheet (or the bedrock elevation b for the grounded ice sheet). It is also ensured that thinning due to grounding-line retreat does not exceed the maximum permissible rate, using theoretical knowledge of maximum possible stresses at the grounding line that is called the "maximum strain check". Similar to Ritz et al. (2015), tensile stresses at the grounding line are ensured to not exceed those from buttressing by water alone, i.e. the free-water tensile stress, and calculate the maximum corresponding strain rate, expressed as a maximum thinning rate.

Calving and sub-shelf pinning
Ice-front calving is obtained from the large-scale stress field (Pollard and DeConto, 2012a), based on the horizontal divergence of the ice shelf velocities and which is similar to parametrizations used elsewhere Winkelmann et al., 2011;Levermann et al., 2012). The calving rate C r is defined as where w c = min(1, h e /200) is a weight factor and h e is the sub-grid ice thickness within a fraction of the ice edge grid cell that is occupied by ice (Pollard and DeConto, 2012a), defined by where a minimum ice thickness of 30 m avoids too thin ice shelves. The value of h max is defined as the maximum ice thickness of the surrounding grid cells (grounded or floating) that are not adjacent to the ocean (Pollard and DeConto, 2012a). The calving rate C r is then subtracted from the basal melt rate M in Eq. (22). Given the relatively low spatial resolution of a large-scale ice sheet model, small pinning points underneath ice shelves due to small bathymetric rises scraping the bottom of the ice and exerting an extra back pressure on the ice shelf Favier et al., 2016) are not taken into account. To overcome this a simple parametrization based on the standard deviation of observed bathymetry within each model cell was accounted for to introduce a given amount of basal friction of the ice shelf (Pollard and DeConto, 2012a). The fractional area f g of ice in contact with sub-grid bathymetric high is defined as (modified from Pollard and DeConto, 2012a) where h w is the thickness of the water column underneath the ice shelf and σ b is the standard deviation of the bedrock variability (see above). This factor f g is multiplied with β 2 in the basal friction. For the grounded ice sheet, f g = 1; for the floating ice shelf in deeper waters, f g = 0, so that the ice shelf does not experience any friction.

Ice sheet temperature
The diffusion-advection equation for an ice sheet is given by (Huybrechts, 1992) where κ = K/ρ i c p is the thermal diffusivity of ice, K is the thermal conductivity, c p is the heat capacity of ice, θ is the ice temperature, and The basal boundary condition is given by where G is the geothermal heat flux and the second term represents frictional heating at the base. The last term in Eq. (28) represents strain heating. Given the two-dimensional nature of the model, the temperature field employs shape functions for vertical profiles of deformational velocity v d , its vertical gradient, and the vertical velocity, based on SIA (Hindmarsh, 1999). Equation (27) is then solved in scaled vertical coordinates ζ = (h s − z)/ h, with ζ = 0 at the surface and ζ = 1 at the bottom of the ice sheet. The use of shape function allows for a faster calculation of the thermodynamic model. However, since this is an approximation compared to fully solv-ing Eq. (27), the EISMINT-I benchmark experiments (Huybrechts et al., 1996) were performed and results are given in Appendix A.

Ice shelf temperature
In ice shelves, a simple temperature model is adopted, considering the accumulation at the surface balanced by basal melting underneath an ice shelf and with only vertical diffusion and advection in play (Holland and Jenkins, 1999): where β 1 =ȧζ h/κ, β 2 =ȧh/κ and θ s b is the ocean temperature at the base of the ice shelf, corrected for ice shelf depth, i.e. θ s b = T oc = −1.7 − 0.12 × 10 −3 h b (Maris et al., 2014).

Thermomechanical coupling
The mean column temperature T is obtained by integrating θ from the base of the ice sheet to a given height in the ice column. Since most of the ice deformation is in the bottom layers of the ice sheet, the temperature closest to the bottom determines to a large extent the deformational properties. Compared to full thermomechanically coupled ice sheet models, satisfactory results where obtained by considering a mean column temperature for the lower most 10-40 % of the ice column. This fraction can also be regarded as an extra tuning parameter in an ensemble run, especially given the large uncertainties pertaining to geothermal heat flow underneath major ice sheets. The flow parameter A and its temperature dependence on temperature are specified as in Huybrechts (1992) and Pollard and DeConto (2012a): where T = T − T m is the homologous temperature, with T m = −8.66 × 10 −4 (1 − ζ )h the pressure melting correction and R the gas constant. Units of A are Pa −3 yr −1 corresponding to n = 3. The enhancement factor E f is set to 1 for the main ice sheet model, and to E f = 0.5 for ice shelves. The ratio of enhancement factors represent differences in fabric anisotropy between grounded and ice shelf ice (Ma et al., 2010). Verification of the thermomechanical coupling scheme using a vertical mean value of A follows the EISMINT-II benchmark experiments (Payne et al., 2000) and is detailed in Appendix B.

Bedrock deformation
The response of the bedrock to changing ice and ocean loads is solved through a combined time-lagged asthenospheric relaxation and elastic lithospheric response due to the applied load (Huybrechts and de Wolde, 1999;Pollard and DeConto, 2012a). The deflection of the lithosphere is given by where D is the flexural rigidity of the lithosphere and ρ b is the bedrock density. The load is then defined by where h w is the ocean column thickness and h eq and h eq w are the values of ice thickness and ocean column thickness in equilibrium, respectively, taken from modern observed fields. Equation (32) is solved by Green's function (Huybrechts and de Wolde, 1999). The response to a point load P w (q b × area) versus distance from the point load l is then given by where kei is a Kelvin function of zeroth order (defined as the imaginary part of a modified Bessel function of the second kind) and L w = (D/ρ b g) 1/4 ≈ 132 km is the flexural length scale. For any load, the different values of the point loads w p are summed over all grid cells to yield w b (x, y). Finally, the actual rate of change in bedrock elevation is given by a simple relaxation scheme: where b is the actual bedrock elevation, b eq is the elevation in equilibrium (taken from modern observed fields) and τ w = 3000 years (Pollard and DeConto, 2012a).

Numerical grid and solution
The ice sheet-shelf model uses a finite-difference staggered grid, where horizontal velocities (u, v) are calculated on two separate staggered Arakawa C grids, as is usual for vector fields (Rommelaere and Ritz, 1996), while diffusion coefficients for the ice sheet equation d are calculated on an Arakawa B grid, staggered in both the x and y direction since these are scalar quantities (Fig. 2)

once in one matrix
where N x and N y are the number of grid points in the x and y direction, respectively. The submatrices A ux and A vx contain the coefficients for the solution in the x direction for u and v, respectively. A uy and A vy are defined in a similar way. Due to the dependence of the effective viscosity η on u, v, the solution requires a few iterations to reach convergence. A similar solution approach is taken for solving the continuity equation for ice thickness (Payne and Dongelmans, 1997), which was favoured over an Alternating Direct Implicit scheme used in several ice sheet models (Huybrechts, 1992;Pollard and DeConto, 2012a). The f.ETISh model is implemented in MATLAB ® . Computational improvements involved the omission of all "for" loops by using circular shifts (with exception of the time loop), thereby optimizing the use of matrix operations. The bulk of computational time is devoted to the solution of the sparse matrix systems, which are natively optimized in MATLAB ® using multi-threading. A preconditioned conjugate gradient method is used for solving the ice sheet-ice shelf continuity equation. The velocity field in the hybrid model is solved using a stabilized bi-conjugate gradients method, which is also preconditioned and further initialized by the velocity field solution from the previous time step. Both numerical solvers are iterative and the preconditioning limits the number of iterations to reach convergence. They are considerably faster compared to the direct solution.
The f.ETISh model is compared to other ice sheet models via a series of benchmarks, such as the EISMINT-I benchmark for isothermal ice sheet models (Huybrechts et al., 1996, Appendix A), the EISMINT-II benchmark for thermomechanically coupled ice sheet models (Payne et al., 2000, Appendix B), and the MISMIP experiments for marine ice sheet models (Pattyn et al., 2012, Appendix C). Results show that the f.ETISh model is in close agreement with all of the benchmark experiments.
3 Input and climate forcing

Input datasets
For modelling the Antarctic ice sheet, the bedrock topography is based on the Bedmap2 data (Fretwell et al., 2013), from which ice thickness, present-day surface topography and grounding-line position are derived. Surface mass balance and temperatures are obtained from Van Wessem et al. (2014), based on the output of the regional atmospheric climate model RACMO2 for the period 1979-2011 and evaluated using 3234 in situ mass balance observations and icebalance velocities.
For geothermal heat flux we employ a recent update of Fox-Maule et al. (2005) by Purucker (2013). It is based on low-resolution magnetic observations acquired by the CHAMP satellite between 2000 and 2010 and produced from the MF-6 model following the same technique as described in Fox-Maule et al. (2005).
All datasets are resampled on the spatial resolution used for the experiments. The experiments shown in this paper employ a grid spacing of 25 (and in a few cases 40 or 16) km.

Atmospheric and ocean forcing
Atmospheric forcing is applied in a parametrized way, based on the observed fields of precipitation (accumulation rate) and surface temperature. For a change in background (forcing) temperature T , corresponding fields of precipitation P and atmospheric temperature T s are defined by (Huybrechts et al., 1998;Pollard and DeConto, 2012a) where γ = 0.008 • C m −1 is the lapse rate and δT is 10 • C (Pollard and DeConto, 2012a). The subscript "obs" refers to the present-day observed value. Any forcing (increase) in background then leads to an overall increase in surface temperature corrected for elevation changes according to the environmental lapse rate γ . The parametrizations of T s and P can easily be replaced by values that stem from global circulation models, with appropriate corrections for surface elevation (e.g. de Boer et al., 2015).
Surface melt is parametrized using a PDD model (Huybrechts and de Wolde, 1999). The total amount of PDDs is obtained as where σ is taken as 5 • C (Reeh, 1989) and T is the mean annual temperature. The annual number of PDDs represents a melt potential, used to melt snow and (superimposed) ice. This is determined by applying a seasonal cycle to the atmospheric temperatures with a double amplitude of 20 • C, linearly increasing to 30 • C at an elevation of 3000 m, and kept at 30 • C at higher elevations (Pollard and DeConto, 2012a). The PDD melt potential is related to surface melt through a coefficient of 0.005 m of melt per degree day (Pollard and DeConto, 2012a). Although more complex schemes are often used, taking into account refreezing of percolating meltwater in the snow pack and melting of superimposed ice with different melt coefficients (Huybrechts and de Wolde, 1999), which is also confirmed by recent observations (Machguth et al., 2016), surface melt is rather limited for the present-day Antarctic ice sheet. Surface mass balance is then the sum of the different components, i.e.ȧ = P − S, where S = 0.005× PDD is the surface melt rate. Melting underneath the floating ice shelves is often based on parametrizations that relate sub-shelf melting to ocean temperature and ice shelf depth (Beckmann and Goosse, 2003;Holland et al., 2008), either in a linear or a quadratic way Pollard and DeConto, 2012a;de Boer et al., 2015;DeConto and Pollard, 2016). This leads to higher melt rates close to the grounding line, as the ice shelf bottom is the lowest. While the adaptation by Holland et al. (2008) and Pollard and DeConto (2012a) is implemented in f.ETISh, only constant values of ice shelf melt were used in the experiments for this paper, scaled by a melt factor F melt . This factor distinguishes protected ice shelves (Ross and Ronne-Filchner; Fig. 3), with a melt scaling factor of F melt = 0.125, from all other ice shelves that have a scaling factor of F melt = 1. A similar approach has been taken by many other ice sheet models cited in de Boer et al. (2015).
4 Present-day Antarctic ice sheet simulation

Initialization
Model initialization to the modern Antarctic ice sheet geometry is based on the method by Pollard and DeConto (2012b) by optimizing basal sliding coefficients in an iterative fashion. This nudging scheme is applied to both the Weertman type power law and the Coulomb friction law so that it can be used in conjunction with the two types of grounding-line flux conditions. The model (with grounding lines and floating ice constrained as described above) is run forward in time, starting from modern observed bed and ice surface elevations and driven by the observed climatology (surface mass balance and temperature). Full thermomechanical coupling and temperature evolution, isostatic bedrock adjustment, calving and sub-grid ice shelf pinning are equally considered. For the Weertman sliding law, basal sliding coefficients A b (x, y) are initialized with a constant value (A b = 3 × 10 −9 m a −1 Pa −2 ) for the grounded ice sheet and a higher value (A b = 10 −5 m a −1 Pa −2 ) underneath ice shelves and the ocean to account for slippery saturated marine sediments in case of re-grounding. At intervals of t inv years, at each grounded ice grid point, the local basal sliding coefficients A b (x, y) in Eq. (9) are adjusted by a multiplicative factor (Pollard and DeConto, 2012b): where z = max −1.5, min 1.5, is the observed ice surface elevation and h inv s is a scaling constant. During the inversion procedure, basal temperature is still allowed to influence sliding. Adjusted A b (x, y) values are also not allowed to exceed 10 −5 m a −1 Pa −2 , representing the slipperiest deformable sediment. At the grounding line, observed surface velocities  are used to define the buttressing factors at the grounding line in the grounding-line flux condition. Values for A b are only updated when r > 0 in Eq. (10), so that they are kept unchanged when ice is frozen to the bedrock.
In addition to Pollard and DeConto (2012b) we also introduce a regularization term that essentially smooths highfrequency noise in the basal sliding coefficients by using a Savitzky-Golay filter of degree 3, with a span of 160 km (surrounding influence matrix). The advantage of such a filter is that it keeps lower-frequency variability intact while removing high-frequency noise. The filter is only applied for marine areas (b − z sl < 0) as it improves the fit in these areas compared to the non-regularized case and guarantees a smooth transition between the inland bed and the more slippery ocean beds under present-day ice shelves.
For the Coulomb friction law, optimization starts with a constant field of O b = 1. Equation (40) then transforms to Values of O b are limited between 0.01 and 5 in order to keep tan φ between physically plausible values.
Optimized basal sliding coefficients (Fig. 4) for the Antarctic ice sheet on a spatial resolution of 25 km were obtained after a forward integration of 80 000 years with h inv = 2000 and t inv = 1000 year. This results in a small difference (within 100 m) between the observed and the steadystate modelled topographic surface of the interior ice sheet (Fig. 4). The highest sliding coefficients are found in the marginal areas, especially in the Siple Coast sector, as well as under Pine Island (PIG) and Thwaites (TWG) glaciers. Higher values are also encountered in the centre of the ice sheet, which is also obvious in other studies (Pollard and DeConto, 2012b;Bernales et al., 2017). These areas also show larger misfits (Fig. 4) and may be attributed to the poor knowledge of bedrock topography and so uncertainties are translated into a basal friction anomaly. The obtained patterns are in general agreement with the results from Pollard and DeConto (2012a, b); i.e. the largest errors are found around the major mountain ranges (e.g. Transantarctic Mountains), since outlet glaciers protruding through these mountain ranges are not well represented on coarser grid cells. However, this fit has been improved by including bedrock variability in determining basal sliding coefficients A b in Eq. (10) to allow for basal sliding of smaller outlet glaciers across mountain ranges.
The lower row of Fig. 4 displays the result for the Coulomb friction law, in combination with the grounding-line flux condition of Tsai et al. (2015). The pattern of optimized friction parameters is similar to the one obtained for Weertman sliding (but inverse, since it displays friction instead of sliding). The optimization results in a slightly larger misfit (especially near the Wilkes Basin in East Antarctica), and this may be attributed to the rather coarse approach taken here to account for the spatial distribution of subglacial water pressure and till friction angle.
The basal temperature fields (Fig. 4) for both optimizations are quite similar and in general agreement with basal temperature fields from other Antarctic modelling studies. Differences can easily be attributed to the use of geothermal heat flow datasets, which has the largest impact on basal temperature distribution (Pattyn, 2010).

Model validation
Modelled velocities form an independent check of the model performance, since the optimized basal sliding coefficients are obtained solely from the observed surface topography. The modelled flow field of the Antarctic ice sheet (Fig. 5) compares well to observations of surface velocities due to , such as the delineation of the different drainage basins and major ice streams discharging into the ice shelves. Some disagreement is found on glaciers discharging through the Transantarctic Mountains in the Ross ice shelf as well as glaciers near the Ellsworth Mountains discharging in the Ronne ice shelf. Those mismatches can be traced back to the difficulty in resolving those feature during the initialization process.
A direct comparison between the present-day velocity field  and modelled velocities is shown in Fig. 5. The scatterplot shows a qualitatively good one-to-one fit for both the grounded ice sheet and the floating ice shelves. Quantitative error analysis shows a mean misfit of 11 m a −1 with a standard deviation of 190 m a −1 for the grounded ice flow and a mean misfit of 97 m a −1 with a standard deviation of 1572 m a −1 for the floating ice shelves. The histogram comparison demonstrates a good overall fit of observed and modelled velocity magnitudes and the result is in line with other model studies (e.g. Martin et al., 2011).

Sensitivity to ice shelf de-buttressing
Ice shelves are the prime gatekeepers of Antarctic continental ice discharge. The breakup of the Larsen B ice shelf (Fig. 3) and the subsequent speed-up of outlet glaciers that previously discharged into the ice shelf witness this important instability mechanism (Scambos et al., 2000(Scambos et al., , 2004. In West Antarctica, observational evidence  as well as modelling studies (Favier et al., 2014;Joughin et al., 2014;Seroussi et al., 2014) show that the reduction in buttressing of ice shelves in the Amundsen Sea embayment may lead to significant inland ice mass loss and that unstoppable retreat of the grounding line of Thwaites Glacier may already be on its way . Since ice shelf buttressing is a key element in the stability of the Antarctic ice sheet, a useful experiment to understand underlying model buttressing physics is the sudden removal of all floating ice shelves, starting from the initialized model state, and to let the model evolve over time. Over this period ice shelves were not allowed to regrow, which is equivalent to removing all floating ice at each time step. This experiment is carried out for three cases: (i) power-law sliding with the flux condition according to Schoof (2007a) (SGL), (ii) Coulomb friction with flux condition according to Tsai et al. (2015) (TGL) and (iii) power-law sliding with the TGL condition (TGL-1). All experiments result in a sudden ice-mass loss and grounding-line retreat, whereby the West Antarctic ice sheet collapses entirely in less than 200 years according to SGL and less than 100 years according to TGL (Fig. 6). Both TGL experiments lead to a similar mass loss (both in terms of timing and volume). Therefore, the decisive factor governing mass changes is the grounding-line flux condition and not the sliding/friction law that is employed for the grounded ice sheet.
For all experiments, grounding-line retreat starts in the marine sections discharging in the Ronne and Ross ice shelves. For the SGL experiment, the retreat from Ellsworth Land leads to thinning in the inland sectors of the Pine Island basin, which after > 50 years triggers grounding-line retreat from Pine Island Glacier and subsequently Thwaites Glacier.
Grounding-line retreat then spreads rapidly towards the Ross sector of the West Antarctic ice sheet, leading to a complete disintegration of the ice sheet within 150 years. However, for both TGL experiments, initial grounding-line retreat also occurs in the Amundsen Sea sector, whereby the retreat is much faster and the ice sheet collapses within less than 100 years. Another major difference between SGL and TGL experiments is that the total SLR contribution for TGL is three times as large compared to SGL, i.e. ∼ 16 m for TGL compared to ∼ 5 m for SGL after 500 years. The extra mass loss is essentially located in the East Antarctic ice sheet, i.e. Wilkes and Aurora basins (Wilkes Land; Fig. 3), both losing substantial amounts of ice. Despite the presence of a sill at the outlet of Wilkes subglacial basin, grounding-line retreat occurs without invoking any other physical mechanism than the flux condition at the grounding line in combination with complete ice shelf collapse. These results contrast with Mengel and Levermann (2014) who require the removal of a specific coastal ice volume equivalent to 80 mm of SLR in order to provoke an unstable grounding-line retreat within Wilkes basin.
The higher TGL grounding-line sensitivity must be sought in its underlying physics: at the grounding line the basal shear stress vanishes in a smooth way to reach zero exactly at the grounding line. As shown by Tsai et al. (2015), this is not the case for the SGL algorithm, where a sharp contrast between the inland non-zero basal shear stress and the ocean exists. This boundary becomes smoother with larger sliding velocities, leading to a larger transition zone (Pattyn et al., 2006;Gladstone et al., 2012;Feldmann et al., 2014), but the transition jump does not vanish. The SGL condition at the grounding line is therefore a function of the friction coefficient A b , while the TGL condition is related to a single parameter in the friction law, i.e. the till friction angle. The latter is also limited in its range, contrary to A b ranging across several orders of magnitude (from saturated till to nearly frozen bedrock). Furthermore, the TGL condition is a function of ice thickness h to a higher power compared to SGL. Since the TGL ice flux is larger than the SGL flux for similar conditions, the surface gradient at the grounding line is generally higher, hence leading to higher driving stresses. These steeper surface slopes make the grounding line to retreat (and advance) more rapidly than with the power-law condition (SGL). The higher sensitivity for TGL is also demonstrated in the modified MISMIP experiments (Appendix C). Additionally, I carried out a series of sensitivity tests by fixing the value of the till friction angle at the grounding line φ, ranging from 10 to 60 • . Only for φ ≥ 50 • did the sensitivity decrease, but the amount of mass loss was still significantly higher than with the SGL condition.

Sensitivity to sub-shelf melt
Antarctic ice sheet sensitivity to sub-shelf melting is investigated with a multi-parameter/multi-resolution forcing ensemble over a period of 500 years. Atmospheric forcing includes changes in background temperature T , ranging from 0 to +8.5 • C, affecting both surface temperature, Eq. (37), and surface mass balance, Eq. (38), through the mass balance-elevation feedback. Surface melt is calculated with the PDD model (Eq. 39). Ocean forcing is based on constant forcing values of sub-shelf melting M, ranging from 0 to 50 m a −1 underneath the freely floating ice shelves surrounding the Antarctic ice sheet and between 0 and 6.25 m a −1 for the Ronne-Filchner and Ross ice shelves (factor 8 less compared to the freely floating ice shelves). Melting is only applied to fully floating grid cells, without taking into account the fractional area of grounded grid points that are actually afloat, as done in a few studies (Feldmann et al., 2014;Golledge et al., 2015). All forcings are applied as a sudden change in temperature/melt rate starting from the initialized model. A background run (without applying the forcing anomaly) is also performed to determine the model drift on the different timescales. in background temperature/basal melting rate underneath the ice shelves on a grid size of = 25 km (as well as on a = 40 km grid to test grid-size dependence). A few runs are performed on a grid size of = 16 km for comparison.
Sea-level contribution according to the forcing experiments and rate of change of sea level for the = 25, 40 km spatial resolutions are shown in Fig. 7. These are determined from the change in ice volume above floatation . SLR according to the forcings ranges between −0.5 and 7 m after 500 years. It increases with increasing sub-shelf melt rates and slightly decreases with increasing atmospheric temperature forcing. The latter is due to higher precipitation rates in a warmer climate, leading to an increase in grounded ice mass. However, for larger atmospheric forcing (+8.5 • C), mass loss is generally enhanced due to the dominance of surface melt and/or increase in ice flux with increased precipitation rates. The different curves in Fig. 7 are clustered according to sub-shelf melt rate, which is the most decisive process governing mass loss. Atmospheric forcing, however, has only a limited effect, probably because the timescale considered (500 years) is too short to relax the ice sheet to the imposed temperature and precipitation changes and because weakening of ice shelves through hydrofracturing is not taken into account. Model drift (zero forcing anomaly) is between 60 and 75 cm of sea level lowering over a period of 500 years, or 0.2-0.3 % of the total Antarctic ice sheet volume per century. This is comparable to other Antarctic model studies (e.g.  and shows that the initialization is rather stable and close to steady state. The major differences in sea-level response are due to the treatment of grounding-line fluxes. As shown above, the TGL flux condition systematically leads to significant higher mass losses, making grounding-line migration a more sensitive process (Sect. 5.1). The higher sensitivity leads to a rate of change in sea level of up to 30 mm a −1 . These high values correspond to periods of MISI. Note, however, that such rates are still significantly lower than those obtained during the ice shelf removal experiment. For the SGL flux condition, these values are half as much, and major MISIs occur generally at a later stage during the model run. Compared to other studies Ritz et al., 2015;DeConto and Pollard, 2016), the TGL flux conditions puts sea-level contributions at the high end of the spectrum.
Only the higher melt-rate scenarios (> 10 m a −1 ) produce significant MISIs over this time period. They first occur in the West Antarctic ice sheet (WAIS), starting from either Pine Island or Thwaites Glacier, progressing inland. Other MISI-prone areas are the Bellingshausen Sea (WAIS) and Wilkes basin (East Antarctic ice sheet -EAIS). Contrary to the de-buttressing experiment in Sect. 5.1, MISIs are not initially triggered in the Siple Coast area nor through Ellsworth Land. This is probably due to the lower imposed melt rates, so that both Ronne and Ross ice shelves remain buttressed for a longer period of time.
The effect of spatial resolution on model result is summarized in Fig. 8 in addition to the data presented in Fig. 7. Coarser resolutions (40 km) give comparable results to the 25 km grid, especially for zero melt forcing and the highest melt forcing according to Coulomb TGL. Those cases correspond to either absence of MISIs (low SLR) or complete disintegration of WAIS (high SLR). The main reason for this relatively good fit must be sought in the grounding-line flux conditions (SGL and TGL) that make the model resolution independent. Models that are not based on such heuristics have to resolve grounding-line migration at sub-kilometre resolutions . Differences in response (medium scenarios) are due to the precise timing of MISIs, which seems to be resolution dependent; some of the MISIs are not completed after 500 years (Fig. 7). However, a spatial resolution of 40 km generally remains to coarse, and results are much improved at 25 km. This is demonstrated by the comparison of 16 to 25 km resolution for which obtained SLR is almost the same (crosses in Fig. 8), even for the medium scenarios. Nevertheless, it is expected that at very high spatial resolutions (< 5 km), grounding-line retreat is influenced by bedrock irregularities as well as the presence of ice shelf pinning points that are not always properly resolved at coarser resolutions. The parametrization of sub-grid processes, such as basal sliding in mountainous areas and sub-shelf pinning at sub-grid level, have to some extent reduced this dependency in the model, but differences remain.
In order to validate this claim, two more experiments were carried out to make comparison with an existing experimental result at high resolution possible . Here, sub-shelf melting is taken as a function of ice thickness , i.e. 1 M = max min 4 7 (H − 100), 400 , 0 .
It limits the melt rate between zero (for ice shelves thinner than 100 m) and 400 m a −1 (for ice shelves thicker than 800 m). Results are shown in Fig. 9. The total contribution to SLR after 500 years in the SGL experiment (3.9 m) is comparable to the finest mesh experiment in Cornford et al. (2016). As expected, the TGL experiment gives a much higher mass loss due to its inherent physics. Differences between the model response are sought in the timing of grounding-line retreat within particular drainage basins. For instance, the grounding line in the SGL experiment starts to retreat in the Siple Coast, Ellsworth Land and PIG (as in Cornford et al., 2016), while TWG kicks in at a later time. However, for the TGL experiment, both PIG and TWG retreat at about the same time at the start of the model run. Such differences in response are to be expected, since both experiments are run on a much coarser resolution (25 km) than in Cornford et al. (2016) and hence have a different basal topography.

Discussion
In terms of model complexity, the f.ETISh model is comparable to the Pollard and DeConto (2012a) model. The major difference lies in a number of simplifications that makes the f.ETISh model two dimensional. This is obtained by approximating the temperature coupling by relating the mean icecolumn temperature to the velocity field via the commonly used Arrhenius relationship (Cuffey and Paterson, 2010). Another major difference pertains to the marine boundary, with the implementation of the grounding-line flux condition according to Tsai et al. (2015), based on a Coulomb friction law (TGL), further extended with a Coulomb friction law for the interior ice sheet. Finally, model initialization based on Pollard and DeConto (2012b) has been further extended with a regularization term that essentially smooths the basal friction field across marine basins and makes the results independent of spatial resolution, since regularization is made a function of horizontal distance instead of number of grid cells. Moreover, the optimization does not involve an optimization of ice shelf basal mass balance, since observed ice shelf velocities are used to determine the amount of buttressing at the grounding line. The resulting initialization is characterized by a small drift once the grounding line is allowed to relax, of the order of 0.2-0.3 % of the ice sheet volume in 100 years. Other marine elements such as hydrofracture and cliff failure (Pollard et al., 2015;DeConto and Pollard, 2016) are not taken into account.
Given the differences in approach with continental-scale ice sheet models, such as AISM-VUB (Huybrechts, 1990(Huybrechts, , 2002, ANICE (de Boer et al., 2013), GRISLI (Ritz et al., 2015), ISSM (Larour et al., 2012), PISM (Bueler and Brown, 2009), PISM-PIK Winkelmann et al., 2011;Golledge et al., 2015), PSU-ISM (Pollard and De-Conto, 2012a), RIMBAY (Thoma et al., 2014) or SICOPO-LIS (Sato and Greve, 2012), verification of the f.ETISh model requires a detailed comparison with existing benchmarks. These are generally based on results of the models cited above. The EISMINT-I benchmark (Huybrechts et al., 1996) shows that the ice-dynamical characteristics of f.ETISh are in very close agreement with the benchmarks shown in Appendix A despite a different numerical solution scheme. The basal temperature field is also in close agreement. The results of thermomechanical coupling of ice sheet flow are also in good agreement with the EISMINT-II benchmark (Payne et al., 2000), although the range of uncertainty between the different participating models on which the benchmark is based is also much larger.
An important experiment for marine ice sheet models is a test of steady-state grounding-line positions in absence of buttressing (Pattyn et al., 2012). Boundary layer theory indeed predicts that unique grounding-line positions exist on a downward-sloping bed, while no stable solutions are found on reversed bed slopes (Schoof, 2007a), unless buttressing is significant . While the experiments are designed for flow-line models, they can be extended to two dimensions to evaluate the behaviour in a qualitative way. Here, the f.ETISh model successfully passes the test independent of model resolution, as grounding-line migration is governed through a heuristic based on the abovementioned boundary layer theory DeConto, 2009, 2012a) and is extended with a heuristic based on Tsai et al. (2015) that qualitatively gives the same results.
The main advantage of using a grounding-line flux parametrization based on a heuristic rule (Sect. 2.1.6) is that the model can be run at lower spatial resolutions, which is be resolved with sufficient detail (Schoof, 2007a), which requires the use of sub-kilometre grid sizes (Pattyn et al., 2012), unless sub-grid grounding-line parametrizations are used that may allow for larger grid sizes (Feldmann et al., 2014;Cornford et al., 2016). The main disadvantage of the heuristic rule is that its parametrization is derived from a steady-state solution based on the SSA model. It can therefore be questioned whether the formulation still holds for transients. It also overrules the hybrid model at this particular location. A major finding in this paper is the increased sensitivity of the grounding line based on a Coulomb friction law (Tsai et al., 2015), compared to a power-law sliding condition at the grounding line. Power-law sliding mechanisms near grounding lines have been extensively discussed, since they lead to sudden jumps in basal drag at the grounding line, especially at relatively low sliding speeds (such as in the MISMIP and MISMIP3d experiments; Pattyn et al., 2012. However, sliding velocities in the Antarctic experiments are not preconditioned by a specific sliding coefficient at the grounding line but determined from the optimization procedure. Therefore, the type of boundary is controlled by the model physics itself. The Coulomb friction condition at the grounding line is consistent with observations, as the ice sheet profiles "taper off" towards a flattening upper surface, contrary to the power-law case, and basal stresses vanish at the grounding line (Tsai et al., 2015). Moreover, the grounding-line ice flux according to Coulomb friction also depends more strongly on floatation ice thickness, implying higher sensitivity to atmospheric and ocean forcing. Furthermore, grounding is facilitated in shallower water compared to the power-law case, so that smaller perturbations may push the grounding line more easily into regions with a retrograde slope, provoking a grounding-line instability (Tsai et al., 2015). As a result of the higher sensitivity, Antarctic sea-level contribution to a given perturbation is also more than twice as high and rates of sea-level change three times as fast compared to a power-law sliding case.
Direct comparison is not possible with recent studies of Antarctic ice mass loss that are forced by atmosphere-ocean models following so-called RCPs (Representative Concentration Pathways). Direct comparison with the SeaRISE experiments  is also hampered due to the lower melt rates applied to the Ross and Ronne-Filchner ice shelves. This differentiation was deliberately chosen, as the de-buttressing experiments show that the highest buttressing stems from those large ice shelves. However, their grounding lines are also farthest from the continental shelf break, hampering the intrusion of warmer waters compared to the smaller ice shelves that are closer to the edge. However, considering the f.ETISh model with the SGL condition comparable to the PSU-ISM model DeConto, 2009, 2012a), some comparison on sensitivity can be made. For the SeaRISE experiments, the PSU-ISM model predicts a sea-level contribution after 500 years according to a 2 × A1B scenario (with-out sub-shelf melting) of ∼ 0.45 m, while the f.ETISh SGL model results in ∼ 0.4 m for similar forcing conditions. One has to note, however, that the initialization of both models is different (spin-up versus optimization).
However, the TGL model is less sensitive than the PSU-ISM model including cliff failure and hydrofracturing (De-Conto and Pollard, 2016). These processes potentially lead to a sea-level contribution of 12-13 m after 500 years under a RCP8.5 scenario forced by atmosphere-ocean models. This result corresponds well with the results of the f.ETISh TGL model under complete de-buttressing (without ice shelf growth), with complete collapse of the West Antarctic ice sheet and major ice loss in the Wilkes and Aurora basins (Fig. 6).
Finally, computational time of f.ETISh largely depends on the spatial resolution, which also governs time steps needed under the CFL condition. A hybrid model 5000-year run with a grid size of 40 km and a time step of 0.2 years takes approximately 10 000 CPU seconds on a single AMD Opteron 2378 2.4 GHz core of the Hydra cluster (VUB-ULB) and 20 000 CPU seconds for a 500-year run with a grid size of 16 km and time step of 0.02 years on a multicore. Future developments will focus on improving the numerical solution schemes in order to reduce the calculation time (larger time steps), especially at higher spatial resolutions.

Conclusions
I developed a new marine ice sheet model, based on common descriptions of ice physics (combined shallow-ice and shallow-shelf approximation) and novel implementation of parametrizations of thermodynamics and grounding-line migration. The model has been extensively tested against existing benchmarks and has been shown to be scale independent, with the exception of grounding zones with small-scale bedrock variability, where grounding-line response to atmospheric and oceanic forcing is sensitive to spatial resolution. This makes the model extremely attractive to couple within Earth system models.
The model has been initialized to the present-day Antarctic ice sheet conditions in order to obtain initial steady-state conditions as close as possible to the observed ice sheet. Independent validation has been obtained through comparison with observed surface velocities that are not utilised during the optimization phase.
Two forcing experiments over a period of 500 years are carried out, one during which all floating ice shelves are removed and one during which sudden atmospheric and oceanic forcing is applied. Both experiments show a very high sensitivity to grounding-line conditions, as Coulomb friction in the grounding-line transition zone leads to significantly higher mass loss in both West and East Antarctica, compared to commonly used power-law sliding laws (such as Weertman type). For the ice shelf removal experiment this leads to 5 and 16 m SLR for the power-law basal sliding and Coulomb friction conditions at the grounding line, respectively. This high-end response is of the same order of magnitude as obtained by DeConto and Pollard (2016) using ice shelf debuttressing caused by hydrofracture and cliff failure.
The atmospheric-oceanic forcing experiments clearly show the dominance of ocean forcing in sea-level response, where significant MISIs occur under relatively mild sub-shelf melt scenarios over centennial timescales (500 years).
Data availability. All datasets used in this paper are publicly available, such as Bedmap2 (Fretwell et al., 2013) and geothermal heat flow data (Purucker, 2013). Results of the RACMO2 model were kindly provided by Melchior Van Wessem. The EISMINT-I benchmark is the first series of ice sheet model intercomparisons aiming at benchmarking large-scale ice sheet models under idealized and controlled conditions (Huybrechts et al., 1996). The first (fixed margin) experiment considers a square grid of 1500 × 1500 km with a flat bed at zero elevation. Grid spacing is taken as = 50 km, leading to 31 × 31 regularly spaced grid points. Starting from zero ice thickness, the model is forced with a constant surface mass balance of 0.3 m a −1 and surface temperature according to T s = 239 K +(8 × 10 −8 )d 3 summit , where d summit is defined as max(|x −x summit |, |y −y summit |), expressed in km. Further boundary conditions for the model are zero ice thickness at the edges of the domain and a constant geothermal heat flux of G = 0.042 W m −2 . The ice temperature is not coupled to the ice flow field and a constant value for the flow parameter of 10 16 Pa −n a −1 is considered.
The f.ETISh model is a 3D Type I model according to the classification scheme in EISMINT-I; i.e. diffusion coefficients for the grounded ice sheet are calculated on a staggered Arakawa B grid. Table A1 lists the comparison with data from other 3D Type I models. Both ice thickness and flux compare very well within error bounds of the sample range (limited to only two to three models in the EISMINT-I benchmark, unfortunately). Also the basal temperature at the divide and along the profile is within the limits given by the EISMINT-I benchmark. Differences can be attributed to the use of the shape functions for the velocity field as well as to the use of a staggered grid for the temperature field, whereby the temperature at the divide and along the profile are interpolated values along the central line.

A2 Moving margin experiment
The moving margin experiment includes ice ablation, hence the presence of an equilibrium line on the ice sheet. This is obtained by defining the climatic conditions bẏ a = min{0.5, h s (R el −d summit )} and T s = 270−0.01h, where d summit is here defined as the radial distance from the centre (in km), and s and R el are 10 −2 m a −1 km −1 and 450 km, respectively (Huybrechts et al., 1996). The steady-state ice sheet according to this experiment does not reach the edge of the domain but is circular in shape. Note that, contrary to the fixed margin experiment, surface temperature is a function of surface elevation and not of the geometrical characteristics of the domain. Surface mass balance, however, remains a function of the distance to the centre of the domain.
Basic characteristics of the experiment are listed in Table A1, and simulated values of ice thickness (h summit ) and basal temperature at the divide (T b summit ), as well as ice flux between divide and margin, are in good agreement with the benchmark. Also the basal temperature profile agrees well Table A1. Comparison of f.ETISh with the EISMINT-I fixed (FM) and moving margin (MM) experiment benchmark based on an ensemble of two to three models (Huybrechts et al., 1996) Figure A1. Homologous basal temperatures along the central line according to the EISMINT-I experiment calculated with f.ETISh (circles) and according to the EISMINT-I benchmark (crosses) for the fixed margin (blue) and moving margin (red) experiment.
with the benchmark and differences can be attributed to the factors listed in Appendix A1.

A3 Transient experiment
Temporal changes in ice thickness/volume and basal temperature are analysed with a forcing experiment, where the surface temperature and mass balance perturbations are defined as follows (Huybrechts et al., 1996): R el = 100 sin 2π t T for moving margin .
The model run starts from the steady-state ice sheet obtained in the previous section and the forcing is applied for a period of 200 ka, with a periodicity of T = 20 and 40 ka, respectively. Results are depicted in Fig. A2 for the fixed margin and in Fig. A3 for the moving margin experiment. Table A2  Table A2. Comparison of f.ETISh with the EISMINT-I fixed (FM) and moving margin (MM) experiment benchmark based on an ensemble of two to three models (Huybrechts et al., 1996)  lists the main characteristics of ice thickness and basal temperature amplitude variations, as well as ice thickness at the divide at the end of the experiment (200 ka). All ice thickness changes (amplitude and phase) as well as the phase in temperature according to the two forcing scenarios are in close agreement with the benchmark. However, amplitude differences for the basal temperatures deviate, but the EISMINT I data sample is rather limited for comparison. The phase of the basal temperature response is in agreement with the benchmark. All other parameters are within the bounds of the benchmark (Table A2).

Appendix B: EISMINT-II benchmark
The EISMINT-II benchmark (Payne et al., 2000) is based on the moving margin experiment of Huybrechts et al. (1996) but includes thermomechanical coupling of the ice flow to the temperature field. Contrary to the EISMINT-I benchmark, inter-model differences are considerably larger, especially with respect to the area of the ice sheet that reaches pressure melting point at the base. The standard experiment consists of a flat bed of the same size as the EISMINT-I benchmark, but with a spatial resolution of 25 km, leading to 61 × 61 grid points. The basic experiment (A in Payne et al., 2000) runs the ice sheet in equilibrium starting from zero ice thickness on the domain and with u b = 0. The climatic conditions are defined as where d summit is defined as in the moving margin experiment as the radial distance from the centre (in km), s and R el are taken as in the moving margin experiment (10 −2 m a −1 km −1 and 450 km, respectively), andȧ max , T min and s T are defined as 0.5 m a −1 , 238.15 K and 1.67 × 10 −2 K km −1 , respectively. Contrary to the moving margin experiment, climatic conditions are independent of ice sheet surface elevation, and hence the mass-balance elevation feedback is excluded.
Six further experiments were carried out: B, C, D, F, G and H (in Payne et al., 2000). They consist in a stepwise change in surface temperature, T min = 243.15 K (B), a stepwise change in surface mass balanceȧ max = 0.25, R el = 425 km (C) and a stepwise shift in equilibrium-line altitude R el = 425 km. Experiments B, C and D start from the steadystate solution of A. Experiment F is similar to A, but starting with a value of T min = 223.15 K (model run starting without ice). Experiment G incorporates basal slip according to a linear sliding law (m = 1 and A b = 10 −3 m a −1 Pa −1 ) with a similar set-up as A. Finally, experiment H is similar to G, but 1872 F. Pattyn: Antarctic sensitivity to sub-shelf melting where sliding is limited to areas that are at pressure melting at the base.
Results for experiments A-H are summarized in Table B1. The majority of parameters are within the bounds of the benchmark, but major differences are related to the basal temperature at the divide. All experiments exhibit a radial pattern in basal temperatures that are at pressure melting point for the outer part of the ice sheet, with a cold spike in the centre of the ice sheet. In all experiments, our temperature spike is slightly less cold than the one given by the benchmark. However, despite this difference, the size of the basal area at pressure melting point is in accord with the benchmark. Again, the main reason for this difference is that temperatures in f.ETISh are calculated on a staggered Arakawa B grid and not exactly at the ice divide. Despite these differences in temperature, ice volume and area coverage are totally in agreement with the benchmark mean.
The emblematic experiments F and H in Payne et al. (2000) displayed an irregular pattern in the basal temperatures of the benchmark for all participating models, leading to cold spikes reaching to the edge of the ice sheet. The pattern was shown to be model dependent and further investigations traced its origin to an interaction between vertical advection (cooling down the base) and strain heating  (Hulton and Mineter, 2000). The pattern was found to be highly dependent on spatial grid resolution due to the lack of membrane stresses in the shallow-ice approximation (Hindmarsh, 2006(Hindmarsh, , 2009). Also f.ETISh produces a similar patterning for this particular experiment despite the approximations in the thermomechanical coupling (using a vertically integrated temperature) and the use of shape functions (Fig. B1).

Appendix C: Modified MISMIP experiments
The capacity of an ice sheet model to cope with the marine boundary, and more specifically migration of the grounding line, is essential in Antarctic ice sheet modelling. Since grounding-line dynamics were elucidated mathematically based on boundary layer theory (Schoof, 2007a(Schoof, , b, 2011, two intercomparison exercises were established. The first one tested grounding-line migration and stability on downwardsloping beds and instability on retrograde slopes for flow-line models (Pattyn et al., 2012), and the second tested the effect of buttressing for two-and three-dimensional ice sheet models . Given that marine ice sheet instability is a crucial feedback process in marine ice sheet behaviour, we performed the flow-line experiments for a planview model set-up. Experiments were carried out for both grounding-line flux conditions SGL and TGL. Ice shelves are included, but without exerting any buttressing strength, i.e. τ xx = τ f . The first experiment is an ice sheet on a seawardsloping bedrock, which in plan view results in a conic bed, defined by (Pattyn et al., 2012) B = 720 − 778.5 750 d summit , where d summit (km) is the radial distance from the centre of the domain. The second experiment consists of an overdeepened section in the bedrock profile, hence the presence of a Table B1. Comparison of f.ETISh with the EISMINT-II experiments (Payne et al., 2000).
The initial ice sheet is obtained for a constant value of the flow parameter A of 10 −16 Pa −n a −1 and a constant surface mass balance ofȧ = 0.3 m a −1 . A grid-size spacing of = 50 km is employed. All other parameters are listed in Tables 1-3. Subsequently, the flow-rate parameter A is altered to a new value to obtain a new steady state, where lower/higher values of A leads to groundingline advance/retreat, respectively. According to theory, a given set of boundary conditions leads to unique steady-state grounding-line positions on a downward-sloping bedrock, while the grounding line never reaches a steady-state position on an upward-sloping bedrock, which is depicted in Fig. C1. For the overdeepened bed, this leads to hysteresis, i.e. multivalued grounding-line positions and ice sheet profiles for the same set of boundary conditions (Figs. C1 and C2). The numerical error was estimated by determining the position of each grounding-line grid cell compared to its radial distance from the centre of the ice sheet (both experiments results in radial ice caps). The mean position of the grounding line and the standard deviation corresponding to each steady state are shown in Fig. C2. Interpolation of the exact position within a grid cell was not considered. All errors are smaller than the nominal grid size of 50 km. The lowest numerical error corresponds to the grounding-line treatment according to the power-law sliding law without the presence of ice shelves (σ ∼ 20 km). Including ice shelves makes the ice sheet more rapidly advance across the unstable section, since ice shelf thickness increases for lower values of A. Associated errors are also larger. Finally, the flux condition for Coulomb friction (Tsai et al., 2015) results in a generally smaller ice sheet, as the ice flux across the grounding line is higher than in the previous case. The ice sheet is also more sensitive to changes in A, i.e. small changes make the grounding-line advance and  Figure C2. Median (upper panel) and standard deviation (centre panel) of steady-state grounding-line positions according to the MISMIP experiments for a circular ice sheet as a function of the flow parameter A (Pa −n a −1 ) for the overdeepened bedrock experiment according to different flux conditions at the grounding line and inclusion or exclusion of ice shelves. Solid lines represent advance and dashed lines represent retreat experiments. The lower panel displays the difference in grounding-line position for the steadystate ice sheets obtained during advance with those obtained during retreat (same parameter values) for the linear-sloping bed (linear) and the overdeepened bed. The latter has two sections with stable steady-state solutions: a small ice sheet (sloping 1) and a large one (sloping 2). retreat more rapidly. Associated errors are smaller for the noshelf experiment but significantly larger for the ice shelf experiment. Given the larger sensitivity, the numerical solution is also less stable compared to the power-law flux condition SGL of Schoof (2007a) and the use of smaller time steps could probably improve the results.
Errors on the advance and retreat grounding-line positions are displayed in the bottom panel of Fig. C2. In all cases, the difference in grounding-line position between advance and retreat is less than 10 km (one-fifth of the spatial resolution of the model). In some cases the error is exactly zero, meaning that the steady-state ice sheets (the one obtained during advance compared to the one obtained after retreat) are exactly the same.