Articles | Volume 12, issue 6
Model experiment description paper
12 Jun 2019
Model experiment description paper |  | 12 Jun 2019

Assessment of sub-shelf melting parameterisations using the ocean–ice-sheet coupled model NEMO(v3.6)–Elmer/Ice(v8.3)

Lionel Favier, Nicolas C. Jourdain, Adrian Jenkins, Nacho Merino, Gaël Durand, Olivier Gagliardini, Fabien Gillet-Chaulet, and Pierre Mathiot

Oceanic melting beneath ice shelves is the main driver of the current mass loss of the Antarctic ice sheet and is mostly parameterised in stand-alone ice-sheet modelling. Parameterisations are crude representations of reality, and their response to ocean warming has not been compared to 3-D ocean–ice-sheet coupled models. Here, we assess various melting parameterisations ranging from simple scalings with far-field thermal driving to emulators of box and plume models, using a new coupling framework combining the ocean model NEMO and the ice-sheet model Elmer/Ice. We define six idealised one-century scenarios for the far-field ocean ranging from cold to warm, and representative of potential futures for typical Antarctic ice shelves. The scenarios are used to constrain an idealised geometry of the Pine Island glacier representative of a relatively small cavity. Melt rates and sea-level contributions obtained with the parameterised stand-alone ice-sheet model are compared to the coupled model results. The plume parameterisations give good results for cold scenarios but fail and underestimate sea level contribution by tens of percent for warm(ing) scenarios, which may be improved by adapting its empirical scaling. The box parameterisation with five boxes compares fairly well to the coupled results for almost all scenarios, but further work is needed to grasp the correct number of boxes. For simple scalings, the comparison to the coupled framework shows that a quadratic as opposed to linear dependency on thermal forcing is required. In addition, the quadratic dependency is improved when melting depends on both local and non-local, i.e. averaged over the ice shelf, thermal forcing. The results of both the box and the two quadratic parameterisations fall within or close to the coupled model uncertainty. All parameterisations overestimate melting for thin ice shelves while underestimating melting in deep water near the grounding line. Further work is therefore needed to assess the validity of these melting parameteriations in more realistic set-ups.

1 Introduction

The majority of grounded ice in Antarctica is drained through its floating extensions advancing in the Southern Ocean. The increase in ice-mass loss since the 1990s has been mostly driven by ice-shelf thinning in the western part of the ice sheet (Paolo et al.2015; Shepherd2018). In the Amundsen and Bellingshausen seas, ice-shelf thinning is due to incursions of Circumpolar Deep Water (CDW) beneath the ice-shelf base all the way to the line boundary between the grounded and floating part of the ice sheet, i.e. the grounding line. These incursions episodically increase the ocean–ice heat flux and drive sub-shelf melting and ice-shelf thinning (Jacobs et al.2011; Dutrieux et al.2014; Jenkins et al.2018 for West Antarctica and Gwyther et al.2018 for East Antarctica). The thinning of floating ice decreases the back force restraining the upstream ice, leading to ice-sheet acceleration (Mouginot et al.2014), ice-surface lowering (Konrad et al.2017), retreating grounding lines (Rignot et al.2014; Konrad et al.2018), and eventually increased sea level rise.

West Antarctic grounding lines often rest on retrograde bed up-sloping towards the ocean (Fretwell et al.2013). This makes the glaciers vulnerable to the marine ice-sheet instability (MISI), which states that an ice sheet starting to retreat over a retrograde bed slope keeps retreating until the slope becomes prograde (Mercer1978; Thomas and Bentley1978; Weertman1974; Schoof2007; Durand et al.2009). Confined ice shelves resist horizontal shearing and potentially stabilise an ice sheet undergoing MISI (Gudmundsson et al.2012; Gudmundsson2013; Haseloff and Sergienko2018). Ice-sheet modelling results suggest that the Pine Island and Thwaites glaciers may have started an unstable retreat (Favier et al.2014; Joughin et al.2014), but the tipping point beyond which MISI occurs is not clearly identified (Pattyn et al.2018).

Ocean warming is currently the main driver of the West Antarctic ice-sheet retreat, and can potentially trigger further MISI (Favier et al.2014; Joughin et al.2014). Using realistic ice-shelf basal melt rates in ice-sheet simulations is therefore crucial. The most comprehensive way to do so consists of using an ocean model that solves the 3-D Navier–Stokes equations in ice-shelf cavities and represents ocean–ice heat exchanges (Losch2008). The existence of strong feedbacks between the cavity geometry, melt rates, and the ocean circulation (De Rydt et al.2014; Donat-Magnin et al.2017) has motivated the development of coupled ocean–ice-sheet models presenting a moving ocean–ice boundary. To date, this kind of coupled model has been used in idealised configurations (e.g. De Rydt and Gudmundsson2016; Asay-Davis et al.2016; Jordan et al.2018; Goldberg et al.2018) or with more realistic configurations representing a single ice shelf (Thoma et al.2015; Seroussi et al.2017). However, the required numerical developments and the relatively high computational cost of the ocean component strongly limit the use of ocean–ice coupled models for long-term simulations of the Antarctic ice sheet.

A much simpler approach to account for oceanic forcing in stand-alone ice-sheet models is to prescribe melting by plugging off-line ocean model outputs (e.g. Seroussi et al.2014). The melt rates cannot evolve with cavity geometry changes. Mengel and Levermann (2014) improved the method by correcting the dependency of the freezing point to a changing ice draft, but it is still unable to account for the dependency on far-field temperature and salinity stratification, and for circulation changes driven by the evolution of the cavity geometry (Donat-Magnin et al.2017). This approach also requires the choice of empirical ad hoc melt rates underneath newly floating ice wherever the grounding line is retreating during the prognostic simulations. To circumvent this issue, Cornford et al. (2015) and Nias et al. (2016) consider the ice-mass flux near and away from the grounding line to build a sound initial melting pattern that depends on the distance to the grounding line and adapts to its further migration. By construction, the melt rates are much larger at the grounding line and decrease exponentially away from it. Spatially and temporally varying melt rates (anomalies) taken from ocean models are added to these initial melt rates to predict future sea level contribution. This latter approach is also empirical and does not account for potential change in oceanic circulation (e.g. due to feedbacks with ice dynamical changes).

The melt rates can also be parameterised using two main approaches, being either an explicit function of depth or a function depending on far-field ocean temperature and salinity. In the first approach (followed by, for example, Favier et al.2014; Joughin et al.2014, with more examples given in Asay-Davis et al.2017), they are computed by a piecewise linear function of depth, and an initial calibration is done to match current observations on average (e.g. using datasets from Rignot et al.2013b; Depoorter et al.2013). The oversimplicity of the depth dependence not only makes the initial pattern very different from the observed pattern, but also leads to a significant overestimation of the grounding-line retreat compared to ocean–ice-sheet coupled models (Seroussi et al.2017; Jordan et al.2018; De Rydt and Gudmundsson2016).

The second approach parameterises the melt rates as a function of ocean temperature and salinity profiles. The simplest parameterisations are mere functions of the difference between the temperature and the melting–freezing point at the ice–ocean boundary, the thermal forcing, using a linear (e.g. Beckmann and Goosse2003; Favier et al.2016) or a quadratic dependency (e.g. DeConto and Pollard2016). More complexity is accounted for in the box model proposed by Reese et al. (2018a) and based on the 1-D ocean-box model from Olbers and Hellmer (2010), and also in the 2-D emulation of a 1-D plume model (Jenkins1991) proposed by Lazeroms et al. (2018).

Assessing these last parameterisations with regard to melt rates computed by a stand-alone ocean model would enable the patterns differences in a static cavity geometry to be investigated. However, the melt-rate pattern also has an effect on the ice-sheet response. The study of Gagliardini et al. (2010) highlights configurations where less melting leads to a grounding line relatively further upstream, or where the same average melting leads to two different ice-sheet responses and grounding-line positions. An ice-sheet model is therefore needed to carry out a meaningful comparison between parameterised and simulated melt rates.

In this paper, we assess several flavours of the aforementioned ocean temperature- and salinity-dependent parameterisations with regard to ocean–ice-sheet coupled simulations. We include the uncertainties arising from the ocean model by considering an ensemble of four ocean–ice coupled configurations. Following an initial calibration that allows further comparisons between parameterised and coupled simulations, we use six one-century far-field ocean temperature and salinity scenarios, which we apply to drive the melting parameterisations in stand-alone ice-sheet simulations and force the members of the ocean ensemble in ocean–ice-sheet coupled simulations. Overall, the MISOMIP (Asay-Davis et al.2016) framework is used to perform 138 one-century simulations (19 sub-shelf melt parameterisations + 4 coupled members × 6 scenarios).

The paper is organised as follows. The second section describes the models: the ice-sheet model Elmer/Ice, the ocean model NEMO and the framework for coupling those two models. The section also describes the sub-shelf melt-rate parameterisations and the members of the ocean–ice ensemble. The third section describes the experiments, including the reference set-up of the ocean–ice-sheet system, the initial calibration of the parameterised and coupled simulations, and the set of far-field ocean temperature and salinity scenarios. Then in the fourth section, we detail the results with regard to sea-level contribution and sub-shelf melting evolution, and in the fifth section, we discuss the use of sub-shelf melt parameterisations in stand-alone ice-sheet modelling at a regional or a global scale.

2 Models

2.1 The ice-sheet model, Elmer/Ice

We perform the ice-sheet simulations with the finite-element ice-sheet model Elmer/Ice (Gagliardini et al.2013). The ice rheology is non-linear and controlled by Glen's flow law (Appendix A), enabling the deviatoric stress tensor to be linked with the strain rate tensor from which ice velocities are retrieved. The version of the ice-sheet model used solves the SSA* solution, a variant of the L1L2 solution of Schoof and Hindmarsh (2010), solving the shallow shelf approximation of the Stokes equations and accounting for vertical shearing in the effective strain rate. The SSA* approximation was recently implemented in Elmer/Ice following the work of Cornford et al. (2015).

To calculate the basal friction, the grounding line position is calculated from hydrostatic equilibrium and can thus be located anywhere within an element. We use a sub-element parameterisation to affect basal friction to the part of the element that is grounded by increasing its number of integration points (equivalent to the SEP3 method in Seroussi et al.2014). The basal friction is computed by a Schoof-like friction law based on the theoretical work of Schoof (2005) applied to a linear ice rheology, and which was extended to a non-linear rheology by Gagliardini et al. (2007). The Schoof friction law (Appendix A) depends on the effective pressure, the difference between the ice overburden pressure and the basal water pressure, here approximated by the ocean pressure. This friction law therefore exhibits two asymptotic behaviours, behaving as a non-linear power law away from the grounding line and as a Coulomb friction law near the grounding line, and thus ensuring a smooth transition of stress state near and at the grounding line. The Schoof friction law was recently compared to various other types of friction laws commonly used in ice-sheet modelling, for an idealised framework (Brondex et al.2017) and a real drainage basin (Brondex et al.2018).

Melting is applied to floating nodes but not to grounded nodes, meaning that the first floating element (partially or not) may be affected by melting. The mesh grid is unstructured and made of triangles, the size of which is about 500 m in the vicinity of the grounding line and up to 4 km away. The Elmer/Ice configuration is identical for parameterised and coupled simulations.

