Simulations of firn processes over the Greenland and Antarctic ice sheets: 1980–2021

Conversion of altimetry-derived ice-sheet volume change to mass requires an understanding of the evolution of the combined ice and air content within the firn column. In the absence of suitable techniques to observe the changes to the firn column across the entirety of an ice sheet, the firn column processes are typically modeled. Here, we present new simulations of firn processes over the Greenland and Antarctic ice sheets (GrIS and AIS) using the Community Firn Model and atmospheric reanalysis variables for more than four decades. A data set of more than 250 measured depth–density profiles from both ice sheets provides the basis of the calibration of the dry-snow densification scheme. The resulting scheme results in a reduction in the rate of densification, relative to a commonly used semi-empirical model, through a decreased dependence on the accumulation rate, a proxy for overburden stress. The 1980–2020 modeled firn column runoff, when combined with atmospheric variables from MERRA-2, generates realistic mean integrated surface mass balance values for the Greenland (+390 Gt yr⁻¹) and Antarctic (+2612 Gt yr⁻¹) ice sheets when compared to published model-ensemble means. We find that seasonal volume changes associated with firn air content are on average approximately 2.5 times larger than those associated with mass fluxes from surface processes for the AIS and 1.5 times larger for the GrIS; however, when averaged over multiple years, ice and air-volume fluctuations within the firn column are of comparable magnitudes. Between 1996 and 2019, the Greenland Ice Sheet lost nearly 5 % of its firn air content, indicating a reduction in the total meltwater retention capability. Nearly all (94 %) of the meltwater produced over the Antarctic Ice Sheet is retained within the firn column through infiltration and refreezing.

1 Introduction

One of the most robust methods for measuring ice-sheet mass balance uses satellite altimetry (Shepherd et al., 2012, 2018) to measure changes in surface height through time and ultimately provide ice-sheet-wide volume change estimates (Helm et al., 2014; Paolo et al., 2015; Pritchard et al., 2009; Zwally et al., 2005, 2015). Interpretation of volume changes, however, requires ancillary information because there are several processes that generate height changes observable by satellite altimeters (Ligtenberg et al., 2011; Smith et al., 2020). The measured surface height change is a combination of signals, which reflect processes that involve ice or solid-earth mass change, or even no mass change at all. Even if we remove the solid-earth processes, partitioning the remaining ice-sheet-volume change to the appropriate material densities remains a challenge. Specifically, volume change due to ice dynamics represents a change at the density of ice (917 kg m⁻³), whereas surface processes (snowfall, sublimation, melt) typically (but not always) represent change under much lower densities (200–600 kg m⁻³) (Zwally et al., 2015). Additionally, the role of surface processes in observed volume change varies substantially in space and time, yet it remains largely unmeasured. Here, we present techniques that use modeling to constrain surface mass balance and firn processes over both the Greenland and Antarctic ice sheets (GrIS and AIS, respectively) for improved mass balance studies. Specifically, we provide details on a new approach to densification model calibration, an investigation of relevant spatial and temporal scales, uncertainty quantification, and a model of initial density.

In our modeling, we divide the ice sheets into two vertical layers of different material density, referred to hereinafter as the firn and ice columns. Typically extending tens of meters to over 100 m down from the surface (Ligtenberg et al., 2011), the firn column represents snow that has fallen, was subsequently buried, and is undergoing densification yet remains less dense than ice. The rate at which firn compacts varies and is dependent on its age, the weight of snow pressing down on it from above, temperature, and meltwater infiltration and refreezing. The ice column begins at a depth where material density becomes approximately constant (917 kg m⁻³) and terminates at the bed. In a constant climate, the annually averaged upward vertical velocity of the surface due to snow accumulation is perfectly balanced by ablation, compaction of the firn column, transformation to solid ice, and finally divergence of the underlying ice column (Zwally and Li, 2002), and the thicknesses of the firn and ice columns remain constant. In this scenario, the height change is zero.

The firn column is constantly evolving due to a changing climate, across all timescales, and the deviations in snow accumulation, meltwater production, and temperature from steady-state conditions drive changes in the firn layer thickness. The goal of this work is to simulate these changes in the firn column over the past 40+ years (1980–2021) using a firn densification model and atmospheric reanalysis climate forcing to determine its manifestation in altimetry-derived ice-sheet height change and the subsequent height change correction for mass balance studies.

1.1 Ice-sheet height and mass change

Changes in ice-sheet surface height reflect the integrated signal of several processes, some of which are related to ice or solid-earth mass change and others that reflect no mass change at all. Thus, we must decompose the full signal into various components in order to derive the quantity of interest; here, we are focusing on ice-mass change.

Observed height change ($\mathrm{d}h/\mathrm{d}t$) is defined as

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}h}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{h}_{\mathrm{f}}}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{i}}}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{e}}}{\mathrm{d}t}},\end{array}$$

where f, i, and e represent the component of $\mathrm{d}h/\mathrm{d}t$ resulting from changes in firn processes, solid-ice processes, and solid-earth movement, respectively. Here, $\mathrm{d}h/\mathrm{d}t$ is the surface height change; however, this value is not synonymous with actual height fluctuations of the full-ice-sheet column change ($\mathrm{d}{h}_{\mathrm{I}}/\mathrm{d}t$). Solid-earth uplift or subsidence impacts measured height changes yet reflect changes in bedrock elevation in response to current and past ice-mass changes rather than ice-thickness changes alone. This signal must be removed in order to isolate the height change due to combined firn and ice processes, $\mathrm{d}{h}_{\mathrm{I}}/\mathrm{d}t$:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}{h}_{\mathrm{I}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}h}{\mathrm{d}t}}-{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{e}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{h}_{\mathrm{f}}}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{i}}}{\mathrm{d}t}}.\end{array}$$

Height changes that manifest from solid-ice processes ($\mathrm{d}{h}_{\mathrm{i}}/\mathrm{d}t$) result from ice dynamical change over grounded ice, but over floating ice, there is an additional component due to sub-ice-shelf melt. These processes are difficult to observe or quantify; thus, we can approximate the solid-ice changes by further reworking Eq. (2) to remove the firn column height change signal ($\mathrm{d}{h}_{\mathrm{f}}/\mathrm{d}t$) from the total ice-sheet column change ($\mathrm{d}{h}_{\mathrm{I}}/\mathrm{d}t$):

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}{h}_{\mathrm{i}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{h}_{\mathrm{I}}}{\mathrm{d}t}}-{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{f}}}{\mathrm{d}t}},\end{array}$$

which provides the groundwork for determining ice-sheet mass balance. If the role of firn processes in ice-sheet height change is adequately modeled, we can isolate the contribution due to ice dynamical changes, which are easily converted to mass because the material is assumed to be solid ice (917 kg m⁻³). Ice-sheet mass balance estimates remain highly sensitive to small errors in the height change measurements and the modeled firn thickness signal.

Firn column changes, however, have a complicated relationship with mass change. Height changes due to variable rates of compaction of the firn column do not reflect a change in mass but impact the observed ice-sheet height variations through changes in volume and density. Meltwater production is more ambiguous: when it can infiltrate the firn and refreeze, there is no resulting mass change, but when infiltration is impeded and meltwater runs off, there is mass change. The effect of net snow accumulation always reflects a change in mass and can be positive or negative. As a result, the conversion between height, volume, and ultimately mass change requires understanding the material density of each component, which is neither constant in time nor space.

Rather than partition firn column changes by its individual components (see above), we divide total firn column height change into changes in the thickness of ice and the air thickness: surface mass balance (SMB) and firn air content (FAC), respectively. Specifically, we define $\mathrm{d}{h}_{\mathrm{f}}/\mathrm{d}t$ as

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}{h}_{\mathrm{f}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{h}_{\mathrm{SMB}}}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{FAC}}}{\mathrm{d}t}},\end{array}$$

where $\mathrm{d}{h}_{\mathrm{SMB}}/\mathrm{d}t$ and $\mathrm{d}{h}_{\mathrm{FAC}}/\mathrm{d}t$ represent height change fluctuations due to SMB and FAC. These components are defined below. These two components are not independent of one another: snow accumulates at the surface as a mixture of ice and air. We elect to partition firn height change into ice and air components for two reasons: (1) to better support ice-sheet altimetry studies and allow for removal of non-ice-mass change from the observed volume changes and (2) to partially isolate the firn modeling effort presented here from the reanalysis-generated surface mass balance variables used as forcing. Apart from surface runoff, the latter ensures that we take the SMB signal directly from the reanalysis model without modification, so the focus of the modeling work presented is almost entirely on $\mathrm{d}{h}_{\mathrm{FAC}}/\mathrm{d}t$; however, we do provide analysis of $\mathrm{d}{h}_{\mathrm{SMB}}/\mathrm{d}t$ for completeness.

1.1.1 Surface mass balance

The SMB is the summation of mass fluxes at the surface, including precipitation (solid and liquid), evaporation/sublimation, and runoff (Lenaerts et al., 2019). Here, we do not account for blowing snow processes that likely impact local-scale SMB; however, these processes comprise an overall small percentage of total SMB (Van Wessem et al., 2018). Specifically,

$$\begin{array}{}\text{(5)}& \mathrm{SMB}=\mathrm{Sn}+\mathrm{Ra}-\mathrm{Ev}-\mathrm{Ru},\end{array}$$

where Sn is snowfall, Ra is rainfall, Ev is evaporation/sublimation, and Ru is runoff. All are in units of meter ice equivalent (i.e.) per year.

1.1.2 Firn air content

The FAC or depth-integrated porosity represents the integrated volume of air within the entire firn column and is defined as

$$\begin{array}{}\text{(6)}& \mathrm{FAC}={\int }_{\mathrm{0}}^{z_{{\mathit{\rho }}_{\mathrm{i}}}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}-\mathit{\rho }\left(z\right)}{{\mathit{\rho }}_{\mathrm{i}}}}\mathrm{d}z,\end{array}$$

where ρ_i is the density of ice, $z_{{\mathit{\rho }}_{\mathrm{i}}}$ is the depth in meters at which the density of ice is reached, and ρ(z) is the density at a given depth. The FAC is in units of meters of air.

2 Materials and methods

We simulated firn column processes over both the GrIS and AIS using the Community Firn Model (CFM) framework (Stevens et al., 2020), forced by reanalysis climate variables. These simulations are referred to as the Goddard Space Flight Center firn densification model (GSFC-FDM, v1.2.1). First, we provide specifics relating to the CFM as well as our methodology for calibration, spin-up, and implementation. We then describe our selected climate forcing from NASA's Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) used in our simulations. Third, we discuss the differences between GSFC-FDMv1.2 and its earlier versions, v1 and v0, the latter of which was used in Smith et al. (2020) and Adusumilli et al. (2020). Finally, we provide details regarding our uncertainty assessment as well as our SMB evaluation approach.

2.1 Firn densification modeling: GSFC-FDMv1.2.1

2.1.1 The Community Firn Model

The Community Firn Model was built as a resource to the glaciology community, consisting of a modular, open-source framework for Lagrangian modeling of several firn- and firn-air-related processes (Stevens et al., 2020). The CFM allows the user to select the processes and/or physics of each simulation. The core CFM modules contain physics for firn density and temperature evolution; however, there are several modules for additional processes that the user can implement. For the GSFC-FDMv1.2.1 simulations, we use modules for grain-size evolution, meltwater percolation and refreezing, and sublimation. Grain-size evolution is simulated for testing purposes and not considered realistic. The user also has several options of firn densification physics from which to choose. Several of the models are calibrated using climate forcing from a regional climate model (RCM), atmospheric reanalysis, or even satellite-derived products, which means that any biases in these climate variables will bias the calibration coefficients in the firn densification model. Thus, it is necessary to have consistent climate forcing between the calibration and actual model runs, so we perform our own densification model calibration (Sect. 2.1.3). Finally, we use a simple bucket scheme for simulating meltwater percolation and refreezing; while the CFM contains a choice of physics of varying complexity, recent work by Verjans et al. (2019) suggests there is currently no evidence that the higher-order models perform better. Here, we use CFM v1.1.6 (Stevens et al., 2020, 2021).

2.1.2 Model spin-up

To ensure that we do not impose any unwanted transients in our simulations, we must have a sufficiently long spin-up interval during which most of the firn column is refreshed. Due to variable snow accumulation rates across the ice sheets, the time required to fully refresh the firn column can vary significantly. Thus, we impose a variable spin-up time that is dependent on the long-term mean climate. Specifically, we use the Herron and Langway (1980) densification model to approximate the depth to the bottom of the firn column (delineated at a density of 910 kg m⁻³) using the long-term reference snow accumulation, temperature, and surface density (see Sect. 2.1.5). This depth is divided by a burial rate (snowfall – sublimation – melt) to estimate the time needed to refresh the firn column for a given site: this first-order approximation of the age of the firn near the transition is an overestimate, which ensures that we refresh the entire column. Spin-up intervals typically span 300 to 7000 years in the Antarctic and 200 to 1500 years for Greenland. In regions with no net accumulation (snowfall < sublimation + melt), no spin-up is implemented. Rather, the simulations begin with a solid-ice column allowing the model to simulate seasonal snowfall, snowmelt, and runoff.

