A new model of dry firn-densification constrained by continuous strain measurements near South Pole

Abstract

Converting measurements of ice-sheet surface elevation change to mass change requires measurements of accumulation and knowledge of the evolution of the density profile in the firn. Most firn-densification models are tuned using measured depth–density profiles, a method which is based on an assumption that the density profile in the firn is invariant through time. Here we present continuous measurements of firn-compaction rates in 12 boreholes near the South Pole over a 2 year period. To our knowledge, these are the first continuous measurements of firn compaction on the Antarctic plateau. We use the data to derive a new firn-densification algorithm framed as a constitutive relationship. We also compare our measurements to compaction rates predicted by several existing firn-densification models. Results indicate that an activation energy of 60 kJ mol⁻¹, a value within the range used by current models, best predicts the seasonal cycle in compaction rates on the Antarctic plateau. Our results suggest models can predict firn-compaction rates with at best 7% uncertainty and cumulative firn compaction on a 2 year timescale with at best 8% uncertainty.

1. Introduction

Snow that accumulates on the polar ice sheets transitions to ice over the course of hundreds to thousands of years through an intermediate stage known as firn. Knowledge of the rate of firn compaction is needed both (1) to extract climate records from ice cores, because the rate of firn compaction affects the age difference between gasses trapped in the ice and the ice itself; and (2) to determine the rate of ice-sheet mass loss from repeat satellite altimetry, where elevation-change (volume) measurements need to be converted to mass change. These applications use firn-densification models to estimate the depth–age and depth–density profiles of the firn and how they change through time. Modern advances in ice-core science provide annual to sub-annual records of climate (Sjolte and others, 2020; Jones and others, 2023), and altimetry techniques provide centimeter-scale measurements of ice-sheet elevation (Shepherd and others, 2019; Smith and others, 2020). Taking full advantage of these high-resolution data requires models that can predict firn evolution accurately on sub-annual timescales and in variable and changing climates. Further, these applications require quantitative understanding of the uncertainties in firn-model calculations. Improving firn models relies on high-resolution measurements of firn processes, including compaction rates and microstructure evolution.

1.1 Firn densification

In most models, firn is divided into three stages based on the dominant physical processes driving densification. Stage 1 extends from the surface, where the density is often assumed to be between 300 and 400 kg m⁻³, to the depth of the 550 kg m⁻³ density horizon; densification in this stage is thought to be primarily due to grain boundary sliding (Alley, 1987). Stage 2 extends from 550 kg m⁻³ to the bubble close off (BCO) density, which varies by site but is commonly assumed to be near 830 kg m⁻³. Densification in stage 2 is thought to be due to sintering processes including power-law creep and lattice diffusion (Wilkinson, 1988). At the BCO density, air becomes trapped in bubbles in the ice. Beyond the BCO density, densification continues due to overburden pressure deforming ice around the pores, but increasing pressure in the bubbles slows the densification significantly (Goujon and others, 2003).