2.2 Ocean melting from a 3-D ocean–ice-sheet coupled model

The melt rates beneath the ice shelf are either parameterised or computed through the coupling of NEMO and Elmer/Ice. Here we describe the ocean model and the ocean–ice-sheet coupling framework.

2.2.1 The ocean model, NEMO

We make use of the 3-D primitive-equation ocean model NEMO-3.6 (Nucleus for European Modelling of the Ocean; Madec and NEMO-team2016). NEMO solves the prognostic equations for the ocean temperature, salinity, and velocities and includes ice-shelf cavities (Mathiot et al.2017). The sub-shelf melting is parameterised through the so-called “three equations” representing (1) the heat balance at the ice–ocean interface accounting for phase change, turbulent exchange in water, and diffusion in the ice; (2) the salt balance accounting for freezing, melting, and turbulent exchange; and (3) the pressure and salinity dependence of the potential temperature at which seawater freezes (Hellmer and Olbers1989; Holland and Jenkins1999; Losch2008; Jenkins et al.2010). In this parameterisation, we assume a constant top-boundary-layer (TBL) thickness along the ice-shelf draft (Mathiot et al.2017), and we use a velocity-dependent formulation in which the heat exchange velocity is defined as follows:

(1) γ T = Γ T C d ( u TBL 2 + u tide 2 ) ,

where uTBL the TBL-averaged velocity resolved by NEMO, ΓT is the non-dimensional heat exchange coefficient, Cd the non-dimensional drag coefficient and utide is a uniform background velocity representing the main effect of tides on ice-shelf melting (Jourdain et al.2018). The values of ΓT, Cd, and utide are given in Table 1.

The ocean configuration used in this study is very similar to the ISOMIP+ configuration described by Asay-Davis et al. (2016): we use a linearised equation of state and the only lateral boundary condition is a temperature and salinity restoration along the vertical boundary representing offshore conditions; neither sea ice nor atmospheric forcing nor tides are represented. The only differences with the general MISOMIP protocol is that we use different temperature and salinity restoration and initial conditions (Sect. 3.3). We use a variety of resolutions and parameters for NEMO to build an ensemble of NEMO-Elmer/Ice coupled simulations as described in Sect. 2.2.3.

2.2.2 The ocean–ice-sheet coupled model framework

We couple NEMO and Elmer/Ice, meaning that Elmer/Ice sees sub-shelf melt rates calculated by NEMO, while NEMO sees the ice-shelf geometry resulting from the ice dynamics resolved by Elmer/Ice. A given coupling period (typically of few months) is first covered by the ocean model with the cavity geometry from the end of the previous coupling period; then, the period is covered by the ice-sheet model forced by the oceanic melt rates averaged over this coupling period in order to conserve mass as much as possible (Fig. 1).

Figure 1NEMO-Elmer/Ice coupling framework. T and S stand for temperature and salinity.


As the respective grids of the two models differ, some interpolation is required for each exchange. Following each NEMO run, Elmer/Ice restarts from its previous time step (ice geometry and velocities). The melt rates provided by NEMO are bi-linearly interpolated onto Elmer/Ice's unstructured grid. A multiplicative correction factor computed over the entire ice shelf ensures that the same mass flux is seen by the two models (this factor is very close to one in our case). In case Elmer/Ice has a floating element but the water column is too thin to be captured by NEMO (a minimum thickness of 20 m allows NEMO to have a minimum of two vertical cells under the partial cell conditions, Mathiot et al.2017), the melt rate seen by Elmer/Ice is set to zero.

Every coupling period, NEMO restarts with temperature, salinity, and velocities from its previous time step using the updated geometry from Elmer/Ice. If new ocean cells appear (previously masked ice cells), temperature and salinity are an average of the four closest wet cells (horizontally if possible, vertically extrapolated otherwise), and ocean velocities are set to zero. To avoid the generation of spurious barotropic waves as a result of sudden changes in water column thickness, we impose a conservation of barotropic velocities across the step change in the ice-shelf geometry. We also conserve the sea surface height (SSH) value for all the water columns, and if a new water column is created, SSH is an average of the four closest wet cells.

We use the same initial state for Elmer/Ice as in MISOMIP (Asay-Davis et al.2016), i.e. a steady state obtained with zero melt, and NEMO is spun up for 5 years with this initial ice-shelf geometry before being coupled to Elmer/Ice. The respective time steps of Elmer/Ice and NEMO are 1 month and 200 s, and the coupling period ranges between 2 and 6 months, depending on the configuration. We performed a sensitivity study following the MISOMIP protocol (Asay-Davis et al.2016), which indicates very little sensitivity to coupling periods between 1 month and 1 year, with less than 3 % difference in sea-level contribution after 100 years (Appendix B).

2.2.3 The ensemble of ocean configurations within the coupled framework

While the NEMO ocean model is much more representative of the ocean physics than any sub-shelf melting parameterisation, there are still processes like turbulence and convection that need to be parameterised. The model is also sensitive to both the horizontal and vertical resolutions. To account for the consequent ocean model uncertainty, we consider four NEMO configurations with the varying parameters listed in Table 1. For each coupled configuration, the ΓT parameter is adjusted following the exact ISOMIP+ calibration protocol after 4 years of ocean spin-up with a steady ice-shelf draft (more details of the protocol relevant to our study are given in Sect. 3.2, and the protocol is fully described in Asay-Davis et al.2016; Sect. 3.2.1).

Table 1Ocean parameters used for the four NEMO-Elmer/Ice coupled simulations. Δx is the horizontal resolution, TCPL is the ocean–ice-sheet coupling period and Δz is the nominal vertical resolution. The actual resolution near the sea floor or ice shelf draft can be smaller due to the use of partial steps, but the TBL thickness is always equal to Δz (i.e. TBL quantities are averaged over several levels in the case of partial steps). ΓT and utide are defined in Eq. (1), and the salt exchange coefficient ΓS is taken as ΓT∕35. Also defined in Eq. (1) is the drag coefficient Cd=2.5×10-3 . The stable vertical diffusivity and viscosity coefficients (Kstab and νstab respectively) are either constant, at the same values as in Asay-Davis et al. (2016), or calculated through the TKE scheme with the same parameter values as in Treguier et al. (2014). Convection is parameterised through enhanced diffusivity and viscosity (Kunstab and νunstab respectively) in case of static instability (0.1 m2 s−1 as Asay-Davis et al.2016 and 10 m2 s−1 as Treguier et al.2014). The remaining parameters are exactly the same as in the common ISOMIP+ configuration described in Asay-Davis et al. (2016).

Download Print Version | Download XLSX

2.3 Ocean melting from ocean-dependent sub-shelf parameterisations

All the parameterisations are linked to ambient temperature and salinity vertical profiles in the far-field ocean. The stand-alone ice-sheet simulations start from the same initial state as for the ocean–ice-sheet coupled simulations. The parameterisations respond instantaneously to changes in ambient temperatures and salinities; i.e. they do not account for ocean circulation timescales (e.g. water residence time in ice-shelf cavities, Holland2017). None of the parameterisations account for the Coriolis effect or for bathymetric features (e.g. sills, channels). To avoid areas of very thin ice that would affect the stability of the ice-sheet model, melting is not permitted wherever the ice base is shallower than 10 m depth.

2.3.1 Simple functions of thermal forcing

The following three parameterisations are based on an expression for the ice–ocean heat transfer that is analogous to the one used in more complex ocean circulation models (Grosfeld et al.1997). However, they make the simplifying assumption that the thermal forcing across the ice–ocean boundary layer can be determined directly from far-field ocean conditions. Thus, cooling of the water as it is advected from the far field into the cavity and then mixed into the ice–ocean boundary layer is accounted for simply through the choice of an effective heat transfer coefficient.

The linear, local dependency on thermal forcing assumes a balance between vertical diffusive heat flux across the ocean cavity top boundary layer and latent heat due to melting and freezing. Its formulation is based on Beckmann and Goosse (2003) and written as follows:

(2) M lin = γ T ρ sw c po ρ i L i ( T o - T f ) ,

with γT the heat exchange velocity (aimed at being calibrated; see Sect. 3.2), ρsw and ρi the respective densities of ocean water and ice, cpo the specific heat capacity of the ocean mixed layer, and Li the latent heat of fusion of ice (Table 2). The melting–freezing point Tf at the interface between the ocean and the ice-shelf basal surface is defined as follows:

(3) T f = λ 1 S o + λ 2 + λ 3 z b .

The practical salinity So and the potential temperature To are taken from the far-field ocean as detailed below in this section; zb is the ice base elevation, which is negative below sea level; and the coefficients λ1, λ2, and λ3 are respectively the liquidus slope, intercept, and pressure coefficient.

The linear formulation with a constant exchange velocity assumes a circulation in the ice-shelf cavity that is independent from the ocean temperature. This assumption is neither supported by modelling (Holland et al.2008; Donat-Magnin et al.2017) nor by observational studies (Jenkins et al.2018) that suggest a more vigorous circulation in response to a warmer ocean, subsequently increasing melt rates.

The quadratic, local dependency on thermal forcing accounts for this positive feedback between the sub-shelf melting and the circulation in the cavity (Holland et al.2008), using a heat exchange velocity linearly depending on local thermal forcing. The formulation is written as follows:

(4) M quad = γ T ρ sw c po ρ i L i 2 ( T o - T f ) 2 .

These last two parameterisations were used in numerous studies (e.g. review in Asay-Davis et al.2017). As the ocean properties used to calculate melting for every draft point are taken at the very same point, they are tagged as local.

The quadratic, local and non-local dependency on thermal forcing is a new parameterisation assuming that the local circulation (at a draft point) is not only affected by local thermal forcing, but also by its average over the ice basal surface, which is written as follows:

(5) M + = γ T ρ sw c po ρ i L i 2 ( T o - T f ) T o - T f .

This formulation is inspired by Jourdain et al. (2017), who showed an overturning circulation proportional to total melt rates. It is equivalent to assuming that melting is first generated by local thermal forcing, and that this first-guess melting generates a circulation at the scale of the ice-shelf cavity that feeds back on melt rates. In other words, this formulation reflects the three equations with a uniform exchange velocity that is proportional to the cavity-average thermal forcing.

In Eqs. (2), (4), and (5), the values of To and So are either depth-dependent or taken from a constant depth in the far field (Sect. 3.3 details the different far-field ocean temperature and salinity vertical profiles). The former situation (for which To=To(z) and So=So(z)) assumes a horizontal circulation between the far-field ocean and the ice draft that would transport constant ocean properties. This can be viewed as an asymptotic case where the circulation in the cavity is driven by tides rather than melt-induced buoyancy forces, which is equivalent to the aforementioned three equations with a constant and uniform velocity along the ice base. Alternatively, in the latter situation, To and So are taken at either 500 m or 700 m depths, i.e. near the sea floor. This assumes that ocean water is advected into the cavity along the sea floor up to the grounding line, then upward along the ice base with constant ocean temperature and salinity.

The value of Tf is therefore calculated with either So(z) in the first option, or So(500) or So(700) in the second option (in a consistent way with To), but with the local ice base depth. For each far-field ocean temperature and salinity profile, we thus run three Elmer/Ice simulations for each simple function of the thermal forcing.

