Ice flow models typically assume that ice is incompressible, a reasonable assumption because ice density changes are indeed small and have a correspondingly small effect on the overall mass balance of glaciers and ice sheets. Given the immense volume of the ice sheets, however, even relatively small changes may influence global mean sea level to a degree that severely impacts humanity (Hauer and others, Reference Hauer2020). Here, we quantify the role of gravitational compression and thermal contraction in estimating ice sheet mass.

We describe gravitational compression by considering an ice column in vertical hydrostatic equilibrium, ∂σ_z/∂z = −ρg with vertical Cauchy stress σ_z, gravitational acceleration g = 9.81 m s⁻², ice density ρ and vertical coordinate z. Although the effect of compressibility is most commonly described through the equation of mass conservation (e.g., Lipovsky and Dunham, Reference Lipovsky and Dunham2015), here we take a different analytical approach following Cathles (Reference Cathles2015) that examines the long-term behavior of a linear compressible Maxwell viscoelastic material following, $\hat \sigma _{z} = [ \lambda ( \omega ) + 2 \mu ( \omega ) ] \hat \epsilon _{z}$. Here, hats denote Fourier transforms and λ and μ are complex valued, frequency ω dependent constitutive parameters that correspond to volumetric and non-volumetric deformations, respectively. At times much greater than the Maxwell time (Reeh and others, Reference Reeh, Christensen, Mayer and Olesen2003, ~1 d for ice) that are relevant to long-term glacier evolution, the model describes ice as a compressible fluid with μ(ω) → 0 and λ(ω) → K = 8.9 GPa the bulk modulus (Schulson and Duval, Reference Schulson and Duval2009). These limits have the interpretation that long-term non-volumetric stresses relax but volumetric stresses do not. This explains why compression persists over glacial evolution timescales.

The total surface displacement due to gravitational compression w may then be calculated by combining the constitutive relation with the condition of hydrostatic equilibrium. The result is,

(1)$$w = {\rho_0\, g H^2\over 2 K}$$

with ice thickness H and reference density ρ₀ = 917 kg m⁻³. From Eqn 1 we may confirm that strains are small, e.g. ρ₀gH/K ~ 10⁻³ for 1 km thick ice.

Thermal contraction also changes the density of glacier ice because the vast majority of the interior of the ice sheets is at a temperature well below the freezing point. If the temperature of this ice were raised to the melting point it would undergo thermal expansion equal to αΔT(z) with thermal expansion coefficient α = 5.3 × 10⁻⁵ °C⁻¹ and temperature difference from the pressure melting point ΔT (Schulson and Duval, Reference Schulson and Duval2009). In order to put a lower bound on thermal strains, we assume a linear temperature profile with depth that reaches 0°C at the bed and is equal to the mean annual surface temperature at z = H. Although simplistic, this approach certainly results in a lower bound because many temperature profiles are steeper than linear (i.e. colder) due to the effect of advection and because much of the bed is colder than the melting point temperature (Cuffey and Paterson, Reference Cuffey and Paterson2010, Ch. 9). Using published ice temperature (Ettema and others, Reference Ettema2020; Le Brocq and others, Reference Le Brocq, Payne and Vieli2010) and geometry datasets (Morlighem and others, Reference Morlighem2017, Reference Morlighem2020), we calculate gravitational compression and thermal contraction in the ice sheets.

We calculate that gravitational compression deforms the surface of the Greenland and Antarctic ice sheets by up to 5.8 and 11.3 m, respectively. For thermal contraction, the corresponding values are 2.9 and 6.2 m, respectively. The corresponding ice sheet-averaged values are 0.5 and 0.7 m for thermal contraction and 0.8 and 0.7 m for gravitational compression. Integrated over the entire ice sheet, the total mass bias from the gravitational effect is equal to 3000 Gt in Greenland and 30 200 Gt in Antarctica. Although these long-term biases are only $\sim 0.2\percnt$ of the total ice-sheet mass, they are significantly larger than recent annual ice mass losses (~318 Gt a^-1, Smith and others, Reference Smith2020). Perhaps more important than this long-term bias are potential shorter-term biases.

Estimates of total glacier mass from satellite altimetry commonly account for compression of the solid Earth but not compression of the ice column (Smith and others, Reference Smith2020). To gauge the importance of seasonal ice compression, we calculate the ice column compression caused by surface mass-balance anomalies in Antarctica as described by Medley and others (Reference Medley, Neumann, Zwally and Smith2020). The resulting ice-sheet surface vertical velocities are as large as 1–2 mm a⁻¹ near the coasts and are therefore of comparable or slightly smaller magnitude as solid Earth uplift rates (e.g. Whitehouse and others, Reference Whitehouse, Bentley, Milne, King and Thomas2012; Ivins and others, Reference Ivins2013). Future ice-sheet mass-balance studies using satellite altimetry should therefore take into account the effect of ice-sheet compressibility. Along with improved firn modeling, seasonal compression could plausibly explain some of the greater variability of altimetric versus gravimetric ice-sheet mass estimates (Shepherd and others, Reference Shepherd2018). We note that seasonal compression could be directly measured with optical fiber strain sensing (Hartog, Reference Hartog2017).