Effect of horizontal divergence on estimates of firn-air content

Abstract

Ice-sheet mass-balance estimates derived from repeat satellite-altimetry observations require accurate calculation of spatiotemporal variability in firn-air content (FAC). However, firn-compaction models remain a large source of uncertainty within mass-balance estimates. In this study, we investigate one process that is neglected in FAC estimates derived from firn-compaction models: enhanced layer thinning due to horizontal divergence. We incorporate a layer-thinning scheme into the Community Firn Model. At every time step, firn layers first densify according to a firn-compaction model and then thin further due to an imposed horizontal divergence rate without additional density changes. We find that horizontal divergence on Thwaites (THW) and Pine Island Glaciers can reduce local FAC by up to 41% and 18%, respectively. We also assess the impact of temporal variability of horizontal divergence on FAC. We find a 15% decrease in FAC between 2007 and 2016 due to horizontal divergence at a location that is characteristic of lower THW. This decrease accounts for 16% of the observed surface lowering, whereas climate variability alone causes negligible changes in FAC at this location. Omitting transient horizontal divergence in estimates of FAC leads to an overestimation of ice loss via satellite-altimetry methods in regions of dynamic ice flow.

Introduction

Many outlet glaciers of the polar ice sheets have accelerated and thinned markedly in the last 25 years (Joughin and others, 2012; Mouginot and others, 2014; Smith and others, 2020). While horizontal divergence is low in the ice-sheet interior, outlet glaciers often have substantial spatially and temporally evolving horizontal divergence rates. Investigating many fundamental glaciological problems depends on accurately estimating changes in firn-air content (FAC; the volume of air in a firn column of unit cross-sectional area), including the determination of mass loss from thinning due to marine ice-sheet instability, an important component of the land-ice contribution to sea-level rise (Shepherd and others, 2012; Depoorter and others, 2013; Shepherd and others, 2019; Smith and others, 2020). Despite its importance, most firn-compaction models lack a method to account for horizontal divergence in estimates of FAC in regions of fast flow. Here, we formulate a simple kinematic model scheme that is an accessible and easily applicable alternative to a material-specific constitutive relation (Gagliardini and Meyssonnier, 1997; Lüthi and Funk, 2000). We use this scheme to account for horizontal divergence within the firn to show the importance of the effect of horizontal divergence on FAC estimates in regions of the ice sheet with high and rapidly changing horizontal divergence rates.

Background

Firn compaction occurs through a variety of microphysical mechanisms, such as grain-boundary sliding, sintering and bubble compression. These mechanisms respond to overburden stresses and temperature gradients (Maeno and Ebinuma, 1983; Burr and others, 2019), to processes related to melting and refreezing (Reeh and others, 2005), and to ice-flow stresses (Alley and Bentley, 1988). These processes vary spatially and temporally, resulting in variable firn-density profiles. Model estimates of firn-column thickness, FAC and the rate at which these evolve often have substantial uncertainty, partly because not all physical processes are included or accurately captured within the current generation of firn-compaction models. For example, thinning in the firn column due to horizontal divergence in the underlying ice is often assumed to be negligible even in regions with high horizontal divergence rates (Kuipers Munneke and others, 2015), despite observational evidence of its impact on firn-density structure (Christianson and others, 2014; Morris and others, 2017; Riverman and others, 2019).

Firn-air content

Estimating mass change of the ice sheets with repeat satellite-altimetry observations requires model-derived estimates of changes in FAC (Shepherd and others, 2012; Depoorter and others, 2013; Shepherd and others, 2019). FAC, also known as depth-integrated porosity or DIP, is defined as the porosity integrated over depth z, from the surface to the depth z _i where ice density ρ _i is attained:

(1)$${\rm FAC} = \int_{0}^{z_{\rm i}} \left({\rho_{\rm i} - \rho\lpar z\rpar \over \rho_{\rm i}} \right)\, {\rm d}z.$$

The change in mass of the ice sheet Δm can be calculated in terms of the observed surface height change Δh, the change in FAC ΔFAC, ice density ρ _i and area A:

(2)$$\Delta m = \left(\Delta h - \Delta {\rm FAC} \right)\rho_{{\rm i}}A.$$

Improving firn-compaction models, and specifically producing more accurate estimates of FAC, is an essential step in reducing uncertainty in altimetry-derived mass-balance products. Currently, model estimates of FAC can have large uncertainty partly because firn-compaction models have been calibrated only to, and therefore are appropriate only for, a limited range of climate and ice-dynamic settings (Lundin and others, 2017). Most firn-compaction models also are compatible with the suggestions by Robin (1958) under steady-state conditions. Robin (1958) suggested that, in steady-state conditions, the change in density with depth dρ(z)/dz (and thus FAC) is proportional to the change in overburden stress. This requires that no horizontal divergence occurs within the firn column (Morris and others, 2017):

(3)$${{\rm d} \rho\lpar z\rpar \over {\rm d}z} = \rho \lpar z\rpar \left(\rho_{\rm i} - \rho\lpar z\rpar \right).$$

