Estimating ice discharge of the Antarctic Peninsula using different ice-thickness datasets

Abstract

The Antarctic Peninsula Ice Sheet (APIS) has become a significant contributor to sea-level rise over recent decades. Accurately estimating the ice discharge from the outlet glaciers of the APIS is crucial to quantify the mass balance of the Antarctic Peninsula. We here compute the ice discharge from the outlet glaciers of the APIS north of 70${^\circ }$S for the five most widely used ice-thickness reconstructions, using a common surface velocity field and a common set of flux gates, so the differences in ice discharge can be solely attributed to the differences in ice thickness at the flux gates. The total volumetric ice discharge for 2015–2017 ranges within 45–141 km³ a⁻¹, depending on the ice-thickness model, with a mean of 87 ± 44 km³ a⁻¹. The substantial differences between the ice-discharge results, and a multi-model normalized root-mean-squared deviation of 0.91 for the whole data set, reveal large differences and inconsistencies between the ice-thickness models, giving an indication of the large uncertainty in the current ice-discharge estimates for the APIS. This manifests a fundamental problem of the region: the scarcity of appropriate ice-thickness measurements and the difficulty of the current models to reconstruct the ice-thickness distribution in this complex region.

1. Introduction

Over the last decades, the Antarctic Peninsula Ice Sheet (APIS) has become a significant contributor to sea-level rise (SLR). Its rate of mass change has increased almost linearly from −6 ± 13 Gt a⁻¹ for 1997–2002 to −35 ± 17 Gt a⁻¹ for 2007-2012, stabilizing afterwards at −33 ± 16 Gt a⁻¹ for 2012–2017 (Shepherd and others, 2018). This substantially augmented mass loss since about 2000 is mostly due to increased ice discharge, associated to the acceleration of the glaciers feeding the ice shelves following ice-shelf disintegration events like those in 1995 (Larsen A) and 2002 (Larsen B), with a partial offset of the losses by increased snowfall (Fox-Kemper and others, 2021). Though the accelerated flow of tributary glaciers was ongoing well after the collapses (e.g. Helm and others, 2014), the rate of mass loss has decreased in the 20 years since the speed-up following the major ice-shelf disintegrations in 1995 and 2002 (e.g. Meredith and others, 2019, and references therein). Such ice-shelf collapses have been attributed to hydrofracturing of the ice shelf resulting from strong regional atmospheric warming (Scambos and others, 2009; Banwell and others, 2013). The particular recent history in this area becomes even more intriguing because of the alternating periods of warming and cooling. Specifically, the Antarctic Peninsula (AP) showed one of the strongest warming trends on Earth during the second half of the 20th century (Turner and othets, 2005), but this was followed by a cooling period between the end of the 20th century and the mid-2010s (Turner and others, 2016), mostly focused on the northern AP and the South Shetland Islands (Oliva and others, 2017), in turn followed by a return to warming conditions (Carrasco and others, 2021).

Glacier speed-up has also been observed in the southern part of the western coast of the AP and attributed to ice-shelf and tributary glaciers thinning, as well as marine-terminating glacier front retreat by increased calving, both induced by upwelling of warm and saline circumpolar depth water (Pritchard and others, 2012; Cook and others, 2016; Hogg and others, 2017; Friedl and others, 2018).

The projected contribution to SLR to the end of the 21st century from the APIS is, however, expected to be small (Edwards and others, 2021). Excluding the possible effects of tributary glacier acceleration due to ice-shelf collapse, the modeled response of the APIS to climate warming shows little detectable temperature dependence. The reason is that model projections predict a virtual balance between increased ice discharge and increased snowfall. For the APIS, this is even more marked under high-emission scenarios, in which snowfall accumulation increases more than mass loss, which is predominantly ocean-induced (Edwards and others, 2021). If ice-shelf collapse scenarios are considered, these dominate the APIS response. Even so, the effect is modest, with a median contribution of 1 cm of sea-level equivalent (SLE) for 2015–2100 (Edwards and others, 2021). The effect is small because surface meltwater is not projected to be sufficient to trigger ice-shelf collapses until the second half of the century (Seroussi and others, 2020). Curiously, the situation for the peripheral glaciers (mostly located on the surrounding islands) is opposite, with limited contribution to SLR at present (Hugonnet and others, 2021; Shahateet and others, 2021) but a much larger projected contribution to the end of the 21st century (Edwards and others, 2021).

Accurately estimating the ice discharge from either tributary glaciers feeding the ice shelves (or the embayments resulting from their collapse) or from tidewater glaciers is therefore crucial to quantify the mass losses from the APIS. Ice discharge is computed as the flow of glacier ice through flux gates defined at convenient locations (either close to the grounding zone of an ice shelf tributary glacier, or of a glacier with floating tongue, or close to the calving front of a tidewater glacier). Such flux is given by the integral, over the flux-gate section, of the velocity field times the cross-sectional area (Cogley and others, 2011). While the velocity field is fairly well known for glaciers on the AP (e.g. Seehaus and others, 2018, and references therein), knowing the cross-section at the flux gate requires the knowledge of the ice-thickness. In the case of the AP, such ice thickness has been only measured at particular locations covered by airborne radar campaigns, mostly from NASA's operation IceBridge (MacGregor and others, 2021). Where no direct ice-thickness measurements are available, the ice thickness has to be inferred, using inverse methods, from other available data such as surface topography, surface velocity, surface elevation change rates or surface mass balance.

