Assessment of sub-shelf melting parameterisations using the ocean-ice-sheet coupled model

. 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-ﬁeld 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 deﬁne six idealised one-century scenarios for the far-ﬁeld 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 ﬁve 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.


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;Shepherd, 2018). 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 oceanice 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 , ice-surface lowering (Konrad et al., 2017), retreating grounding lines 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 insta-Published by Copernicus Publications on behalf of the European Geosciences Union. bility (MISI), which states that an ice sheet starting to retreat over a retrograde bed slope keeps retreating until the slope becomes prograde (Mercer, 1978;Thomas and Bentley, 1978;Weertman, 1974;Schoof, 2007;Durand et al., 2009). Confined ice shelves resist horizontal shearing and potentially stabilise an ice sheet undergoing MISI Haseloff and Sergienko, 2018). 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 (Losch, 2008). The existence of strong feedbacks between the cavity geometry, melt rates, and the ocean circulation 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 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 Gudmundsson, 2016).
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 Goosse, 2003;Favier et al., 2016) or a quadratic dependency (e.g. DeConto and Pollard, 2016). 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 (Jenkins, 1991) 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-icesheet coupled simulations. Overall, the MISOMIP (Asay- Davis et al., 2016) framework is used to perform 138 onecentury 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.

The ice-sheet model, Elmer/Ice
We perform the ice-sheet simulations with the finite-element ice-sheet model Elmer/Ice . 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 Schooflike 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.

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.

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-team, 2016). NEMO solves the prognostic equations for the ocean temperature, salinity, and velocities and includes ice-shelf cavities . The sub-shelf melting is parameterised through the so-called "three equations" representing (1) the heat balance at the iceocean 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 Olbers, 1989;Holland and Jenkins, 1999;Losch, 2008;Jenkins et al., 2010). In this parameterisation, we assume a constant top-boundary-layer (TBL) thickness along the ice-shelf draft , and we use a velocity-dependent formulation in which the heat exchange velocity is defined as follows: where u TBL the TBL-averaged velocity resolved by NEMO, T is the non-dimensional heat exchange coefficient, C d the non-dimensional drag coefficient and u tide is a uniform background velocity representing the main effect of tides on iceshelf melting (Jourdain et al., 2018). The values of T , C d , and u tide 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.

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).
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).

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).

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 standalone 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, Holland, 2017). 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.

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 iceocean boundary layer can be determined directly from farfield 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: 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, c po the specific heat capacity of the x is the horizontal resolution, T CPL is the ocean-icesheet 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 u tide are defined in Eq. (1), and the salt exchange coefficient S is taken as T /35. Also defined in Eq. (1) is the drag ID Name  Table 2). The melting-freezing point T f at the interface between the ocean and the ice-shelf basal surface is defined as follows: The practical salinity S o and the potential temperature T o are taken from the far-field ocean as detailed below in this section; z b 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 ) 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: 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: 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 T o and S o 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 T o = T o (z) and S o = S o (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, T o and S o 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 T f is therefore calculated with either S o (z) in the first option, or S o (500) or S o (700) in the second option (in a consistent way with T o ), 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.

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: 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 (Jenkins, 2014). The melt rates are given by the following: where PME means plume model emulator; M o is an overall scaling parameter; g(θ ) is a function of the ice-shelf basal Elmer/Ice grid resolution 500 m at the grounding line to 4 km away NEMO grid resolution 1 or 2 km in the horizontal, 10 or 20 m in the vertical (Table 1) Coupling period between 2 and 6 months ( Table 1) Plume parameterisation PME 1 published implementation Appendix D α = 0.75 (Lazeroms et al., 2018) PME 2 alternative implementation (Appendix B in the discussion paper) Appendix D α = 0.53 PME 3 simple implementation Appendix D α = 0.32 PME 4 asymmetric implementation Appendix D α = 0.63 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 (Jenkins, 2014). α 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 pa-rameterisation 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).

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 icesheet 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.

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. 3; Asay-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 m t 8.5± m a −1 . These four spin-ups will thus be used as initial states for subsequent coupled simulations.
Then, m t is used as a target for stand-alone icesheet 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 MI-SOMIP 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 onecentury scenarios described in Sect. 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: -Warm 0 resembles the present-day typical Amundsen Sea conditions . There is no temporal change of temperature and salinity profiles. -Warm 1 starts from the Warm 0 profile and then the temperature uniformly increases by 1 • C per century. The salinity profile is constant in time.
-Warm 2 is similar to Warm 1 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.
-Warm 3 starts from the Warm 0 profile and undergoes a 200 m uplift of both the thermocline and the halocline.
-Cold 0 resembles a cold cavity such as beneath the Ronne-Filchner ice shelves. There is no temporal change of temperature and salinity profiles.
-Cold 1 starts from the Cold 0 profile and then warms to reach a warm cavity state within a century. The salinity is also increased.
These profiles are slightly more realistic than in MIS-OMIP. They all include a thermocline, because its importance in ice-shelf melting has been pointed out by previous studies (e.g. De . The Warm 0 profile corresponds to a linear representation of the average hydrographic profiles measured in front of Pine Island glacier . By contrast, the Cold 0 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 Warm 1 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 Warm 2 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 Warm 3 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 Cold 1 scenario is an idealised representation of the ocean tipping point described by Hellmer et al. (2012Hellmer et al. ( , 2017, in which the Ronne-Filchner cavities switch from a cold to a warm state. The salinity profile is unchanged throughout Warm 0 , Warm 1 , and Warm 2 and is sufficiently stratified to keep a stable density profile. In the Warm 3 scenario, the halocline is lifted together with the thermocline to mimic an Ekmandriven uplift of the pycnocline, and in Cold 1 , 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 standalone 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.

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 param- 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 M quad 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. M lin , M quad , 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 M quad parameterisation, thinner ice being less melted and thicker ice being more melted, compared to M lin and M + . The M lin and M + patterns are similar by construction because the melting average is driven by the (T o − T f ) term, which appears only once in the two respective formulations. However, the respective calibrations are different (Table 3) because of the term T o − T f 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 PME 1 , PME 2 , and PME 3 on the one hand and PME 4 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 Geosci. Model Dev., 12, 2255-2283, 2019 www.geosci-model-dev.net/12/2255/2019/ line, and not from the sides where the basal surface is higher than the draft point (PME 1 in Fig. 4). Farther away, PME 1 and PME 2 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 PME 3 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 PME 1 and PME 2 for which the plumes can also come from less deep parts of the cavity and mitigate the melt rates.
Similarly to the M lin , M quad , 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.

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 Cold 1 scenario and only shows for the coupled simulations and not for the parameterised simulation for the Cold 0 scenario. For the Warm i 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. In the Cold 0 scenario, almost all parameterised and coupled simulations produce constant melt rates. Only the M lin parameterisations produce very high melt rates at the start and decrease monotonically afterwards. In the other constant scenario, which is Warm 0 , 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. M lin _700 forced by Warm 3 ). Finally, the Warm i scenarios end up with between 40 and 175 Gt a −1 of melting and the Cold 1 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 Warm i scenar-ios, 2 to 4 mm for the Cold 1 scenario and 0.5 to 3 mm for the Cold 0 scenario (Figs. 7, 8).
For the Warm i 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 Cold i 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 Cold 0 and Warm 2 scenarios, where it is about ± 20 %, respectively. A larger spread of about ± 30 % for the Warm i scenarios, and about ± 50 % and ± 100 % for the Cold 1 and Cold 0 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 M lin parameterisations tend to largely overestimate the melt rates for the Cold i scenarios and underestimate them for the Warm i 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 M quad parameterisations give melting in fair agreement with coupled results for the Cold i scenarios. For the Warm i scenarios, the tendency is a slight underestimation of SLC using the M quad and M quad _700 parameterisations, and a larger underestimation using M quad _500. Compared to the M lin parameterisations, it behaves much better and for a larger range of scenarios. All the M quad parameterisations behave quite well when confronted with a rise in the thermocline (Warm 3 scenario), apart from M quad _700, which slightly underestimates SLC.
The M + parameterisation results are almost as close to the coupled simulations as the M quad parameterisations for the Cold i scenarios, and closest for the Warm i 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 M quad _500 and M quad _700 parameterisations, thus slightly underestimating SLC.
Forcing a parameterisation by the far-field depthdependent 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 PME i 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 Warm 0 and Warm 2 scenarios.
The quality of the PME i 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 (Cold 0 ) 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 uni-form 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 PME i 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 Warm i scenarios compared to the Cold i 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 Cold i scenarios is between two and five, while for the Warm i 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 (Warm 3 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 iceshelf melting, and explains why the box model is closer to the coupled model when ocean properties are taken at 700 m depth.

Discussion
Parameterising sub-shelf melt rates in ice-sheet modelling is currently the only way to account for melting in largeensemble or multi-millennium simulations of the Antarctic ice sheet (DeConto and Pollard, 2016), 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 oceaninduced 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. 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  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 sup-plied 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 Warm 1 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., , 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., 2012Pattyn et al., , 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 farfield 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 Antarc-tic realistic ocean-ice-sheet systems in order to improve sea level projections.

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 icesheet 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 nonlocal 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 oceanice-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 https://github.com/ ElmerCSC/elmerfem (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 https://doi.org/10.5281/zenodo.2562838 (Jourdain and Favier, 2019); the NEMO set-up version Feb-2019, available at https://doi.org/10.5281/zenodo.2562731 (Jourdain, 2019a); and the Elmer/Ice set-up version 1.2, available at https://doi.org/10.5281/zenodo.2563156 (Jourdain, 2019b The Glen's flow law relates deviatoric stresses τ ij to strain ratesε ij as follows: 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: with τ b the basal friction, C s a friction parameter, u b the basal velocity, C max 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. The plume parameterisation is derived empirically from the results of a full plume model (Jenkins, 2014) 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 (Jenkins, 2014). 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: 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 Figure C1. Similar to Fig. 9 but to evaluate the use of the local gradient in PME 3 (which gives PME 5 ) and PME 4 (which gives PME 6 ) 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 PME 5 and PME 6 is done with α = 0.34 and α = 0.65, respectively (Sect. 3.2). The grey shading is only to ease the comparisons between the parameterisations.
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). An equivalent solution to the problem of non-uniform farfield 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 groundingline depth and the local slope in that direction were used to scale the melt rates. In PME 1 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 (PME 2 ) decreases the checkerboard noise in the PME 1 melt rates (Fig. 4) but does not greatly influence the results. In PME 3 and PME 4 , 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 PME 1 (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 z gl (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 re-called here below: z gl (i, j ) = 1 N ij valid n z n (i, j ), tan[θ (i, j )] = 1 N ij valid n s n (i, j ), with N ij the number of valid directions, z n (i, j ) the grounding line depth, and s n (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).

D2 PME 2
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 groundingline 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.

D3 PME 3
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).

D4 PME 4
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.

D5 PME 5
This implementation is similar to PME 3 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.

D6 PME 6
This implementation is similar to PME 4 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 D1. Computing the effective grounding-line depth and angle for the four plume parameterisation implementations. Panel (a) illustrates both the published implementation PME 1 (Lazeroms et al., 2018) and the one appearing in the corresponding discussion article PME 2 (Lazeroms et al., 2018), here with 12 directions (vs. 64 triangles used in the present paper). PME 1 and PME 2 differ by the calculation of the effective slope, which is not depicted here. Panel (b) illustrates the simple implementation PME 3 and panel (c) illustrates the asymmetric implementation PME 4 .