Antarctic sub-shelf melt rates via PICO

Ocean-induced melting below ice shelves is one of the dominant drivers for mass loss from the Antarctic Ice Sheet at present. An appropriate representation of sub-shelf melt rates is therefore essential for model simulations of marine-based ice sheet evolution. Continental-scale ice sheet models often rely on simple melt-parameterizations, in particular for long-term simulations, when fully coupled ice–ocean interaction becomes computationally too expensive. Such parameterizations can account for the influence of the local depth of the ice-shelf draft or its slope on melting. However, they do not capture the effect of ocean circulation underneath the ice shelf. Here we present the Potsdam Ice-shelf Cavity mOdel (PICO), which simulates the vertical overturning circulation in ice-shelf cavities and thus enables the computation of sub-shelf melt rates consistent with this circulation. PICO is based on an ocean box model that coarsely resolves ice shelf cavities and uses a boundary layer melt formulation. We implement it as a module of the Parallel Ice Sheet Model (PISM) and evaluate its performance under present-day conditions of the Southern Ocean. We identify a set of parameters that yield two-dimensional melt rate fields that qualitatively reproduce the typical pattern of comparably high melting near the grounding line and lower melting or refreezing towards the calving front. PICO captures the wide range of melt rates observed for Antarctic ice shelves, with an average of about 0.1 m a⁻¹ for cold sub-shelf cavities, for example, underneath Ross or Ronne ice shelves, to 16 m a⁻¹ for warm cavities such as in the Amundsen Sea region. This makes PICO a computationally feasible and more physical alternative to melt parameterizations purely based on ice draft geometry.

1 Introduction

Dynamic ice discharge across the grounding lines into floating ice shelves is the main mass loss process of the Antarctic Ice Sheet. Surrounding most of Antarctica's coastlines, the ice shelves themselves lose mass by ocean-induced melting from below or calving of icebergs (Depoorter et al., 2013; Liu et al., 2015). Observations show that many Antarctic ice shelves are thinning at present, driven by enhanced sub-shelf melting (Paolo et al., 2015; Pritchard et al., 2012). Thinning reduces the ice shelves' buttressing potential, i.e., the restraining force at the grounding line provided by the ice shelves (Dupont and Alley, 2005; Gudmundsson et al., 2012; Thomas, 1979), and can thereby accelerate upstream glacier flow. The observed acceleration of tributary glaciers is seen as the major contributor to the current mass loss in the West Antarctic Ice Sheet (Pritchard et al., 2012). In particular, the recent dynamic ice loss in the Amundsen Sea sector (MacGregor et al., 2012; Mouginot et al., 2014) is associated with high melt rates that result from inflow of relatively warm circumpolar deep water (CDW) in the ice-shelf cavities (D. M. Holland et al., 2008; Hellmer et al., 2017; Jacobs et al., 2011; Pritchard et al., 2012; Schmidtko et al., 2014; Thoma et al., 2008). Also in East Antarctica, particularly at Totten glacier, as well as along the Southern Antarctic Peninsula, glacier thinning seems to be linked to CDW reaching the deep grounding lines (Greenbaum et al., 2015; Wouters et al., 2015). An appropriate representation of melt rates at the ice–ocean interface is hence crucial for simulating the dynamics of the Antarctic Ice Sheet. Melting in ice-shelf cavities can occur in different modes that depend on the ocean properties in the proximity of the ice shelf, the topography of the ocean bed and the ice-shelf subsurface (Jacobs et al., 1992). Antarctica's ice-shelf cavities can be classified into “cold” and “warm” with typical mean melt rates ranging from 𝒪(0.1–1.0) m a⁻¹ in “cold” cavities as for the Filchner-Ronne Ice-Shelf and 𝒪(10) m a⁻¹ in “warm” cavities like the one adjacent to Pine Island Glacier (Joughin et al., 2012). For the “cold” cavities of the large Ross, Filchner-Ronne and Amery ice shelves, freezing to the shelf base is observed in the shallower areas near the center of the ice shelf and towards the calving front (Moholdt et al., 2014; Rignot et al., 2013).

Since the stability of the ice sheet is strongly linked to the dynamics of the buttressing ice shelves, it is essential to adequately represent their mass balance. A number of parameterizations with different levels of complexity have been developed to capture the effect of sub-shelf melting. Simplistic parameterizations that depend on the local ocean and ice-shelf properties have been applied in long-term and large-scale ice sheet simulations (Favier et al., 2014; Joughin et al., 2014; Martin et al., 2011; Pollard and DeConto, 2012). These parameterizations make melt rates piece-wise linear functions of the depth of the ice-shelf draft (Beckmann and Goosse, 2003) or of the slope of the ice-shelf base (Little et al., 2012). Other models describe the evolution of melt-water plumes forming at the ice-shelf base (Jenkins, 1991). Plumes evolve depending on the ice-shelf draft and slope, sub-glacial discharge and entrainment of ambient ocean water. This approach has been applied to models with characteristic conditions for Antarctic ice shelves (Holland et al., 2007; Lazeroms et al., 2018; Payne et al., 2007) and for Greenland outlet glaciers and fjord systems (Beckmann et al., 2018; Carroll et al., 2015; Jenkins, 2011). Interactively coupled ice–ocean models that resolve both the ice flow and the water circulation below ice shelves are now becoming available (De Rydt and Gudmundsson, 2016; Goldberg et al., 2012; Seroussi et al., 2017; Thoma et al., 2015). There is a community effort to better understand effects of ice–ocean interaction in such coupled models (Asay-Davis et al., 2016). However, as ocean models have many more degrees of freedom than ice sheet models and require for much shorter time steps, coupled simulations are currently limited to short time scales (on the order of decades to centuries).

Here, we present the Potsdam Ice-shelf Cavity mOdel (PICO), which provides sub-shelf melt rates in a computationally efficient manner and accounts for the basic vertical overturning circulation in ice-shelf cavities driven by the ice pump (Lewis and Perkin, 1986). It is based on the earlier work of Olbers and Hellmer (2010) and is implemented as a module in the Parallel Ice Sheet Model (PISM: Bueler and Brown, 2009; Winkelmann et al., 2011).¹ Ocean temperature and salinity at the depth of the bathymetry in the continental shelf region serve as input data. PICO allows for long-term simulations (on centennial to millennial time scales) and for large ensembles of simulations which makes it applicable, for example, in paleo-climate studies or as a coupling module between ice-sheet and Earth System models.

In this paper, we give a brief overview of the cavity circulation and melt physics and describe the ocean box geometry in PICO and implementation in PISM in Sect. 2. In Sect. 3, we derive a valid parameter range for present-day Antarctica and compare the resulting sub-shelf melt rates to observational data. This is followed by a discussion of the applicability and limitations of the model (Sect. 4) and conclusions (Sect. 5).

2 Model description

PICO is developed from the ocean box model of Olbers and Hellmer (2010), henceforth referred to as OH10. The OH10 model is designed to capture the basic overturning circulation in ice-shelf cavities which is driven by the “ice pump” mechanism: melting at the ice-shelf base near the grounding line reduces salinity and the ambient ocean water becomes buoyant, rising along the ice-shelf base towards the calving front. Since the ocean temperatures on the Antarctic continental shelf are generally close to the local freezing point, density variations are primarily controlled by salinity changes. Melting at the ice-shelf base hence reduces the density of ambient water masses, resulting in a haline driven circulation. Buoyant water rising along the shelf base draws in ocean water at depth, which flows across the continental shelf towards the deep grounding lines of the ice shelves. The warmer these water masses are, the stronger the melting-induced ice pump will be. The OH10 box model describes the relevant physical processes and captures this vertical overturning circulation by defining consecutive boxes following the flow within the ice-shelf cavity. The strength of the overturning flux q is determined from the density difference between the incoming water masses on the continental shelf and the buoyant water masses near the deep grounding lines of the ice shelf.

Table 1 PICO parameters and typical values.