However, observations collected over the last 60 years indicate that large deviatoric stresses in the underlying solid ice affect the firn by increasing the rate of firn-density change with depth (Zumberge, 1960; Crary and Charles, 1961; Gow, 1968; Kirchner and others, 1979; Alley and Bentley, 1988; Christianson and others, 2014; Vallelonga and others, 2014; Riverman and others, 2019; Morris and others, 2017).

Previous observations and modeling work

Active-source seismic and radar surveys across the Northeast Greenland Ice Stream show that the firn column is 30 m thinner in the shear margins than outside the margins (Christianson and others, 2014; Vallelonga and others, 2014; Riverman and others, 2019). Riverman and others (2019) inferred that this spatial variability of firn density is due to both horizontal divergence and strain softening (i.e., an increased firn-compaction rate due to the acceleration of time-dependent microphysical processes from increased ice-flow stresses). Morris and others (2017) accounted for horizontal divergence in the firn along the iSTAR traverse on Pine Island Glacier by using a layer-thinning scheme similar to the one we propose here (see Section 2 of the Supplementary Material). Morris and others (2017) showed that the negative ratio of the vertical densification rate to the density-corrected volumetric strain rate is not equal to the mean-annual accumulation rate in some cases; they attributed this to a non-negligible horizontal divergence, and illustrated that the steady-state suggestion by Robin (1958) does not hold in those cases.

Some modeling studies have attempted to incorporate the impact of horizontal ice-flow stresses on firn within a constitutive formulation using a generalized form of Glen's Law (Gagliardini and Meyssonnier, 1997; Lüthi and Funk, 2000). This approach has not been widely adopted by the firn community because it is difficult to integrate into commonly used firn-compaction models. To our knowledge, no current firn-compaction model used for inferring ice-sheet mass balance from repeat-altimetry observations (e.g., Ligtenberg and others, 2011; Li and Zwally, 2015, etc.) accounts for horizontal divergence.

In this study, we incorporate a layer-thinning scheme to account for horizontal divergence in the Community Firn Model (CFM) (Stevens and others, 2020), and systematically investigate the effects of horizontal divergence on FAC. We first describe the specifics of our implementation of the layer-thinning scheme in the CFM. Then, we test the scheme through a series of idealized conceptual runs for climate conditions representative of West Antarctica. We subsequently apply the layer-thinning scheme to a firn column along two flowlines on Thwaites Glacier and Pine Island Glacier, Antarctica, where dynamic ice-sheet thinning due to accelerating ice flow from marine ice-sheet instability is occurring (Joughin and others, 2014; Rignot and others, 2014). We then quantify how much of recent observed thinning results from FAC changes due to horizontal divergence for lower Thwaites Glacier. Finally, we map where horizontal divergence should be considered in FAC estimates across the entire Antarctic Ice Sheet

Methods

Horizontal divergence in the Community Firn Model

The CFM is an open-source, modular model framework that is designed to simulate evolution of firn properties, including density, compaction rate and temperature (Stevens and others, 2020). It utilizes a suite of 13 published snow- and firn-compaction models. The CFM uses a one-dimensional Lagrangian framework to track the properties of firn parcels as they advect from the surface into the underlying ice sheet. Users stipulate the surface-boundary conditions, including accumulation rate, surface temperature, surface-snow density and other parameters necessary for the chosen firn- or snow-compaction model.

To simulate layer-thickness changes within the CFM due to horizontal divergence, we adopt a kinematic two-part layer-thinning scheme for firn compaction (Fig. 1). During each time step, the firn first densifies via the equations of the user-specified firn-compaction model:

(4)$$\lambda_{{\rm part}1} = \lambda_{{\rm old}}\lpar 1 + \dot{\epsilon}_{zz}\Delta{t}\rpar.$$

Fig. 1. Our layer-thinning scheme that accounts for horizontal divergence in the CFM. At each time step, the firn first compresses vertically and densifies (Eqn (4)) following the equations of the user-specified firn-compaction model (part one). Then the firn stretches horizontally without further density change, as determined by the prescribed horizontal divergence rate $\dot {\epsilon }_{{\rm h}}$ in Eqn (6) (part two).

where λ_old is the firn-parcel thickness at the previous time step, λ_part1 is the firn-parcel thickness after densification from the firn-compaction model, λ_part2 is the firn-parcel thickness after thinning from horizontal divergence, and Δt is the time step. $\dot {\epsilon }_{zz}$ is the vertical strain rate due to densification of the firn, provided by the chosen model physics, and commonly implemented in models as:

(5)$$\dot{\epsilon}_{zz} = {1\over \rho\lpar z\rpar } {{\rm d}\rho\lpar z\rpar \over {\rm d}t}.$$

Equation (4) (part one) is identical to the procedure in a conventional firn-compaction model that neglects horizontal divergence.

