Characteristics of the 1979–2020 Antarctic firn layer simulated with IMAU-FDM v1.2A

Firn simulations are essential for understanding Antarctic ice sheet mass change, as they enable us to convert satellite altimetry observed volume changes to mass changes and column thickness to ice thickness and to quantify the meltwater buffering capacity of firn. Here, we present and evaluate a simulation of the contemporary Antarctic firn layer using the updated semi-empirical IMAU Firn Densification Model (IMAU-FDM) for the period 1979–2020. We have improved previous fresh-snow density and firn compaction parameterizations and used updated atmospheric forcing. In addition, the model has been calibrated and evaluated using 112 firn core density observations across the ice sheet. We found that 62 % of the seasonal and 67 % of the decadal surface height variability are due to variations in firn air content rather than firn mass. Comparison of simulated surface elevation change with a previously published multi-mission altimetry product for the period 2003–2015 shows that performance of the updated model has improved, notably in Dronning Maud Land and Wilkes Land. However, a substantial trend difference (>10 cm yr⁻¹) remains in the Antarctic Peninsula and Ellsworth Land, mainly caused by uncertainties in the spin-up forcing. By estimating previous climatic conditions from ice core data, these trend differences can be reduced by 38 %.

1 Introduction

The Antarctic ice sheet (AIS) is the largest freshwater reservoir on Earth and has been losing mass since at least 2002 (Shepherd et al., 2018; Rignot et al., 2019), thereby contributing to around 10 % of global average sea level rise since 1993 (Oppenheimer et al., 2019). About 99 % of the AIS is covered by a firn layer (Winther et al., 2001), which represents the transition of fresh snow to glacier ice. The firn layer thickness typically ranges from 40 m in the coastal zone and can exceed 100 m in the interior (Van den Broeke, 2008) and also fluctuates in time due to changes in accumulation, melt and firn compaction. A detailed understanding of the firn layer and its variability is important for constraining current and future AIS mass change for three key reasons. Firstly, firn depth and density estimates are required to convert altimetry observed volume changes to mass changes, which remains a major source of uncertainty in mass balance studies (Morris and Wingham, 2015; Verjans et al., 2021). Secondly, firn depth and density estimates are required to convert column thickness to ice mass to calculate solid-ice discharge over the grounding line (Rignot et al., 2019). Thirdly, firn provides pore space in which nearly all (>94 %) the surface meltwater refreezes or is retained, thereby buffering surface meltwater and preventing mass loss by runoff (Medley et al., 2022).

From both observational and modelling studies, firn layer depth and density are known to vary greatly across the AIS, resulting from the wide range of (near)-surface climatic conditions (e.g. Keenan et al., 2021), i.e. temperature, wind speed, accumulation and melt. Low temperatures in the interior (Fig. 1) cause slow densification rates, resulting in a firn layer thickness (defined here as the depth of the 830 kg m⁻³ density level) greater than 100 m. On the other hand, on warm and dry ice shelves, densification and melt are enhanced, resulting in denser firn layers with a depth less than 50 m (Van den Broeke, 2008).

Figure 1 Map of Antarctica showing annual average surface temperature from the Regional Atmospheric Climate Model (RACMO2.3p2; 1979–2020) and the locations of the in situ density observations. The yellow circles indicate observational locations of the depths of the critical density level of ρ=550 kg m⁻³ (z₅₅₀), the white circles indicate observational locations of the depths of the critical density levels of both ρ=550 kg m⁻³ (z₅₅₀) and ρ=830 kg m⁻³ (z₈₃₀), the purple circles indicate observational locations of the depths of the critical density levels of both ρ=550 kg m⁻³ (z₅₅₀) and ρ=830 kg m⁻³ (z₈₃₀) and of the firn air content (FAC), and the black crosses indicate observational locations of the surface snow densities (top 0.5 m of the firn column). The numbers indicate regions mentioned throughout the paper.

Changes of AIS surface elevation are an expression of changes in the firn layer and of dynamical change in the underlying ice and bedrock. The difficulty is that these seasonal cycles and decadal changes in elevation are not solely caused by fluctuations in mass but also largely by fluctuations in the firn density (Arthern and Wingham, 1998; Ligtenberg et al., 2012; Medley et al., 2022). In altimetry studies, it is therefore important to accurately separate these changes in mass- and density-change components. The mass component can subsequently be separated into ice-dynamical, bedrock and surface mass balance components (Willen et al., 2021).

Another important characteristic of the firn layer is its capability to retain meltwater in pore spaces by capillary forces and refreezing. In a warmer future climate, enhanced firn compaction, melt and refreezing can all potentially lead to firn air depletion, thereby limiting the meltwater storage capacity of the firn. This is especially important over ice shelves, where meltwater accumulation can lead to hydrofracturing-induced ice shelf collapse (Kuipers Munneke et al., 2014; Datta et al., 2019), thereby accelerating future AIS mass loss and sea level rise (Gilbert and Kittel, 2021). Accurate information about firn conditions is thus essential both to understand the current AIS mass balance and to predict its future.

Firn models can fill spatial and temporal gaps between the relatively sparse observations from firn cores and snow pits. While remote sensing products provide high temporal and spatial resolution information about surface height changes, and indirectly, firn properties of the top few metres (e.g. Davis and Poznyak, 1993; Medley et al., 2015; Schröder et al., 2019a; Shepherd et al., 2019), they do not provide information about the deeper firn column. Firn models are, therefore, also important to obtain vertical variations in density and can be used to perform future simulations of firn depth and density in response to climate change scenarios (Ligtenberg et al., 2012; Kuipers Munneke et al., 2014).

Firn models can roughly be divided in two classes: physics-based and semi-empirical models. Physics-based models describe densification using a constitutive relationship between stress and strain for snow, while semi-empirical models use physics-based densification equations in combination with parameters that are calibrated with observed density profiles. Semi-empirical models neglect some processes such as drifting-snow compaction (e.g. Sommer et al., 2018) and are calibrated with observations representative of the past or present climate, but these conditions might not be representative of a future climate. Nevertheless, semi-empirical models require a smaller number of poorly constrained parameters, such as roughness length and snow grain shape and size, and computational demands are lower (Keenan et al., 2021). This is especially advantageous for data-scarce regions and large-scale studies. Therefore, semi-empirical firn models have been used widely for correction of satellite altimetry products (e.g. Adusumilli et al., 2018; Smith et al., 2020; Willen et al., 2021) and to assess the meltwater buffering capacity of the Antarctic firn layer (Ligtenberg et al., 2014; Kuipers Munneke et al., 2014). Several studies have found that the performance of semi-empirical models is comparable to that of physics-based models (Steger et al., 2017; Vandecrux et al., 2020; Keenan et al., 2021).

Differences in the formulations of e.g. firn densification and fresh-snow density can lead to a substantial spread in modelled firn thickness and air content (Lundin et al., 2017; Verjans et al., 2021). Using statistical emulation of firn densification models, Verjans et al. (2021) quantified how different sources of uncertainty contribute to the total uncertainty in modelled firn thickness change of the East Antarctic ice sheet. They show that, at the basin scale, the ensemble uncertainty in firn thickness changes ranges from 0.2 to 1.0 cm yr⁻¹, which amounts to a 15 %–300 % relative uncertainty. While differences in climate forcing have the largest influence on the spread, in basins with high snowfall and a large spatial variability of climatic conditions, firn densification formulation also has a substantial contribution to the spread (up to 46 %). In basins with recent increases in snowfall rates the contribution of fresh-snow density to the total uncertainty is also substantial (up to 28 %).

In this study, we have improved the semi-empirical IMAU Firn Densification Model (IMAU-FDM) (Ligtenberg et al., 2011) by updating (i) the fresh-snow density and (ii) the firn densification parameterizations and (iii) by using updated atmospheric forcing fields. In addition, the model has been calibrated and evaluated using 112 in situ firn core density measurements across the AIS. This version, 1.2A, of IMAU-FDM supersedes the previous version (Ligtenberg et al., 2011). Here, we present the simulated contemporary characteristics of the Antarctic firn layer for the period 1979–2020. We focus on spatial, seasonal and decadal variability in firn thickness and firn air content (FAC). To evaluate the decadal surface elevation changes and seasonal cycles we compare the simulations to a multi-mission remote sensing surface height change product for the period 1992–2015 (Schröder et al., 2019a). Finally, we test how sensitive the modelled firn layer is to uncertainties in model variables and to the spin-up procedure.

2 Materials and methods

2.1 IMAU-FDM