Table 2Physical parameters, model grid resolutions, and coupling period.

Download Print Version | Download XLSX

(Reese et al.2018a)(Lazeroms et al.2018)

Table 3Parameterisations used to compute melting in stand-alone ice-sheet simulations. The last column lists the calibrated γT obtained from the WARM profile, except for the plume parameterisation where a multiplicative coefficient α is used instead.

Download Print Version | Download XLSX

2.3.2 More complex functions of thermal forcing

The following two parameterisations attempt to improve on the above by including a representation of some of the processes that determine the temperature within the ice–ocean boundary layer. Cooling of the water as it is advected into the cavity is still neglected, so that the waters incorporated into the boundary layer have far-field properties. However, cooling of the boundary layer by melting at depth, the rise of the waters along the ice shelf base, and the change in the freezing point with depth are all considered with different levels of detail. Critically, including such processes enables these parameterisations to simulate regions of basal freezing, something that the simple functions of far-field temperature cannot reproduce.

The box parameterisation was developed by Reese et al. (2018a) based on the analytical steady-state solution of the box model of Olbers and Hellmer (2010). The latter, initially developed for a 2-D cavity, represents the buoyancy-driven advection of ambient ocean water into the ice-shelf cavity at depth up to the grounding line, then upward along the ice draft in consecutive boxes. The melt rates are given by the following:

(6) BM = γ T ρ sw c po ρ i L i ( T k - T f , k ) ,

where the k subscript indicates properties evaluated in each box. Those properties account for the transformation of ocean temperature and salinity in consecutive boxes through heat and salt turbulent exchange across the ocean boundary layer underneath ice shelves. Hence, the box model is entirely driven by ocean temperature and salinity near the sea floor. Unlike plume models, the box model does not entrain deep water all along the upward transport, it advects deep water from the open ocean to the grounding zone then transports it upward. Therefore, this parameterisation produces maximum melt rates near the grounding line.

A key assumption is that the overturning circulation (i.e. volume transport through the boxes) is taken proportionally to the density difference between the ambient ocean (open ocean seaward of the ice shelf) and the deepest box including an ocean–ice interface. Similarly to the simple parameterisations, the box model assumes constant heat and salt exchange velocities.

In their implementation, Reese et al. (2018a) calibrated both the heat exchange and overturning coefficients to obtain realistic melt rates for both the Pine Island and Ronne-Filchner ice shelves. Here, we keep the overturning coefficient used by Reese et al. (2018a), and we calibrate the effective heat exchange velocity in the same way as the other parameterisations (Sect. 3.2).

In our implementation of the box model, the calving front position that is used to build the boxes' positions is considered to be at either x=640 km or defined by the 10 m depth contour, the limit below which no melting is permitted for the ice-sheet model. In the Reese et al. (2018a), the dependence of sub-shelf melting on the local pressure due to the vertical ice column induces a lack of energy conservation. We thus decided not to implement this dependence, resulting in a uniform melting within each box.

For each temperature and salinity scenario, we run six Elmer/Ice simulations using the box parameterisation, with either 2, 5, or 10 boxes, and with ocean temperature and salinity taken at constant depths of either 500 or 700 m.

The plume parameterisation developed by Lazeroms et al. (2018) emulates the 2-D behaviour of the 1-D plume model proposed by Jenkins (1991). This model describes the evolution of a buoyant plume originating from the grounding line with zero thickness and velocity, and temperature and salinity taken from the ambient ocean. Away from the grounding line, the thickness, velocity, temperature, and salinity of the plume evolve through advection, turbulent exchange across the ocean boundary layer underneath the ice shelf, and entrainment of deep water. Among the melt formulations presented in this paper, the plume parameterisation is the only one to include velocity-dependent heat and salt exchange velocity. No background or tidal velocity is prescribed, so turbulent exchanges and melt rates are zero right at the grounding line.

The plume model can be scaled with external parameters and applied to 1-D ice drafts of any slope, ambient temperature, and salinity (Jenkins2014). The melt rates are given by the following:

(7) PME = α M o g ( θ ) ( T o - T f , gl ) 2 M ^ ( X ^ ) ,

where PME means plume model emulator; Mo is an overall scaling parameter; g(θ) is a function of the ice-shelf basal slope θ, but also of physical constants (heat exchange coefficient, drag coefficient, and entrainment); the f, gl subscript indicates the freezing temperature at the depth of the grounding line; and the final term gives the scaled melt rate, M^, as a universal function of scaled distance, X^, that was derived from empirical fitting of results generated by the full plume model on idealised geometries (Jenkins2014). α is a multiplicative coefficient that will be used for calibrating purposes in our study (see further details in this section). The far-field temperature used here is taken at the depth of the grounding line, as in the box model, and enters the parameterisation explicitly because the subsequent evolution of the ice–ocean boundary layer temperature through entrainment of the far-field ocean, melting, and freezing is captured through the slope-dependent scaling and the universal function. The non-linear dependence on temperature arises because the melt rates depend on the product of plume temperature and plume speed. The latter is a function of the plume buoyancy, which is itself linearly dependent on plume temperature. The physical basis for the scaling is discussed further in Appendix C, but we note here that when the ice-shelf basal slope and far-field conditions are non-uniform, there is no longer a unique choice for those variables in the parameterisation, and choices other than the ones used in this study are equally valid.

Another major issue with the plume parameterisation is the transition from a 1-D to a 2-D ice draft. It is indeed difficult to identify the pathway from a given location of the ice draft to the grounding line point where the plume has emerged, which is enhanced by the fact that several plumes may end up at a given location. To define effective pathways, we apply 4 empirical methods that are all based on different calculations of effective values for the grounding line depth and the basal slope. The first method was originally published in Lazeroms et al. (2018) and applied to a structured grid. The second method was proposed in the corresponding discussion paper but finally discarded to simplify the publication. The last two methods propose simpler ways to calculate the effective grounding-line depth and basal slope. All the methods and their adaptation to unstructured grids are described in Appendix D.

The plume parameterisation from Lazeroms et al. (2018) includes a heat exchange coefficient that is a function of the plume velocity along the ice-shelf base, which is similar to the ocean model but not to the other parameterisations. The complexity of this parameterisation motivated us to calibrate it by adding a multiplicative coefficient α (Table 3) to the melt expression (Eq. 7) rather than calibrating physical parameters.

For each temperature and salinity scenario, we run four Elmer/Ice simulations with the plume parameterisation, using the four aforementioned methods to calculate the effective plume pathway (Appendix D) and with ambient temperature and salinity taken at the effective grounding line depth (as defined in  Lazeroms et al.2018).

3 Experiments

3.1 Initial geometry and set-up

We simulate the evolution of an ideal ice-sheet inspired by the Pine Island glacier in West Antarctica. The domain is the same as the MISOMIP domain for the coupled simulations and as the MISMIP+ domain for the stand-alone ice-sheet simulations (Asay-Davis et al.2016). The ice sheet is marine based and its grounding line rests on a retrograde bed sloping upward towards the ocean. The entire domain, including the ice sheet and the ocean, is 800 km long and 80 km wide (Fig. 2). The ice-sheet calving front is located at x=640 km, while the remaining domain, up to x=800 km, and also the cavity beneath the ice shelf are filled with ocean water. The ice sheet is in equilibrium state with an accumulation rate of 0.3 m a−1 and no sub-shelf melting, as required by MISMIP+, using the ice-sheet configuration detailed in Sect. 2.1. The initial grounding line central position is x=450 km.

Figure 2Initial ice-sheet in equilibrium calculated by Elmer/Ice with an accumulation rate of 0.3 m a−1 and no sub-shelf melting as required by the MISMIP+ protocol (Asay-Davis et al.2016). (a) Side-view geometry in the central flow line, also indicating the position of the ocean restoration used by the ocean model, and the velocity magnitude along the central flow line shown in panel (b). (b) Velocity magnitude seen from above. The black solid line indicates the grounding line. (c) Cross section of the ice sheet at x=480 km.


3.2 Initial state and calibration

The initial calibration purpose is to assess whether the parameterisations represent the response of melt rates to changing ocean temperature and salinity. We thus make sure that all the parameterised and coupled configurations produce the same melting average for the WARM profile of MISOMIP (Fig. 3Asay-Davis et al.2016).

This average is obtained through a spin-up of the ocean model applied to the initial ice-shelf draft (Fig. 2 and Sect. 3.1) and performed before further coupled simulation (Fig. 1). We follow the ISOMIP+ protocol (Asay-Davis et al.2016) to achieve the required sub-shelf melt rate average of 30±2 m a−1 below 300 m depth after 4 years of ocean spin-up. The value of ΓT, which is not known with accuracy and is usually calibrated in ocean models (Asay-Davis et al.2016; Jourdain et al.2017), is therefore adjusted to achieve these melt rates (Table 1). The remaining steps of our calibration, described here below, differ from the ISOMIP+ protocol and are specific to our study. We compute the melting average of all four configurations over the ice draft (excluding parts shallower than 10 m for which no melt is applied), which gives mt8.5± m a−1. These four spin-ups will thus be used as initial states for subsequent coupled simulations.

Then, mt is used as a target for stand-alone ice-sheet simulations forced by the WARM profile from the ISOMIP+/MISOMIP protocol. For the parameterisations in which γT is constant (Eq. 1), we achieve the target by adjusting γT. For the plume parameterisation, which accounts for a top boundary layer velocity, we adjust the value of the multiplicative coefficient α (calibrated values shown in Table 3) to achieve the same target (see Sect. 2.3).

The reason why we did not calibrate the parameterisations to reproduce the average melt rates below 300 m as done in MISOMIP is because all of them produce substantial melt rates underneath the shallowest parts of the ice shelf, as opposed to the ocean models. To emphasise this point, we also performed the simulations with the calibration done as in MISOMIP below 300 m depth, the results of which are given in Appendix F.

The WARM profile was put forward in MISOMIP because it enables a short spin-up of the ocean model, which is useful for calibration purposes as here. After this calibration phase, we keep the calibration reported in Table 2 for all the one-century scenarios described in Sect. 3.3.

3.3 The set of ocean temperature and salinity scenarios

We consider the following six scenarios over a century (Fig. 3), the first two being kept constant, and the other four linearly evolving in time:

  • Warm0 resembles the present-day typical Amundsen Sea conditions (Dutrieux et al.2014). There is no temporal change of temperature and salinity profiles.

  • Warm1 starts from the Warm0 profile and then the temperature uniformly increases by 1 C per century. The salinity profile is constant in time.

  • Warm2 is similar to Warm1 but the warming rate increases with depth, from zero in the surface layer to 1 C per century below the deep thermocline. The salinity profile is constant in time.

  • Warm3 starts from the Warm0 profile and undergoes a 200 m uplift of both the thermocline and the halocline.

  • Cold0 resembles a cold cavity such as beneath the Ronne-Filchner ice shelves. There is no temporal change of temperature and salinity profiles.

  • Cold1 starts from the Cold0 profile and then warms to reach a warm cavity state within a century. The salinity is also increased.