In the case of the AP, five main ice-thickness reconstructions are currently available. They are based on spatial interpolation of observations (Fretwell and others, 2013), in cases additionally supported by velocity observations (Morlighem and others, 2020), or are mass-conservation-based inverse modeling approaches calibrated with observations (Huss and Farinotti, 2014), or models based on a perfect plasticity assumption (Carrivick and others, 2018), or deep-learning-based approaches (Leong and Horgan, 2020). We note that the study by Morlighem and others (2020) uses, in general, a mass-conservation-based approach. However, in cases where the flow velocity is small, it uses instead an anisotropic interpolation scheme based on streamline diffusion. The latter has been the approach mostly used by Morlighem and others (2020) for reconstructing the ice thickness of the AP north of 70°S.

The aim of this paper is to compute the ice discharge from the outlet glaciers of the APIS north of 70°S (Fig. 1) for the five above-mentioned ice-thickness reconstructions, using a common surface velocity field and set of flux gates, so the differences in ice discharge can solely be attributed to the differences in ice thickness at the flux gates available for such reconstructions. In this way, the differences in computed ice discharge can be considered as an indicator of the current uncertainties in the ice discharge estimate and hence of the contribution of the APIS to sea-level rise from this source.

Figure 1. Study site, surface velocity field and flux gates on the Antarctic Peninsula. The velocity field is the mean value for the period 2015–2017 of the velocity values provided by the ESA Antarctic Ice Sheet Climate Change Initiative project.

2. Data and methods

To evaluate the differences implied by the use of different ice-thickness models in the calculation of the ice discharge from the APIS, we kept all other intervening variables (surface velocity field and flux gate definition) fixed and only the employed ice-thickness map varies from case to case.

Ice discharge is calculated as the flux of ice through a flux gate close to the grounding line or calving front (Fig. 2). Therefore, it requires ice velocity field and ice-thickness data at the flux gate (Cuffey and Paterson, 2010):

(1)$$\phi = \int_S \rho {\bf v} .{\bf dS} = \sum_i{\rho L_i H_i f v_i \cos{\alpha_i}}$$

Figure 2. Scheme of the ice discharge calculation using flux gates. Ice discharge is calculated as the ice mass that crosses the flux gate (with area A) per time unit. The surface velocity field across the flux gate is represented, as well as a vertical profile of it at the central line.

The flux calculation is discretized by dividing the flux gate into smaller evenly spaced bins, so the flux is expressed as a finite sum, as shown in Eqn. (1) (Sánchez-Gámez and Navarro, 2018). ρ is ice density, L _i is the width of the individual segments (we used L _i = 200 m to make it comparable to the resolution of the velocity grid), H _i is the ice-thickness of the segment i, f is the ratio of the surface to depth-averaged velocities, v _i represents the speed of the individual segment i, and α is the angle between the velocity vector and the vector normal to the flux gate surface. For the factor f we took a value of 1.0, which means that all motion is due to sliding rather than internal deformation. This is a common assumption for ice-sheet outlet glaciers next to the grounding line (Gardner and others, 2018; Mankoff and others, 2020). To facilitate the comparison with the results of other authors without being conditioned by the choice of a particular value of ice density, we will henceforth refer to volumetric ice discharge φ (same as Eqn. (1) but excluding ρ).

We used velocity data available at https://cryoportal.enveo.at/data/, processed and provided by ENVEO (Environmental Earth Observation IT GmbH), as part of the ESA Antarctic Ice Sheet Climate Change Initiative project (Fig. 1). The data consist of monthly Antarctic ice velocity mosaics derived from Sentinel-1 Synthetic Aperture Radar data using feature-tracking techniques for the period 01/Jan/2014-31/Dec/2021. The data are provided in stereographic projection (EPSG:3031) with 200 m resolution. The horizontal velocity components are provided in true metres per day, towards the x and y directions of the grid. From these data we calculated a mean velocity field for the period 01/Jan/2015-31/Dec/2017 for the Antarctic Peninsula region, north of 70°S, converting the grid to true geographic directions (North and East).