IMAU-FDM is a semi-empirical 1D firn densification model that simulates the time evolution of firn density, temperature, liquid water content and surface height changes due to firn and surface mass balance (SMB) processes. The model employs up to 3000 layers of 3 to 15 cm thickness, which represent the firn properties in a Lagrangian way. Abbreviations of all versions of IMAU-FDM used in this study are listed in Table 1. The model was originally developed by Helsen et al. (2008) and updated by Ligtenberg et al. (2011) to a version abbreviated here to FDM v1.1. Brils et al. (2022) recently improved IMAU-FDM applied to the Greenland ice sheet (GrIS) by including a refined parameterization of the thermal conductivity to version FDM v1.2G. The various versions of IMAU-FDM have been extensively evaluated against firn density and temperature observations from both Greenland and Antarctica (Ligtenberg et al., 2011; Kuipers Munneke et al., 2015; Brils et al., 2022). We further improved the model to FDM v1.2A by updating (i) the fresh-snow density and (ii) the firn densification parameterizations, as discussed below. The model is applied to the AIS, and we use this term to broadly include the East Antarctic, West Antarctic and Antarctic Peninsula ice sheets, as well as the ice shelves and disconnected glaciers and ice caps on the continent.

Kaspers et al. (2004)Lenaerts et al. (2012)Kaspers et al. (2004)Kaspers et al. (2004)Fausto et al. (2018)

Table 1 Abbreviations and characteristics of IMAU-FDM versions used in this study.

2.1.1 Fresh-snow density

The fresh-snow density (ρ_s; kg m⁻³) is an important but poorly constrained upper boundary condition of firn models. FDM v1.1 used a fresh-snow density parameterization based on Kaspers et al. (2004) with a correction by Helsen et al. (2008), which has been calibrated over Antarctica and yields density values that typically represent the first 0.5 m of the snowpack. It varies as a function of annual average surface temperature (${\overline{T}}_{\mathrm{s}}$; K), 10 m wind speed (${\overline{V}}_{\mathrm{10}}$; m s⁻¹) and accumulation ($\dot{b}$; $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$):

$$\begin{array}{}\text{(1)}& {\mathit{\rho }}_{\mathrm{s}}=A+B{\overline{T}}_{\mathrm{s}}+C{\overline{V}}_{\mathrm{10}}+E\dot{b},\end{array}$$

in which A, B, C and E are fit coefficients. Snow crystal type and degree of riming are indeed temperature dependent, and previous work shows that fresh-snow density increases with increasing temperature (Judson and Doesken, 2000). Also, strong winds during snowfall cause crystal breaking, thereby reducing the snow effective grain size (Sato et al., 2008), which results in more efficient packing and increased fresh-snow densities. However, at least two concerns can be raised against Eq. (1). First, due to the partly similar spatial patterns of accumulation and temperature in Antarctica, a statistical relationship between fresh-snow density and accumulation may exist; however this relationship has no obvious physical grounds. Other studies have not found a dependency between annual accumulation and fresh-snow density in Antarctica, Greenland and the Alps (Lehning et al., 2002; Fausto et al., 2018; Lenaerts et al., 2012). Second, Eq. (1) neglects the impact of the meteorological conditions at the time of deposition, which can be highly variable in time. Therefore, we tested the more recent parameterization of Lenaerts et al. (2012), which has been calibrated over Antarctica, in which ρ_s is a function of instantaneous surface temperature (T_s; K) and 10 m wind speed (V₁₀; m s⁻¹):

$$\begin{array}{}\text{(2)}& {\mathit{\rho }}_{\mathrm{s}}=A+B{T}_{\mathrm{s}}+C{V}_{\mathrm{10}},\end{array}$$

in which A, B and C are fit coefficients. The resulting density values are representative of the fresh-snow-layer density, which typically represents the upper few centimetres of the snowpack. Within days after deposition, the crystal structure of freshly deposited snow breaks down due to wind and redistribution of drifting snow, which increases the surface snow density (Groot Zwaaftink et al., 2013). IMAU-FDM does not calculate this densification by wind packing and redistribution of drifting snow. The fit coefficients A, B and C in Eq. (2) for the fresh-snow density are recalibrated to improve the fit of the simulated with observed surface snow densities, defined as the top 0.5 m of the firn column, as this best matches the thickness of the sampled layer. These fit coefficients are then used to simulate the fresh-snow density, i.e. the density of the fresh snow that is added on top of the firn column. This recalibration is described in Sect. 3.1. As we only use 10 additional measurements in the evaluation dataset compared to the calibration dataset, we also perform a 10-fold cross evaluation.

2.1.2 Dry-snow densification rate

Over time, the freshly deposited snow is buried and becomes denser. The rate of the dry-firn densification ($\mathrm{d}\mathit{\rho }/\mathrm{d}t$) is calculated using the semi-empirical equation of Arthern et al. (2010):

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\text{d}\mathit{\rho }}{\text{d}t}}=D\dot{b}g({\mathit{\rho }}_{\mathrm{i}}-\mathit{\rho })e^{(-\frac{E_{\mathrm{c}}}{RT}+\frac{E_{\mathrm{g}}}{R{T}_{\text{ave}}})},\end{array}$$

where D is a constant, $\dot{b}$ is the long-term average accumulation rate ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), g is the gravitational acceleration, ρ_i is the density of bubble-free ice (917 kg m⁻³), R is the gas constant, ρ is the layer density (kg m⁻³), T is the instantaneous layer temperature (K), T_ave is the long-term average surface skin temperature (K), and E_c and E_g are the activation energies for creep (60 kJ mol⁻¹) and grain growth (42.4 kJ mol⁻¹), respectively. $\dot{b}$ is used as a proxy of the overburden pressure, thereby approximating the accumulation rate as constant. Using FDM v1.2G for Greenland, Brils et al. (2022) show that assuming a constant $\dot{b}$ introduces a minor error in the load experienced by a layer of firn (e.g. <3.2 % at Summit and <1.9 % at Dye-2). The constant D has different values above (0.03) and below (0.07) the critical density level ρ=550 kg m⁻³ to represent the distinct densification mechanisms: for ρ<550 kg m⁻³, densification mainly occurs by settling and sliding of grains, and for ρ>550 kg m⁻³, it mainly occurs by deformation, recrystallization and molecular diffusion (Herron and Langway, 1980).

By comparing simulated and observed depths of the 550 and 830 kg m⁻³ density levels (z₅₅₀ and z₈₃₀, respectively), Ligtenberg et al. (2011) found that Eq. (3) requires dimensionless correction terms MO₅₅₀ for ρ<550 kg m⁻³ and MO_830* for $\mathrm{550}<\mathit{\rho }<\mathrm{830}$ kg m⁻³, which are defined as the ratio of modelled and observed values of z₅₅₀ and z_830*, where $z_{\mathrm{830}*}=z_{\mathrm{830}}-z_{\mathrm{550}}$. The correction terms MO₅₅₀ and MO_830* are added as a multiplier to Eq. (3); MO values below 1 reduce the densification rate, and values above 1 enhance the densification rate. MO₅₅₀ and MO_830* are chosen as functions of the long-term mean accumulation rate. Ligtenberg et al. (2011) and Brils et al. (2022) used logarithmic correction functions, thus

$$\begin{array}{}\text{(4)}& \text{MO}=\mathit{\alpha }-\mathit{\beta }\ln \left(\dot{b}\right),\end{array}$$

in which MO is either MO₅₅₀ or MO_830*, and α and β are fit coefficients which differ for MO₅₅₀ and MO_830*. In FDM v1.1 a minimum value of 0.25 is imposed on MO₅₅₀ and MO_830*. For Greenland and Antarctica, different values for α and β have been used. Here, we also tested a power-law function for MO_830*:

$$\begin{array}{}\text{(5)}& {\text{MO}}_{\mathrm{830}*}=\mathit{\delta }{\dot{b}}^{-\mathit{ϵ}}+\mathit{\varphi },\end{array}$$

in which δ, ϵ and ϕ are fit coefficients. The fit coefficients in Eqs. (4) and (5) are calibrated, which is described in Sect. 3.2. As we only use 10 additional measurements in the evaluation dataset compared to the calibration dataset, we also perform a 10-fold cross evaluation.

2.1.3 Heat conduction, meltwater percolation and refreezing

The conduction of heat is simulated by using a 1D heat transfer equation, which couples vertical heat conduction to temperature gradients through the thermal conductivity of firn. In FDM v1.1 the thermal conductivity is a function of firn density only (Anderson, 1976), and in FDM v1.2G this has been extended to a function of density and temperature (Calonne et al., 2019), which is also included in FDM v1.2A. The lower boundary condition assumes a constant heat flux across the lowest layer, where the deep temperature is allowed to change along with long-term changes in surface temperature and/or internal heat release. Meltwater percolation is simulated using the bucket method, whereby each firn layer has a maximum irreducible water content that decreases with increasing density (Coléou et al., 1999). The meltwater can percolate through all layers in a single time step. The meltwater (partly) refreezes when it reaches a layer with a temperature below the freezing point. To capture the percolation and refreezing of meltwater more accurately, the model time step is reduced from 3600 to 300 s for locations with melt. Further details and evaluation of heat conduction, percolation and refreezing in IMAU-FDM are provided by Brils et al. (2022).

2.2 Model strategy