The CFM has the option to impose a dry-snow spin-up; however, this solution would build a firn column that is in dynamic equilibrium under dry conditions only. If melt were then imposed, meltwater processes would create large negative $\mathrm{d}{h}_{\mathrm{f}}/\mathrm{d}t$ and $\mathrm{d}{h}_{\mathrm{FAC}}/\mathrm{d}t$ that are not realistic. Instead, we only apply a 30-year spin-up to build a dry-firn column. We next repeat a baseline reference climate interval (RCI) time series the number of times required to match the estimated spin-up time described above. For example, if a location needs an 800-year spin-up and the RCI is 40 years, the latter is repeated 20 times. If the spin-up time required is not divisible by the RCI interval, we round up to the next integer to exceed the required spin-up time. The CFM is then run using this synthetic time series to generate a firn column that is in dynamic equilibrium with the climate under both dry and wet conditions over the RCI. Because the firn column consists of snow that has fallen years to decades to even centuries ago, its density evolution at present is still responding to atmospheric conditions from the recent and distant past. Thus, an RCI consisting of true atmospheric forcing over these longer timescales would provide the ideal RCI for firn column spin-up; however, these records do not exist over the ice sheets but rather span only a few decades, which means we must make assumptions regarding model spin-up and the RCI. First, we use a baseline RCI for GrIS of 1 January 1980–31 December 1995, which we assume is representative of a longer-term mean climate state. We then assume that the firn column is in equilibrium with the atmospheric conditions over that RCI, which imposes no changes of firn conditions over that interval, and the firn column can evolve freely beginning in 1996. The GrIS underwent a significant increase in temperatures and meltwater production after 1995 (see Sect. 3.2.1). For Antarctica, we define the baseline RCI as 1 January 1980–31 December 2019 because there were no appreciable shifts in climate during that time and to remain consistent with prior GSFC-FDM simulations. Like the GrIS, this assumption of RCI allows firn conditions to evolve in time; however, they are constrained by our steady-state assumption in our spin-up, which requires no net change in firn conditions over the entire RCI. We discuss the selection of RCI for both ice sheets in the Sect. 3.2 and explore the limitations of the approach in the Sect. 4.

2.1.3 Densification model and calibration

We use a subset of 256 published firn depth–density profiles from both the GrIS and AIS as the basis of our calibration procedure and perform a single calibration that is representative of both ice sheets. The density profile data set is described in Appendix A. The Arthern et al. (2010) dry-snow densification model provides the physical basis for our GSFC-FDMv1.2.1 simulations. Specifically, modeled dry-snow densification rates are separated into two stages during which the parcels experience different compaction processes and that are defined by the density of the parcel:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}t}}=c_{\mathrm{0}}\left({\mathit{\rho }}_{\mathrm{i}}-\mathit{\rho }\right);\mathit{\rho }\le \mathrm{550}\mathrm{kg}{\mathrm{m}}^{-}\mathrm{3},\text{(8)}& {\displaystyle \frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}t}}=c_{\mathrm{1}}\left({\mathit{\rho }}_{\mathrm{i}}-\mathit{\rho }\right);\mathit{\rho }>\mathrm{550}\mathrm{kg}{\mathrm{m}}^{-}\mathrm{3},\end{array}$$

where the densification rate coefficients for stage 1 and stage 2 (c₀ and c₁) are defined as a function of the cumulative accumulation above a given firn parcel ($\dot{b}$: defined as the mean accumulation rate in ice equivalence, m i.e. yr⁻¹, experienced since that parcel was deposited), the temperature of the parcel in kelvin, T, and the mean annual temperature, $\stackrel{\mathrm{‾}}{T}$:

$$\begin{array}{}\text{(9)}& {\displaystyle }{c}_{\mathrm{0}}=\mathrm{0.07}{\dot{b}}^{{\mathit{\alpha }}_{\mathrm{0}}}g\exp \left({\displaystyle \frac{-E_{{\mathrm{c}}_{\mathrm{0}}}}{RT}}+{\displaystyle \frac{E_{\mathrm{g}}}{R\stackrel{\mathrm{‾}}{T}}}\right),\text{(10)}& {\displaystyle }{c}_{\mathrm{1}}=\mathrm{0.03}{\dot{b}}^{{\mathit{\alpha }}_{\mathrm{1}}}g\exp \left({\displaystyle \frac{-E_{{\mathrm{c}}_{\mathrm{1}}}}{RT}}+{\displaystyle \frac{E_{\mathrm{g}}}{R\stackrel{\mathrm{‾}}{T}}}\right),\end{array}$$

where g is the gravitational acceleration constant (9.8 m s⁻²), the activation energy for lattice diffusion commonly used for ice is $E_{{\mathrm{c}}_{\mathrm{0}}}=E_{{\mathrm{c}}_{\mathrm{1}}}=\mathrm{60}$ kJ mol⁻¹ (Cuffey and Paterson, 2010), E_g is the activation energy for grain growth (42.4 kJ mol⁻¹), R is the gas constant (8.314 J K⁻¹ mol⁻¹), and the exponential dependence of overburden is ${\mathit{\alpha }}_{\mathrm{0}}={\mathit{\alpha }}_{\mathrm{1}}=\mathrm{1}$ (Arthern et al., 2010). Thus, the dry densification rate experienced by a given firn parcel varies in time and is based on ρ, $\dot{b}$, and T within a single stage of densification.

Figure 1 Modeled time-invariant initial density for the (a) Greenland and (b) Antarctic ice sheets. These results are based on MERRA-2 mean surface climate conditions. The solid black circles indicate locations that were used to train and test the initial density model as well as used in stage 1 calibration (Sect. 2.1.5), and the red open circles indicate sites used in stage 2 calibration. Note the differences in color scale for Greenland and Antarctica.

To begin the calibration procedure, we first run the model in its original form at 226 calibration sites across Greenland and Antarctica (Fig. 1). The number of model calibration runs is less than the actual number of observations (256) as some fall within the same grid cell (e.g., several observations from the vicinity of Summit, Greenland). All 256 observations are used. Unlike other calibration efforts (Kuipers Munneke et al., 2015; Li and Zwally, 2004, 2011; Ligtenberg et al., 2011), the calibration procedure presented here treats dry-firn densification from both ice sheets together, forming a single calibration parameterization, which benefits from a much wider range of climate conditions than if each ice sheet was treated individually.

The logarithm of the firn density profile with depth is approximately linear, largely for stage 2 (Herron and Langway, 1980). More discussion on the use of a logarithmic density profile is in Appendix B. For each calibration site, we compare the slopes of the logarithmic density versus depth for the two stages of densification between observations ($C_{\mathrm{0}}^{\mathrm{O}}$, $C_{\mathrm{1}}^{\mathrm{O}}$) and the equivalent model output ($C_{\mathrm{0}}^{\mathrm{M}}$, $C_{\mathrm{1}}^{\mathrm{M}}$) using the original Arthern et al. (2010) model configuration (Eqs. 7–10) forced by the RCI. Both the density measurements and model output are binned into half-meter depth increments to obtain similar sampling intervals. After binning, the slopes are estimated. Because density measurements are noisy, we determine the slopes in an iterative fashion, removing individual density measurements with residuals to the linear model larger than 3σ, recalculating the linear model, and repeating until all residuals are less than 3σ (i.e., an iterative 3σ edit). Calibration sites were not used in a given stage if (1) they did not contain more than seven data points for that stage prior to the 3σ edit, (2) they did not span more than 5 m in depth, (3) the final linear model produced a slope that was not significant (p>0.01), or (4) they encountered significant melt (mean annual surface melt exceeds 1 % of the mean annual snow accumulation). The latter ensures that we are only calibrating to dry-snow densification. Our final calibration data set contains 141 depth–density profiles spanning stage 1 and 76 spanning stage 2. Note that there are fewer profiles for stage 2 because not every density profile extends to stage 2 densities. There are a limited number of sites used in the calibration from Greenland because most of them cannot fully reflect dry-snow conditions (i.e., they do not meet the aforementioned melt criterion) (Fig. 1). The ratios (R₀, R₁) of the observed slopes ($C_{\mathrm{0}}^{\mathrm{O}}$, $C_{\mathrm{1}}^{\mathrm{O}}$) to the modeled slopes ($C_{\mathrm{0}}^{\mathrm{M}}$, $C_{\mathrm{1}}^{\mathrm{M}}$) provide the necessary correction (or calibration coefficient) for each site as described below.

Rather than develop a new physical form for calibration, we optimize two parameters within the Arthern et al. (2010) model: the exponential dependence on the mean annual accumulation rate since the parcel was deposited and the activation energy for creep. Arthern et al. (2010) found evidence that the activation energy is not well constrained for the sites investigated, suggesting that the physical processes at play under various conditions are not fully understood. Similarly, Ligtenberg et al. (2011) and Kuipers Munneke et al. (2015) found that the Arthern et al. (2010) model required additional dependence on snow accumulation to best fit observations. Thus, we elect to calibrate the parameters relating to variations in snow accumulation (α₀, α₁) and temperature ($E_{{\mathrm{c}}_{\mathrm{0}}}$, $E_{{\mathrm{c}}_{\mathrm{1}}}$) for each stage of densification. This choice of calibration parameters is also important because the climate forcing contains unknown biases, which can be partially overcome through calibration. Thus, the calibrated model presented here (along with all others) is only relevant when used with the same climate forcing (see Sect. 2.2).

We define our calibration coefficients for the two stages of densification (R₀, R₁) as a function of the mean accumulation rate ($\stackrel{\mathrm{‾}}{b}$) and temperature ($\stackrel{\mathrm{‾}}{T}$):

$$\begin{array}{}\text{(11)}& {\displaystyle }{R}_{\mathrm{0}}={\displaystyle \frac{C_{\mathrm{0}}^{\mathrm{O}}}{C_{\mathrm{0}}^{\mathrm{M}}}}={\stackrel{\mathrm{‾}}{b}}^{{\mathit{\beta }}_{\mathrm{0}}}\exp \left({\displaystyle \frac{-E_{\mathrm{0}}}{R\stackrel{\mathrm{‾}}{T_{\mathrm{0}}}}}\right),\text{(12)}& {\displaystyle }{R}_{\mathrm{1}}={\displaystyle \frac{C_{\mathrm{1}}^{\mathrm{O}}}{C_{\mathrm{1}}^{\mathrm{M}}}}={\stackrel{\mathrm{‾}}{b}}^{{\mathit{\beta }}_{\mathrm{1}}}\exp \left({\displaystyle \frac{-E_{\mathrm{1}}}{R\stackrel{\mathrm{‾}}{T_{\mathrm{1}}}}}\right).\end{array}$$

We solve for β and E using a least-squares fit regression model (intercept =0) with the climate forcing, $\stackrel{\mathrm{‾}}{b}$ and $\stackrel{\mathrm{‾}}{T}$, as predictor variables and our calibration coefficient, R, as the response variable. We force the intercept to zero to minimize overdetermination and allow the changes in the Arthern et al. (2010) functional form to be linked to a physical control (e.g., overburden, temperature) rather than a bulk bias shift. We first linearize Eqs. (11) and (12):

$$\begin{array}{}\text{(13)}& {\displaystyle }\ln \left(R_{\mathrm{0}}\right)={\mathit{\beta }}_{\mathrm{0}}\ln \left(\stackrel{\mathrm{‾}}{b}\right)-E_{\mathrm{0}}\left({\displaystyle \frac{\mathrm{1}}{R\stackrel{\mathrm{‾}}{T_{\mathrm{0}}}}}\right),\text{(14)}& {\displaystyle }\ln \left(R_{\mathrm{1}}\right)={\mathit{\beta }}_{\mathrm{1}}\ln \left(\stackrel{\mathrm{‾}}{b}\right)-E_{\mathrm{1}}\left({\displaystyle \frac{\mathrm{1}}{R\stackrel{\mathrm{‾}}{T_{\mathrm{1}}}}}\right).\end{array}$$

To generate $\stackrel{\mathrm{‾}}{b}$, we calculate the mean accumulation rate for each parcel since its deposition and take the average for each stage of densification. For stage 2, the firn column is effectively isothermal, so we substitute in Eq. (14) the mean temperature of all firn parcels with a density greater than 550 kg m⁻³, $\stackrel{\mathrm{‾}}{T_{\mathrm{1}}}$. Parcels undergoing stage 1 densification incur much larger fluctuations in temperature, especially near the surface. Prior versions of GSFC-FDM used a mean effective temperature within stage 1 to capture the non-linear relationship between temperature and compaction rates; however, that practice was abandoned after further evaluation against daily simulations suggest more refinement to the CFM is required to use an effective temperature (Appendix C). Thus, as is done for stage 2, we use the mean temperature of all firn parcels with a density less than 550 kg m⁻³, $\stackrel{\mathrm{‾}}{T_{\mathrm{0}}}$, in Eq. (13) to calibrate stage 1 densification.

Figure 2 The dry-snow densification calibration coefficients for (a) stage 1, R₀, and (b) stage 2, R₁ (Sect. 2.1.3; Eqs. 11–12, 15), over a range of mean annual temperatures and snow accumulation rates provide the background color contours. The coefficients derived for each of the calibration sites are plotted as closed circles, colored by their scale factor (i.e., calibration coefficients). The background color contours are derived directly from the calibrated R₀ and R₁ equations. The circles each represent the calibration coefficient derived directly from the depth–density profile at each calibration site. The black contour separates the region of enhanced densification (blue) from the region of reduced densification (red).