Figure 3Far-field ocean temperature (a) and salinity (b) profiles scenarios in front of the cavity. The WARM profile is used for calibrating the initial state of parameterised and coupled simulations. The Warm0 and Cold0 scenarios are constant in time, while the others evolve linearly in time following the arrows. The Warm1, Warm2, and Warm3 scenarios start with the Warm0 profile and end up after a century in their respective profiles, while the Cold1 scenario starts from the Cold0 profile. In (b), profiles Warm0, Warm1, and Warm2 are equal. Thermal forcing is calculated from the far-field temperature and salinity and applied to the ice-shelf draft, (c) assuming horizontal circulation between the far-field ocean and the cavity or assuming that the circulation is driven by oceanic properties at (d) 500 m and (e) 700 m depths. Profiles from the (c) panel are superimposed to panels (d) and (e) as a watermark for comparison purposes. The Warm1 and Warm2 profiles are equal in panel (e).


These profiles are slightly more realistic than in MISOMIP. They all include a thermocline, because its importance in ice-shelf melting has been pointed out by previous studies (e.g. De Rydt et al.2014). The Warm0 profile corresponds to a linear representation of the average hydrographic profiles measured in front of Pine Island glacier (Dutrieux et al.2014). By contrast, the Cold0 profile represents typical cold-cavity conditions in which deep ocean convection associated with sea ice formation prevents the stratification (e.g. for the Ronne-Filchner and Ross ice shelves). The Warm1 scenario leads to 1 C warming at all depths after 100 years, which corresponds to the upper 80th to 90th percentile of ocean warming projected in the Amundsen Sea by 33 CMIP5 models (Appendix E). The Warm2 scenario is more conceptual and assumes that the sea ice cover will persist over 100 years, i.e. that the ocean surface remains close to the freezing point while the subsurface gets warmer. The Warm3 scenario is inspired by the study of Spence et al. (2014) suggesting that poleward shifting winds over the 21st century will uplift the coastal thermocline due to decreased Ekman downwelling. Last, the Cold1 scenario is an idealised representation of the ocean tipping point described by Hellmer et al. (2012, 2017), in which the Ronne-Filchner cavities switch from a cold to a warm state.

The salinity profile is unchanged throughout Warm0, Warm1, and Warm2 and is sufficiently stratified to keep a stable density profile. In the Warm3 scenario, the halocline is lifted together with the thermocline to mimic an Ekman-driven uplift of the pycnocline, and in Cold1, the stratification in salinity is increased linearly in time to keep a stable stratification when the cavity switches from cold to warm states. Note that none of the temperature profiles account for a salinity compensation (as opposed to the MISOMIP protocol), so the density profile is different in each scenario.

Figure 3c–e show the thermal forcings applied to stand-alone ice-sheet simulations for the different hypotheses for temperature and salinity inputs (Sect. 2.3), while Table 3 summarises the ensemble of sub-shelf melting parameterisations.

4 Results

4.1 Melting patterns resulting from the initial calibration

The calibrated parameters are given in Table 3 and the melting patterns are shown in Fig. 4 (not all the patterns are shown). The patterns obtained from the coupled and parameterised simulations are quite different, even though all of them result in similar cavity melt rates. The coupled simulations give the most melting below approximately 300 m depth and almost no melting near the ocean surface, which also highlights why the calibration was performed below 300 m depth in ISOMIP+ (Asay-Davis et al.2016). The parameterised simulations give significant melt rates at all depths.

Figure 4Diagnostic sub-shelf melt rates obtained through the calibration process by forcing the coupled and the parameterised models with the WARM profile from Asay-Davis et al. (2016). All the ocean members are represented (last column) but not all the parameterisations (first three columns). The average melting for every parameterisation equals 8.5 m a−1, while being in the range 8.5± 1 m a−1 for the ocean members. In the PME1 panel the 200, 300, and 400 m draft contours are shown. The grounded ice is coloured in grey.


Near the grounding line, melt rates higher than 50 m a−1 are predicted by all coupled simulations, while this value is only and hardly reached by the Mquad parameterisation and never reached in the other cases. Away from the grounding line, where the ice shelf is also thinner, melt rates are close to zero for the coupled simulations while they mostly remain above 10 m a−1 when parameterised. Such differences in melt rate patterns are expected to induce diverging responses from the ice sheet (Gagliardini et al.2010; Reese et al.2018b).

While the patterns in the coupled simulations are quite similar to each other, the parameterised patterns differ to various extents. The parameterisations that have a simple dependence on thermal forcing (i.e. Mlin, Mquad, and M+) compute the highest melt rates at depth, which also falls close to the grounding line in the central flow line. They also result in a rather uniform pattern when the basal surface is closer to the sea surface, which occurs away from the grounding line in the central flow line but also close to the grounding line on the sides of the ice shelf, where two bits (or horns) of grounded ice penetrate seaward. The range of melt rates is wider for the Mquad parameterisation, thinner ice being less melted and thicker ice being more melted, compared to Mlin and M+. The Mlin and M+ patterns are similar by construction because the melting average is driven by the (ToTf) term, which appears only once in the two respective formulations. However, the respective calibrations are different (Table 3) because of the term ToTf appearing in M+ only, and the sensitivity to ocean warming will therefore be different.

The implementations of the 2-D plume emulator produce quite different patterns between PME1, PME2, and PME3 on the one hand and PME4 on the other hand, mostly because the latter is highly asymmetric. In the first three implementations, the different approaches adopted to calculate the effective depth and angle (Lazeroms et al.2018) all result in very similar patterns. They all induce zero to small melt rates near the central grounding line because the valid directions are associated with low basal slopes. However, along the sides of the main trunk, on the inner side of the horns, the melt rates get higher at the grounding line because the plumes mostly emerge from the central, much deeper part of the grounding line, and not from the sides where the basal surface is higher than the draft point (PME1 in Fig. 4). Farther away, PME1 and PME2 produce a slight decrease in melting near the calving front, which reflects the empirical scaling with the distance to the grounding line made in Lazeroms et al. (2018) and may not be adapted to our relatively small ice shelf. In the PME3 parameterisation, the plume arises only from the deepest grounding line, whatever the position in the ice draft. On the external sides of the domain, it induces strong melting compared to PME1 and PME2 for which the plumes can also come from less deep parts of the cavity and mitigate the melt rates.

Similarly to the Mlin, Mquad, and M+ parameterisations, the box parameterisation produces its highest melt rates near the grounding line. Away from the grounding line, the melt rates get lower to end up with the lowest values close to the calving front. The larger the number of boxes, the larger the melt rates near the grounding line, and the smaller the melt rates near the calving front.

4.2 Ice-mass loss and sub-shelf melt rates

The initial ice sheet is built within the framework of MISMIP+ (Asay-Davis et al.2016), requiring no sub-shelf melting, and is thus in equilibrium under such conditions. The simulations thus all start with an initial dynamical adjustment of the ice-sheet geometry to new ocean conditions (Fig. 4), which generates a melting pulse despite the 5 years of ocean spin-up. The adjustment is larger for relatively warmer scenarios (Figs. 5, 6). The pulse is therefore much lower and hardly visible for the Cold1 scenario and only shows for the coupled simulations and not for the parameterised simulation for the Cold0 scenario. For the Warmi scenarios, the peak of the pulse yields similar melting of up to 130 Gt a−1 for parameterised and coupled simulations. However, it lasts longer for the former, about 20 a, than for the latter, about 5 a. The pulse in coupled simulations quickly adds a lot of fresh water in the cavity, which further decreases melting. Such feedback is either not or poorly accounted for in the parameterisations, thus increasing the duration of the pulse compared to the coupled simulations.

Figure 5Total melt rates for simple parameterisations (Sect. 2.3.1). The coupled simulations are shown in solid light grey. The coloured lines correspond to parameterised simulations. The black solid lines correspond to a 50 % underestimation or overestimation compared to the average of coupled runs members.


Figure 6Similar to Fig. 5 but for more complex parameterisations (Sect. 2.3.2).


In the Cold0 scenario, almost all parameterised and coupled simulations produce constant melt rates. Only the Mlin parameterisations produce very high melt rates at the start and decrease monotonically afterwards. In the other constant scenario, which is Warm0, the pulse is followed by a decrease in melting, which becomes constant after tens of years for most parameterisations, as opposed to the coupled simulations where the melt rates slightly increase up to the end. In the other scenarios, which are all warming in some way, the pulse is always followed by a melting minimum, after which almost all the parameterised melt rates slightly increase up to the end (there are few exceptions where they are more constant, e.g. Mlin_700 forced by Warm3). Finally, the Warmi scenarios end up with between 40 and 175 Gt a−1 of melting and the Cold1 scenario with between 50 and 100 Gt a−1 of melting. This means that the ice-sheet is contributing 4 to 12 mm to the sea level equivalent mass for the Warmi scenarios, 2 to 4 mm for the Cold1 scenario and 0.5 to 3 mm for the Cold0 scenario (Figs. 7, 8).

Figure 7Sea level contribution (SLC) for simple parameterisations (Sect. 2.3.1). The coupled simulations are shown in solid light grey and their envelope in grey shading. The coloured lines correspond to parameterised simulations. The black solid lines correspond to a 50 % underestimation or overestimation compared to the average of coupled run members.


Figure 8Similar to Fig. 7 but for more complex parameterisations (Sect. 2.3.2).


For the Warmi scenarios, the parameterisations in general tend to overestimate the melting close to the sea surface and underestimate it at depth. This results in initially melting a large part of thinner ice, which makes overall melting higher compared to coupled simulations. Along with the disappearance of thinner ice, the overall melting becomes progressively lower than for coupled simulations. In the end, this results in lower sea level contribution (SLC) from the parameterised simulations, apart from few exceptions. In the Coldi scenarios, melting is never high enough to completely remove thin ice and the SLC from parameterised simulations is more in agreement with the coupled simulation on average.

The uncertainties linked to the ocean model are emphasised by the spread of SLC calculated from the coupled model. The spread is about ± 10 % around the average for all the scenarios except for the Cold0 and Warm2 scenarios, where it is about ± 20 %, respectively. A larger spread of about ± 30 % for the Warmi scenarios, and about ± 50 % and ± 100 % for the Cold1 and Cold0 scenarios, respectively, is obtained from the parameterisations, which reflects the wide variety of approaches and indicates that it makes sense to inter-compare parameterisations with respect to the coupled model.

Whatever the type of hypothesis for the depth at which the far-field ocean temperature and salinity profiles are taken (Sect. 2.3), the Mlin parameterisations tend to largely overestimate the melt rates for the Coldi scenarios and underestimate them for the Warmi scenarios, leading to respective overestimation and underestimation of SLC. This reflects a poor representation of melting by these parameterisations when the change in ocean forcing is too large.

The Mquad parameterisations give melting in fair agreement with coupled results for the Coldi scenarios. For the Warmi scenarios, the tendency is a slight underestimation of SLC using the Mquad and Mquad_700 parameterisations, and a larger underestimation using Mquad_500. Compared to the Mlin parameterisations, it behaves much better and for a larger range of scenarios. All the Mquad parameterisations behave quite well when confronted with a rise in the thermocline (Warm3 scenario), apart from Mquad_700, which slightly underestimates SLC.

