A comparison of contemporaneous airborne altimetry and ice-thickness measurements of Antarctic ice shelves

Abstract

Estimates of ice shelf mass loss are typically based on surface height measurements, assuming hydrostatic equilibrium and estimated firn thickness. Recent investigations, however, challenge the assumption that ice shelves are freely floating, particularly in proximity to narrow structures such as basal channels and shear margins. We compare contemporaneous measurements of Antarctic ice shelf thickness, from ice-penetrating radar, to freeboard height, from laser altimetry, acquired during multiple airborne surveys. On average, the hydrostatic thickness differs from observed thickness by at least ~17 ± 98 m, but this difference varies well beyond the propagated error within and among ice shelves, and depends on the corrections applied. We find that uncertainty in firn thickness can account for most, but not all, of the imbalance. Overall, errors in hydrostatic thickness do not significantly impact estimated basal melt rates. Our results indicate that localized approaches to estimating ice shelf thickness and rates of change are not applicable at large scales, and vice versa, and point to the need for more abundant and accurate firn and ice thickness measurements to improve estimates and predictions of ice shelf mass loss.

1. Introduction

Estimates of ice shelf mass loss are strongly dependent on the assumption of hydrostatic equilibrium, especially when ice thickness measurements are unavailable. The assumption that the ice shelf is freely floating allows estimates of ice thickness from abundant surface height measurements combined with estimates of firn thickness. Recent investigations, however, provide evidence that some areas of ice shelves are not freely floating, particularly in regions associated with steep gradients in ice thickness, such as basal channels and shear margins (e.g. Le Brocq and others, 2013; Drews, 2015; Drews and others, 2016; Alley and others, 2019; Chartrand and Howat, 2020; Dow and others, 2021; Wearing and others, 2021). However, the spatial scales of this imbalance and their impacts on ice shelf thickness and mass balance estimates are not fully understood.

The validity of the hydrostatic assumption has been investigated several times, though rarely with contemporaneous surface height and ice thickness data. The hydrostatic assumption was used to compare ice surface heights from the European Remote Sensing Satellite, ERS-1 (launched in 1991), to surface heights derived from Ross Ice Shelf Geophysical and Glaciological Survey (RIGGS) and Scott Polar Research Institute ice penetrating radar (IPR) surveys from the 1970s (Bamber and Bentley, 1994). They found agreement between the two datasets within the combined errors of the measurements across the Ross Ice Shelf, with some exceptions near grounding lines and in the vicinity of flow stripes. Later, ERS-1 data from 1994–1995 were supplemented by ICESat laser altimetry to derive a 1-km resolution gridded ice thickness dataset for all Antarctic ice shelves, and these results were compared to several independent airborne IPR-derived thickness datasets, with varying levels of agreement depending on proximity to the grounding line, data gaps and unknown marine ice density and thickness (Griggs and Bamber, 2011). Chuter and Bamber (2015) compared 1 km-resolution gridded ice thickness estimates derived from Cryosat-2 radar altimetry from 2011–2014 to IPR thickness measurements from 2001 and earlier on the Amery Ice Shelf, showing a mean difference in thickness of ~3% between the two estimates. Similarly, Fricker and others (2001) showed general agreement between hydrostatic thickness and IPR thickness measurements except where marine ice was expected, which they found accounts for about 9% of the ice shelf volume.

Observed and estimated ice shelf thickness from contemporaneous, or nearly contemporaneous, IPR and surface height measurements have been compared on smaller spatial scales, particularly in investigations of basal channels (e.g. Chartrand and Howat, 2020; Dow and others, 2021), showing that the hydrostatic assumption underestimates variability in ice thickness over distances <1 km. Similarly, simulations of stress fields on ice shelves have shown that the freeboard above basal channels is maintained at a higher height than expected based on ice thickness and the hydrostatic assumption, likely due to bridging stresses (Le Brocq and others, 2013; Drews, 2015).

Accurate measurements and/or estimates of ice shelf thickness, including near grounding zones, are crucial for estimates of mass balance. Disagreement between observed and hydrostatic thickness has consistently been identified near grounding lines, where ice is generally thinner than expected under the hydrostatic assumption, due to tidal flexure in the grounding zone (Bindschadler and others, 2011; Griggs and Bamber, 2011; Chuter and Bamber, 2015), leading to uncertainties in ice flux and mass balance estimates, and ice shelf cavity and ice sheet models. Bamber and Bentley (1994) also found that mismatch between hydrostatic and measured surface heights near the grounding line on the Ross Ice Shelf were associated with high densities of ice draining fast-flowing East Antarctic glaciers. Furthermore, disagreements between observed and hydrostatic thickness on sub-kilometer scales and near ice fronts (e.g. in regions of accreted marine ice) introduce inaccuracies in estimates of basal mass change (Bamber and Bentley, 1994; Griggs and Bamber, 2011; Chuter and Bamber, 2015) and complicate understanding of the impact of small-scale features like basal channels on ice shelf stability (Drews, 2015).

Over a decade of ice shelf thickness measurements from airborne IPR collected by the NASA Operation IceBridge (OIB) and pre-OIB and NSF Investigating the Cryospheric Evolution of the Central Antarctic Plate (ICECAP) programs provide an extensive dataset with which to examine possible deviations from hydrostatic equilibrium on Antarctic ice shelves (MacGregor and others, 2021). Here, we use airborne laser altimeters (OIB Airborne Topographic Mapper (ATM) and ICECAP Riegl Laser Altimeter (RLA)) to estimate hydrostatic thicknesses, and compare these to measured thicknesses from IPRs that were flown simultaneously (OIB Multichannel Coherent Radar Depth Sounder (MCoRDS) and ICECAP High Capability Radar Sounder (HiCARS)), on sub-kilometer to ice-shelf scales. We also test the sensitivity of the hydrostatic thickness to the use of different firn and mean dynamic topography (MDT) corrections to elucidate the implications of the hydrostatic assumption.

2. Methods

2.1 Study area

Contemporaneous surface and thickness data are binned by discrete ice shelves or ice shelf systems. This results in 20 ice shelf systems, including the Ronne-Filchner Ice Shelf, two Antarctic Peninsula ice shelves, six West Antarctic ice shelves, the Western Ross/McMurdo Ice Shelf and ten East Antarctic ice shelves (Fig. 1). Henceforth, the Ronne-Filchner Ice Shelf, Antarctic Peninsula ice shelves and West Antarctic ice shelves will be collectively referred to as West Antarctica, and the Western Ross/McMurdo and East Antarctic ice shelves will be collectively referred to as East Antarctica.