The coefficients in the equation of state (EOS), the turbulent exchange velocities for heat and salt, are taken from Olbers and Hellmer (2010). We linearized the potential freezing temperature equation with a least-squares fit with salinity values over a range of 20–40 PSU and pressure values of 0–10⁷ Pa using Gibbs SeaWater Oceanographic Package of TEOS-10 (McDougall and Barker, 2011). All values are kept constant, except for ${\mathit{\gamma }}_{\mathrm{T}}^{*}$ and C, which vary between experiments. The values of these two parameters are the best-fit from Sect. 3.1.

As PICO is implemented in an ice sheet model with characteristic time scales much slower than typical response times of the ocean, we assume steady state ocean conditions and hence reduce the complexity of the governing equations of the OH10 model. Using this assumption facilitates adaptive box adjustment to grounding line migration, especially since PICO transfers the box model approach into two horizontal dimensions. We assume stable vertical stratification; OH10 found that a circulation state for an unstable vertical water column, which would imply a high (parametrized) diffusive transport between boxes, only occurs transiently (Olbers and Hellmer, 2010). This motivates neglecting the diffusive heat and salt transport between boxes which is small under these conditions. Without diffusive transport between the boxes, some of the original ocean boxes from OH10 become passive and can be incorporated into the governing equations of the set of boxes used in PICO. We explicitly model a single open ocean box which provides the boundary conditions for the boxes adjacent to the ice-shelf base following the overturning circulation, as shown in Fig. 1. In order to better resolve the complex melt patterns, PICO adapts the number of boxes based on the evolving geometry of the ice shelf. These simplifying assumptions allow us to analytically solve the system of governing equations, which is presented in the following two sections. A detailed derivation of the analytic solutions is given in Appendix A. In Sect. 2.3, we describe how the ice-model grid relates to the ocean box geometry of PICO. The system of equations is solved locally on the ice-model grid, as described in Sect. 2.4. Table 1 summarizes the model parameters and typical values.

Figure 1 Schematic view of the PICO model. The model mimics the overturning circulation in ice-shelf cavities: ocean water from box B₀ enters the ice-shelf cavity at the depth of the sea floor and is advected to the grounding line box B₁. Freshwater influx from melting at the ice-shelf base makes the water buoyant, causing it to rise. The cavity is divided into n boxes along the ice-shelf base. Generally, the highest melt rates can be found near the grounding line, with lower melt rates or refreezing towards the calving front.

2.1 Physics of the overturning circulation in ice-shelf cavities

PICO solves for the transport of heat and salt between the ocean boxes as depicted in Fig. 1. Although box B₀, which is located at depth between the ice-shelf front and the edge of the continental shelf, does not extend into the shelf cavity, its properties are transported unchanged from box B₀ to box B₁ near the grounding line. The heat and salt balances for all boxes in contact with the ice-shelf base (boxes B_k for $k\in \mathit{\{}\mathrm{1},\mathrm{\dots }n\mathit{\}}$) can be written as

$$\begin{array}{ll}{\displaystyle }{V}_k\dot{T_k}=& {\displaystyle }q{T}_{k-\mathrm{1}}-q{T}_k+A_k{m}_k{T}_{bk}-A_k{m}_k{T}_k\\ \text{(1)}& {\displaystyle }& {\displaystyle }+A_k{\mathit{\gamma }}_{\mathrm{T}}\left(T_{bk}-T_k\right),{\displaystyle }{V}_k\dot{S_k}=& {\displaystyle }q{S}_{k-\mathrm{1}}-q{S}_k+A_k{m}_k{S}_{bk}-A_k{m}_k{S}_k\\ \text{(2)}& {\displaystyle }& {\displaystyle }+A_k{\mathit{\gamma }}_{\mathrm{S}}\left(S_{bk}-S_k\right).\end{array}$$

The local application of these equations for each ice model cell is described in Sect. 2.4. Since we assume steady circulation, the terms on the left-hand side are neglected. For the box B_k with volume V_k, heat or salt content change due to advection from the adjacent box B_k−1 with overturning flux q (first term on the right-hand side) and due to advection to the neighboring box B_k+1 (or the open ocean for k=n; second term). Vertical melt flux into the box B_k across the ice–ocean interface with area A_k (third term) and out of the box (fourth term) play a minor role and are neglected in the analytic solution of the equation system employed in PICO (a detailed discussion of these terms is given in Jenkins et al., 2001). The melt rate m_k is negative if ambient water freezes to the shelf base. The last term represents heat and salt changes via turbulent, vertical diffusion across the boundary layer beneath the ice–ocean interface. The parameters γ_T and γ_S are the turbulent heat and salt exchange velocities which we assume, following OH10, to be constant.

The overturning flux q > 0 is assumed to be driven by the density difference between the ocean reservoir box B₀ and the grounding line box B₁. This is parametrized as in OH10 as

$$\begin{array}{}\text{(3)}& {\displaystyle }q=C\left({\mathit{\rho }}_{\mathrm{0}}-{\mathit{\rho }}_{\mathrm{1}}\right),\end{array}$$

where C is a constant overturning coefficient that captures effects of friction, rotation and bottom form stress.² The circulation strength in PICO is hence determined by density changes through sub-shelf melting in the grounding zone box B₁. From there, water follows the ice-shelf base towards the open ocean assuming the overturning flux q to be the same for all subsequent boxes. Ocean water densities are computed assuming a linear approximation of the equation of state:

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathit{\rho }={\mathit{\rho }}_{*}\left(\mathrm{1}-\mathit{\alpha }\left(T-T_{*}\right)+\mathit{\beta }\left(S-S_{*}\right)\right),\end{array}$$

where $T_{*}=\mathrm{0}$ ^∘C, $S_{*}=\mathrm{34}$ PSU and α, β and ρ_* are constants with values given in Table 1.

2.2 Melting physics

Melting physics are derived from the widely used 3-equation model (Hellmer and Olbers, 1989; Holland and Jenkins, 1999), which assumes the presence of a boundary layer below the ice–ocean interface. The temperature at this interface in box B_k is assumed to be at the in situ freezing point T_bk, which is linearly approximated by

$$\begin{array}{}\text{(5)}& {\displaystyle }{T}_{bk}=a{S}_{bk}+b-c{p}_k,\end{array}$$

where p_k is the overburden pressure, here calculated as static-fluid pressure given by the weight of the ice on top. At the ice–ocean interface, the heat flux from the ambient ocean across the boundary layer due to turbulent mixing, $Q_{\mathrm{T}}={\mathit{\rho }}_{\mathrm{w}}{c}_{\mathrm{p}}{\mathit{\gamma }}_{\mathrm{T}}\left(T_{bk}-T_k\right)$, equals the heat flux due to melting or freezing $Q_{\mathrm{Tb}}=-{\mathit{\rho }}_{\mathrm{i}}L{m}_k$. Neglecting heat flux into the ice, the heat balance equation thus reads

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{T}}\left(T_{bk}-T_k\right)=-\mathit{\nu }\mathit{\lambda }{m}_k,\end{array}$$

where $\mathit{\nu }={\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{w}}\sim \mathrm{0.89}$, $\mathit{\lambda }=L/c_{\mathrm{p}}\sim \mathrm{84}$ ^∘C. We obtain the salt flux boundary condition as the balance between turbulent salt transfer across the boundary layer, $Q_{\mathrm{S}}={\mathit{\rho }}_{\mathrm{w}}{\mathit{\gamma }}_{\mathrm{S}}\left(S_{bk}-S_k\right)$, and reduced salinity due to melt water input, $Q_{\mathrm{Sb}}=-{\mathit{\rho }}_{\mathrm{i}}{S}_{bk}{m}_k$,

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{S}}\left(S_{bk}-S_k\right)=-\mathit{\nu }{S}_{bk}{m}_k.\end{array}$$

To compute melt rates, we apply a simplified version of the 3-equations model (Holland and Jenkins, 1999; McPhee, 1992, 1999) which allows for a simple, analytic solution of the system of governing equations. It has been shown that this formulation yields realistic heat fluxes (McPhee, 1992, 1999). This simplification is used only for melt rates, we nevertheless solve for the boundary layer salinity which is central to the solution of the system of equations as detailed in Appendix A. Melt rates are given by