An initial firn layer is obtained by spinning up the model over a reference period until the firn layer is in equilibrium with the surface climate. The equilibrium is roughly reached when the entire firn column, defined here as the depth where the density ρ reaches 830 kg m⁻³, is fully refreshed by accumulation. The required spin-up time in years is therefore calculated as the total mass of the modelled firn layer divided by the annual accumulation in w.e. (Kuipers Munneke et al., 2015). As no obvious strong long-term trends have been detected in Antarctica's surface climate and SMB during the period 1979–2020 (e.g. Favier et al., 2017; Mottram et al., 2021) and no reliable surface climate fields are available prior to 1979, the spin-up forcing is obtained by looping over the 1979–2020 forcing data. Thereby, we make the rather strong assumption that this 42-year period is representative of the past 100–1000 years and that the firn layer is in equilibrium with that climatic period. In the actual simulation following the spin-up, a minor trend (<0.6 mm yr⁻¹ averaged over the AIS) in total FAC remains, because at the bottom of the column, ice with a density between 830–917 kg m⁻³ slowly replaces bubble-free ice (ρ=917 kg m⁻³), which originates from the initialization of the firn column prior to the spin-up. This trend is removed before further analysis of the results.

2.3 Surface elevation change

IMAU-FDM quantifies the impact of firn and SMB processes on the depth of the firn layer. The resulting vertical velocity of the firn top surface (v_tot) is defined as the sum of these components:

$$\begin{array}{}\text{(6)}& v_{\text{tot}}=v_{\text{snow}}+v_{\text{sub}}+v_{\text{snd}}+v_{\text{melt}}+v_{\text{fc}}+v_{\text{by}}+v_{\text{ice}},\end{array}$$

where v_snow represents the vertical velocity component as a result of snowfall, v_sub of sublimation/riming, v_snd of snowdrift erosion/deposition, v_melt of snowmelt and v_fc of firn compaction. v_by represents the vertical motion associated with the changing ice shelf draft when the mass of the firn layer changes, and v_ice represents the downward/upward movement of the surface by the local divergence/convergence of the ice flow, driven by the long-term vertical mass flux through the lower boundary of the firn column. In a steady-state firn layer, this equals the mass flux at the upper boundary. v_ice is therefore equal to the long-term average annual SMB ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) over the reference period divided by ρ_i but of opposite sign. All components are derived from the run-time output of IMAU-FDM. The minor residual trend in $\mathrm{d}H/\mathrm{d}t$ due to increasing air content in the ice below the firn layer is added to v_ice. v_by, only relevant over ice shelves, is equal to the negative change in firn mass divided by the density of sea water. The contributions of ice-dynamical imbalance and bedrock motion to changes in the surface elevation are not included in the current model.

2.4 RACMO2.3p2 forcing

IMAU-FDM is forced at the upper boundary with 3-hourly fields of snowfall, sublimation, snowdrift erosion, 10 m wind speed, surface temperature, snowmelt and rainfall. In this study, FDM v1.2A is driven with ERA5 reanalysis data (Hersbach et al., 2020), dynamically downscaled with the regional atmospheric climate model RACMO2.3p2 (Van Wessem et al., 2018). RACMO2.3p2 is run with a horizontal resolution of 27 km×27 km, and the simulation covers the period 1979–2020. RACMO2.3p2 has been thoroughly evaluated for the Antarctic and Greenland ice sheets by Van Wessem et al. (2018) and Noël et al. (2018), respectively. In comparison to RACMO2.3p1, RACMO2.3p2 employs upper-air relaxation, has updated topography, and has a retuned cloud scheme and modified snow properties. Upper-air relaxation is the indiscriminate nudging applied to the upper part of the atmosphere only. The most important effects of these changes are (i) increased snowfall in the ice sheet interior, (ii) reduced snowdrift sublimation and (iii) increased snowmelt (Van Wessem et al., 2018). Several recent studies use output from FDM v1.1 forced with RACMO2.3p2 over the AIS (e.g. Schröder et al., 2019a; Keenan et al., 2021), but the impact of the improved forcing has not yet been described. As snowfall and surface melt are defining drivers of the state of the firn, we also assess the impact of the updated forcing in this paper. In addition to the updates of RACMO2.3p2, previous FDM simulations for Antarctica used dynamically downscaled ERA-Interim reanalysis data, which have been replaced by ERA5 in 2019. This results in an improved forcing. For example, surface mass balance anomalies in Dronning Maud Land from 2006 to 2018 are better captured in ERA5 compared to ERA-Interim (Gossart et al., 2019). We expect that forcing differences have a larger impact on the temporal evolution of the surface height compared to adjustments in the fresh-snow density and dry-snow densification parameterizations (Verjans et al., 2021).

2.5 Observational data

To calibrate and evaluate the firn model, we collected published firn density profiles from widely varying locations across the AIS (Fig. 1). We used 112 density profiles from firn cores and 8 density profiles from neutron density probe measurements, combining multiple datasets (Van den Broeke, 2008; Schwanck et al., 2016; Bréant et al., 2017; Fernandoy et al., 2010; Montgomery et al., 2018; Fourteau et al., 2019; Olmi et al., 2021; Winstrup et al., 2019). The firn cores and neutron density probe measurements are mostly obtained in summer months between 1980 and 2020. For calibrating the fresh-snow density parameterization, observations of the density of surface snow, defined as the top 0.5 m, from 61 firn cores could be used. Eight additional surface snow density values from neutron density probe measurements and two from firn cores (Montgomery et al., 2018) were added after calibration and are thus only used for evaluation. We used densities of the top 0.5 m; as this depth is representative of many in situ surface snow density measurements, it generally consists of snow from multiple snowfall events, and it captures the densification by wind and redistribution of drifting snow (Brun et al., 1997). For calibrating the MO fits, only dry firn cores are used, as Eq. (3) describes dry-snow densification only. A firn core is considered dry if its location experiences on average less melt than 3 % of the average annual accumulation. For the MO fits, 93 dry firn cores could be used, which improves upon Ligtenberg et al. (2011), who used data from 48 dry firn cores. In addition to the 93 dry cores used for calibration of the MO ratios, another 11 wet cores and 4 additional dry cores that were later added were used to evaluate the z₅₅₀ and z₈₃₀ model output. For evaluating FAC output, 31 firn cores could be used. Table S1 in the Supplement lists all measurements that have been used, the corresponding coordinates, methods and citations.

In order to verify whether the modelled surface elevation trends and seasonal variability align with observations, we used a multi-mission satellite AIS altimetry product covering the period from January 1992 to December 2015 as presented by Schröder et al. (2019a). In this product, data from the ERS-1, ERS-2, Envisat, ICESat and CryoSat-2 missions are combined into a monthly time series on a 10 km×10 km grid, with the polar gap, floating ice shelves and steep terrain being excluded. Validation of this dataset with in situ and airborne observations showed that the products are successfully combined (Schröder et al., 2019a). The altimetry product has been resampled to the FDM grid using bilinear resampling. The measurement precision over flat terrain ranges from 5 to 15 cm, depending on the satellite. However, as the precision decreases with increasing surface slope, the uncertainty can reach several metres in the coastal regions with steep topography. Offsets due to varying penetration depth of the radar signals (thus relevant for ERS-1, ERS-2, Envisat and CryoSat-2), caused by variations in electromagnetic properties of the ice sheet, are not included in this uncertainty estimate. These additional offsets can reach up to several centimetres per year and can for instance increase the seasonal amplitude at continent-wide scales by up to several centimetres (Nilsson et al., 2022). While corrections are applied to account for variations in penetration depths, these corrections are in many cases unable to fully correct the artificial signals, which can e.g. be seen from the intermission variation in seasonal amplitude (Schröder et al., 2019a; Nilsson et al., 2022). On the other hand, signal penetration is negligible for laser measurements (such as ICESat, ICESat-2, Operation IceBridge). However, ICESat and Operation IceBridge datasets have relatively poor spatial and temporal coverage, and we believe that the record of ICESat-2 is currently too short (only 2 years and 2 months overlap with FDM v1.2A and no overlap with FDM v1.1p1) to provide a viable evaluation. Therefore, we used the multi-mission satellite altimetry product as presented by Schröder et al. (2019a), which has a good spatial and temporal coverage. Our firn simulations only provide elevation change due to anomalies in firn thickness caused by variations in e.g. snow accumulation, surface temperature, firn densification and melt. The difference between the observed and simulated elevation change therefore represents a combination of elevation changes due to ice-dynamical imbalance, bedrock motion, errors in the firn model and its climate forcing, and errors in the altimetry observations. For this comparison, we therefore exclude regions with known ice-dynamical imbalance, by using the ice-dynamical imbalance mask of Shepherd et al. (2019), which mainly consists of Pine Island, Thwaites, Totten and Getz Glacier basins and the Kamb Ice Stream.

2.6 Sensitivity analysis

To study the model sensitivity to input and parameter uncertainty we performed a sensitivity analysis for the 106 locations shown in Fig. S1a in the Supplement. These locations include 95 locations from firn core sites (the firn cores initially used for MO calibration and z₅₅₀ and z₈₃₀ evaluation are located at 95 unique grid points). To improve the representativeness of the climatic conditions over the AIS, we used 11 additional locations in this analysis that are located in high-accumulation or low-temperature regions, as these areas are underrepresented in the 95 firn core sites (Fig. S1b and c).