Figure 1. Map of Antarctica showing the ice shelf system boundaries (boxes) colored by the mean hydrostatic residual for the case in which steady state FDM firn corrections and MDT corrections are applied. Also shown are the IPR ground track coordinates (gray points represent all IPR data; white points are those used in the hydrostatic residual analysis). Base map is the REMA DEM hillshade image, and the black curve shows the 2007–09 InSAR grounding line.

2.2 Estimation of hydrostatic thickness

The hydrostatic ice shelf thickness, H_E, is estimated from ice shelf freeboard height (h) as:

(1)$$\matrix{ {H_E = h\displaystyle{{\rho _s} \over {\rho _s-\rho _i}}-H_a\displaystyle{{\rho _i-\rho _a} \over {\rho _s-\rho _i}}} \cr } $$

where ρ_s is seawater density (1027 kg m⁻³), ρ_i is meteoric ice density (918 kg m⁻³), ρ_a is the firn-air column density (2 kg m⁻³), and H_a is the thickness of the firn-air column within the freeboard (specifically defined as the length of the change in firn thickness resulting from compressing the firn column to ice density (Ligtenberg and others, 2011)), and the subscript E denotes that H_E is an estimate of ice thickness (Fig. 2).

Figure 2. Schematic showing relevant quantities for a column of ice floating in seawater. The ice below sea level is discontinuous to exaggerate the vertical scale. Quantities on the left of the ice column represent observed values for an ice column that is not necessarily in hydrostatic equilibrium. H, total observed ice thickness; h, observed freeboard thickness; H_a, thickness of the firn air column; Quantities on the right represent H_E, total hydrostatic ice thickness for an ice column with the observed freeboard, h, in hydrostatic equilibrium, calculated using Eqn (1); R is equal to the difference between H_E and H.

Surface elevation measurements, which are corrected to freeboard heights were collected by laser altimeters from the OIB and ICECAP programs. For West Antarctica, we use surface elevations from the ATM L1B (i.e. geolocated ice elevation) datasets for both pre-OIB campaigns, which cover 26 November 2002–29 November 2004 (Studinger, 2012), and OIB campaigns, which cover 16 October 2009–16 November 2018 (Studinger, 2013). For East Antarctica, we use surface elevations from the OIB/ICECAP RLA (a LD90-3800-HiP-LR distance meter) L2 (i.e. geolocated ice elevation) dataset for 13 January 2009–21 January 2013 (Blankenship and others, 2012b).

Only OIB and ICECAP campaigns in which ATM/RLA surface elevation data were collected simultaneously with MCoRDS/HiCARS thickness data between 26 November 2002 and 16 November 2018 are used in this study. This enables a direct comparison between hydrostatic thicknesses derived from ATM/RLA surface elevation data and thicknesses measured by the IPRs. As such, ATM and RLA point cloud surface elevation data are interpolated to contemporaneous MCoRDS and HiCARS ground track coordinates using natural neighbor triangulation with no extrapolation. This method is chosen to reduce unconstrained extrapolation on the edges of the point cloud, and because ATM and RLA point clouds are spaced similarly to or more densely than the MCoRDS and HiCARS point clouds. To estimate hydrostatic thickness, surface elevation data (z) must be corrected for tides, MDT and referenced to the geoid to obtain freeboard heights (h), given h = z – tide – geoid – MDT. The Eigen-6C4 geoid (Förste and others, 2014; Morlighem and others, 2020) is bilinearly interpolated to the ground track coordinates, and tide corrections are obtained from the CATS2008b tide model (Padman and others, 2018). However, MDT values do not extend into most ice shelf cavities in Antarctica (Andersen and Knudsen, 2009), and model results must either be extrapolated, or MDT must be accounted for as an uncertainty. We estimate hydrostatic thickness both with and without extrapolated MDT values (setting MDT to zero in the nonextrapolated case). MDT values are obtained from the DTU22 model (Knudsen and others, 2021); these gridded data are bilinearly interpolated to the ground track coordinates where both datasets overlap. Because MDT data doesn't extend into several of the ice shelf cavities, we use nearest neighbor extrapolation to fill in missing values along-track with the nearest interpolated value.

We quality-control the freeboard heights by comparison with the Reference Elevation Model of Antarctica (REMA) 200 m Digital Elevation Model (DEM) mosaic (Howat and others, 2019; Howat and others, 2022). REMA freeboard heights are bilinearly interpolated to the ground track coordinates for each ice shelf, and points at which the absolute value of the difference between ATM/Riegl and REMA freeboard heights falls outside the 95% confidence interval are removed. This window is used to avoid exclusion of airborne observations that differ from REMA due to advection of surface features, but to exclude erroneous observations due to clouds or measurement errors. All freeboard heights less than 20 m in magnitude are removed, as these data likely reflect open ocean or sea ice.

Gridded firn air column thickness values, H_a, were obtained from both a steady-state firn densification model (FDM) and a time-evolving FDM. The steady-state FDM (henceforth termed sFDM) is the Institute for Marine and Atmospheric research Utrecht steady-state FDM (Ligtenberg and others, 2011), forced at the surface by output of the regional climate model RACMO2.3p2 (van Wessem and others, 2018). These H_a values are bilinearly interpolated to the MCoRDS/HiCARS ground track coordinates. This FDM output was selected because it is included in BedMachine Antarctica, Version 2 (Morlighem and others, 2020). The time-evolving FDM (henceforth termed tFDM) is the NASA Goddard Space Flight Center FDM, version 1.2.1 at 25 km resolution (Medley and others, 2022a). We bilinearly interpolate the 6-day firn air content closest in time to each airborne campaign to the IPR ground track coordinates. The data are filtered further by removing all ground track points where firn air content, H_a, exceeds freeboard height, as these may produce negative thicknesses (Griggs and Bamber, 2011), although it is possible for firn to extend below the freeboard height (Cook and others, 2018). Freeboard heights and firn air column thicknesses are then used to obtain hydrostatic thicknesses using Eqn (1).

2.3. Estimation of hydrostatic residual

The hydrostatic residual (R) is defined as the difference between the estimated H_E and the observed H:

(2)$$\matrix{ {R = H_E\;-\;H.} \cr } $$

Thus, a positive value for R indicates that the actual ice thickness is less than hydrostatic thickness, or that the freeboard is elevated relative to sea level (Fig. 2), with the opposite for a negative R. We expect R to depend on measurement errors, assumptions in ice density and firn thickness and ice dynamics, such as transfer of vertical stresses.