$$\begin{array}{}\text{(8)}& {\displaystyle }{m}_k=-{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{T}}^{*}}{\mathit{\nu }\mathit{\lambda }}}\left(a{S}_k+b-c{p}_k-T_k\right),\end{array}$$

with ambient ocean temperature T_k and salinity S_k in box B_k. Here, we use the effective turbulent heat exchange coefficient ${\mathit{\gamma }}_{\mathrm{T}}^{*}$. The relation between γ_T and ${\mathit{\gamma }}_{\mathrm{T}}^{*}$ is discussed in the Appendix A.

2.3 PICO ocean box geometry

PICO is implemented as a module in the three-dimensional ice sheet model PISM as described in Sect. 2.4. Since the original system of box-model equations is formulated for only one horizontal and one vertical dimension, it needed to be extended for the use in the three-dimensional ice sheet model. The system of governing equations as described in the previous two sections is solved for each ice shelf independently. PICO adapts the ocean boxes to the evolving ice shelves at every time step.

For any shelf D, we determine the number of ocean boxes n_D by interpolating between 1 and n_max depending on its size and geometry such that larger ice shelves are resolved with more boxes. The maximum number of boxes n_max is a model parameter; a value of 5 is suitable for the Antarctic setup, as discussed further in Sect. 3.2. We determine the number of boxes n_D for each individual ice shelf D with

$$\begin{array}{}\text{(9)}& {\displaystyle }{n}_D=\mathrm{1}+\text{rd}\left(\sqrt{d_{\mathrm{GL}}\left(D\right)/d_{max}}\left(n_{max}-\mathrm{1}\right)\right),\end{array}$$

where rd( ) rounds to the nearest integer. Here, d_GL(x,y) is the local distance to the grounding line from an ice-model grid cell with horizontal coordinates (x, y), d_GL(D) is the maximum distance within ice shelf D and d_max is the maximum distance to the grounding line within the entire computational domain.

Knowing the maximum number of boxes n_D for an ice shelf D, we next define the ocean boxes underneath it. The extent of boxes $B_{\mathrm{1}},\mathrm{\dots },B_{n_D}$ is determined using the distance to the grounding line and the shelf front. The non-dimensional relative distance to the grounding line r is defined as

$$\begin{array}{}\text{(10)}& {\displaystyle }r\left(x,y\right)=d_{\mathrm{GL}}\left(x,y\right)/\left(d_{\mathrm{GL}}\left(x,y\right)+d_{\mathrm{IF}}\left(x,y\right)\right),\end{array}$$

with d_IF(x,y) the horizontal distance to the ice front. We assign all ice cells with horizontal coordinates $\left(x,y\right)\in D$ to box B_k if the following condition is met:

$$\begin{array}{}\text{(11)}& {\displaystyle }\mathrm{1}-\sqrt{\left(n_D-k+\mathrm{1}\right)/n_D}\le r\left(x,y\right)\le \mathrm{1}-\sqrt{\left(n_D-k\right)/n_D}.\end{array}$$

This leads to comparable areas for the different boxes within a shelf, which is motivated in Appendix B. Thus, for example, the box B₁ adjacent to the grounding line interacts with all ice shelf grid cells with $\mathrm{0}\le r\le \mathrm{1}-\sqrt{\left(n_D-\mathrm{1}\right)/n_D}$. Figure 3 shows an example of the ocean box areas for Antarctica. PICO does currently not account for melting along vertical ice cliffs, as, for example, the termini of some Greenland outlet fjords.

2.4 Implementation in the Parallel Ice Sheet Model

PICO is implemented in the open-source Parallel Ice Sheet Model (Bueler and Brown, 2009; Winkelmann et al., 2011). In the three-dimensional, thermo-mechanically coupled, finite-difference model, ice velocities are computed through a superposition of the shallow approximations for the slow, shear-dominated flow in ice sheets (Hutter, 1983) and the fast, membrane-like flow in ice streams and ice shelves (Morland, 1987). In PISM, the grounding lines (diagnosed via the flotation criterion) and ice fronts evolve freely. Grounding line movement has been evaluated in the model intercomparison project MISMIP3d (Feldmann et al., 2014; Pattyn et al., 2013).

PICO is synchronously coupled to the ice-sheet model, i.e., they use the same adaptive time steps. The cavity model provides sub-shelf melt rates and temperatures at the ice–ocean boundary to PISM, with temperatures being at the in situ freezing point. PISM supplies the evolving ice-shelf geometry to PICO, which in turn adjusts in each time step the ocean box geometry to the ice-shelf geometry as described in Sect. 2.3.

PICO computes the melt rates progressively over the ocean boxes, independently for each ice shelf. Since the ice-sheet model has a much higher resolution, each ocean box interacts with a number of ice-shelf grid cells. PICO applies the analytic solutions of the system of governing equations summarized in Sect. 2.1 and 2.2 locally to the ice model grid as detailed below. Model parameters that are varied between the experiments are the effective turbulent heat exchange velocity ${\mathit{\gamma }}_{\mathrm{T}}^{*}$ from the melt parametrization described in Sect. 2.2 and the overturning coefficient C described in Sect. 2.1. The two parameter values are applied to the entire ice sheet.

Input for PICO in the ocean reservoir box B₀ is data from observations or large-scale ocean models in front of the ice shelves. Temperature T₀ and salinity S₀ are averaged at the depth of the bathymetry in the continental shelf region. In box B₁ adjacent to the grounding line, PICO solves the system of governing equations in each ice grid cell (x, y) to attain the overturning flux q(x,y), temperature T₁(x,y), salinity S₁(x,y), and the melting m₁(x,y) at its ice–ocean interface (given by the local solution of Eqs. 3, A12, A8 and 8, respectively). The model proceeds progressively from box B_k to box B_k+1 to solve for sub-shelf melt rate $m_{k+\mathrm{1}}\left(x,y\right)$, ambient ocean temperature $T_{k+\mathrm{1}}\left(x,y\right)$, and salinity $S_{k+\mathrm{1}}\left(x,y\right)$ (given by the local solution of Eqs. 13, A13 and A8, respectively) based on the previous solutions S_k and T_k in box B_k and conditions at the ice–ocean interface. PICO provides the boundary conditions T_k and S_k to box B_k+1 as the average over the ice-grid cells within box B_k, i.e.,

$$\begin{array}{}\text{(12)}& {\displaystyle }{T}_k=\langle T_k\left(x,y\right)\text{with}\left(x,y\right)\text{in}{B}_k\rangle \end{array}$$

and analogously for S_k, where 〈 〉 denotes the average.

The overturning is solved in Box B₁ and given by $q=\langle q\left(x,y\right)\text{with}\left(x,y\right)\text{in}{B}_{\mathrm{1}}\rangle$. Melt rates in box B_k are computed using the local overburden pressure p_k(x,y) in each ice-shelf grid cell that is given by the weight of the ice column provided by PISM, i.e.,

$$\begin{array}{}\text{(13)}& {\displaystyle }{m}_k\left(x,y\right)=-{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{T}}^{*}}{\mathit{\nu }\mathit{\lambda }}}\left(a{S}_k\left(x,y\right)+b-c{p}_k\left(x,y\right)-T_k\left(x,y\right)\right).\end{array}$$

This reflects the pressure dependence of heat available for melting and leads to a depth-dependent melt rate pattern within each box. The implications for energy and mass conservation are discussed in Sects. 3.2 and 4.

Figure 2 PICO input for Antarctic basins. The ice sheet, ice shelves and the surrounding Southern Ocean are split into 19 basins that are based on Zwally et al. (2012) and indicated by black contour lines and labels. For each ice shelf, the governing equations are solved separately with the respective oceanic boundary conditions. For ice shelves that cross basin boundaries, the input is averaged, weighted with the fractional area of the shelf within the corresponding basin. Numbers show the temperature and salinity input in each basin, obtained by averaging observed properties of the ocean water in front of the ice-shelf cavities at depth of the continental shelf (Schmidtko et al., 2014), indicated by the color shading. Grey lines show Antarctic grounding lines and ice-shelf fronts (Fretwell et al., 2013).

3 Results for present-day Antarctica