All the performed experiments are listed in Table 2. Verjans et al. (2021) and Lundin et al. (2017) showed that differences in accumulation, temperature, firn densification formulation and fresh-snow density can lead to a substantial spread in modelled firn thickness and air content. For each of these components, we separately performed a range of model sensitivity tests, which includes the spin-up as well as the main run. The fresh-snow density was uniformly varied with the RMSE from the evaluation with in situ measurements (±30 kg m⁻³; Sect. 3.1). To test the sensitivity to uncertainty of the dry-snow densification rate, we perform a simulation where we use the 95 % confidence interval boundaries of the MO fits (Eqs. 4 and 5, Fig. 3). To test the sensitivity to uncertainty in accumulation and temperature forcing, we use the standard deviation of snowfall and temperature among regional climate models and reanalysis products (Carter et al., 2022) (Fig. S2 in the Supplement). This spread varies spatially but is assumed to be constant through time. The accumulation is varied with the ratio of the ensemble standard deviation to the ensemble mean and the temperature with the ensemble standard deviation.

Carter et al. (2022)Carter et al. (2022)Stenni et al. (2017)Thomas et al. (2017)Stenni et al. (2017)Thomas et al. (2017)

Table 2 Overview of the sensitivity experiments. The inputs and parameters are adjusted by subtracting and adding the uncertainties either during the entire run or only during the spin-up. The values used are based on the references provided or on analysis within this work.

^∗ The values on the left side of the slash indicate the variation for the Antarctic Peninsula and Ellsworth Land and the values on the right side of the slash for the remainder of the AIS. Temperature and accumulation over the Antarctic Peninsula and Ellsworth Land are the second-to-last loop of spin-up reduced by 6.66 % and 0.66 K and in the last loop with 3.33 % and 0.33 K; see main text.

Another source of forcing uncertainty is the climate forcing during the spin-up period. As explained in Sect. 2.2, the spin-up forcing is obtained by looping over the 1979–2020 forcing data. However, firn core (Thomas et al., 2015, 2017; Medley and Thomas, 2019) and isotope (Stenni et al., 2017) studies show that in the Antarctic Peninsula and Ellsworth Land, the accumulation and temperature were typically about 10 % and 1 K lower during the last few centuries. In addition, over the remainder of the AIS, the accumulation and temperature were approximately 5 % and 0.5 K higher or lower during the last few centuries (Thomas et al., 2017; Stenni et al., 2017). To investigate the typical effect this has on the IMAU-FDM results, we include schematic sensitivity tests in which we adjusted the accumulation and temperature forcing with these values only during the spin-up period. To mimic a gradual increase in precipitation over the Antarctic Peninsula and Ellsworth Land, we performed a spin-up experiment in which the precipitation up until the third-to-last loop (the 42-year-long reference period) is decreased by 10 %, in the second to last by 6.66 % and in the last by 3.33 % with respect to the mean precipitation of 1979–2020. To mimic a gradual increase in temperature over the Antarctic Peninsula and Ellsworth Land, we create a spin-up experiment in which the temperature up until the third to last is decreased by 1 K, in the second to last by 0.66 K and in the last by 0.33 K. Over the other regions of the AIS, accumulation and temperature were both decreased and increased with 5 % and 0.5 K during the spin-up period. As the uncertainties between temperature and accumulation during the spin-up period are dependent, we also performed two tests in which we simultaneously adjusted the temperature and precipitation.

For all our experiments, except for the experiments where we perturbed the coefficients of the MO fits, we reran our model calibration procedure to get the optimal MO fits. We then ran the model with these updated MO fits and compared the outputs to our reference FDM v1.2A run. We investigate the impact of the uncertainties on the FAC and surface elevation change. As the model set-up of FDM v1.2A assumes a firn layer in a steady state, there is no net surface elevation trend in these simulations over the period 1979–2020, except for sensitivity tests in which the accumulation and/or temperature during the spin-up period are adjusted. In the remainder of the sensitivity experiments, we therefore compare surface elevation change over the period 2015–2020, as most locations experienced clear trends in surface elevation over that period. By correlating the sensitivities of FAC and surface elevation change of each experiment to climatic variables at the sample locations, we can expand the uncertainty estimates over the entire ice sheet. The empirical relations for each sensitivity experiment are shown in Fig. S3 in the Supplement. The uncertainties of each test were combined quadratically to obtain a value for the total uncertainty, where for the spin-up experiments only the combined temperature and accumulation experiments are included.

3 Calibration and model performance

3.1 Fresh-snow density

To calibrate and evaluate the fresh-snow density parametrization, we performed simulations with FDM v1.2A, using the parameterizations of Kaspers et al. (2004) (Eq. 1), Lenaerts et al. (2012) (Eq. 2), the latter with updated constants, and FDM v1.1p1. These simulations are abbreviated as FDM FS-K, FDM FS-L, FDM v1.2A and FDM v1.1p1. The fit coefficients and statistics of the simulations are listed in Table 3. Figure 2a shows the modelled against observed surface snow densities. FDM FS-K somewhat overestimates the surface snow densities (bias=17.6 kg m⁻³), especially for high surface density, and FDM v1.1p1 to a lesser extent (bias=10.3 kg m⁻³). It is the update in the atmospheric forcing that causes a poorer performance of FDM FS-K compared to FDM v1.1p1. FDM FS-L substantially underestimates surface snow densities ($\text{bias}=-\mathrm{45.1}$ kg m⁻³) but has a higher correlation with the observations than FDM FS-K. The degree of underestimation increases with increasing wind speed, which suggests that a higher value of fit coefficient C is needed to properly capture the impact of wind compaction. This aligns with the fact that IMAU-FDM does not include densification by wind packing. In the optimized fit coefficients (FDM v1.2A, Table 3) the value of C is 2.6 times higher than the value reported by Lenaerts et al. (2012). The wind dependency is also 2.4 times stronger than proposed by Kaspers et al. (2004), which is in agreement with Sugiyama et al. (2012) who found a factor 2.8. The RMSE of the surface snow density of FDM v1.2A is, respectively, 28 % and 23 % lower compared to FDM FS-K and FDM v1.1p1, due to both minimizing the bias and improving the representation of the variability. Compared to FDM FS-K, the surface snow density is reduced by 18 kg m⁻³ on average, especially in the high-accumulation margins, and to a lesser extent in windy escarpment regions (Fig. 2b). With the 10-fold cross evaluation for FDM v1.2A, we found that the RMSE for the surface snow density slightly increases to 30.4 kg m⁻³ (compared to 28.8 kg m⁻³ for FDM v1.2A). The fresh-snow density parameterization for Greenland used in FDM v1.2G is only a function of annual average temperature (Brils et al., 2022), owing to a lack of co-located surface snow density, temperature and wind speed observations on the GrIS. For all sensitivity tests presented here, we calculated the associated MO fits. However, the impact of the different MO fits on the surface snow densities is small (<0.1 %).

Table 3 Fit coefficients and evaluation of the fresh-snow density parameterizations (Eqs. 1 and 2). The fit coefficients of FDM FS-K and FDM FS-L are taken from Kaspers et al. (2004) and Lenaerts et al. (2012), respectively.

^a Using annual average accumulation, temperature and wind speed. ^b Using instantaneous temperature and wind speed.

Figure 2 (a) Scatter plot of simulated against observed surface snow density using FDM v1.2A, FDM FS-K (Kaspers et al., 2004), FDM FS-L (Lenaerts et al., 2012) and FDM v1.1p1. (b) Map of average difference between the simulated surface snow densities of FDM v1.2A and FDM FS-K for the locations with observations.

3.2 Dry-snow densification rate

To calibrate the dry-snow densification rate we first performed a simulation of FDM v1.2A without MO corrections, thus only using Eq. (3). The resulting MO fits and statistics are listed in Table 4 and shown in Fig. 3a and b. The optimal MO₅₅₀ fit is less steep with respect to accumulation (β=0.117 vs. 0.151) than in FDM v1.1p1. The fresh-snow density is now independent of accumulation and generally lower, and hence more densification is needed at high-accumulation sites to match the density measurements. This might explain the slightly deteriorated correlation for MO₅₅₀. For the MO_830* regression, we found that a power fit (FDM v1.2A) compared to a logarithmic fit (FDM v1.2A-log) improves correlation (R²=0.88 compared to 0.83). Furthermore, a minimum value is now explicitly derived with the power fit so that a rather arbitrary minimum value no longer needs to be prescribed. When the locations are rerun using the new MO values with the power fit for MO_830*, the resulting RMSE of the modelled z_830*, z₈₃₀ (firn thickness) and FAC are, respectively, 41 %, 39 % and 48 % lower compared to the logarithmic fit (Fig. 3c and d). The update in atmospheric forcing without updated MO calibration does not result in an improvement in FAC (shown by the similar RMSE of FDM v1.1p1 and FDM v1.1p2, which is FDM v1.1 forced with the updated forcing RACMO2.3p2). Also, the RMSE of FAC simulated with FDM v1.2A is, respectively, 18 % and 20 % lower compared to FDM v1.1p1 and FDM v1.1p2, which indicates a smaller improvement compared to FDM v1.2A-log (48 %). The poor performance of FDM 1.2A-log shows that the additional observations necessitated the introduction of the power fit. When we only include the seven observations that were used to calibrate FDM v1.1, and two wet firn cores, the RMSE of FAC simulated by FDM v1.2a is 2.79 m, which is 19.5 % and 9.0 % lower compared to FDM v1.1p1 and FDM v.1.1p2 for these observations.