The M+ parameterisation results are almost as close to the coupled simulations as the Mquad parameterisations for the Coldi scenarios, and closest for the Warmi scenarios. Regarding all the scenarios, this makes this parameterisation the best among simple parameterisations. When the far-field ocean temperature and salinity profiles are taken at depth, the results are comparable to the Mquad_500 and Mquad_700 parameterisations, thus slightly underestimating SLC.

Forcing a parameterisation by the far-field depth-dependent or the constant depth ocean properties changes the thermal forcing at the ice–ocean interface (Fig. 3) but also the initial calibration (Table 3). Considering a constant depth for instance, the deeper the considered depth, the larger the thermal forcing, but also the lower the calibrated parameter (γT or α for the PMEi parameterisations), which affects the further evolution of melt rates in a complicated way. For example, the thermal forcing for a given constant depth of 700 m is at all depths higher than the depth-dependent thermal forcing but results in less SLC for all scenarios apart from the Warm0 and Warm2 scenarios.

The quality of the PMEi parameterisation results, with regard to the coupled simulations, is linked to the degree of warming. The higher the thermal forcing, the poorer are the results. The SLC is systematically underestimated except for the coldest (Cold0) scenario, for which the SLC prediction is in agreement with the coupled results. In terms of melt rates, this parameterisation computes a different pattern compared to the other parameterisations. The melt rates are very low near the central grounding line and almost uniform downstream. This could explain why, compared to the other parameterisations, the prior pulse that they undergo is shorter in time and why after this pulse the melt rates drop down to much lower melt rates compared to others. After this pulse, the ice shelf is mostly composed of thick ice, and the low melt rates near the grounding line, where the ice is thicker, hamper the impact of melting on buttressing relative to the coupled and parameterised simulations. Surprisingly, the PMEi parameterisations are quite close to one another, regardless of the approach used to define effective grounding line and angle.

The box parameterisations are forced by the ocean properties at a constant depth, being either 500 or 700 m depths. Whatever the depth, the higher the number of boxes, the larger both the overall melting and the SLC in our experiments, which is enhanced for the Warmi scenarios compared to the Coldi scenarios. Note that during the melt pulse in the beginning, the order seems to be reversed and total melting decreases with the number of boxes. The optimal number of boxes for the Coldi scenarios is between two and five, while for the Warmi scenarios using five boxes results in a good agreement with the coupled simulations and seems to be the best trade-off within the box model, regardless of the given forcing depth. Note that using 700 m for the forcing depth gives pretty good results whatever the number of boxes, while using 500 m ends up in a larger spread in our experiments.

A rise of the thermocline (Warm3 scenario) does not affect the coupled simulations, likely because sea-floor ocean properties remain unchanged in this experiment. This emphasises the importance of sea-floor ocean properties for ice-shelf melting, and explains why the box model is closer to the coupled model when ocean properties are taken at 700 m depth.

5 Discussion

Parameterising sub-shelf melt rates in ice-sheet modelling is currently the only way to account for melting in large-ensemble or multi-millennium simulations of the Antarctic ice sheet (DeConto and Pollard2016), and even shorter term simulations applied to single Antarctic basins have been done on very few occasions and only very recently (Thoma et al.2015; Seroussi et al.2017). Our study suggests that parameterisations should be chosen with caution. To assess the capacity of the parameterisations to reproduce the ocean-induced melting and its effect on ice-sheet dynamics under a wide range of scenarios, we set-up a performance indicator (Fig. 9). We define it as the root-mean-square deviation (RMSD) in SLC of every parameterisation with respect to the average of coupled simulations on a given year. We choose to calculate this performance indicator at the 50th year of the simulations, for a significant part of the ice shelf is melted out by the parameterisations after this year.

Figure 9Performance of parameterisations compared to coupled simulations, calculated at the 50th year of simulations. (a) Root mean square deviation (RMSD) in SLC of every parameterised simulation with respect to the average of coupled simulations. (b) Difference between SLCs from parameterisations and coupled simulations for all the experiments. The grey shading is only to ease the comparisons between the parameterisations.


While the plume parameterisation is in pretty good agreement with coupled simulations for the cold forcings, it consistently underestimates both the melt rates and subsequent SLC for the warm forcings. Lazeroms et al. (2018) show melt rate patterns in good agreement with observations for the large ice shelves such as Ronne-Filchner and Ross. However, for smaller ice shelves such as the Pine Island and Thwaites glaciers, the patterns exhibit very strong melting near the calving front, and quite a uniform melting in the entire cavity. This is contradictory to observation-based estimates (Rignot et al.2013a; Dutrieux et al.2013) and to high-resolution ocean simulations (Dutrieux et al.2014) showing large melt rates near the grounding line that drop abruptly from a few kilometres downstream to almost zero near the calving front. In our plume parameterisation configuration, the melt rates are zero at the grounding line, close to zero nearby (not seen in the Lazeroms et al.2018, paper because the resolution is too coarse) and the strongest at the calving front. We suspect this is due to the empirical relationship used in Lazeroms et al. (2018) that relates melt rates to the depth difference between the effective grounding line point and the ice draft, which may wrongly place the melting–accretion point for small ice shelves as opposed to large ice shelves. The fact that the same ice-sheet response occurs regardless of the type of implementation supports this point.

The box parameterisation tends to give relatively good results regardless of the number of boxes or the near sea-floor depth at which the ocean properties are taken. Using five boxes seems to yield the best results. Reese et al. (2018a) found that increasing the number of boxes in a static cavity would converge to almost constant average melt rates above five boxes. In our study, increasing the number of boxes leads neither to convergence of the calibrated parameter nor to converging SLC during the prognostic simulations. The melting pattern has an effect on the ice-sheet dynamics, so even though convergence could be expected from the work of Reese et al. (2018a) for a static cavity, the ice-sheet response to the different patterns related to the various number of boxes could have suppressed the initial convergence.

A key issue in our implementations of the 1-D plume parameterisation might be in the use of deep ocean temperatures, which will lead to an overestimate of melting near the ice front. Our calibration procedure then scales back the melting near the grounding line and leads to an underestimate of the reduction in buttressing. The box model also uses the deep temperatures, but in that parameterisation heat is supplied to the overturning circulation in the grounding zone only, beyond which melt rates must fall as a result of the extraction of latent heat and the rise in the freezing point. Calibrating the heat transfer coefficient alters the balance between heat used to melt in the grounding zone and that advected downstream to melt elsewhere. Hence, the calibration redistributes the melting rather than just scaling a fixed melt pattern, and that may be the reason that the results compare quite well with those from the coupled model, especially when the parameterisation is used with five boxes.

Among the simple functions of thermal forcing, the two quadratic (local and non-local) functions are in good agreement with the coupled simulations. A non-local dependency leads to slightly better results. Taking the ocean properties at a varying depth gives better results. In that case, these two parameterisations are the only ones to capture the increased melting in coupled simulations after the initial adjustment phase in the Warm1 scenario. When these simple functions depend on constant depth ocean properties, deeper temperature and salinity inputs result in better agreement with coupled simulations.

We chose to calibrate the parameterisations using the same far-field ocean temperature and salinity constant profiles, which is different from the temperature and salinity scenarios used in the rest of the study. Such an approach is actually very selective but enables to distinguish between parameterisations that could be applied to real cases, because they adapt well to a change in ocean properties, from those that either need to be improved or discarded with regard to changing ocean conditions.

All parameterisations yield too large melt rates in thin ice areas and too small melt rates near the deepest parts around the grounding line. Even though our geometrical set-up is ideal, the distribution of thicknesses within the ice shelf are not far from reality, meaning that applying these parameterisations to real ice shelves would also induce too much thinning of initially thin floating ice. The studies of Jenkins (2016) and Jenkins et al. (2018) suggest that the basal slope of the ice shelf influences the mixing across the thermocline. Accounting for this effect in simple functions of thermal forcing may allow to redistribute more melting over the steep areas near grounding line and less melting over flat areas near calving fronts, thus decreasing the overmelting of thin floating ice.

The choice of a parameterisation for real applications may account for the local circulation in the ice-shelf cavity. Whether the circulation is horizontal or vertical may guide the choice of the dependence on thermal forcing being either a function of varying depth or taken at a constant depth. For instance, the circulation in the Amundsen sea embayment appears to be a mix between vertical overturning fed by incursions of CDW and horizontal barotropic flow generated by tides (Jourdain et al.2017, 2018). It should be noted that our study does not account for sea ice, which tends to limit the Ekman pumping due to wind stress and vertical mixing, nor for tides.

The spatial distribution of melt rates affect ice-shelf buttressing in a complicated way. Similar total melt rates distributed differently beneath the ice shelf is likely to induce distinct responses of the ice sheet (Reese et al.2018b; Gagliardini et al.2010). Conversely, different melting patterns can induce similar responses of the ice sheet if the integrated losses in buttressing happen to be well balanced from one another. This is illustrated in our simulations, for instance by the two types of quadratic functions of the thermal forcing that exhibit different patterns but lead to similar SLC. The study of Reese et al. (2018b) attributes an equal effect of bits of ice-shelf removal on ice-sheet dynamics in places where ice thicknesses can be very different. Removing floating ice near the deepest grounding lines or near ice rises can remove the same amount of buttressing and lead to similar SLC. Ice rises are generally found in shallow waters, and thus a parameterisation that computes overly large melt rates near this sensitive area may remove too much buttressing restraining the upstream ice sheet compared to coupled simulations.

An ocean–ice-sheet coupled model is needed as a reference to assess the melting parameterisations. Only an ocean model can convey the complexity of ocean physics to melting at the ice-shelf base, as opposed to parameterisations, and only an ice-sheet model can respond to a change in ice-shelf buttressing induced by changing melt rates. On the one hand, the ocean model NEMO was used to calculate the melt rates in the coupled framework. On the other hand, the ice sheet was simulated by the Elmer/Ice model using the SSA* approximation of the Stokes equations and a Schoof friction law at the ice–bed interface. Over the last decade, many ice-sheet and ocean models were developed, which motivated various model intercomparison projects to evaluate the caveats and assets of models and their physics with regard to ideal simulations (the MISMIP and MISMIP3D projects in Pattyn et al.2012, 2013, for ice-sheet models; the ISOMIP project in Holland et al.2003, for ice-shelf–ocean models; and the MISMIP+, ISOMIP+, and MISOMIP1 projects in Asay-Davis et al.2016, for ice sheet, ice-shelf–ocean, and ocean–ice-sheet coupled models). These intercomparison projects have highlighted differences between models that have not been accounted for in our study, even though we included an ensemble of coupled configurations to quantify uncertainties in the ocean model grid and physics. This present study will need to be pursued using other types of models and physics to further assess the robustness of our results.

Our study highlights the assets and caveat of sub-shelf melt parameterisations that can be constrained by the far-field ocean, some of which are used over a decade without thorough assessment. This work was performed with an idealised representation of a relatively small outlet glacier in West Antarctica and now needs to be extended to Antarctic realistic ocean–ice-sheet systems in order to improve sea level projections.

6 Conclusions