We finally iteratively solve for β₀, β₁, E₀, and E₁; however, only a single iteration was sufficient for both stages. To determine the uncertainties in our parameterization, we use the Monte Carlo method to explore the impact of uncertainties in the predictors. Specifically, we perform n=10 000 least-squares-fit regression models (intercept =0) using randomly perturbed predictors and predictands. The uncertainty in the modeled predictors ($\stackrel{\mathrm{‾}}{T_{\mathrm{0}}}$, $\stackrel{\mathrm{‾}}{T_{\mathrm{1}}}$, $\stackrel{\mathrm{‾}}{b}$, $C_{\mathrm{0}}^{\mathrm{M}}$, $C_{\mathrm{1}}^{\mathrm{M}}$) is derived from their variability in time (i.e., are randomly sampled in time), and the uncertainty in the observed predictands ($C_{\mathrm{0}}^{\mathrm{O}}$, $C_{\mathrm{1}}^{\mathrm{O}}$) is derived from the uncertainty in the logarithmic linear fit, which represents a Gaussian spread. Our final parameters are the mean of all 10 000 regression models, and their uncertainties are equal to the 2σ deviations. We find the optimal parameters for Eqs. (11) and (12) are

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& {\mathit{\beta }}_{\mathrm{0}}=-\mathrm{0.09}\pm \mathrm{0.03},{\mathit{\beta }}_{\mathrm{1}}=-\mathrm{0.356}\pm \mathrm{0.017},\\ & E_{\mathrm{0}}=-\mathrm{500}\pm \mathrm{300}\mathrm{J}{\mathrm{mol}}^{-\mathrm{1}},\\ & E_{\mathrm{1}}=-\mathrm{3130}\pm \mathrm{100}\mathrm{J}{\mathrm{mol}}^{-\mathrm{1}}.\end{array}\end{array}$$

These calibrated parameters, when plugged into Eqs. (11)–(12), provide the calibration coefficients for the two stages of densification, which scale the densification rate provided by the original Arthern et al. (2010) rate model (Fig. 2). The calibration largely finds reduced rates of densification during stage 1, especially at higher accumulation rates. For stage 2, the modeled compaction at the coldest and driest sites will increase, while compaction at sites experiencing moderate to high accumulation rates ($>\sim \mathrm{100}$ mm i.e. yr⁻¹) will decrease, largely as a function of the accumulation rate.

Combining Eqs. (9)–(12), (15), we define our densification rate coefficients as

$$\begin{array}{}\text{(16)}& {\displaystyle }{c}_{\mathrm{0}}=R_{\mathrm{0}}\mathrm{0.07}\dot{b}g\exp \left({\displaystyle \frac{-\mathrm{60}\mathrm{000}}{RT}}+{\displaystyle \frac{\mathrm{42}\mathrm{400}}{R\stackrel{\mathrm{‾}}{T}}}\right),\text{(17)}& {\displaystyle }{c}_{\mathrm{1}}=R_{\mathrm{1}}\mathrm{0.03}\dot{b}g\exp \left({\displaystyle \frac{-\mathrm{60}\mathrm{000}}{RT}}+{\displaystyle \frac{\mathrm{42}\mathrm{400}}{R\stackrel{\mathrm{‾}}{T}}}\right),\end{array}$$

which requires certain assumptions. Specifically, we assume that $\stackrel{\mathrm{‾}}{b}\approx \dot{b}$ and $\stackrel{\mathrm{‾}}{T_{\mathrm{0}}}\approx \stackrel{\mathrm{‾}}{T_{\mathrm{1}}}\approx T$. Concerning the former, because the CFM defines $\dot{b}$ as the mean accumulation rate after deposition of a parcel, $\dot{b}$ approaches $\stackrel{\mathrm{‾}}{b}$ with depth. Near the surface, $\dot{b}$ of a parcel can differ from $\stackrel{\mathrm{‾}}{b}$; however, the assumption is that, integrated across all parcels, the deviation is negligible. The same is true for $\stackrel{\mathrm{‾}}{T_{\mathrm{0}}}$ and $\stackrel{\mathrm{‾}}{T_{\mathrm{1}}}$: the firn pack reaches thermal equilibrium with depth, so the temperature of a parcel will deviate from the mean closer to the surface, but with increasing depth, T approaches $\stackrel{\mathrm{‾}}{T_{\mathrm{1}}}$. While these assumptions are valid for deeper firn, they are practical simplifications within the upper part where deviations in the integrated accumulation rate and temperature from the mean exist. The assumption is that in a column-integrated sense, the impact is minimized. Therefore, Eqs. (11)–(12), (15), (16)–(17), and the aforementioned assumptions produced newly calibrated parameters for use with Eqs. (9)–(10):

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& {\mathit{\alpha }}_{\mathrm{0}}=\mathrm{0.91},{\mathit{\alpha }}_{\mathrm{1}}=\mathrm{0.644},E_{{\mathrm{c}}_{\mathrm{0}}}=\mathrm{59}\mathrm{500}\mathrm{J}{\mathrm{mol}}^{-\mathrm{1}},\\ & E_{{\mathrm{c}}_{\mathrm{1}}}=\mathrm{56}\mathrm{870}\mathrm{J}{\mathrm{mol}}^{-\mathrm{1}}.\end{array}\end{array}$$

The dry compaction model used in the GSFC-FDMv1.2.1 simulations presented here is summarized by Eqs. (7)–(10) and (18). We note that the new parameters in Eq. (18) are similar to those developed by Verjans et al. (2020) despite substantial differences in the techniques used to complete the calibration. Model performance is discussed in Sect. 2.4.

2.1.4 Spatial domain

For Greenland, we define the ice boundary using the Greenland Mapping Project (GIMP) ice mask posted at 90 m spatial resolution (Howat et al., 2014). We identified approximately 13 200 of the 12.5 km GSFC-FDMv1.2.1 pixels as ice if any of the GIMP pixels that fell within were flagged as ice. For integrated SMB determination, we scale each pixel by the area of ice within based on the GIMP ice mask; the total ice-sheet area along with the peripheral ice not connected to the main ice sheet is 1.78×10⁶ km². The grid cells with positive net accumulation (i.e., snowfall – sublimation – meltwater >0), a condition required to build a firn column, amount to just over 9000. About 4100 grid cells did not meet the requirements to sustain a firn column, and their seasonal snowfall, melt, and runoff are simulated as described in Sect. 2.1.2.

For Antarctica, we use the drainage basins at 1 km resolution defined by Zwally et al. (2012). We identified any of its 12.5 km pixels that contain an ice-flagged pixel from Zwally et al. (2012) as ice, resulting in just over 88 300 ice-covered pixels. We assume all the pixels are 100 % ice-covered, which is equivalent in area to 13.6×10⁶ km² (grounded-ice-sheet area: 12.1×10⁶ km²). Most meet the positive net accumulation condition to sustain a firn column (87 800). To improve efficiency, we do not simulate firn column processes for each grid cell. Rather, we investigate the similarities in atmospheric forcing between neighboring pixels to eliminate redundant simulations. If a cell has a neighbor where (1) its mean annual temperature is consistent within 0.75 K, (2) the root mean square difference (RMSD) in snowfall minus sublimation is less than 10 % of the mean annual snowfall minus sublimation, (3) the RMSD in skin temperature is less than 0.25 K, and (4) the RMSD in meltwater production is less than 5 % of the mean annual meltwater production, then we do not run a simulation for that grid cell. These selection criteria reduce the number of simulations to 38 200. With these criteria, the fine spatial resolution is preserved in coastal regions where climate gradients are strong and is coarsened in the interior where correlation length scales are quite large (Fig. 3). Once the subset of simulations was complete, we linearly interpolated the runoff and FAC time series to fill the ∼88 300 ice-covered cells.

Figure 3 The Antarctic GSFC-FDM simulation locations colored by the representative size of their neighborhood. Darker colors (larger neighborhoods) with more white space (redundant simulations) indicate that the gradients in mean annual climate variables do not vary significantly over short length scales. Paler colors suggest stronger gradients with fewer redundant simulations.

2.1.5 Initial (surface) density

Because of the low accumulation rates over the ice sheets and the coarse (5 d) time resolution of our simulations, we anticipate significant reworking of the initial, low-density surface snowpack. Ideally, the imposed initial density would vary in time based on the ambient climate conditions; however, there are few studies that focus on the temporal evolution of freshly fallen snow over the ice sheets (e.g., Groot Zwaaftink et al., 2013). Thus, we focus rather on improving the bulk (or time-invariant) initial density for each grid cell based on the mean annual climate conditions as done by Helsen et al. (2008) and Kuipers Munneke et al. (2015). This approach means that, on average, we will approximate the surface density well, but we accept that there might be significant deviations from this bulk density over shorter timescales.

To build a model of initial density (ρ₀), we estimate initial densities from 233 depth–density profiles (stage 1) by finding the surface intercept of a linear fit to the logarithm of density versus depth (Fig. 1). These represent the best fit of the initial density to the observed density profile and include sites that are both dry and wet, resulting in a larger number of useable stage 1 profiles than for the dry-snow densification calibration (Sect. 2.1.3; n=141). We then trained a Gaussian process regression model to predict the observed initial densities using the mean annual MERRA-2 surface climate (snow accumulation, air temperature, total wind speed, and specific humidity), of which air temperature had the largest impact on prediction. The 233 initial densities were split into a training (n=187) and testing (n=46) partition, the latter of which provides an assessment of model performance. The model results are shown in Fig. 4: while we capture nearly 50 % of the variability within the testing partition (r²=0.46), predicted densities remain too high at the lowest densities. Specifically, the RMSE of all observations (training and testing; n=233) is 16.6 kg m⁻³, whereas the RMSE for observed initial densities less than 330 kg m⁻³ is 30.9 kg m⁻³ (n=20). The ρ₀ values used in GSFC-FDMv1.2.1 are displayed in Fig. 1. The upper and lower 5 % of the initial densities span 327–387 kg m⁻³ for the GrIS and 350–417 kg m⁻³ for the AIS with respective median values of 369 and 382 kg m⁻³. For Greenland, we find higher densities around the periphery, which is in line with other studies (Machguth et al., 2016a; Fausto et al., 2018; Kuipers Munneke et al., 2015); however, the lower end of our distribution is biased high, which will have implications on the modeled firn air content.

Figure 4 Comparison of observed and modeled initial densities. Solid circles indicate that they were not used in the final model development and represent an independent testing data set, whereas open circles represent the training partition used to build the Gaussian process regression model (see Sect. 2.1.5). The statistics in black are in reference to the solid circles only (testing partition), and those in gray are in reference to the open circles (training partition).

2.2 MERRA-2

MERRA-2 is a global atmospheric reanalysis developed at the Global Modeling and Assimilation Office (GMAO) at the NASA Goddard Space Flight Center (Gelaro et al., 2017). Atmospheric variables are provided at $\mathrm{0.625}{}^{\circ }\mathrm{longitude}\times \mathrm{0.5}{}^{\circ }$ latitude resolution and span the satellite era (1980–present). Here, we use the MERRA-2 snowfall, total precipitation, evaporation, 2 m air temperature, and skin temperature at hourly resolution and land ice runoff flux at monthly resolution covering 1 January 1980 to 30 September 2021 (GMAO, 2015a, b, c). At the $\pm \mathrm{70}{}^{\circ }$ latitude bands, the model has a resolution of 24 km×56 km, which is too coarse to resolve steep coastal topography such as the Antarctic Peninsula or the GrIS ablation zone. Thus, we rely on offline, 12.5 km “replay” MERRA-2 runs over both the GrIS and AIS to improve representation of regions of steeply sloping topography.

MERRA-2 employs the incremental analysis update (IAU) scheme of Bloom et al. (1996). The IAU uses predictor and corrector model forward integrations where differences with observations are first computed in the predictor segment and then added as an additional forcing term in the corrector run. It may be noted that an entirely different global model may employ the IAU scheme to correct to MERRA-2 innovation variables every 6 h, a process referred to as “replay” (e.g., Mapes and Bacmeister, 2012). The MERRA-2 12 km replay integration (M2R12K) was produced as part of the NASA Downscaling Project (Tian et al., 2017) and covers the period December 1999 to November 2015. A non-hydrostatic version of the Goddard Earth Observing System (GEOS) model was used in the replay integration with an output grid spacing of $\mathrm{1}/\mathrm{8}{}^{\circ }$ by $\mathrm{1}/\mathrm{8}{}^{\circ }$, but with the same vertical resolution as the original MERRA-2. The atmospheric model was modified to repartition large-scale and convective processes, and the analysis increment was filtered to allow for features of a higher resolution than resolved in the original MERRA-2 analysis grid.

Figure 5 Mean annual net accumulation (snowfall minus sublimation) for the Greenland (a, b) and Antarctic (c, d) ice sheets from MERRA-2 (a, c) and M2R12K (b, d) over their contemporaneous time span (2000–2014). Note the differences in color scale for Greenland and Antarctica.