We apply PICO to compute sub-shelf melt rates for all Antarctic ice shelves under present-day conditions. Oceanic input for each ice shelf is given by observations of temperature (converted to potential temperature) and salinity (converted to practical salinity) of the water masses occupying the sea floor on the continental shelf (Schmidtko et al., 2014), averaged over the time period 1975 to 2012. Water masses within an ice-shelf cavity originate from source regions: neglecting ocean dynamics, we approximate these by averaging ocean properties on the depth of the continental shelf within regions that are chosen to encompass areas of similar, large-scale ocean conditions. Oceanic input is given for 19 basins of the Antarctic Ice Sheet, which are based on Zwally et al. (2012) and extended to the attached ice shelves and the surrounding Southern Ocean (Fig. 2). For each ice shelf, temperature T₀ and salinity S₀ in box B₀ are obtained by averaging the basin input weighted with the fractional area of the shelf within the corresponding basin. Figure 2 shows the basin-mean ocean temperature (shadings and numbers) and salinity (numbers) used.³

Here, we use $n_{max}=\mathrm{5}$ from which PICO determines the number of ocean boxes in each shelf via Eq. (9). Figure 3 displays the resulting extent of the ocean boxes for Antarctica, ordered in elongated bands beneath the ice shelves. For the large ice-shelf cavities of Filchner-Ronne and Ross, the ice–ocean boundary is divided into five ocean boxes while smaller ice shelves have two to four boxes (see Table 2). Introducing more than five ocean boxes has a negligible effect on the melt rates, as discussed in Sect. 3.2.

To validate our model, we run diagnostic simulations with PISM+PICO based on bed topography and ice thickness from BEDMAP2 (Fretwell et al., 2013), mapped to a grid with 5 km horizontal resolution. Diagnostic simulations allow us to assess the sensitivity of the model to the parameters C and ${\mathit{\gamma }}_{\mathrm{T}}^{*}$ and to the number of boxes n_max as well as the ice model resolution. Transient behavior is exemplified in a simulation starting from an equilibrium state of the Antarctic Ice Sheet forced with ocean temperature changes, see Video S1 in the Supplement and Sect. 3.3.

Figure 3 Extent of PICO ocean boxes for Antarctic ice shelves. Most ice shelves are split into two, three, or four ocean boxes interacting with the ice cells on a much higher resolution. The largest ice shelves, Filchner-Ronne and Ross, have five ocean boxes. One ocean box typically corresponds to many ice-shelf grid cells.

Table 2 Results from the reference simulation as displayed in Fig. 5.

The number of boxes for each is ice-shelf is given by b_n, T₀ (S₀) is the temperature (salinity) in ocean box B₀, T_n (S_n) the temperature (salinity) averaged over the ocean box at the ice-shelf front, $\mathrm{\Delta }T=T_n-T_{\mathrm{0}}$ and $\mathrm{\Delta }S=S_n-S_{\mathrm{0}}$. m is the average sub-shelf melt rate, m_observed the observed melt rates from Rignot et al. (2013). q is the overturning flux. Unit of temperatures is ^∘C, salinity is given in PSU, melt rates in m a⁻¹ and overturning flux in Sv.

3.1 Sensitivity to model parameters C and ${\mathit{\gamma }}_{\mathrm{T}}^{*}$

We test the sensitivity of sub-shelf melt rates to the model parameters for overturning strength $C\in [\mathrm{0.1},\mathrm{9}]$ Sv m³ kg⁻¹ and the effective turbulent heat exchange velocity ${\mathit{\gamma }}_{\mathrm{T}}^{*}\in [\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}},\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}]$ m s⁻¹. These ranges encompass the values identified in OH10, discussed further in Appendix A. The same parameters for C and ${\mathit{\gamma }}_{\mathrm{T}}^{*}$ are applied to all shelves. We constrain the results by the following qualitative criteria (1) and (2), as well as the quantitative constraints (3) and (4), summarized in Fig. 4:

Figure 4 Sensitivity of PICO sub-shelf melt rates to the overturning coefficient C and the effective turbulent heat exchange coefficient ${\mathit{\gamma }}_{\mathrm{T}}^{*}$. Black contour indicates the valid range of parameters, all other parameter combinations are excluded by one of the following criteria: no freezing may occur in the first ocean box (a), mean basal melt rates must decrease between the first and second ocean box (b). Green contour indicates the valid range of parameters where the following quantitative constraints are additionally met: mean basal melt rates for Filchner-Ronne Ice Shelf should be within the range of 0.01 to 1.0 m a⁻¹ (c), mean basal melt rates for Pine Island Glacier Ice Shelf should be within the range of 10.0 to 25.0 m a⁻¹ (d). The best-fit parameters (brown contour) minimize the root-mean-square error of mean melt rates to observations for both ice shelves.

Figure 5 Sub-shelf melt rates for present-day Antarctica computed with PISM+PICO. The mean melt rate per ice shelf (upper numbers) is compared to the observed melt rates (lower numbers with uncertainty ranges) from Rignot et al. (2013). In the model, the same parameters ${\mathit{\gamma }}_{\mathrm{T}}^{*}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ and C=1 Sv m³ kg⁻¹ are applied to all ice shelves around Antarctica. The respective oceanic boundary condition are shown in Fig. 2. Ice geometry and bedrock topography are from the BEDMAP2 data set on 5 km resolution (Fretwell et al., 2013). Refreezing occurs in some parts of the larger shelves like Filchner-Ronne and Ross.

Figure 6 Sensitivity of PICO sub-shelf melt rates to ocean temperature changes for the entirety of Antarctica (black), Filchner-Ronne Ice Shelf (blue), and Pine Island Glacier Ice Shelf (red). Ocean input temperatures are varied by 0.1 up to 2 ^∘C. Melting depends quadratically on temperature for “cold” cavities like the one adjacent to Filchner-Ronne, and linearly for “warm” cavities like the ones in the Amundsen Region.

Criterion (1). Freezing must not occur in the first box B₁ of any ice shelf, i.e., the ocean box closest to the grounding line. Freezing in box B₁ would increase ambient salinity, and since the overturning circulation in ice-shelf cavities is mainly haline driven, the circulation would shut down, violating the model assumption q>0 (see Sect. 2). As shown in Fig. 4a, the condition is not met for a combination of relatively high turbulent heat exchange and relatively low overturning parameters. In such cases, freezing near the grounding line occurs because of the strong heat exchange between the ambient ocean and the ice–ocean boundary in box B₁ that cannot be balanced by the resupply of heat from the open ocean through overturning.

Criterion (2). Sub-shelf melt rates must decrease between the first and second box for each ice shelf. This condition is based on general observations of melt-rate patterns and on the assumption that ocean water masses move consecutively through the boxes and cool down along the way, as long as melting in these boxes outweighs freezing. As shown in Fig. 4b, this condition is violated for either high overturning and low turbulent heat exchange or, vice versa, low overturning and high turbulent heat exchange. An appropriate balance between the strength of these values is hence necessary for a realistic melt rate pattern.

If criterion (1) or (2) fails, basic assumptions of PICO are violated. Thus, we choose the model parameters ${\mathit{\gamma }}_{\mathrm{T}}^{*}$ and C such that both criteria are strictly met. The following quantitative, observational constraints (3) and (4) compare modeled average melt rates with observations and thus depend on our choice of valid ranges. We choose Filchner-Ronne Ice Shelf and the ice shelf adjacent to Pine Island Glacier to further constrain our model parameters for Antarctica. These shelves represent two different types regarding both the mode of melting and the ice-shelf size.

Observational constraint (3). Average melt rates in Filchner-Ronne Ice Shelf comply with the classification of a “cold” cavity and lie between 0.01 and 1.0 m a⁻¹ (Fig. 4c).

Observational constraint (4). For Pine Island Glacier, with “warm” ocean conditions, average melt rates lie between 10 and 25 m a⁻¹ (Fig. 4d).