We compared a wide variety of sub-shelf melting parameterisations depending on oceanic properties to an ensemble of ocean–ice-sheet coupled simulations, using a new coupled model combining the ocean model NEMO and the ice-sheet model Elmer/Ice. Among the complex parameterisations that we assessed, representing melting through a 2-D emulation of a 1-D plume model gives good results for cold conditions (e.g. in the Ronne-Filchner cavity) but underestimates the melt rates and sea level contribution for warm conditions (e.g. in Pine Island glacier cavity). Given the high degree of complexity in the physics represented in the plume model, it is possible that calibrating more parameters could improve the validity of the scaling across multiple ice-shelf sizes. More work may also improve the way to extend the 1-D plume model to a realistic ice draft. The box parameterisation representing the vertical overturning in the cavity gives results relatively close to the coupled simulations, especially when used with five boxes. We showed that a linear parameterisation of thermal forcing is not able to represent ocean-induced melting beneath an ice shelf. Instead, a quadratic parameterisation of thermal forcing gives much better results, which are even improved for a local and non-local approach, as opposed to a fully local approach. Studies aiming at projecting the future contribution of Antarctica to sea level should take care about the choice of the melting parameterisation before providing predictions. We recommend validating the chosen parameterisation with regard to ocean–ice-sheet model coupled simulations within the specific environmental conditions and ice physics, although our results have to be taken carefully, until assessment based upon other models are produced.

Code availability

We used Elmer/Ice Version 8.3 at revision 6be9699, which is available at (last access: 18 April 2019), and NEMO-3.6 at revision 6402. The experimental protocol is composed of the coupling framework version 1.1, available at (Jourdain and Favier2019); the NEMO set-up version Feb-2019, available at (Jourdain2019a); and the Elmer/Ice set-up version 1.2, available at (Jourdain2019b).

Appendix A: Schoof friction law

The Glen's flow law relates deviatoric stresses τij to strain rates ε˙ij as follows:

(A1) τ i j = A - 1 / n ε ˙ e ( 1 - n ) / n ε ˙ i j ,

with A the fluidity parameter, ε˙e the second invariant of strain rates, and n Glen's exponent.

The Schoof friction law is written as in Brondex et al. (2017) and Brondex et al. (2018) as follows:

(A2) τ b = C s u b m 1 + C s C max N 1 / m u b m ,

with τb the basal friction, Cs a friction parameter, ub the basal velocity, Cmax Iken's bound parameter, N the effective pressure, and m the basal friction exponent.

The values of the parameters accounted for in Eqs. (A1) and (A2) are given in Table A1.

Table A1Parameters of Glen's flow law and Schoof friction law. n/a – not applicable

Download Print Version | Download XLSX

Appendix B: Sensitivity to the coupling period


Figure B1Mean cavity melt rate seen by Elmer/Ice for various coupling periods (a). Global mean sea level rise equivalent to the ice-mass loss simulated by Elmer/Ice for various coupling periods (b). The four simulations correspond to the IceOcean1r experiment of the standard MISOMIP protocol (Asay-Davis et al.2016).


Appendix C: Physical basis for the plume parameterisation empirical scaling

The plume parameterisation is derived empirically from the results of a full plume model (Jenkins2014) applied to a range of simple ice-shelf geometries and water properties. If the ice-shelf base is linear and the far-field ocean uniform, results for a wide range of ocean temperatures, ice-shelf basal slopes, and grounding-line depths, when appropriately scaled, collapse (within ±20 %) onto a universal melt rate curve (Jenkins2014). The plume parameterisation of Lazeroms et al. (2018) was created by fitting an 11th-order polynomial function to the universal curve.

When applying the parameterisation in practice, there are a number of issues to deal with: the ice-shelf basal slope will vary; the far-field ocean will be non-uniform; and for 2-D ice-shelf geometries, there is no unique grounding-line point. The first two are generic problems that arise from the simplifications that are required to allow the derivation of a universal melt rate curve. The latter arises when the 1-D parameterisation is implemented in 2-D.

The ice shelf basal slope θ enters the parameterisation through the following function:


The last of these terms scales the thermal driving in the plume as a fraction of the far-field thermal driving, while the first two scale the plume speed based on the balance between buoyancy and friction (first) and the dependence of the buoyancy on far-field thermal driving (second). Since the inertia of the plume is small, its speed rapidly adjusts to changing slope, and the first term of the above expression therefore represents a local balance between the upslope buoyancy force and frictional drag. The latter two terms, on the other hand, reflect the balance between entrainment and melting over the path of the plume, and so cannot be directly related to the local slope, if the slope is non-uniform. However, for low slopes the turbulent transfer of heat and momentum at the ice base tends to dominate over entrainment, giving the following:

(C2) g ( θ ) = sin θ C d 1 / 2 E 0 sin θ C d Γ TS .

Hence, for low slopes the thermal driving evolves along the plume path with a simple sin θ scaling that effectively makes it a function of the depth change between the grounding line and the point of interest. It does not matter if that path is short and steep with rapid entrainment, or long and gentle with slow entrainment; the net result is the same. The first term remains a local scaling, so when the parameterisation is applied to 1-D problems with varying slope, using the local slope to estimate g(θ) gives good results (Lazeroms et al.2018).

Figure C1Similar to Fig. 9 but to evaluate the use of the local gradient in PME3 (which gives PME5) and PME4 (which gives PME6) to calculate the effective angle instead of using the slope between the ice draft and the grounding line from which the plume starts. The calibration of PME5 and PME6 is done with α=0.34 and α=0.65, respectively (Sect. 3.2). The grey shading is only to ease the comparisons between the parameterisations.


An equivalent solution to the problem of non-uniform far-field properties is less obvious. Taking a depth-average to reflect the range of properties entrained into the plume is the option used in Lazeroms et al. (2018). However, entrainment is strongest where the basal slope is steepest, so we might expect deeper waters to contribute more to the plume. In this study, we use the temperature at the depth of the grounding line. Since the temperature profiles all have a thick isothermal layer at the seabed, taking an average over the depth range from which waters are entrained into the plume would probably yield similar results, because the isothermal layer would dominate.

Implementation of the plume parameterisation in 2-D is a more complex problem, to which there are many possible solutions. The procedure implemented by Lazeroms et al. (2018) was effectively an average of 1-D implementations along whichever of 16 prescribed directions represented valid plume paths. For each valid direction the grounding-line depth and the local slope in that direction were used to scale the melt rates. In PME1 we implemented that procedure as closely as possible, given the unstructured model grid (Appendix D). However, following the above reasoning, it could be argued that the path followed by the plume should not matter. Only the depth change from the grounding line to the point of interest should influence the results, and the plume should flow locally up the steepest slope, i.e. parallel to the gradient vector. Using the magnitude of the local gradient vector as the slope scale (PME2) decreases the checkerboard noise in the PME1 melt rates (Fig. 4) but does not greatly influence the results. In PME3 and PME4, we adopted two procedures for picking a unique grounding line point for each grid point, rather than using an average of many. In each case, we scaled the melt rate using the depth of the grounding-line point and the mean slope along a straight line connecting the grounding line point and the grid point. Following the earlier reasoning, using the magnitude of the local gradient vector is equally valid. Results using that alternative are shown in Fig. C1 but differ little from those presented in Fig. 9.

Appendix D: Implementations of the plume parameterisation in Elmer/Ice

The plume parameterisation was originally implemented for regular grids. We adapt the method used to calculate the effective grounding line depth and effective angle as defined in Lazeroms et al. (2018) to the unstructured grids used in Elmer/Ice. Here, we describe the adapted method and the alternative implementations also discussed in the paper.

D1 PME1 (plume model emulator 1)

In the original algorithm published in Lazeroms et al. (2018), the melt rates at a draft point are calculated by considering the effective grounding line depth, which is calculated by searching in the 16 grid directions equally distributed around the draft point, and starting from it, insofar as those directions are valid. The 16 directions follow the grid points as shown in Fig. 3 of Lazeroms et al. (2018). A direction is valid if (i) the local slope in this direction is negative and (ii) the first grounded point met in this direction is deeper than the draft point in this direction. The effective grounding-line depth of the draft point (i,j) of the regular grid zgl(i,j) is calculated using Eq. (13a) of Lazeroms et al. (2018), and its effective angle θ(i,j) is calculated using Eq. (13b) by considering the local slopes in the valid directions. These equations are recalled here below:


with Nij the number of valid directions, zn(i,j) the grounding line depth, and sn(i,j) the local slope in the direction of the grounding line.

Instead of grid directions, we consider directional triangles that are angularly equally distributed around the draft point. Analogously to the original criterion to find valid directions, the criterion to make a cone valid is based on (i) the average of the local angles of all the directions connecting the draft point to the grounding line points included in the cone and (ii) the average of these grounding line point depths (Fig. D1a). The simulations were all done using 64 triangles around the draft point, which enables a rather smooth melting pattern compared to using 16 triangles (analogously to Lazeroms et al.2018, for directions).


In the algorithm published in Appendix 2 of Lazeroms et al. (2018), i.e. the discussion version of Lazeroms et al. (2018), the criterion to make a direction valid is the same as the first algorithm, but the computation of the effective grounding-line depth and angle is slightly different, taking for instance the local gradient instead of the local slope in each direction to calculate the effective angle.

We calculate the effective grounding line depth and angle as it is in Eqs. (B1) and (B2) of Lazeroms et al. (2018) but using the search for valid directions as explained in Sect. D1 (Fig. D1a), also using 64 triangles.


In this algorithm, we simply take the deepest grounding line (which is located in the central flow line) to calculate the effective grounding line depth, and the effective angle is the slope between this grounding-line point and the draft point (Fig. D1b).


This algorithm accounts for the asymmetry resulting from the Coriolis effect, although in a very crude manner. The effective grounding-line depth is found by starting from the closest grounding-line point and looking for a deeper contiguous grounding-line point in the anti-clockwise direction as long as the grounding line deepens. Two examples of this algorithm are shown in (Fig. D1c) for two draft points in the left and the right side of the cavity, respectively. The effective angle is calculated as in Appendix D3.


This implementation is similar to PME3 but the effective angle is calculated from the ice-draft local gradient. The results are not given in the main article but compared to the other plume parameterisations in Fig. C1.


This implementation is similar to PME4 but the effective angle is calculated from the ice-draft local gradient. The results are not given in the main article but compared to the other plume parameterisations in Fig. C1.

Figure D1Computing the effective grounding-line depth and angle for the four plume parameterisation implementations. Panel (a) illustrates both the published implementation PME1 (Lazeroms et al.2018) and the one appearing in the corresponding discussion article PME2 (Lazeroms et al.2018), here with 12 directions (vs. 64 triangles used in the present paper). PME1 and PME2 differ by the calculation of the effective slope, which is not depicted here. Panel (b) illustrates the simple implementation PME3 and panel (c) illustrates the asymmetric implementation PME4.


Appendix E: CMIP5 temperature anomalies in the Amundsen Sea under the RCP8.5 scenario

Figure E1Temperature anomaly (2080–2100 mean minus 1989–2009 mean) in the Amundsen Sea (76–69 S, 128–90 W) from 33 CMIP5 models in the RCP85 scenario. Continental shelf temperatures (a) are averaged over the area where the sea floor is shallower than 1500 m, while offshore temperatures (b) are averaged over the rest of the domain. The numbers for individual CMIP5 models and the multi-model mean (MMM) indicate the mean ocean warming in the 500–800 m layer.