Figure 6 Mean annual skin temperature for the Greenland (a, b) and Antarctic (c, d) ice sheets from MERRA-2 (a, c) and M2R12K (b, d) over their contemporaneous time span (2000–2014). Note the differences in color scale for Greenland and Antarctica.

The high-resolution M2R12K only spans 15 years, so it cannot be used as direct forcing of the firn densification model. Rather, we retain the seasonal magnitudes in the atmospheric variables from the M2R12K to provide hybridized MERRA-2 output. First, the MERRA-2 output is oversampled to the M2R12K grid. We then determine the 2000–2014 monthly means in MERRA-2 and remove them from the full MERRA-2 record (1980–2021). The 2000–2014 M2R12K monthly means are then added to the MERRA-2 residuals to form the hybridized MERRA-2 atmospheric variables. In such a manner, the magnitude of the gradients in precipitation and temperature from the high-resolution M2R12K are transferred to the coarse MERRA-2 output. Figures 5 and 6 show the mean annual net accumulation (snowfall minus sublimation) and skin temperature, respectively, for the GrIS and AIS. For simplicity, we hereinafter refer to the hybridized MERRA-2 as MERRA-2.

While the variables are provided at hourly resolution, to maximize computational efficiency, we perform the firn simulations at a resolution of 5 d. The 5 d MERRA-2 time series are built by averaging the hourly data over 5 d intervals. Although MERRA-2 includes meltwater processes, only net runoff is retained. Thus, we use a degree-day approach to build gridded meltwater time series, which is described in Sect. 2.2.1.

2.2.1 Degree-day model

For both ice sheets, we used a simple model to generate meltwater fluxes for input into the CFM. Specifically, meltwater production (m) was estimated using a calibrated degree-day model (e.g., van den Broeke et al., 2010):

$$\begin{array}{}\text{(19)}& m=\mathrm{DDF}\times \sum _t(T_{\mathrm{2}\mathrm{m}}-T_{\mathrm{0}})\mathrm{\Delta }t;T_{\mathrm{2}\mathrm{m}}>T_{\mathrm{0}}.\end{array}$$

Melt was activated when the 2 m air temperature (T_2 m) exceeded a calibrated temperature threshold (T₀); the exceedance is then scaled by the calibrated degree-day factor (DDF: kg m⁻² h⁻¹ K⁻¹) to generate the magnitude of melt. Here, we used hourly temperatures (Δt=1 h) to estimate 5 d (t=5 d) meltwater production. While degree-day models traditionally use Δt=1 d, we used a finer temporal resolution to ensure more realistic meltwater production, but ultimately, melt was accumulated over a 5 d window.

Figure 7 Degree-day model evaluation for (a) Antarctica and (b) Greenland for different temperature thresholds, T₀. The orange and black lines represent the mean r² and RMSE of every grid cell for a given temperature threshold. The orange and black dots indicate the thresholds that maximize r² and minimize RMSE, respectively. The dashed line represents the threshold selected for each ice sheet as it maximizes the normalized distance between the two curves and yields the best model according to these evaluators. Note the differences in vertical scales.

Figure 8 Mean degree-day factors over the Greenland Ice Sheet binned at 250 m elevation intervals for a temperature threshold, T₀, of 270.25 K. We assume that degree-day factors above the vertical dashed red line at 1500 m are non-physical (i.e., too large). For all elevations greater than 1500 m, we apply a maximum degree-day factor equal to the mean within the 1250–1500 m elevation bin (0.13 kg m⁻² h⁻¹ K⁻¹).

We calibrated our melt model for Antarctica using a calibration data set of surface meltwater fluxes (Trusel et al., 2013) that span the 1999 to 2009 melt seasons, which are linearly interpolated to our 12.5 km grid. Rather than calibrate our model to 5 d meltwater fluxes, we optimized correspondence of annual meltwater production between the model and calibration data and set t=1 year. For each grid cell, we quantified the DDF that best relates the annual accumulated exceedance of T_2 m over a predetermined threshold, T₀, which does not vary in space. To evaluate which temperature threshold yielded the best model, we calculated DDF under a wide range of T₀ (265–273 K) at quarter-degree intervals. To eliminate unrealistic DDF, we set all DDF in the upper 1 % to the 99th percentile factor. We evaluated the performance of these models in reproducing the annual time series of Antarctic-wide meltwater production as compared to our calibration data set (Trusel et al., 2013). Specifically, we compared on a grid-cell basis their ability to reproduce interannual variability (r²) and to minimize mismatch (RMSE) with the observations. Giving equal weight to the aforementioned, we found the ice-sheet-wide mean r² and RMSE for each threshold and selected the T₀ that maximizes the normalized distance between the curves of the two evaluators (r² and RMSE) in Fig. 7a. This approach selected a threshold that lies in between the threshold if determined by one evaluator alone. For Antarctica, we used T₀=270.25 K.

We estimated a temperature threshold over the GrIS using a similar approach. While we used an observation-based calibration data set over Antarctica, a similar data set does not exist for Greenland, so we instead used independent model output as the basis of our calibration. Specifically, we used the 1980–2014 annual meltwater rates from the MARv3.5.2 regional climate model (RCM) (Fettweis et al., 2017). Although this product provides sub-annual resolution, we opted to calibrate to annual meltwater production once more. In such a manner, the short-timescale meltwater fluxes were driven by MERRA-2, but the calibration to annual RCM output ensured that the simple model remains aligned with realistic annual magnitudes from MAR. For Greenland, we found a threshold, identical to the Antarctic, of T₀=270.25 K (Fig. 7b).

Figure 9 Mean annual meltwater fluxes for the (a) Greenland and (b) Antarctic ice sheets based on a degree-day approach. Note the differences in color scale for Greenland and Antarctica.

For both ice sheets, the temperature threshold is below freezing, which suggests either (1) a cold bias in MERRA-2 or (2) too strong a melt within the calibration data sets. The former has been found over Greenland (Hearty III et al., 2018) and Antarctica (Gossart et al., 2019; Huai et al., 2019) for summer months, but we cannot eliminate the latter as a contributor to the sub-freezing threshold as well, which we discuss more in Sect. 4. We assess the realism of the calibrated GrIS DDF by plotting the mean values over 250 m elevation bins. Moving into the interior, we would expect lower DDFs as the surface is typically bright snow, whereas lower elevations are more likely to exhibit bare ice and lower albedos, which would yield higher DDFs. For Greenland, the relationship between elevation and DDF exhibits high values at lower elevations, which drop off to a near-stable value around 1500 m, above which the values rapidly increase (Fig. 8). We assume that the stable values around 0.13 kg m⁻² h⁻¹ K⁻¹ are likely more representative of expected values moving upward into the dry-snow zone than the values obtained by allowing the DDFs to unrealistically rise. Above 1500 m, DDFs are capped at 0.13 kg m⁻² h⁻¹ K⁻¹ while calibrated values below that cap are untouched. The lower and upper 5 % DDF bounds over the GrIS are 0.06 and 0.21 kg m⁻² h⁻¹ K⁻¹. For Antarctica, we cannot take a similar approach as its geometry and the presence of large floating ice shelves complicate the interpretation of the relationship between the DDF and elevation. Thus, since the majority of melting occurs over the ice shelves, we use typical values from the ice shelves to limit DDF over the higher elevations of the grounded ice sheet. Specifically, we found the 95 % bounds of DDFs over the ice shelves. If a DDF is less than the lower bound, we set it to zero, and if it is larger than the upper bound, we cap it at that upper bound. The lower and upper bounds for T₀=270.25 K are 0.01 and 0.18 kg m⁻² h⁻¹ K⁻¹ with a mean of 0.06 kg m⁻² h⁻¹ K⁻¹. After capping the lower and upper bounds, the range of DDFs remains effectively the same: 0.01–0.18 kg m⁻² h⁻¹ K⁻¹ (95th percentile), but the mean goes up to 0.09 kg m⁻² h⁻¹ K⁻¹. These modified DDFs are then used to generate 5 d meltwater production using Eq. (19) with a temperature threshold of T₀=270.25 K (Fig. 9). The calibrated DDFs for Antarctica are typically smaller than those from Greenland, which is logical given that a significant portion of the GrIS has exposed bare ice during the summer months, enhancing meltwater production. The lower bound was capped for Antarctica to exclude unrealistically low melt factors; however, follow-on analysis should involve studying the impact of this assumption on the final meltwater fluxes. The meltwater production model implemented is a source of substantial uncertainty within our results; development of a surface energy balance scheme within the CFM is underway and will provide a more robust representation of meltwater fluxes in the future.

2.3 Improvement from GSFC-FDMv0 and v1

The results presented here build off prior simulations, GSFC-FDMv0 and v1, detailed in a previous publication (Smith et al., 2020; Medley et al., 2020). We have since incorporated major improvements to the GSFC-FDMv1.2.1, which we outline below. Version 1.2 is obsolete as there was a bug in the CFM that excluded time steps with net sublimation. The CFM bug was fixed for v1.2.1 runs.

1. GSFC-FDMv1.2.1 includes a spatially variable initial density (ρ₀; see Sect. 2.1.5), whereas v0 used a constant 350 kg m⁻³. While v1 also used a spatially varying initial density, the formulation was not physically realistic as it only used northward wind and not eastward winds as predictors. For v1.2.1, we also include more observations, even those where wet firn processes occur, whereas in v0 and v1, the observations used were limited to largely dry-firn conditions, which generated a larger mismatch in modeled initial densities around the periphery of the GrIS with observations.

2. GSFC-FDMv1.2.1 includes the calibration of the dry-snow/firn compaction model that limits the inclusion of observations based on the ratio of mean annual meltwater production to snowfall (see Sect. 2.1.3). The calibration approach for v0 did not discard observations based on their exposure to liquid water processes. This change in v1 and v1.2.1 should lead to an improvement in the representation of dry compaction, but we note that this calibrated dry-snow/firn compaction model is still used in regions of meltwater percolation.

3. GSFC-FDMv1.2.1 includes a more robust approach to handling mass fluxes at the surface. The CFM underwent a significant update between v0 and v1, including allowing the explicit removal of mass via sublimation and also inclusion of rainfall. For v0, sublimation was handled by aggregating the accumulation from neighboring time steps until positive, thereby still accounting for sublimation but at the cost of smoothing out the accumulation signal. Rainfall was not included in v0. For v1 and v1.2.1, mass via rainfall can now be added to the total liquid volume present and become subject to liquid water processes.

4. GSFC-FDMv1.2.1 includes an improved meltwater model. The degree-day approach for both v0 and v1 is the same; however, the v0 model was built using skin temperature, which cannot exceed 273.15 K and will not capture the large temperature deviations above freezing, especially in Greenland. For v1, we use 2 m air temperature (see Sect. 2.2.1), which is a more robust approach; however, extreme DDF values, largely in the interior, resulted in unrealistic melt rates. Thus, for v1.2.1, we capped DDFs based on realistic dry-snow values, which should improve meltwater fluxes in the nearly dry interior of the GrIS.

5. GSFC-FDMv1.2.1 includes runoff as an output. The older CFM version used for v0 did allow for melt, percolation, and refreezing but did not provide runoff as an output. Thus, we are now able to calculate surface mass balance using v1 and v1.2.1.

6. GSFC-FDMv1.2.1 includes an uncertainty analysis of the dry-snow calibration coefficients, which was not completed in v0 or v1. This exercise provides part of the basis for estimating total uncertainty in FAC and its evolution in time as well as total height and volume change.

7. GSFC-FDMv1.2.1 includes a time resolution of 5 d for both the GrIS and AIS. The prior versions (v0 and v1) ran subsets of the AIS at 5, 10, and 20 d, depending on their mean climate. Within v1.2.1, the entire AIS is run at 5 d resolution.

Figure 10 Comparison of observed and modeled slopes of the logarithmic density with depth for (a) stage 1 and (b) stage 2. Green circles reflect comparisons of observed slopes with those from the original densification model (Arthern et al., 2010), and the orange circles compare the same observations after calibration. Circles without an edge color are sites from Antarctica, and those with a black edge are from Greenland. Summary statistics are also color coded, where r² is the coefficient of determination and MAE is the mean absolute error. The dashed black line is the 1:1 line. Note the differences in scale.

2.4 Model performance

To evaluate the model improvement through our calibration procedure (Sect. 2.1.3), we evaluate the uncalibrated and calibrated model abilities to capture the slopes of the logarithmic density versus depth for both stages against the calibration data set. We found that the mean absolute error (MAE) in modeled slopes for both stages was reduced by nearly one-half after calibration, and the explained variance between observed and modeled was significantly increased in stage 2 (Fig. 10). The mean observed slope is 0.067 m⁻¹ for stage 1 and 0.030 m⁻¹ for stage 2. After calibration, the MAE in modeled slope reduces from 0.021 to 0.013 m⁻¹ for stage 1 and from 0.009 to 0.005 m⁻¹ for stage 2 (Fig. 10). The calibration relies heavily on modification to the accumulation rate (i.e., overburden) component of densification for both stages. Modification to the temperature dependence is necessary for stage 2 and of very minor importance for stage 1. For both stages, the sensitivity of densification rates to increasing accumulation is reduced, although the changes are more pronounced for stage 2. Densification due to temperature fluctuations during stage 2 is increased, especially at colder temperatures. Ligtenberg et al. (2011) and Kuipers Munneke et al. (2015) similarly found that the semi-empirical Arthern et al. (2010) model mostly overestimated the rate of densification and found an empirical link with the accumulation rate.