Table 4 Fit coefficients and statistics of the MO fits (Eqs. 4 and 5). The fit coefficients and statistics of FDM v1.1p1 are taken from Ligtenberg et al. (2011). Please note that R² for v1.1p1 applies to the calibration dataset used in that study. The FDM v1.2A-log uses the settings of FDM v1.2A; only Eq. (3) is optimized using a logarithmic MO_830* fit. The root mean squared deviation (RMSD) quantifies the deviation with respect to the fitted function.

Figure 3 Scatter plots defining the MO ratios for (a) z₅₅₀ and (b) z₈₃₀* as a function of the annual average accumulation including various MO fits. The blue-shaded areas indicate the 95 % confidence intervals. (c) Scatter plots of the simulated against observed z_830* (z₈₃₀ minus z₅₅₀) for FDM v1.2A and FDM v1.2A-log. (d) Scatter plots of the simulated against observed FAC for FDM v1.2A, FDM v1.2A-log, FDM v1.1p1 and FDM v1.1p2. (e) Map of the difference between simulated firn thickness (z₈₃₀) of the two tests (FDM v1.2A vs. FDM v1.2A-log) for locations with observations. The statistics and fits of the MO fits are listed in Table 4.

Using a power compared to logarithmic fit results in higher MO_830* values for low (<400 mm yr⁻¹) and high (>1400 mm yr⁻¹) accumulation regimes and lower MO_830* values for intermediate accumulation regimes (400–1400 mm yr⁻¹) (Fig. 3b). Therefore, the densification rate is increased in high- and low-accumulation regions and reduced in intermediate-accumulation regions. This results in decreased firn thickness in high- and low-accumulation regions and in increased firn thickness in intermediate-accumulation regions (Fig. 3d). Especially large differences are found in high-accumulation areas where the MO fit curve decreases less rapidly in FDM v1.2A-power compared to FDM v1.2A-log. With the 10-fold cross evaluation for FDM v1.2A, we found that the RMSE for z₅₅₀ slightly increases to 2.33 m (compared to 2.26 m for FDM v1.2A) and the RMSE of z_830* to 8.91 m (compared to 8.33 m for FDM v1.2A). For the 11 wet firn cores, we found that the RMSE of modelled z₅₅₀, z_830* and FAC amounts to 3.17, 14.34 and 4.40 m, respectively, and thus larger than the RMSEs of the entire firn core dataset (2.26, 8.33 and 2.74 m, respectively).

In FDM v1.2G, the version for the Greenland ice sheet, a logarithmic MO₈₃₀ fit was applied. Dry firn cores from the GrIS cover a smaller range of average annual accumulation, 80–680 mm yr⁻¹, compared to 20–960 mm yr⁻¹ in Antarctica. If we only include Antarctic cores from the 80–680 mm yr⁻¹ accumulation range, we find similar model fits to the data for the logarithmic model (R²=0.67) and the power model (R²=0.68). Another difference between FDM v1.2G and FDM v1.2A is that in FDM v1.2G an almost constant value of 0.67 for the MO₅₅₀ fit was found, which implies a linear correlation between the densification rate and the accumulation rate. In FDM v1.2A we do find a moderate (R²=0.37, β=0.12) correlation of MO₅₅₀ with accumulation; i.e. the AIS densification rate depends more strongly on accumulation than temperature compared to the GrIS. Again, if we only include cores from the 80–680 mm yr⁻¹ accumulation range this correlation weakens (R²=0.08, β=0.07).

4 Spatial and temporal variability of AIS firn depth and density

4.1 Spatial density variability

Using the improved firn model FDM v1.2A, we simulate the spatial and temporal patterns of the AIS firn characteristics at 27 km×27 km horizontal resolution. The spatial patterns of the densities, including observed values, are shown in Fig. 4. Figure 4a shows the surface snow density (top 0.5 m), which is on average 376 kg m⁻³. The patterns are primarily a reflection of the variation in fresh-snow density (Eq. 2), with the lowest values on the cold and calm East Antarctic plateau and higher values in the windy escarpment regions and along the milder coastal margins. In addition to this, the surface snow density further increases over the ice shelves, where the average surface snow density amounts to 414 kg m⁻³, and over the coastal margins, where intermittent melting and refreezing of meltwater increases the near-surface density.

Figure 4 Maps from FDM v1.2A of the average (a) surface snow density (defined as the top 0.5 m), depths of critical density levels of (b) ρ=550 kg m⁻³ (z₅₅₀) and of (c) ρ=830 kg m⁻³ (z₈₃₀) (the firn thickness), and (d) firn age at the critical density level of (c) ρ=830 kg m⁻³. Circles indicate the values found by field measurements.

In line with observations (Winther et al., 2001), FDM v1.2A predicts that nearly all (>99 %) of the AIS is covered by a layer of firn. This represents the ice sheet fraction where average SMB is positive. The spatial pattern of the depths of the critical density levels z₅₅₀ and z₈₃₀ are shown in Fig. 4b and c, which are on average 14.1 and 73.8 m. The patterns are roughly the inversed pattern of the surface snow density, which constitutes the upper boundary of the depth–density profile. For simplicity, henceforth we define the firn thickness as z₈₃₀. The patterns vary spatially across climatic regions with temperature, accumulation and surface melt as drivers. The highest z₅₅₀ and z₈₃₀ values (30.6 and 114.2 m, respectively) are found in the interior, where the densification rate is low due to low temperatures, and the surface snow density is low. Generally, z₅₅₀ and z₈₃₀ are small along the warm coastal margins due to faster densification, higher surface snow densities and melt. However, the z₈₃₀ values can be relatively high in high-accumulation regions along the coast, particularly in West Antarctica and the Antarctic Peninsula. The high accumulation rates in these regions result in younger, thus less densified, snow with depth and thus thick firn layers despite relatively high temperatures. The z₅₅₀ and z₈₃₀ values are the lowest in relatively warm and dry regions on and around Ross, Filchner–Ronne and Amery ice shelves, where fresh snow is buried more slowly, and the densification rates are high due to relatively high temperatures.

Figure 4d shows the age of the firn at the pore close-off depth, here assumed to equal z₈₃₀, which is on average 754 years. The combination of strongly increasing firn depth and decreasing accumulation towards the interior leads to firn ages at the pore close-off depth of up to 3240 years in central East Antarctica, whereas the pore close-off firn age in warm and wet coastal margins can be as low as 20 years. When comparing this to firn age inferred from δ¹⁵N measurements for 15 locations (Bréant et al., 2017), we find an RMSE=231 years and R²=0.99. On average, the absolute errors are 25 % of the observed values. In comparison, the absolute errors of z₈₃₀ are 15 % of the observed values on average, and R²=0.73.

4.2 Firn air content

Firn air content (FAC), the total (vertically integrated) pore space, is an important output of firn models, as it is used to convert observed volume changes to mass changes and column thickness to ice mass and as it is an indicator of the meltwater buffering capacity of the firn. FAC is expressed in metres and represents the change in depth that would occur if the firn layer would be compressed to the density of glacier ice. The average FAC of the AIS according to IMAU-FDM is 22.8 m, and it varies spatially between 0 to 40 m across climatic regions (Fig. 5a). The pattern resembles the z₈₃₀ pattern (Fig. 4c); thick firn layers contain most air. High values are found in the cold interior and along parts of the coast, particularly in West Antarctica, and in areas with high accumulation. Low values are found along the warmer coastal margins, with ice shelves having an average FAC of 15.9 m. These average FAC estimates are 5 % and 7 % lower compared to values simulated with the Community Firn Model forced with MERRA-2 climate data of 24.0 and 17.0 m, respectively (Medley et al., 2022).

Figure 5 Maps from FDM v1.2A of average (a) firn air content (FAC), (b) the difference in FAC (FDM v1.2A vs. FDM v1.1p1), (c) peak-to-peak FAC difference and (d) relative peak-to-peak FAC difference compared to the average FAC from (a). The peak to peak is defined as the difference between the highest and the lowest FAC value occurring during the study period (1979–2020).

