We diagnose the fluxes at the grounding line with the finite-element ice-flow model Úa (Gudmundsson, 2013). The flow model Úa has been used to calculate the ice-flow for various geometries involving ice-shelf buttressing (De Rydt and Gudmundsson, 2016; Gudmundsson et al., 2017; Royston and Gudmundsson, 2016), and results obtained by the model submitted to a number of model intercomparison experiments (Pattyn et al., 2012) and (Pattyn et al., 2013). The model employs an unstructured grid and hence allows for resolving the grounding-line zone locally with high resolution. Simultaneous inversion for the ice rate factor (A) and the basal slipperiness (C) can be done either over nodal or over element values, and using either Bayesian- or Tikhonov-type regularisation. The gradient of the resulting objective function is calculated using the adjoint method.

Here we use Úa to solve the shallow ice-stream equations (MacAyeal, 1989; Morland, 1987) in a diagnostic mode using a Weertman-type sliding law (see Eq. 7) and Glen's flow law (see Eq. 8). In the glaciological literature the shallow ice-stream equations are also referred to as the shallow shelf approximation or shelfy stream approximation and often abbreviated as SSA. In the two horizontal dimensional situation (2HD) the SSA momentum equations are

where

and R is the tensor of resistive stresses given by Eq. (15), h is the ice thickness, s is the ice surface elevation, ρ_i is the vertically averaged ice density, and τ_bh is the horizontal part of the bed-tangential basal traction τ_b. Where the ice is floating τ_bh=0. In the SSA the flotation criterion has the form h<h_f with

where S is the ocean surface, B the bedrock, and ρ_w is the ocean density. The flotation criterion in Úa is evaluated at each integration point of the elements of the finite element mesh and the basal drag term is evaluated accordingly through a standard finite-element procedure involving element-wise integration.

2.1 Methodology

Using the ice flow model Úa, we calculate ice velocities for the entire Antarctic Ice Sheet, including all ice shelves. The SSA equations are solved throughout the computational domain. Stress boundary conditions (i.e. Neumann boundary conditions) are applied at the margins of the computational domain. Since the modelling domain covers the whole of the AIS, no inflow or outflow boundary conditions (i.e. Dirichlet boundary conditions) need to be applied at any sections of the boundary.

Two different computational meshes were generated and the sensitivity of the results was evaluated using linear (3-node), quadratic (6-node) and cubic (10-node) triangular elements. All results presented here were obtained using a very high-resolution mesh generated with the finite-element mesh generator Gmsh (Geuzaine and Remacle, 2009) with 1 360 894 triangular linear elements and 689 042 nodes. Within 5 km distance to the grounding line, the mesh was refined such that element sizes decrease towards the grounding line to a maximum size of 250 m directly at the grounding line. Overall, the elements have a maximal size of 179 307 m in the interior of the continent and a minimal size of 56 m along the grounding line. The mean element size is 1596 m and the median is 480 m. A regional example of the mesh is given in Fig. S1 in the Supplement. The robustness of the results was also tested based on the mesh used in Reese et al. (2018), as discussed in Appendix B.

Ice thickness and bed geometry input is based on the Bedmap2 estimates (Fretwell et al., 2013). Vertically averaged ice densities were calculated using firn thickness fields from RACMO2 (Lenaerts et al., 2012) and assuming a constant ice density of 910 kg m⁻³ and a firn density of 500 kg m⁻³. Resulting densities range from 770 to 910 kg m⁻³ and the horizontal gradients in vertically averaged densities are hence small; see Fig. S2. In a few places the bathymetry around the grounding lines was vertically modified to improve its alignment with Bindschadler et al. (2011), with vertical adjustments of maximally 50 m being allowed.

For the entire Antarctic set-up we inverted for basal slipperiness C (see Eq. 7) and ice softness fields A (see Eq. 8) to match observed 2015–2016 velocities derived from Landsat 8 imagery (Gardner et al., 2018). The stress exponent of Glen's flow law was set to n=3 and we repeated the inversion for a whole sequence of sliding law exponents m=1, 2, 3, 4, 5, 7, 9, 11. We inverted for A and C over the computational nodes using Tikhonov-type regularisation. The inversion procedure minimizes the function