When the firn parcels thin due to horizontal divergence (Eqn (6); part two), during the same time step as densification due to Eqn (4), the firn is horizontally stretched using a prescribed horizontal divergence rate $\dot {\epsilon }_{{\rm h}}$:

(6)$$\lambda_{{\rm part}2} = \lambda_{{\rm part}1}\lpar 1 - \dot{\epsilon}_{{\rm h}}\Delta{t}\rpar.$$

Individual firn parcels stretch horizontally as a result of thinning due to horizontal divergence. Density and FAC calculations of the individual firn parcels in the CFM are produced per unit cross-sectional area of the firn parcels. Eqn (6) does not consider density changes of individual firn parcels. This is because the scheme is solely kinematic, and it is assumed that the material properties of a given parcel of firn do not change with horizontal stretching. However, Eqn (6) reduces the depth at which any density appears in the firn column.

Choice of firn-compaction models

In this study, we run the CFM using the firn-compaction equations from the Herron and Langway (1980) model (HL) and the Ligtenberg and others (2011) model (LIG). Multiple firn-compaction models have been developed in the last 40 years. We use HL because it is still seen as a benchmark firn-compaction model, and most models of polar firn compaction are based on its general framework and assumptions (Li and Zwally, 2011; Ligtenberg and others, 2011; Morris and Wingham, 2014). HL used Antarctic and Greenlandic firn depth-density profiles to derive empirical equations describing the firn-compaction rate in stage one and stage two of the firn column. Stage one is where density ρ < 550 kg m ⁻³, and represents the shallowest portion of the firn column. Stage two occurs deeper and extends to bubble close-off, where densities range from 550 kg m ⁻³ < ρ < 830 kg m ⁻³. Robin (1958) assumed that the change in firn-pore volume is proportional to the change in the overburden pressure, a steady-state assumption. This assumption allows HL to parametrize the overburden pressure using mean annual snow accumulation. In addition, the HL densification rate includes an Arrhenius-type temperature dependence, which represents densification via temperature-dependent microphysical mechanisms, such as grain growth (Gow, 1969).

We also use LIG because it is the firn-compaction model included in subsurface processes of the regional climate model RACMO (Van Wessem and others, 2018), one of the most common reanalysis products used in Antarctic mass-balance calculations. LIG has been used to estimate FAC changes in multiple ice-sheet mass-balance studies (Shepherd and others, 2012; Gardner and others, 2013; McMillan and others, 2016; Shepherd and others, 2019). LIG used 48 depth-density profiles from Antarctica to tune the firn-compaction model by Arthern and others (2010) for general applicability in Antarctica. Arthern and others (2010) measured vertical strain rates in the firn at four sites in Antarctica and used those data to derive a semi-empirical compaction model based on rate equations of Nabarro-Herring creep and grain growth (Coble, 1970; Gow, 1969), with a form similar to that of HL.

Model inputs

Total horizontal divergence rates from ice velocities

We use the Mouginot and others (2019) Antarctic ice-velocity map to compute mean horizontal divergence rates from 1996 to 2018 at a spatial resolution of 450 m. We also use the Mouginot and others (2017) ice-velocity time series to compute the annual horizontal divergence rates from 2007 to 2016 in the Amundsen Sea Embayment with a spatial resolution of 1 km. To compute the total horizontal divergence rate from the ice velocities, we implement the logarithmic strain-rate formulation from Alley and others (2018) to produce continent-wide horizontal divergence rates. Designating u and v as the x and y components of the velocity field, respectively, in a polar stereographic coordinate system (EPSG: 3031), the two-dimensional strain-rate tensor is:

(7)$$\dot{\epsilon} = \left[ \ \ \ \ \matrix{{\partial u\over \partial x} \ \ \ \ \ \ \ \ {1\over 2}\lpar {\partial v\over \partial x} + {\partial u\over \partial y}\rpar \cr \hskip-2.7pc {1\over 2}\lpar {\partial v\over \partial x} + {\partial u\over \partial y}\rpar \ \ \ \ \ \ \ {\partial v\over \partial y} }\right].$$

If the strain-rate tensor is rotated with the local ice-flow direction, the total horizontal divergence rate, $\dot {\epsilon }_{\rm h}$, can be computed as the trace of $\dot {\epsilon }$, or the sum of the longitudinal and transverse strain rates, respectively:

(8)$$\dot{\epsilon}_{\rm h} = {\partial u\over \partial x} + {\partial v\over \partial y}.$$

If the strain-rate tensor is not with the local ice-flow direction, the longitudinal and transverse strain rates can be reoriented to calculate the horizontal divergence rate $\dot {\epsilon }_{\rm h}$ (see Alley and others (2018)). Errors can arise when generating the strain-rate tensor from a satellite-derived velocity field in areas of high strain when a nominal strain formulation is used (Alley and others, 2018). Therefore, we apply the logarithmic formulation from Alley and others (2018), which compares the change in length with the previous length, and not the original length, to account for the history of strain that is experienced by that region of the ice sheet.