We would ideally prefer to perform an evaluation of modeled firn densification rates, but a substantial number of published observations are lacking. Here, we further evaluate the ability of GSFC-FDMv1.2 to reproduce the observed densities in our full data set of sites that are in both dry and wet conditions. Most of these observations were used in the calibration; however, those with significant melt were excluded (see Sect. 2.1.3). Thus, we break out our evaluation into sites exhibiting zero, moderate, and high melt rates, quantified by their ratio to net snowfall, and there are at least two observations within a stage. Specifically, these are respectively defined as 0 %, less than 10 %, and more than 10 % of the mean annual snowfall, and we evaluate the modeled mean absolute error in reproducing depth–density observations (Fig. 11). The error increases with larger melt fractions, especially for stage 1, where the impact of melt is stronger. Because most of these observations are included in the calibration, we report them as interquartile ranges and assume the upper bounds are more representative of a realistic error for each group. For stage 1, we expect density errors of 15.2 to 29.9 kg m⁻³ for dry snow/firn, 26.0 to 41.3 kg m⁻³ for moderate melt fractions, and 49.4 to 72.5 kg m⁻³ for high melt fractions. For stage 2, we expect density errors of 10.7 to 25.9 kg m⁻³ for dry snow/firn, 18.1 to 38.2 kg m⁻³ for moderate melt fractions, and 32.5 to 78.5 kg m⁻³ for high melt fractions. We note here that we assume that each observation was taken on 1 January 1980 for comparison with the model, which likely introduces additional error.

If we evaluate the bias in our model-derived density profiles for each stage, we find that with increasing melt, the modeled profiles exhibit a more positive bias (bias = model−observation). Specifically, the median stage 1 bias under melt scenarios of 0 %, less than 10 %, and more than 10 % of the mean annual snowfall are 7.4, 31.4, and 61.5 kg m⁻³, respectively. The respective biases for stage 2 are 11.2, 23.1, and 40.3 kg m⁻³. These biases suggest that the model likely underestimates the FAC, i.e., overestimates the density in regions of strong melt.

2.5 Uncertainty analysis

2.5.1 Firn air content, surface mass balance, and height change

We estimated the uncertainty in the total FAC and its variability through time through ensemble perturbation runs of the CFM at select locations over each ice sheet. Specifically, we completed principal component analysis (PCA) on the 5 d climate time series of variables of critical importance to our simulations: SMB and temperature. We then found the principal components that account for 95 % of the variability for both SMB and temperature. This selection yielded 41 principal components (PCs) for SMB and 4 for temperature for AIS and 14 and 4 for GrIS, respectively. We then correlated each individual PC time series with the equivalent time series at every grid cell over the respective ice sheet. The grid cell with the largest correlation with the PC was selected as a perturbation site. As such, we had 45 sites for AIS and 18 sites for GrIS. We used these locations because they are the most representative of the forcing time series across the entire ice sheet. PCA analysis of melt was not performed because it is determined by the temperature (Sect. 2.2.1).

Figure 11 Boxplot of the mean absolute errors in density for stages (a) 1 and (b) 2 of densification in relation to the total melt experienced. The number of observations that make up each distribution is listed above each bar. The black rectangle represents the interquartile range of the mean absolute error, while the horizontal line within represents the median. The whiskers show the maximum and minimum, with outliers shown as open circles. Melt fraction is defined as the ratio of mean annual melt to mean annual snow accumulation expressed as a percentage.

Table 1 The 1σ (standard deviation) uncertainties in various CFM parameters and atmospheric forcing provide the basis of our uncertainty analysis (Sect. 2.5). We also uniformly sample two assumptions regarding model setup: the thermal conductivity and RCI end year. The perturbation developed samples randomly from a Gaussian distribution with a mean of zero and a standard deviation provided within the table. The values used are either based on analysis within this work or based on other references provided.

^* Choice of seven parameterizations within the Community Firn Model.

For each of the calibration sites, we ran the CFM 100 times, each time applying 11 perturbations to the climate forcing variables, CFM parameters, and the reference climate interval. Each of the perturbations sampled from a Gaussian distribution with a mean of zero and a standard deviation based on observations or model performance with specifics found in the second column of Table 1. The choice of reference climate interval and the parameterization for the thermal conductivity of ice were exceptions: we assumed a uniform distribution of each of the various scenarios in Table 1, which details each perturbation, their sampling window, and any references. For each of the 100 perturbations, we sampled each of the aforementioned Gaussian distribution of uncertainties for the modeled initial density (ρ₀) and the calibration parameters (α₀, α₁, $E_{{\mathrm{c}}_{\mathrm{0}}}$, $E_{{\mathrm{c}}_{\mathrm{1}}}$). We also sampled approximated Gaussian uncertainties in the mean snow accumulation rate (Sn–Ev), rainfall (Ra), melt (Me), and skin temperature and then scaled the original MERRA-2 time series to modify the climate forcing used. Finally, we selected our choice of the parameterization of the thermal conductivity of firn from seven different models within the CFM and our choice of the end of the reference climate interval assuming a uniform distribution. Each calibration-site perturbation was sampled independently of others, resulting in 4500 and 1800 unique CFM runs for the AIS and GrIS. We note that this is a simplified approach into uncertainties in firn column evolution, considering atmospheric variables and CFM parameters are undoubtedly correlated, which we do not consider but could be implemented in future model versions with some modification.

Figure 12 Comparisons of the estimated 2σ uncertainty in FAC for both GrIS (solid circles) and AIS (open circles) and the model prediction of the 2σ uncertainty derived from Eq. (20), which only applies to sites with a liquid-to-solid ratio (LSR) between 0 and 1 (Sect. 2.5.1). The uncertainty for those sites that fall outside of those LSR bounds was defined as their standard deviation through time (Eq. 21). Performance statistics provided apply to the entire population (i.e., GrIS and AIS combined).

We assessed uncertainties by taking the standard deviation of the mean FAC for each of the 100 perturbations over the entire time series for a given site. We next used mean annual climate parameters (snow accumulation, rain, melt, and temperature) for each site (the original, non-perturbed MERRA-2 mean values) to predict the standard deviations in FAC. We broke the regression into two groups based on the ratio of the mean annual liquid water content (melt + rain) divided by the mean annual snowfall. This ratio is defined as the liquid-to-solid ratio (LSR). We created two uncertainty models in FAC: one for LSR <0 or LSR ≥1 (n=55) and another for 0≤ LSR <1 (n=8). The latter approximates the existence of a firn column where the mass of snowfall received outweighs the combination of meltwater production and rainfall. The former is indicative of conditions that are not suitable for firn development: a negative LSR suggests net sublimation (i.e., no solid accumulation) and an LSR greater than 1 reflects conditions where liquid processes outweigh the solid limiting formation of a firn column. Therefore, locations with LSR <0 or LSR ≥1 conditions experience only transient snow or firn pack, so we estimated the uncertainty in the mean FAC by simply taking the standard deviation of the FAC time series. Combining results from both AIS and GrIS where 0≤ LSR <1, we developed a linear regression model to approximate uncertainty in the mean annual FAC. We found that mean snow accumulation, $\stackrel{\mathrm{‾}}{b}$, and skin temperature, $\stackrel{\mathrm{‾}}{T}$, provide a robust prediction (Fig. 12) of the 1σ uncertainties in mean FAC, ${\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{\mathrm{FAC}}}:$

$$\begin{array}{}\text{(20)}& {\displaystyle }{\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{\mathrm{FAC}}}=\mathrm{15.1}+\mathrm{0.78}\stackrel{\mathrm{‾}}{b}-\mathrm{0.055}\stackrel{\mathrm{‾}}{T},\mathrm{0}\le \mathrm{LSR}<\mathrm{1},\text{(21)}& {\displaystyle }{\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{\mathrm{FAC}}}={\mathit{\sigma }}_{\mathrm{FAC}},\mathrm{LSR}<\mathrm{0}\text{or\hspace{0.17em} LSR}\ge \mathrm{1}.\end{array}$$

Using Eqs. (20)–(21), we estimated the 2σ uncertainty in the mean FAC ($\mathrm{2}{\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{\mathrm{FAC}}}$) for the GrIS and AIS, which yielded typical values ranging from 0.2 to 3.9 m for the GrIS and from 2.8 to 6.4 m for the AIS (lower and upper 5 % bounds). Colder temperatures and higher accumulation rates produce larger uncertainties in FAC. Melt, rainfall, and the LSR were not significant predictors of ${\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{\mathrm{FAC}}}$ where 0≤ LSR <1, so they were excluded from the prediction.

To quantify the uncertainty in FAC variability through time, we used the same set of perturbations and estimated the standard deviation in FAC change for each of the 100 perturbation runs over every 5 d time step, producing a time series of standard deviations. We then scaled the standard deviations in 5 d FAC change by dividing them by the absolute value of the mean 5 d FAC change, yielding a time series of standard deviations relative to the absolute value of the mean FAC change. Finally, we calculated the median scaled standard deviation over the entire time series to approximate the typical uncertainties in FAC change, which was done for each of the perturbation sites. We were unable to quantify a relationship between the relative error in FAC change and the mean climate forcing even when separating between sites that experience melt and those that do not. Rather, the relative uncertainty in 5 d FAC change (${\mathit{\sigma }}_{\mathrm{d}{h}_{\mathrm{FAC}}/\mathrm{d}t}$) did not largely change between sites, so we use the mean relative error for all sites:

$$\begin{array}{}\text{(22)}& {\mathit{\sigma }}_{\mathrm{d}{h}_{\mathrm{FAC}}/\mathrm{d}t}=\mathrm{0.134}\left|\mathrm{d}{h}_{\mathrm{FAC}}/\mathrm{d}t\right|,\end{array}$$

where $\mathrm{d}{h}_{\mathrm{FAC}}/\mathrm{d}t$ is the firn thickness change due to changes in FAC in units of meters per unit time. We use the results of our SMB evaluation to assess the uncertainty in total thickness change due to SMB (Sect. 3.4). We found the median absolute bias when comparing our mean annual SMB to a series of observations for each ice sheet (Sect. 2.5.2). Specifically, we found a 1σ uncertainty of 14 % and 23 % for GrIS and AIS, respectively:

$$\begin{array}{}\text{(23)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{d}{h}_{\mathrm{SMB}}/\mathrm{d}t}=\mathrm{0.14}\left|\mathrm{d}{h}_{\mathrm{SMB}}/\mathrm{d}t\right|,\mathrm{GrIS},\text{(24)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{d}{h}_{\mathrm{SMB}}/\mathrm{d}t}=\mathrm{0.23}\left|\mathrm{d}{h}_{\mathrm{SMB}}/\mathrm{d}t\right|,\mathrm{AIS},\end{array}$$

where $\mathrm{d}{h}_{\mathrm{SMB}}/\mathrm{d}t$ is the SMB-induced height change in units of meters per time. Future work would likely involve developing a more comprehensive assessment of SMB against observations to quantify SMB uncertainties. For instance, SMB over the AIS is largely biased at lower accumulation rates, so uncertainty development in the future could explore more complex relationships under different climate conditions or even explore spatial biases. All uncertainties listed in the publication are expressed as their 2σ equivalent.

2.5.2 Surface mass balance

We evaluated our SMB estimates through comparison with in situ measurements from across both ice sheets. For the AIS, we attempted to replicate the analyses as presented by Mottram et al. (2021) to ease comparisons of our performance against a suite of state-of-the-art SMB models. We used a new compilation of SMB observations from Wang et al. (2021), excluding those from Dattler et al. (2019) and Medley et al. (2013). The former study generated SMB using airborne shallow radar; however, because of the lack of an age constraint of the observed radar horizons, the layers were dated in a way to allow the derived SMB estimates to match the large-scale MERRA-2 mean. Thus, the Dattler et al. (2019) data set is dependent on the MERRA-2 SMB and is excluded. The Medley et al. (2013) data set was not used because it was excluded in the Mottram et al. (2021) evaluation, which cited the challenge in evaluating a coarsely resolved SMB data set against finely resolved radar-derived measurements. We performed a separate analysis that includes the Medley et al. (2013) data set.

After filtering the observations as described in Mottram et al. (2021) by limiting observations to the 1950–2018 interval, we arrived at a total number used in the evaluation of 16 427. We used a reference interval of 1987–2015 to match Mottram et al. (2021). For SMB observations that fall entirely within the reference interval, we compared the observation against the model mean SMB over the contemporaneous period. For the observations that cover years outside of the reference interval, we used those that span more than 5 years and compare the mean against the mean SMB over the reference interval. We also used the same aggregation approach by (1) interpolating the modeled SMB values to the location of the SMB observation and (2) averaging all the interpolated model values and observations that fall within the same grid cell. We do not do the comparison on the same common grid as Mottram et al. (2021) but rather use the 12.5 km grid used in this analysis. The final number of aggregated observations for comparison against modeled SMB was 1037 as many of the observations fall within the same grid cell (1207 if the Medley et al., 2013, data set is included).