$${\displaystyle }{\displaystyle }J(u,p)=I\left(u\right)+R\left(p\right),$$

with respect to p, where p stands for model parameters to be determined (i.e. A and C), u are modelled surface velocities, I is the data misfit function, and R is the regularisation term. The misfit function I has the form

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }I={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathcal{A}}}\int {\left({\mathit{v}}_{\mathrm{modelled}}-{\mathit{v}}_{\mathrm{observed}}\right)}^{\mathrm{2}}/e^{\mathrm{2}}\mathrm{d}A,\end{array}$$

where $\mathcal{A}=\int \mathrm{d}A$ is the total area, v_modelled and v_observed are modelled and observed velocities, respectively, and e are the data errors. The regularisation function R has the form

$$\begin{array}{ll}{\displaystyle }R& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathcal{A}}}\int \left({\mathit{\gamma }}_{\mathrm{s}}^{\mathrm{2}}{\left(\mathrm{\nabla }\left({\log }_{\mathrm{10}}\left(p\right)-{\log }_{\mathrm{10}}\left(\widehat{p}\right)\right)\right)}^{\mathrm{2}}\right.\\ {\displaystyle }& {\displaystyle }\left.+{\mathit{\gamma }}_{\mathrm{a}}^{\mathrm{2}}{\left({\log }_{\mathrm{10}}\left(p\right)-{\log }_{\mathrm{10}}\left(\widehat{p}\right)\right)}^{\mathrm{2}}\right)\mathrm{d}A\\ \text{(5)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathcal{A}}}\int \left({\mathit{\gamma }}_{\mathrm{s}}^{\mathrm{2}}{\left(\mathrm{\nabla }{\log }_{\mathrm{10}}(p/\widehat{p})\right)}^{\mathrm{2}}+{\mathit{\gamma }}_{\mathrm{a}}^{\mathrm{2}}{\left({\log }_{\mathrm{10}}(p/\widehat{p})\right)}^{\mathrm{2}}\right)\mathrm{d}A,\end{array}$$

where γ_a and γ_s  are regularisation parameters, and $\widehat{p}$ are the a priori values for model parameters. Inversions were done for a wide range of γ_s and γ_a and optimal values were determined from an L-curve analysis. In the results shown here, we use γ_a=1 and γ_s=10 000 m. However, our results are insensitive to the particular values chosen.

For γ_s=10 000 m, γ_a=1 and the sliding exponent m=3, the corresponding basal slipperiness C and the ice rate factor A distributions are shown in Figs. S3 and S4. The average difference between modelled and observed ice speed is 29 m per year with a median of 13 m per year and a root mean square error of 103 m per year. The measured and modelled velocity fields for the region of Institute Ice Stream are displayed in the upper panels of Fig. 1. They agree well in this area, as the residual histogram for this region shows in the lower-left panel, but also for the entire continent; see Fig. S5. As a consequence of our inverse methodology, modelled ice velocities are in close agreement with measurements.

From the modelled stresses obtained with our ice-flow model we calculate the buttressing parameter θ as defined in Sect. 3. We do this for each of the three different definitions for θ (see Eqs. 11, 12, and 13). We then calculate the analytical fluxes predicted by the flux formula, i.e. Eq. (6). Note that we refer to these fluxes as analytical fluxes, although their calculation involves the use of our numerical ice-flow model for estimating the buttressing number θ.

We also calculate modelled grounding-line fluxes from modelled ice velocities. Since our modelled velocities are in good agreement with observed velocities, these modelled grounding-flux estimates will be in an equally good agreement with fluxes estimated directly from observed velocities. The analytical and the modelled flux estimates are then compared and analysed.