We specify the accumulation rate, surface temperature and surface-snow density as the boundary conditions of the CFM. We force the CFM with a modified MERRA-2 climate reanalysis product (Smith and others, 2020, Brooke Medley, personal communication, 3 March 2020) at 5-day temporal resolution and 12.5 km spatial resolution. All accumulation rates we use in the model runs are in ice-equivalent units. Additionally, for all runs, we prescribe a constant surface-snow density of 400 kg m⁻³. Fausto and others (2018) discuss uncertainty arising from surface-snow density for firn-density calculations. While variable surface-snow density affects FAC estimates, we suspect it does not affect the relative change in FAC estimates from different horizontal divergence rates enough to alter our conclusions.

Results

Idealized runs

To provide quantitative insight into the impact of horizontal divergence on FAC, we conduct idealized simulations using the HL and LIG models with constant climate forcing under six horizontal-divergence-rate scenarios that span the range of horizontal divergence rates observed on the ice sheets (Experiment 1; Table 1 and Fig. 2). We choose surface-boundary conditions representative of central West Antarctica: a constant accumulation rate of 0.30 m a⁻¹, and a constant temperature of $-20^\circ$C (Kaspari and others, 2004; Steig and others, 2005; Medley and others, 2013; Fudge and others, 2016).

Fig. 2. Model response of FAC (b) to a step change in the horizontal divergence rate (a) from 0 to 1×10⁻⁴ a⁻¹ (solid light gray line), to 0 to 1×10⁻² a⁻¹ (dotted black line), using the LIG firn-compaction model. The model was forced with an accumulation rate of 0.30 m ice eq. a⁻¹, surface temperature of -20$^\circ$C and surface-snow density of 400 kg m⁻³. Evolution of the depth-density profile for a step-change in horizontal divergence rate of 5 × 10⁻³ a^-1 is shown in (c) as an example. Black line indicates the BCO depth. Contour interval is 50 kg m⁻³.

Table 1. Steps in horizontal divergence rates used in the layer-thinning scheme.

These encompass the range of horizontal divergence rates commonly observed on the ice sheets. The model was forced with an accumulation rate of 0.30 m ice eq. a⁻¹, surface temperature of $-20^\circ$C and surface-snow density of 400 kg m⁻³. Percent decrease values are shown using LIG; values from HL are shown in parentheses.

We spin the model up for 600 years to steady state under this constant climate with no horizontal divergence rate. Then, we run the model for 600 additional years, starting at time t = 0 years, using the same constant climate. We apply a step-change in horizontal divergence rate to the steady-state firn column at t = 100 years and track how the simulated FAC evolves for the next 500 years. We run this routine using HL and LIG for six horizontal-divergence-rate scenarios (Table 1 and Fig. 2).

Figure 2 shows the evolution of FAC through time predicted by LIG, including the evolution of the depth-density profile and bubble close-off (BCO) depth (i.e. density horizon of 830 kg m⁻³). After an initial adjustment period, the simulated firn column reaches a new steady state approximately 500 years following the onset of the horizontal divergence rate for all runs. In this new steady state, the firn parcels have thinned, and the FAC has correspondingly decreased (Fig. 3 and Table 1).

Fig. 3. Estimated firn-air content (FAC) using the layer-thinning scheme to account for horizontal divergence with the HL and LIG firn-compaction models. The greater the step-change in horizontal divergence rate, the greater the decrease in the FAC after the step change.

Adding a horizontal divergence rate of 10⁻³ a⁻¹ with the LIG (HL) model results in a FAC that is 4% (6%) less than for a model with no horizontal divergence (Experiment 1; Table 1 and Fig. 3). Imposing a horizontal divergence rate of 10⁻² a⁻¹ reduces the FAC by 31% (36%) compared to the no-horizontal-divergence scenario. HL estimates larger decreases in FAC than LIG for all idealized runs, and the difference between the models is larger for higher horizontal divergence rates (Fig. 3).

Spatial variability: flowline runs

To determine the impact of horizontal divergence on FAC in realistic climate and ice-flow conditions, we apply the layer-thinning scheme to two flowlines on Thwaites (THW – Experiment 2) and Pine Island (PIG – Experiment 3) Glaciers, West Antarctica (Fig. 4). We choose these flowlines because ice-surface speeds have been monotonically increasing in this area during the satellite record (since 1992; Wingham and others, 1998; Mouginot and others, 2014), and the region is thinning rapidly (Schröder and others, 2019). This makes mass-balance estimates from the region sensitive to the treatment of horizontal divergence in firn-compaction models. The surface temperature and accumulation rate increase non-linearly along both flowlines, as they approach the lower-elevation coast. Along the THW flowline, the surface temperature increases from $-27 ^\circ {\rm C}$ to $-18 ^\circ {\rm C}$, and the accumulation rate increases from 0.5 to 0.9 m a⁻¹. Along the PIG flowline, the surface temperature increases from $-26^\circ {\rm C}$ to $-17^\circ {\rm C}$, and the accumulation rate increases from 0.4 to 0.9 m a⁻¹. The ice speed along each flowline increases as the ice enters streaming flow, from speeds <10 m a⁻¹ to speeds >1000 m a⁻¹, with corresponding increases in along-flow horizontal divergence rates (from 0 to > 10⁻³ a⁻¹).