For the GrIS, we performed a similar analysis as with the AIS using ice core observations of SMB compiled by Fettweis et al. (2020) and PROMICE (v2020) SMB observations compiled by Machguth et al. (2016b), filtering the latter to observations of greater than 3 months with a start date after 1980. We also used an ensemble mean of 13 SMB models (GrSMBMIP) to add context to the evaluation (Fettweis et al., 2020). For each observation, we linearly interpolated the model SMB to the observation location, repeating for both the GSFC and GrSMBMIP models. To minimize bias imparted by poor spatial sampling, we averaged all the observations and their associated model values into the 12.5 km grid used in this study, as done in Mottram et al. (2021). We compared the observations to the models in three ways. First, we determined the mean GSFC SMB over the exact observation interval. Second, to ease comparison with GrSMBMIP, we calculated the annual GSFC mean SMB and took the mean GSFC annual SMB of the years the observation interval covered (referred to as GSFC/ANN). Third, we perform the same as the latter with the GrSMBMIP ensemble (referred to as GrSMBMIP/ANN). The final number of aggregated observations for comparison against GSFC, GSFC/ANN, and GrSMBMIP/ANN was 312. Results from the SMB evaluation follow in Sect. 3.4 and provide the basis of our SMB uncertainty analysis in Sect. 2.5.1.

Figure 13 The mean firn air content (FAC) for the (a) Greenland Ice Sheet and (b) the Antarctic Ice Sheet over their respective reference climate intervals and (c, d) the respective 2σ uncertainty. Note the differences in color scale for Greenland and Antarctica and their uncertainties.

3 Results

3.1 Firn air content

During the RCI, the average firn air content over the GrIS was 15.7 m (the mean 2σ FAC uncertainty was 3.0 m), but it varied quite substantially in space (Fig. 13a) from 0.1 to 23.9 m (lower and upper 5 % bounds). The 2σ FAC uncertainty varied from 0.2 to 3.9 m (Fig. 13c). The peripheral ice contained less FAC with an average of 1.6 m (the mean 2σ FAC uncertainty was 0.7 m), yet, like the GrIS, there was a substantial range (0.1–12.4 m; 2σ FAC uncertainty: 0.1–2.6 m). Between 1 September 1996 and 1 September 2021, the mean loss of FAC over the GrIS was 4.8 %; however, local losses up to 100 % exist while the majority ranged between a loss of 19.1 % and a gain of 1.6 %. These change estimates were based on locations where the mean annual RCI SMB was greater than zero (i.e., a firn column exists). We note that our surface density model likely overpredicts the initial density value at the lowest density values (Fig. 4), which suggests that the model might underpredict total FAC where the modeled initial densities are the lowest (Fig. 1). We attempted to account for this bias within our uncertainty analysis by perturbing the initial density (Sect. 2.5.1).

Because of the much colder conditions, the AIS firn column contains, on average, substantially more air than the GrIS. The average FAC during the RCI for the AIS was 24.0 m (the mean 2σ FAC uncertainty was 4.7 m), which typically ranges in space between 14.2 and 36.6 m (2σ FAC uncertainty: 2.9–6.4 m) (Fig. 13b, d). Floating ice has a lower average FAC (17.0 m) than the grounded ice (24.9 m) because of higher temperatures and increased meltwater production.

3.2 Surface mass balance

The net mass flux at the surface of an ice sheet is referred to as the surface mass balance (SMB; Eq. 5) and is typically presented in units of mass per unit time. Here, we use gigatons per year (Gt yr⁻¹) to refer to area-integrated values and meters of ice equivalence per year (m i.e. yr⁻¹) for local values (i.e., grid cell). We also present total meltwater production, Me. The excess Ra + Me over Ru is retained within the firn column in either a solid or liquid state. Ru is taken directly from the CFM output and not from MERRA-2.

Figure 14 (a) The mean annual surface mass balance of the Greenland Ice Sheet and its peripheral ice over the reference climate interval (RCI: 1980–1995). (b) Time series of annual SMB (black) and its components (see Eq. 5) for the Greenland Ice Sheet only. Net accumulation (blue) and runoff (red) are combined to make SMB. Meltwater production is in yellow. The annual values are calculated from October through September and are defined by the year of the latter.

3.2.1 Greenland Ice Sheet and peripheral ice

Over the RCI (1980–1995), the mean annual SMB of the GrIS was 406±103 Gt yr⁻¹ (±1 standard deviation), which was comprised of 617±62 Gt yr⁻¹ in net accumulation (Sn + Ev), 25±5 Gt yr⁻¹ in rainfall, and 237±59 Gt yr⁻¹ in runoff (Fig. 14). Total meltwater production averaged 361±68 Gt yr⁻¹, suggesting that the firn column accommodated 39 % of all liquid water at the surface (Ra + Me). The average local SMB was 0.25 m i.e. yr⁻¹; however, it typically ranged from −0.68 to +0.93 m i.e. yr⁻¹ (lower and upper 5 % bounds) where approximately 10 % of the ice sheet by area experienced SMB <0. The largest positive SMB (+3.9 m i.e. yr⁻¹) was found in the snowfall-rich southeastern GrIS, while the largest negative SMB (−5.6 m i.e. yr⁻¹) was found along the most coastal portion of the southwestern GrIS. Such large magnitudes, however, are extremely atypical. Our choice of RCI remains an assumption, and we chose to select the MERRA-2 interval that is likely similar in state to conditions of prior decades over the GrIS. We find that for the GrIS neither SMB nor any of its components nor skin temperatures experienced a significant trend over our chosen RCI (p values >0.3; Fig. 14b). We also used a two-sample t test to evaluate whether the variables from the RCI are sampled from a population with different means than after the RCI (1996–2021). We found no significant difference in annual means for SMB and Sn + Ev between the intervals during and after the RCI; however, rainfall, meltwater production, and skin temperatures are significantly elevated post-RCI (p values <0.05). Because our spin-up involves repeating the RCI until the entire column is refreshed, our choice of RCI (1980–1995) should not generate transients associated with the initialization process in our simulation, and the firn column at the beginning of the transient simulation is in steady state with atmospheric conditions over the RCI (1980–1995), after which the firn will evolve freely in response to post-RCI conditions.

After 2003, the mean annual SMB for the GrIS was 347±118 Gt yr⁻¹, a reduction of 59 Gt yr⁻¹ as compared to the RCI. Insignificant increases in solid and liquid precipitation (21 Gt yr⁻¹) were outweighed by a strong increase in meltwater production (107 Gt yr⁻¹) and ultimately runoff (79 Gt yr⁻¹). The firn column only accommodated 37 % of liquid water present at the surface, suggesting decreased firn air storage. The ablation zone grew in area by 36 %, covering 13 % of the entire GrIS. After the major melt event of 2012 when GrIS experienced its second-lowest SMB (+133 Gt yr⁻¹) over the 40-year interval, sharp reductions in runoff coupled with above-normal net precipitation allowed the SMB to recover between 2013 and 2018. In 2019, however, the GrIS incurred its lowest annual SMB (+131 Gt yr⁻¹) due to a combination of well-below-average precipitation and well-above-average melt.

The SMB of Greenland peripheral ice was never positive over the entire 1980–2021 period with a mean of $-\mathrm{55}\pm \mathrm{22}$ and $-\mathrm{74}\pm \mathrm{25}$ Gt yr⁻¹ during the RCI and after 2003, respectively. After 2003, like the GrIS, the peripheral ice bodies experienced minimal precipitation gains (5 Gt yr⁻¹) in conjunction with moderate increases in melt and runoff (both 20 Gt yr⁻¹). Over the entire 40-year record, the firn only accommodated 18 % of all liquid water, indicating that the majority of Greenland's peripheral ice is bare ice. Local SMB over the RCI ranges from −3.6 to +1.0 m i.e. yr⁻¹ with a mean of −0.9 m i.e. yr⁻¹. As expected, 73 % of the peripheral ice experienced SMB <0.

Figure 15 (a) The mean annual surface mass balance of the Antarctic Ice Sheet (AIS; floating and grounded ice) over the reference climate interval (RCI: 1980–2019). (b) Time series of annual SMB and its components (see Eq. 5) for grounded and floating AIS. SMB (black) and net accumulation (blue) use the left axis and runoff (red) and meltwater production (yellow) use the right axis. Note that the two axes span the same range (800 Gt yr⁻¹) but are shifted in magnitude. SMB is presented as a dashed line because of overlap with net accumulation. The annual values are calculated from April through March and are defined by the year of the latter.

3.2.2 Antarctic Ice Sheet

The SMB of the AIS is nearly entirely controlled by snowfall (Fig. 15). Of the 2605±145 Gt yr⁻¹ annual mass gain over the RCI (1980–2019), net accumulation (Sn + Ev) accounts for 2605±146 Gt yr⁻¹, whereas rainfall contributed a mere 6±3 Gt yr⁻¹ and runoff removed only 6±4 Gt yr⁻¹. Meltwater production does exist (96±30 Gt yr⁻¹); however, the majority (94 %) is retained within the firn column. Local SMB is predominantly positive with a mean of +0.21 m i.e. yr⁻¹, and values commonly span +0.04 to +0.71 m i.e. yr⁻¹ (lower and upper 5 % bounds). Approximately 0.4 % of the ice sheet by area exhibited mean annual SMB <0. The maximum SMB of +6.11 m i.e. yr⁻¹ was found along the spine of the western Antarctic Peninsula, whereas the minimum of −0.38 m i.e. yr⁻¹ was found at the northwestern corner of the Ross Ice Shelf. Net snow accumulation, rainfall, runoff, and skin temperatures did not experience significant trends (p values >0.4) over the RCI; thus, we assume the full 40-year record (1980–2019) is a realistic guess regarding atmospheric conditions before 1980 in the absence of longer-term atmospheric models. Meltwater production exhibited a significant negative trend (−1.1 Gt yr⁻²; p value =0.01); however, because of the extremely small contribution relative to net accumulation (Fig. 15b) and its highly localized spatial distribution, the firn column initialized over the RCI spin-up should be in equilibrium with steady-state climate conditions.

Most mass gains over the AIS occur in the form of net accumulation over the grounded ice sheet (2148±127 Gt yr⁻¹), whereas floating ice accumulates 457±27 Gt yr⁻¹. Although not substantial, meltwater production over floating ice (66±19 Gt yr⁻¹) was on average double that over grounded ice (31±11 Gt yr⁻¹). Nearly all this meltwater is retained within the firn column as runoff averages 1 and 5 Gt yr⁻¹ for grounded and floating ice, respectively. We note that the area of grounded (12.1×10⁶ km²) ice is an order of magnitude larger than floating (1.5×10⁶ km²) ice.

3.3 Height and volume change

The combined fluctuations in SMB and FAC drive the total ice-sheet volume changes due to surface processes, yet only the former constitutes an actual mass change. We evaluate the relative contributions of mass (SMB) and air (FAC) at seasonal and multi-annual timescales. When propagating errors, we account for the variable correlation in time and space.

Figure 16 Height and volume change of the Greenland Ice Sheet. (a) Seasonal changes separated into surface mass balance (SMB; blue) and firn air content (FAC; green) components. The thin lines represent each year, and the thick, dotted line represents the mean over the RCI (1980–1996). The solid and dashed lines represent the 2012 and 2019 mass balance years, respectively, where the black lines represent FAC and the red lines represent SMB. (b) A 40-year time series of volume and height change due to SMB and FAC. (c, d) The FAC and SMB combined height and volume change. The shaded regions represent the cumulative 2σ uncertainty in the volume changes, which accounts for spatial and temporal correlation in the time series. Note the differences in scale.

3.3.1 Greenland Ice Sheet

The seasonal amplitudes of the SMB and FAC components of ice-sheet-wide volume change averaged over the RCI are 143 and 236 km³, respectively (Fig. 16a), indicating that changes in the FAC are more than 1.5× larger than SMB at sub-annual timescales. When combined, this volume change translates into an ice-sheet-wide average height change of 23 cm due to seasonal variability of surface processes. During the RCI, volume increases until May, when it typically reaches its maximum, and rapidly decreases to its minimum in August, bringing the ice sheet effectively back in balance (i.e., net zero change) as by design. After the RCI, the GrIS seasonal amplitudes of the SMB and FAC components increased, respectively, to 218 and 308 km³; however, that increase is driven largely by two extreme years in 2012 and 2019 when the GrIS lost in total 585 and 594 km³, respectively (Fig. 16c). Although FAC exhibits a larger seasonal cycle, the contribution to ice-sheet-wide volume change over longer timescales (i.e., several years) is smaller for the FAC than SMB (Fig. 16b). Between 1 September 2003 and 1 September 2021, SMB anomalies and FAC changes contributed to a decrease in GrIS volume of 1591±245 km³: 368±139 km³ due to FAC (Fig. 17) and 1224±119 km³ due to SMB.

Figure 17 Rate of height change resulting from changes in firn air content over the Greenland Ice Sheet between 1 September 2003 and 1 September 2021. Note the asymmetric color bar.