Generally, an increase in overturning strength C will supply more heat and thus yield higher melt rates, especially for the large and “cold” ice shelves like Filchner-Ronne. Larger C leads to higher melt rates also in the smaller and “warm” Pine-Island Glacier Ice Shelf but the effect is less pronounced. In contrast, the turbulent heat exchange alters melting, particularly in small ice shelves, while it might decrease melt rates in large ice shelves with “cold” cavities. Hence, modeled melting in the Filchner-Ronne Ice Shelf is dominated by overturning while in the Amundsen region melting is dominated by turbulent exchange across the ice–ocean boundary layer. For three different parameter combinations, the resulting spatial patterns of melt rates in the Filchner-Ronne and Pine Island regions are displayed in Fig. S1 in the Supplement.

All of the above criteria restrict the parameter space to a bounded set with lower and upper limits as depicted by the green contour lines in Fig. 4. We determine a best-fit pair of parameters which minimizes the root-mean-square error of average melt rates for both ice shelves. The valid range of model parameters around the best-fit parameters with C=1 Sv m³ kg⁻¹ and ${\mathit{\gamma }}_{\mathrm{T}}^{*}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$ m a⁻¹ compares well with parameters found in OH10 and Holland and Jenkins (1999), discussed further in Appendix A.

3.2 Diagnostic melt rates for present-day Antarctica

Using the best-fit values C=1 Sv m³ kg⁻¹ and ${\mathit{\gamma }}_{\mathrm{T}}^{*}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ found in Sect. 3.1, we apply PICO to present-day Antarctica, solving for sub-shelf melt rates and water properties in the ocean boxes. This model simulation is referred to as “reference simulation” hereafter.

The average melt rates computed with PICO range from 0.06 m a⁻¹ under the Ross Ice Shelf to 16.15 m a⁻¹ for the Amundsen Region (Fig. 5). Consistent with the model assumptions, melt rates are highest in the vicinity of the grounding line and decrease towards the calving front. In some regions of the large ice shelves, refreezing occurs, e.g., towards the center of Filchner-Ronne or Amery ice shelves. The melt pattern depends on the local pressure melting temperature. Thus, melt rates are highest where the shelf is thickest, i.e., near the grounding lines within box B₁. Furthermore, freezing and melting can occur in the same box determined by local ice-shelf thickness. For the vast majority of ice shelves, the modeled average melt rates compare well with the observed ranges derived from Table 1 in Rignot et al. (2013). Exceptions are the Wilkins and Stange ice shelves (in basin 16) with average modeled melt rates of 9.50 m a⁻¹, which is much higher than the observed range of 1.46 ± 1.0 m a⁻¹. This is most likely due to the ocean temperature input for this basin (1.17 ^∘C, see Fig. 2) which is higher than for the basin containing Pine Island located nearby (0.47 ^∘C, basin 14), explaining why the melt rates are of the same order of magnitude in these basins. Modification of water masses flowing into the shelf cavities, not captured by PICO, might explain the low observed melt rates in this area despite the relatively high ocean temperatures.

For all ice shelves, ocean temperatures and salinities consistently decrease in overturning direction, i.e., from the ocean reservoir box B₀ to the last box adjacent to the ice front B_n, as shown in Table 2. Most shelves contain small areas with accretion with a maximum rate of −1.22 m a⁻¹ for the Amery Ice Shelf, see Table S1 in the Supplement. No freezing occurs at the Western Antarctic Peninsula nor in the Amundsen and Bellingshausen Seas. A detailed map of sub-shelf melt rates in this region as well as for Filchner-Ronne Ice Shelf can be found in the middle panels of Fig. S1. For the Filchner-Ronne Ice Shelf melt rates vary between −0.67 and 1.76 m a⁻¹ and for Pine Island Glacier, melt rates range from 12.39 to 21.01 m a⁻¹.

PICO tends to smooth out melt rate patterns, see Fig. S4: for Filchner-Ronne and Ross ice shelves the deviations in observed melt rates (Moholdt et al., 2016) from the average melting are larger than the deviations modeled in PICO. The box-wide averages compare well and are at the same order of magnitude for both ice shelves, except for modeled accretion in the later boxes towards the shelf front which is not reflected in the observations. This disagreement might be explained by the fact that the overturning circulation in the model cannot reach neutral density and detach from the ice-shelf base while flowing towards the shelf front.

Aggregated over all Antarctic ice shelves, the total melt flux is 1718 Gt a⁻¹, close to the observed estimate of 1500 ± 237 Gt a⁻¹ (Rignot et al., 2013). Overturning fluxes in our reference simulation range from 0.02 Sv for the small Drygalski Ice Shelf to 0.32 Sv for Wilkings, Stange, Bach and George VI ice shelves. Because these are treated in PICO as one ice shelf and they have a high input temperature, these ice shelves reach together this high overturning value. The second highest value of 0.23 Sv is found for Getz Ice Shelf in the Bellingshausen Sea. These overturning fluxes compare well with the estimates in OH10.

PICO solves the system of governing equations locally in each ice-model grid cell and calculates input for each ocean box as an average over the previous box as described in Sect. 2.4. Due to these model assumptions and because temperature, salinity, overturning and melting are non-linear functions of pressure in the first box (Eq. A12), mass and energy are a priori not perfectly conserved. In Table S1, we compare (within individual shelves) heat fluxes into the ice-shelf cavities with the heat flux out of the cavities into the ocean and the latent heat flux for melting. For the whole of Antarctica, the deviation in heat flux is 403.63 G J s⁻¹ which is equivalent to 2.2 % of the latent heat flux for melting. The per basin deviations are generally low (< 5 %), except around the large Filchner-Ronne and Ross ice shelves (68 and 43 %). In PICO we assume q to be constant, neglecting changes due to melt water input along the shelf base. This melt water input amounts to 3.06 % or less of the overturning flux within ice shelves, and 0.62 % for the entire continent, discussed in Sect. 2.1.

Melt rates are strongly affected by changes in the ambient ocean temperatures, see Fig. 6. The dependence is approximately linear for high and quadratic for lower ambient ocean temperatures. This relationship is similar to the one observed in OH10 and as expected from the governing equations. In Pine Island Glacier, melt rates increase by approximately 6 m a⁻¹ per degree of warming. Changes in the ice-sheet model resolution have little effect on the resulting melt rates (Fig. S2). For increasing the maximum number of boxes n_max, average melt rates converge to almost constant values for $n_{max}\ge \mathrm{5}$ within all basins, compare Fig. S3.

3.3 Transient evolution of PICO boxes and melt rates

PICO is capable of adjusting to changing ice-shelf geometries and migrating grounding lines. We demonstrate its behavior as a module in PISM in a transient simulation. Based on an Antarctic equilibrium state at 8 km resolution comparable to the state submitted to initMIP (Nowicki et al., 2016), we run PISM+PICO over time with a simple temperature forcing applied: starting from equilibrium conditions, ocean temperatures increase linearly over 50 years until an ocean-wide warming of 1 ^∘C is reached. It is then held constant over the next 250 modeled years. The Video S1 shows the temporal evolution of the ocean temperature input for the ice shelf adjacent to Pine Island Glacier as well as the Filchner-Ronne Ice Shelf. The ocean temperature increase enhances the sub-shelf melting for both ice shelves, especially in the first box. Ice-shelf thinning reduces buttressing and causes the grounding lines to retreat with the ocean boxes adjusting accordingly.

4 Discussion

PICO models the dominant vertical overturning circulation in ice-shelf cavities and translates ocean conditions in front of the ice shelves, either from observations or large-scale ocean models, into physically based sub-shelf melt rates. For present-day ocean fields and ice-shelf cavity geometries, PICO as an ocean module in PISM reproduces average melt rates of the same order of magnitude as observations for most Antarctic ice shelves. With a single combination of overturning parameter C and effective turbulent heat exchange parameter ${\mathit{\gamma }}_{\mathrm{T}}^{*}$ applied to all shelves, a wide range of melt rates for the different ice shelves is obtained, reproducing the large-scale patterns observed in Antarctica. The results are consistent across different ice-sheet and cavity model resolutions. Additionally, PICO reproduces the common pattern of maximum melt close to the grounding line and decreasing melt rates towards the ice-shelf front, eventually with re-freezing in the shallow parts of the large ice shelves. The governing model equations are solved for individual grid cells of the ice sheet model (and not for each ocean box with representative depth value), which allows spatial variability in the resulting melt rate field at comparatively smaller scale. PICO can adapt to evolving cavities and is applicable to ice shelves in two horizontal dimensions. It is hence suited as a sub-shelf melt module for ice-sheet models.