Using the above velocity field, together with glacier catchment delineations (Silva and others, 2020), optical satellite images and available grounding line information (from the Antarctic Digital Database), we manually defined 325 flux gates across our study area. Following Seehaus and others (2018), we defined the flux gates as across-glacier profiles for each glacier basin. For marine-terminating glaciers, such profiles were defined close to the terminus while taking into account the maximum frontal retreat state, based on the multi-temporal inventory by Silva and others (2020). For glaciers with known grounding line position, in particular for those in the Larsen-C embayment, the flux gates were defined as across-glacier profiles slightly up-glacier from the grounding line. Whenever possible, a straight line, perpendicular to the main ice flow direction, was defined. However, due to the complex geometries of various glacier basins (e.g. coalescence of multiple glacier branches towards the terminus), kinked lines were used as flux gates to better account for the complex ice flow pattern. Figure 3 illustrates the position of the flux gates, together with the radar flight lines supplying the ice-thickness observations. The flux gate positions are provided as supplementary material to this paper. We note that not all ice-thickness models use ice-thickness data as input, and that models using them as input (often for model calibration) employ such data in different ways (this will be dealt with in the Results and Discussion section). Moreover, not all the models using ice-thickness data as input used all the radar data shown, as such models were developed at different dates, implying slightly different sets of radar data available.

Figure 3. Map of the Antarctic Peninsula north of 70${^\circ }$S showing: (a) the location of three sectors (red rectangles) shown in the other panels. (b–d) Flux gates used in our study (orange lines) and Operation Ice Bridge radar flight lines (black lines) used as input.

As pointed out in previous sections, our aim is to assess the differences in ice discharge calculation for the APIS arising from the use of different ice-thickness models. In particular, we used the following five models with the given spatial resolutions: Carrivick and others (2018) and Huss and Farinotti (2014) with a resolution of 100 m, Bedmap2 (Fretwell and others, 2013) with 1 km resolution, DeepBedMap (Leong and Horgan, 2020) with 250 m resolution and Bedmachine v2 (Morlighem and others, 2020) with 500 m resolution. In the case of DeepBedMap, a digital elevation model (DEM) of the bed topography is originally provided. To convert into ice thickness, we calculated the difference between DeepBedMap and the REMA DEM (Howat and others, 2019), which was used as input for DeepBadMap, and resampled at the same resolution as DeepBedMap (250 m). In cases where the resulting ice thickness presented negative values, we masked them out.

To have better spatial consistency, we resampled the ice-thickness datasets to the same resolution as the velocity field (200 m), using bilinear interpolation.

In the following, we briefly describe the various methodologies used by Carrivick and others (2018), Huss and Farinotti (2014), Fretwell and others (2013), Leong and Horgan (2020) and Morlighem and others (2020) to generate the ice-thickness distribution for the AP (north of 70°S) from the available data. These models will be the focus of our analysis.

2.1 Ice thickness datasets

2.1.1 Carrivick

Carrivick and others (2018) implemented an ice-thickness model of the Antarctic Peninsula, between 60°S and 74°S, which does not require bed information. They included all of the mainland ice sheet and all glaciers on the surrounding islands at a resolution of 100 m. We note that we will use only their data north of 70°S and exclude the islands surrounding the peninsula. Carrivick and others (2018) used the VOLTA model described in James and Carrivick (2016), which is a perfect-plasticity approach to calculate ice thickness, bed elevation and glacier volume. The VOLTA model solves, along glacier centerlines, the equation

(2)$$h = {\tau_{\rm b}\over f \rho g \tan \alpha}$$

where h is the desired ice thickness, τ _b is the basal shear stress, f is a shape factor that accounts for drag at the glacier valley sides, ρ is the ice density, g is the acceleration of gravity and α is the surface slope. Following Driedger and Kennard (1986), Carrivick and others (2018) calculated τ _b by means of the empirical equation

(3)$$\tau_{\rm b} = 2.7 \cdot 10^4 \sum_{i = 1}^n\left({A_i\over \cos \alpha_i}\right)^{0.106}$$

where τ _b is given in Pa, A _i is the elevation band area given in m² and α _i is the surface slope. In addition to a value for f, Eqns (2) and (3) only require glacier outlines, surface elevation and the definition of the centerlines. Carrivick and others (2018) used the DEM provided by Cook and others (2012). The glacier outlines were a combination of Davies and others (2012), Bliss and others (2013) and Cook and others (2014). The centerlines were manually defined using the velocity field of Rignot and others (2011) as a reference. Once ice thickness along centerlines was calculated, distributed ice thickness was interpolated across each glacier, at 100 m resolution, using an interpolation routine by Hutchinson (1989), designed to generate concave-shaped beds that mimic the typical parabolic shape of idealized glacier beds. More detailed information can be found in Carrivick and others (2018) and James and Carrivick (2016).

2.1.2 Huss and Farinotti (2014)

Huss and Farinotti (2014) provided an ice-thickness map of the AP region (north of 70°S) at a resolution of 100 m. They inferred the local ice thickness assuming the shallow ice approximation (SIA), using the DEM from Cook and others (2012), surface mass balance from RACMO, ice-thickness measurements mostly from Operation IceBridge and velocity fields from Rignot and others (2011). Assuming conservation of mass, they calculated an apparent mass balance (Eqn (4)) using the information on the surface mass balance of downscaled RACMO output, together with surface elevation change.