Appendix F: MISOMIP original calibration below 300 m depth

Figure F1Same as Fig. 9 but for an initial calibration based on averaged melt rates below 300 m depth.


Author contributions

LF and NCJ designed the experiments and implemented the melt rate parameterisations in Elmer/Ice. FGC implemented the SSA* in Elmer/Ice. NM developed the coupling interface. PM developed the ice-shelf routines for evolving geometries. LF led the writing, and all authors contributed to the writing and discussion of ideas.

Competing interests

The authors declare that they have no conflict of interest.


The simulations were performed using HPC resources from GENCI-CINES (Grant 2018, project A0040106066). We thank two anonymous reviewers for their valuable and constructive comments, which helped to improve the clarity of the paper.

Financial support

This work was funded by the French National Research Agency (ANR) through the TROIS-AS (ANR-15-CE01-0005-01) and the ANR-15-IDEX-02 projects.

Review statement

This paper was edited by Alexander Robel and reviewed by two anonymous referees.


Asay-Davis, X. S., Cornford, S. L., Durand, G., Galton-Fenzi, B. K., Gladstone, R. M., Gudmundsson, G. H., Hattermann, T., Holland, D. M., Holland, D., Holland, P. R., Martin, D. F., Mathiot, P., Pattyn, F., and Seroussi, H.: Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects: MISMIP v. 3 (MISMIP +), ISOMIP v. 2 (ISOMIP +) and MISOMIP v. 1 (MISOMIP1), Geosci. Model Dev., 9, 2471–2497,, 2016. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s

Asay-Davis, X. S., Jourdain, N. C., and Nakayama, Y.: Developments in Simulating and Parameterizing Interactions Between the Southern Ocean and the Antarctic Ice Sheet, Curr. Clim. Chang. Reports, 3, 316–329,, 2017. a, b

Beckmann, A. and Goosse, H.: A parameterization of ice shelf-ocean interaction for climate models, Ocean Model., 5, 157–170,, 2003. a, b

Brondex, J., Gagliardini, O., Gillet-Chaulet, F., and Durand, G.: Sensitivity of grounding line dynamics to the choice of the friction law, J. Glaciol., 63, 854–866,, 2017. a, b

Brondex, J., Gillet-chaulet, F., and Gagliardini, O.: Sensitivity of centennial mass loss projections of the Amundsen basin to the friction law, The Cryosphere, 13, 177–195,, 2019. a, b

Cornford, S. L., Martin, D. F., Payne, A. J., Ng, E. G., Le Brocq, A. M., Gladstone, R. M., Edwards, T. L., Shannon, S. R., Agosta, C., van den Broeke, M. R., Hellmer, H. H., Krinner, G., Ligtenberg, S. R. M., Timmermann, R., and Vaughan, D. G.: Century-scale simulations of the response of the West Antarctic Ice Sheet to a warming climate, The Cryosphere, 9, 1579–1600,, 2015. a, b

DeConto, R. M. and Pollard, D.: Contribution of Antarctica to past and future sea-level rise, Nature, 531, 591–597,, 2016. a, b

Depoorter, M. A., Bamber, J. L., Griggs, J. A., Lenaerts, J. T., Ligtenberg, S. R., Van Den Broeke, M. R., and Moholdt, G.: Calving fluxes and basal melt rates of Antarctic ice shelves, Nature, 502, 89–92,, 2013. a

De Rydt, J. and Gudmundsson, G. H.: Coupled ice shelf-ocean modeling and complex grounding line retreat from a seabed ridge, J. Geophys. Res.-Earth, 121, 865–880,, 2016. a, b

De Rydt, J., Holland, P. R., Dutrieux, P., and Jenkins, A.: Geometric and oceanographic controls on melting beneath Pine Island Glacier, J. Geophys. Res.-Oceans, 119, 2420–2438,, 2014. a, b

Donat-Magnin, M., Jourdain, N. C., Spence, P., Le Sommer, J., Gallée, H., and Durand, G.: Ice-Shelf Melt Response to Changing Winds and Glacier Dynamics in the Amundsen Sea Sector, Antarctica, J. Geophys. Res.-Oceans, 122, 2090–2107,, 2017. a, b, c

Durand, G., Gagliardini, O., De Fleurian, B., Zwinger, T., and Le Meur, E.: Marine ice sheet dynamics: Hysteresis and neutral equilibrium, J. Geophys. Res.-Sol. Ea., 114, 1–10,, 2009. a

Dutrieux, P., Vaughan, D. G., Corr, H. F. J., Jenkins, A., Holland, P. R., Joughin, I., and Fleming, A. H.: Pine Island glacier ice shelf melt distributed at kilometre scales, The Cryosphere, 7, 1543–1555,, 2013. a

Dutrieux, P., De Rydt, J., Jenkins, A., Holland, P. R., Ha, H. K., Lee, S. H., Steig, E. J., Ding, Q., Abrahamsen, E. P., and Schröder, M.: Strong sensitivity of pine Island ice-shelf melting to climatic variability, Science, 343, 174–178,, 2014. a, b, c, d

Favier, L., Durand, G., Cornford, S. L., Gudmundsson, G. H., Gagliardini, O., Gillet-Chaulet, F., Zwinger, T., Payne, A. J., and Le Brocq, A. M.: Retreat of Pine Island Glacier controlled by marine ice-sheet instability, Nat. Clim. Change, 4, 117–121,, 2014. a, b, c

Favier, L., Pattyn, F., Berger, S., and Drews, R.: Dynamic influence of pinning points on marine ice-sheet stability: a numerical study in Dronning Maud Land, East Antarctica, The Cryosphere, 10, 2623–2635,, 2016. a

Fretwell, P., Pritchard, H. D., Vaughan, D. G., Bamber, J. L., Barrand, N. E., Bell, R., Bianchi, C., Bingham, R. G., Blankenship, D. D., Casassa, G., Catania, G., Callens, D., Conway, H., Cook, A. J., Corr, H. F. J., Damaske, D., Damm, V., Ferraccioli, F., Forsberg, R., Fujita, S., Gim, Y., Gogineni, P., Griggs, J. A., Hindmarsh, R. C. A., Holmlund, P., Holt, J. W., Jacobel, R. W., Jenkins, A., Jokat, W., Jordan, T., King, E. C., Kohler, J., Krabill, W., Riger-Kusk, M., Langley, K. A., Leitchenkov, G., Leuschen, C., Luyendyk, B. P., Matsuoka, K., Mouginot, J., Nitsche, F. O., Nogi, Y., Nost, O. A., Popov, S. V., Rignot, E., Rippin, D. M., Rivera, A., Roberts, J., Ross, N., Siegert, M. J., Smith, A. M., Steinhage, D., Studinger, M., Sun, B., Tinto, B. K., Welch, B. C., Wilson, D., Young, D. A., Xiangbin, C., and Zirizzotti, A.: Bedmap2: improved ice bed, surface and thickness datasets for Antarctica, The Cryosphere, 7, 375–393,, 2013. a

Gagliardini, O., Cohen, D., Råback, P., and Zwinger, T.: Finite-element modeling of subglacial cavities and related friction law, J. Geophys. Res.-Earth, 112, 1–11,, 2007. a

Gagliardini, O., Durand, G., Zwinger, T., Hindmarsh, R. C., and Le Meur, E.: Coupling of ice-shelf melting and buttressing is a key process in ice-sheets dynamics, Geophys. Res. Lett., 37, 1–5,, 2010. a, b, c

Gagliardini, O., Zwinger, T., Gillet-Chaulet, F., Durand, G., Favier, L., de Fleurian, B., Greve, R., Malinen, M., Martín, C., Råback, P., Ruokolainen, J., Sacchettini, M., Schäfer, M., Seddik, H., and Thies, J.: Capabilities and performance of Elmer/Ice, a new-generation ice sheet model, Geosci. Model Dev., 6, 1299–1318,, 2013. a

Goldberg, D. N., Snow, K., Holland, P., Jordan, J. R., Campin, J. M., Heimbach, P., Arthern, R., and Jenkins, A.: Representing grounding line migration in synchronous coupling between a marine ice sheet model and a z-coordinate ocean model, Ocean Model., 125, 45–60,, 2018. a

Grosfeld, K., Gerdes, R., and Determann, J.: Thermohaline circulation and interaction between ice shelf cavities and the adjacent open ocean, J. Geophys. Res., 102, 595–610, 1997. a

Gudmundsson, G. H.: Ice-shelf buttressing and the stability of marine ice sheets, The Cryosphere, 7, 647–655,, 2013. a

Gudmundsson, G. H., Krug, J., Durand, G., Favier, L., and Gagliardini, O.: The stability of grounding lines on retrograde slopes, The Cryosphere, 6, 1497–1505,, 2012. a

Gwyther, D. E., O'Kane, T. J., Galton-Fenzi, B. K., Monselesan, D. P., and Greenbaum, J. S.: Intrinsic processes drive variability in basal melting of the Totten Glacier Ice Shelf, Nat. Commun., 9, 3141,, 2018. a

Haseloff, M. and Sergienko, O. V.: The effect of buttressing on grounding line dynamics, J. Glaciol., 64, 417–431,, 2018. a

Hellmer, H. and Olbers, D.: A two-dimensional model for the thermohaline circulation under an ice shelf, Antarct. Sci., 1, 325–336,, 1989. a

Hellmer, H. H., Kauker, F., Timmermann, R., Determann, J., and Rae, J.: Twenty-first-century warming of a large Antarctic ice-shelf cavity by a redirected coastal current, Nature, 485, 225–228,, 2012. a

Hellmer, H. H., Kauker, F., Timmermann, R., and Hattermann, T.: The fate of the Southern Weddell sea continental shelf in a warming climate, J. Climate, 30, 4337–4350,, 2017. a

Holland, D., Hunter, S. J., Grosfeld, K., Hellmer, H. H., Jenkins, A., Morales-Maqueda, M. A., Hemer, M., Williams, M., Klinck, J. M., and Dinniman, M. S.: The Ice Shelf -– Ocean Model Intercomparison Project (ISOMIP),, Eos Trans. AGU, 84, Abstr. C41A–05, Fall Meet. Suppl., 2003. a

Holland, D. M. and Jenkins, A.: Modeling Thermodynamic Ice–Ocean Interactions at the Base of an Ice Shelf, J. Phys. Oceanogr., 29, 1787–1800,<1787:MTIOIA>2.0.CO;2, 1999. a

Holland, P. R.: The Transient Response of Ice Shelf Melting to Ocean Change, J. Phys. Oceanogr., 47, 2101–2114,, 2017. a

Holland, P. R., Jenkins, A., and Holland, D. M.: The response of Ice shelf basal melting to variations in ocean temperature, J. Climate, 21, 2558–2572,, 2008. a, b

Jacobs, S. S., Jenkins, A., Giulivi, C. F., and Dutrieux, P.: Stronger ocean circulation and increased melting under Pine Island Glacier ice shelf, Nat. Geosci., 4, 519–523,, 2011. a

Jenkins, A.: A one-dimensional model of ice shelf-ocean interaction, J. Geophys. Res., 96, 671–677, 1991. a, b