Yet PICO is a coarse model designed as an ocean coupler for large-scale ice-sheet models. It is based on the OH10 model and hence shares some simplifying assumptions with that model: PICO does not resolve ocean dynamics and it parameterizes the vertical overturning circulation in the ice-shelf cavities which is given by one value for each ice shelf. The underlying box model equations of PICO do not resolve horizontal ocean circulation in the ice-shelf cavities driven by the Coriolis force nor seasonal melt variation due to intrusion of warm water from the calving front during Austral summer. Hence, we do not expect that horizontal variations or small scale patterns of basal melt are accurately captured in PICO. Due to the box model formulation, maximum melting in PICO is found directly at the grounding line and not slightly downstream as seen in the high-resolution coupled ice–ocean simulation by, e.g., De Rydt and Gudmundsson (2016). We find that PICO tends to produce smoother melt rate patterns than observed, though the box-wide averages are in good agreement with observations. The effect on ice dynamics of small-scale differences in melt patterns in relation to the large-scale mean melt fluxes is not well established yet and needs further investigation. Following OH10, meltwater does not contribute to the volume flux in the cavity in PICO, introducing a minor error regarding mass conservation. The relative error regarding mass conservation is however small and below 0.7 % of the total overturning strength in our reference run.

A necessary condition for the box model to work is further the assumption that melting outweighs accretion in box B₁ which is consistent with the majority of available observations. In PICO, melt rates show a quadratic dependency on ocean temperature input for lower temperatures, e.g., in the Filchner-Ronne basin, and a rather linear dependency for higher temperatures, e.g., in the Amundsen basin. This is consistent with the results from OH10 and the implemented melting physics assuming a constant coefficient for turbulent heat exchange. In contrast, P. R. Holland et al. (2008) employ a dependency of the turbulent heat exchange coefficient on the velocity of the overturning circulation, suggesting melt rates to respond quadratically to warming of the ambient ocean water. Here we follow the approach taken in OH10.

Differently from the OH10 model and relying on much longer timescales of ice dynamics in comparison to ocean dynamics, PICO assumes the overturning circulation to be in steady state. In their analysis, OH10 find unstable vertical water columns to occur only transiently, and hence for PICO a stable stratification of the vertical water column is assumed. Under these conditions, diffusive transport between the boxes is generally small in OH10 and it is hence omitted in PICO. Because PICO also does not consider vertical variations in ambient ocean density, under-ice flow is prevented from reaching neutral density and detaching from the ice-shelf base on its way towards the shelf front. The spatial pattern of melting closer to the calving front of cold ice shelves may in such cases be not represented well.

PICO input is determined by averaging bottom temperatures and salinities over the continental shelf, this is done for 19 different basins. Hence PICO – as a coarse model – misses the nuances of how ocean currents transport and modify CDW over the regions being averaged. The procedure to determine melt rates in PICO is based on the assumption that ocean water that is present on the continental shelf can access the ice-shelf cavities and reach their grounding lines. This implies for example, that barriers like sills that may prevent intrusion of warm CDW are not accounted for and might explain why PICO melting is too high for the ice shelves located along the Southern Antarctic Peninsula. Such phenomena could be tested by varying the ocean input of PICO by evolving the temperature and salinity outside of the cavity over time. Because of the dependence of sub-shelf melting on the local pressure of the ice column above, the model is not fully energy conserving. For the estimated heat fluxes, the relative error is lower than 2.2 % of the latent heat flux due to sub-shelf melting. Hence, we consider our choice of model simplifications as justified, as it introduces small errors in the heat and mass balances in our reference simulation.

PICO is computationally fast, as it uses analytic solutions of the equations of motion with a small number of boxes along the ice shelf. As boundary conditions for PICO are aggregated based on predefined regional basins, the model can act as an efficient coupler of large-scale ice-sheet and ocean models. For this purpose, heat flux into the ice should be added to the boundary layer melt formulation.

5 Conclusions

The Antarctic Ice Sheet plays a vital role in modulating global sea level. The ice grounded below sea level in its marine basins is susceptible to ocean forcing and might respond nonlinearly to changes in ocean boundary conditions (Mercer, 1978; Schoof, 2007). We therefore need carefully estimated conditions at the ice–ocean boundary to better constrain the dynamics of the Antarctic ice and its contribution to sea-level rise for the past and the future.

The PICO model presented here provides a physics-based yet efficient approach for estimating the ocean circulation below ice shelves and the heat provided for ice-shelf melt. The model extends the one-horizontal dimensional ocean box model by OH10 to realistic ice-shelf geometries following the shape of the grounding line and calving front. PICO is a comparably simple alternative to full ocean models, but goes beyond local melt parameterizations, which do not account for the circulation below ice shelves. We find a set of possible parameters for present-day ocean conditions and ice geometries which yield PICO melt rates in agreement with average melt rate observations. PICO qualitatively reproduces the general pattern of ice-shelf melt, with high melting at the grounding line and low melting or refreezing towards the calving front. Its sensitivity to changes in input ocean temperatures and model parameters is comparable to earlier estimates (Holland and Jenkins, 1999; Olbers and Hellmer, 2010). The model accurately captures the large variety of observed Antarctic melt rates using only two calibrated parameters applied to all ice shelves.

The ocean models that are part of the large Earth system and global circulation models do not yet resolve the circulation below ice shelves. PICO is able to fill this gap and can be used as an intermediary between global circulation models and ice sheet models. We expect that PICO will be useful for providing ocean forcing to ice sheet models with the standardized input from climate model intercomparison projects like CMIP5 and CMIP6 (Eyring et al., 2016; Meehl et al., 2014; Taylor et al., 2012). Since PICO can deal with evolving ice-shelf geometries in a computationally efficient way, it is in particular suitable for modeling the ice sheet evolution on paleo-climate timescales as well as for future projections.

PICO is implemented as a module in the open-source Parallel Ice Sheet Model (PISM). The source code is fully accessible and documented as we want to encourage improvements and implementation in other ice sheet models. This includes the adaption to other ice sheets than present-day Antarctica.

Appendix A: Derivation of the analytic solutions

Here, we derive the analytic solutions of the equations system describing the overturning circulation (see Sect. 2.1) and the melting at the ice–ocean interface (see Sect. 2.2).

For box B_k with k>1 we solve progressively for melt rate m_k, temperature T_k and salinity S_k in box B_k, dependent on the local pressure p_k, the area of box adjacent to the ice-shelf base A_k and the temperature T_k−1 and salinity S_k−1 of the upstream box B_k−1. For box B₁, we additionally solve for the overturning q as explained below. These derivations advance the ideas presented in the appendix of OH10. Assuming steady state conditions, the balance Eqs. (1) and (2) for box B_k from Sect. 2.1 are

$$\begin{array}{ll}{\displaystyle }\mathrm{0}& {\displaystyle }=q\left(T_{k-\mathrm{1}}-T_k\right)+A_k{\mathit{\gamma }}_{\mathrm{T}}\left(T_{bk}-T_k\right)+A_k{m}_k\left(T_{bk}-T_k\right),\\ \text{(A1)}& {\displaystyle }\mathrm{0}& {\displaystyle }=q\left(S_{k-\mathrm{1}}-S_k\right)+A_k{\mathit{\gamma }}_{\mathrm{S}}\left(S_{bk}-S_k\right)+A_k{m}_k\left(S_{bk}-S_k\right).\end{array}$$

The heat fluxes balance at the boundary layer interface, i.e., the heat flux across the boundary layer due to turbulent mixing $Q_{\mathrm{T}}={\mathit{\rho }}_{\mathrm{w}}{c}_{\mathrm{p}}{\mathit{\gamma }}_{\mathrm{T}}\left(T_{bk}-T_k\right)$ equals the heat flux due to melting or freezing $Q_{\mathrm{Tb}}=-{\mathit{\rho }}_{\mathrm{i}}L{m}_k$, omitting the heat flux into the ice. This yields

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{T}}\left(T_{bk}-T_k\right)=-\mathit{\nu }\mathit{\lambda }{m}_k,\end{array}$$

