## Introduction

Many outlet glaciers of the polar ice sheets have accelerated and thinned markedly in the last 25 years (Joughin and others, Reference Joughin2012; Mouginot and others, Reference Mouginot, Rignot and Scheuchl2014; Smith and others, Reference Smith2020). 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, Reference Shepherd2012; Depoorter and others, Reference Depoorter2013; Shepherd and others, Reference Shepherd2019; Smith and others, Reference Smith2020). 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, Reference Gagliardini and Meyssonnier1997; Lüthi and Funk, Reference Lüthi and Funk2000). 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, Reference Maeno and Ebinuma1983; Burr and others, Reference Burr, Lhuissier, Martin and Philip2019), to processes related to melting and refreezing (Reeh and others, Reference Reeh, Fisher, Koerner and Clausen2005), and to ice-flow stresses (Alley and Bentley, Reference Alley and Bentley1988). 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, Reference Kuipers Munneke2015), despite observational evidence of its impact on firn-density structure (Christianson and others, Reference Christianson2014; Morris and others, Reference Morris2017; Riverman and others, Reference Riverman2019).

### 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, Reference Shepherd2012; Depoorter and others, Reference Depoorter2013; Shepherd and others, Reference Shepherd2019). 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:

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*:

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, Reference Lundin2017). Most firn-compaction models also are compatible with the suggestions by Robin (Reference Robin1958) under steady-state conditions. Robin (Reference Robin1958) suggested that, in steady-state conditions, the change in density with depth d*ρ*(*z*)/d*z* (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, Reference Morris2017):

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, Reference Zumberge1960; Crary and Charles, Reference Crary and Charles1961; Gow, Reference Gow1968; Kirchner and others, Reference Kirchner, Bentley and Robertson1979; Alley and Bentley, Reference Alley and Bentley1988; Christianson and others, Reference Christianson2014; Vallelonga and others, Reference Vallelonga2014; Riverman and others, Reference Riverman2019; Morris and others, Reference Morris2017).

### 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, Reference Christianson2014; Vallelonga and others, Reference Vallelonga2014; Riverman and others, Reference Riverman2019). Riverman and others (Reference Riverman2019) 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 (Reference Morris2017) 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 (Reference Morris2017) 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 (Reference Robin1958) 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, Reference Gagliardini and Meyssonnier1997; Lüthi and Funk, Reference Lüthi and Funk2000). 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, Reference Ligtenberg, Helsen and Van den Broeke2011; Li and Zwally, Reference Li and Zwally2015, 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, Reference Stevens2020), 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, Reference Joughin, Smith and Medley2014; Rignot and others, Reference Rignot, Mouginot, Morlighem, Seroussi and Scheuchl2014). 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, Reference Stevens2020). 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:

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:

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}}$:

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 (Reference Herron and Langway1980) model (HL) and the Ligtenberg and others (Reference Ligtenberg, Helsen and Van den Broeke2011) 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, Reference Li and Zwally2011; Ligtenberg and others, Reference Ligtenberg, Helsen and Van den Broeke2011; Morris and Wingham, Reference Morris and Wingham2014). 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* ^{−3}, 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* ^{−3} < *ρ* < 830 kg *m* ^{−3}. Robin (Reference Robin1958) 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, Reference Gow1969).

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, Reference Van Wessem2018), 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, Reference Shepherd2012; Gardner and others, Reference Gardner2013; McMillan and others, Reference McMillan2016; Shepherd and others, Reference Shepherd2019). LIG used 48 depth-density profiles from Antarctica to tune the firn-compaction model by Arthern and others (Reference Arthern, Vaughan, Rankin, Mulvaney and Thomas2010) for general applicability in Antarctica. Arthern and others (Reference Arthern, Vaughan, Rankin, Mulvaney and Thomas2010) 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, Reference Coble1970; Gow, Reference Gow1969), with a form similar to that of HL.

### Model inputs

#### Total horizontal divergence rates from ice velocities