The average FAC over the AIS simulated with FDM v1.2A is 0.6 m higher than that simulated with FDM v1.1. The spatial FAC differences between FDM v1.2A and FDM v1.1 are shown in Fig. 5b. In the highly elevated and low-accumulation regions the updated FAC is lower, whereas in the intermediate-accumulation regions the updated FAC is higher. This implies that the range of FAC is smaller and less spatially variable than previously modelled. The differences in the FAC pattern can be explained by differences in (i) the MO_830* fit, (ii) the fresh-snow density, (iii) snowfall and (iv) snowmelt. FDM v1.2A has higher MO_830* values for low- and high-accumulation regions compared to FDM v1.1p1, which results in a faster densification and therefore lower FAC values. In intermediate-accumulation regions the MO_830* values are lower, resulting in a slower densification and higher FAC values. The fresh-snow densities of FDM v1.2A are on average lower (19.1 kg m⁻³, Fig. 2), which contributes to FAC values being higher on average. On top of that, increased accumulation in the interior and increased snowmelt along the coastal margins in RACMO2.3p2 compared to RACMO2.3p1 (Van Wessem et al., 2018) also contribute to the general pattern of higher FAC in the interior and lower values along the coastal margins. However, in the highly elevated and lowest-accumulation regions, we find a reduction in FAC, which implies that the MO_830* update outweighs the effects of higher accumulation and lower fresh-snow density in these regions.

Figure 5c and d focus on temporal variability in FAC. Figure 5c shows the peak-to-peak FAC difference, which represents the difference between the highest and the lowest FAC values occurring during the study period (1979–2020). The peak-to-peak FAC difference is low in the interior, less than 0.1 m, despite the high FAC. Along the coastal margins, the peak-to-peak FAC difference is higher, up to 5 m, even though the FAC is generally lower. Therefore, the relative peak-to-peak FAC difference is low in the interior (<1 %) and high near the ice sheet margins (>5 %) (Fig. 5d). The average relative peak-to-peak FAC difference for the ice shelves, excluding the two largest ones, the Ross and Filchner–Ronne ice shelves, amounts to 13 %. This high peak-to-peak variability is caused by a combination of substantial temporal variability in accumulation and melt. The temporal FAC variability and the driving processes are discussed further in Sects. 4.3 and 4.4 below.

4.3 Seasonal firn air content and firn depth variability

In order to better understand the FAC temporal variability and its drivers, we decompose the signal into seasonal and interannual contributions. Figure 6a shows the average seasonal cycle of the vertical firn surface velocity (v_tot) including its components from Eq. (6). Four main processes drive the seasonality of the firn thickness: snowfall (v_snow), firn densification (v_fc), melt (v_melt) and sublimation (v_sub). Averaged over the ice sheet, snowfall is the highest in winter, while firn densification, melt and sublimation peak in summer. The net result is that the firn thickness and FAC steadily increase over winter due to enhanced accumulation and reduced firn densification, melt, and sublimation (Fig. 6b). This is followed by a sharp decrease over summer due to reduced accumulation and increased firn densification, melt, and sublimation. The firn height change by ice flux divergence (v_ice) is constant by definition, while the AIS average impacts of snowdrift (v_snd) and buoyancy (v_by) are small (<1 %). The vertical bars in Fig. 6a indicate the interannual variability, which is driven solely by snowfall except for the two summer months, which also have contributions from sublimation and melt.

Figure 6 Ice-sheet-averaged seasonal cycle from FDM v1.2A of the (a) vertical firn surface velocity, including its separate components from Eq. (6), and of (b) surface elevation and firn air content (FAC). The bars and coloured areas indicate the interannual standard deviations. Maps from FDM v1.2A (c) of the average seasonal firn thickness amplitude and (d) of the average day of the year of the maximum firn height.

The seasonal amplitudes of the ice-sheet-wide firn thickness and FAC are 3.1 and 2.1 cm, which amounts to Antarctic-wide integrated volumes of 428 and 307 km³, respectively. This agrees with the modelled firn thickness amplitude of 3 cm found by Medley et al. (2022). The averages of all local seasonal amplitudes are considerably higher, 6.8 and 4.2 cm, amounting to volume changes of 892 and 562 km³, respectively (Fig. 6c). The reason is that there is a spatial variation in the timing of the seasonal maximum and minimum (Fig. 6d). The difference between the firn thickness and FAC seasonal cycles (6.8 vs. 4.2 cm) suggests that around 62 % of the seasonal surface elevation fluctuations are caused by a change in air content rather than actual mass change. This result implies that not including the seasonal cycles of FAC in altimetry studies can generate significant errors, especially when investigating changes in subannual periods and in regions with large seasonal SMB cycles.

Figure 6c shows that the seasonal firn thickness amplitude can be up to 3 m in the western Antarctic Peninsula, caused by a rare combination of high snowfall and strong melt. Maps of mean annual precipitation, melt and sublimation are shown in Fig. S9 in the Supplement. In regions with little melt, such as coastal West Antarctica, the values are typically 1 m. In interior East Antarctica, values of less than 2.5 cm are found, as the annual average accumulation is low, and melt is absent. The average day of the year of the maximum firn height generally occurs during late winter or spring, before the summer decrease in accumulation and increase in temperature (Fig. 6d). In some regions accumulation and firn depth peak in late summer or early winter, e.g. along the Weddell Sea coast, in Enderby Land, in Adélie Land and in parts of the interior.

4.4 Interannual firn air content and firn depth variability

Figure 7 shows time series of integrated FAC, the cumulative surface temperature anomalies (against the long-term mean) and the vertical firn surface velocity anomalies (against the long-term mean) for the entire AIS (Fig. 7a) and above/below 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Fig. 7b and c). In addition, the surface velocity anomalies are broken down into their components (Eq. 6). The figure confirms that the seasonal firn thickness and FAC variability are driven by snowfall, densification, melt and sublimation, resulting in a firn thickness peak in late winter. The figure also reveals decadal firn thickness and FAC variability, mainly resulting from fluctuations in snowfall (46 % of the total cumulative dH anomaly), and that firn densification, despite the slower response timescale, reduces these snowfall-induced fluctuations by about 20 % (calculated from the difference in standard deviation (SD) between the cumulative anomalies of v_snow and v_snow+fc). We also find interannual variability in melt, but snowfall variability dominates. Below 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ the variability in firn thickness and FAC is substantially larger, which can be attributed to the higher accumulation (70 % of the total accumulation falls on only 40 % of the ice sheet), in combination with larger temporal variability in snowmelt and firn densification. Above 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ the relative contribution of firn densification to the total cumulative dH anomaly is larger (32 %) than below 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, where this is only 14 %. This difference can partly be explained by the larger interannual temperature variability above 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (SD in annual means = 0.78 K compared to 0.48 K), the high sensitivity of firn compaction to temperature variability for firn at low temperatures and the absence of melt. The interannual variabilities in snowfall and temperature are related to low- and mid-latitude climate oscillations such as the El Niño–Southern Oscillation and the Southern Annular Mode (Kwok and Comiso, 2002; Marshall et al., 2017). Figure 7 shows that the magnitude of the decadal FAC variability (SD=1.8 cm) constitutes a large fraction (67 %) of the total surface elevation change variability (SD=2.7 cm). This reiterates the necessity of the common practice to remove decadal variability in FAC from altimetry signals to avoid errors in mass change estimates.

Figure 7 Time series from FDM v1.2A of FAC, the cumulative anomalies (against long-term mean) of surface temperature, the vertical firn surface velocity and its separate components from Eq. (6). Time series are shown for (a) the entire ice sheet, (b) the part of the ice sheet situated above 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ and (c) the part of the ice sheet situated below 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$

5 Comparison with altimetry and FDM v1.1 surface elevation change

In Sects. 3.1, 3.2 and 4 we presented and evaluated the spatial variation in firn thickness and density. To also evaluate the temporal variation, we compare our simulated elevation trends and seasonal variability to multi-mission satellite altimetry observations reported by Schröder et al. (2019a) for the period 1992–2015.

5.1 Seasonal amplitude

The altimetry observations prior to 2003 exhibit a larger short-term variability as a result of the lower measurement precision of the satellites during that period (Schröder et al., 2019a; Nilsson et al., 2022) and are excluded from the seasonal amplitude comparison here. To enable comparison with the FDM v1.1p1 simulation, we focus on the period 2003–2015. To reduce the impact of the long-term trends we remove a linear trend from all datasets. After doing so, we find that the seasonal amplitude averaged over the ice sheet of the altimetry observations amounts to 5.2 cm, of FDM v1.2A to 3.1 cm (60 % of the observed amplitude) and of FDM v1.1p1 to 2.5 cm (48 % of the observed amplitude). The higher seasonal amplitude of FDM v1.2A results from increased snowmelt along the coastal margins in RACMO2.3p2 compared to RACMO2.3p1 in combination with a lower fresh-snow density (Fig. 2). Using updated corrections and a normalization of seasonal amplitudes from different missions, Nilsson et al. (2022) found that seasonal amplitudes during 1992–2016 are 47 % smaller than the ones from Schröder et al. (2019a). This indicates that our AIS-wide seasonal amplitude lies within the range of observational uncertainty.

5.2 Decadal variability