where $\mathit{\nu }={\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{w}}\sim \mathrm{0.89}$, $\mathit{\lambda }=L/c_{\mathrm{p}}\sim \mathrm{84}$ ^∘C. Regarding the salt flux balance in the boundary layer, with $Q_{\mathrm{S}}={\mathit{\rho }}_{\mathrm{w}}{\mathit{\gamma }}_{\mathrm{S}}\left(S_{bk}-S_k\right)$ at the lower interface of the boundary layer and “virtual” salt flux due to meltwater input $Q_{\mathrm{Sb}}=-{\mathit{\rho }}_{\mathrm{i}}{S}_{bk}{m}_k$, we obtain

$$\begin{array}{}\text{(A3)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{S}}\left(S_{bk}-S_k\right)=-\mathit{\nu }{S}_{bk}{m}_k.\end{array}$$

Inserting Eqs. (A2) and (A3) into Eqs. (A1) yields

$$\begin{array}{ll}{\displaystyle }\mathrm{0}& {\displaystyle }=q\left(T_{k-\mathrm{1}}-T_k\right)-A_k{m}_k\mathit{\nu }\mathit{\lambda }+A_k{m}_k\left(T_{bk}-T_k\right),\\ {\displaystyle }\mathrm{0}& {\displaystyle }=q\left(S_{k-\mathrm{1}}-S_k\right)-A_k{m}_k\mathit{\nu }{S}_{bk}+A_k{m}_k\left(S_{bk}-S_k\right).\end{array}$$

Comparing $\left(T_{bk}-T_k\right)\ll \mathit{\nu }\mathit{\lambda }\approx \mathrm{75}$ ^∘C, allows us to neglect the last term in the temperature equation. Considering the last two terms of the salinity equation, we find that $S_k>\left(\mathrm{1}-\mathit{\nu }\right)S_{bk}\approx \mathrm{0.1}{S}_{bk}$, allowing us to neglect the terms containing S_bk, which simplifies the equations to

$$\begin{array}{ll}{\displaystyle }\mathrm{0}& {\displaystyle }=q\left(T_{k-\mathrm{1}}-T_k\right)-A_k\mathit{\nu }\mathit{\lambda }{m}_k\\ \text{(A4)}& {\displaystyle }\mathrm{0}& {\displaystyle }=q\left(S_{k-\mathrm{1}}-S_k\right)-A_k{m}_k{S}_k.\end{array}$$

We use a simplified version of the melt law described by McPhee (1992) and detailed in Sect. 2.2, which makes use of Eqs. (6) and (5) in which the salinity in the boundary layer S_bk is replaced by salinity of the ambient ocean water:

$$\begin{array}{}\text{(A5)}& {\displaystyle }{m}_k=-{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{T}}^{*}}{\mathit{\nu }\mathit{\lambda }}}\left(a{S}_k+b-c{p}_k-T_k\right).\end{array}$$

Holland and Jenkins (1999) suggest that this simplification requires ${\mathit{\gamma }}_{\mathrm{T}}^{*}$ to be a factor of 1.35 to 1.6 smaller than γ_T in the 3-equation formulation for the constant values of γ_T ranging from $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ to $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ used in OH10. This implies that ${\mathit{\gamma }}_{\mathrm{T}}^{*}$ ranges from $\mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{5}}$ to $\mathrm{3.2}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹, which is consistent with the parameter range we derive in Sect. 3.1. We apply this assumption in the computation of melt rates. For the solution of the transport Eqs. (A1), it is essential to take the salinity of the boundary layer S_bk into account, since otherwise the salinity transport equation would reduce to $S_k=S_{k-\mathrm{1}}$ and the overturning circulation, which is predominantly haline driven, would be reduced. Inserting the simplified melt law in Eq. (A4) yields

$$\begin{array}{ll}{\displaystyle }\mathrm{0}& {\displaystyle }=q\left(T_{k-\mathrm{1}}-T_k\right)+A_k{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{T}}^{*}}{\mathit{\nu }\mathit{\lambda }}}\left(a{S}_k+b-c{p}_k-T_k\right)\mathit{\nu }\mathit{\lambda },\\ {\displaystyle }\mathrm{0}& {\displaystyle }=q\left(S_{k-\mathrm{1}}-S_k\right)+A_k{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{T}}^{*}}{\mathit{\nu }\mathit{\lambda }}}\left(a{S}_k+b-c{p}_k-T_k\right)S_k.\end{array}$$

Replacing $x=T_{k-\mathrm{1}}-T_k$, $y=S_{k-\mathrm{1}}-S_k$, $T^{*}=a{S}_{k-\mathrm{1}}+b-c{p}_k-T_{k-\mathrm{1}}$, $g_{\mathrm{1}}=A_k{\mathit{\gamma }}_{\mathrm{T}}^{*}$ and $g_{\mathrm{2}}=\frac{g_{\mathrm{1}}}{\mathit{\nu }\mathit{\lambda }}$, we obtain

$$\begin{array}{}\text{(A6)}& {\displaystyle }\mathrm{0}& {\displaystyle }=qx+g_{\mathrm{1}}\left(T^{*}+x-ay\right),\text{(A7)}& {\displaystyle }\mathrm{0}& {\displaystyle }=qy+g_{\mathrm{2}}\left(S_{k-\mathrm{1}}-y\right)\left(T^{*}+x-ay\right).\end{array}$$

We simplify the previous equations as follows. Rewriting Eq. (A6) as

$$\left(T^{*}+x-ay\right)={\displaystyle \frac{-qx}{g_{\mathrm{1}}}}$$

and inserting it into Eq. (A7), we obtain

$$\begin{array}{ll}{\displaystyle }\mathrm{0}& {\displaystyle }=qy+g_{\mathrm{2}}(S_{k-\mathrm{1}}-y)\left(-{\displaystyle \frac{qx}{g_{\mathrm{1}}}}\right)=qy-qx{\displaystyle \frac{S_{k-\mathrm{1}}-y}{\mathit{\nu }\mathit{\lambda }}}\\ {\displaystyle }⟺\mathrm{0}& {\displaystyle }=\mathit{\nu }\mathit{\lambda }y-S_{k-\mathrm{1}}x+xy\\ {\displaystyle }⟺\mathrm{0}& {\displaystyle }=\left(\mathit{\nu }\mathit{\lambda }+x\right)y-S_{k-\mathrm{1}}x\\ {\displaystyle }\Rightarrow y& {\displaystyle }={\displaystyle \frac{S_{k-\mathrm{1}}x}{\mathit{\nu }\mathit{\lambda }+x}}.\end{array}$$

Note that we can divide the first line by q since, by the model assumptions, q>0. Because $x=T_{k-\mathrm{1}}-T_k\ll \mathit{\nu }\mathit{\lambda }\approx \mathrm{75}$ ^∘C, we may approximate

$$\begin{array}{}\text{(A8)}& {\displaystyle }y\approx {\displaystyle \frac{S_{k-\mathrm{1}}x}{\mathit{\nu }\mathit{\lambda }}}.\end{array}$$

Using this approximation, we may proceed to solve the system of equations. Since we also need to solve for the overturning q in box B₁, which is adjacent to the grounding line, a slightly different approach is needed than for the other boxes, as discussed in the next section.

Solution for box B₁

The overturning flux q is parameterized as

$$\begin{array}{}\text{(A9)}& {\displaystyle }q=C{\mathit{\rho }}_{*}\left(\mathit{\beta }\left(S_{\mathrm{0}}-S_{\mathrm{1}}\right)-\mathit{\alpha }\left(T_{\mathrm{0}}-T_{\mathrm{1}}\right)\right),\end{array}$$

in the model, see Sect. 2.1. Substituting this equation into Eqs. (A6) and (A7), we obtain