Fig. 4. Location of Experiment 2 (THW) and Experiment 3 (PIG) on a map of mean thinning rate for 1978–2018 (Schröder and others, 2019). The black star represents the location on lower Thwaites used in Experiment 4. Map is superimposed on Reference Elevation Model of Antarctica (REMA) ice-sheet surface elevation (Howat and others, 2019). Inset shows location of figure domain in Antarctica. The projection is polar stereographic (EPSG: 3031).

The starting points of the flowlines were chosen so that the flowlines would extend through the main trunks of the glaciers. We again spin up the CFM for 600 years with the mean 1980–2019 climate (accumulation rate and temperature) for the head of the flowline, and no horizontal divergence rate. We run the model twice for each flowline, with and without an imposed horizontal divergence rate. The firn column advects through the flowline based on the temporally static, but spatially variable 1996–2018 mean ice velocities from Mouginot and others (2019). We map the temporally evolving firn column to its position along-flow, and plot the associated FAC for each position along the flowline in Figures 5 and 6.

Fig. 5. Results from the layer-thinning scheme for the flowline on Thwaites Glacier using the LIG firn-compaction model (Experiment 2). (a) Horizontal divergence rates for the flowline. Horizontal divergence rates were derived from Mouginot and others (2019) following the approach of Alley and others (2018), and exclude compression. (b) The firn depth-density profile along the flowline for the model that accounts for horizontal divergence. Black line indicates the BCO depth. Contour interval is 50 kg m⁻³. (c) FAC results from model runs including the horizontal divergence rates shown in (a) (dotted line) and from a model without the horizontal divergence rates (dashed line).

Fig. 6. Results from the layer-thinning scheme for a flowline on Pine Island Glacier using the Ligtenberg and others (2011) firn-compaction model (Experiment 3). (a) Horizontal divergence rates for the flowline. Horizontal divergence rates were derived from Mouginot and others (2019) following the approach of Alley and others (2018), and exclude compression. (b) The firn depth-density profile along the flowline for the model that accounts for horizontal divergence. Black line indicates the BCO depth. Contour interval is 50 kg m⁻³. (c) FAC results from model runs with the horizontal divergence rates shown in (a) (dotted line) and from a model without horizontal divergence rates (dashed line).

Figure 5 shows model results with and without including horizontal divergence using LIG for the THW flowline. Horizontal divergence along the THW flowline results in a mean 3.7% (4.0%) difference for the entire flowline between the no-divergence and divergence runs using LIG (HL). However, the effect of horizontal divergence on modeled FAC increases toward the terminus of the glacier and influences FAC most within the last 150–200 km. At the end of the flowline, horizontal divergence causes the FAC to be up to 9.6 m (9.2 m) less than with a method that does not account for horizontal divergence (41% (40%) difference). For the PIG flowline (Fig. 6), horizontal divergence results in an average 0.68% (0.81%) difference over the entire flowline. Like THW, horizontal divergence has the greatest impact on FAC estimates for the last 150–200 km of the flowline, where the FAC is 4.4 m (4.8 m) less while using a model with horizontal divergence compared to a conventional firn-compaction model (18% (19%) difference).

Experiments 2 and 3 indicate that horizontal divergence becomes most important to consider in estimating the FAC in the most coastal 150–200 km of THW and PIG. Horizontal divergence rates increase at approximately 150–200 km before the end of the THW and PIG flowlines, and the firn column begins to thin there (Experiments 2 and 3; Figs 5 and 6). The thickness of the entire firn layer along the flowline can be seen in Figures 5b and 6b, where the black line shows the 830 kg m⁻³ density horizon (BCO depth). The firn thins non-uniformly along the flowline as a result of the variability of the horizontal divergence rates, and thinning or thickening induced by the spatially variable accumulation rate and temperature. Note that higher accumulation rates increase the FAC, whereas higher temperatures decrease the FAC. Thus, higher accumulation rates will offset the effect of horizontal divergence on the FAC, whereas higher temperatures will reinforce the net thinning effect. Results from using HL for both flowlines on THW and PIG show qualitatively similar results and therefore are shown in Section 5 of the Supplementary Material.

Temporal variability: static location on lower Thwaites Glacier

To investigate the effect of temporal variability of horizontal divergence on FAC estimates during the satellite record, we next consider the time evolution of FAC (ΔFAC) from 2007 to 2016 for a fixed location on lower THW (Experiment 4; black star in Fig. 4). We choose this location because the time series of ice speed is annually continuous, and mean thinning rates here are characteristic of a large portion of lower THW and of PIG. For spin up of the CFM, we use a constant horizontal divergence rate of 0.015 a⁻¹, based on observations from the beginning of the time series. We spin up the CFM for 600 years using a climate forcing randomly generated from the normal distribution of the 1980–2007 mean climate.