Jenkins, A.: Scaling laws for the melt rate and overturning circulation beneath ice shelves derived from simple plume theory, EGU Gen. Assem. Conf. Abstr. 16, 13755, 2014. a, b, c, d

Jenkins, A.: A Simple Model of the Ice Shelf – Ocean Boundary Layer and Current, J. Phys. Oceanogr., 46, 1785–1803,, 2016. a

Jenkins, A., Dutrieux, P., Jacobs, S. S., McPhail, S. D., Perrett, J. R., Webb, A. T., and White, D.: Observations beneath Pine Island Glacier in West-Antarctica and implications for its retreat, Nat. Geosci., 3, 468–472,, 2010. a

Jenkins, A., Shoosmith, D., Dutrieux, P., Jacobs, S., Kim, T. W., Lee, S. H., Ha, H. K., and Stammerjohn, S.: West Antarctic Ice Sheet retreat in the Amundsen Sea driven by decadal oceanic variability, Nat. Geosci., 11, 733-738,, 2018. a, b, c

Jordan, J. R., Holland, P. R., Goldberg, D., Snow, K., Arthern, R., Campin, J. M., Heimbach, P., and Jenkins, A.: Ocean-Forced Ice-Shelf Thinning in a Synchronously Coupled Ice-Ocean Model, J. Geophys. Res.-Oceans, 123, 864–882,, 2018. a, b

Joughin, I., Smith, B. E., and Medley, B.: Marine Ice Sheet Collapse Potentially Under Way for the Thwaites Glacier Basin, West Antarctica, Science, 344, 735–738,, 2014. a, b, c

Jourdain, N. C.: nicojourdain/NEMO_PARAMS_SIMU: NEMO material (including Favier et al. 2019) (Version Feb-2019), Zenodo,, 2019a. a

Jourdain, N. C.: nicojourdain/ElmerIce_Experiments: Initial Version (Version 1.2), Zenodo,, 2019b. a

Jourdain, N. C. and Favier, L.: nicojourdain/CPL_NEMO_ELMER: NEMO – Elmer/Ice coupler (v1.1) (Version v1.1), Zenodo,, 2019. a

Jourdain, N. C., Mathiot, P., Merino, N., Durand, G., Le Sommer, J., Spence, P., Dutrieux, P., and Madec, G.: Ocean circulation and sea-ice thinning induced bymelting ice shelves in the Amundsen Sea, J. Geophys. Res.-Oceans, 122, 2550–2573,, 2017. a, b, c

Jourdain, N. C., Molines, J.-M., Le Sommer, J., Mathiot, P., Chanut, J., de Lavergne, C., and Madec, G.: Simulating or prescribing the influence of tides on the Amundsen Sea ice shelves, Ocean Model., 33, 44–55,, 2018. a, b

Konrad, H., Gilbert, L., Cornford, S. L., Payne, A., Hogg, A., Muir, A., and Shepherd, A.: Uneven onset and pace of ice-dynamical imbalance in the Amundsen Sea Embayment, West Antarctica, Geophys. Res. Lett., 44, 910–918,, 2017. a

Konrad, H., Shepherd, A., Gilbert, L., Hogg, A. E., McMillan, M., Muir, A., and Slater, T.: Net retreat of Antarctic glacier grounding lines, Nat. Geosci., 11, 258–262,, 2018. a

Lazeroms, W. M. J., Jenkins, A., Gudmundsson, G. H., and van de Wal, R. S. W.: Modelling present-day basal melt rates for Antarctic ice shelvesusing a parametrization of buoyant meltwater plumes, The Cryosphere Discuss.,, 2017. a, b, c

Lazeroms, W. M. J., Jenkins, A., Gudmundsson, G. H., and van de Wal, R. S. W.: Modelling present-day basal melt rates for Antarctic ice shelves using a parametrization of buoyant meltwater plumes, The Cryosphere, 12, 49–70,, 2018. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v

Losch, M.: Modeling ice shelf cavities in a z coordinate ocean general circulation model, J. Geophys. Res. Ocean., 113, 1–15,, 2008. a, b

Madec, G. and NEMO-team: NEMO ocean engine, version 3.6 stable, Note du Pôle de modélisation de l'Institut Pierre-Simon Laplace, SSN No 1288-1619, Technical Report, IPSL, France, 2016. a

Mathiot, P., Jenkins, A., Harris, C., and Madec, G.: Explicit representation and parametrised impacts of under ice shelf seas in the z coordinate ocean model NEMO 3.6, Geosci. Model Dev., 10, 2849–2874,, 2017 a, b, c

Mengel, M. and Levermann, A.: Ice plug prevents irreversible discharge from east Antarctica, Nat. Clim. Change, 4, 451–455,, 2014. a

Mercer, J. H.: West Antarctic ice sheet and CO2 greenhouse effect: a threat of disaster, Nature, 271, 321–325,, 1978. a

Mouginot, J., Rignot, E., and Scheuchl, B.: Sustained increase in ice discharge from the Amundsen Sea Embayment, West Antarctica, from 1973 to 2013, 41, 1576–1584,, 2014. a

Nias, I. J., Cornford, S. L., and Payne, A. J.: Contrasting the Modelled sensitivity of the Amundsen Sea Embayment ice streams, J. Glaciol., 62, 552–562,, 2016. a

Olbers, D. and Hellmer, H.: A box model of circulation and melting in ice shelf caverns, Ocean Dynam., 60, 141–153,, 2010. a, b

Paolo, F. S., Fricker, H. A., and Padman, L.: Volume loss from Antarctic ice shelves is accelerating, Science, 348, 327–331,, 2015. a

Pattyn, F., Schoof, C., Perichon, L., Hindmarsh, R. C. A., Bueler, E., de Fleurian, B., Durand, G., Gagliardini, O., Gladstone, R., Goldberg, D., Gudmundsson, G. H., Huybrechts, P., Lee, V., Nick, F. M., Payne, A. J., Pollard, D., Rybak, O., Saito, F., and Vieli, A.: Results of the Marine Ice Sheet Model Intercomparison Project, MISMIP, The Cryosphere, 6, 573–588,, 2012. a

Pattyn, F., Perichon, L., Durand, G., Favier, L., Gagliardini, O., Hindmarsh, R. C., Zwinger, T., Albrecht, T., Cornford, S., Docquier, D., Fürst, J. J., Goldberg, D., Gudmundsson, G. H., Humbert, A., Hütten, M., Huybrechts, P., Jouvet, G., Kleiner, T., Larour, E., Martin, D., Morlighem, M., Payne, A. J., Pollard, D., Rückamp, M., Rybak, O., Seroussi, H., Thoma, M., and Wilkens, N.: Grounding-line migration in plan-view marine ice-sheet models: Results of the ice2sea MISMIP3d intercomparison, J. Glaciol., 59, 410–422,, 2013. a

Pattyn, F., Ritz, C., Hanna, E., Asay-Davis, X., DeConto, R., Durand, G., Favier, L., Fettweis, X., Goelzer, H., Golledge, N. R., Munneke, P. K., Lenaerts, J. T. M., Nowicki, S., Payne, A. J., Robinson, A., Seroussi, H., and Tr, M.: The Greenland and Antarctic ice sheets under 1.5 C global warming, Nat. Clim. Change, 8, 1053–1061,, 2018. a

Reese, R., Albrecht, T., Mengel, M., Asay-Davis, X., and Winkelmann, R.: Antarctic sub-shelf melt rates via PICO, The Cryosphere, 12, 1969–1985,, 2018a. a, b, c, d, e, f, g, h

Reese, R., Gudmundsson, G. H., Levermann, A., and Winkelmann, R.: The far reach of ice-shelf thinning in Antarctica, Nat. Clim. Change, 8, 53–57,, 2018b. a, b, c

Rignot, E., Jacobs, S., Mouginot, J., and Scheuchl, B.: Ice-shelf melting around antarctica, Science, 341, 266–270,, 2013a. a

Rignot, E., Jacobs, S., Mouginot, J., and Scheuchl, B.: Ice-shelf melting around antarctica, Science, 341, 266–270,, 2013b. a

Rignot, E., Mouginot, J., Morlighem, M., Seroussi, H., and Scheuchl, B.: Widespread, rapid grounding line retreat of Pine Island, Thwaites, Smith, and Kohler glaciers, West Antarctica, from 1992 to 2011, Geophys. Res. Lett., 41, 3502–3509,, 2014. a

Schoof, C.: The effect of cavitation on glacier sliding, Proc. R. Soc. A, 461, 609–627,, 2005. a

Schoof, C.: Ice sheet grounding line dynamics: Steady states, stability, and hysteresis, J. Geophys. Res.-Earth, 112, 1–19,, 2007. a

Schoof, C. and Hindmarsh, R. C.: Thin-film flows with wall slip: An asymptotic analysis of higher order glacier flow models, Q. J. Mech. Appl. Math., 63, 73–114,, 2010. a

Seroussi, H., Morlighem, M., Larour, E., Rignot, E., and Khazendar, A.: Hydrostatic grounding line parameterization in ice sheet models, The Cryosphere, 8, 2075–2087,, 2014. a, b

Seroussi, H., Nakayama, Y., Larour, E., Menemenlis, D., Morlighem, M., Rignot, E., and Khazendar, A.: Continued retreat of Thwaites Glacier, West Antarctica, controlled by bed topography and ocean circulation, Geophys. Res. Lett., 44, 6191–6199,, 2017. a, b, c

Shepherd, A.: Mass balance of the Antarctic ice sheet, Nature, 364, 1627–1635,, 2018. a

Spence, P., Griffies, S. M., England, M. H., Hogg, A. M., Saenko, O. A., and Jourdain, N. C.: Rapid subsurface warming and circulation changes of Antarctic coastal waters by poleward shifting winds, Geophys. Res. Lett., 41, 4601–4610, 2014. a

Thoma, M., Determann, J., Grosfeld, K., Goeller, S., and Hellmer, H. H.: Future sea-level rise due to projected ocean warming beneath the Filchner Ronne Ice Shelf: A coupled model study, Earth Planet. Sc. Lett., 431, 217–224,, 2015. a, b

Thomas, R. H. and Bentley, C. R.: A model for Holocene retreat of the West Antarctic Ice Sheet, Quaternary Res., 10, 150–170,, 1978. a

Treguier, A. M., Deshayes, J., Le Sommer, J., Lique, C., Madec, G., Penduff, T., Molines, J.-M., Barnier, B., Bourdalle-Badie, R., and Talandier, C.: Meridional transport of salt in the global ocean from an eddy-resolving model, Ocean Sci., 10, 243–255,, 2014.  a, b

Weertman, J.: Stability of the Junction of an Ice Sheet and an Ice Shelf, J. Glaciol., 13, 3–11,, 1974. a

Short summary
The melting at the base of floating ice shelves is the main driver of the Antarctic ice sheet current retreat. Here, we use an ideal set-up to assess a wide range of melting parameterisations depending on oceanic properties with regard to a new ocean–ice-sheet coupled model, published here for the first time. A parameterisation that depends quadratically on thermal forcing in both a local and a non-local way yields the best results and needs to be further assessed with more realistic set-ups.