(4)$$\tilde{b} = b \cdot f_{{\rm RACMO}}-\rho \cdot {\partial h\over \partial t}$$

In Eqn (4), ρ is the ice density and f _RACMO is a calibration parameter to account for the potential bias in the RACMO model. They parameterized the surface elevation change (∂h/∂t) from 0 m a⁻¹ at 700 m a.s.l. to −10 m a⁻¹ at the calving front for glaciers flowing into Larsen A and B. For the rest of the glaciers, they defined −1 m a⁻¹ on the ice front. By integrating $\tilde {b}$ over the entire glacier, they obtained the ice volume flux q _tot for each glacier.

Huss and Farinotti (2014) then applied Glen's flow law to calculate the ice thickness h _i for every elevation band i:

(5)$$h_{i} = \sqrt{[ ( n + 2) ] {{q_{{\rm d}, i}\over 2 A_{{\rm f}}} \cdot {n + 2\over \left(F_{{\rm s}, i} \rho g \sin \overline{\alpha_{i}}\right)^{n}}}}$$

where q _d,i is the deformation component of the total ice flux normalized by the glacier width (w), n = 3, g is the acceleration of gravity, $\overline {\alpha _{i}}$ is the average slope of the surface of each elevation and F _s,i is the valley shape factor F _s,i = w _i/(2 h _i + w _i). F _s is solved iteratively because it is a function of h. The flow rate factor (A _f) is another calibration parameter. q _d is the deformation component of the total ice flux (q _tot) for each elevation band i. It is calculated as

(6)$$q_{{\rm d}, i} = q_{{\rm tot}, i}-q_{{\rm b}, i} = q_{{\rm tot}, i} \cdot\left(1-{\,f_{{\rm sl}, i}\over ( 1-r) \cdot f_{{\rm sl}, i} + r}\right)$$

where f _sl is the fraction of basal sliding velocity to surface velocity, q _b is the ice flux due to basal sliding and r = (n + 1)/(n + 2) = 0.8 is the relationship between the average speed of the deformational flow and the surface velocity of an ice column.

2.1.3 Bedmap2

Bedmap2 (Fretwell and others, 2013) provides an ice thickness and bed elevation distribution of the entire Antarctica with a resolution of 1 km. It is a refinement (they used 10 times more ice-thickness measurements) of Lythe and Vaughan (2001), which provided a resolution of 5 km. Fretwell and others (2013) used different methods depending on the availability of ice-thickness measurements. For the Antarctic Peninsula, three methods were applied: direct ice-thickness measurements, Griggs and Bamber (2011) shelf thickness reconstruction and synthetic ice-thickness data.