Figure 18 Height and volume change of the Antarctic Ice Sheet (grounded and floating ice). (a) Seasonal changes separated into surface mass balance (SMB; blue) and firn air content (FAC; green) components. The thin lines represent each year, and the thick, dotted line represents the mean over the RCI (1980–2019). (b) A 40-year time series of volume and height change due to SMB and FAC. (c, d) The FAC and SMB combined height and volume change. The shaded regions represent the cumulative 2σ uncertainty in the volume changes, which accounts for spatial and temporal correlation in the time series. Note the differences in scale.

3.3.2 Antarctic Ice Sheet

The seasonal amplitude of height change due to surface processes alone averages to 6 cm over the entire AIS (Fig. 18c), which is one-fourth that of the GrIS (23 cm). Due to its large area, however, the seasonal volume change amounts to 808 km³. The change in FAC is 2.5 times larger than SMB and dominates the seasonal signal, amounting to 576 km³ in seasonal change (Fig. 18a), which is larger than the seasonal signal of 340 km³ from Ligtenberg et al. (2012). While the maximum and minimum volume changes due to SMB occur in October and February, respectively, those due to FAC variability occur 1 month later (November and March). Between 31 March 2003 and 31 March 2021, the AIS has grown in volume by 1526±505 km³ from surface processes alone, of which 1050±252 km³ resulted from FAC changes (Fig. 19) and 477±248 km³ from SMB. In sum, surface processes contributed $+\mathrm{85}\pm \mathrm{28}$ km³ yr⁻¹ to the volume of the AIS since 2003, a number that is vastly overshadowed by the seasonal cycle. Because the RCI encompassed the entire 1980–2019 interval, the height and volume changes in our model experiments begin and return to zero at the end of 2019 (i.e., no height change over the entire RCI).

3.4 Surface mass balance evaluation

To contextualize the SMB values derived here from MERRA-2 and the CFM, we perform SMB evaluation against observations inspired by two recent SMB model intercomparison exercises for the AIS (Mottram et al., 2021) and GrIS (Fettweis et al., 2020).

Table 2 Modeled GrIS SMB performance statistics against N=312 observations (see Sect. 2.5.2 for description), which are plotted in Fig. 20. The final column presents a comparison (N=1 688 416) between the ensemble mean annual SMB from GrSMBMIP and the GSFC mean annual SMB interpolated onto the GrSMBMIP grid that is mapped in Fig. 21. N is the total number of observations used in the model comparison, μ is the mean of the observations, σ is the standard deviation, MB is the mean bias (GSFC minus observation), RMSE is the root mean square error, and r is the correlation coefficient.

3.4.1 Greenland Ice Sheet

The comparison between observations and modeled SMB for the GrIS indicates that our model performs similar to several of the models within the GrSMBMIP exercise. Figure 20 shows the performance of the GSFC/ANN comparison against the GrSMBMIP/ANN, and Table 2 provides the statistical comparison with the observations. We note here that the GrSMBMIP ensemble mean resolved SMB better than any individual model within the ensemble, so we expect the GSFC model to have lower performance metrics than the ensemble mean. The GSFC model reproduces observed SMB under near-equal performance as GrSMBMIP for observations with SMB $>\sim -\mathrm{2}$ kg m⁻² yr⁻¹ but experiences more spread from the observations at higher melt rates. Table 2 indicates that while the net bias of the GSFC/ANN model is comparable to the GrSMBMIP/ANN, the GSFC model experiences higher spread from the observations (RMSE =0.35 kg m⁻² yr⁻¹), which indicates partly diminished performance in capturing the spatial variability.

Figure 19 Rate of height change resulting from changes in firn air content over the Antarctic Ice Sheet between 31 March 2003 and 31 March 2021.

Figure 20 Comparison of SMB observations and the GSFC/ANN (purple) and GrSMBMIP/ANN (green) data sets for the GrIS (see Sect. 2.5.2 for the method).

Using the n=312 observation–model comparison pairs (Fig. 20 and Table 2), we approximate the uncertainty in the GSFC modeled SMB in a relative sense. Specifically, we found the absolute bias for each pair, $\mathrm{bias}=\left|\left(\text{model\hspace{0.17em}–\hspace{0.17em}observation}\right)/\mathrm{model}\right|$, and assigned an uncertainty in modeled SMB equal to the median absolute bias, which is less sensitive to outliers than the mean. The typical relative bias for GrIS is 14 %, which we employ as the 1σ uncertainty in SMB (Sect. 2.5.1; Eq. 22).

Figure 21 The difference in mean annual SMB between GSFC and the GrSMBMIP ensemble mean.

Figure 22 Comparison of annual melt from the MERRA-2 degree-day model (Sect. 2.2.1) and MARv3.5.2, the model used to train the degree-day model (Fettweis et al., 2017).

We also directly compare the GrSMBMIP ensemble mean annual SMB with our GSFC results in Fig. 21 over the common 1980–2012 interval, interpolating our model results onto the GrSMBMIP grid. The GSFC model exhibits elevated SMB over the interior relative to the GrSMBMIP ensemble mean with variable differences in sign around the periphery (i.e., exhibits positive and negative differences). The statistical summary in Table 2 suggests that integrated over the entire ice sheet, the GSFC model has a slightly higher SMB. The annual mean SMB from the GrSMBMIP of 347 Gt yr⁻¹ is smaller than the GSFC mean of 383 Gt yr⁻¹. The GrSMBMIP ensemble mean SMB trend is −7.2 Gt yr⁻², whereas our GSFC results have a slightly less negative trend of −4.5 Gt yr⁻², which falls within the entire ensemble spread (−3.1 to −12.9 Gt yr⁻²). We also compare runoff values between the GrSMBMIP ensemble mean (328 Gt yr⁻¹) and GSFC (304 Gt yr⁻¹), which suggests our runoff estimates are more muted than some models, and the GSFC trend (5.4 Gt yr⁻¹) is less positive than GrSMBMIP ensemble (8.0 Gt yr⁻¹) but still falls within the ensemble spread (4.0 to 13.4 Gt yr⁻¹). These findings suggest that GSFC SMB is on average larger than the ensemble mean because of larger snow accumulation and less runoff, which is also evidenced by the differences across the interior in Fig. 21. The difference in SMB trend from the ensemble average is largely sourced from a difference in runoff trends as snowfall exhibits no trend in both.

Finally, we compare our degree-day model annual melt rates with those used to train our model (i.e., MARv3.5.2; Fig. 22). The time series have a high correlation (r²=0.94); however, their agreement in magnitude differs between the RCI and post-RCI. The MARv3.5.2 produces a stronger increase in melt than our degree-day model. This difference could stem from multiple sources including (1) a weaker increase in temperature within the MERRA-2 model, (2) our capping of melt factors above 1500 m (see Sect. 2.2.1), and (3) our final selection of the temperature threshold. Over the contemporaneous interval (1980–2014), we find that MARv3.5.2, MERRA-2, and our GSFC GrIS runoff values average 258, 279, and 274 Gt yr⁻¹. While the GSFC melt values derived in this study are lower than the training data set, runoff values are on average larger. We note that the older GSFC-FDMv1 model, along with a newer version of MAR, showed poor performance when compared with ICESat-2-derived surface height changes in the low-melt, high-elevation portions of the ice sheet over the summer melt seasons of 2019 and 2020 (Smith et al., 2022). The same study found that our new degree-day model parameterization with reduced runoff presented here performed better than v1, which averaged 307 Gt yr⁻¹ of runoff from 1980 through 2014. Thus, recent data suggest our modifications to the degree-day model better replicated observations; however, meltwater flux and its ultimate fate is at present the largest discrepancy between the SMB and FAC models and is the largest source of uncertainty in our results.

Table 3 A breakdown of GSFC SMB performance against observations separated into different elevation bands to match Mottram et al. (2021) analysis. N is the total number of observations, L is the number used for the logarithmic analysis, μ is the mean of the observations, σ is the standard deviation, MB is the mean bias (model minus observation), RMSE is the root mean square error, r is the correlation coefficient, and rlog is the correlation coefficient of the logarithmic values. Values not in parentheses are our best attempt to match the Mottram et al. (2021) model–observation comparison method. Values within the parentheses also include SMB values from Medley et al. (2013), which were excluded from the Mottram et al. (2021) analysis.

3.4.2 Antarctic Ice Sheet

Replicating the analysis within Mottram et al. (2021) was more straightforward, so we present analysis that allows for direct comparison with their results. Figure 23 compares all SMB observations with the GSFC modeled SMB, and statistics of the evaluation are presented in Table 3, broken down into different categories as done by Mottram et al. (2021). Considering the AIS as a whole, GSFC SMB has a very small positive mean bias (6 kg m⁻² yr⁻¹) as compared to larger, negative biases from the ensemble of models in Mottram et al. (2021). Otherwise, the performance is very similar. Table 3 suggests that the GSFC SMB over ice shelves is remarkably good as compared to the Mottram et al. (2021) ensemble that suggests most models underestimate SMB. Notable differences between the GSFC SMB and the ensemble from Mottram et al. (2021) include (1) a smaller SMB bias at lower elevations than the ensemble, (2) similar performances over mid-elevations, and (3) a larger, positive bias in GSFC SMB at the highest elevations (<2800 m) where snowfall is the lowest. We observe this bias in Fig. 23 as well where the GSFC SMB values fall above the 1:1 line for the lowest observed SMB values. Thus, we find that the GSFC SMB performs well over the ice shelves and coastal grounded ice sheet but likely overestimates SMB in the dry interior (Fig. 24). As done with the GrIS, we assigned an uncertainty in modeled SMB equal to the median absolute bias between the observation pairs (n=1201), which yielded a relative uncertainty of 23 % for the AIS, providing the 1σ uncertainty in SMB within our uncertainty analysis (Sect. 2.5.1; Eq. 23).

Figure 23 Comparison of N=1201 observations of SMB with GSFC modeled SMB over the AIS. Note the logarithmic scale.

Figure 24 Difference in GSFC SMB and observations over the AIS, including the results from Medley et al. (2013) (see Sect. 2.5.2 for method).

The 1980–2010 mean annual GSFC SMB is 2620 Gt yr⁻¹, which is larger than the Mottram et al. (2021) ensemble mean of 2483 Gt yr⁻¹ but remains within the ensemble spread (2023–2752 Gt yr⁻¹). For the grounded ice sheet, the mean GSFC SMB is 2170 Gt yr⁻¹, which is also within the model spread (1743–2323 Gt yr⁻¹) and similar to the ensemble mean of 2073 Gt yr⁻¹. The fact that both the grounded and total AIS SMB are larger in the GSFC model than the ensemble average is not surprising given that (1) most of the ensemble models have a negative bias over ice shelves and (2) the GSFC model has a positive bias over the interior of the ice sheet. In fact, the GSFC SMB over the ice shelves is 450 Gt yr⁻¹ – a value that is only exceeded by two models within the ensemble of nine models. The integrated GSFC SMB did not exhibit any trends through time, which is also evident in the ensemble of models from Mottram et al. (2021). We note that we use a different grid than in Mottram et al. (2021), which could have a large impact on integrated SMB (Hansen et al., 2022).

Figure 25 Comparison of annual melt from the MERRA-2 degree-day model (Sect. 2.2.1) and Trusel et al. (2013) QuikSCAT-derived surface meltwater fluxes used to calibrate the degree-day model, as well as two regional climate models, RACMO2.3p2 (Van Wessem et al., 2018) and MARv3.6.4 (Agosta et al., 2019).

Finally, we compare our degree-day model annual melt rates with those used to train our model (i.e., Trusel et al., 2013; Fig. 25). We also compare our annual melt fluxes against two regional climate models (Van Wessem et al., 2018; Agosta et al., 2019) to provide a longer context because the QuikSCAT observations cover only a decade. By design, our degree-day model best matched the magnitude of the observations from Trusel et al. (2013). The contemporaneous (1981–2016) mean annual melt rates from our degree-day model, RACMO2.3p2, and MARv3.6.4 are 99, 107, and 83 Gt yr⁻¹. We note that the annual means are accumulated over each melt season, so the degree-day model begins in 1981, which spans 1 July 1980 to 30 June 1981. All melt fluxes are calculated in the same fashion. The degree-day model annual melt corresponds closest with RACMO2.3p2 (r²=0.72), followed by MARv3.6.4 (r²=0.56). The time series from the two RCMs show a similar correspondence (r²=0.73). These three models similarly agree on very low runoff amounts (6, 1, and 2 Gt yr⁻¹, respectively); however, the MERRA-2 land ice runoff is nearly an order of magnitude larger (68 Gt yr⁻¹). The annual runoff from our degree-day model and MERRA-2 significantly correlate in time (r²=0.65). Thus, there is a discrepancy between the firn and regional climate modeling runoff and the reanalysis-derived runoff over the Antarctic Ice Sheet. Without meltwater fluxes directly from MERRA-2, we cannot determine whether this is related to the snow model within the MERRA-2 framework or whether MERRA-2 predicts larger melt fluxes than our degree-day model, leading to more runoff.

4 Discussion and conclusion