We use the Mouginot and others (Reference Mouginot, Rignot and Scheuchl2019) 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 (Reference Mouginot, Rignot, Scheuchl and Millan2017) 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 (Reference Alley2018) 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:

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:

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 (Reference Alley2018)). 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, Reference Alley2018). Therefore, we apply the logarithmic formulation from Alley and others (Reference Alley2018), 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, Reference Smith2020, 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^{−3}. Fausto and others (Reference Fausto2018) 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^{−1}, and a constant temperature of $-20^\circ$C (Kaspari and others, Reference Kaspari2004; Steig and others, Reference Steig2005; Medley and others, Reference Medley2013; Fudge and others, Reference Fudge2016).

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^{−3}). 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).

Adding a horizontal divergence rate of 10^{−3} a^{−1} 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^{−2} a^{−1} 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, Reference Wingham, Ridout, Scharroo, Arthern and Shum1998; Mouginot and others, Reference Mouginot, Rignot and Scheuchl2014), and the region is thinning rapidly (Schröder and others, Reference Schröder2019). 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^{−1}. 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^{−1}. The ice speed along each flowline increases as the ice enters streaming flow, from speeds <10 m a^{−1} to speeds >1000 m a^{−1}, with corresponding increases in along-flow horizontal divergence rates (from 0 to > 10^{−3} a^{−1}).

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 (Reference Mouginot, Rignot and Scheuchl2019). 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.

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^{−3} 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^{−1}, 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, Reference Smith2020, 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) 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) A run with climate from (1), and a constant horizontal divergence rate of 0.015 a

^{−1}through the entire spin up and model run.(3) A run with climate from (1); constant horizontal divergence rate of 0.015 a

^{−1}through the spin up until 2007; then the horizontal divergence rate evolves based on the 2007–2016 ice-velocity time series.(4) Run with climate from (1); constant horizontal divergence rate of 0.015 a

^{−1}through the spin up until 1997; linear ramp up to a horizontal divergence rate of 0.04 a^{−1}at 2007; then the horizontal divergence rate evolves based on the 2007–2016 ice-velocity time series.

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.

## Discussion

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

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

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

(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^{−4} a^{−1} 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, Reference Joughin2012; Mouginot and others, Reference Mouginot, Rignot and Scheuchl2014; Smith and others, Reference Smith2020).

### 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, Reference Schröder2019), 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 (Reference Herron and Langway1980), and the total horizontal divergence rate $\dot \epsilon _{h}$, as calculated from velocity data from Rignot and others (2017) following Alley and others (Reference Alley2018):

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.

### 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, Reference Maeno and Ebinuma1983; Alley and Bentley, Reference Alley and Bentley1988). 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, Reference Riverman2019). 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, Reference Chawla and Deng2005; Alley and Fitzpatrick, Reference Alley and Fitzpatrick1999). 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, Reference Duddu and Waisman2012). 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, Reference Larson2009; Gutmann and others, Reference Gutmann, Larson, Williams, Nievinski and Zavorotny2012) and autonomous phase sensitive radio echo sounding (Corr and others, Reference Corr, Jenkins, Nicholls and Doake2002; Jenkins and others, Reference Jenkins, Corr, Nicholls, Stewart and Doake2006; Nicholls and others, Reference Nicholls2015), 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, Reference Adolph and Albert2014; Gregory and others, Reference Gregory, Albert and Baker2014; Keegan and others, Reference Keegan, Albert, McConnell and Baker2019). 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, Reference Stevens2020). 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^{−2} a^{−1} 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^{−4} a^{−1} 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, Reference Mouginot, Rignot and Scheuchl2014). 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.

## Supplementary material

The supplementary material for this article can be found at https://doi.org/10.1017/jog.2020.105.

## Code and data availability

The CFM code is publicly available at https://github.com/UWGlaciology/CommunityFirnModel. Documentation for the CFM is online at https://communityfirnmodel.readthedocs.io/TS5. Model output and scripts used to make the figures will be available on the University of Washington ResearchWorks Archive.

## Acknowledgments

We thank Brooke Medley for her modified MERRA-2 climate reanalysis data. We also thank Tyler J. Fudge and Michelle Koutnik for their comments on earlier versions of the manuscript. This work was supported by NASA grant NNX16AM01G and NSF grant 0968391.

## Author contributions

ANH, KC, EDW and CMS designed the study. ANH implemented the layer-thinning procedure, ran the model experiments and led writing of the paper. CMS aided implementation of the layer-thinning scheme in the CFM and performed some model runs. NH and KC contributed to analysis of ice-velocity data. All authors contributed to paper writing and editing.