2.1 Plume model

The parametrization used in this study is based on the plume model developed by Jenkins (1991). Here we summarize the key assumptions and physics behind this model. The ice-shelf cavity is modelled by a two-dimensional geometry (Fig. 1), in which the ice-shelf base has a (local) slope given by the angle α. This geometry is assumed to be uniform in the direction perpendicular to the plane and constant in time and can be seen as a vertical cross section along a flow line of the ice shelf. We can define a coordinate X along the ice-shelf base with slope grounding line (X = 0) and moving up along the ice-shelf base due to positive buoyancy with respect to the ambient ocean water.

The situation depicted in Fig. 1 essentially yields a two-layer system of the meltwater plume with varying thickness D, velocity U, temperature T and salinity S lying above the ambient ocean with temperature T_a and salinity S_a. As explained in Jenkins (1991), the typically small values of the slope angle α allow us to consider conservation of mass, momentum, heat and salt within the plume in a depth-averaged sense. Moreover, as the plume travels upward in the direction of X, it is affected by entrainment (at rate $\dot{e}$) of ambient ocean water, as well as the fluxes of meltwater (with melt rate $\dot{m}$) and heat at the ice–ocean interface (with temperature T_b and salinity S_b). These considerations yield the following quasi-one-dimensional system of equations for (D, U, T, S) as a function of the coordinate X along the shelf base, denoting the balance of mass, momentum, heat and salt within the plume:

where g is the gravitational acceleration, C_d is the (constant) drag coefficient, Δρ = ρ_a − ρ is the difference in density between plume and ambient ocean, and $C_{\mathrm{d}}^{\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_{\mathrm{T}}$ and $C_{\mathrm{d}}^{\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_{\mathrm{S}}$ are the turbulent exchange coefficients (Stanton numbers) of heat and salinity at the ice–ocean interface. The above formulation makes explicit the linear dependence of the turbulent exchange velocities on the ocean current (γ_T = $C_{\mathrm{d}}^{\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_{\mathrm{T}}U$, γ_S = $C_{\mathrm{d}}^{\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_{\mathrm{S}}U$). The system of Eq. (2) is closed using suitable expressions for the entrainment rate $\dot{e}$, an equation of state ρ = ρ(T, S), the balance of heat and salt at the ice–ocean interface, and the liquidus condition. The expression for the entrainment rate is assumed to have the following form (Bo Pedersen, 1980):

with E₀ a dimensionless constant. Hence, the entrainment rate increases linearly with the plume velocity, is zero for a horizontal ice-shelf base and grows with increasing slope angle. Furthermore, a linearized equation of state yields

where β_S is the haline contraction coefficient and β_T the thermal expansion coefficient. The boundary conditions at the ice–ocean interface are given by

i.e. the first equation balances the turbulent exchange of heat with heat conduction and latent heat of fusion L in the ice, where c_w and c_i are the specific heat capacities of ocean water and ice, respectively, and T_i is the ice temperature. Similarly, Eq. (5b) is a balance between the turbulent exchange of salt and diffusion into the ice. Equation (5c) is the (linearized) liquidus condition that puts the interface temperature equal to the pressure freezing point at the local depth z_b of the ice-shelf base, which is equivalent to Eq. (1b).

Equations (2)–(5) form a closed set that can be solved to obtain the prognostic variables (D, U, T, S) of the plume as a function of the plume path X, given the ice-shelf draught z_b(X) with slope angle α(X), the ambient ocean properties T_a(z) and S_a(z), and the ice properties T_i and S_i. Of particular interest for the current work, however, are the ice–ocean interface conditions (Eq. 5), which essentially determine the melt rate $\dot{m}$, the key quantity of this study. In other words, the melt rate is determined by the fluxes of heat and salt at the interface, which in turn are linked to the development of the plume. Note that these boundary conditions can be simplified (McPhee, 1992; McPhee et al., 1999) to only two equations containing the freezing temperature T_f of the plume, rather than the interface properties T_b and S_b:

where $C_{\mathrm{d}}^{\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_{\mathrm{TS}}$ is an effective heat exchange coefficient. This simplified formulation can be used together with the prognostic Eq. (2) by substituting T_b with T_f in Eq. (2c) (note that T_b and T_f are not necessarily equal), whereas S_b disappears from the problem by substituting Eq. (5b) in Eq. (2d). Strictly speaking, Eq. (6) is only valid after assuming a constant ratio Γ_T∕Γ_S of the exchange coefficients, as explained by Jenkins et al. (2010), who also show that both Eqs. (5) and (6) give similar results when used to describe basal melt rates under Ronne Ice Shelf. Also note the similarity between Eq. (6) and the simple melt model described by Eq. (1), with the difference being the inclusion of heat conduction and the parametrization γ_T = $C_{\mathrm{d}}^{\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_{\mathrm{TS}}U$ as well as the plume variables T and S instead of ambient ocean properties. Hence, the turbulent exchange in this model is directly determined by the plume velocity that appears as a prognostic variable.

Without giving further details, we mention that the plume model described above can be evaluated for different ice-shelf geometries (i.e. vertical cross sections along flow lines) and different vertical temperature and salinity profiles of the ambient ocean (Jenkins, 2011, 2014). In this model, the general physical mechanism governing the development of the plume is the addition of meltwater at the ice–ocean interface, which increases its buoyancy. Changes in buoyancy affect plume speed, and that, combined with its temperature and salinity, determines the subsequent input of meltwater.

2.2 Basal melt parametrization along a flow line

Evaluating the aforementioned plume model for different geometries and ocean properties leads to a wide variety of solutions for the basal melt rates. The question arises whether there exists an appropriate scaling with external parameters that combines these results into a universal melt pattern. Here we will summarize how such a scaling can be found, leading to the basal melt parametrization of Jenkins (2014) for the quasi-one-dimensional geometries along flow lines described in the previous section; more details can be found in Appendix A. It is important to note that the following derivation is based on simple geometries with a constant basal slope and constant ambient ocean properties, though the resulting parametrization can easily be applied to more general cases, as shown in Sect. 3.1. Section 2.3 will discuss the extension of this parametrization to more realistic two-dimensional geometries.

Figure 2Dimensionless melt curve $\widehat{M}\left(\widehat{X}\right)$ used in the basal melt parametrization. Higher melt rates typically occur close to the grounding line with a maximum at $\widehat{X}$ ≈ 0.2. A transition from melting (red) to refreezing (blue) may occur further away from the grounding line, depending on the position of the ice front. Note that the value of $\widehat{X}$ depends on the distance to the grounding line, as well as the temperature difference T_a − T_f and the local slope α (see Appendix A). In other words, $\widehat{X}$ = 0 corresponds to the grounding line, but the dimensionless position of the ice-shelf front depends on the length scale and is not necessarily equal to $\widehat{X}$ = 1.

Download

The basal melt parametrization used in this study consists of a general expression for a dimensionless melt rate $\widehat{M}$ as a function of the dimensionless coordinate $\widehat{X}$ measured from the grounding line (Fig. 2). This dimensionless coordinate is essentially the vertical distance of the ice-shelf base from the grounding line, scaled by a temperature- and geometry-dependent length scale l:

where z_gl is the grounding-line depth and f(α) a slope-dependent factor. Hence, $\widehat{X}$ = 0 corresponds to the grounding line and any shelf point downstream from the grounding line corresponds to a value 0 < $\widehat{X}$ < 1 depending on T_a, S_a, z_gl and α. This scaling also implies that the ice-shelf front is not necessarily located at $\widehat{X}$ = 1, but its location is highly dependent on the input variables. Similarly, the melt rate is scaled as follows:

with a different slope-dependent factor g(α) and a constant parameter M₀. The dimensionless curve $\widehat{M}\left(\widehat{X}\right)$ in Fig. 2 is now defined by polynomial coefficients that were found empirically from the plume model results (Jenkins, 2014; Appendix A). In summary, to obtain the basal melt rate $\dot{m}$ at any point beneath the ice shelf, one requires the local depth z_b, local slope α, grounding-line depth z_gl, and ambient ocean properties T_a and S_a to calculate $\widehat{X}$ and find the corresponding value on the dimensionless curve $\widehat{M}\left(\widehat{X}\right)$, which then has to be multiplied by the physical scale given in Eq. (8) (see Appendix A for details). The physical quantities and constant parameters required for evaluating the parametrization are summarized in Table 1.

Table 1Physical quantities and constant parameters serving as input for the basal melt parametrization.

Download Print Version | Download XLSX

Although the scaling defined by Eqs. (7) and (8) is found in a purely empirical way, it is possible to derive the various factors analytically, as sketched in Appendix A. The empirical procedure and the physical meaning of the different factors are outlined in the following. A general solution to the problem is challenging to find as there are at least four length scales that determine the plume evolution (Jenkins, 2011). The first governing length scale is associated with the pressure dependence of the freezing point that imposes an external control on the relationship between plume temperature, plume salinity and the melt rate. Lane-Serff (1995) discussed how this length scale, (T_a − T_f)∕λ₃, approximately determines the distribution of melting and freezing beneath an ice shelf. Jenkins (2014) extended the analysis of Lane-Serff (1995) by making the transition point between melting and freezing dependent on the ice-shelf basal slope, resulting in the length scale Eq. (7) with slope factor f(α).

The second of these four length scales is associated with the ambient stratification, which determines how far the plume can rise before reaching a level of neutral buoyancy. Magorrian and Wells (2016) discuss the plume behaviour and resulting melt rates when this length scale dominates. Critically, with the pressure dependence of the freezing point assumed to be negligible, as required in the analysis of Magorrian and Wells (2016), no freezing can occur. The third length scale can be formulated by comparing the input of buoyancy from freshwater outflow at the grounding line with the input of buoyancy by melting at the ice–ocean interface (Jenkins, 2011). This length scale indicates the size of the zone next to the grounding line where the impact of ice-shelf melting on plume buoyancy can be ignored and conventional plume theory (Morton et al., 1956; Ellison and Turner, 1959) applied, and it is generally small compared with typical ice-shelf dimensions. The final length scale is that at which the Coriolis force takes over from friction as the primary force balancing the plume buoyancy in the momentum budget. Jenkins (2011) discussed these length scales in the context of which would take over as the dominant control on plume behaviour beyond the initial zone near the grounding line where the initial source of buoyancy dominates and showed the length scale associated with the pressure dependence of the freezing point, (T_a − T_f)∕λ₃, to be most important for typical ice-shelf conditions.

Hence, we obtain the second factor of the length scale l in Eq. (7) used in the parametrization. However, this length scale contains two more ingredients. First, as discussed by Jenkins (2011), the entrainment rate in the mass conservation Eq. (2a) explicitly depends on the slope α, whereas the melt rate is only affected indirectly, so there is a geometrical factor that scales the elevation of the plume temperature above the local freezing point:

This factor gives rise to the slope dependence f(α) in l, which is essentially an empirically derived scaling of the transition point between melting and freezing (Appendix A). The second ingredient is related to the coefficient Γ_TS, which appears in f(α) through the simplified interface conditions (Eq. 6). Jenkins (2014) retained the more complex melt formulation (Eq. 5) in the plume model while seeking empirical scalings based on an effective Γ_TS. As discussed by Holland and Jenkins (1999), the factor relating Γ_T and Γ_TS is itself a function of the plume temperature, so Jenkins (2014) expressed the effective Γ_TS as an empirical function of Γ_T, T_a − T_f and Eq. (9) including a constant initial value ${\mathrm{\Gamma }}_{{\mathrm{TS}}_{\mathrm{0}}}$ (see Appendix A). When distance along the plume path is scaled with this slightly more complex factor (see Eq. A10), the melt rates produced by the plume model conform to a universal form – first rising to a peak at $\widehat{X}$ ≈ 0.2 before falling and transitioning to freezing at $\widehat{X}$ ≈ 0.56 (Fig. 2).

With the distance along the plume path appropriately scaled, all that remains is to scale the amplitude of the melt-rate curves produced by the plume model and find the melt-rate scale M in Eq. (8). As in Jenkins (2011) the appropriate physical scales are (1) the temperature of the ambient ocean water relative to the freezing point; (2) the factor in Eq. (9) scaling the temperature elevation of the plume above freezing; (3) a factor that scales the plume speed, given by the ratio of plume buoyancy to frictional drag:

The second term in parenthesis is the factor that scales the plume temperature relative to the ambient temperature and thus controls plume buoyancy. It replaces the initial buoyancy flux at the grounding line used in the scaling of Jenkins (2011). The final expression includes factors and powers that are derived empirically (though some theoretical arguments can be applied, see Appendix A), giving rise to the form of M with slope factor g(α) in Eq. (8). In summary, the result of this scaling procedure is an approximately universal melt-rate curve, which can then be represented by a single polynomial expression that is accurate to about 20 % for melt rates ranging over many orders of magnitude (Jenkins, 2014).

Table 2Definition of the ice mask. The ice-shelf criterion is that for uniform ice with density ρ_i floating on ocean water with density ρ_w. The minimum ice thickness used here is H_i,min = 2 m.

Download Print Version | Download XLSX

2.3 Basal melt parametrization in 2-D: effective plume path

As explained in the previous section, an important feature of the basal melt parametrization is its dependence on non-local quantities, in particular the grounding-line depth z_gl from which the plume originated. Therefore, in order to apply the parametrization to realistic geometries, one needs to know for each ice-shelf point the corresponding grounding-line point(s) serving as the origin of the plume(s) reaching that particular shelf point. For the quasi-one-dimensional settings considered so far, this is not an issue, since the plume can only travel in one direction. However, for general ice-shelf cavities, an arbitrary shelf point can be reached by plumes from multiple directions, corresponding not only to different values for z_gl, but also to different slope angles α. This means that the plume parametrization cannot be directly applied to such geometries. An algorithm is needed to determine effective values for z_gl and α. The development of this algorithm is the main focus of the current work and discussed below.

Figure 3Schematic of the algorithm for finding the average grounding-line depth and associated slope angle used by the basal melt parametrization. (a) Top view of an ice shelf on a horizontal grid. The algorithm searches in 16 directions on the grid from the shelf point (i, j). Possible grounded points found in this way are marked by ×. (b) Vertical slice along the nth direction (e.g. the red dotted line in a). If the grounded point is higher than the previous shelf point, the grounding-line depth z_n is found by interpolation along the bed (z_n = $\frac{\mathrm{1}}{\mathrm{2}}(H_{\mathrm{bed},\mathrm{1}}$ + H_bed,2)). (c) Interpolation along the ice base if the grounded point in the nth direction is deeper than the previous shelf point (z_n = $\frac{\mathrm{1}}{\mathrm{2}}(z_{\mathrm{b},\mathrm{1}}$ + z_b,2)).

Download

As a starting point, we consider the usual topographic data in terms of two-dimensional fields for the ice thickness H_i, bedrock/seabed elevation H_bed and elevation of the upper ice surface H_s used by ice-dynamical models. The following algorithm is valid for any topographic data on a rectangular grid with any resolution Δx × Δy. First of all, the topographic data are used to define an ice mask based on the criterion for floating uniform ice, as shown in Table 2. Furthermore, the depth of the ice base is determined to be

In order to apply the basal melt parametrization to these two-dimensional data, effective values for z_gl and α must be determined for every ice-shelf point (i, j) with basal depth z_b(i, j), where the indices i and j denote the position on the grid. This is done by first searching for “valid” grounding-line points in 16 directions on the grid, starting from any shelf point (i,  j), as depicted in Fig. 3a. Note that we can calculate a local basal slope s_n(i,j) at the point (i,j) in the nth direction as follows:

where (i_n, j_n) denotes a direction vector on the grid, i.e. (i_n, j_n) = (1, 0) denotes right, (i_n, j_n) = (0, 1) denotes up, etc., and Δx and Δy denote the horizontal grid size in the x and y direction, respectively. To determine whether a grounding-line point found in 1 of the 16 directions is valid for the calculation of the basal melt, the following two criteria are applied:

1. Assuming that a buoyant meltwater plume can only reach the point (i, j) from the nth direction if the basal slope in that direction is positive, the algorithm only searches in directions for which s_n(i, j) > 0.

2. If the first criterion is met for the nth direction, the algorithm searches in this direction for the nearest ice-sheet point. More precisely, the associated direction vector (i_n, j_n) is added to the grid indices and the mask value in the resulting point is checked. This process is repeated until either an ice-sheet point, an ocean point or the domain boundary is encountered. An ice-sheet point found in this way is only considered to be a valid grounding-line point if it lies deeper than the original ice-shelf point at (i, j), assuming again that a buoyant meltwater plume from the grounding line can only go up. The second criterion then becomes z_n(i, j) < z_b(i, j), where z_n(i, j) is the grounding-line depth in the nth direction.

Note, however, that in determining the second criterion the depth difference between the encountered sheet point and the adjacent shelf point can be considerable, especially for coarser resolutions. In such cases, the algorithm tries to obtain a better estimate of the true grounding-line depth in this direction, z_n(i, j), by interpolating along either the bed or the ice base, as shown in Fig. 3b and c. The two cases shown in these figures account for either a positive or a negative basal slope beyond the grounding line. One should note that this additional step assumes the grounding line to be located halfway between the sheet and shelf points, which could be improved by more sophisticated interpolation techniques.

Following the above procedure yields for each ice-shelf point (i, j) a set of grounding-line depths z_n and local slopes s_n in the directions that are valid according to the aforementioned two criteria. Keep in mind that not all directions may yield a (valid) grounding-line point – in particular those towards the open ocean. Now, in order to determine the effective grounding-line depth z_gl(i, j) and effective slope angle α(i, j) necessary for calculating the basal melt in the shelf point (i, j), we simply take the average of the values found for z_n and s_n:

where N_ij denotes the number of valid directions found for the shelf point (i, j). On the other hand, if no valid values for z_n and s_n are found for a particular shelf point, we take z_gl = z_b and α = 0, leading to zero basal melt in that point (see Appendix A).

In summary, the method described above yields two-dimensional fields for the effective grounding-line depth z_gl and effective slope tan(α), given topographic data in terms of H_i, H_s, and H_bed and a suitable ice mask, such as the one defined in Table 2. These fields, in turn, serve as input for the basal melt parametrization described in the previous section, together with appropriate data for the ocean temperature T_a and salinity S_a (discussed in Sect. 3.2). We thus obtain a complete method for calculating the basal melt for all Antarctic ice shelves, given the topography and ocean properties, which can also be used in conjunction with ice-dynamical models. In the following, however, we use the Bedmap2 dataset (Fretwell et al., 2013) to define the present-day topography of Antarctica and disregard the ice dynamics. More specifically, the original Bedmap2 data are remapped to a rectangular grid with grid size Δx = Δy = 20 km, using the mapping package OBLIMAP 2.0 (Reerink et al., 2016). The resulting topographic data can be used as input for the algorithm described here, leading to the fields for z_gl and tan(α) shown in Fig. 4, which are used for the basal melt calculations discussed in Sect. 3. Note that the mask in Fig. 4a does not exactly match the Bedmap2 mask because a constant ρ_i was used in formulation of Table 2 as is common in many ice-sheet models. This might cause discrepancies in the position of the grounding line, which, however, are likely compensated for by the rather coarse resolution. In Fig. 4b one can see that the lowest values of z_gl are obtained towards the inland regions of Filchner–Ronne Ice Shelf (FRIS) and Amery Ice Shelf. The values for the local slope are typically high near the grounding line and in some places also near the ice front, as shown in Fig. 4c.

One should note that, although we attempt to directly translate the concept of a quasi-1-D plume to a multitude of plumes in two dimensions, there are important physical effects not taken into account by this approach. Most importantly, a realistic two-dimensional plume has an additional degree of freedom because it also develops in the cross-flow direction, causing the width to be a dynamic variable in addition to the thickness D. This can have significant consequences for the mass budget currently described by Eq. (2a). Hattermann (2012) explored the possibility of adding a variable plume width to the original plume model and Hattermann et al. (2014) showed that such a 2-D formulation improves the prediction of melt rates for a realistic ice-shelf geometry compared to the 1-D model. Although this appears to be an important extension of the plume model that should be taken into account, the aim of the current work is to explore the capabilities of the original 1-D plume parametrization in predicting melt rates around Antarctica. The current approach is meant to be a simple method to parametrize the net circulation within an ice-shelf cavity as the average effect of multiple plumes, in order to be applied around the entire ice sheet. Further extensions for obtaining a 2-D plume model are beyond the scope of this work.

3.1 Comparison of basal melt parametrizations along flow lines

Topographic data along flow lines for both Filchner–Ronne Ice Shelf and Ross Ice Shelf are taken from Bombosch and Jenkins (1995) and Shabtaie and Bentley (1987), respectively. These data can be used to determine the quantities z_b, α and z_gl necessary for calculating the basal melt with the parametrization of Sect. 2.2. Furthermore, we define a uniform ambient ocean temperature T_a = −1.9 ^∘C + ΔT, where ΔT is varied between runs, and a constant ambient ocean salinity S_a = 34.65 psu. The results of these calculations are shown in Fig. 5 and compared with those of the full plume model described in Sect. 2.1. Moreover, we compare with two simple basal melt parametrizations based on Eq. (1), namely the linear (i.e. in T_a − T_f) parametrization by Beckmann and Goosse (2003) with constant γ_T and the quadratic parametrization by DeConto and Pollard (2016) with γ_T = κ_T|T_a − T_f|. Apart from the values listed in Table 1, additional model parameters used for these calculations are given in Table 3.

Table 3Additional model parameters used for evaluating the plume model and the simple parametrizations described in Sect. 3.1. BG2003 refers to Beckmann and Goosse (2003) and DCP2016 refers to DeConto and Pollard (2016).

Download Print Version | Download XLSX

Figure 4Effective plume paths under the Antarctic ice shelves as calculated by the algorithm of Sect. 2.3 using the Bedmap2 topographic data remapped on a 20 km by 20 km grid. (a) Ice mask according to Table 2. (b) The effective grounding-line depth z_gl. (c) The effective slope tan(α). (d) The difference between local ice-base depth and associated grounding-line depth, z_b − z_gl.

Download

Figure 5Comparison of the plume model (Sect. 2.1) with the 1-D basal melt parametrization (Sect. 2.2), as well as the parametrizations of Beckmann and Goosse (2003) (BG2003) and DeConto and Pollard (2016) (DCP2016), for flow lines along Filchner–Ronne Ice Shelf (a, c, e, g) and Ross Ice Shelf (b, d, f, h), both with uniform ocean temperature T_a = −1.9 ^∘C + ΔT and constant salinity S_a = 34.65 psu. (a, b) Geometry of the ice-shelf base. (c, d) Melt pattern for ΔT = 0 ^∘C. (e, f) Melt pattern for ΔT = 0.8 ^∘C. (g, h) Melt-rate average along the flow line as a function of ΔT. Note that the black curve is nearly identical to the green curve and might appear below it. Also note the difference in vertical scale between the left and right columns. The flow-line locations are indicated in Fig. 6.

Download

Figure 5 shows that both the current parametrization and the original plume model yield approximately the same melt-rate patterns as a function of the horizontal distance from the grounding line. These patterns roughly correspond to the dimensionless melt curve in Fig. 2, i.e. maximum melt near the grounding line and possibly refreezing further away along the flow line. This is most apparent in Fig. 5c, which shows a transition from melting to freezing, since the relatively deep draught of FRIS allows higher values of the dimensionless coordinate $\widehat{X}$. On the other hand, Fig. 5d does not show refreezing because the draught of Ross Ice Shelf is much shallower. Increasing the ocean temperature (through ΔT) can significantly enhance basal melt and remove the area of refreezing, as shown in Fig. 5e and f. In these cases, additional melt peaks occur in regions of high basal slope. Moreover, although the general agreement is good, the discrepancies between the current parametrization and the plume model are largest when the basal slope changes rapidly, because the parametrization responds immediately to the change while the full model has an inherent lag as the plume adjusts to the new conditions. On the whole, we see that the melt patterns given by the plume parametrization can be quite complex, while the two simple parametrizations give nearly constant curves (i.e. independent of the position with respect to the grounding line).

It is interesting to investigate the temperature sensitivity of the four models in terms of the horizontally averaged melt rate as a function of ΔT, as shown in Fig. 5g and h. In the case of FRIS, the plume model and parametrization are much more sensitive to the ocean temperature than the two simpler models. However, the average melt rates for Ross Ice Shelf are rather similar for all four models and all values of ΔT. Hence, the difference in the temperature sensitivity depends significantly on the ice-shelf geometry, where the plume parametrization appears to have a larger potential for capturing diverse melt values than the simpler models. Note that in both cases the temperature dependence of the plume parametrization is slightly non-linear, similar to the DeConto and Pollard (2016) parametrization, while the Beckmann and Goosse (2003) parametrization has a linear temperature dependence. Following the discussion of Holland et al. (2008), the temperature dependence of the plume parametrization should therefore be more realistic than the one of Beckmann and Goosse (2003). However, the quadratic parametrization of DeConto and Pollard (2016) tends to significantly underestimate the melt rates as well, despite its non-linearity. It appears that the geometry dependence of the plume parametrization is an important factor for the temperature sensitivity of the calculated basal melt rates. In Sect.3.3 we show that these geometrical effects are indeed crucial for obtaining realistic melt rates with the 2-D parametrization, but first we discuss the matter of determining a suitable input field for the ocean temperature.

3.2 Effective ocean temperature

The previous section dealt with the 1-D basal melt parametrization along a flow line using a uniform ambient ocean temperature for the entire ice-shelf cavity. While a uniform temperature might appear a reasonable first approximation for a single ice shelf, it is far from realistic to apply a single ocean temperature for multiple ice shelves around the entire Antarctic continent. Therefore, in order to apply the parametrization to the 2-D geometry defined by Fig. 4, a suitable 2-D field for the ocean temperature T_a is required. In principle, the same is true for the salinity S_a, but we will assume that the horizontal variations in ocean salinity around Antarctica are so small that the pressure freezing point T_f is only affected by variations in depth. In the following, we will therefore take a uniform salinity S_a = 34.6 psu. One should realize that vertical variations in S_a, which are not accounted for in the current parametrization, would be important in reality, as discussed in Sect. 4.

Two problems arise when considering a 2-D ocean temperature field for forcing the parametrization. First of all, such a field should ideally be based on observational data, but ocean temperature measurements in the Antarctic ice-shelf cavities are sparse. A more feasible approach would be to compute an interpolated field based on ocean temperature data in the surrounding ocean, which inevitably contains artefacts resulting from the non-uniform and predominantly summertime sampling. Secondly, even if a complete dataset of ocean temperatures were available, it would not be immediately clear which temperatures (i.e. at which depth) are characteristic of the ocean water reaching the grounding lines (e.g. Jenkins et al., 2010). In principle, detailed knowledge of the bottom topography and the ocean circulation would be required for this, which goes beyond the scope of the current modelling approach.

In view of these issues, we construct an effective ocean temperature field with which the current plume parametrization yields melt rates that are as close as possible to present-day observations, averaged over entire ice shelves. In other words, this can be regarded as the inverse problem of computing the unknown ocean temperatures by assuming that the model output for the melt rates matches the (averaged) observations. For this purpose, we use the results of Rignot et al. (2013), who calculated the area-averaged melt rates for each Antarctic ice shelf, based on a combination of observational data and regional climate model output for the different terms in the local ice-shelf mass balance. Other datasets for recent Antarctic basal melt rates exist (e.g. Depoorter et al., 2013), as well as more recent data for ice-shelf thinning (Paolo et al., 2015) from which the basal melt rates can be calculated when combined with the other terms in the mass balance (e.g. velocity and surface melt rates). These alternative datasets for the (area-averaged) basal melt rates are expected to be at least of the same order of magnitude, which we deem sufficient for the purpose of the current study. Since it is impossible to resolve each individual ice shelf from the Rignot et al. (2013) dataset with the currently used 20 km resolution (Fig. 4), we consider a set of 13 ice-shelf groups and determine the area-averaged basal melt for each group from the data of Rignot et al. (2013). The definition of these groups and the calculated average melt rates are shown in Fig. 6. Note that the shelves have been grouped based on their geographical location but also for more practical reasons such as the possibility of distinguishing their boundaries on the 20 km grid.

Figure 6The 13 groups of ice shelves used for constructing the effective ocean temperature field. Average melt rates and error estimates (1 SD – 1 standard deviation) for each group are calculated from the data of Rignot et al. (2013) for individual ice shelves. Green lines indicate the approximate positions of the flow lines used in Fig. 5.

Download

As a starting point for constructing the effective ocean temperature, we consider the observational data of the World Ocean Atlas 2013 (WOA13, Locarnini et al., 2013), which contains a global dataset of (annual mean) ocean temperatures within a range of depths (0–5500 m). Restricting ourselves to the temperature data for latitudes south of 60^∘ S, we average the ocean temperatures over depth intervals [z₁, z₂], where z₁ is the level of the bed (i.e. the deepest level for which data are available) with the additional constraint z₁ ≥ −1000 m and z₂ = min{0, z₁ + 400 m}. This results in a relatively smooth 2-D temperature field containing an inherent dependence on the bottom topography, which can be considered a first estimate for the ocean water flowing into the ice-shelf cavities. The depth-averaged temperature field is now remapped on the same 20 km grid as the topography data (see Sect. 2.3 and Fig. 4) and interpolated using natural-neighbour interpolation (i.e. a weighted version of nearest-neighbour interpolation, giving smoother results) to obtain data in the entire domain of interest. The resulting temperature field, called T₀, is shown in Fig. 7a. One should note that both the depth-averaging and the interpolation procedures introduce biases in the resulting field. In particular, the rather simple interpolation technique also interpolates ocean temperatures between ice-shelf cavities separated by the continent or grounded ice, which is not realistic as it propagates temperatures into cavities that the corresponding ocean water cannot reach. Using the natural-neighbour interpolation method appears to limit these effects. However, the details of the resulting field T₀ are somewhat arbitrary as it needs to be adjusted in order to obtain melt rates that agree with the data of Rignot et al. (2013).

Figure 7(a) Depth-averaged and interpolated ocean temperature, T₀, calculated from annual mean WOA13 data. (b) Effective ocean temperature T_eff = max{T₀ + ΔT, −1.9} constructed from T₀ as described in Sect. 3.2. The circles indicate the positions of the sample points in which the values of ΔT are imposed. The colour of each circle corresponds to the imposed value of ΔT (same colour scale), ranging from −1.4 to 0.8 ^∘C. The full ΔT field is obtained by linearly interpolating these values.

Download

The aim is now to modify this depth-averaged, interpolated temperature field T₀ in such a way that the basal melt parametrization yields melt rates close to those shown in Fig. 6 for each ice-shelf group. As explained earlier, this modification is necessary for eliminating biases in T₀ caused by the sparse observations and numerical interpolation and also because the flow dynamics of the ocean are not resolved. The field T₀ is now modified by adding a 2-D field of temperature differences (ΔT), which, in turn, is the result of linearly interpolating individual values of ΔT in 29 carefully chosen sample points, with ΔT = 0 on the domain boundary. The sample points and values of ΔT have been determined by trial and error and are certainly not a unique nor optimal configuration. The points are mainly located in regions that are most affected by interpolation between strictly separated cavities (e.g. grounding line of FRIS) or extrapolation of warm open-ocean temperatures into cavities (e.g. Dronning Maud Land, shelf groups 2 and 3 in Fig. 6). The resulting effective temperature field, T_eff = T₀ + ΔT, is shown in Fig. 7b, which also indicates the positions of the aforementioned sample points along with the used values of ΔT in these points. Note that, for technical reasons explained in Appendix A, we have applied a lower limit to the effective temperature equal to the pressure freezing point at surface level. With the current choice S_a = 34.6 psu, this implies T_eff ≥ −1.9 ^∘C. Comparing Fig. 7a and b, we see that the main effect of ΔT is a decrease in the ocean temperature over most of the continental shelf and most ice-shelf cavities (in particular for Ross and Amery ice shelves) and a slight increase in the ocean temperature in West Antarctica and some regions in East Antarctica (e.g. shelf group 6 in Fig. 6). Again, note that the details in the procedure for calculating T₀ and ΔT are somewhat arbitrary, since increasing one term would require decreasing the other term in order to obtain similar values for T_eff with similar basal melt rates.

Figure 8 shows the basal melt rates computed by the parametrization using the effective temperature T_eff of Fig. 7b as forcing. An area-averaged value is obtained for each of the 13 ice-shelf groups in Fig. 6 and compared with the observational values from the Rignot et al. (2013) data. By construction, the modelled basal melt rates correspond closely to the observational values and fall within the error estimates. A notable exception is the value for Filchner–Ronne Ice Shelf, which is 0.32 ± 0.08 m yr⁻¹ according to the observations, whereas the parametrization gives a value just above 0.5 m yr⁻¹. This discrepancy is caused by the lower bound of −1.9 ^∘C imposed on the effective temperature, whereas in reality the temperatures can reach values below −2.0 ^∘C (e.g. Nicholls et al., 2009). As we can see in Fig. 7b, the ocean water below FRIS is almost entirely at this minimum temperature, making it impossible to further improve the basal melt rate without using unfeasibly low values for T_eff. This rather technical constraint might be relaxed in various ways, as briefly discussed in Appendix A, possibly improving the melt rates in very cold cavities.

Figure 8Area-averaged basal melt rates for each ice-shelf group in Fig. 6 obtained with the plume parametrization and the effective temperature field of Fig. 7b. The modelled melt rates are plotted against the averaged observational values given in Fig. 6. For four important shelf groups, the data points are explicitly labelled along with the corresponding group number in Fig. 6. The horizontal error bar is 1 standard deviation uncertainty in the observations.

Download

Nevertheless, the plume parametrization in conjunction with the constructed effective temperature field appears to yield realistic present-day melt rates for all shelf groups. By construction, the effective temperature shown in Fig. 7b contains an inherent dependence on the bottom topography, with typically lower temperatures above the continental shelves (and thus in the ice-shelf cavities), while still retaining the spatial variation in temperature of the surrounding deep ocean (e.g. higher temperatures for West Antarctica, leading to higher melt rates for ice-shelf groups 11 and 12 as defined in Fig. 6).

Figure 9Basal melt rates in metres per year with the Bedmap2 topographic data and the effective temperature field of Fig. 7b as obtained from (a) the plume parametrization with additional input parameters from Fig. 4 and (b) the quadratic parametrization of DeConto and Pollard (2016).

Download

3.3 Comparison of 2-D melt-rate patterns

The effective grounding-line depth and effective slope in Fig. 4, the effective ocean temperature in Fig. 7b and the assumption S_a = 34.6 psu constitute the full set of input parameters necessary for evaluating the plume parametrization on the entire 2-D geometry. The resulting 2-D field of basal melt rates under all Antarctic ice shelves is shown in Fig. 9a (note that these are the same data used for the area-averaged melt rates in Fig. 8 but now plotted as a spatial field rather than averaged values over the ice shelves). A general pattern that can be observed, especially on the bigger ice shelves, consists of regions of higher melt close to the grounding line and lower melt or patches of refreezing closer to the ice front, with the latter being most apparent at the ice fronts of shelf groups 1, 2 and 9. This pattern is a consequence of the underlying plume model, as shown in Sect. 3.1 for data along a flow line. Moreover, the highest melt rates occur in West Antarctica (shelf groups 11 and 12) and some specific shelves in East Antarctica (shelf groups 6 and 7), where the constructed effective temperature is significantly higher than elsewhere. The general melt patterns within individual cavities appear to be in line with observations, e.g. Rignot et al. (2013). However, one should note that the Rignot et al. (2013) melt pattern shows a greater spatial variability, with more patches of (stronger) refreezing occurring between patches of melting (Fig. 11a). Especially beneath FRIS and Ross Ice Shelf, the melt pattern appears quite complex and local deviations from the general pattern can be considerable (Fig. 11b). These discrepancies in the current parametrization might have different reasons, such as the coarse resolution or the fact that we disregard the details of the ocean circulation within the ice-shelf cavities, as well as effects due to the Coriolis force and both seasonal and vertical variability in the temperature and salinity fields.

Figure 10As Fig. 9a, but with a logarithmic colour scale (negative and zero values shown white) and zoomed in on (a) Filcher–Ronne Ice Shelf (group 1), (b) West Antarctica including Pine Island and Thwaites (group 11) and (c) Ross Ice Shelf (group 9).

Download

Figure 11Basal melt rates in metres per year extracted from the Rignot et al. (2013) observational dataset (courtesy of Dr. Jeremie Mouginot): (a) raw data plotted together with the currently used mask; (b) difference between the plume parametrization (Fig. 9a) and the observations interpolated on the 20 km grid.

Download

Furthermore, Fig. 9 shows the melt-rate patterns of the plume parametrization zoomed in on three regions, giving more insight into the orders of magnitude of the highest melt rates. The high near-grounding-line melt rates for FRIS have values between 1 and 10 m yr⁻¹, while those for Ross Ice Shelf appear 1 order of magnitude smaller. On the other hand, the West Antarctic melt rates shown in Fig. 10b have values around 10 m yr⁻¹ or more due to the higher ocean temperatures here. It should be noted, however, that the latter values are still lower than those observed in the Rignot et al. (2013) data, where local melt rates close the grounding line can reach 100 m yr⁻¹, while the average melt rates over the full area of Pine Island and Thwaites are 16.2 and 17.7 m yr⁻¹, respectively.

For comparison, we also evaluate the quadratic parametrization of DeConto and Pollard (2016), described in Sect. 3.1, using the same geometric data and the effective temperature field of Fig. 7b as input. The resulting basal melt-rate pattern is shown in Fig. 9b. Comparing this figure to Fig. 9a shows that the quadratic parametrization yields significantly lower melt rates than the plume parametrization, at least with the current effective temperature as input. The only visible patches of basal melt are located in the aforementioned regions where the ocean temperature is high, as well as near the grounding line of Filchner–Ronne Ice Shelf. Therefore, if the effective temperature in Fig. 7b is indeed characteristic of the true temperatures in the ice-shelf cavities, the quadratic parametrization would require significant tuning in order to obtain a similar agreement with observed melt rates as currently found with the plume parametrization. For completeness, we mention that the linear parametrization of Beckmann and Goosse (2003) yields even lower melt rates due to its low temperature sensitivity, as discussed in Sect. 3.1.

To further clarify the differences between the two parametrizations in Fig. 9, we have repeated the steps outlined in Sect. 3.2 and constructed a second effective temperature field based on the quadratic parametrization by DeConto and Pollard (2016) instead of the plume parametrization. The resulting temperature field is shown in Fig. 12a. Note that the difference between this field and the one in Fig. 7b only lies in the values chosen for ΔT and not in the underlying interpolated observations (T₀). For simplicity, the ΔT values have been imposed on the same sample points as used for Fig. 7b. Comparing the two effective temperature fields in Figs. 7b and 12a shows that much higher ocean temperatures are required for the quadratic parametrization to give realistic area-averaged melt rates. The ΔT values imposed on the sample points indicated in Fig. 12a range from −0.5 to 5.4 ^∘C, while those used for Fig. 7b range from −1.4 to 0.8 ^∘C. Furthermore, we can calculate the root mean square values of T_eff − T₀ over the entire domain (disregarding the continental points), yielding 0.3 ^∘C for Fig. 7b and 1.1 ^∘C for Fig. 12a. Hence, the effective temperature in Fig. 7b lies closer to the underlying observational data T₀ than the field in Fig. 12a.

Figure 12(a) Effective temperature field constructed in a similar way as Fig. 7b, but with different values for ΔT (indicated by the circles and ranging from −0.5 to 5.4 ^∘C), chosen in order to match the melt rates of the quadratic parametrization of DeConto and Pollard (2016) with the data of Rignot et al. (2013). (b) Basal melt rates obtained with the quadratic parametrization of DeConto and Pollard (2016) using the Bedmap2 topographic data and the effective temperature in (a) as input.

Download

The basal melt rates resulting from the quadratic parametrization and the new effective temperature field are shown in Fig. 12b. Clearly, the higher ocean temperatures cause significantly higher melt rates than those shown in Fig. 9b. However, compared with the plume parametrization in Fig. 9a, the spatial distribution of these melt rates is more uniform, showing less prominent melt peaks near grounding lines and no patches of refreezing. It appears that the quadratic temperature dependence together with the (slight) depth dependence through the pressure freezing point T_f (Eq. 1b) is not sufficient for obtaining realistic melt rates without significantly increasing the input ocean temperature, which can be considered equivalent to using different tuning factors for different ice shelves. On the other hand, the plume parametrization, containing an additional geometry dependence through the grounding-line depth and local slope, appears to yield the required melt rates rather naturally with ocean temperatures constructed in a plausible way, and it results in a more realistic spatial pattern with highest basal melt rates near the grounding line as well as areas of refreezing.