Figure 1Observed (Gardner et al., 2018) and modelled (b) ice speed in the region of Institute Ice Stream. The inset displays the location of the plotted area in Antarctica. Grounding lines are shown as black lines and streamlines are displayed in blue. Panel (c) shows a normalised bivariate histogram of the velocity residuals which are the differences between modelled and observed velocities within this area, that is, $\mathrm{\Delta }u=u_{\mathrm{modelled}}-u_{\mathrm{observed}}$ and $\mathrm{\Delta }v=v_{\mathrm{modelled}}-v_{\mathrm{observed}}$, and u and v are the horizontal components of the surface velocity vector, respectively. Panel (d) shows an ice-speed profile along the central line of Institute Ice Stream that is indicated in green in (a).

Download

When calculating grounding-line fluxes we interpolate nodal quantities of the computational mesh onto the (calculated) grounding line. The grounding line does not, as such, enter the numerical calculations made by our numerical ice flow model. As described in Sect. 2, it is the flotation mask – evaluated at the integration points – that determines the impact of the basal drag term. However, in a post-processing step we determine the positions of the grounding lines from the flotation mask. Our approximation of the grounding line is a piecewise linear curve, with each linear segment representing the grounding line within a given computational element (see Figs. S1 and B1). We then interpolate nodal values onto the central point of each linear segment. The same procedure is employed when calculating both analytical and modelled fluxes.

In Schoof (2007a, b), an expression for the grounding-line flux (q) of marine ice sheets is derived. While the analysis is primarily focused on a flow-line configuration where ice-shelf buttressing plays no role, Schoof (2007a) also estimates how the flux might be affected by a reduction θ in axial stress at the grounding line due to ice-shelf buttressing. The resulting analytical flux expression is

where q is the ice flux across the grounding line, h is the ice thickness, ρ_i the ice density, ρ_w the density of ocean water and g the gravitational acceleration (please note that in the related Eq. 17 of Gudmundsson, 2013 for the flux q there is a typo in the exponent of the basal slipperiness C). For grounded ice, the tangential component of the basal traction (τ_b) is related to the basal velocity (v_b) through the Weertman-type sliding law

where C is the basal slipperiness, and m is the stress exponent, while deviatoric stresses τ_ij and strain rates ${\dot{\mathit{ϵ}}}_{ij}$ in ice flow are linked via Glen's flow law

with $\mathit{\tau }=\sqrt{{\mathit{\tau }}_{ij}{\mathit{\tau }}_{ij}/\mathrm{2}}$ the second invariant of the deviatoric stress tensor, exponent n (often set to 3) and rate factor A. Here τ_ij denote the components of the deviatoric stress tensor and ${\dot{\mathit{ϵ}}}_{ij}$ the components of the strain rate tensor.

Figure 2Buttressing ratio θ₁ along the grounding lines of Filchner–Ronne Ice Shelf (a) and Ross Ice Shelf (b). Insets indicate the ice shelves' locations in Antarctica. Regions where the grounding line is over-buttressed, that is, θ≤0, are displayed in black. Modelled speed is plotted in grey ranging up to 1500 m a⁻¹. Grounding line and ice front locations are indicated in black. IS denotes ice streams; IR denotes ice rises or rumples.

Download

As mentioned above, θ is a scalar quantity that describes the deviation in deviatoric axial stress at the grounding line from the unbuttressed situation. For an unbuttressed grounding line in one horizontal dimension (i.e. no variations in any quantities transverse to the flow direction) and assuming that the x axis of the coordinate system is aligned with the flow, we have τ_xx=τ_f where

which can be derived from the stress boundary condition at the calving front (see Appendix A). In the buttressed case, τ_xx is, however, no longer necessarily equal to τ_f, and θ is defined as

Here we have used the superscript 1HD to indicate that this definition of θ is only unambiguous in the one horizontal dimensional situation (1HD). In the more general two horizontal dimensional situation (2HD), where the flow direction is not necessarily aligned with the (horizontal) normal to the grounding line, several different definitions of θ are possible, and in the literature at least three different definitions of θ have been suggested. In the following we denote these by θ₁, θ₂, and θ₃, with

