Analysis of Antarctic Peninsula glacier frontal ablation rates with respect to iceberg melt-inferred variability in ocean conditions

Abstract

Marine-terminating glaciers on the Antarctic Peninsula (AP) have retreated, accelerated and thinned in response to climate change in recent decades. Ocean warming has been implicated as a trigger for these changes in glacier dynamics, yet little data exist near glacier termini to assess the role of ocean warming here. We use remotely-sensed iceberg melt rates seaward of two glaciers on the eastern and six glaciers on the western AP from 2013 to 2019 to explore connections between variations in ocean conditions and glacier frontal ablation. We find iceberg melt rates follow regional ocean temperature variations, with the highest melt rates (mean ≈ 10 cm d⁻¹) at Cadman and Widdowson glaciers in the west and the lowest melt rates (mean ≈ 0.5 cm d⁻¹) at Crane Glacier in the east. Near-coincident glacier frontal ablation rates from 2014 to 2018 vary from ~450 m a⁻¹ at Edgeworth and Blanchard glaciers to ~3000 m a⁻¹ at Seller Glacier, former Wordie Ice Shelf tributary. Variations in iceberg melt rates and glacier frontal ablation rates are significantly positively correlated around the AP (Spearman's ρ = 0.71, p-value = 0.003). We interpret this correlation as support for previous research suggesting submarine melting of glacier termini exerts control on glacier frontal dynamics around the AP.

1. Introduction

Many marine-terminating glaciers and ice shelves along the Antarctic Peninsula (AP) have undergone significant thinning and retreat in the last three decades (e.g. Liu and others, 2015; Cook and others, 2016; Rott and others, 2018). Following the collapse of ice shelves along the Antarctic Peninsula (AP), ice discharge increased from the marine-terminating glaciers that were formerly buttressed by the ice shelves (i.e. dynamic mass loss) (Rignot and others, 2004a, 2004b). The observed retreat and associated increase in dynamic mass loss have been coincident with regional oceanic warming, suggesting that this warming may have played a strong role in the ice-shelf collapse, glacier grounding line retreat and glacier area changes along the AP (Cook and Vaughan, 2010; Pritchard and others, 2012; Cook and others, 2016). However, the timing and magnitude of glacier terminus retreat and acceleration have varied around the AP (Cook and others, 2005; Pritchard and Vaughan, 2007; Wouters and others, 2015; Friedl and others, 2018), potentially in part due to local variations in ocean forcing. For example, following the collapse of the Wordie Ice Shelf between the 1960s and 1990s, the grounding line of Fleming Glacier (a former ice-shelf tributary) retreated ~6–9 km through 2011 and its flow increased by 27% between 2008 and 2011 (Friedl and others, 2018). This retreat has been attributed to enhanced submarine melt caused by increased upwelling of Circumpolar Deep Water (CDW), which was then sustained through internal feedbacks associated with the retrograde bed geometry (Friedl and others, 2018). At Seller and Prospect glaciers in the Marguerite Bay region, the >100 m a⁻¹ glacier acceleration from 2008 to 2015 (Gardner and others, 2018) has been linked with ocean warming related to large-scale changes in the Southern Oscillation Index and El Niño Southern Oscillation (Walker and Gardner, 2017). Meanwhile, on the eastern Antarctic Peninsula (EAP), oceanic warming may have been a factor in preconditioning the Larsen A and B ice shelves to their rapid surface meltwater-driven collapse in 1995 and 2002, respectively (Shepherd and others, 2003; Glasser and Scambos, 2008; McGrath and others, 2012). Taken together, these observations suggest that variations in ocean forcing of glaciers and ice shelves across the AP may strongly influence glacier dynamics in this region.

Despite their morphological differences, glaciers that feed ice shelves and those that terminate at grounded marine cliffs are both sensitive to ocean change through its influence on submarine melting. Submarine melting varies with both ocean water temperature and the flux of buoyant subglacial discharge emitted from the grounding line (Jenkins, 2011). For ice shelves, enhanced submarine melting thins the ice. If this thinning reduces the ability of the ice shelf to buttress flow from its tributary glaciers, this thinning will lead to acceleration and longitudinal extension, which can increase crevassing and trigger ice-shelf collapse (Glasser and Scambos, 2008; Pritchard and others, 2012). As observed at the glaciers that formerly fed the Larsen A and B ice shelves, ice-shelf collapse results in widespread acceleration, thinning and dynamic mass loss (Scambos and others, 2004; Hulbe and others, 2008). These dynamic changes can persist for years to decades, with glacier speeds on the EAP relatively unchanged from 2000 to 2015 but still in excess of pre-collapse speeds (Gardner and others, 2018). For glaciers that lack floating ice tongues, submarine melting occurs along the glaciers' near-vertical submerged margins. The combination of relatively cold subsurface water temperatures observed around the AP (Moffat and Meredith, 2018) and small meltwater runoff fluxes from AP glaciers (Van Wessem and others, 2016) suggest that submarine melt rates for marine-terminating termini may be much lower than the meters per day melt rates that have been inferred across the Arctic (Bartholomaus and others, 2013; Rignot and others, 2015; Holmes and others, 2019). The correlation between the magnitude of recent marine-terminating glacier retreat and subsurface water temperatures on the nearby continental shelves suggests that ocean-driven changes in submarine melting may exert a strong control on marine-terminating glacier dynamics along the AP (Pritchard and others, 2012; Cook and others, 2016).

Despite the potential importance of spatio-temporal variations in ocean conditions on dynamic mass loss to the AP, the present understanding of the importance of controls on glacier dynamics in this region is limited by sparse ocean observations proximal to glacier margins and/or estimates of submarine melt rates. Argo floats provide valuable information about spatio-temporal variations in regional ocean properties (i.e. temperature and salinity) in the Southern Ocean (e.g. Gille, 2008), but these far afield properties may differ from ocean conditions near glacier termini. As such, these data are not sufficient to tease-out the relative importance of submarine melt variability and other climate controls on dynamic mass loss from these glaciers.

Here, we explore the use of remotely-sensed iceberg melt rates near marine-terminating glaciers around the AP as a means to infer variability in ocean forcing of glacier terminus position change (Fig. 1). We show that iceberg melt rates generally follow variations in regional ocean temperatures and share a positive, statistically significant relationship with glacier frontal ablation rates (i.e. produced by the combined effects of submarine melting and iceberg calving). Although we do not directly quantify glacier submarine melt rates, as we are unable to with the data we have available, our analysis supports the important role of the ocean as a control on glacier dynamics around the AP.

Fig. 1. Location of study sites around the Antarctic Peninsula. The name for the glacier at each site is listed in the legend, followed by its location with respect to the spine of the peninsula (W = west, E = east). The size of the symbol corresponds to the mean iceberg melt rate proximal to the glacier terminus. Symbols are overlain on the 250 m-resolution Landsat Image Mosaic of Antarctica (LIMA). The inset in the lower left corner shows the full LIMA, with the red box outlining the study region.

2. Data and methods

2.1. Study sites

Eight study sites were selected proximal to glaciers along the AP: six along the west side near Leonardo Glacier and Blanchard Glacier on the Danco Coast, Cadman Glacier on the Graham Coast, Widdowson Glacier on the Loubet Coast, and Heim Glacier and Seller Glacier in the north and south Fallières Coast; and two from the Larsen A (Edgeworth Glacier) and Larsen B (Crane Glacier) embayments on the east side (Fig. 1). Sites were selected based on WorldView imagery availability and the repeat presence of icebergs in the 17 km-by-17 km image footprint. Persistent cloud cover along the western side of the peninsula, combined with the rapid evacuation of icebergs from the fjords in austral summer, restricted the intra- and inter-annual coverage of iceberg data for our six study sites in this region (Table 1). Along the eastern peninsula, there were a larger number of repeat iceberg observations but intra-annual variations in iceberg melting could not be resolved using the method described below (Table 1). Melt rates were calculated for 14 separate observation periods during 2013–2019 on the EAP and 2014–2019 on the western Antarctic Peninsula (WAP) for no more than 20 icebergs at each site. Melt rates were produced over approximate annual time periods on the EAP, with the exception of one observation for the Crane Glacier study site. WAP melt rates were calculated over time intervals ≤6 months, with the exception of the Seller Glacier study site.

Table 1. Iceberg melt data for all study sites and date pairs

Columns 3–8 are the number of shallow-drafted and deep-drafted icebergs, mean melt rates and corresponding std dev., median melt rates and the interquartile range for melt rates, respectively. Mean melt rates are the slope of the linear polynomials describing the relationship between submerged area and meltwater flux, with one std dev. confidence bounds. Mean and median melt rates with asterisks are reported without the two anomalously high iceberg meltwater fluxes circled in Figure 4.

2.2. Iceberg melt rates

Changes in iceberg surface elevation (i.e. iceberg freeboard) over time were used to estimate iceberg meltwater fluxes, and then converted to melt rates following the methods of Enderlin and Hamilton (2014). Briefly, we used very high-resolution stereo satellite images collected by the WorldView satellite constellation from 2013 to 2019 to create 2 m-resolution, high-precision (~3 m vertical uncertainty) iceberg digital elevation models (DEMs) using the NASA Ames Stereo Pipeline (ASP) software package (Shean and others, 2016). Vertical biases in iceberg elevations due to satellite orbital parameter uncertainties and tidal change between sequential DEMs were quantified using open ocean or thin sea ice elevations in close proximity to the icebergs. These biases were used to vertically co-register the iceberg DEMs so that elevations were calculated with respect to a constant reference (i.e. sea level at zero meters elevation). Changes in iceberg freeboard were then manually extracted from the DEMs for the largest icebergs that were easily recognizable from one time period to the next, typically 10 000–500 000 m² in surface area. Differences in elevation over time were used to estimate volume change, ΔV, as

(1)$$\Delta V = A_{{\rm surf}}\left[{\left({\displaystyle{{\rho_{\rm w}} \over {\rho_{\rm w}-\rho_{\rm i}}}} \right)\Delta h} \right]= A_{{\rm surf}}\Delta H\comma \;$$

where A _surf is the mean iceberg surface area, ρ _w and ρ _i are the density of the sea water and iceberg, respectively, Δh is the observed median change in freeboard and ΔH is the median thickness change.

An ocean water density of 1027.5 kg m⁻³ (Rackow and others, 2017) and uncertainty of 2 kg m⁻³ were used in all calculations. Iceberg density values varied both within and between study sites due to variations in iceberg shape and firn density. To estimate iceberg densities, we extracted firn density profiles from the pixel closest to the icebergs in the Institute for Marine and Atmospheric Research (IMAU) firn densification model (Ligtenberg and others, 2011). Exponential curves were fit to these discrete density profiles such that the density at any depth (z) could be estimated from 917 − (917 − a)^−z/b, where 917 kg m⁻³ is the density of bubble-free ice and a and b are the free parameters. These curves were used to estimate dry firn and ice densities (i.e. density of tabular icebergs above the waterline). The density of water-saturated firn/ice was estimated from these curves under the assumptions that (1) all pore space was filled, (2) pore close-off occurred at a density of 830 kg m⁻³, and (3) the maximum possible density was equal to sea water. Uncertainties in firn density were estimated from the ±1σ (i.e. 68% or ±1 std dev.) confidence interval for the density profiles and an ice density uncertainty of 10 kg m⁻³. For non-tabular icebergs, there was considerably larger uncertainty in the iceberg density as it was uncertain whether the iceberg consists solely of water-saturated firn, solely of bubble-free ice or some combination of water-saturated firn and ice. The mean and range of iceberg densities assuming water-saturated firn, as well as pure ice, were used for our elevation to volume conversions. The median density along the WAP varied from ~850 kg m⁻³ for Seller and Cadman glaciers to 877 kg m⁻³ for Leonardo Glacier, and the median density for icebergs calved from Edgeworth and Crane glaciers along the EAP was 897 and 886 kg m⁻³, respectively, reflecting the predominance of over-turned icebergs in our analysis.

Changes in volume estimated using Eqn (1) are the result of surface meltwater runoff, creep thinning and submarine melting (Enderlin and Hamilton, 2014). Structural disintegration of icebergs is not accounted for in this equation since the icebergs selected for our analysis (1) were surrounded by sea ice for the majority of the observation period, limiting the wave erosion that strongly influences the iceberg ‘footloose’ decay mechanism (Wagner and others, 2014) and (2) lacked obvious changes in surface area between repeat images. The contribution of surface meltwater runoff to the observed volume change was estimated from 5.5 km-resolution Regional Atmospheric Climate Model (RACMO2.3) outputs for the AP (e.g. Van Wessem and others, 2016) and subtracted from the observed volume change. Thinning due to viscous creep of the iceberg's unconfined lateral margins (Thomas, 1973) was also subtracted from the observed volume change. The creep calculation necessitates an estimate of the temperature-dependent rate factor, which is used to describe the ice viscosity. Icebergs are unlikely to have reached equilibrium with coastal air temperatures given the fast flow-velocities of AP glaciers, but are also unlikely to have maintained the same temperature as the snow that falls on the AP's high elevation accumulation zone (~−17°C annual average across the AP spine). In the absence of direct observations of iceberg temperatures, we estimated that iceberg internal temperatures were 5°C colder than the RACMO2.3 annual average air temperatures at the iceberg locations. The residual volume change estimates were converted to meltwater equivalents and divided by the time between imagery dates to produce submarine meltwater flux estimates (ΔV/Δt). To normalize these data for comparisons between study sites, we divided meltwater flux by the submerged ice area to produce an iceberg melt rate per unit submerged area (cm d⁻¹). Icebergs were assumed to have cylindrical submerged geometries that were constrained by the surface area, perimeter at the waterline and median thickness under the assumption of hydrostatic equilibrium. Uncertainties and potential biases in meltwater fluxes and melt rates estimated using this method are described in detail in Enderlin and Hamilton (2014) and Section 2.4 below.

2.3. Frontal ablation

Frontal ablation (FA), defined as glacier mass loss due to calving and submarine melt (though neither was directly measured here), was calculated as the difference between surface speed and the terminus position change rate along each glacier centerline (Fig. 2). The terminus change and speed datasets are described below.

Fig. 2. Schematic showing glacier frontal ablation, where ΔL/Δt is the change in terminus position over time, U is the velocity and FA is the frontal ablation. The green terminus position shows a scenario of glacier advance, light blue shows a scenario of no terminus position change and dark blue shows a scenario of glacier retreat. Not to scale.

Terminus position time series were constructed for the primary glaciers that contribute icebergs to each study site. Satellite images from Landsat-4,5,7,8 and Sentinel-1 and 2 acquired between 2012 and 2019 were used to map terminus position change using the Google Earth Engine Digitisation Tool (GEEDiT) (Lea, 2018). Terminus positions were recorded monthly during the austral summer (October to March) in all available imagery for which a terminus could be clearly delineated. The spatial and temporal coverage of imagery varied, with more dense coverage after the launch of Landsat-8 in 2013 and the Sentinel satellites in 2014/2015. We manually traced the centerline of each glacier and extracted terminus position change time series along the centerline. We chose the centerline method, rather than the more complex methods (Lea and others, 2014), because inter-annual changes in terminus position observed over the study period were relatively uniform across the glaciers. Additionally, Walsh and others (2012) found that the centerline and box methods resulted in nearly identical estimates of terminus position change for glaciers spanning a wide range of geometries, supporting our use of the centerline method for quantification of terminus position change and FA rates.

Glacier speed time series were extracted from the nearest neighboring pixel ~300 m inland of the most-retreated terminus centerline position. For all study sites, glacier speeds were obtained from the GoLIVE dataset. GoLIVE velocities were constructed using feature-tracking techniques applied to Landsat optical images. GoLIVE velocity maps with 300 m spatial resolution were available from 2013 to 2018 (Fahnestock and others, 2016). Inspection of the GoLIVE data revealed that although the data have sufficient temporal coverage for our analysis, velocity time series extracted near the glacier termini show that clouds introduce considerable scatter into the WAP velocity time series. In order to filter the erroneous data, we added velocity observations from the Landsat feature-tracking-based ITS_LIVE dataset (Gardner and others, 2018) as well as from two TerraSAR interferometry-based datasets which extended back in time to as early as 1995 (Wuite and others, 2015; Rott and others, 2018). The ITS_LIVE data have a slightly higher spatial resolution (240 m-resolution) but have an approximate annual temporal resolution. The 50 m-resolution TerraSAR-X velocities have a comparable temporal resolution to the GoLIVE datasets (~monthly) but are restricted to Edgeworth Glacier and Crane Glacier study sites from 2013 and 2016. To filter the data, we use an annual moving window centered about the median ±3 × 1.48 times the median of the absolute deviation (MAD) of the entire ITS_LIVE speed time series. For normally distributed data, this filtering envelop is analogous to the mean ±3 std dev. For each year, we first calculate the median annual speed using only the ITS_LIVE and TerraSAR speeds. If the GoLIVE annual median speed falls within the median ±3MAD envelope, we recalculate the median using the GoLIVE data to ensure that seasonality not resolved by ITS_LIVE is accounted for in our filtering. This process effectively removes unrealistic outliers while preserving inter-annual temporal variations in speed. Potential speed (and FA) biases associated with differences in spatial resolution of the various velocity datasets and the use of a fixed velocity sampling location are described in Section 2.4 and summarized in Table 2.

Table 2. Potential speed biases due to differences in spatial resolution of the velocity datasets and the use of velocities from a fixed inland location

Positive values in columns 3 and 4 indicate ITS_LIVE and TerraSAR speeds, respectively, were less than GoLIVE speeds from overlapping time periods. The GoLIVE velocity dataset is used as a reference due to its extensive temporal coverage. The average along-flow change in speed and per cent acceleration between the fixed inland location used for our frontal ablation calculations and the glacier termini is included in the last column.

We constructed the time series of FA rates from around 2014 to 2018 for all study sites. Although the terminus position record extends for decades, we restricted our FA estimates to time periods with sufficient coverage to resolve variations in FA at the same temporal resolution as our melt rate estimates. Terminus change rates were calculated by dividing terminus change by the time period between the first and last terminus observations with approximately the same temporal resolution of our iceberg melt rate estimates: 3 months for the WAP and 1 year for the EAP. Assuming ice flow is primarily due to basal sliding, surface speeds were treated as depth-averaged speeds such that FA was calculated as the difference between the average surface speed minus the corresponding terminus change rate.

2.4. Uncertainties and potential biases

The analysis of spatio-temporal variations in iceberg melting and FA, as well as comparisons between these variables, is limited by uncertainties and potential biases in our estimates. Uncertainties and biases in iceberg melt and FA estimates are presented below.

Quantifiable sources of uncertainty in our meltwater flux (ΔV/Δt) estimates include random errors in surface elevations, uncertainty in iceberg mapping and uncertainty in iceberg density estimates. WorldView DEMs have a random uncertainty of ~3 m (Enderlin and Hamilton, 2014; Noh and Howat, 2015; Shean and others, 2016). Uncertainty in the manual mapping of icebergs introduces uncertainty through both its influence on the estimates of surface elevation change and area of the icebergs. Each iceberg was manually delineated ten times and interpreter uncertainty was quantified from the std dev. in elevation change estimates constructed from the repeat delineations (Enderlin and Hamilton, 2014). Changes in iceberg surface area over time also introduce uncertainty into our estimates, as these can result from differences in iceberg delineation accuracy as well as real changes in iceberg extent. The temporal difference in the average surface area of the repeat delineations was used to constrain area change uncertainty. Uncertainties in iceberg densities were computed as described above. These quantifiable uncertainty terms were summed in quadrature to estimate meltwater flux uncertainties.

The median meltwater flux uncertainty for the WAP was 0.025 m³ s⁻¹, while this median uncertainty on the EAP was 0.001 m³ s⁻¹. Interpreter and random errors dominated the meltwater flux uncertainties. For the WAP, interpreter and random errors typically contributed 46 and 7%, respectively. For the EAP, interpreter and random errors typically contributed 65 and 27%, respectively. Interpreter and random errors were higher on the EAP likely because little change in elevation occurred from one time period to the next in that region, increasing sensitivity to small variations in iceberg delineation accuracy. Uncertainties associated with temporal changes in the surface area of the iceberg accounted for 5% of WAP uncertainty and were negligible (0.2%) for the EAP. Surface area changes were larger for the WAP, likely due to mass loss via iceberg fragmentation following the seasonal loss of sea ice. Uncertainties related to density variations were negligible (median values of <0.1% for the EAP and 0.1% for the WAP).

There are several sources of potential biases in our iceberg meltwater flux and melt rate estimates that were less easily quantified. RACMO2.3 estimated no surface meltwater runoff at all sites, so the 30% uncertainty quoted for RACMO does not contribute to meltwater flux uncertainty. Given the observations of surface melting on glaciers and ice shelves around the AP (Scambos and others, 2000; Tuckett and others, 2019), it is likely that the coarse resolution of RACMO and the steep topographic relief in this region result in under-estimates of surface melting. However, we did not attempt to downscale the RACMO outputs to better estimate mass loss as we did not have in situ data to constrain any statistical downscaling approach. The potential biases introduced to our iceberg melt estimates as a result of the apparent lack of surface melting vary across the AP. Along the WAP, there were no visible surface melt features in our austral spring satellite images and we assume that the estimates of zero surface melting (from RACMO outputs) have a negligible influence on our iceberg submarine melt estimates. Surface meltwater runoff likely occurs along the EAP, as supported by our interpretation of subglacial meltwater plume-enhanced iceberg melt rates described below, but there were no in situ observations with which to constrain its contribution to the observed iceberg volume change. Creep thinning also constitutes a potential source of bias. We estimated a median creep thinning of 0.07 m (0.12 m) for icebergs along the WAP (EAP), constituting ~1% (~34%) of observed iceberg thickness change for icebergs at those sites. Although we could modify the temperature of the ice to see how creep thinning is influenced by temperature uncertainty, there were no in situ constraints on ice temperature on which to base these temperature variations. Thus, we cannot confidently assess whether our creep thinning is accurate or represents a potential over- or under-estimate of its contribution to thickness change, and therefore volume change. Finally, we used simple cylindrical geometries for our iceberg area estimates. Iceberg shapes almost certainly deviate from the assumed geometries, as submerged ice area is generally more complex than a cylinder, but we cannot quantify these variations from surface observations. Even if we used a variety of simple geometries, such as cones or ‘skirted’ icebergs, we cannot quantify surface complexity/roughness and potential impacts on melt rates. Based on the complex scalloped geometries of capsized icebergs, we recommend that our melt rates are treated as maximum estimates.

Speed uncertainties were provided with each velocity dataset. GoLIVE velocity uncertainties vary from ~0.5 to 1.0 m d⁻¹ and we conservatively estimated uncertainties of ±1 m d⁻¹ for GoLIVE velocities (Fahnestock and others, 2016). ITS_LIVE uncertainties were provided for each pixel and can vary widely, but are generally less than the GoLIVE uncertainties. ENVEO Cryoportal (http://cryoportal.enveo.at/data/) velocities were estimated as ±0.05 m d⁻¹ (Rott and others, 2018). The additional TerraSAR-X velocity datasets used to expand on the Crane and Edgeworth time series have an average uncertainty of ±0.08 m d⁻¹. We compared nearly-coincident velocities from the datasets to quantify potential biases introduced in our FA estimates as the result of variations in temporal coverage of the three velocity datasets. All datasets were compared to GoLIVE as its temporal coverage overlaps the other datasets. Biases vary by site, with median offsets between GoLIVE and the other speed datasets in Table 2.

To quantify potential FA biases introduced by the use of speeds from a fixed location rather than speeds from the terminus, we extracted along-flow changes in speed from the centerline between our fixed inland location and the seaward-most observation for each velocity map. The use of speeds from a fixed inland location results in an average under-estimation of speed by 242 m a⁻¹ (33.7%) for the EAP and 40 m a⁻¹ (4.4%) for the WAP (Table 2).

Uncertainties in glacier speed and terminus position, which were assumed to be ±15 m (one pixel, falling below image resolution) following previous studies (e.g. Burgess and Sharp, 2004; Carr and others, 2013), were summed in quadrature to obtain the estimates of FA uncertainty. We did not apply adjustments to the FA data to account for the biases presented in Table 2 since these biases can vary in space and time; instead, we discuss the effects of these biases below.

3. Results

3.1. Iceberg melt rates

Iceberg melt rates vary spatially (Fig. 1; Fig. 3; Table 1), with higher melt rates on the western side of the peninsula and lower melt rates on the eastern side of the peninsula. In line with iceberg melt rate estimates for Greenland (Enderlin and others, 2018; Moon and others, 2018), we find that iceberg melt rates generally increase with the draft. To minimize the influence of iceberg size on our interpretation, we focus our analysis on iceberg melt rates estimated as the slope of linear polynomials fit the submerged area and meltwater flux estimates for each observation period (Enderlin and others, 2018), hereafter referred to as the mean melt rate (Table 1). The 68% (i.e. ±1 std dev.) confidence interval for these polynomials, as well as the median melt rates and interquartile ranges in melt rates, are reported in Table 1. The lowest mean melt rate of all observation periods for deep-drafted icebergs was 0.08 cm d⁻¹ at Crane Glacier in 2016–2017. The highest mean melt rate was 10.57 cm d⁻¹ at Cadman Glacier in 2018–2019. We found no consistent patterns in iceberg melt rates with respect to proximity to terminus position. We do, however, note anomalously large melt rates for two icebergs proximal to Edgeworth Glacier on two different dates and we exclude these outliers from the melt rate metrics in Table 1.

Fig. 3. Individual iceberg (a) meltwater fluxes and (b) melt rates at study sites along the Antarctic Peninsula from 2013 to 2018 (see Table 1). Vertical error bars account for uncertainty in elevation, surface area, density and surface melting, summed together in quadrature. The horizontal error bars represent uncertainties associated with surface elevation and area change over time.

We do not have sufficient coverage for detailed analysis of temporal variations in melt rates, with the exception of Edgeworth and Crane glaciers on the EAP (Fig. 4a). The size distributions for the icebergs at these sites vary in space and time, however, making a direct comparison of melt rates difficult. Differences in mean melt rate that exceed the bounds of the uncertainty envelopes in Figures 4b, c are considered significant. For Edgeworth Glacier, the mean iceberg melt rate from 2013 (austral spring) to 2014 (austral winter) was significantly higher than the mean melt rates from the other three observation periods, with a value of 1.85 cm d⁻¹ excluding outliers (Fig. 4b). The interpretation of these temporal differences is highly sensitive to the small sample size of the dataset: if we remove the next two largest icebergs from the 2013–2014 dataset, we find the average melt rate for this observation period drops to 1.0 cm d⁻¹ and is no longer significantly different than the other observation periods. Potential drivers of these anomalous melt rates are discussed in detail below. For Crane Glacier, we found that iceberg melt rates do not significantly differ in time, with most observations falling within the one-sigma confidence interval for the 2015 (austral spring) to 2016 (austral summer) observations.

Fig. 4. Iceberg meltwater fluxes plotted against the estimated submerged area for eastern Antarctic Peninsula study sites. (a) All meltwater flux observations for Edgeworth (triangles) and Crane (circles). Symbol colors distinguish observation year. (b) Meltwater fluxes for the Edgeworth study site. Shaded regions outline the ±1 std dev. confidence interval for the linear trendlines. The two circled, anomalously large meltwater fluxes for Edgeworth are excluded. (c) The same as (b) except at the Crane study site. For all panels, vertical error bars represent uncertainties associated with elevation and surface area observations as well as density, surface meltwater runoff and creep thinning estimates. Horizontal error bars represent uncertainties associated with surface elevation and area change over time. Note that each panel has a different scale.

3.2. Frontal ablation results

3.2.1. Terminus change

Figure 5 shows annual terminus positions from the austral summer for each study site. Terminus geometry remained fairly constant for each study site, such that the centerline terminus position change rate time series in Figure 6 can be considered representative of the width-averaged terminus retreat or advance. Temporal patterns in terminus change varied between study sites. Crane Glacier and Cadman Glacier termini advanced >3 km from 2013 to 2019. Edgeworth Glacier advanced until 2017, then retreated back to its 2014 position. Heim gradually retreated throughout the study period, with a total retreat of ~0.6 km. Widdowson, Blanchard, Leonardo and Seller Glaciers maintained relatively stable terminus positions (<0.5 km change) in recent years. All glaciers were characterized by terminus position change rates at the same rate as ice flow throughout the majority of the year and periods of rapid retreat during the austral summer (Fig. 6). Maximum retreat rates of >200 m d⁻¹ were observed at Blanchard and Seller glaciers. Short-term rates of terminus retreat are highly sensitive to the temporal sampling of the satellite images, however, and we suspect that comparably high terminus position change rates are possible at the other study sites given the prevalence of large (>0.5 km-wide) icebergs in the image record.

Fig. 5. Summer glacier terminus positions for all study sites (see legend). Maps are displayed on the same spatial scale. The white X's indicate iceberg locations, with marker sizes scaled according to melt rate for all icebergs measured from 2013 to 2019.

Fig. 6. Speed (left axes) and terminus position change rate (right axes, squares) time series for the major outlet glacier at each study site. Markers for speeds distinguish data sources: stars = TerraSAR, diamonds = GoLIVE, triangles = ITS_LIVE. Marker face colors distinguish location (map inset) and the outline colors distinguish the (black) speed and (pink) terminus position change rate datasets. WAP sites are shown in the left column and EAP sites are shown on the right column. Vertical error bars represent uncertainties in the speed and terminus position change rate datasets and horizontal error bars indicate the time period over which speeds and terminus position change rates were calculated. The map inset illustrates spatial variations in near-terminus speed, with marker sizes scaled according to the median speeds in legend.

3.2.2. Glacier speeds

Glacier speeds at all study sites remained relatively stable (within uncertainty bounds) since 2014 (Fig. 6). However, the speed time series suggest Cadman and Widdowson glaciers along the WAP accelerated by ~1 m d⁻¹ (from ~5 to 6 m d⁻¹) from 2014 through 2018. No multi-year trends in speed are apparent for Leonardo, Blanchard and Heim glaciers, also along the WAP. Seller Glacier along the WAP and Edgeworth Glacier along the EAP slowed by ~0.6 m d⁻¹ from 2013 to 2018. The most pronounced deceleration (1.1 m d⁻¹) occurred from 2013 through 2017 at Crane Glacier along the EAP, followed by an apparent stabilization in speed.

3.2.3. Frontal ablation

In line with the temporal resolution of our iceberg melt observations, we calculated FA rates over 3-month time periods for study sites along the WAP – referred to here as ‘seasonal FA’ – and over annual time periods for study sites along the EAP – ‘annual FA’. We observe spatially-variable FA rates over the study period (Fig. 7; Table 2). Seller and Widdowson Glaciers have the highest FA rates, with mean rates of 2958 and 2008 m a⁻¹, respectively. Blanchard and Edgeworth Glaciers have the lowest FA rates, with mean rates of 442 and 543 m a⁻¹, respectively. The rest of the sites on the WAP have mean FA rates ranging from 585 to 1531 m a⁻¹.

Fig. 7. Frontal ablation rates at the major outlet glacier at each study site. Three-month and 1-year frontal ablation rates are plotted for the WAP sites in the left column and EAP sites in the right column, respectively. Vertical error bars represent uncertainty and horizontal bars indicate the time period over which rates were averaged. The map inset illustrates spatial variations in the frontal ablation rate, with marker sizes scaled according to the median frontal ablation rates in the legend.

Across the study period, there is a clear relationship between FA rate and speed: average FA rates are larger for faster-flowing glaciers (Fig. 6; Table 2). In line with the multi-year stability or decrease in speed observed at most glaciers, FA rates remained relatively consistent or decreased between 2014 and late 2018 across the AP (Fig. 7). Leonardo, Blanchard and Seller reached their maximum FA rates in 2014, followed by Heim Glacier which reached peak rates in mid-2015, while Cadman and Widdowson Glaciers reached their maximum FA rates from mid-2016 to mid-2017; however, none of these variations exceed the uncertainties in FA rates. On the EAP, maximum annual FA rates occurred from 2004 to 2005.

4. Discussion

4.1. Ocean conditions along the peninsula and iceberg melt

At regional scales, spatial variations in mean and median iceberg melt rates are generally in agreement with shelf temperature variations (Schmidtko and others, 2014; Van Caspel and others, 2015; Moffat and Meredith, 2018). Water temperatures along the EAP and northern WAP are relatively cold as they are sourced from the Weddell Sea, which has been cooling gradually since the 1970s and is characterized by very cold (~−1.8°C) surface waters (above 500 m) on the shelf (Schmidtko and others, 2014; Van Caspel and others, 2015). The cold water temperatures on the EAP explain the relatively homogeneous, low iceberg melt rates at the Edgeworth Glacier and Crane Glacier sites, with the exception of the very large melt rates from icebergs proximal to the Edgeworth Glacier margin (Figs 4a, 5f), where we hypothesize glacial meltwater plumes exist. Meltwater plumes could enhance turbulent melting, driving the approximate order of magnitude difference in melting between the large icebergs proximal to the glacier terminus and those farther afield (Xu and others, 2012; Moon and others, 2018). Ignoring these anomalously large melt rates, there is no temporal variation in iceberg melt rates that exceeds our uncertainties on the EAP.

Along the WAP, our observed iceberg melt rates reflect regional variations in the vertical hydrographic structure, as discussed in Moffat and Meredith (2018, Fig. 3). Just east of Cape Jeremy and near Adelaide Island, the thermal structure of the water is strongly stratified with warm (~0.75°C) Antarctic Surface Water (AASW) at the surface, very cold Winter Water (WW) (~−1.25°C) from ~50 to 150 m, water temperatures between 0 and 1°C from ~150 to 225 m and 1.25°C below ~225 m depth. Even for the simplified iceberg geometries used in our melt rate estimates, which may under-estimate the true draft of the icebergs, we find that Seller Glacier icebergs are sufficiently deep-drafted during the 2014–2016 and 2019 observation periods (with median iceberg drafts of 121 and 167 m and maximum drafts of 192 and 244 m in 2014–16 and 2019, respectively) to contact with the warmer subsurface waters, explaining the relatively high mean melt rate at that site. Lower melt rates at Seller during 2016–2017 correspond with the inclusion of shallower icebergs in our analysis: the median iceberg draft of 68 m (maximum draft = 104 m) does not penetrate the lower boundary of the WW layer. This stratification also likely contributes to the relatively low iceberg melt rates at Heim Glacier, given that the icebergs sampled have a maximum draft of only ~73 m and may not penetrate through the colder WW layer. Moving north along the WAP, spatial patterns in iceberg melt rates also align with variations in the temperature structure on the shelf: iceberg melt rates are highest for the relatively deep icebergs (up to 173 m draft) observed at Widdowson and Cadman glaciers, corresponding to a region where AASW expands deeper into the water column (to ~50 m depth) and WW is compressed. Approximately half of the icebergs included at Widdowson Glacier are for icebergs with keel depths between 70 and 150 m, where a rapid transition occurs between −1°C WW and 0.75°C upper CDW. Although no icebergs appear to penetrate the warm (up to 1.75°C) CDW found in this region, these are the warmest waters on the shelf. Iceberg melt rates near Cadman Glacier are comparable to those estimated near Widdowson as a result of the compressed AASW and slightly warmer (~0.25 to −0.75°C) but expanded WW layer thickness. Moving north, iceberg melt rates are exceptionally low at Blanchard likely due to the shallow depths of the icebergs measured at this site (~40–60 m), and more enclosed geometry of the Blanchard Glacier fjord. Iceberg melt rates at neighboring Leonardo are comparably high as the melt rates nearby Cadman and Widdowson glaciers, likely due to the inclusion of deeper-drafted icebergs (up to 84 m) and the expansion of AASW and relatively warm (~0–0.25°C) WW.

While there is general agreement between our iceberg melt rates and regional ocean conditions, some local deviations warrant further discussion. The shallow draft of icebergs near Heim Glacier may at least partially explain their low melt rates. However, these melt rates are still surprisingly low given their proximity to the deeper bathymetry of the Marguerite Trough and observations of anomalously warm waters in this region as recently as 2012 (Venables and others, 2017; Walker and Gardner, 2017; Moffat and Meredith, 2018). We hypothesize that these relatively low iceberg melt rates reflect the deflection of warm CDW by bathymetric obstructions seaward of the Heim Glacier study site (e.g. Moffat and Meredith, 2018). While Heim Glacier is located proximal to Marguerite Trough, we infer that it terminates in an area of shallower bathymetry based on a review of the mapped bathymetry in the region (Arndt and others, 2013) and the presence of Blaiklock and Pourquoi Pas Islands nearby, which likely prevent warmer sub-surface waters from reaching the glacier. These results suggest that we cannot solely use far afield (i.e. shelf) observations to infer water mass properties at glacier termini.

4.2. Glacier frontal ablation comparisons to iceberg melt rates

We find that there is a moderately strong positive relationship (Spearman's correlation coefficient of 0.71, p-value = 0.003) between near-coincident mean iceberg melt rates and glacier FA rates (Fig. 8). Although our sample size is limited and there is considerable scatter in the relationship, we interpret the positive relationship between iceberg melt rates and FA rates as an indicator that ocean conditions influence FA of marine-terminating glaciers around the AP. Below we discuss potential reasons for the observed scatter, including controls on iceberg melt rates that differ from the controls on FA and the potential influence of temporal and/or spatial sampling limitations on our analysis.

Fig. 8. Scatterplot of iceberg melt rates and frontal ablation rates for nearby glaciers over near-coincident time periods. Symbols mark the median iceberg melt rate with horizontal bars indicating the interquartile range (Table 1). Vertical bars indicate uncertainty in the frontal ablation rate estimates. Symbol colors distinguish study sites (legend). The equation for the linear polynomial fit to the data (black line) is indicated in the top left, where m _iceberg is the iceberg melt rate in cm d⁻¹ and FA_glacier is the frontal ablation rate in m d⁻¹.

We first address the controls on iceberg melting and glacier FA rates. Our analysis selected the largest icebergs that did not noticeably fragment and remained within ~15 km of the glacier termini between observation periods. These criteria maximize the likelihood that spatio-temporal patterns in iceberg volume change are primarily controlled by subsurface oceanographic conditions. However, it remains possible that the icebergs and glaciers are exposed to different water masses, particularly for the smaller icebergs included at the Leonardo and Blanchard study sites. Additionally, it is possible that the primary controls on iceberg submarine melting and glacier FA rates differ. Both iceberg and glacier submarine melting are strongly controlled by the velocity of the water mass relative to the ice (Xu and others, 2012; Moon and others, 2018), but shear past icebergs can drive high rates of submarine melting even when buoyancy-driven shear along the glacier terminus is near zero (Moon and others, 2018).

Glacier FA rates are also controlled by processes other than submarine melting, including factors that change the resistive forces acting on the glacier terminus. Daily to decadal variations in buttressing provided by sea ice, ice mélange and floating ice shelves, as well as surface meltwater runoff-driven changes in basal sliding, can influence glacier FA rates. At Seller and Crane glaciers, for example, FA rates are likely influenced by the reduction in ice flow buttressing that accompanied the collapses of the Wordie Ice Shelf between 1966 and 1989 and Larsen B ice shelf in 2002, respectively. Following the ice-shelf collapses, tributary glaciers for both sites accelerated (Doake and Vaughan, 1991; Rignot and others, 2004b, 2005; Scambos and others, 2004; Wendt and others, 2010; Friedl and others, 2018; Gardner and others, 2018). For Seller Glacier, we observe relatively stable, elevated speeds from 2014 to 2018, indicating that the dynamic response to ice shelf persisted throughout our observation period. Crane Glacier underwent prolonged retreat from 2002 to 2007 before stabilizing at sustained high velocities (Gardner and others, 2018; Seehaus and others, 2018), which appear to drive FA rates today. Thus, we suggest that the FA rates at Seller and Crane glaciers are also a reflection of the glaciers' continued long-term dynamic adjustment to ice-shelf collapse. Glacier FA rates at our other study sites may also vary in part due to controls beyond submarine melting, and variations in such controls over the time scales of our observations also likely influence our iceberg melt rate estimates due to the highly-coupled nature of the glacier–ocean system. For example, variations in subglacial discharge will influence fjord circulation, in turn altering iceberg melt rates relative to periods when runoff is absent. Thus, although some of the scatter in the observed relationship between iceberg melting and glacier FA rates may be caused by long-term changes in glacier dynamics, it is highly unlikely that the positive relationship between iceberg melt rates and glacier FA rates is purely coincidental.

We next address the issues in temporal and spatial sampling. Our focus on FA rates and nearly coincident with submarine iceberg melt rates over ~3- to 12-month intervals minimizes potential aliasing due to large individual calving events and synoptic weather phenomena. Although the icebergs included in our analysis were mostly surrounded by sea ice in the satellite images used to construct our iceberg DEMs, the difference in iceberg melt and glacier FA rates for the WAP and EAP is potentially enhanced by the differences in temporal resolution (approximately seasonal vs annual) of the regions, but this cannot be assessed with the limited data available at present.

The spatial sampling of the iceberg melt estimates, varying spatial resolution of the glacier velocity data and spatial sampling of the glacier velocity datasets potentially introduce bias. The influence of iceberg location on our analysis is discussed with regard to controls on iceberg melting and glacier FA rates above. The potential biases introduced by using speed data with the varying spatial resolution are presented in Table 2. Our comparison of nearly-coincident speeds from the different datasets suggests that the speeds differ by up to ~10%. Unrealistic velocity data near the glacier termini also limited the along-flow location from which we extracted speeds, potentially resulting in slight under-estimation of FA rates (~1–10%) for the glaciers along the WAP. The average along-flow acceleration estimates in Table 2 suggest that FA rates may be underestimated by ~30% for the EAP. If the FA estimates are adjusted to account for the mean along-flow speed change along the centerline (Table 2), the relationship between glacier FA rates and iceberg melt rates weakens (Spearman's correlation coefficient of 0.61, p-value = 0.015), and the relationship – while still positive – becomes less so.

We interpret the observed relationship between iceberg melt rates and FA rates for nearby glaciers as support for the influence of ocean conditions on submarine melting of glacier termini and, in turn, glacier dynamics. We hypothesize that the statistical relationship between these variables would be improved with an increased number of coincident observations of glacier FA and iceberg melt rates. The positive correlation between iceberg melt rates and glacier FA rates, as well as the general agreement between spatial patterns in iceberg melting and water temperatures on the nearby shelf, provides additional support for the connection between ocean conditions and glacier dynamics along the AP presented in Cook and others (2016). Our data also suggest that the rates of submarine melting of icebergs are strongly controlled by subglacial discharge plumes: in locations where iceberg melt rates are likely enhanced by subglacial discharge plumes, such as for deep-drafted icebergs adjacent to the Edgeworth terminus, these plume-enhanced iceberg melt rates increase by an order of magnitude to maxima of ~5.5 cm d⁻¹, respectively. Glacier submarine melting is also sensitive to the presence of subglacial discharge, with the potential for considerable enhancement of terminus-wide submarine melt undercutting driven by local meltwater plumes (Slater and others, 2018). Thus, we interpret our results as an indicator of the influence of regional shelf-sourced and local glacier-forced variations in submarine melt conditions on AP FA rates and glacier dynamics, but suggest that more research still needs to be done to better quantify the effects of the ocean on AP glaciers.

Unfortunately, we do not presently have a long enough record of iceberg melt rates for a robust review of temporal variations in ocean conditions over time. The western Weddell Sea has cooled by ~0.05°C per decade since the 1970s (Meredith and King, 2005; Schmidtko and others, 2014; Cook and others, 2016) but this small decrease in ocean temperature over our 5-year study period is insufficient to drive change in iceberg melt rates at Edgeworth and Crane Glaciers that exceed our uncertainties. Despite the relatively stable iceberg melt rates (and FA rates) observed during our study period, continued comparisons between ocean conditions as inferred through iceberg melt rates and local glacier dynamics may shed light on the extent to which dynamic mass loss is forced by changing ocean conditions (Depoorter and others, 2013; Rignot and others, 2013; Luckman and others, 2015). Therefore, we recommend continued observations of iceberg melt rates and glacier dynamics so that the full potential of iceberg melt rates as a tool to infer changes in ice–ocean interactions can be explored using longer time series.

5. Conclusions

Correlations between spatio-temporal patterns in ocean conditions and glacier dynamics along the AP suggest that the marine-terminating glaciers in this region are sensitive to ocean conditions yet there are limited in situ ocean observations available to fully assess the role of the ocean as a control on dynamic mass loss. Using repeat high-resolution DEMs from 2013 to 2019 for icebergs around the AP, we show that the magnitude of submarine iceberg melt rates generally follows regional patterns in ocean conditions on the eastern and western sides of the AP, with exceptions where we infer the presence of shallow bathymetry. Mean FA rates were 543 and 1445 m a⁻¹ for the Edgeworth and Crane glaciers on the EAP, and ranged from 442 m a⁻¹ at Blanchard Glacier to 2959 m a⁻¹ at Seller Glacier on the WAP. Iceberg melt rates are positively correlated with nearly-coincident rates of FA at nearby glaciers, providing support for the hypothesis that glacier dynamics along the AP are linked to ocean conditions. Based on this analysis, we recommend that remotely-sensed iceberg melt rates continue to be used as a proxy for ocean conditions in remote areas where in situ measurements are logistically difficult and/or prohibitively expensive to collect. The generation of long-term records of iceberg melting will potentially enable a better understanding of the sensitivity of marine-terminating glaciers to ocean forcing, which is ever important for predicting the contributions of such glaciers to global sea level rise under climate change.