We perform four model runs, all of which use the 2007–2016 temperature and accumulation-rate fields from MERRA-2, with annual time steps (Smith and others, 2020, Medley, personal communication, 3 March 2020). The accumulation rate, temperature and divergence-rate boundary conditions are shown in Figure 7. Descriptions of each run are as follows:

1. (1) A baseline conventional firn-compaction-model run, which entails running the CFM with the evolving temperature and surface-accumulation rate but no horizontal divergence rate.

2. (2) A run with climate from (1), and a constant horizontal divergence rate of 0.015 a⁻¹ through the entire spin up and model run.

3. (3) A run with climate from (1); constant horizontal divergence rate of 0.015 a⁻¹ through the spin up until 2007; then the horizontal divergence rate evolves based on the 2007–2016 ice-velocity time series.

4. (4) Run with climate from (1); constant horizontal divergence rate of 0.015 a⁻¹ through the spin up until 1997; linear ramp up to a horizontal divergence rate of 0.04 a⁻¹ at 2007; then the horizontal divergence rate evolves based on the 2007–2016 ice-velocity time series.

Fig. 7. Surface boundary conditions, horizontal divergence rates and estimated FAC using the layer-thinning scheme with the LIG firn-compaction model for a location on lower Thwaites Glacier (Experiment 4). The model spin up from 1980 to 2007 is shown. Run 1 represents a conventional firn-compaction model run with no horizontal divergence. A constant horizontal divergence rate of 0.015 a⁻¹ is used in run 2. For runs 3 and 4, after spin up with a constant divergence rate of 0.015 a⁻¹, the model is run from 2007 to 2016 with temporally variable horizontal divergence rates derived from the Mouginot and others (2017) velocity time series. Run 4 also includes a linear ramp between horizontal divergence rates from the 1997 to 2007 values.

We choose these runs to demonstrate (a) the impact of horizontal divergence on FAC estimates through time, and (b) the effects of initializing the runs using firn columns with different initial conditions, which will result in different FACs due to the initial-state dependence of the firn-compaction rate (Fig. 7). The lower panel in Figure 7 shows the model-predicted FAC for the four runs. Horizontal divergence applied over longer histories has higher influence on the FAC. Run 1 (no imposed horizontal divergence rate) consistently predicts a FAC that is 7–8 m greater than the other runs. Runs 3 and 4 address the influence of temporal variability of horizontal divergence on FAC estimates, which is important for assessing ice-sheet mass balance from repeat satellite-altimetry observations. These runs indicate a marked decrease in FAC associated with greater horizontal divergence rates.

Temporally increasing horizontal divergence rates reduce the FAC substantially (Table 2 and Fig. 7). In Run 1, the FAC slightly increases by 0.20 m (0.17 m), or 0.77% (0.66%), almost a negligible change, from 2007 to 2016. In Run 2, the FAC also slightly increases by 0.08 m (0.12 m) or 0.42% (0.68%). In contrast, in Run 3, the FAC is reduced by 2.66 m (2.71 m), or 15.30% (15.55%). In Run 4, the FAC is reduced by 1.83 m (1.91 m), or 11.66% (12.08%). Even the short-term application of greater and time-variable horizontal divergence rates initialized in 2007 (Runs 3 and 4) produces a substantial difference in predicted FAC compared to the constant horizontal-divergence case (Run 2). FAC in 2016 is estimated to be 3.04 m (3.14 m) less in Run 3 than in Run 2 (18.75% (19.29%) difference) and is 3.86 m (4.00 m) less in Run 4 than in Run 2 (24.38% (25.21%) difference). FAC estimated in Run 3 and Run 4 begin to differ when the Run-4 horizontal divergence rate ramps up in 1997, ultimately resulting in a slightly lower FAC in 2016 (0.81 m (0.86 m); 5.70% (6.0%) difference) because the FAC in 1997 is already less in Run 4 than Run 3.

Table 2. Summary of results from temporally varying the horizontal divergence rate for the location on Thwaites Glacier from 2007 to 2016. Results using the LIG firn-compaction model in the layer-thinning scheme are shown, with results using the HL firn-compaction model in parentheses.

Discussion

In the following sections, we address the central questions motivating this study:

1. (1) What is the importance of horizontal divergence in controlling FAC?

2. (2) How do estimates of the time-evolution of FAC (ΔFAC) estimates change by including horizontal divergence in the calculations?

3. (3) Where does horizontal divergence matter for estimating FAC on the Antarctic Ice Sheet?