We present simulations of GrIS and AIS firn processes using the CFM forced by MERRA-2 atmospheric reanalysis data spanning more than 40 years. Specifically, we calibrate the Arthern et al. (2010) firn densification model through modification of its dependence on overburden and temperature. The resulting model reduces the rates of densification, largely due to reduced sensitivity to increasing overburden, which is approximated by the mean accumulation rate. Modification to the temperature dependence was necessary for the second stage of densification, which is in line with other studies that found the accumulation rate as a key parameter in model calibration (Kuipers Munneke et al., 2015; Ligtenberg et al., 2011). Our calibration differs and is comparable to the approach by Verjans et al. (2020), as we derive the form of our calibration using the original form of the Arthern et al. (2010) densification equation, which provides adjusted model parameters that best fit observed depth–density profiles and the MERRA-2 climate conditions. Additionally, we calibrate the model using observations from both ice sheets, resulting in one set of adjusted parameters. It is important to note that the adjustments to the densification model parameters reflect missing physical processes as well as persistent biases within the climate forcing (e.g., if the forcing exhibited a cold bias). Thus, application of these adjustments when using a different climate forcing is not recommended. Future work will investigate use of alternative calibration equations to assess its impact on the resulting volume changes.

The surface density parameterization is also dependent on the mean annual climate conditions derived from MERRA-2, so any biases will manifest in the derived coefficients. We note that while the model does a satisfactory job of reproducing moderate to high surface densities (325–415 kg m⁻³), it appears insufficient at capturing the lowest observed densities. Thus, our model potentially overestimates the initial density, predominantly over GrIS, which leads to an underestimation of the FAC. We do perturb the initial density in our uncertainty CFM runs, so we expect our FAC errors to reflect this lack of constraint. More exploration into the density of new snow accumulations and their subsequent evolution over short timescales (hours to days) and across several locations is necessary to improve this simple density model. While new snow accumulation is often very low density, these values cannot be directly applied to the firn densification model, which models density evolution over coarse time steps (5 d) during which the snow can undergo rapid densification. Thus, the GSFC-FDMv1.2.1 does not account for sub-time-step surface density evolution and requires a bulk density representative of snowfall that has been exposed at the surface for several days. We note that the Arthern et al. (2010) densification model was not developed for densification at very low densities, so even with a more realistic fresh snow density, the model as presented here and in Arthern et al. (2010) would not adequately reproduce densification of freshly fallen snow. Future improvements in the time resolution of the simulations as well as observations of the rapid evolution of new snow accumulation should provide important future improvements to the model presented. Furthermore, the modeled surface density does not evolve in time, which is likely an oversimplification, but future work will evaluate the potential to capture seasonal initial density in future versions of the GSFC-FDM.

We next review other limitations of the work we have presented, which will be the focus of future work. The choice of running the model at 5 d time steps was a subjective choice, based on the need for computational efficiency. The firn is subject to diurnal changes in temperature and melt that our model is not capable of resolving; however, we attempt to capture much of the signal at 5 d windows through accumulating fluxes at hourly resolution such as melt and snow accumulation. In the prior simulations, we used an effective mean temperature to try to capture the non-linear impact of the large diurnal fluctuations in temperature and their resulting impact on the densification rate. We abandoned that effort (see Sect. 2.1.3 and Appendix C) given its degraded performance when compared against simulations performed at 1 d time steps. Future work preserving both the physical and effective temperature means through time will help us better understand if we can adequately capture the sub-time-step temperature impact on densification moving forward.

While we indicate that our choice of RCI was our best attempt at capturing the long-term conditions, it remains a partly subjective choice that does have an ultimate impact on our results and interpretation. The challenge for all firn modeling efforts is that the firn column was built of tens to thousands of years of snow accumulation, yet we only have a spatiotemporally complete understanding of polar climate conditions arguably since the beginning of the satellite era (1979 and onwards). Thus, we make assumptions regarding how that firn column will respond to modern conditions without knowledge of the past prevailing conditions. Studies suggest variable spatial trends in both snow accumulation rates (Medley and Thomas, 2019; Thomas et al., 2017) and air temperatures (Steig et al., 2009; Nicolas and Bromwich, 2014; Bromwich et al., 2013) over the AIS, which are not considered in this work. Similarly, the reconstructed SMB from a 20th century reanalysis found significant trends over GrIS since 1870, which this work does not capture (Hanna et al., 2011). Thus, any deviation of the RCI atmospheric conditions from reality will bias the trends in firn column evolution. Future work investigating the impact of these reconstructed trends would help to quantify the resulting uncertainty in height changes due to long-term climate change. Deviations from observed height changes thus reflect both errors in firn modeling efforts and unknown trends due to a lack of constraint on recent climate, impacting results over both ice sheets.

Meltwater fluxes as well as their ultimate fate remain the largest source of uncertainty in our firn modeling effort. Our simple degree-day model of melt was employed due to the absence of MERRA-2 meltwater flux output. At present, the CFM does not have an energy balance model subroutine, although it is in preparation, so future versions of GSFC-FDM will use a physically based melt model. Comparisons against the degree-day model training data, as well as other RCM results, suggest that we are capturing a significant portion of the annual signal (Figs. 22 and 25). The total magnitude of melt is less than the training data set for the GrIS, which might be due to (1) an overestimation of melt within the RCM used to train our model, (2) a cold bias in the MERRA-2 air temperatures, (3) capping DDFs above 1500 m, or (4) a combination of the aforementioned. Thus, the runoff produced by GSFC-FDMv1.2.1 is on the lower end of several existing SMB models for the GrIS and exhibits a smaller increase in runoff through time. We note that a recent study by Smith et al. (2022) found that the older melt model used for GSFC-FDMv1.1, as well as a more recent version of MAR than used in this study (i.e., MARv3.11.5; Amory et al., 2021), systematically overpredicted the height changes within the high-elevation parts of the ice sheet, particularly in association with melt events. After capping the unrealistic melt factors above 1500 m, the melt model in v1.2.1 yields a better match of the firn height changes with satellite altimetry. Because this comparison only covers two melt seasons, the evaluation suggests improvement, but comparison against more melt event/seasons is necessary to fully evaluate this improvement, rule out possible compensating errors, and highlight other potential future improvements.

While not the focus of the work, one important output from the GSFC-FDMv1.2.1 simulation is surface runoff, which allows us to estimate ice-sheet SMB. The mean annual GrIS SMB is comparable to SMB estimates from an ensemble of models of varying complexity (Fettweis et al., 2020). Our estimates of the 1980–2012 GrIS mean annual SMB (383±111 Gt yr⁻¹) and runoff (304±86 Gt yr⁻¹) are similar to the ensemble averages (347±111 and 328±101 Gt yr⁻¹, respectively). The lower runoff derived in this study along with slightly larger snow accumulation rates accounts for the larger SMB. Comparison with an accompanying Antarctic model ensemble suggests that our AIS SMB estimate for grounded and floating ice is larger than most models: Mottram et al. (2021) found the ensemble mean of AIS SMB of 2483 Gt yr⁻¹ (range: 2023–2752 Gt yr⁻¹), which is less than our estimate of 2620 Gt yr⁻¹. We note, however, that the evaluation in Mottram et al. (2021) of each model against observations suggests they each contain a negative bias (i.e., the modeled SMB is typically less than the observed).

Deviations in SMB from its mean over the RCI result in ice-sheet height and volume fluctuations; however, these SMB deviations along with changes in temperature also modulate the total air content within the firn column, amplifying the mass-related height and volume fluctuations. Thus, the SMB impact on height change is twofold: imposing both a change in mass and a change in air (i.e., fresh snowfall is a matrix of ice and air), which means that the fluctuations in SMB and FAC change are strongly correlated. We keep height changes due to mass separated from those due to air because of the relevance to interpretation of satellite-derived height changes. While the SMB and FAC contributions to total firn volume change over multiannual timescales are somewhat comparable, the seasonal signal is dominated by FAC for both ice sheets. This difference suggests that 62 % for the GrIS and 71 % for the AIS of sub-annual volume fluctuations are in response to a change in the air content rather than actual mass change. Thus, determination of seasonal mass change using satellite altimetry requires a substantial FAC correction, highlighting the importance of firn densification and the atmospheric models that force the FDMs, especially when investigating shorter intervals of change as not being mindful of the seasonal cycles of SMB and FAC can generate large biases.

Finally, we briefly note the differences between the GrIS GSFC results for v1.1 and v1.2.1 (differences were negligible over the AIS). The largest difference is the muted FAC change through time integrated over the GrIS (Fig. 16b). This change is partly due to the improved surface density model that yields lower densities over the interior and higher densities around the periphery, which led to a larger increase in firn air over the interior in response to additional snowfall and a smaller decrease in firn air in the percolation zone (Fig. 17). Other factors include the modification of the melt regime at high elevations, which acted to reduce total meltwater fluxes in the interior, reducing the FAC losses due to melt. Finally, the overestimation of density (or underestimation of FAC) at sites with high melt would potentially generate FAC change biased low as there is less air to lose when melt occurs. Thus, substantial FAC loss occurs along the periphery of the GrIS, but those losses are partly balanced by gains in the interior. Small changes in the surface density and liquid water processes yield measurable changes in FAC and SMB, and their uncertainty limits our ability to constrain mass balance estimates from satellite altimetry. Thus, future work constraining melt, its routing, and the initial density and their spatiotemporal evolution is necessary and should be a priority.

The time series of firn height and volume change, split into its respective SMB (ice) and FAC (air) components, provide the data necessary to isolate the ice-dynamical change from the changes observed using airborne and satellite altimeters. Future work improving the representation of the near-surface climate, initial density, and especially liquid water processes within the firn column should improve future iterations of GSFC-FDM-modeled firn volume changes. Because of the challenges in measuring firn processes in situ, future evaluations of firn densification model representation will likely rely on direct comparisons with altimetry-derived volume changes.

Appendix A: Density data

The calibration depth–density data were compiled through combination of the SUMup data sets (Koenig and Montgomery, 2018; Montgomery et al., 2018) and other compiled sources that are listed in Table A1.

Table A1 Locations and sources of the depth–density profiles used in model calibration.

Appendix B: Linear fit to the logarithmic density profile

We compared the fit statistics when making a linear fit to the logarithmic density profile versus a linear fit to the actual density profile for each stage of densification. Table B1 summarizes the results taken from the n=141 stage 1 observations and n=76 stage 2 observations (Sect. 2.1.3). For each observation, we use all the depth–density measurements in depth to calculate the corresponding RMSE and r² of the n=141 stage 1 fits and n=76 stage 2 fits. The performances are nearly identical for stage 1; however, the fit to logarithmic density profile is significantly better than using the actual density data based on a two-sample t test (p<0.01).

Table B1 Fit statistics. Bold values indicate where the fit to logarithmic density profile is significantly better than using the actual density data based on a two-sample t test (p<0.01).

Appendix C: Discontinued use of the effective mean

We tested how the model results are affected by the surface-temperature averaging scheme, which is needed to upscale the forcing data from its native 1 h resolution to the desired 5 d resolution for the CFM runs.

To do so, we performed three types of model runs. In the first, we ran the CFM with 1 d time steps, using the 1 h MERRA-2 fields (labeled in Fig. C1 as “1 d”). In the second, we ran the CFM with 5 d time steps, and the surface temperature was calculated by taking the mean temperature for each 5 d period (labeled as “15 d, mean T”). In the third, we also ran the CFM with 5 d time steps, but we calculated the 5 d “effective” mean temperature, given by

$$\begin{array}{}\text{(C1)}& T_{\mathrm{eff}}={\displaystyle \frac{-Q}{R\times \ln \left(K\right)}},\end{array}$$

with

$$\begin{array}{}\text{(C2)}& K={\displaystyle \frac{\mathrm{1}}{n}}\sum _{j=\mathrm{1}}^n{e}^{-\frac{Q}{R{T}_j}},\end{array}$$

where n is the number of days to average over (here, n=5), Q=59.5 kJ mol⁻¹ is the activation energy, R=8.314 J mol K⁻¹, and T_j represents the temperatures (K) of each day of the resampling interval. The value for Q used was based on the calibrated activation energy for the prior GSFC model v1.

We ran the CFM with the three types of runs for two different sites (South Pole and Summit, Greenland). Figure C1 shows the firn air content (FAC) change from 1980 to 2021 predicted for the two sites for each of the three model run types. Table C1 shows, for each site, the mean FAC for the entirety of each model run (mean FAC row), the change in FAC from the start of the model run to the final time step (FAC change), and the mean modeled FAC in 2020 minus the mean modeled FAC in 1980.

In both cases, the effective mean runs produce a lower total FAC than the 1 and 5 d mean runs. The FAC change using the 5 d mean setting gives a FAC change that is closer to the 1 d value, whereas the effective mean runs predict a smaller FAC change than the 1 d runs. Thus, the use of an effective mean was abandoned; however, future work on the CFM might allow for tracking of both effective mean and physical mean of the firn parcels, which might resolve these discrepancies.

Table C1 The mean FAC, change in FAC, and 2020 mean FAC minus 1980 mean FAC predicted for each of the three model run types, for each site. The 5 d mean T results are closer to the 1 d results than the 5 d effective T method.

Figure C1 The change in FAC for the duration of the model run for the South Pole (a) and Summit (b) for each of the three model run types.