Firn densification is driven by processes acting at the scale of individual grains (order of a mm; Blackford, 2007); these processes can be simulated in a model using a microstructural approach, a macroscale approach or a hybrid of the two (Morris and Wingham, 2014). The microscale approach models the evolution of microstructural properties directly and predicts the densification rate based on those properties. For example, Alley (1987) proposed a grain-boundary sliding model for stage 1 by considering the size and contact area of idealized spherical grains. Alternatively, macroscale models relate densification to site-specific properties such as temperature and accumulation rate using empirically derived coefficients; these coefficients can be interpreted as parameterizations of the microstructural processes. Macroscale models often use dated depth–density profiles to convert the change in density ρ with depth z, ∂ρ/∂z, to change in density with time t, ∂ρ/∂t, i.e. the densification rate. This method assumes the firn depth–density profile is invariant through time, i.e. that it is in steady state (Sorge's law; Bader, 1954). The benchmark Herron and Langway (1980) model is an example of a macroscale firn-densification model. Hybrid models use a combination of microscale and macroscale approaches; for example, Arthern and others (2010) measured firn-compaction rates in Antarctica and derived a model based on previously derived expressions for creep of ice and porous media, and they integrated a grain-growth equation into their model formulation. They then used this information to determine best-fit coefficients in a macroscale firn-densification expression similar to that of Herron and Langway (1980).

1.2 Measuring firn compaction

Efforts to validate the accuracy of firn-densification models are hampered by the fact that there have been relatively few measurements of firn-compaction rates (i.e. vertical strain rates). These measurements provide temporal information, whereas depth–density profiles only provide a snapshot of the firn at one moment in time. Comparing a modeled density profile to a measured depth–density profile does not provide information about how well a model simulates compaction in a transient climate or in a climate with natural variability (i.e. variability in temperature and surface mass balance on sub-daily to inter-annual timescales). Measuring firn compaction directly provides opportunities to develop firn-compaction models using measured (rather than inferred) strain rates and to compare measured compaction rates against those predicted by firn densification models. This is particularly important for altimetry products that utilize simulations of firn-thickness change at daily to weekly timescales to separate the portions of surface elevation variations due to new accumulation, ice dynamics, and firn-air content (FAC) change (e.g. Khan and others, 2022; Smith and others, 2023).

Several methods have been used to measure firn compaction. Hulbe and Whillans (1994) developed a method to measure firn compaction by placing poles at the bottom of firn boreholes near South Pole and measuring the surface elevation relative to the poles; Hamilton and others (1998) further described this method and termed it the ‘coffee-can’ method because the poles were frozen into coffee cans at the bottom of the boreholes. These measurements were manual, so firn compaction could be measured only when an observer was present. Arthern and others (2010) automated this method using a temporally continuous data-logging system attached to a draw-wire sensor. The FirnCover project (MacFerrin and others, 2022) used similar instrumentation to measure firn compaction in 48 boreholes across eight sites in Greenland. Other methods have measured firn compaction by measuring the distance between pairs of firn layers at different times. Morris and Wingham (2014) used a neutron probe to make repeat, high-resolution density logs in boreholes in Greenland; they derived compaction rates by determining how the distance between adjacent layers changed between logs. Hawley and Waddington (2011) used borehole optical stratigraphy to determine firn-compaction rates at Summit, Greenland; this technique used a camera to track relative position changes of optical features in a borehole. Hubbard and others (2020) also used this technique to measure compaction at Derwael Ice Rise in Antarctica, using a high-resolution camera that allowed for measurements with millimeter-scale resolution. Both the neutron probe and optical logging methods provide high-spatial-resolution data, but their temporal resolution is limited to the frequency of visits to the borehole sites.

Additional methods circumvent the need for boreholes, but require an estimate of the firn-density profile to determine the wave speed of radar through the firn. Medley and others (2015) derived firn-compaction rates from repeat airborne-radar surveys that were constrained by overflying ITASE boreholes with measured density-age profiles. This method provides compaction-rate data over a large spatial area, but its temporal resolution is limited by the frequency of radar flights and its vertical resolution is limited by the radar's resolution (at best 10 cm, depending on radar frequency and bandwidth). Case and Kingslake (2021) used autonomous phase-sensitive radar to derive firn-compaction rates by tracking relative positions of internal reflectors in the firn through time.

Here, we made continuous firn-compaction and temperature measurements over 2 years near South Pole. To our knowledge these are the first continuous measurements of firn compaction on the Antarctic plateau. We use these measurements to derive a new firn-densification equation using a continuum-mechanics framework. We then compare the measurements to predictions from several previous firn-densification equations.

2. Methods

2.1 Field measurements

We measured firn compaction continuously for 2 years using the ‘coffee-can’ method in 12 boreholes at USP50 (89.54$^\circ$S, 137.04$^\circ$E), a site 50 km upstream from the South Pole on the flowline (Lilien and others, 2018). The depths of the boreholes were 4.4, 4.42, 9.65, 9.75, 14.65, 15.1, 19.53, 24.85, 29.82, 40.16, 80.0 and 106.0 m. The layout of the site is shown in Figure 1b. For simplicity, we label and refer to the boreholes and instruments installed in them by rounding their depths to the nearest integer meter. Duplicate-depth boreholes and instruments are referred to using letters a and b after the integer.

Figure 1. (a) Cartoon showing how the ‘coffee-can’ type instruments measure firn compaction. (b) Map-view schematic of the experiment site. Black dots show locations of boreholes, which are labeled by the borehole depth rounded to the nearest integer meter. (c) Photographs of the experiment setup.

We measured compaction rates in the firn by installing draw-wire string potentiometers in the 12 boreholes. Figures 1a, c show a schematic of the instrument design and photographs of the setup, respectively. The instruments were enclosed in wooden boxes, which were attached to wooden platforms with an area of ~0.4 m². The platforms were installed 20–25 cm below the surface (Fig. 1 photographs). A length of braided Vectran string connected each potentiometer to an anchor (0.56 kg lead weight) at the bottom of its boreholes. As the firn compacted, the potentiometers reeled in string to maintain constant tension. The potentiometers operate as variable resistors, and thus a measurement of voltage drop across the instrument can be converted to a measurement of cable length extending out the potentiometer.

We measured density on the core recovered from the 106 m borehole by sampling a ~10 cm portion of every ~100 cm of core; each section's length, diameter and mass were recorded. We also measured near-surface density and firn temperature in a 2 m deep snow pit. An additional 126 m core was recovered and returned to the NSF Ice Core Facility (ICF), where we measured density of each core section and determined the depth–age relationship with electrical conductivity measurements (Lilien and others, 2018) and major ions (Winski and others, 2019). We also attempted to measure a high-resolution density profile in the 126 m open borehole using a neutron probe. The probe failed, which we hypothesize was due to the cold temperatures in the borehole, although neutron probes have been used successfully on the Antarctic plateau, notably on a traverse through Queen Maud Land in 1967–68 (Kane, 1969).

We also measured firn temperature using a 40 m thermistor string with 23 sensors in a backfilled borehole. The depths of the thermistors at time of installation were 0, 0.25, 0.5, 0.75, 1.00, 1.25, 1.5, 2.0, 2.5, 3, 4, 5, 6, 8, 10, 12, 14, 17, 20, 25, 30, 35 and 40 m.

Campbell CR1000 dataloggers recorded the change in length of each borehole and the temperature of each thermistor every 6 h from 9 January 2017 to 24 December 2018. Power was provided by a battery bank in an insulated box connected to a solar-panel array.

2.2 Data processing

The change in cable extension between two measurement epochs is equal to the change in borehole length ΔL during that time. We define borehole length as the distance from the instrument platform to the borehole bottom. For each borehole, we calculate borehole length through time, L(t), by subtracting the cumulative borehole length change from the initial length L ₀. We discard the first month of data to remove any effects of instrument settling and cable strain (MacFerrin and others, 2022), and then we filter the displacement time series using a rolling Gaussian window (window size = 30 observations (7.5 d), std dev. = 10 observations (2.5 d)). We define the compaction rate as ΔL divided by the time Δt between measurements.

We calculate the strain $\epsilon$ at each measurement time using the definition of logarithmic strain:

(1)$$\epsilon( t) = \ln{\left({L( t) \over L_{0}}\right)},\; $$

and take the time derivative of $\epsilon ( t)$ to calculate the strain rate $\dot {\epsilon }( t)$:

(2)$$\dot{\epsilon}( t) = {1\over L( t) } {{\rm d}L( t) \over {\rm d}t}.$$

We assume that strain due to horizontal divergence (Horlings and others, 2021) is negligible, so the vertical strain rate is equal to the volumetric strain rate. We also use the depth–density profile to calculate the overburden stress (with the convention of negative stress being compressive) at depth z:

(3)$$\sigma( z) = - g \int_{0}^{z} \rho( z') \, {\rm d}z',\; $$

where g is the gravity and ρ is the firn density. Instrument 10a failed during the first year of measurements and was replaced on 1 January 2018. We estimate the length of borehole 10a at the time of replacement installation by assuming that the compaction in borehole 10a was the same as in borehole 10b during the first year of the experiment.

3. Results

3.1 Firn core measurements

Figure 2 shows the depth–age, depth–density and annual accumulation rate data from the 126 m firn core. The firn at the bottom of the core had a density ~800 kg m⁻³ and its age was 1021 years. The mean accumulation for the 1021 years was 69.3 kg m⁻² a⁻¹ (0.076 m ice eq.). Consistent with Fudge and others (2020), our data show that the accumulation rate has increased in the past 150 years: the mean accumulation for the past 150 years (calendar years 1868–2017) was 79.4 kg m⁻² a⁻¹ (0.087 m ice eq.), and the mean for the 871 years previous to that (calendar years 996–1867) was 67.6 kg m⁻² a⁻¹ (0.074 m ice eq.).

Figure 2. Firn core data from USP50: (a) depth–age profile; (b) depth–density profile and (c) derived accumulation rate for the past 1021 years. The orange curve shows the annual record, and the black curve is the 10 year running mean.

3.2 Firn compaction

Figure 3 shows compaction rates measured in each of the 12 boreholes from February 2017 to December 2018. The density at 106 m depth was ~800 kg m⁻³, which is slightly less than the commonly cited bubble close-off density of 830 kg m⁻³ (e.g. Blunier and Schwander, 2000). Firn cores from South Pole have shown that the BCO depth is near 120 m (Mayewski and Dixon, 2005; Winski and others, 2019). For our analyses, we consider observations in the 106 m borehole to be representative of the full firn column depth because there is very little additional compaction that occurs beyond 106  m. Compaction in this deepest borehole was 0.262 m over 680 d, which is an average rate of 0.141 m a⁻¹. For comparison, the steady-state compaction rate between a 300 kg m⁻³ surface and the 800 kg m⁻³ density horizon, assuming an accumulation rate of 69.31 kg m⁻² a⁻¹, is 0.144 m a⁻¹.

Figure 3. Compaction rates measured in each of the 12 boreholes from February 2017 to December 2018. Blue shading indicates winter months, and red shading indicates summer months. Note that the y-axis scales differ among the panels.

The compaction signal is dominated by compaction in the upper firn. Total compaction in the two 4 m deep boreholes over the same period was 0.083 and 0.088 m, or about a third of that in the entire firn column. The total compaction in the two 15 m deep boreholes was 0.155 and 0.149 m, almost 60% of the full-column compaction.

The data have a clear seasonal signal. The highest compaction rates occur in summer (late December through mid March; highlighted in red in Figure 3), and the lowest rates occur in late winter (blue) into early spring. The minimum compaction rate in the 106 m borehole was 0.093 m a⁻¹, which occurred in August 2017. The compaction rate more than doubles in the summer: the maximum observed was 0.254 m a⁻¹ in February 2017, and its maximum is during the second summer of observation was 0.221 m a⁻¹ in January 2018. In general, the compaction rates during the second year of measurement are expected to be lower because the firn spanned by the strain meter is slightly denser in the second year (firn viscosity increases with density); however, this effect will be altered by temperature differences between the 2 years (firn viscosity decreases with increasing temperature and vice versa).

The seasonal increase and decrease in compaction rate is not symmetrical; the compaction rate increases quickly in the spring but decreases more slowly in the fall. This effect is well-captured by a firn model (Section 5.3). We do not fully analyze here how a symmetrical surface temperature signal is transformed into an asymmetric compaction rate, but we speculate that it is due to the compaction rate having a non-linear dependence on temperature and density (Section 5.1). In addition, the firn spanned by the instrument is denser at the end of the summer than the beginning, which could also contribute to this asymmetry.

The variability in compaction rates in the boreholes was also highly correlated (Fig. 4); for example, correlating the filtered compaction-rate time series for instruments 80 and 40 gave r = 0.97, which indicates that the dominant factor driving compaction on weekly and longer timescales is the same for these boreholes.

Figure 4. (a) Temperatures measured in the top 10 m of firn from February 2017 to December 2018. (b) Compaction rates for the 40 m (gray) and 80 m (brown) boreholes for the first 2 months of observation (February–April 2017). Despite their spatial separation, the compaction rates are highly correlated with each other (r = 0.97 for the filtered data over the full 2 years of data).

3.3 Firn temperature

Figure 4 shows the temperature measurements from the thermistor string for the upper 10 m of firn (we do not show measurements deeper than 10 m because the firn is nearly isothermal at those depths). The uppermost thermistor, which began at the surface but was slowly buried throughout the experiment, measured a maximum temperature of $-23.10^\circ$C and a minimum temperature of $-64.45^\circ$C. The thermistor that was installed at 10 m depth had a maximum temperature of $-50.55^\circ$C and a minimum temperature of $-51.25^\circ$C. The deepest four thermistors (25, 30, 35, 40 m) had respective mean temperatures of $-51.13$, $-51.19$, $-51.05$ and $-51.25^\circ$C, and the temperatures measured on these four thermistors varied by at most 0.25$^\circ$C.

3.4 Measurement uncertainty

We duplicated our instrument setup for the 4, 10 and 15 m boreholes to assess the uncertainty of our method. As previously stated, one of the instruments in a 10 m borehole failed during the first year, but the second year provided useful data for comparison. We measured 0.0845 and 0.0897 m of total shortening in boreholes 4a and 4b, respectively, during the experiment; this is approximately a 6% difference between the two. The total shortening in boreholes 10a and 10b during 2018 was 0.0700 and 0.0755 m ($\sim \!7.5\%$ difference). Total shortening in boreholes 15a and 15b for the entirety of the experiment was 0.1576 and 0.1508 m ($\sim \!4.4\%$ difference). Averaging these, we estimate a measurement uncertainty of 6%.

The agreement between measurements over common depth ranges gives us confidence that our instruments were reliably measuring the firn-compaction signal. We attribute the differences between measurements to two factors: instrument precision and spatial variability. Our measurements do not allow us to differentiate between these two factors. The instrument uncertainty includes the potentiometers themselves as well as the larger setup; e.g. one platform could settle differently than another, which does not directly reflect natural firn compaction. Regardless of whether the uncertainty is due to instrument precision or spatial variability, our data indicate that the ‘coffee-can’ method is a robust method of measuring firn compaction.

4. Firn viscosity

Our high-temporal-resolution measurements allow for a more detailed examination of the mechanical properties of the firn than discrete measurements of depth–density profiles allow. We follow the work of Bader (1960), who defined a one-dimensional ‘compactive viscosity factor’ relating local strain rate of a layer of firn to the overburden stress. However, because our strain measurements span many meters of firn over which the firn properties, including density and stress, vary, we define a ‘parcel viscosity factor’ η _p (hereafter referred to as parcel viscosity) relating the stress σ at a depth z to the vertical strain rate $\dot {\epsilon }$:

(4)$$\dot{\epsilon} = {1\over 2\eta_{\rm p}}\sigma.$$

We calculate the ‘parcel viscosity’ of a parcel of firn using the bulk density of that parcel and the stress at the bottom of the parcel. The parcel viscosity is thus a ‘bulk’ viscosity value which describes the macroscale deformation behavior of a parcel of firn subjected to a load, and it differs from the viscosity of a small firn sample of the same bulk density with nearly homogeneous properties (i.e. a representative elementary volume, or REV).

Our goal is to find an expression to predict η _p as a function of density ρ, temperature T and potentially other variables, so that η _p can then be used in Eqn (4) to predict the strain rate in the firn. We note that it is possible that the parcel viscosity may also be a function of σ as assumed in previous firn research (e.g. Alley, 1987; Barnola and others, 1991), but in the present work we assume the viscosity is independent of stress.

We rearrange Eqn (4) and use our calculated strain rates (Eqn (2)) and the stresses at the bottoms of the boreholes to calculate the parcel viscosity of the firn spanning ‘depth intervals’. We use the term ‘depth interval X–Y’ to refer to the firn between the boreholes with bottom depths of X and Y m. We find the parcel viscosity for depth intervals between borehole bottoms by subtracting the compaction measured in boreholes of different depths and following the steps described above. For example, to find the amount of compaction that occurred between the depths of X = 4 and Y = 10 m (depth interval 4–10), we subtract the compaction measured in the shallower hole X from the compaction measured in the deeper hole Y. We use that compaction to calculate the strain rate for depth interval X–Y, and then we use that strain rate and the stress at the bottom of hole Y to find the parcel viscosity of the firn in depth interval X–Y. This method assumes that the compaction measured in the shallower borehole X is equal to the compaction in the upper portion of the deeper borehole Y; compaction measurements in boreholes of similar depths indicate that this assumption is valid. We calculate the strain rates and parcel viscosities for each pair of boreholes with adjacent depths (i.e. 4–10, 10–15 m, etc.). To do so, we down-sample our data from 6 h to weekly resolution by taking the median value of the parcel viscosity in each depth interval for each week of the observation period.

Figure 5 shows the firn parcel viscosity calculated for each depth interval. Each blue dot is the parcel viscosity calculated for 1 week of observation, and the blue vertical bars are the median of those weekly parcel viscosities for each depth interval. There are multiple vertical bars in several of the depth intervals because the boreholes with similar depths allow us to calculate the parcel viscosity for multiple pairs of boreholes spanning those depth intervals (e.g. boreholes 4a–10a, 4a–10b, 4b–10a, 4b–10b).

Figure 5. Firn parcel viscosities derived from compaction-rate measurements and steady-state-derived viscosities for various depth intervals. Blue dots show the derived parcel viscosities for each week of observation, and the blue vertical lines show the median values of those weekly values. The orange curve shows the viscosity derived from depth–density data using a steady-state assumption (Eqn (A.5)), and the orange vertical lines are the median steady-state viscosities in the depth intervals. The meter-scale variability in the ‘steady-state’ profile is a result of the variability in the depth–density profile.

We can also estimate the theoretical steady-state firn viscosity (Appendix A) implied by a depth–density profile. The orange curve in Figure 5 shows the steady-state firn viscosity profile predicted by Eqn (A.5) using our measured depth–density data and an accumulation rate of 69.37 kg m⁻² a⁻¹ (the mean accumulation rate for the 1021 years recorded in the 126 m core). The curve fits our measured parcel viscosities qualitatively well, suggesting that (1) the parcel viscosities are a good approximation of the firn's compactive viscosity and (2) the firn near South Pole is in or near steady state. The higher accumulation rate at USP50 in the last 150 years compared to the years prior may explain part of the misfit between our observations and the steady-state viscosity approximation in the deep firn.

5. Firn-densification modeling

5.1 A new firn-densification equation

Our measurements allow development of a firn-densification equation framed as a constitutive relationship, in contrast to many previous firn-densification equations that have been developed by converting depth–density–age data to densification rates. The advantage of a constitutive-relation model form is that it should, in theory, be able to reliably predict firn densification on various timescales in variable and transient climates, whereas an equation based on fitting to depth–density profiles may have limited applicability in these applications (Lundin and others, 2017). Note that while we are using a constitutive relationship form tuned with strain-rate data, we are using an empirical, macroscale approach for our proposed equation.

We use our strain-rate and density data, along with our derived weekly parcel viscosities, to derive an expression describing firn densification at USP50. We follow previous work (e.g. Morris and Wingham, 2014) and assume that the viscosity increases as the porosity decreases (i.e. η ∝ 1/(ρ _i − ρ)), that the densification rate has an Arrhenius dependence on temperature T (unit: K) with activation energy Q, and that the strain rate is linearly related to the stress. It is possible that the firn strain rate has a non-linear stress dependence, as ice visco-plasticity has, but we do not explore that possibility here. We begin by assuming that Q = 60 kJ mol⁻¹ (Arthern and others, 2010), but subsequently test this assumption with other values for Q. We also assume that the viscosity increases as the firn ages due to processes operating at the grain scale, such as normal grain growth (Gow, 1969). In absence of data that would allow us to determine and model those processes specifically, we introduce a factor f(τ) that is a function of the firn's age τ (years). The equation describing densification then has the form:

(5)$$\dot{\epsilon} = -{1\over \rho}{{\rm D}\rho\over {\rm D}t} = {1\over 2\eta}\sigma = {( \rho_{i}-\rho) {\rm e}^{( -Q/RT) }\over K( \rho) f( \tau) }\sigma,\; $$

where K(ρ) is a density-dependent prefactor (unit: kg² m⁻⁴ s⁻²).

Several previously developed firn-densification models have similarly used a constitutive relation to predict the densification rate (e.g. Alley, 1987; Arnaud and others, 2000); Arthern and others, 2010; Morris and Wingham, 2014). Arthern and others (2010) and Morris and Wingham (2014) both used strain-rate measurements to define firn-densification equations, though they did not explicitly define the viscosity. Arthern and others (2010) used a grain growth term in the implied viscosity: as grains grow through time, the firn becomes more viscous. Morris and Wingham (2014) introduced a ‘temperature-history’ term, which implies the evolution of the viscosity depends on past temperature in the firn. Though the formulations of these models are different from each other, both imply that the firn's viscosity increases with time. Our proposed equation is also based on the assumption that firn viscosity increases through time. However, we do not have reliable data that could be used to identify the specific physical mechanisms that stiffen the firn as it ages (or to develop an expression to describe those mechanisms in a numerical model). Therefore, we chose to simply let f(τ) = τ. That is, we use the firn's age explicitly as a model parameter because it is a measurable quantity.

Equation (5) implies that the viscosity of the firn is:

(6)$$\eta = {K( \rho) \tau\over 2 ( \rho_{i}-\rho) {\rm e}^{( -Q/RT) }}.$$

The quantities σ, ρ and τ are known from our data. They did increase throughout the period of observation, but we assume those changes were small enough over the duration of our experiment to consider them constant. However, Eqns (5) and (6) include the viscosity, whereas our data allowed us to calculate the parcel viscosity. For the development of our densification equation, we assume that the parcel viscosity for each depth interval and each week is equal to the viscosity of a firn parcel with homogenous properties. Specifically, we assign these theoretical parcels the mean density, maximum age and weekly mean temperature for each depth interval. We then use Eqn (6) to calculate values of K(ρ) for each depth interval for each week of the measurements. Those values are shown as blue dots in Figure 6. We exclude values outside the interquartile range of K(ρ) in order to exclude outliers, which are likely due to spatial variability, as the parcel viscosity in depth intervals is based on differencing compaction rates in boreholes several meters apart.

Figure 6. Values of the prefactor K (blue dots) derived from parcel viscosity data in Figure 5 and modeled fit to those K values (red line).

Our use of a prefactor K(ρ) follows previous work (e.g. Herron and Langway, 1980; Arthern and others, 2010) that used coefficients to tune densification models to observations. Although Herron and Langway (1980) and Arthern and others (2010) used constant values of K for each densification stage (thereby implying an abrupt transition from stage 1 to stage 2), we use the suggestion from Morris (2018) that there should be a smooth transition from stage 1 to stage 2 densification. We assume that K is a logistic function of density with the following form:

(7)$$K( \rho) = {L\over 1 + {\rm e}^{-a( \rho-\rho_{\rm c}) }} + b.$$

We use the curve_fit function in SciPy's Optimize package to fit Eqn (7) to the K(ρ) values derived from our parcel viscosity data. The optimal values for the variables are L = 9.52 × 10⁻⁷ kg² m⁻⁴ s⁻², a = 4.11 × 10⁻² m³ kg⁻¹, b = 2.82 × 10⁻⁷ kg² m⁻⁴ s⁻² and ρ _c = 515.6 kg m⁻³. The red line in Figure 6 shows the modeled fit to K(ρ).

We perform the same tuning procedure using different values for Q to test the assumption that Q = 60 kJ mol⁻¹ (Section 5.4).

5.2 Simulations using the Community Firn Model

We used the Community Firn Model (CFM; Stevens and others, 2020) to simulate firn compaction using our new equation. Additionally, we ran the CFM using five previously published firn-densification equations (with abbreviations): Herron and Langway (1980, HL), Ligtenberg and others (2011, LIG), Li and Zwally (2015, LZ), CROCUS (Vionnet and others, 2012, CRO) and Medley and others (2022, GSFC). Note that each of these equations use globally, rather than locally, tuned parameters. Each model run simulated firn compaction from February 2017 to December 2018 with daily time steps. We used our measured depth–density profile as an initial condition and our daily observations of the firn–temperature profile to set the firn temperature in the model at each time step.

We did not measure accumulation rate continuously at our site. Our only accumulation measurements occurred during site visits when we measured the previous year's total accumulation. The CFM requires accumulation at each time step (daily in this case) as a model input. We estimated this by scaling the accumulation measured at South Pole station to match our annual snow-accumulation observations at USP50. This scaled accumulation rate likely has significant uncertainty. However, this uncertainty does not significantly affect the model results because the firn-densification equations we used are dependent on either the stress or the long-term accumulation rate, neither of which vary significantly on a 2 year timescale.

The CFM uses a Lagrangian (layer-following) grid scheme, which mimics our experimental setup exactly. In the model, we track the compaction predicted between a layer near the surface (the platform) and another layer at some depth (the anchor) as that system is advected downward. We calculated the firn-compaction rate in the model by tracking the distance between the model grid layers that were closest to our field measurements; for example, to compare model predictions with measurements from instrument 4a, we found the model grid layers nearest 0.25 and 4.4 m depth and tracked the distance between those layers for the duration of the model run.

5.3 Model outputs and compaction-rate data

We compare our measurements of compaction rates against values predicted by the CFM. We consider both the variation of the compaction rates throughout the year and the cumulative compaction over the course of the measurements.

5.3.1 Weekly to seasonal compaction rates

Figure 7 shows compaction rates predicted by the six firn-densification equations and the measured compaction rates for each borehole.

Figure 7. Measured and modeled compaction rates for each borehole.

The models show a range of responses and vary in their ability to predict the compaction rates. We calculated the root mean square deviation (RMSD) and normalized RMSD (NRMSD; calculated by dividing the RMSD by the mean compaction rate in each borehole over the period of observation) between each borehole's observed compaction rate and the models’ predictions. The RMSD and NRMSD for each model and each borehole are shown in Table 1. The RMSD for each model is roughly constant regardless of borehole depth, which means that the NRMSD decreases with increasing borehole depth (because compaction rates are greater in those boreholes). This suggests that the greatest model uncertainty is in predicting compaction rates in the near surface where there is higher density variability due to variable surface properties and processes. These results can be used to inform future ‘coffee-can’ experiments: installing multiple instruments to shallow depths would improve our ability to quantify the spatial variability and predict the near-surface compaction rate. Fortunately, shallow installations are logistically easier.

Table 1. RMSD (m a⁻¹) (NRMSD, %) for each of the models compared to the compaction data

Our new firn-densification equation has the lowest NRMSD (~20 and $\sim \!10\%$ in the 4 and 10 m boreholes respectively, and ~7–8% in the 40, 80 and 106 m boreholes), which is expected because it was tuned to the data to which we are comparing it. In developing the equation we found that no single set of parameters could fit the data perfectly. For example, adjusting parameters to improve the fit to data from one borehole had the tradeoff of worsening the fit to data from another borehole. We note also that the equation tuning is independent of the CFM runs (the equation is tuned entirely using the field measurements and then used in the CFM to predict compaction).

Among the previously published models, GSFC has the lowest NRMSD (~20% in the 4 and 10 m boreholes and $\sim \!10\%$ in the deeper boreholes). LIG has slightly higher (30%) NRMSD in the shallow holes, but NRMSD is lower than GSFC in the deepest boreholes. In the 106 m borehole, the model misfits occur during opposite times of the year: GFSC predicts the summer compaction rate well but predicts too-fast compaction in the winter; LIG predicts the winter rate well but has too little compaction in the summer. The other three models (HL, CRO, LZ) have larger RMSDs: LZ and HL do not predict the temperature-driven seasonal cycle. These findings mirror those of Arthern and others (2010), who also compared measured compaction rates with outputs from HL and an earlier generation of LZ (Li and Zwally, 2004). CRO predicts too little compaction at greater depths and does not reproduce the seasonal cycle. This is consistent with results from Stevens and others (2020), who showed that CRO predicts compaction rates that are too slow at Summit, Greenland, possibly due to the fact that CRO was originally designed for mountain snowpacks and not deep firn (Brun and others, 1992).

5.4 Activation energy

The temperature sensitivity of firn has been debated in the literature. Most firn-densification models (including all in this study except for CRO) use an Arrhenius relationship to represent the temperature sensitivity. Herron and Langway (1980) proposed activation energies in the Arrhenius term of 10.16 and 21.4 kJ mol⁻¹ in stages 1 and 2, respectively. Arthern and others (2010), which the LIG and GSFC models are based on, found optimum activation energies of 70, 80 and 120 kJ mol⁻¹ for three different sites. For their firn-densification equation, they used two activation energies based on theoretical values representing two different processes: 60 kJ mol⁻¹ for self-diffusion of water molecules through the ice lattice, and 42.2 kJ mol⁻¹ for the grain-growth activation energy. The Arthern and others (2010) equation is formulated so that the self-diffusion activation energy is sensitive to the firn temperature and the grain-growth energy is sensitive to the mean annual site temperature. This dual-energy method highlights the fact that densification is likely an amalgam of microstructural processes, each of which likely has its own temperature sensitivity. The dominant mechanism (and thus activation energy) may change with temperature and/or microstructure, which means that the correct activation energy to use in a macroscale firn-densification equation could change based on temperature or location. The LZ model uses a temperature-dependent activation energy based on ice crystal growth and deformation measurements, which is 27.7 kJ mol⁻¹ at $-50^\circ$C. Morris and Wingham (2014) used near surface, summer densification observations in Greenland to derive a densification equation (not included in this study) and found an optimal activation energy of 110 kJ mol⁻¹.

The amplitude of the compaction-rate seasonality we measured is best replicated by the equations with activation energy near 60 kJ mol⁻¹ (LIG, GSFC and our proposed equation). HL and LZ, with their lower activation energies, show only a small compaction-rate increase during the warmer months. CRO does show a compaction-rate increase in the warmer months, but its amplitude is lower than the measurements.

We investigated the assumption that the activation energy in the densification equation (Eqn (5)) is 60 kJ mol⁻¹ by repeating our model tuning procedure (Section 5.1) with assumed activation energies of 40, 50, 55, 65 and 70 kJ mol⁻¹. Each activation energy yielded its own set of optimized parameters L,  a,  b and ρ _c for Eqn (7), which we used to run the CFM. Figure 8 shows model results compared with the observations, and Table 2 shows the NRMSD using each of the activation energies. For all of the boreholes, the seasonal amplitude predicted with Q = 70 kJ mol⁻¹ is too large, and it is too little with Q = 40 and Q = 50 kJ mol⁻¹. The NRMSD is similar for Q = 55,  60 and 65 kJ mol⁻¹, and no one of those values predicts the best fit at all borehole depths. The best fit for the 106 m borehole uses Q = 60 kJ mol⁻¹, and as such we choose it to be our optimal value.

Figure 8. Compaction rates predicted by Eqn (5) using a range of values for the Arrhenius activation energy Q compared to the compaction-rate data. Each instance of Eqn (5) with the various values of Q is tuned to have its own prefactor K(ρ).

Table 2. NRMSD (%) of compaction rates predicted by Eqn (5) using a range of values for the Arrhenius activation energy Q compared to the compaction-rate data

Each instance of Eqn (5) with the various values of Q is tuned to have its own prefactor K(ρ).

5.4.1 Cumulative compaction

It is also informative to consider the cumulative compaction for 2 years of observation. Figure 9 shows this for the 106 m borehole. The gray shaded area shows the estimated ±6% measurement uncertainty (Section 3.4). Our proposed model matches the observations to within a millimeter, because it was tuned to the data. Of the previously published models, LZ and GFSC over-predict, and HL and LIG under-predict, the cumulative compaction. The misfits are on the order of 1–2 cm, and fall within or nearly within the ±6% uncertainty bounds. CRO significantly underestimates the cumulative compaction. The effects of the temperature sensitivity can also be seen in this figure: HL and LZ are nearly linear, whereas the other models have a distinct curve. This results in a seasonality in model-data misfit.

Figure 9. Cumulative measured and modeled shortening of the 106 m borehole from February 2017 to December 2018. The gray shaded region is a ±6% uncertainty bound on the measurements. The legend includes the value at the end of observations/model runs.

6. Discussion

6.1 Firn data and model development

Firn-compaction measurements are key for developing improved firn-densification models and for validating existing models. As altimetry measurements have increased in resolution and surface meltwater fluxes have increased on the ice sheets, the need for accurate firn models and quantified estimates of their uncertainties has also increased. Many of the efforts to develop new firn-densification equations have followed the same method as Herron and Langway (1980), which leverages depth–density profiles from firn cores to tune a firn-densification equation. Doing so assumes steady state (invoking Sorge's law to convert depth to time), but these equations are used for transient model simulations. Steady-state calibration has historically been a necessity due to a lack of measurements that enable analyses of the firn's transient evolution (Ligtenberg and others, 2014). Similarly, validation of model predictions have often tested a model's ability to predict a depth–density profile, which is testing the model performance at one point in time. Continuous measurements of firn evolution, including compaction measurements like those we present in this paper, are essential to allow steady-state assumptions in firn models to be relaxed.

The spatial resolution of our data (i.e. measuring strain over meters to tens of meters of firn) required us to calculate the parcel viscosity of a thick firn layer rather than the local viscosity of a small firn layer. We assumed that this parcel viscosity is a reasonable approximation of the local viscosity, which is supported by the comparison between the parcel viscosity and steady-state-derived viscosity. However, the parcel viscosity of a meters-long portion of a firn column with density $\bar {\rho }$ is likely greater than the local viscosity of a thin layer of firn with density $\bar {\rho }$. Although this overestimation is likely small relative to the orders-of-magnitude variability in firn viscosity, it is prudent for future studies using similar methods to carefully consider the differences between the local, REV-scale viscosity (which might be measured in a lab) and the parcel viscosity of a firn layer that spans many meters and has varying properties. Field measurements of smaller and more homogeneous firn parcels can be expected to be a better approximation of the REV-scale viscosity. Despite the limitations from resolution in determining firn viscosity, our measurements are a step toward a widely applicable constitutive relation for firn.

We chose to use an empirical, macroscale framework for our proposed equation, similar to other models (e.g. Arthern and others, 2010; Ligtenberg and others, 2011; Morris and Wingham, 2014; Medley and others, 2020). Like Morris and Wingham (2014), our approach is to frame the densification equation as a constitutive relation in which the firn's strain rate is assumed to be linearly related to the stress by a viscosity. Empirical, macroscale models have the advantage that they are simply tuned to observations, whereas physically based models require knowledge of the evolution of numerous state variables (e.g. grain size, grain shape, coordination number) which may not be well constrained for polar firn. Additionally, empirical firn models have an advantage that they are computationally fast to run at ice-sheet scale, as opposed to a physically based model that may require a suite of equations to be solved dynamically to determine a number of state variables at each time step. The snow model SNOWPACK (Bartelt and Lehning, 2002; Lehning and others, 2002a, 2002b) holds promise as a physically based model for firn evolution, but its ability to simulate deeper firn remains uncertain (Keenan and others, 2021).

Conversely, empirical, macroscale models have several shortfalls. Notably, they are tuned to observations in past and/or present climates, which may not be representative of future climates (Ligtenberg and others, 2014; Keenan and others, 2021). However, this deficiency may be overcome by sampling the full range of modern glacierized climates, including ice sheets, ice caps and mountain glaciers. These climate regimes likely represent most or all firn-producing climatic conditions that would exist in a future warmer climate. Application of these models to past, colder climates with no modern analogue presents a challenge. Development of firn models for these climates likely needs to rely on lab measurements, but scaling those measurements to simulate transient firn evolution is not trivial.

Lack of descriptions of physical processes in empirical models can lead to model features that are not realistic. For example, in our densification equation (Eqn (5)), the densification rate is zero at the surface where there is no stress. In reality, grain-scale processes do densify the firn at the surface absent stress, but any model with a strain rate based solely on the stress will have this deficiency (Kingslake and others, 2022). One means of overcoming this is to develop a semi-empirical model that includes both tuning factors and descriptions of physical processes (e.g. Arthern and others, 2010). We do not attempt to make a correction for this issue, but one solution could be coupling the firn model to a surface snow model.

Our new results include an equation that predicts both the seasonal cycle and longer-term compaction rates well at USP50, which is expected because we used USP50 data to tune the model. However, the equation's specific tuning (i.e. the prefactor K) should be used with caution outside of this limited experiment. Our goal in developing the equation was not to develop a universally applicable equation (which is not possible with data from just one site), but rather to suggest a general functional form of an equation based on previous firn research and a constitutive relationship. An advantage of our approach compared to previous empirical firn models is that it uses a viscosity, which can be integrated in a constitutive-relation model. In our proposed form, the viscosity increases with age, but future research is likely to elucidate how specific firn microstructural processes affect the firn's viscosity, which could prove to be time independent. As our proposed equation is empirical, we urge caution in interpreting its form in terms of the underlying physical processes. Previous research suggests that the density and temperature dependence of our equation is robust and widely applicable, but it is possible that a universal firn constitutive relation could have an entirely different form (e.g. Meyer and others, 2020).

Future work could build a more universal firn-densification equation using firn-compaction measurements from additional sites. We acknowledge that a comprehensive dataset comprising compaction-rate and microstructure data will take many years to assemble. In the shorter term, our firn-densification equation can be evaluated against firn-core data from other geographic locations to test its applicability across various polar climates. Further, it can be tuned using other extant compaction-rate data and depth–density data to broaden its applicability.

6.2 Model-data comparison

Our modeling results show that different models will give different results over different timescales. At USP50, LZ predicts the cumulative compaction best of the previously published models. However, the low temperature sensitivities of HL, LZ and CRO fail to capture the seasonal cycle in compaction rates, rendering them inappropriate for use for corrections of altimetric measurements with sub-annual repeat intervals (e.g. ICESat-2; Markus and others, 2017; Smith and others, 2020). Additionally, incorrect temperature sensitivity leads to temporally varying model-data misfits. LIG and GSFC predict the amplitude of the seasonal compaction cycle reasonably well: they may be more appropriate to correct sub-annual surface-elevation measurements. However, our results show that modeling the seasonal cycle correctly can still lead to a long-term bias in the cumulative compaction.

The high correlations between compaction rates in pairs of boreholes indicate that the observed higher-frequency (daily to weekly) variability in the compaction rate is a real signal. None of the firn-densification models was able to simulate that variability; we hypothesize that it is a result of one or more microstructural processes that have a higher temperature sensitivity than the 60 kJ mol⁻¹ used in the macroscale models. The activation energy of 110  kJ mol⁻¹ determined by Morris and Wingham (2014) for the near-surface firn in Greenland in summer corroborates this idea, and these results underscore the need to understand the temperature sensitivity of individual microstructural processes and to add those processes to models.

There are several caveats that come with our data-model comparison that limit the conclusions that should be drawn. First, our measurements and model results are for a single site, and we cannot assess how model predictions would compare to continuous observations at other sites. South Pole has a relatively high accumulation rate considering its low temperature (i.e. that would be predicted using a Clausius–Clapeyron relationship), and the models may behave differently in different climate regimes (similar temperature with lower accumulation, e.g. Dome C, or warmer and higher accumulation, e.g. WAIS Divide). Additional field campaigns are required to identify how the models perform across a range of climate conditions.

The fact that the ‘best’ model (using RMSD as a metric) generated by our tuning process is different for the various boreholes suggests that there is spatial variability in the processes driving densification. Densification is driven at the grain scale, and grain structures are likely to vary on the ~10 m scale of our study site. For example, Proksch and others (2015) found significant density variability in the top meter of snow along a 46 m transect near Kohnen station, Antarctica. Additionally, there can be large variability in local-scale accumulation due to sastrugi and wind crusts (Frezzotti and others, 2005). This variability is likely to affect compaction rates on these local spatial scales.

In addition, some of the model/data mismatch is likely due to the initial conditions prescribed in the model runs rather than inaccuracies in how the densification equations represent physical processes. That is, we initialized the model runs using density data from a snow pit and a firn core, but these data have their own uncertainty, and errors in the initial density will propagate in the model results. We highlight the uncertainty in the initial density specifically because research has shown its importance (Fausto and others, 2018), and the density may have been different in our boreholes than in the snow pit due to spatial variability. However, uncertainty from a density initial condition cannot be avoided in model runs simulating firn evolution over entire ice sheets, which rely on initial density profiles generated by a model spin-up routine (e.g. Ligtenberg and others, 2011).

The model outputs may also be biased due to the tuning of firn models to uncertain climate forcings. For all six model runs, we forced the CFM using measured firn temperatures and a scaled accumulation rate. These forcings represent our best estimate of the actual site conditions. However, some of the firn-densification equations were developed by tuning firn-model outputs to match depth–density profiles when forced by outputs from a regional climate model (RCM). For example, LIG was developed using outputs from the RCM RACMO2 (Lenaerts and others, 2012); in this case, any biases in the RACMO outputs will manifest in the firn-densification model. That is, the RACMO-predicted climate is different than the actual conditions at our site, and it is possible that LIG-predicted compaction rates would match the measurements better if forced by RACMO. Quantifying the uncertainty due to these boundary-condition effects is outside the scope of the present work but is worthy of future study.

Within the scope our model-data comparison, the previously published firn-densification equations can predict compaction rates to at best 9% on sub-annual timescales, and our site-specific equation predicts the compaction rate to at best 7%. We interpret this as the minimum uncertainty introduced by an empirically tuned firn-densification equation into estimates of firn compaction rates. Given that our data are from a single site spanning 2 years, it is difficult to partition this uncertainty between the model initial/boundary conditions and the model formulation itself, but studies incorporating compaction-rate data from additional sites may provide further insight. On annual timescales or longer, the models can predict cumulative compaction to ~5% under the ideal forcing conditions. This estimate assumes that the observed compaction rate has zero uncertainty; adding the uncertainty in our compaction-rate measurements (6%) gives an uncertainty in the modeled cumulative compaction on annual and longer timescales of at best 8%. These results do not negate the possibility that any of the previously published equations might have smaller model-data misfit if it were forced with optimal initial and boundary conditions and/or tuned to specific site. However, firn model simulations on large spatial and temporal scales rely on accumulation histories and surface temperatures (and/or energy fluxes) from RCMs. The uncertainty in these forcings and in model tuning is likely to add to the uncertainty in modeled firn compaction on ice-sheet scales. The uncertainty in firn-compaction rate is one component of the uncertainty in modeling the change in FAC, but since changes in FAC are also affected by the accumulation rate, the uncertainty in modeled FAC change at ice-sheet scales is likely higher than we find the present modeling exercise.

7. Summary and conclusions

We measured firn compaction continuously for 2 years in 12 boreholes near the South Pole. Simultaneous measurements in holes of similar depth demonstrated that the continuous ‘coffee-can’ method can successfully be used to measure dry firn compaction to within ~5%. The measurements showed that the parcel viscosity of the firn is close to that which would be predicted using a steady-state assumption. This suggests that the firn near South Pole is near steady state. However, the steady-state viscosity masks significant variability in viscosity on seasonal and shorter timescales.

We developed a new firn-densification equation using a continuum-mechanics framework. Our equation reproduced the firn-compaction measurements favorably using a viscosity determined by the firn's age, density and a tuning parameter. Models of firn-microstructure evolution and its effect on viscosity could be added to this parameterized approach. However, at present lack of adequate microstructure data necessitates the continued use of empirical models. Model simulations of our study site with our new equation capture the compaction seasonality and amplitude, while simulations with five previously published firn-densification equations all failed to capture either the observed seasonality or amplitude. Our results indicate that an activation energy of 60 kJ mol⁻¹ in the Arrhenius term, as suggested by Arthern and others (2010), yields the best fit to the data. In general, uncertainty in the modeled firn compaction is at least 7%.