In Fig. 8 we compare observed (Fig. 8a) and simulated 2003–2015 surface elevation trends. To obtain these numbers, a linear regression is fitted to monthly values for the period 2003–2015. The FDM v1.2A surface elevation trend over this period is positive in Dronning Maud Land and Enderby Land and negative in most other regions (Fig. 8b). FDM v1.1p1 (Fig. 8c) has roughly similar patterns as FDM v1.2A, but the regional magnitudes differ. Figure 8d shows the residual surface elevation trend of FDM v1.2A, which is calculated by subtracting the FDM v1.2A elevation trend from the altimetry trend, and thus combines unresolved trends and the errors in the firn model and altimetry. In Fig. 8e we show the difference in the altimetry agreement between FDM v1.2A and FDM v1.1p1, which is calculated by subtracting the absolute residual surface elevation trend compared to the altimetry (Fig. 8d) of FDM v1.1p1 from the similarly derived absolute residual of FDM v1.2A. The blue areas indicate an improvement in the altimetry agreement of FDM v1.2A compared to FDM v1.1p1. The average absolute residual trend of FDM v1.2A compared to altimetry has been reduced by 17 % compared to FDM v1.1p1 (from 2.6 to 2.1 cm yr⁻¹). The improvement is most notable in Dronning Maud Land, Wilkes Land and Adélie Land. In Dronning Maud Land the FDM v1.2A trend has become more positive, and in Wilkes Land and Adélie Land the trend has become less negative (Fig. 8b and c; for region names see Fig. 1). In the Antarctic Peninsula and Ellsworth Land the altimetry agreement of FDM v1.2A has reduced compared to FDM v1.1p1. In these regions, a substantial residual positive trend remains (>10 cm yr⁻¹) (Fig. 8d), which could be due to a long-term increase in accumulation and/or ice-dynamical thickening. Smith et al. (2020) used ICESat data to show that the mass gain is mainly located along the steep slopes of the Antarctic Peninsula and decreases with distance from the ocean, which is indicative of an increase in snow accumulation. The increase in accumulation on centennial timescales in these regions (e.g. Thomas et al., 2015, 2017; Medley and Thomas, 2019) implies that the actual firn column is not in balance with the 1979–2020 climate, as assumed in the FDM spin-up procedure. As a result, the vertical downward ice flow (Eq. 6) is overestimated, which results in underestimated surface elevation change. This is further discussed in the next paragraph.

Figure 8 Maps of trends in surface elevation change over the period 2003–2015 from (a) multi-mission altimetry (Schröder et al., 2019a), (b) FDM v1.2A and (c) FDM v1.1p1 obtained using linear regression. (d) Residual trend of altimetry minus FDM v1.2A. (e) Difference in the altimetry agreement between FDM v1.2A and FDM v1.1p1 (absolute residual of FDM v1.1p1 compared to altimetry minus absolute residual of FDM v1.2A compared to altimetry). Regions within the ice-dynamical imbalance mask reported by Shepherd et al. (2019) are excluded in (d) and (e) and indicated by the grey-shaded area.

In Fig. 9 we compare time series of altimetry observed elevation change over the period 1992–2015 to simulated elevation change for 11 locations across the AIS. These locations were selected as they cover the main distinct patterns from Fig. 8 and have continuous observations. The time series have been smoothed using a 6-monthly moving average window in order to reduce the impact of seasonal variations in radar penetration depth. The red-shaded area is the aggregated uncertainty range from the sensitivity analysis, which is primarily determined by uncertainty in the climate forcing during the spin-up period (see the next section for details of the sensitivity analysis). The agreement of FDM v1.2A compared to FDM v1.1p1 has improved (R²=0.58 compared to 0.36), which is mainly due to changes in the RACMO2.3p2 forcing, given the strong similarity between the elevation change of FDM v1.2A and FDM v1.1p2. The altimetry observations prior to 2003 exhibit stronger short-term variability (+9 % SD compared to after 2003), in line with the smaller measurement precision (Schröder et al., 2019a). The agreement with FDM v1.2A increases after 2003 (R²=0.93 compared to 0.64, RMSE reduced by 49 %); however the altimetry variability remains higher than the simulated variability (+13 % SD compared to FDM v1.2A), as these more recent observations still contain noise (<11 cm over flat terrain) and have variations in radar penetration depth (Schröder et al., 2019a). The uncertainty bands in Fig. 9 show that a large part (73 %) of the difference between FDM v1.2A and altimetry falls within the uncertainty range of the firn model. As in Fig. 8d, we find substantial positive residual altimetry trends (∼14 and 5 cm yr⁻¹) at locations 5 (Antarctic Peninsula) and 6 (Ellsworth Land). The upper bound of the uncertainty for those two locations is determined by the sensitivity tests in which the climate during the spin-up period was colder and drier (spin-up T_s, $\dot{b}$ experiment), indicated by the dotted red lines in Fig. 9. With this assumption, the trend difference for location 5 greatly reduces and largely disappears for location 6.

Figure 9 Time series of the altimetry observations and modelled elevation change for 11 locations over the period 1992–2015. The time series have been smoothed using a 6-monthly moving average window. The modelled elevation change includes FDM v1.2A, FDM v1.1p1 and FDM v1.1p2. The number of each panel corresponds to the numbered locations on the maps. The red-shaded areas indicate the total uncertainty from the sensitivity analysis. The dotted red lines in panels (5) and (6) indicate the surface elevation change from the spin-up T_s, $\dot{b}$ experiment.

6 Sensitivity analysis

We investigated the impacts of uncertainties in the model densification rate expression, fresh-snow density, climate forcing and climate forcing during the spin-up period on the final FAC and surface elevation change results for (i) the 106 sample locations and (ii) expanded to the entire ice sheet (Table 5). For the sample locations, the sensitivities with updated MO fits are lower than with the FDM v1.2A fits (i.e. without updated MO fits), especially the FAC uncertainties, which are on average 41 % lower. The results for the sample locations with updated MO fits are then expanded to the entire ice sheet. Some differences between the ice-sheet-wide and sample locations exist, which can be explained by the sample locations being biased towards high accumulation, which we used to obtain optimal empirical relations for expanding the sample location uncertainties (Fig. S1). The average simulated FAC, for the sample locations and entire ice sheet (both with updated MO fits), is the most sensitive to the uncertainties in the MO fits (5.2 % and 5.1 %), temperature forcing (2.3 % and 1.8 %) and accumulation forcing (3.7 % and 4.0 %) and the least sensitive to uncertainties in fresh-snow density (0.7 % and 1.2 %) and the spin-up forcing (ranging from 0.5 % to 0.7 %). The simulated surface elevation change, for the sample locations and the entire ice sheet, is the most sensitive to uncertainties in the accumulation during the spin-up (22.4 and 19.0 mm yr⁻¹), followed by the accumulation forcing (7.8 and 5.6 mm yr⁻¹) and temperature during the spin-up forcing (5.2 and 4.4 mm yr⁻¹). As the uncertainties in temperature and accumulation during the spin-up period have opposite effects, the tests in which accumulation and temperature were adjusted simultaneously yield smaller surface elevation trend differences of 17.4 and 14.9 mm yr⁻¹. The simulated surface elevation change is the least sensitive to uncertainties in the MO fits (0.7 and 0.5 mm yr⁻¹), fresh-snow density (2.1 and 1.5 mm yr⁻¹) and temperature forcing (0.7 mm yr⁻¹).

Table 5 Overview of the relative change in annual average firn air content (FAC) and surface elevation change ($\mathrm{d}H/\mathrm{d}t$) of the ice-sheet-wide estimates and sample location experiments. The sample location experiments also include the experiments without adjusted MO fits. The FAC uncertainty for all experiments and surface elevation change uncertainty for the spin-up experiments are calculated over 1979–2020, and the surface elevation change uncertainty for the MO fits, ρ_s, T_s and $\dot{b}$ experiments are calculated over 2015–2020. An overview of the experiments including the prescribed uncertainties are listed in Table 2.

^a The values between brackets indicate the uncertainties of the experiments without updated MO fits. ^b For the spin-up experiments only the combined temperature and accumulation experiments are included.

The spatial patterns of the total FAC and surface elevation change uncertainties are shown in Fig. 10 and for all experiments individually in Figs. S4 and S5 in the Supplement. The total FAC uncertainty, dominated by the MO fits, accumulation and temperature uncertainties, is high (>10 %) along the high-accumulation coastal margins and in a part of the interior, especially near the Transantarctic Mountains. The surface elevation change uncertainty is dominated by the spin-up and accumulation uncertainty, and mainly follows the accumulation pattern, with high values (>25 mm yr⁻¹) along the high-accumulation coastal margins, especially in the Antarctic Peninsula, and in the Transantarctic Mountains.

Figure 10 Maps from FDM v1.2A of the estimated ice-sheet-wide total uncertainty in (a) firn air content (FAC) and (b) surface elevation change ($\mathrm{d}H/\mathrm{d}t$) from the sensitivity analysis. The FAC uncertainty for all experiments and surface elevation change uncertainty for the spin-up experiments are calculated over 1979–2020, and the surface elevation change uncertainty for the MO fits, ρ_s, T_s and $\dot{b}$ experiments are calculated over 2015–2020.