where n₁ is the normal to the grounding line, pointing horizontally outwards from the grounded ice into the ice shelf, and

and

where n₂ is the direction of ice flow at the grounding line and

is the (horizontal) deviatoric stress tensor, and

is the tensor of resistive stresses. In the 1HD unbuttressed case where n₁=n₂, ${\mathit{\tau }}_{xx}={\mathit{\rho }}_{\mathrm{i}}gh(\mathrm{1}-{\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{w}})/\mathrm{4}$, and ${\mathit{\tau }}_{yy}={\mathit{\tau }}_{xy}=\mathrm{0}$, all these three definitions of θ result in ${\mathit{\theta }}_{\mathrm{1}}={\mathit{\theta }}_{\mathrm{2}}={\mathit{\theta }}_{\mathrm{3}}=\mathrm{1}$. The first definition (i.e. θ₁) has, for example, been used by Gudmundsson (2013) to diagnose buttressing at the grounding line of an idealised set-up. The second definition has, for example, been used by Pollard and DeConto (2012), Thoma et al. (2014), and Pattyn (2017) as a flux condition, and the third one has been used by Fürst et al. (2016) to diagnose “flow buttressing” within Antarctic ice shelves. Note, however, that Pollard and DeConto (2012), Thoma et al. (2014, see Sect. 4.4), Fürst et al. (2016, see Supplementary Eq. 2), and Pattyn (2017; see Eq. 20) appear to use a different expression for τ_f, with ${\mathit{\tau }}_{\mathrm{f}}={\mathit{\rho }}_{\mathrm{i}}gh(\mathrm{1}-{\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{w}})/\mathrm{2}$, in which case $\mathit{\theta }=\mathrm{1}/\mathrm{2}$ in the unbuttressed case and θ in Eq. (6) must be replaced by 2θ.

Figure 3Buttressing ratio and differences in grounding-line flux for Institute Ice Stream draining into the Filchner–Ronne Ice Shelf (location shown in inset). (a) Buttressing values θ₁ are displayed along the grounding line and principle deviatoric stresses are shown with compression in red and extension in blue. The length of the vectors indicate the magnitude of each principle stress. (b) Differences between analytical and modelled fluxes and observed ice velocities ranging up to 500 m a⁻¹ (Gardner et al., 2018). Analytical fluxes are set to 0 where θ₁<0. Grounding-line positions are indicated in black.

Download

Figure 4Comparison of fluxes calculated with Úa (blue) and analytical fluxes (black) along the grounding lines of four major ice streams draining into the Filchner–Ronne Ice Shelf. Locations where the flux formula provides unphysical results are marked in grey. Plotted grounding-line segments are located as displayed in the inset with western margins indicated by a yellow dot.

Download

The definition of θ₁ is motivated by the form of the boundary condition at the calving front in the shallow ice-stream approximation (see Appendix A). For θ₁=1 the normal traction at the grounding line equals that of a calving front. In the general 2HD situation, this same interpretation does not hold for the definitions of θ₂ and θ₃. If θ₁>1 the ice shelf can be considered to be “pulling” the ice at the grounding line, while θ₁<1 implies that the ice shelf causes a reduction in normal traction at the grounding line; i.e. the ice shelf “holds the ice back”. Note that, for all these three different definitions, it is possible for θ to become negative. If, however, a negative θ value is inserted into Eq. (6), the resulting value for the flux q is a negative or even a complex number for most combinations of n and m – a clear indication that the analytical flux formula fails in such situations. Only the specific combinations of n and m, such that $nm/(m+\mathrm{1})=\mathrm{2}k$ for k∈ℕ (for instance the combination n=3 and m=2), “fix” the flux to a positive real number; however they introduce a non-substantiated dependency between the flow law and the sliding law. Furthermore, for these combinations and θ<0, enhanced buttressing inconsistently yields an increase in ice flux. Physically, θ₁<0 corresponds to a situation where the traction vector at the grounding line points in upstream direction. One possible situation giving rise to θ₁<0 would be τ_xx<0 while τ_yy=0, with x being the flow direction and the grounding line aligned with the y axis. In this case, the ice at the grounding line experiences compression in the along-flow direction and, hence, longitudinal strain rates are negative and ice velocities become smaller as the grounding line is approached from upstream direction. Another situation giving rise to θ₁<0 is that of equal transversal compression and vertical extension of the ice column at the grounding line; i.e. ${\mathit{\tau }}_{yy}=-{\mathit{\tau }}_{zz}<\mathrm{0}$ while τ_xx=0.