$$\begin{array}{}\text{(A10)}& {\displaystyle }\mathrm{0}& {\displaystyle }=\mathit{\alpha }{x}^{\mathrm{2}}-\mathit{\beta }xy-{\displaystyle \frac{g_{\mathrm{1}}}{C{\mathit{\rho }}_{*}}}\left(T^{*}+x-ay\right),\text{(A11)}& {\displaystyle }\mathrm{0}& {\displaystyle }=-\mathit{\beta }{y}^{\mathrm{2}}+\mathit{\alpha }xy-{\displaystyle \frac{g_{\mathrm{2}}}{C{\mathit{\rho }}_{*}}}\left(S_{\mathrm{0}}-y\right)\left(T^{*}+x-ay\right).\end{array}$$

Inserting the approximation for y from Eq. (A8) into the Eq. (A10), we obtain a quadratic equation for x,

$$\left(\mathit{\beta }s-\mathit{\alpha }\right)x^{\mathrm{2}}+{\displaystyle \frac{g_{\mathrm{1}}}{C{\mathit{\rho }}_{*}}}\left(T^{*}+x\left(\mathrm{1}-as\right)\right)=\mathrm{0},$$

with $s=S_{\mathrm{0}}/\mathit{\nu }\mathit{\lambda }$. Since $as=-\mathrm{0.057}\times S_{\mathrm{0}}/\mathrm{74.76}=-\mathrm{0.000762}\times S_{\mathrm{0}}\ll \mathrm{1}$, we can omit the last part of the last term,

$$\left(\mathit{\beta }s-\mathit{\alpha }\right)x^{\mathrm{2}}+{\displaystyle \frac{g_{\mathrm{1}}}{C{\mathit{\rho }}_{*}}}\left(T^{*}+x\right)=\mathrm{0}.$$

Rearranging (assuming that $\mathit{\beta }s-\mathit{\alpha }>\mathrm{0}$, which we demonstrate below), we obtain

$$x^{\mathrm{2}}+{\displaystyle \frac{g_{\mathrm{1}}}{C{\mathit{\rho }}_{*}\left(\mathit{\beta }s-\mathit{\alpha }\right)}}x+{\displaystyle \frac{g_{\mathrm{1}}{T}^{*}}{C{\mathit{\rho }}_{*}\left(\mathit{\beta }s-\mathit{\alpha }\right)}}=\mathrm{0},$$

and hence we obtain the solution

$$\begin{array}{ll}{\displaystyle }x=& {\displaystyle }-{\displaystyle \frac{g_{\mathrm{1}}}{\mathrm{2}C{\mathit{\rho }}_{*}\left(\mathit{\beta }s-\mathit{\alpha }\right)}}\\ \text{(A12)}& {\displaystyle }& {\displaystyle }\pm \sqrt{{\left({\displaystyle \frac{g_{\mathrm{1}}}{\mathrm{2}C{\mathit{\rho }}_{*}\left(\mathit{\beta }s-\mathit{\alpha }\right)}}\right)}^{\mathrm{2}}-{\displaystyle \frac{g_{\mathrm{1}}{T}^{*}}{C{\mathit{\rho }}_{*}\left(\mathit{\beta }s-\mathit{\alpha }\right)}}}.\end{array}$$

The temperature in the box B₁ near the grounding line is supposed to be smaller than in the ocean box B₀, since, in general, melting will occur in box B₁ and hence T₁<T₀, or equivalently $x=T_{\mathrm{0}}-T_{\mathrm{1}}>\mathrm{0}$. Furthermore, we know that $g_{\mathrm{1}}/\left(C{\mathit{\rho }}_{*}\right)=A_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{T}}^{*}/\left(C{\mathit{\rho }}_{*}\right)$ is positive, as all factors are positive. Since $\mathit{\alpha }=\mathrm{7.5}\times {\mathrm{10}}^{-\mathrm{5}}$, $\mathit{\beta }=\mathrm{7.7}\times {\mathrm{10}}^{-\mathrm{4}}$ and $s=S_{\mathrm{0}}/\left(\mathit{\nu }\mathit{\lambda }\right)=S_{\mathrm{0}}/\mathrm{74.76}\ge \mathrm{0.4}$, it follows that βs>α. This means that the first summand of Eq. (A12) is negative and the second (negative) solution can be excluded. From here, we use $T_{\mathrm{1}}=T_{\mathrm{0}}+x$ and $y=x{S}_{\mathrm{0}}/\left(\mathit{\nu }\mathit{\lambda }\right)$ to solve for T₁, S₁, m₁ and q.

Solution for box B_k, k>1

Now, we give the solution for the other boxes B_k with k>1. By inserting the approximation for y in Eq. (A8) into Eq. (A6), we can solve for x as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{0}=qx+g_{\mathrm{1}}\left(T^{*}+x-a{\displaystyle \frac{S_{k-\mathrm{1}}x}{\mathit{\nu }\mathit{\lambda }}}\right)\\ {\displaystyle }⟺& {\displaystyle }\mathrm{0}=qx+g_{\mathrm{1}}{T}^{*}+g_{\mathrm{1}}x-g_{\mathrm{2}}a{S}_{k-\mathrm{1}}x\\ {\displaystyle }⟺& {\displaystyle }-g_{\mathrm{1}}{T}^{*}=x\left(q+g_{\mathrm{1}}-g_{\mathrm{2}}{S}_{k-\mathrm{1}}a\right)\\ \text{(A13)}& {\displaystyle }⟺& {\displaystyle }x={\displaystyle \frac{-g_{\mathrm{1}}{T}^{*}}{q+g_{\mathrm{1}}-g_{\mathrm{2}}a{S}_{k-\mathrm{1}}}}.\end{array}$$

The denominator is positive, as all terms are positive, and the sign of the numerator depends on T^*. The equation can now be solved for T_k, and then Eq. (A8) for S_k and Eq. (13) for m_k.

Appendix B: Motivation for geometric rule

Here, we want to motivate the rule that determines the extent of the boxes under each ice shelf. The rule aims at equal areas for all boxes. Assuming a half-circle with radius r=1, we want to split it into a fixed number n of equal-area rings. Generalized to the individual shapes of ice shelves, we will define the “radius” of an ice shelf as $r=\mathrm{1}-d_{\mathrm{GL}}/\left(d_{\mathrm{GL}}+d_{\mathrm{IF}}\right)$. We define r₁=1 the outer (grounding-line ward) radius of the half-circle ring covering an area A₁ and corresponding to box B₁ adjacent to the grounding line, r₂ as the outer radius of second outer-most half-ring, etc. The box B_k is then given by all shelf cells with horizontal coordinates (x, y) such that $r_{k+\mathrm{1}}\le r\left(x,y\right)\le r_k$ where $r_{n+\mathrm{1}}=\mathrm{0}$ is the center point of the circle. We can use these to determine the areas $A_n=\mathrm{0.5}\mathit{\pi }{r}_n^{\mathrm{2}},A_{n-\mathrm{1}}=\mathrm{0.5}\mathit{\pi }\left(r_{n-\mathrm{1}}^{\mathrm{2}}-r_n^{\mathrm{2}}\right),\mathrm{\dots },A_{n-k}=\mathrm{0.5}\mathit{\pi }\left(r_{n-k}^{\mathrm{2}}-r_{n-k+\mathrm{1}}^{\mathrm{2}}\right)$. If we require that $A_{\mathrm{1}}=A_{\mathrm{2}}=\mathrm{\dots }=A_n$, then, solving progressively, $r_{n-k}=\sqrt{k+\mathrm{1}}{r}_n$. By our assumption is r₁=1, hence $\mathrm{1}=r_{n-\left(n-\mathrm{1}\right)}=\sqrt{n}{r}_n$. This implies that $r_n=\mathrm{1}/\sqrt{n}$ and thus $r_{n-k}=\sqrt{\frac{k+\mathrm{1}}{n}}$. Hence, the box B_k for $k=\mathrm{1},\mathrm{\dots },n$ is defined as $\mathrm{1}-\sqrt{\left(n-k+\mathrm{1}\right)/n}\le d_{\mathrm{GL}}/\left(d_{\mathrm{GL}}+d_{\mathrm{IF}}\right)\le \mathrm{1}-\sqrt{\left(n-k\right)/n}$.