We show time series of the uncertainty in surface elevation change and FAC at the sample locations in Fig. S7a in the Supplement for all experiments and in Fig. S7b for the non-spin-up experiments, i.e. the MO fits, fresh-snow density, temperature and accumulation experiments. First, we discuss uncertainties in surface elevation change for non-spin-up and spin-up experiments. Next, we look at uncertainties in FAC. The uncertainty of the non-spin-up experiments (Fig. S7b) is only 4 % of the total surface elevation change uncertainty (Fig. S7a). This uncertainty is smaller than the 40 % that can be derived from Table 5. This is because an imposed uncertainty has opposite effects across the AIS and in time. For example, lowering the fresh-snow density increases the rate of a surface elevation increase/decrease. Table 5 lists the absolute uncertainties in surface elevation change over a short 5-year period. However, spatial and temporal variability in the sign of elevation change tend to cancel out the uncertainties. In contrast, the uncertainty in the spin-up experiments gives a spatially and temporally uniform response. This enables us to expand ice-sheet-wide estimates of uncertainty of the spin-up experiments into time series (Fig. S7c) by adding a uniform uncertainty for each time step.

Uncertainty in ice-sheet-wide FAC change is only 19 % of the uncertainty in ice-sheet-wide surface elevation change, similar to the 23 % for the suite of sample location experiments (Fig. S7a and c). The spatial pattern in FAC change uncertainty for the spin-up experiments is in general similar to the surface elevation change uncertainty, apart from a decreasing uncertainty in high-accumulation regions (>600 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$) (Fig. S6 in the Supplement). We attribute this to the fact that in high-accumulation regions, a large part of the upper firn layer does not originate from the spin-up period. Here, the main impact is on deeper high-density firn, meaning that the impact on a volume-to-mass conversion is limited compared to the large uncertainty in surface elevation change. However, this large uncertainty in the surface elevation change compromises our ability to separate ice-dynamical from SMB and firn-related surface elevation changes. In contrast, in the non-spin-up experiments, the magnitude of the FAC uncertainty is similar (∼85 %) to the uncertainty of surface elevation change.

For the Antarctic Peninsula and Ellsworth Land, accumulation and temperature have increased during the 20th century (Thomas et al., 2017). To evaluate our results from the spin-up sensitivity tests in these two regions, we used surface elevation trends from altimetry observations (Schröder et al., 2019a) (Fig. S8 in the Supplement). We found that the absolute residual trend of altimetry minus the spin-up sensitivity test was reduced by 25 % compared to the reference FDM v1.2A (from 6.1 to 4.5 cm yr⁻¹). When we leave out pixels with a z₈₃₀ age<42 years (6 % of the pixels), where we assume the underlying ice to have responded to the changing climatic conditions, the residual trend was reduced by 38 % (from 5.6 to 3.5 cm yr⁻¹). This suggests that the uncertainty for those locations (z₈₃₀ age<42 years) is overestimated, as v_ice can not be assumed constant. Over the remainder of the AIS, some locations experienced an increase in accumulation and temperature, whereas others experienced a decrease in accumulation and temperature over the last few centuries (Thomas et al., 2017). This has opposite effects on the surface elevation and FAC change, which suggests that the surface elevation and FAC change uncertainties averaged over the AIS are likely lower than the results of the spin-up temperature and accumulation uncertainty experiments (Fig. S7c).

7 Implications, remaining limitations and outlook

In this study, the fresh-snow density and firn compaction parameterizations of IMAU-FDM have been improved, and an updated atmospheric forcing has been used over the AIS. The resulting differences in firn changes over time impact e.g. the conversion of volume changes to mass changes in altimetry studies. Over the period 2003–2019, Smith et al. (2020) used volume changes from altimetry to estimate that the total mass change of East Antarctica amounts to 195 Gt yr⁻¹ and of West Antarctica to −245 Gt yr⁻¹. When FDM v1.2A is used instead of FDM v1.1p1, the mass change over the period 2003–2015 for East Antarctica is 7.3 Gt yr⁻¹ lower (187.7 Gt yr⁻¹) and for West Antarctica this is 1.2 Gt yr⁻¹ higher (−243.8 Gt yr⁻¹). These temporal firn differences are predominantly caused by the updated climatic forcing. In addition, updates in the fresh-snow density and firn densification parameterizations have a substantial impact on the time-averaged firn density profiles. These outputs are required to assess the meltwater buffering capacity of the firn layer and to convert column thickness to ice mass in order to calculate solid-ice discharge over the grounding line (Rignot et al., 2019).

The comparison with satellite altimetry measurements is used to assess the performance of FDM. However, especially in the interior, both the long-term trend and seasonal variability are small compared to the error range of the altimetry observations, which explains part of the remaining discrepancy that exists there, but also hampers the comparison (Verjans et al., 2021). Accurate in situ measurements of surface height evolution and firn densification are sparse. Future evaluation of firn models will therefore likely continue to depend on satellite altimetry comparison. The ICESat-2 ATL15 (Smith et al., 2021) is a precise Antarctic-wide 3-month height change laser product that has been available from May 2018 onwards, thereby providing an opportunity to assess seasonal and interannual variations, especially in the coming years when more data will become available. Smith et al. (2023) already used ICESat-2 measurements to evaluate regional climate and firn densification models in Greenland.

The sensitivity analysis shows that updating the MO fits reduces the average FAC sensitivity to uncertainties in fresh-snow density and climatic forcing by 41 % (Table 5), which implies that the MO factors partly correct for biases in these boundary conditions. Thus, the application of FDM v1.2A calibrated MO fits when using a different fresh-snow density or climate forcing is not recommended. The fresh-snow density parameterization is also dependent on the climate conditions derived from RACMO2.3p2, so any biases will impact the updates of the coefficients. In addition, the simulated fresh-snow density accounts for the wind-driven compaction, and the updated coefficients are thus not recommended for use in combination with a firn model that includes wind-driven densification of the upper firn layers.

Another limitation, which is evident from both the remote sensing comparison and the sensitivity analysis, is the assumption of a steady-state firn layer over the period with reliable climate forcing (1979–present). In our sensitivity analysis we show that correcting for this by estimating previous climatic conditions from ice core data reduces the difference between simulated and altimetry observed surface elevation trends by 38 % over the Antarctic Peninsula and Ellsworth Land. Future work could focus on further accounting for pre-1979 climatic conditions.

The densification scheme used in FDM v1.2A (Eq. 3) is developed for dry-snow densification. As it is used for both wet and dry locations, we assume that the densification rate of dry firn is equal to that of wet firn. Since the presence of liquid water impacts the evolution of grain size and shape, this may also impact the densification rate of firn; however due to a lack of physical understanding and available measurements, this has not been included in FDM v1.2A. The RMSE of the modelled z₅₅₀, z₈₃₀ and FAC for the wet firn cores are, respectively, 3.17, 14.34 and 4.40 m and thus larger than the RMSEs of the entire firn core dataset (2.33, 8.33 and 2.74 m, respectively). Uncertainties in the melt forcing (Carter et al., 2022) impact the simulated density profiles, which may contribute to the lower agreement and hampers this comparison.

8 Summary and conclusions

In this study, we present and evaluate a simulation of the contemporary characteristics of the Antarctic firn layer using the updated IMAU-FDM (IMAU Firn Densification Model), v1.2A, for the period 1979–2020. In this new version, the fresh-snow density and firn compaction parameterizations have been improved, and an updated atmospheric forcing has been used. The RMSEs when compared to observations of the fresh-snow density and firn thickness have been reduced by 28 % and 39 %, respectively, compared to the previous parameterizations. The RMSE of FAC is 18 % lower compared to the previous model version FDM v1.1p1. In general, the modelled fresh-snow densities have decreased, especially in high-accumulation areas, and to a lesser extent in windy escarpment regions. Firn thickness and FAC have decreased in low- and high-accumulation areas, whereas they have increased in intermediate-accumulation areas.

In accordance with observations, the firn thickness and density patterns vary spatially across climatic regions with temperature, accumulation and surface melt as drivers. Along the coastal margins, we found that the firn thickness and FAC vary strongly in time, whereas in the interior the variations are small. The temporal variations can be split into a rather stable seasonal cycle driven by snowfall, compaction, summer melt and sublimation, and more irregular decadal variations mainly driven by slow snowfall variations. Variations in FAC and firn thickness have a similar phase, and 62 % of the seasonal and 67 % of decadal surface height variability are due to variations in FAC rather than in firn mass, which emphasizes the importance of correcting for seasonal and decadal FAC variations in altimetry studies.

Comparison of simulated surface elevation change with altimetry confirms that the performance of the updated model has improved (a 17 % reduction in error), notably in Dronning Maud Land and Wilkes Land, which is predominantly due to updated climate forcing. However, a substantial trend difference (>10 cm yr⁻¹) remains in the Antarctic Peninsula and Ellsworth Land, which is likely caused by increasing accumulation over the past few centuries, violating the model assumption of a steady state over the simulated period (1979–2020). Correcting for this over the Antarctic Peninsula and Ellsworth Land, by estimating previous climatic conditions from ice core data, substantially improves this (an error reduction of 38 %). In general, uncertainties in the model formulation, climate forcing or climate forcing during the spin-up cause rather small changes in simulated FAC (<5.2 %), while uncertainties in the accumulation forcing and in the accumulation and temperature forcing during the spin-up period can lead to substantial differences in simulated surface elevation change (4 to 19 mm yr⁻¹).