From the numerically modelled stress field we calculate the buttressing parameter θ₁ (given by Eq. 11) for all grounding lines of the Antarctic set-up described in Sect. 2.1. While here we focus on the buttressing parameter θ₁, our findings are independent of the exact definition of θ, the choice of the sliding law exponent m, the mesh and the details of the inverse methodology applied (see Appendix B).

We find that the grounding lines of the Filchner–Ronne and Ross ice shelves are, in general, highly buttressed with buttressing values significantly different from unity (see Fig. 2). Typically, θ₁≤0.4, and in many cases θ₁<0. Among the ice streams of these two biggest ice shelves of the AIS, the dormant Kamb Ice Stream is the least buttressed one, with θ₁≈0.4. Over all other ice streams θ values are even smaller. Negative θ values are also found over grounding-line segments located between active ice streams, for example along the grounding line running between the Rutford and Institute ice streams.

An example of an ice stream where θ₁<0 over most of its grounding line is the Institute Ice Stream (see Figs. 3 and 4). Inspection of the velocity field in the vicinity of the grounding line of that ice stream reveals that ice flow velocities decrease with distance as the grounding line is approached from upstream direction (see also Fig. 1d). Consequently, both along-flow strain rates and along-flow deviatoric stresses are negative (compressive). This general feature of ice flow around the grounding line of Institute Ice Stream implies that its grounding line is “over-buttressed” with the traction vector at the grounding line pointing in an inland direction. Hence, independently of our numerical simulation of the stress field, it is clear that for this ice stream a negative value for θ is obtained.

As discussed in Sect. 3 the analytical flux formula (Eq. 6) is clearly not applicable in situations where θ becomes negative. As θ is found to become negative over large sections of the grounding lines of many the ice streams of the two largest Antarctic ice shelves, i.e. Ross and Filchner–Ronne ice shelves, it follows that the formula cannot be used to calculate grounding-line ice fluxes over significant parts of the AIS.

We furthermore compare analytical and numerically modelled grounding-line ice fluxes in all regions where θ₁≥0, i.e. where the application of the analytical flux formula (Eq. 6) results in real-valued ice fluxes. In particular, we compare both the flux values pointwise along all grounding lines (Fig. 5) and the total cumulative fluxes over grounding lines of ice streams and ice shelves (Table 1). When comparing cumulative analytical fluxes, we are forced to assume values for those sections of grounding lines for which θ is negative (and q complex). There we assume q=0, which is equivalent to setting θ=0.

In general, we find significant differences between analytically calculated and numerically modelled flux values. Analytical fluxes are much lower than modelled in many locations of the Filchner–Ronne Ice Shelf, especially along the grounding lines of the Rutford, Institute and Moeller ice streams (Fig. 5). However, cumulative analytical fluxes over all grounding lines of the Filchner–Ronne Ice Shelf are about 30 % larger than modelled for θ₁, and this difference is considerably larger for θ₂ and θ₃ (Table 1). A similar disagreement between analytical and modelled fluxes is found for the Siple Coast ice streams such as Bindschadler and MacAyeal ice streams, and for Byrd Glacier (Fig. 5b). For Ross Ice Shelf the overall difference is only 5 %, but, given the fact that θ₁ is negative over significant sections of its grounding line (where we set the analytical flux values to zero), this agreement appears somewhat fortuitous.