For the Bedmap2 reconstruction, around 24.8 million points of ice-thickness measurements were gathered. Many are freely available for download (http://nsidc.org/data/), while others are presented in summary publications. From these points, after filtering, only about 140 000 were used in the reconstruction. The ice thickness in zones between measurements was linearly interpolated.

In regions covered by ice shelves, Bedmap2 uses ice thickness derived from altimetry satellites, assuming hydrostatic equilibrium. Areas where the ice was grounded were excluded and, to avoid introducing bias from failure of the assumption of hydrostatic equilibrium, data within 5  km of the grounding line were also excluded (Fretwell and others, 2013).

Synthetic ice thickness data were generated to prevent rock outcrops from overly skewing the ice thickness distribution in mountainous areas with few direct measurements. They also used this approach to define major glaciers passing through mountain ranges for which ice-thickness measurements were too sparse to ensure the existence of such glaciers in the results. In the first case, in areas within 10 km of rock outcrops and more than 10 km away from ice-thickness data, they applied a thin-ice model as done by Lythe and Vaughan (2001). This model is based on the assumption of a general correlation between ice thickness (in meters) and the distance (in meters) from rock outcrops (H = 223.98ln(d) − 1108.4). In the second case, the nearest upstream and downstream data points were linearly interpolated to generate the synthetic data.

2.1.4 DeepBedMap

DeepBedMap (Leong and Horgan, 2020) is a super-resolution approach that involves processing a lower resolution raster image (in this case Bedmap2) into a higher resolution one. More specifically, they used an Enhanced Super-Resolution Generative Adversarial Network (ESRGAN) method. This method consists of two neural networks, one working as a generator and the other as a discriminator. The discriminator checks the perceptual quality (the extent to which an image looks like a valid natural image) of the product of the generator, which is quantified mathematically by the discriminator taking into account high-level features of an image like contrast, texture, etc. They also used as input the REMA DEM by Howat and others (2019), the ice velocity field from Mouginot and others (2019) and the snow accumulation from Arthern and others (2006).

The method is applied through a cost function

(7)$$L_{\rm G} = \eta L_1 + \lambda L_{\rm G}^{Ra} + \theta L_{\rm T} + \zeta L_{\rm S}$$

where L ₁ is the content loss, which takes into account the difference between the result of the method and the original ice-thickness measurements. $L_{\rm G}^{Ra}$ is the adversarial loss, where the discriminant neural network evaluates the perceptual quality of the generator. L _T is the topographic loss that accounts for the difference between the elevation values of the resulting map with respect to the original low-resolution image. L _S is the structural loss that analyzes the luminance, contrast and structural information between the predicted and the ice-thickness measurements. Finally, η, λ, θ and ζ are weights for the loss terms mentioned above. More detailed information is available in Leong and Horgan (2020).

2.1.5 Bedmachine Antarctica v2

Although the Bedmachine project (Morlighem and others, 2020) is designed to use the mass-conservation assumption to reconstruct the ice thickness in those regions of Antarctica with surface velocities greater than 50 m  a⁻¹, in the Antarctic Peninsula (even in regions with ice velocity greater than 50 m a⁻¹) it is almost exclusively based on the so-called ice-thickness diffusion equation

(8)$$\left\{\matrix{ \nabla \cdot {\bf D} \nabla H = 0 \hfill & {\rm in }\; \Omega \hfill \cr H = H_{{\rm obs}} \hfill & {\rm on }\; T \hfill }\right.$$

where D is the anisotropic diffusion tensor defined as:

(9)$${\bf D} = {\bf v}_{{\rm s}} \otimes {\bf v}_{{\rm s}} + \tau \;{\bf I}$$

where Ω is the domain of the model, T are the flight tracks where observations are available, v_s is the surface velocity, ⊗ denotes the tensor product, and $\tau = 30^\circ$ m² y⁻², if the ice velocity is greater than 30 m a⁻¹. In the regions where the velocity is undefined or less than 30 m a⁻¹, D becomes a scalar (i.e. isotropic). Equations (8) and (9) are not based on physics, but are a way to interpolate ice thickness between measurements anisotropically (Morlighem and others, 2020). The data are available through Morlighem (2020), and from now on we will refer to this data set as Bedmachine v2.

3. Results and discussion

With the aim of assessing the consistency among the various ice-thickness models, we compare in this section their ice discharge results and evaluate various aspects of the individual models.

Table 1 shows, for each model, the total ice discharge across the whole set of flux gates. Three sets of discharge values become apparent. The largest values correspond to Carrivick and others and Huss and Farinotti models, which provide discharge values nearly twice those of DeepBedMap and Bedmap2 (whose ice discharge estimates differ from each other by just 6%), and three times as large as that of Bedmachine v2. DeepBedMap and Bedmap2, in turn, produce discharge estimates 32 and 40% higher, respectively, than that of Bedmachine v2. One could think that the largest value for Carrivick and others model could be related to the fact that, according to Eqn. (2), ice thickness, h, will tend to infinity as surface slope tends to zero, meaning that h may be overestimated in regions of flatter ice surface. A similar situation is expected, though to a lesser extent, for Huss and Farinotti model (H&F) through Eqn. (5). Note, however, that the flat regions at higher elevations (in particular near the ice divides, where the slopes tend to zero) are of no concern to us, as the ice-discharge calculation is limited to flux gates near the marine termini of tidewater glaciers or the grounding line of ice shelves. In the latter case, the surface slope up-glacier from the grounding line is nonzero. Only close to the terminus of tidewater glacier fronts the surface slope tends to zero. Even so, both Carrivick and others and H&F models set a minimum slope threshold of 1.5${^\circ }$ to avoid this problem. Moreover, H&F model also sets a maximum ice-thickness value at grounding lines in terms of flotation conditions.

Table 1. Total volumetric and mass ice discharge 2015–17 calculated for the individual models

We used an ice density of 900 kg m⁻³ for the volume to mass conversion. The uncertainty in the mean ice discharge is the standard deviation.

With the exception of the Carrivick and others model, which does not use bed data as input, and DeepBedMap, which uses data from other parts of Antarctica for calibration, the three remaining models use basically the same set of ice-thickness data, and thus provide similar discharge results for the flux gates for which thickness observations are available. Although Bedmachine v2 used some additional airborne radar ice-thickness observations in comparison with Bedmap2 and H&F (see e.g. Figure S3 of Morlighem and others, 2020), those for the AP north of 70${^\circ }$S were very few, if we exclude radar flights over the ice shelves. Therefore, the differences among flux calculations for the three latter models mostly come from those flux gates where no ice-thickness observations are available (Fig. 3), so the various models use their respective methods to estimate the ice-thickness. The ice-thickness observations on the northernmost AP are very scarce (Fig. 3, sector 1), and are virtually absent on the Trinity Peninsula (northernmost part of Fig. 3, sector 1). The flight lines on the central and southern sectors (Fig. 3, sectors 2 and 3, respectively) are more abundant, but those on the western coast are generally more distant from the flux gates compared with those on the eastern coast.

Given the striking differences among the discharge estimates, it is worth comparing the results in Table 1 with those of the recent estimate by Rignot and others (2019). The latter authors calculated the mass balance of Antarctica, with detail of regions and subregions, using the input/output method, in which the main output component is ice discharge. For the AP region north of 70°S, they estimated an average ice discharge of 164.33 Gt a⁻¹ over 2015–2017 (their period overlapping with that of our velocity field). To compare with our results in Table 1, we use an ice density of 900 kg m⁻³ to convert from volume to mass ice discharge. These results highlight the proximity between our results using the ice-thickness models of Carrivick and others (2018) and Huss and Farinotti (2014) with those of Rignot and others (2019). The key point here is considering which ice-thickness dataset was used by Rignot and others (2019). In their study, if the available ice-thickness data were not of sufficient quality, they assumed a 1979 ice flux in balance with the average SMB for 1979–2008 and scaled the results based on changes in ice velocity, which is fundamentally a mass-conservation approach, and hence its closest agreement with our result for Huss and Farinotti (2014). In the case of the AP, this resulted in constraining fluxes with ice thickness and time-variable velocity for 78$\%$ of the Antarctic Peninsula. For the rest, they used both Bedmap2 and Bedmachine Antarctica v1.

In our study, the flux gate with the largest ice discharge is that of Seller Glacier catchment (which includes Fleming, Seller, Airy and Rotz glaciers), located on the south-western coast of our study region, with a mean of 12.76 ± 4.62 Gt a⁻¹. The Seller Glacier catchment terminates in Wordie Bay, where a series of ice-shelf disintegration events resulted in an acceleration of the mentioned glaciers (Friedl and others, 2018). Rignot and others (2019) calculated a similar value of 14.44 Gt a⁻¹ (for 2015–2017 in their case). The calculated ice discharge by the models considered here shows considerable spread, differing from Rignot's estimate by amounts between 3% (for H&F model) and 35% (for Carrivick's model).

To analyze how the whole set of models performed, on average, for each individual glacier, we calculated a multi-model normalized root-mean-square deviation (NRMSD, as shortened acronym) between the model results for each individual glacier:

(10)$$NRMSD_i = {1\over \overline{\varphi_{i}}} \sqrt{{\sum_{\,j = 1}^J( \overline{\varphi_{i}}-\varphi_{i, j}) ^2\over J}}$$

Here, φ_i,j is the flux calculated using model j for the flux gate i. $\overline {\varphi _{i}}$ is the average flux calculated using the various models (where j = 1,  2,  3,  4,  5, with J = 5) for an individual flux gate i. The results for each of the individual glaciers are shown in Figure 4, which illustrates that, although the differences between the models are large, they are not equally spatially distributed. We can observe that the various models perform similarly for the glaciers discharging along the eastern coast of the AP, particularly for the glaciers terminating in ice shelves, while on the western coast, there are much larger discrepancies between the results from the different models. We attribute this to the fact that the flux gates in the western coast are more distant from the radar flight lines compared to those in the eastern coast (Fig. 3).

Figure 4. Spatial distribution of the normalized multi-model root-mean square deviations for the ice discharge of the individual flux gates in the Antarctic Peninsula (see Eqn. (10)). Letters from a to f identify glaciers with profile plots in 6.

We also calculated the averaged value $NRMSD = \overline {NRMSD_i}$. The result for the whole set of flux gates and the five models was NRMSD = 0.91. This means that the average NRMSD is 91% of the mean ice discharge of all flux gates in the AP, representing a high average variability in the ice discharge of individual glaciers. Referring back to the differences between the western and eastern coasts, if we only consider the glaciers that terminate in ice shelves, the averaged NRMSD of the whole set is 0.50, while if we only consider the ocean- and land-terminating glaciers, it is 1.03. This difference in the NRMSD values indicates larger discrepancies between models in the estimation of the ice discharge of the ocean- and land-terminating glaciers compared with the ice-shelf-terminating ones (as visually shown in Fig. 4).

To analyze the difference in performance between pairs of models, we introduced

(11)$${\Delta\varphi_{\,jk}} = \overline{\left({\vert \varphi_{i, j}-\varphi_{i, k}\vert \over \overline{\left\lbrace\varphi_{i, 1},\; \varphi_{i, 2},\; \varphi_{i, 3},\; \varphi_{i, 4},\; \varphi_{i, 5}\right\rbrace}}\right)}$$

where the fraction in the right-hand side is the difference in the absolute value of ice discharge between the models j and k for flux gate i normalized by the mean of all models (the term in the denominator, where φ_i,1, φ_i,2, φ_i,3, φ_i,4 and φ_i,5 are the ice discharge calculations for flux gate i using the five different ice-thickness models (1, 2, 3, 4 and 5) and the averaging is understood to be applied over the five quantities enclosed by the braces). The global averaging bar is understood to be applied over the whole set of flux gates (represented by subscript i), so Δφ_jk is the normalized mean difference between models j and k for the whole set of flux gates.

The differences between pairs of ice-discharge calculations (using Eqn. (11)) are summarized in Table 2. It shows that there are large differences between the models. Huss and Farinotti (2014) exhibits the three largest differences, with values of 1.85, 1.73 and 1.62 with respect to Bedmachine v2, DeepBedMap and Bedmap2, respectively. The result of ice discharge using Huss and Farinotti (2014) is in better agreement with the one of Carrivick, with a mean-normalized difference of 0.94. We attribute the greatest differences between models shown by Huss and Farinotti (2014) model to the fact that it uses a mass-conservation approach, which is more physically based and clearly differs from those used by Bedmachine v2, DeepBedMap and Bedmap2. More precisely, as noted by Carrivick and others (2018), Huss and Farinotti (2014) used an ice flow model driven by mass-balance parameters; it is technically not a mass conservation approach due to the correction parameter f _RACMO in Eqn. ((4)). Carrivick and others model is also physically based, but it uses the perfect-plasticity approach instead of SIA (in H&F) and is not a mass-conservation approach. The results using Carrivick and others model present intermediate mean-normalized differences from the other models, with values of 1.30, 1.27 and 1.13 with respect to Bedmachine v2, DeepBedMap and Bedmap2, respectively.

Table 2. Mean of the normalized difference between pairs of models for each flux gate (see Eqn. (11))

The closest agreement is between Bedmap2 and Bedmachine v2, with a value of 0.37. We remind that this value represents the difference relative to the mean. Therefore, the mean difference between Bedmap2 and Bedmachine v2 is 0.37 times the mean of all models.

Although the total ice discharge from the APIS north of 70°S calculated using DeepBedMap agrees well with that obtained using Bedmap2 (see Table 1), comparing the values at the flux-gate level manifests a larger difference (Δφ₂₃ = 0.76). More importantly, focusing now on DeepBedMap, 44% of its ice-thickness data at the AP flux gates are meaningless negative values (see Figs 5e, f). Also, while Leong and Horgan (2020) highlight the sensitivity of DeepBedMap to the training data set, they do not have a training data set for the Antarctic Peninsula, which compromises their results in this region. Furthermore, it is possible to identify artifacts through their ice-thickness map.

Figure 5. Histogram of the distribution of ice-thickness points extracted along the flux gates for each model.

Finally, Fretwell and others (2013) and Morlighem and others (2020) had similar general results, as was expected due to the use of similar approaches (interpolation, although of a different type).

Another important parameter to assess the quality of the results is the number of flux gates producing zero ice-discharge results, as this is an indication of unrealistic thickness at the flux gates. Using Carrivick and others (2018) and Huss and Farinotti (2014) ice-thickness maps, we obtained non-zero ice discharge values for every single flux gate, whereas Fretwell and others (2013), Leong and Horgan (2020) and Morlighem and others (2020) produced, respectively, 88, 86 and 1 zero ice-discharge results at individual flux gates, out of 325 flux gates.

Once again, Figure 5 highlights the similarity between the Bedmap2 and Bedmachine v2 ice-thickness distributions (along the flux gates). As noted earlier, this is expected since they use similar approaches for the Antarctic Peninsula. However, while Bedmachine v2 includes a total of 220 zero ice-thickness points along the flux gates, Bedmap2 includes 1851 out of 7630 (see Table 3). This difference, however, is not as striking when we note that 1810 measurement points along the flux gates in Bedmachine v2 show ice-thickness values of less than 1 m, thus having little impact on the total ice-discharge calculation.

Table 3. Total number of ice thickness points along flux gates with negative, zero and positive values, categorized into values with point velocity greater than 0.5 m d⁻¹, between 0.5 and 0.1 m d⁻¹ and lower than 0.1  m d⁻¹

We also analyzed the individual points of the flux gate discretizations, counting the number of negative, zero and positive ice thicknesses and comparing them with the ice surface velocity (Table 3). As mentioned above, 43.6% of the DeepBedMap points show negative ice thickness values. These meaningless negative values were masked and considered as zero values in our calculation of ice discharge (Fig. 5f shows the ice-thickness distribution for DeepBedMap before masking negative numbers). Also, 24.3% of the Bedmap2 thickness points along the flux gates are zero. In these two cases (DeepBedMap and Bedmap2), we expect a high underestimation of the ice discharge, since the calculations of ice discharge for fast-flowing ice (>0.5 m d⁻¹) produce zero discharge due to the zero ice thickness. In the case of Carrivick and others (2018), Huss and Farinotti (2014) and Morlighem and others (2020), 99.25, 99.5 and 97.1% are points along the flux gates with meaningful positive ice-thickness values.

A fundamental problem of Carrivick and others (2018) model, pointed out by their own authors and also common to Huss and Farinotti (2014) model, is that neither the perfect-plasticity approach (Carrivick) nor the SIA (H&F) is expected to be appropriate for either calving tidewater glaciers or ice shelf tributary glaciers, because they have significant longitudinal stresses resulting in rapidly expanding flow near their terminus (tidewater glaciers) or as they approach the ice shelves. While the mass-conservation approach of Bedmachine does its best when modeling ice-thickness in regions of fast flow (>50 m  a⁻¹), in the case of Bedmachine Antarctica v2 it has rarely been applied to glaciers north of 70${^\circ }$S; instead, the streamline-diffusion interpolation scheme has been applied in most of the region, no matter whether the zones have fast flow or slow flow. This is likely due to the fact, pointed out by Morlighem and others (2020), that the precision of the mass-conservation product is affected by the spacing between ice-thickness measurements, which are used to constrain the calculation. A fundamental problem of DeepBedMap is that it does not use ice-thickness data from the AP to train its neural network algorithm. This could be the reason behind the large number of meaningless negative ice-thickness values produced in this region. Finally, a major limitation of Bedmap2 in the AP region is its coarse nominal spatial resolution of 1 km, which is actually much coarser due to the limited data used for interpolation (Carrivick and others, 2018). This makes it inappropriate for most of our study region, characterized by a complex and rapidly spatially-varying bed geometry. This could be the reason for the large number of zero ice-thickness points generated by this model at the flux gates, even in fast-flowing areas.

Taking into account that most of the ice-thickness reconstruction models used in this study have proved to be more consistent and in more agreement with measurements when applied at other locations (such as the Greenland and central Antarctic ice sheets, and mountainous regions), their deficiencies and inconsistencies in ice-thickness results shown when applied to the AP north of 70 °S, and their associated ice discharge estimates, manifest a fundamental problem of our region under study: the scarcity of appropriate ice-thickness data. The mountainous relief, with many outlet glaciers confined within steep ridges, makes airborne radar data collection extremely difficult. Across-flow profiles close to the calving fronts of tidewater glaciers or next to the grounding line of ice-shelves fed by many of the outlets, are critical for accurate ice-thickness estimates. However, the surrounding steep slopes make the data collection prone to failure. When along-flow profiles are involved, unwanted reflections from the glacier valley walls obscure the radar records and make their interpretation difficult or impossible. The highly crevassed glacier fronts, with their inherent scattering of radar energy, further complicate things. For these reasons, even if the Icebridge radar profiles are in many parts of the AP relatively abundant, when one looks at the actual data recorded a huge amount of data gaps or inconsistent data arises. In regions such as the northernmost AP, where a certain amount of temperate ice is expected, Icebridge data have been collected using radars designed for cold ice and the bed reflection is often simply not visible because of the strong scattering of the radar energy by the temperate ice, preventing data interpretation.

4. Conclusions and outlook

From the above discussion we can conclude that, when calculating the ice discharge of the Antarctic Peninsula, the final result is strongly influenced by the choice of the ice-thickness model. Some inconsistencies were found, such as the huge number of meaningless negative thicknesses at the flux gates generated by Leong and Horgan (2020), its lack of a training set for the AP and the presence of artifacts in the ice-thickness map. Also, Fretwell and others (2013) generated a large number of zero ice-thickness points at the flux gates, even in fast-flowing areas. Although Bedmachine v2 presented few zero ice-thickness points, 24% of them were less than 1 m, presumably implying an underestimation of the ice discharge.

Much of the poor performance of the ice-discharge calculations using all ice-thickness models analyzed, which have shown to be efficient when applied elsewhere, can be attributed to the complex geometry of the region, together with the lack of suitable ice-thickness data in large sectors of the AP. This scarcity has a significant impact on the models due to the insufficiency of constraints in the ice-thickness data, either used as direct input to the models or as calibration and validation data. In particular, DeepBedMap does not use for training of its neural network ice-thickness data from the AP, and Carrivick and others model only uses the ice-thickness data for validation and uncertainty analysis. Of the remaining models, Huss and Farinotti use the data for calibration and Bedmap2 and Bedmachine v2 are, in this region, basically interpolation approaches. Therefore, the three latter models would be those most benefited from the availability of additional ice-thickness observations. Indeed DeepBedMap would benefit from such additional data if they were used for training of its neural network.

This scarcity of suitable ice-thickness data in the region stresses the importance of initiatives such as the RINGS expert group of SCAR (https://www.scar.org/science/rings/home/), basically aimed at improving the knowledge of bed topography at the ice-sheet margin, is critical for accurate estimates of ice discharge. This is to be done through the promotion of internationally coordinated data collection campaigns, primarily focused on the Antarctic regions with severe data scarcity, like the AP, as well as developing a set of protocols to systematically collect, analyze and share comprehensive airborne geophysical data.

Another possible improvement for the ice-discharge calculation in ice-shelf terminating glaciers is using flux gates downstream of the grounding line in combination with a floating criterion. It requires accurate surface elevation information, since the surface elevation error will be multiplied by a factor of 10. Additionally, the ice discharge must also be corrected by the SMB downstream of the grounding line, adding further uncertainty to the calculation.

Finally, it is expected that the use, with the currently available data, of more sophisticated ice-thickness inversion models (e.g. using higher-order dynamical models) could improve the ice discharge results for the AP. In particular, we plan to use in a future paper a different type of model (Fürst and others, 2017). However, we cannot forecast how remarkably the results will improve, since no ice-thickness inversion model can work efficiently without the availability of proper calibration data.