Ice thickness measurements, H, were collected simultaneously with ATM and RLA surface elevations by airborne IPRs from the OIB and ICECAP programs. Thus, we may directly compare H_E derived from ATM/RLA data to H (Eqn (2). For West Antarctica, ice thicknesses are obtained from the MCoRDS L2 (i.e. geolocated ice thickness with ice surface and ice bottom elevation) datasets from both pre-OIB campaigns (26 November 2002–29 November 2004; Paden and others, 2011) and OIB campaigns (16 October 2009–16 November 2018; Paden and others, 2010). For East Antarctica, ice thicknesses are obtained from the HiCARS 1 (a 52.5–67.5 MHz instrument with two 12-bit digitizer channels; 13 January 2009–21 December 2010; Blankenship and others, 2011) and HiCARS 2 (a 52.5–67.5 MHz instrument with two 14-bit digitizer channel; 05 December 2010–21 January 2013; Blankenship and others, 2012a) L2 (i.e. geolocated ice thickness with ice surface and ice bottom elevation) datasets from the ICECAP project, which operated from 2008–2013 and included four OIB campaigns (Blankenship and others, 2011, 2012a). All thickness data less than 20 m in magnitude are removed, as these likely reflect open ocean or sea ice.

The hydrostatic residual is calculated using Eqn (2) for each ice shelf system, excluding all data upstream of the MEaSUREs Antarctic Grounding Line from Differential Satellite Radar Interferometry from the 2007–2009 IPY (Rignot and others, 2013; Mouginot and others, 2017).

3. Uncertainties and errors

Uncertainties in MCoRDS ice thicknesses are estimated to be ±50 m (Medley and others, 2014), and HiCARS ice thicknesses have a reported uncertainty of ±70 m (Blankenship and others, 2011). However, since these are nominal values with potentially different values over ice shelves, we perform a crossover analysis to assess the self-consistency of the data. Crossover points are located by splitting the ground tracks into 5000-point segments for each ice shelf and finding the intersections of all possible segment combinations. This method identifies intersections not only where the ground tracks cross one another at large angles, but also where repeated ground tracks overlap. Where repeated ground tracks overlap, the intersections are frequently only meters apart. We thus ignore repeat-track intersections that are within 1 km to ensure that redundant points are excluded. All measurements that fall within 50 m of an intersection are differenced from each other, showing expected changes in thickness through time. Crossover and repeat-track intersections from the same MCoRDS campaign have a mean absolute difference in H of 3.3 m and a standard deviation of H of 2.6 m (Fig. S1). For HiCARS, these are 3.6 and 4.1 m, respectively. These estimates, however, provide the measurement precision of the instruments, and do not account for biases due to firn penetration or radar attenuation, so we adopt an uncertainty of ±50 m for our propagation of errors based on the literature. Even with this large uncertainty, MCoRDS and HiCARS provide the best, large-scale ice shelf thickness measurements available. We do not consider errors or uncertainties in ρ_s, which varies by <1 kg m⁻³ in the top 1 km of the ocean (Jackett and McDougall, 1997), or ρ_i, which has accepted values ranging from 910–921 kg m⁻³ and often varies by less than ± 5 kg m⁻³ (e.g. Griggs and Bamber, 2011), choosing to keep these values constant throughout our analyses. Errors for other data sets and calculations are reported in Table 1.

Table 1. Errors/Uncertainties for data involved in the calculation of R

We propagate the above errors (σ) and uncertainties in Eqns (1) and (2) as:

(3)$${ \matrix{ {\sigma _R = \sqrt {{( {c_1\sigma_z} ) }^2 + {( {c_1\sigma_{tide}} ) }^2 + {( {c_1\sigma_{geoid}} ) }^2 + {( {c_1\sigma_{MDT}} ) }^2 + {( {c_2\sigma_{Ha}} ) }^2 + \sigma _H^2 } , \;\;\;} \cr } }$$

where c ₁ = ρ_s / (ρ_s – ρ_i) and c ₂ = (ρ_i – ρ_a) / (ρ_s – ρ_i). This gives a combined error of ±84 m for H_E and ±98 m for R with the sFDM, and ±42 m for H_E and ±65 m for R with the tFDM. It's unclear, however, how much this error varies spatially due to its dependence on the firn correction because both the sFDM and tFDM are posted at much lower resolution than the airborne data.

4. Results

4.1. Hydrostatic residual

We perform analyses of R for individual ground tracks, sectors of ice shelf regions, whole ice shelves, and for the complete dataset. To remove outliers, and because R is dependent on H, we bin R by 100 m intervals of H at the corresponding ground track coordinates. We then remove measurements where R < Q1–1.5 × IQR or R > Q3 + 1.5 × IQR, where Q1 and Q3 are the 25 and 75% percentiles of the binned R values, respectively, and IQR = Q3 – Q1 is the interquartile range (Lane and others, 2013).

4.1.1 Comparison of hydrostatic residual between corrections scenarios

We test the impact of applying different H_a and MDT corrections by considering six cases: using a steady state FDM for H_a corrections, using a transient FDM for H_a corrections, and no H_a correction applied, each with and without MDT corrections applied. We find vastly differing mean R values among the cases with different H_a corrections, although for each FDM case, the inclusion of an MDT correction (which ranges from −1.3 to −1.1 around Antarctica, Fig. S3) results in a more positive mean R (Table 2, Fig. 3). This effect holds for individual ice shelves as well as the aggregate results (Table S1). In general, the sFDM produces hydrostatic residuals closest to 0 (mean R = 17 ± a standard deviation of 51 m with MDT applied; mean R = 6 ± 51 m without MDT applied). Notably, the ice shelves into which the MDT model does not extend, and for which MDT was extrapolated (Ronne Filchner, George VI/Wilkins/Stange, Western Ross/McMurdo) have R values closest to 0 with no MDT correction applied (Table S1). The tFDM reports thicker H_a throughout the study period than the sFDM (Text S1.1, Fig. S2), resulting in thinner H_E and more negative R (mean R = −28 ± 64 m with MDT applied, mean R = −39 ± 64 m without MDT applied). As expected, the case with no H_a correction applied results in the most positive and largest magnitude R, because H_E is computed as though the entire ice shelf column is pure ice (ρ _i = 918 kg m⁻³), resulting in a mean R = 122 ± 63 m with MDT applied and 111 ± 64 m without MDT applied. The large standard deviations in R for each case indicate that there is significant spatial variability that is not accounted for by the H_a corrections. Unless otherwise noted, we henceforth report results only for the case with sFDM H_a and MDT corrections applied to further investigate spatial variability, since the different corrections contribute very little to spatial variability due to the models' relatively coarse resolutions.

Table 2. Aggregate hydrostatic residual results for the six cases with different H_a and MDT corrections applied (σ =standard deviation). All units are meters

Figure 3. Histograms of R for each corrections scenario.

4.1.2 Hydrostatic residuals among ice shelves

For the nine West Antarctic ice shelf systems with contemporaneous ATM and MCoRDS measurements, the mean R is ~16 m (7% of measured ice thickness). In other words, the observed ice thickness is 16 m less, on average, than the hydrostatic thickness estimated from freeboard. The hydrostatic residual varies significantly between individual ice shelves both in absolute and relative magnitudes (Table 3). Two ice shelves systems, Dotson/Crosson and Nickerson, have negative mean and median hydrostatic residuals, but these are within 1% of the observed ice shelf thicknesses. Larsen Ice Shelf has both the greatest absolute mean R of 27 m, and the greatest mean percent overestimation at 13%. Although the magnitudes of R for the Abbot and Getz ice shelves are similar to one another, the mean percent overestimation is 12% for the thinner Abbot ice shelf and only 2% for the thicker Getz ice shelf. On average, observed ice thicknesses of West Antarctic ice shelves are 7% less than the hydrostatic thickness predicted from freeboard.

Table 3. Overview of hydrostatic residual (R) and related statistics for all ice shelves in the case with sFDM and MDT corrections applied (σ = standard deviation).

For the 12 East Antarctic ice shelf systems with contemporaneous RLA and HiCARS measurements, the mean R is 26 m, or the observed thickness is 4% thinner than the estimated hydrostatic thickness. The density of observations is much lower in East Antarctica than West Antarctica, and several ice shelves have coverage by only one or two campaigns. Thus, it is difficult to generalize results for most ice shelves, and R values vary more widely than on West Antarctic ice shelves (Table 3). Measurements of the Ross Ice Shelf, while dense, only cover its far western portion near McMurdo and the northernmost glaciers draining the Trans-Antarctic Mountains, so we term this region the Western Ross/McMurdo Ice Shelf system. This ice shelf system and Shackleton Ice Shelf, which has less dense coverage, have R values closest to zero, at 12 and 11 m respectively, both corresponding to a 3% thickness overestimation. Of the East Antarctic ice shelves surveyed by three or more campaigns, the Moscow University and Totten Glacier ice shelves (which were surveyed with similar density as West Antarctic ice shelves) show the largest disagreements in absolute magnitude between hydrostatic and measured thickness, but with low fractional overestimations of ice thickness (3 and 4%). Only the Vincennes Bay/Underwood Ice Shelf system exhibits negative mean R values, but the percent error indicates an overestimation of 11%, indicating that there are negative outliers skewing the mean. The other ice shelf regions with fewer than three campaigns (Cook, Ninnis, Frost/Holmes and West) exhibit positive hydrostatic residuals. The disparities between densely and sparsely surveyed ice shelves indicate that there is high spatial variability within an ice shelf as well as among ice shelves.

Overall, over three-quarters of point estimates of H_E are within 10% of H, and over 60% are within 5%, or 25 m. Histograms of both R and the percent misestimation have a positive skew, although the mode of R values for West Antarctica is positive, while the mode for East Antarctica is negative (Fig. 4). Notably, the mode for the percent difference is negative (between −2 and 0%) for West Antarctica, although it is positive (between 0 and 2%) for East Antarctica and the mean and median for both ice sheets are positive (Figs 1, 4).

Figure 4. a, b: Histograms of R; c, d: histograms of percent difference between hydrostatic and measured ice thickness for all MCoRDS (West Antarctica, a, c) and HiCARS (East Antarctica, b, d) data used in analysis.

4.2. Spatial Variability

4.2.1. Variability on >10 kilometer scales

Although a few patterns emerged, observed and hydrostatic thicknesses vary widely within and among ice shelves. We sample several ground tracks in each ice shelf system to investigate spatial patterns in R (see Supplement). In general, over distances >10 km, the hydrostatic assumption overestimates ice thickness where the ice is relatively thicker, and underestimates ice thickness where the ice is relatively thinner. Alternatively, the freeboard is elevated relative to the predicted flotation level where the ice is thicker and is depressed relative to flotation where the ice is thinner. Furthermore, R tends to increase with distance from the grounding line. Figure 5a shows that after binning ground track coordinates into 25-km increments of the shortest Euclidean distance to the grounding line, both the mean R of points within each bin and the mean R of all points included in the current and all previous bins (cumulative R) increases with distance from the grounding line. Specifically, the mean R in the 0–25 km bin is 14 m and the mean R of all points within 200 km of the grounding line is 17 m. However, only the Ronne-Filchner and Larsen ice shelves have data >200 km from the grounding line (Fig. 1). Notably, we find the opposite pattern within 10 km of the grounding zone, discussed further in the next section.

Figure 5. Left Y axis shows the cumulative (light gray) and bin total (dark gray) number of points within each successive distance from the grounding line (0 km). Right Y axis shows the mean R of all cumulative points (solid curve) and points within each bin (dashed curve) for each successive distance from the grounding line. Panel a shows bins of 25 km; b shows bins of 1 km.

4.2.2. Grounding zones

Figure 5b shows that when we bin ground track coordinates by 1 km increments of the shortest Euclidean distance to the grounding line, both the mean R of points within each bin and the mean R of all points included in the current and all previous bins (cumulative R) decreases with distance from the grounding line. Specifically, the mean R is 36 m for points within 1 km of the grounding line and 15.0 m for all points within 10 km of the grounding line (Fig. 5b), which is close to the mean R value for points within 25 km of the grounding line (Fig. 5a). On several ice shelves, the characteristic surface break-in-slope (Fricker and others, 2009) within 10 km of grounding lines is associated with a highly variable R along-track (Fig. 6). IPR ice thicknesses are generally less than hydrostatic near the grounding line and greater than hydrostatic at the local surface minima (Fig. 6a at 5 km and 122 km, Fig. 6b at ~1 km, Fig. 6c at ~0.5 km) or inflection point (Fig. 6c at 15 km) at the break-in-slope, which is often associated with a local thickness maximum. Where the surface height rebounds further along-track, observed thicknesses drop back below hydrostatic. Beyond the grounding line break-in-slope feature, however, variations in R are not necessarily similar along these ground tracks. For all sampled ground tracks (Figs S4–S45) that intersect the grounding line, 72% show negative values for R coinciding with the break-in-slope.

Figure 6. Selected transects that start and/or end at a grounding line with a break-in-slope feature 1–5 km from the grounding line. Top subpanel of a–c shows freeboard height h (blue curve, left Y axis), IPR thickness H and hydrostatic thickness H_E (orange solid and red dashed curves, right Y axis), while the bottom subpanel shows hydrostatic residual R. Map insets show the location of each transect (a: transect b-b’ downstream of Institute Ice Stream, b: transect f-f’ on Thwaites Ice Shelf and c: transect b-b’ on Cook Ice Shelf), with plotted portions in a-c marked in orange.

4.2.3. Variability on kilometer scales

A pattern common to all ice shelves is that changes in R are generally inversely related to changes in H over distances <10 km, with some exceptions. This indicates that the surface topography is muted relative to the thickness profile, especially where peaks in observed thickness and freeboard height are associated with negative R values, and where local minima in observed thickness and freeboard height are associated with positive R values. However, sampled ground tracks also show that surface peaks and troughs aren't always aligned with variations in the thickness profile, and that there are some regions where the surface topography is exaggerated compared to the thickness profile. Figure 7 shows examples of these patterns along transects from the Foundation ice stream sector of the Ronne-Filchner Ice Shelf, the Getz Ice Shelf and the Totten Ice Shelf.

Figure 7. Selected flow-transverse transects with shading to highlight the relationship between H and H_E ~ different topographic features. Yellow (green) shading highlights where the surface topography is muted (exaggerated) compared to the thickness profile, and blue (red) shading highlights where the freeboard is too low (high) within large surface troughs/thin points (such as basal channels). Top subpanel of a-c shows freeboard height h (blue curve, left Y axis), IPR thickness H and hydrostatic thickness H_E (orange solid and red dashed curves, right Y axis), while the bottom subpanel shows hydrostatic residual R (black curve, left Y axis) and the sum of normal strain rates and the shear strain rates (solid blue and dashed red curves, right Y axis). Map insets show the location of each transect; a: MCoRDS transect b-b’ on Ronne-Filchner ice shelf in the Foundation ice stream sector, b: MCoRDS transect e-e’ on the Getz Ice Shelf, and c: HiCARS transect e-e’ on Totten Ice Shelf.

Two basal channels are intersected by the Foundation sector transect b-b’, at 5–10 km and 20–24 km (Fig. 7a). Both basal channels exhibit a mismatch between surface slope and thickness gradient, leading to thinner ice than hydrostatic on the true right flank and thicker ice on the true left flank, with the mean also thicker than hydrostatic. This pattern is common to other basal channel intersections, such as those intersecting the Getz transect e-e’ at 21–27 km and at 75–80 km (Fig. 7b) and the Totten transect e-e’ at 7–11 km (Fig. 7c). However, some basal channels are thicker than hydrostatic at both flanks and thinner than hydrostatic within the channel, particularly when the surface trough and thickness minimum are aligned, such as those that intersect the Getz transect e-e’ at 42–47 km (Fig. 7b) and the Totten transect e-e’ at 15–21 km and 23–27 km (Fig. 7c). Similar patterns can also be seen in other selected transects, shown and described in the Supplementary Material.

5. Discussion

5.1. Spatial variability in R

The spatial variability shown in our estimates of R is somewhat consistent with other studies, particularly near the grounding zone. Near the grounding zone, we do not expect the ice to be freely floating because it is dynamically linked to the grounded ice and experiences flexure due to variations in sea level (e.g. tides) for several kilometers downstream of the true grounding line rather than simple vertical displacement (Rignot and others, 2011a; Friedl and others, 2020). The grounding line used in this study was identified from differential satellite radar interferometry data acquired in 2007–09, and thus most closely represents the landward limit of tidal flexure (Rignot and others, 2013; Mouginot and others, 2017); much of the airborne thickness and altimetry data included in analysis are likely within the flexure zone, which often extends a few hundred meters to a few kilometers past the break-in-slope or surface minimum (Rignot and others, 2011a). The distance between the grounding line and the first seaward point at which the ice is freely floating depends on ice rheology, surface and basal topography, ice velocities and the thermal forcing of the ocean (Griggs and Bamber, 2011). Changes in ice properties may lead to decoupling between thickness and surface height gradients (Rignot and Jacobs, 2002), leading to high hydrostatic residuals. Griggs and Bamber (2011) showed that IPR thickness measurements were up to 100 m thinner than those obtained from ERS-1 surface heights within 10 km of the grounding line, which was attributed to poor data coverage due to loss of lock by ERS-1 in regions with steep topography (which led to interpolation errors) and/or the breakdown of hydrostatic equilibrium near the grounding line. In contrast, Chuter and Bamber (2015) found the opposite sign in hydrostatic residual near the grounding line, which was attributed to greater data density from CryoSat-2 compared to ERS-1 and ICESat, which reduced interpolation errors and resulted in thinner hydrostatic thicknesses. Both studies found greater absolute hydrostatic residuals and standard deviations near grounding lines than over entire ice shelves, attributed to the breakdown of the hydrostatic assumption near the grounding zone due to vertical stresses associated with elastic bending and to greater uncertainties in firn thickness on the steep slopes within the grounding zone. Our results are more consistent with those of Griggs and Bamber (2011), as the mean R within 10 km of the grounding line is consistently positive (Fig. 5b), although we do find negative R values associated with the break-in-slope of the surface profile within 10 km of the grounding line (Fig. 6). Our more detailed observations show that the ice is possibly freely floating at 6–8 km from the grounding line (Fig. 6), and we concur that the hydrostatic assumption is unreliable within this distance.

Hydrostatic residuals may reflect uncertainties in the parameters used to calculate hydrostatic thickness and/or physical phenomena preventing the ice from floating freely. The flotation of an ice shelf is dependent on its geometry and velocity; stress transfer may bend the ice to be concave or convex, thus raising or lowering the ice column. Furthermore, estimates of hydrostatic thickness rely on the modeled firn air content, H_a, which is highly uncertain, as firn thickness can vary on sub-km scales not captured in FDMs (Medley and others, 2022b). Underestimation of the firn density or thickness would result in an overestimation of hydrostatic ice thickness based on its freeboard, and vice versa. Indeed, the sFDM reported lower H_a values than the tFDM, resulting in more positive hydrostatic residuals (Table 2). Below, we discuss the measurement errors and uncertainties that may contribute to hydrostatic residuals, and we assess their impacts on basal melt rate estimates.

5.2. Confidence in ice penetrating radar thickness measurements

Our crossover analysis shows that radar thickness measurements were highly self-consistent. This indicates that hydrostatic residuals cannot be explained by lack of precision in thickness measurements, but it does not rule out the possibility that the MCoRDS or HiCARS thickness measurements are biased due to radar attenuation. Indeed, HiCARS ice thicknesses are reported to tend to be biased high based on a first return, and biased low based on a nadir return (Blankenship and others, 2011). Outliers likely represent steep thickness gradients near the intersections due to crevassing or other damage to the ice. Furthermore, shear heating in ice sheet shear margins has been associated with radar signal attenuation leading to dimmed basal echoes and absent or low-confidence radar picks (Summers and Schroeder, 2022), however our data show no clear relationship between missing or low-confidence radar picks and high shear strain rates.

5.3. Impact of ice column component thickness and density on hydrostatic imbalance

We do not assess the impact of accreted marine ice on the hydrostatic residual for the ice shelves in this study. Marine ice can have a density of up to 938 kg m⁻³ (Craven and others, 2009), so we expect that failure to consider accreted marine ice would lead to an underestimation of hydrostatic thickness since a denser ice column sits lower in the water column. Griggs and Bamber (2011) found that ice thickness was underestimated by 5% by not including marine ice (thereby underestimating ice density) in the upper-bound case where half of the total thickness is composed of marine ice. The presence of marine ice may also result in low-confidence picks for the ice shelf base due to its higher conductivity and radar wave energy absorption than meteoric ice (Vaňková and others, 2021). The thickness of marine ice has been estimated for several ice shelves, but few of these areas were surveyed in our dataset. On the Ronne-Filchner ice shelf, marine ice exceeding 100 m in thickness is expected north of 80° S (Vaňková and others, 2021), but most of our ground tracks fall south of this latitude. Marine ice up to 80 m thick is also expected along several flowlines on Larsen C ice shelf (Holland and others, 2009; Harrison and others, 2022), but these regions are not associated with anomalous R values (Fig. S14).

Uncertainty in the thickness and density of firn may contribute to hydrostatic residuals. We approximate how these parameters would need to change for the measured H and h to satisfy the hydrostatic assumption for the cases in which the sFDM or tFDM H_a and MDT corrections are applied. When referring to the firn air column thickness and density necessary to satisfy the hydrostatic assumption, we will denote them with the subscript E for consistency with H_E.

Because R is generally positive, the ice must be less dense than assumed for the measured thickness and freeboard to be in hydrostatic equilibrium. This disparity in densities could be a result of uncertainties in the modeled H_a and/or assumed density ρ_a. To independently investigate the thickness of the firn air column needed to account for R (H_aE), we substitute H and H_aE for H_E and H_a in Eqn (1), leaving ρ_a constant, and set the difference between H_E and H equal to R, so that:

(4a)$$\matrix{ {R = ( {H_{aE}-H_a} ) \displaystyle{{\rho _i-\rho _a} \over {\rho _s-\rho _i}}.} \cr } $$

To independently investigate the firn air column density (ρ_aE), needed to account for R, we substitute H and ρ_aE for H_E and ρ_a in Eqn (1), leaving H_a constant, and again take the difference between the equations for H_E and H:

(4b)$$\matrix{ {R = ( {\rho_a-\rho_{aE}} ) \displaystyle{{H_a} \over {\rho _s-\rho _i}}.} \cr } $$

We can then directly solve for H_aE – H_a (assuming ρ_a = 2 kg m⁻³) or ρ_aE – ρ_a (assuming modeled H_a), eliminating the need to explicitly calculate H_aE and ρ_aE.

Equation (4a) shows that when R is positive, the equilibrium H_aE must be proportionally greater than H_a, positive so that air with a density of 2 kg m⁻³ accounts for more of the total thickness of the ice shelf, decreasing the vertically averaged column density to flotation, and vice versa. A thicker firn-air column would account for the higher observed h required for the observed H to satisfy the hydrostatic assumption, because it would lower the density of the observed ice column, forcing it to float higher in the water (i.e. higher freeboard, smaller submerged portion than if the ice column were denser; Fig. 8). In reality, a thicker firn air column, as seen in the tFDM, indicates a deeper firn layer (Ligtenberg and others, 2011). Similarly, if we assume that the modeled firn-air column thickness is correct but that the density is unknown, Eqn (4b) shows that when R is positive, ρ_aE – ρ_a must be negative in order to bring the vertically averaged column density down, and vice versa.

Figure 8. Cartoon graphic showing relevant quantities for a column of ice floating in seawater. The ice below sea level is discontinuous to exaggerate the vertical scale. Quantities represent observed or accepted values as in Fig. 2, with added H_aE, which is the firn air column thickness necessary to bring the observed ice column into hydrostatic equilibrium, and dH_a, which is the difference between H_aE and the modeled firn air column thickness H_a.

Overall, for the case with sFDM H_a and MDT corrections applied, the mean H_aE – H_a is 2 m, and for the case with tFDM Ha and MDT corrections applied, the mean H_aE – H_a is −4 m (Table S3). Thus, the sFDM H_a more closely match the H_aE required for hydrostatic equilibrium than the tFDM H_a. Both mean values are within the nominal uncertainties of both firn models (Table 1), but this uncertainty is poorly spatially constrained. Indeed, a change in H_a of ±10 m would result in an R of ±84 m, and our R values exceed± 84 m in several places (Figs 6, 7), even resulting in negative H_aE over short distances (e.g. Fig. 9). Furthermore, the direct relationship between H_aE – H_a and R means that the firn-air column thickness would vary widely over the same spatial scales as the hydrostatic residual. Although spatial variability in H_a is driven primarily by surface climatic conditions, which have not been modeled on sub-km scales (Ligtenberg and others, 2011; Lenaerts and others, 2014; Ligtenberg and others, 2014), more recent studies have shown that surface accumulation can vary on km scales (Dattler and others, 2019). Our results show that regions like Remnant Larsen B, and near the Bawden ice rise on Larsen C require a > 10 m thicker firn air layer than modeled to satisfy the hydrostatic assumption, despite absent or near zero modeled and observed firn thicknesses in this region (Holland and others, 2011; Ligtenberg and others, 2011).

Figure 9. Thwaites transect d-d’ showing modeled H_a (black curve), and H_aE (gray curve).

The mean change in ρ_a would result in unphysical mean hydrostatic firn air column densities (ρ_aE) for all but three ice shelves in the case with sFDM H_a and MDT corrections applied, indicating that uncertainties in accounting for the firn air column alone cannot explain R (Table S3). However, the case with tFDM and MDT corrections applied resulted in positive ρ_aE for all but four ice shelves. This disparity points to the need for more observations of firn properties and firn densification models of higher confidence. Larsen C has the most negative ρ_aE – ρ_a (ρ_aE = −1633 kg m⁻³) required for balance in the sFDM case and the second most negative in the tFDM case (ρ_aE = −71 kg m⁻³), providing further evidence that the measured freeboard is much too high for ice with the observed thickness and H_a from either FDM to be in hydrostatic equilibrium.

5.4. Relationship between R and strain rates

If hydrostatic balance may partly be due to the transfer of vertical stress (i.e. stress bridging) within the ice shelf, we expect that R will also be related to strain rates (Cuffey and Patterson, 2010). We estimate longitudinal (e_lon), transverse (e_trans) and shear (e_shear) surface strain rates from the NASA MEaSUREs InSAR-derived average velocity map (Rignot and others, 2011b, 2017; Mouginot and others, 2012) following the approach of Bindschadler and others (1996) at each measurement point. The relationships between R (from the case with sFDM H_a and MDT corrections applied) and the median ∇·u = e_lon + e_trans (normal strain rates, where ∇ is the del operator and u is velocity) and absolute value of e_shear within 1 m increments of R, are plotted in Fig. 10. We find that, as expected, low-magnitude R values coincide with low strain rates, and the magnitude of R increases with increasingly positive shear and normal strain rates. A negative R means the ice is thicker than hydrostatic (the freeboard is below flotation), which is consistent with increased vertical stress due to bridging (Le Brocq and others, 2013; Drews, 2015). Normal strain rates increase with both positive and negative R, which may depend on the direction of stress transfer.

Figure 10. a: Median normal strain rates (e _lon + e _shear, black dots) and absolute values of shear strain rates (|e _shear|, gray '+' signs) for points within 1 m bins of R for all IPR points. b, c: Median e _lon + e _trans and |e _shear|, respectively, within 1 m bins of R for West Antarctica (blue dots, + signs), East Antarctica (all shelves, red dots, + signs) and East Antarctica excluding the Western Ross/McMurdo ice shelf system (orange dots, + signs). Bins containing fewer than the 40th percentile of N (1100 points for Panel a) are excluded.

The relationship between R and ∇·u described above (Fig. 10a) is dominated by West Antarctic ice shelves, which have greater data density (Fig. 10b). However, individual ice shelves show significant variability (Text S3.1, Figs S46–S47).

For East Antarctic ice shelves, excluding the western Ross Ice Shelf/McMurdo Ice Shelf system, the median ∇·u is near zero for all values of R, with higher magnitudes of R generally correlating with decreasing ∇·u (Fig. 10b). Overall, we find that smaller hydrostatic imbalances tend to be associated with compression, while larger imbalances, particularly positive R values, are more likely to be associated with extension, when the Western Ross/McMurdo Ice Shelf is excluded. However, this general pattern is not consistent at the scale of individual ice shelves (Figs S46–S47).

Higher shear strain rates are associated with increasing magnitudes of R, and this effect is larger for negative R values for both West and East Antarctica (Fig. 10a, c). We would expect more negative R values in areas of higher shear as shear stresses may be transferred horizontally from the interior of the ice sheet to the margin.

5.5. Impact of R on estimates of basal melt/accretion rates

The over/underestimation of the rate of basal mass change is dependent on the signs of R and the strain rates. Because R and median strain rates for the vast majority of points are near zero (Figs 4,10), we expect that the rate of basal mass change estimated from hydrostatic thickness, M_bE, won't be greatly misestimated. Assuming incompressibility of ice, and following the continuity approach, the basal mass balance is estimated as (e.g. Dutrieux and others, 2013; Berger and others, 2017; Shean and others, 2019; Chartrand and Howat, 2020):

(5)$$\matrix{ {M_b = \left({\displaystyle{{DH} \over {Dt}} + H( {\nabla \cdot u} ) } \right)\displaystyle{{\rho _i} \over {\rho _w}}\;-\;M_s, \;\;} \cr } $$

where M_b is the rate of basal mass loss/gain in m w.e. a⁻¹ (meters of fresh water equivalent per year) and is positive for refreezing and negative for basal melt, M_s is the surface ablation/accumulation rate, which is positive for mass gain, ∇ is the del operator, u is the column-average horizontal velocity of the ice (m a⁻¹), and ρ_w is the density of fresh water, 1000 kg m⁻². The density of ice is assumed to be 918 kg m⁻³. Estimates of basal mass balance from spaceborne freeboard height measurements, such as those from Adusumilli and others (2020) rely on the calculation of H_E, and we will thus refer to these estimates as M_bE.

To explicitly calculate the difference in the rate of basal mass loss/gain estimated from H and H_E, termed R_Mb, we substitute M_bE for M_b and H_E for H in Eqn (5) and subtract M_b from M_bE, assuming that DH/Dt = DH_E/Dt, so that these values and M_s cancel out:

(6)$$\matrix{ {R_{Mb} = M_{bE}-M_b = [ {( {H_E-H} ) ( {\nabla \cdot u} ) } ] \displaystyle{{\rho _i} \over {\rho _w}} = R( {\nabla \cdot u} ) \displaystyle{{\rho _i} \over {\rho _w}}.} \cr } $$

Thus, R_Mb balances the extension or compression of the ice and the hydrostatic residual. We estimate ∇·u as described in Section 5.3. We then compare our results from Eqn (6) with basal mass balance rates obtained from the ICESat and ICESat-2 satellite record (Adusumilli and others, 2020), termed M_bE. Where R is positive (thickness is overestimated) and strain rates are tensile, R_Mb is greater than 0, indicating that M_bE is too positive, and where strain rates are compressive, R_Mb < 0 and M_bE is too negative. Where R is negative (thickness is underestimated), R_Mb < 0 and M_bE is too negative where strain rates are tensile, and R_Mb > 0 and M_bE is too positive where strain rates are compressive. These interpretations are also summarized in Box 1. These relationships hold at each ground track coordinate, but not necessarily for the aggregated ice shelf results (Table S4). We divide the absolute value of R_Mb by the absolute value of M_bE – R_Mb (where M_bE is from the satellite record, bilinearly interpolated to ground track coordinates) and multiply by 100 to obtain a percent error of mass balance estimates (Table S4).

Box 1. Impact of R on basal mass balance estimates

Overall, accounting for R using sFDM H_a and MDT corrections results in a mean of 71% and a median of 3% error in the rate of basal mass change calculated in Adusumilli and others (2020). Since strain rates of ice shelves tend to be on the order of 10⁻³ per year, the median percent error aligns with our expectation that the impact of R on M_b is generally small. We expect that the large mean percent error results from division by very small magnitudes of |M_bE – R_Mb|. However, hydrostatic imbalance may introduce a bias that, when integrated over large areas, may be significant to the total mass balance. Also, the impact may be significant in areas of high strain rates, such as at shear margins, or in areas of high R, such as basal channels.

Overall, the mean and median |R_Mb| for all data points is 0.4 and 0.0 ± 2.1 m w.e. a⁻¹, meaning that on average, the hydrostatic assumption does not dramatically over- or underestimate basal melt rates, but the standard deviation of 2.1 m w.e. a⁻¹ indicates that there is spatial variability, depending on the flow regime. However, the impact of R on basal mass change rate estimates varies between ice shelves, and on local scales (Table S4). The most extreme impacts of R_Mb occur on Thwaites Ice Shelf (mean and median R_Mb = 0.3 and 0.1 ± 6.6 m w.e. a⁻¹), the Ninnis ice shelf (−0.9 and −0.3 ± 2.8 m w.e. a⁻¹), and the Vincennes Bay/Underwood Ice Shelf region (−0.8 and 0.0 ± 12.7 m w.e. a⁻¹), and the Shackleton Ice Shelf (−0.3 and −0.0 ± 2.5 m w.e. a⁻¹). When compared to the melt rates from Adusumilli and others (2020), however, the most extreme relative impacts on basal mass balance were on Thwaites Ice Shelf, where M_bE is misestimated by a median of 11%, Cook Ice Shelf, (9%), Shackleton Ice Shelf (10%) and West Ice Shelf (16%). The R_Mb values of the latter three ice shelves are likely dominated by extreme R values due to the relatively low number of ground tracks in those regions.

Our results showing high magnitudes of R near basal channels and other potentially destabilizing features are consistent with other observations (e.g. Drews, 2015; Chartrand and Howat, 2020; Dow and others, 2021) and point to the need for more detailed measurements near these features to accurately account for them in mass balance estimates. Hydrostatic imbalance has been shown to change over time as ice advects over an actively incising basal channel (Chartrand and Howat, 2020), indicating that repeat freeboard height measurements may yield erroneous basal melt rates. Although temporal analysis of R is not a goal of this study, several ground tracks with repeat coverage show that R changes over time at a variety of ice shelf features (Figs S9, S10, S16, S17, S27–29, S41). Furthermore, analyses on the Roi Baudoin and Nansen Ice shelves have shown that satellite-derived surface velocities and related strain rates may be better suited than the hydrostatic assumption to characterize basal feature morphology (Drews, 2015; Dow and others, 2021). However, these studies used near-contemporaneous surface velocities to test the agreement between strain-rate and IPR-derived morphology, which are not widely available for supplementing hydrostatic calculations of ice thickness prior to the epoch of widespread availability of high-resolution speed and surface elevation data, such as from the GO_LIVE/ITS_LIVE (Fahnestock and others, 2016; Gardner and others, 2018) and REMA projects (Noh and Howat, 2019).

Similarly, short-term and short-spatial-scale freeboard changes are largely unrelated to basal mass balance and, if not accounted for, can lead to magnification of errors in estimating changes in ice thickness (Vaňková and Nicholls, 2022). Our results corroborate the assertion that errors in basal melt rates derived from satellite data (e.g. Adusumilli and others, 2020) are not spatially uniform (Vaňková and Nicholls, 2022), because R is not uniform in time or space, imparting unknown and potentially large errors in basal melt rates estimated from freeboard.

6. Conclusions

We completed the first, large-scale comparison between thickness observed from ice-penetrating radar and the hydrostatic thickness estimated from contemporaneous surface elevation measurements over Antarctic ice shelves. Using MCoRDS/HiCARS IPR and ATM/RLA laser altimetry, we have found that Antarctic ice shelves are, on average, about 17 m (6%) thinner than hydrostatic thickness estimated using a steady state FDM for firn air content correction and with MDT corrections applied. However, the mean hydrostatic residual, or the difference between estimated and observed thickness, R, varies among individual ice shelf systems and can vary by hundreds of meters over sub-kilometer scales, regardless of the choice of corrections. The greatest hydrostatic residuals in West Antarctica are found on the Larsen C ice shelf, where the measured thickness is ~27 m, or 13%, less than hydrostatic. Of the East Antarctic ice shelves with similar data density to West Antarctica, the greatest residuals are found on Moscow University and Totten Ice Shelves (R = 34, or 3% and R = 41, or 4%, respectively), although the ice shelves with three or fewer campaigns also have high magnitude R values, reaching 52 m, or 12%, on the Frost/Holmes ice shelf system. We expect that the sparse coverage on these shelves allows extreme values to dominate the mean hydrostatic residual.

On kilometer scales, few spatial patterns in hydrostatic residual are apparent. Most notably, the break-in-slope feature within 10 km of the grounding line is often associated with negative R values, and the mean R decreases (but remains above zero) with increasing distance from the grounding line up to 10 km. Past 25 km, the mean R increases with increasing distance from the grounding line. We also find that hydrostatic thickness sometimes exaggerates thickness anomalies compared to the measurements, and sometimes mutes thickness anomalies, including for surface and basal crevasses and basal channels.

We assess whether measurement errors, uncertainties in firn thickness and/or density could account for the average hydrostatic residuals. A crossover analysis of same-campaign thickness measurements shows high consistency in both MCoRDS and HiCARS data, and low errors are expected for surface elevation measurements. On average, R can largely be corrected by assuming a lower vertically averaged density for ice shelves when sFDM corrections are used. This can be achieved physically by accounting for a negative bias in the modeled sFDM firn air column thickness. However, the variability in R across sub-kilometer scales cannot be explained by measurement errors or assumed firn properties. We posit that higher spatial resolutions and accuracies in firn column observations and densification models are needed for confidence in estimating hydrostatic thickness.

Furthermore, although most R values and strain rates are near zero, higher shear and normal strain rates are associated with |R| > 50 m, which is consistent with the concept of stress bridging where the hydrostatic thickness is less than the measured thickness (i.e. vertical stress transfer may hold the freeboard below its hydrostatic height). However, on small scales, strain rates do not correlate with R. One of the greatest implications of uncertainties in estimating hydrostatic thickness is that it will lead to uncertainties in estimating basal mass balance. Few studies consider thickness gradients across flow when modeling ice shelf flow and mass balance, yet we show that R has substantial implications for flow-transverse ice shelf dynamics, particularly on small scales. By isolating the impact of hydrostatic residual on basal mass balance, we find that overall, the hydrostatic assumption misestimates the rate of mass gain by a median of 3%, but this varies spatially, depending on strain rates and thickness gradients. Furthermore, sampled repeat ground tracks show that R can change over time (in an Eulerian framework), pointing to the need for greater utilization of available thickness data and future thickness measurements, which will in turn improve estimates of hydrostatic thickness over time as well as spatially.