For other ice shelves, cumulative fluxes are generally underestimated by the flux formula. Analytical fluxes for Pine Island Glacier and Thwaites Glacier, for example, deviate by −33 % and −52 % from the modelled fluxes, respectively. For George VI ice shelf, cumulative analytical fluxes are several times smaller than modelled ones (Table 1).

The analytical flux formula tends to strongly overestimate fluxes over grounding lines where ice flow is approximately tangential to the grounding line. The failure of the flux formula to correctly predict fluxes in such circumstances is not surprising, as the underlying assumptions of the formula are clearly not met in these situations. Nevertheless, this demonstrates the inherent conceptual difficulties in applying the formula to the Antarctic Ice Sheet.

Moreover, the analytical formula produces much higher spatial variability in fluxes than the numerically modelled ones. This can be clearly seen in Fig. 4, where analytical and modelled ice fluxes are plotted along the grounding lines of Rutford, Institute, Foundation and Recovery ice streams. Here, the grey background indicates sections of the respective grounding lines where the flux formula yields unphysical results. Variability in fluxes calculated with Úa occurs when ice flow is nearly aligned with the grounding line. We calculate fluxes within each triangular element using the normal of the piecewise-linear grounding-line curve, which may vary among individual line segments.

Table 1Ice flux integrated along the grounding lines of Antarctic ice shelves. Q_Úa denotes the modelled ice flux with Úa, Q₁ was derived from the analytical flux formula based on θ₁, Q₂ based on θ₂ and Q₃ based on θ₃, respectively. The last column shows the deviation of the analytical flux Q₁ from the modelled Q_Úa.

Download Print Version | Download XLSX

Figure 5Difference between the analytical and the modelled fluxes along the grounding lines of Filchner–Ronne Ice Shelf (a) and Ross Ice Shelf (b). Analytical fluxes are calculated based on θ₁ defined in Eq. (11). In locations where the formula yields unphysical results, fluxes are set to zero. Grey arrows show the modelled ice flow. IS denotes ice streams; IR denotes ice rises or rumples. Grounding line and ice front locations are indicated in black.

Download

We test the sensitivity of our analytical flux calculations to different degrees of regularisation (γ_s and γ_a), different values of the sliding law stress exponent (m) and relaxation of the ice geometry, for which our findings are summarised in Appendix B. Numerically modelled fluxes are, as expected, mostly independent of the value of the sliding law stress exponent m. This can be considered to be a consequence of the inversion procedure, which ensures that modelled velocity fields agree closely with measured data, independently of the value of m. On the other hand, analytically calculated flux values are highly sensitive to the value of m (see Fig. B3). For example, cumulative analytical fluxes for Filchner–Ronne Ice Shelf increase by about a factor of 5 as m is changed from 1 to 7, while numerically modelled fluxes change by less than 10 %. Numerically modelled fluxes are also insensitive to the exact degree of regularisation applied, whereas analytically calculated flux values change significantly (Fig. B5). The dependency of the analytically calculated fluxes on the amount of regularisation used in the numerical model is due to the impact regularisation has on modelled stresses and, therefore, on the value of θ.

We also compare analytical fluxes as calculated using the three different definitions (Eqs. 11, 12 and 13) for θ. While overall spatial variability of θ is similar for these three definitions, with all definitions giving rise to extended areas of negative θ values, the cumulative flux for the alternative definitions θ₂ and θ₃ are generally higher than for θ₁ (see also Fig. B4). Extended negative θ values and large flux deviations are consistently found after a short relaxation period (see Figs. S6 and S7).

