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

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

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:

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

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

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$,

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

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

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

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:

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

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

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

Download

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 4Sensitivity 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.

Download

Figure 5Sub-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.

Download

Figure 6Sensitivity 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.

Download

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.