Through investigating these questions, we find that, firstly, neglecting horizontal divergence where divergence rates exceed 10⁻⁴ a⁻¹ will lead to an overestimate in FAC. This is because a firn column in regions with horizontal divergence is stretched horizontally and is thinner, and thus has less air content per unit volume than a firn column in regions without horizontal divergence. Second, accounting for horizontal divergence in FAC estimates results in a smaller calculated mass loss for regions of increasing horizontal divergence through time. This is because, for a given change in surface elevation Δh, the time-evolution of FAC (ΔFAC) is greater for regions with increasing horizontal divergence through time, which implies that the interpreted Δm is less than estimates that neglect horizontal divergence (Eqn (2)). Lastly, we find that horizontal divergence should be included in FAC estimates in regions entering and within the outlet glaciers and ice shelves of the Antarctic Ice Sheet. Horizontal divergence should be accounted for in estimates of the time-evolution of FAC in these regions, as speeds and thinning rates are currently increasing (Joughin and others, 2012; Mouginot and others, 2014; Smith and others, 2020).

How do FAC estimates needed for altimetry studies change by including time-evolving horizontal divergence?

Results from our transient runs suggest that including horizontal divergence makes a substantial difference in the calculated ΔFAC. Decadal-scale climate variability (Run 1) from 2007 to 2016 results in 0.20 m (0.17 m), a 0.77% (0.66%) increase. By comparison, ΔFAC for Run 2 is 0.08 m (0.12 m), a 0.42% (0.68%) increase. The only difference between Runs 1 and 2 is that Run 2 includes a constant horizontal divergence rate, which makes the total FAC of Run 2 lower than that of Run 1 and means that FAC changes for these two runs are solely a result of the variable climate. Runs 1 and 2 contrast Runs 3 and 4, which include the time-variable horizontal divergence and show substantial decreases in FAC through time (Table 2). These changes in FAC for Runs 3 and 4 constitute 16% and 11% of the observed thinning (27 m; Schröder and others, 2019), respectively, for this location on lower THW from 2007 to 2016. Surface lowering on a glacier is due to both mass loss and firn thinning or thickening; for a given observed elevation change Δh, if there is more firn thinning and thus higher ΔFAC, then there is less mass loss Δm.

Where do horizontal divergence rates matter for estimates of FAC on the Antarctic Ice Sheet?

To identify where horizontal divergence is important to consider in estimates of FAC, we first determine where the horizontal divergence-rate magnitude is comparable to the vertical strain-rate magnitude. We consider horizontal divergence negligible when the ratio of the magnitude of horizontal divergence rate to vertical strain rate is < 0.1.

Figure 8 shows the Antarctic-wide ratio R of the depth-averaged vertical strain rate $\overline {\dot {\epsilon }}_{zz}$ within the firn column, as calculated by the analytic model from Herron and Langway (1980), and the total horizontal divergence rate $\dot \epsilon _{h}$, as calculated from velocity data from Rignot and others (2017) following Alley and others (2018):

(9)$$R = {\dot\epsilon_{h}\over \dot{\epsilon_{zz}}}.$$

Not unexpectedly, R ≥ 0.1 occurs dominantly (1) in regions where ice enters the outlet glaciers along the margins of the Antarctic Ice Sheet, such as Pine Island Glacier, Thwaites Glacier, and other glaciers in the Amundsen Sea sector, and (2) in regions entering the ice shelves. Regions in the interior of the ice sheet have relatively low values of R. Based on this analysis, horizontal divergence can be neglected in broad interior regions of the ice sheet. Some high horizontal divergence rates on ice shelves are associated with ongoing rifting, and may not reflect changes in horizontal divergence rates that are of interest in this study.

Fig. 8. Ratio of the vertical and horizontal divergence rates (R) across the Antarctic Ice Sheet. Higher values of R show where horizontal divergence rates are significant in calculations of firn-air content.

Future work

We show, with a simple kinematic layer-thinning scheme in the CFM, that horizontal divergence must be accounted for in estimates of FAC in regions of dynamic ice flow. However, the effects of several additional processes not included in our model require further study: (1) the role of compressive strain in firn-thickness change, (2) densification during horizontal stretching, and (3) the role of brittle failure in reducing the effect of ductile thinning. To include compressive strain rates in our treatment of layer thinning, either better theoretical treatment of effects of compressional ice-flow stresses on firn compaction or intentional model calibration including high ice-flow-stress sites is needed. Most of the horizontal divergence rates in Experiments 2 and 3 are not compressive because these glaciers are not confined outlet glaciers and do not experience significant lateral compression associated with downstream narrowing of glacier extent. For these reasons, we chose to set any compressional strain rates in our runs to zero. We suspect that compressional stresses would not offset reduction in the FAC due to horizontal divergence, but would have an additive effect to increase firn density and lower FAC, because some microphysical mechanisms (e.g., power-law creep) are dependent on the square of the effective stress (Maeno and Ebinuma, 1983; Alley and Bentley, 1988). A dynamic treatment of horizontal extensional stresses in the firn would also contribute to densification during thinning and result in a further decrease in FAC for similar reasons. Previous studies refer to this effect as strain softening, and identify this process as having a greater impact on the firn-density profile than horizontal divergence in some locations (Riverman and others, 2019). Additionally, theoretical and observational evidence suggests that a decrease in the porosity will occur with an increase in the stress state due to increased damage and strain around bubbles (Chawla and Deng, 2005; Alley and Fitzpatrick, 1999). Finally, we ignore the effects of brittle failure in near-surface crevassing, which must accommodate some of the horizontal divergence in especially high strain-rate areas (Duddu and Waisman, 2012). Including a constitutive relation for density changes in response to horizontal stresses is an obvious next development. However, a dynamic treatment of horizontal stresses should be created along with a dynamic treatment of the vertical forces and vertical compression in firn. Additional experimental and field data are necessary for formulating such a model, and therefore is clearly beyond the scope of this paper.