The analytical grounding-line flux formula Eq. (6) was derived for a flow-line configuration (Schoof, 2007a), and there is no reason to doubt its validity in that particular case. When applied to a flow-line configuration, many current ice-flow models employing the shallow ice-stream approximation (SSA) with Weertman-type sliding law, have demonstrated excellent agreement between modelled and analytical grounding-line fluxes (Pattyn et al., 2012). Ice fluxes and grounding-line positions calculated with the ice flow model Úa also agree closely with those predicted by Eq. (6), where such an agreement is to be expected. The inclusion of the buttressing parameter θ was used by Schoof (2007a) to illustrate the potential impacts of ice-shelf buttressing on ice flux, provided its effects were sufficiently small to not invalidate the basic assumption of a flow-line setting too strongly. However, we find that most of the grounding lines of the AIS are highly buttressed, with θ significantly different from unity. It seems likely that at least part of the reason that the analytical flux formula fails relates to the high degree of buttressing that we find to be characteristic for most Antarctic ice streams.

When applied to the current geometry and the current flow field of the AIS, the flux formula predicts either unphysical or highly inaccurate flux values when compared to modelled ones. While we have done the comparison with numerically modelled fluxes, a comparison with observed fluxes – calculated from measured surface velocities, observed grounding-line positions, and measured ice thicknesses – would not alter our conclusions as, due to our inversion procedure, observed and modelled surface velocities are in good agreement.

The strongest indication that the analytical flux formula fails when applied to the Antarctic Ice Sheet is arguably the fact that it predicts non-real valued fluxes over significant parts of Antarctic grounding lines. This happens whenever θ becomes negative. Although for specific combinations of n and m (such as n=3 and m=2) the resulting exponents in the flux formula are even numbers – in which case the analytical fluxes are always real positive numbers – the flux values are still unphysical (see Sect. 3). As we point out above, even a cursory inspection of the velocity field of the AIS suffices to show that θ is negative for a number of grounding lines (e.g. the Institute Ice Stream grounding line). Hence, the occurrence of negative θ values is not simply a feature of our particular numerical approach, but a general aspect of the current ice-flow regime of the AIS.

As analytical ice fluxes are strongly dependent on ice thickness (h) at the grounding line, they depend somewhat on the specifications of the numerical model: the exact location of the grounding line is influenced by the mesh resolution used by the model. The resulting error is an example of a discretization error that becomes smaller as the mesh is refined. Other numerical models using a different computational mesh may locate the grounding line differently and hence calculate different ice flux values. We tested the dependency of our modelled ice fluxes to grid resolution by using several different meshes – an example of two such meshes is given in Fig. B1 – and found none of our main conclusions to be affected by differences in mesh resolution.

As measured by the buttressing parameter θ₁, almost all grounding lines of the AIS can be considered to be strongly buttressed with, in most cases, θ<0.4. Hence, theoretical concepts based on the assumption of none, or insignificant, ice-shelf buttressing may not apply to present-day Antarctica. One such theoretical prediction of considerable relevance for the possible future of the AIS relates to the stability of its grounding lines. In the absence of ice-shelf buttressing, grounding-line stability is predicted to be related to local bed slope (Schoof, 2007a, b, 2011; Thomas and Bentley, 1978; Weertman, 1974). However, in the presence of ice-shelf buttressing no such simple conclusions can be drawn (Goldberg et al., 2009; Gudmundsson, 2013; Gudmundsson et al., 2012; Pegler, 2016; Schoof et al., 2017). Possibly, rather than being dominated by local bed slope, the stability regime of the Antarctic Ice Sheet is to a leading order dependent on the properties of the ice shelves downstream of its grounding lines (e.g. geometry and structural integrity), as also supported by Pegler et al. (2013) and Haseloff and Sergienko (2018). Further work is needed to address the question of the stability of Antarctica's grounding lines.

The data and code that support the findings of this study are available from the corresponding author upon request.

The supplement related to this article is available online at: https://doi.org/10.5194/tc-12-3229-2018-supplement.

All authors designed the study. RR carried out the analysis based on an Antarctic set-up and using the ice-flow model Úa, both developed by GHG. RR prepared the manuscript with contributions from GHG and RW.

The authors declare that they have no conflict of interest.