Future work should therefore (1) collect more field measurements in regions of high horizontal ice-flow stresses; (2) intentionally calibrate models to high ice-flow-stress sites in an empirical framework; and/or (3) perform lab work investigating how density, viscosity and microstructures change under horizontal ice-flow stresses. Firn-compaction models have been developed using depth-density profiles that are assumed to have negligible thinning due to horizontal ice-flow stresses. This is true for many of those profiles, but the model development processes have not necessarily checked the validity of this assumption for each core; nor have data from high-stress regions been sought out. As a result, firn-compaction models assume there is no dependence on ice-flow stresses due to their functional form. In regions that may indeed have ice-flow stresses, calibrating the models using data from firn density there could lead to error in the calibrated coefficient estimates. However, calibrating additional coefficients in existing models using firn data from identified high ice-flow-stress regions could accommodate some of the impacts of horizontal ice-flow stresses on FAC. Also, new geophysical techniques for measuring time series of firn thickness and density, such as GNSS interferometric reflectometry (Larson and others, 2009; Gutmann and others, 2012) and autonomous phase sensitive radio echo sounding (Corr and others, 2002; Jenkins and others, 2006; Nicholls and others, 2015), are promising low-cost methods for gathering additional firn depth-density data in regions of dynamic ice flow. Work is ongoing to formulate a dynamic, time-dependent expression encompassing the effects of horizontal ice-flow stresses on firn-compaction processes, and new measurements (e.g., micro-CT scans) may provide the necessary data on firn-grain evolution to construct a model based on microstructure evolution. Micro-CT technology is starting to be applied to measure firn properties (Adolph and Albert, 2014; Gregory and others, 2014; Keegan and others, 2019). Further work to characterize the relative roles of ductile and brittle behavior of the firn will also allow for better characterization of firn in high-stress environments.

Conclusions

Estimates of spatially and temporally variable FAC are needed in calculations of ice-sheet mass balance derived from repeat-altimetry observations. Here, we introduced a method that accounts for firn-layer thinning from horizontal divergence into the CFM (Stevens and others, 2020). This scheme consists of (1) densification via an existing firn-compaction model, and (2) thinning of the firn via horizontal stretching due to horizontal divergence.

We assessed the spatial and temporal variability of changes in FAC due to horizontal divergence separately. Horizontal divergence becomes most impactful on FAC estimates in the last 100–150 km of the flowlines on Thwaites and Pine Island Glaciers, where horizontal divergence rates can reach 10⁻² a⁻¹ and higher. At the end of the Thwaites flowline, horizontal divergence causes the FAC to be 41% less than FAC estimates from a conventional firn-compaction model with no horizontal divergence. At the end of the PIG flowline, horizontal divergence leads to a 18% less FAC compared to results from a conventional firn-compaction model. For a representative location on lower Thwaites Glacier, a 15% decrease in FAC occurs from 2007 to 2016 due to horizontal divergence, which corresponds to 16% of the observed surface-elevation change. This contrasts output from a conventional firn-compaction model with no horizontal divergence, which estimates a 0.77% increase in FAC, due to climate variability alone. Horizontal divergence is most important to include in FAC estimates in outlet glaciers in the Amundsen Sea Embayment, and in regions entering and within the ice shelves of Antarctica.

We find that horizontal divergence is important for both (1) estimates of the steady-state FAC and (2) estimates of FAC variability in time. Neglecting horizontal divergence in FAC estimates where horizontal divergence rates exceed 10⁻⁴ a⁻¹ will lead to an overestimate in the steady-state and time-evolving FAC. Improved FAC estimates will produce better calculations of mass change from repeat surface-elevation observations, as well as more accurate estimates of basal-melt rates that depend on a hydrostatic assumption. Including horizontal divergence in FAC estimates will become more important as regions of the Antarctic Ice Sheet, such as the Amundsen Sea Embayment, continue to experience substantial increases in ice speed (Mouginot and others, 2014). Neglecting horizontal divergence within FAC estimates used in altimetry-derived mass-change calculations in these areas will lead to an overestimate in mass loss. Our work highlights the importance of accounting for horizontal divergence in estimates of FAC, and is a first step toward improving firn products that may consider adding the effects of strain thinning to produce better FAC outputs.