Physics-based SNOWPACK model improves representation of near-surface Antarctic snow and firn density

Estimates of snow and firn density are required for satellite altimetry based retrievals of ice sheet mass balance that rely on volume to mass conversions. Therefore, biases and errors in presently used density models confound assessments of ice sheet mass balance, and by extension, ice sheet contribution to sea level rise. Despite this importance, most contemporary firn densification models rely on simplified semi-empirical methods, which are partially reflected by significant modeled density errors when compared to observations. In this study, we present a new drifting snow compaction scheme that we have 5 implemented into SNOWPACK, a physics-based land surface snow model. We show that our new scheme improves over existing versions of SNOWPACK by increasing simulated near-surface (defined as the top 10 m) density to be more in line with observations (near-surface bias reduction from -44.9 to -5.4 kg m−3). Furthermore, we demonstrate high-quality simulation of near-surface Antarctic snow and firn density at 122 observed density profiles across the Antarctic ice sheet, as indicated by reduced model biases throughout most of the near-surface firn column when compared to two semi-empirical firn densi10 fication models (SNOWPACK mean bias = -9.7 kg m−3, IMAU-FDM mean bias = -32.5 kg m−3, GSFC-FDM mean bias = 15.5 kg m−3). Notably, our analysis is restricted to the near-surface where firn density is most variable due to accumulation and compaction variability driven by synoptic weather and seasonal climate variability. Additionally : , the GSFC-FDM exhibits lower mean density bias from 7 10 m (SNOWPACK bias = -22.5 kg m−3, GSFC-FDM bias = 10.6 kg m−3) and throughout the entire near-surface at high accumulation sites (SNOWPACK bias = -31.4 kg m−3, GSFC-FDM bias = -4.7 kg m−3). 15 However, we found that the performance of SNOWPACK did not degrade when applied to sites that were not included in the calibration of semi-empirical models. This suggests that SNOWPACK may possibly better represent firn properties in locations without extensive observations and under future climate scenarios, when firn properties are expected to diverge from their present state.

a constitutive relationship between stress and strain for snow. However, physics-based approaches are hindered by a variety of factors including computational expense, the requirement to estimate unknown parameters such as surface roughness length (controlling turbulent fluxes and friction velocity) and snow activation energy (controlling snow viscocity), as well as the 60 need to calculate intermediate prognostic variables including thermal conductivity, viscosity, and snow grain shape and size.
Despite these drawbacks, the complex material properties of snow, combined with an acute scarcity of snow and firn density measurements lead us to hypothesize that a physics-based modeling approach is preferred.
In order to improve model representation of Antarctic snow and firn properties compared to semi-empirical models, we compare results from the physics-based snow model, SNOWPACK, forced by hourly weather data from MERRA-2 atmospheric 65 reanalysis (section 2.2) to nine automatic weather stations (AWSs), 55 borehole 10 m depth temperatures, and 122 observed near-surface density profiles for a total of 186 locations across the AIS. First, we describe model setup and a new drifting snow routine designed to improve representation of near-surface (depth ≤ 10 m) snow and firn density (section 2.1). Next, we evaluate SNOWPACK's ability to simulate the surface energy balance by comparing with available surface temperature proxies (section 3.1). We then explore the sensitivity of SNOWPACK simulated density profiles to uncertainties in atmospheric forcing 70 and prescribed snow physics (sections 3.3 -3.4). Next, we compare SNOWPACK to observed density profiles and compare the relative performance of SNOWPACK to two other firn densification models (sections 3.5 -3.8). Finally, we conclude with a discussion of SNOWPACK predicted surface density variability and its implications for satellite altimetry based measurements of ice sheet MB (section 3.9).
2 Methods 75 2.1 Physics-based snow model SNOWPACK Models commonly assume a high (> 250 kg m −3 ) new snow density over Antarctica (Ligtenberg et al., 2011;Groot Zwaaftink et al., 2013) despite observational evidence of occasional relatively low (< 100 kg m −3 ) new snow density (Groot Zwaaftink et al., 2013;Sommer et al., 2018). This difference can be explained by accumulation of snowfall without sufficient wind speed required for drifting snow compaction to occur. In order to account for this, we have implemented a new drifting snow com-80 paction routine into the vertical, one-dimensional physics-based land-surface snow model, SNOWPACK (Bartelt and Lehning, 2002;Lehning et al., 2002a, b), which contrasts to most existing firn models (sections 2.5 -2.6) in that it calculates densification using an overburden stress formulation as opposed to an empirical relationship and explicitly determines snow viscosity by calculating the snow microstructure (e.g. grain size and shape) and temperature. Originally designed as an alpine snow cover model capable of describing snow cover properties, including accumulation, densification, temperature, energy balance, 85 and snow microstructure, SNOWPACK has been applied to both the Greenland and Antarctic ice sheets (Groot Zwaaftink et al., 2013;Van Tricht et al., 2016;Steger et al., 2017;Izeboud et al., 2020;Dunmire et al., 2020). These studies have shown that SNOWPACK is capable of capturing important processes in the ice sheet firn layer, namely accumulation in windy environments, surface meltwater production, and subsequent liquid water retention in the firn. However, these studies have not implemented a polar snow accumulation scheme that allows for initial accumulation of low (< 250 kg m −3 ) density snow and subsequent drifting snow compaction. Additionally, no previous study has evaluated SNOWPACK simulated near-surface firn density across the entire AIS.
Our new drifting snow compaction routine combines an existing alpine new snow density parameterization with a polar snow compaction routine. In particular, we impose the key physical restraint that drifting snow is required for near-surface compaction in addition to overburden stress compaction. This mechanism has previously been proposed in the literature (Brun 95 et al., 1997) and recently confirmed in both laboratory experiments (Sommer et al., 2017) as well as field observations on the AIS (Sommer et al., 2018). In our scheme, new snowfall is assigned a typically low (30 -150 kg m −3 ) density ρ new , Eq. (1), according to the default SNOWPACK alpine new snow density parameterization which is a multiple linear regression derived from observations in the Swiss Alps (Lehning et al., 2002a), and varies as a function of 2 m air temperature T 2m ( • C), snow surface temperature T s ( • C), relative humidity RH (0 -1), and 10 m wind speed U 10 m (m s −1 ).

100
ρ new = 70 + 6.5T 2m + 7.5T s + 0.26RH + 13U 10 m − 4.5T 2m T s − 0.65T 2m U 10 m − 0.17RH U 10 m + 0.06T 2m T s RH Densification then occurs via two distinct processes, a) densification due to overburden stress and b) densification due to drifting snow compaction. Densification due to overburden stress is calculated nearly identically as described in Steger et al. (2017), which adapted SNOWPACK for use in the percolation zone of the Greenland Ice Sheet by tuning the viscosity parameters including the snow activation energy Q s and critical exponent β. Here, we find that the parameter tuning proposed 105 by Steger et al. (2017) leads to significantly overestimated densities (bias > 50 kg m −3 ) in the dry snow zone of Antarctica.
Therefore we revert to original SNOWPACK viscosity parameters by setting Q s and β to 67,000 J mol −1 and 0.7, respectively.
The second mechanism driving densification, i.e. drifting snow compaction, occurs when the friction velocity u * (m s −1 ) exceeds the snow threshold friction velocity u * th (m s −1 ), the minimum friction velocity at which surface snow grains are mobilized by the wind (Fig. 1). In our implementation of SNOWPACK, u * is estimated by scaling hourly averaged 10 m 110 wind speeds from the MERRA-2 atmospheric reanalysis (section 2.2) using a logarithmic wind profile and stability corrections (Michlmayr et al., 2008) with a roughness length, z 0 , of 2 mm. Meanwhile, u * th , Eq.
(2), is calculated as a function of of snow microstructural properties internally determined by SNOWPACK, including snow grain sphericity SP (0 -1), radius r g (m), bond radius r b (m), and coordination number N 3 (Lehning and Fierz, 2008).
In Eq. (2), ρ i is the density of ice (917 kg m −3 ), ρ a is the density of air (1.1 kg m −3 ), g is the gravitational acceleration (9.8 m s −2 ), σ is a reference shear strength set to 300 P a, while constants A and B are set to 0.02 and 0.0015 respectively.
When u * exceeds u * th , a saltation mass transport rate Q (kg m −1 s −1 ) is calculated following Lehning and Fierz (2008) and then scaled to a saltation mass flux Φ (kg m −2 s −1 ), Eq. (3), by dividing Q by a characteristic horizontal length scale L.
120 L can be interpreted as a fetch length and characteristic horizontal length scale over which the originally upwind and now mobilized snow particles, which make up the saltation mass flux Φ, have been eroded from the snow surface. We choose a magnitude for L of 10 m, as this is a typical horizontal length scale of wind erosion features including sastrugi and barchan dunes (Filhol and Sturm, 2015). Given a lack of direct observations, L can effectively be considered a tuning parameter, whose magnitude is inversely proportional to Φ. Note that as suspension of drifting snow is not considered in the saltation 125 model (Lehning and Fierz, 2008), the saltation mass transport rate Q may underestimate the total mass flux in the saltation and suspension layers. However, the scaling of the mass transport rate by the poorly constrained fetch length L to obtain the saltation mass flux Φ could be considered a way to include both suspension and saltation in the simulations.
Following erosion, the drifting snow mass flux is redeposited at the snow surface with an updated density ρ redeposit (kg m −3 ), Eq. (4), parameterized according to field measurements of surface snow deposited during drifting snow events (Groot Zwaaftink 130 et al., 2013, equation 1). Note that, in contrast to Groot Zwaaftink et al. (2013), which uses 100-hour average wind speeds, we calculate ρ redeposit using hourly mean MERRA-2 10 m wind speeds. We implement this distinction in order to better resolve the temporal evolution of redeposited snow density during ephemeral (i.e. shorter than 100 hours) drifting snow events. hourly mean 2 m air temperature, relative humidity, 10 m wind speed, incoming shortwave and longwave radiation (ISWR and ILWR), and precipitation rate. At the bottom of the simulated snowpack, we apply the MERRA-2 1980 -2017 mean annual surface temperature as a thermodynamic boundary condition. We choose MERRA-2 because it provides a low release latency (approximately one month) and publicly available description of the state of the atmosphere, whereas regional climate models 140 are not always regularly updated. This advantage, although not strictly necessary for this study as we focus on a past period (1980 -2017), would be advantageous for future applications where timely estimates of snow properties are required, for example rapid interpretation of satellite imagery and altimetry, or determining the status of field assets (e.g. weather stations).
Additionally However, it also smooths out high frequency discrepancies that may be important when evaluating instantaneous simulated 155 density profiles. The AWSs measure temperature, relative humidity, wind speed and direction, air pressure, and the full radiation balance, i.e. incoming and reflected short wave radiation, and incoming and outgoing longwave radiation. Data from these stations are presented in Reijmer (2002) and Jakobs et al. (2020) and have been previously used to evaluate remote sensing retrievals (Trusel et al., 2013), ice core paleoclimate records , and climate models (van Wessem et al., 2018). According to our analysis and consistent with the findings of Lenaerts et al. (2017) and Gossart et al. (2019), MERRA-2 160 well captures observed 2 m air temperature, relative humidity, and wind speed, but significantly underestimates both ISWR and ILWR. We calculated an average MERRA-2 bias across all nine AWSs of -15.1 W m −2 (corresponding to 19.4%) and -16.9 W m −2 for ISWR and ILWR, respectively (Fig. A1). In order to reduce this bias in incoming radiation and thus better

SNOWPACK model initialization
In our scheme, new snow layers are added on top of the snow column when precipitation is prescribed by the atmospheric forcing, in steps of 2 cm. Layers are initialized with a density given by Eq. 1 when they originate from precipitation. When snow layers are eroded and redeposited, the layers originating from drifting snow are initialized with a density given by Eq. 4.
Initial grain size for all newly added layers is 0.2 mm (Groot Zwaaftink et al., 2013). There are two sets of microstructural 175 properties for grain shape (dendricity and sphericity), for high and low wind speed, respectively (Groot Zwaaftink et al., 2013).
Note that precipitation is treated before assessing snow erosion, such that low density snow from precipitation can erode immediately when conditions allow. To reduce computational costs, a sophisticated snow layer merging scheme merges layers with very low ice content due to sublimation or melt and layers with similar properties (density, water content, grain size, and grain shape parameters). The criteria for layer merging are relaxed with depth, to allow for more aggressive layer merging. At 180 10 m depth, typical layer spacing is around 10 -20 cm. Near the surface, layers can be split to maintain a vertical resolution of a few cm in order to numerically represent steep temperature gradients.
In order to ensure a realistic representation of snow and firn properties throughout the entire near-surface, we complete a SNOWPACK model spinup such that simulated snow depth is 10 m at all sites before comparison with observations or other SNOWPACK simulations. We choose 10 m in order to resolve seasonal variability in snow and firn characteristics as well as

IMAU-FDM model
To aid in the evaluation of SNOWPACK modeled near-surface density, we also compare with the widely used Institute for Marine and Atmospheric research Utrecht firn densification model (IMAU-FDM v1.1) (Ligtenberg et al., 2011;Kuipers Munneke et al., 2015;Ligtenberg et al., 2018). The IMAU-FDM is a semi-empirical firn densification model designed to represent snow and firn cover processes including densification, meltwater refreezing and percolation, and surface height change. Gridded

200
IMAU-FDM density profiles are available at 27 km horizontal, 4 cm vertical, and 30 day temporal resolution with atmospheric forcing provided by the regional climate model RACMO 2.3p2 (van Wessem et al., 2018). In contrast to SNOWPACK, which relies on an overburden stress compaction scheme, IMAU-FDM uses a calibrated semi-empirical dry snow densification scheme based on Arthern et al. (2010), which describes densification as a function of density as well as annual average accumulation and temperature. In further contrast to SNOWPACK, the IMAU-FDM parameterizes new snow density as a function 205 of annual average, rather than hourly, meteorology and currently includes no post-deposition mechanism to increase surface snow density due to drifting snow processes.

GSFC-FDM model
In addition to IMAU-FDM we also compare SNOWPACK to the NASA Goddard Space Flight Center firn densification model (GSFC-FDMv1) which provides simulated firn properties over the past 40 years (1980 -2019) for both the Greenland and 210 Antarctic ice sheets . GSFC-FDM uses the Community Firn Model, a modular, open-source framework for Lagrangian modeling of several firn and firn-air related processes within a single column (Stevens et al., 2020). The GSFC-FDM simulations are forced by an enhanced resolution hybridized MERRA-2 that was developed by exploiting a 15-year, 12.5 km resolution offline MERRA-2 reanalysis. The hybridized MERRA-2 forcing retains the spatial gradients in the highresolution reanalysis while maintaining the temporal variations from the original MERRA-2 variables (see Medley et al. (2020) 215 for complete details). The dry snow and firn compaction model, based on Arthern et al. (2010), was calibrated to observed depth-density profiles from both Greenland and Antarctica. A simple initial density scheme was implemented based on mean annual MERRA-2 climate, which provides a spatially variable initial density that does not, in contrast to SNOWPACK, vary in time or vary due to drifting snow processes.

SNOWPACK, IMAU-FDM, and GSFC-FDM comparison with density observations 220
We evaluate skill for all three firn models in representing surface and near-surface snow and firn density, defined here as the average density between depths of 0 -1 m and 0 -10 m, respectively, by comparing with community sourced and publicly available density profiles from the Surface Mass Balance and Snow on Sea Ice Working Group (SUMup) data set (Albert, 2007;Medley et al., 2013;Montgomery et al., 2018). We first identified all available Antarctic SUMup density observations in the near-surface and then retrieved spatially and temporally consistent SNOWPACK, IMAU-FDM, and GSFC-FDM model output.

225
For all models, we retrieve the simulated density profile whose time step is closest to that of the observed profile. SNOWPACK, IMAU-FDM, and GSFC-FDM report simulated density profiles every 1, 30, and 5 days respectively. This classification yielded 122 unique observed profiles (Fig. 2) that are located primarily on the grounded ice sheet, where surface melt is limited (< 50 mm w.e. yr −1 ), absent, or very rare (Trusel et al., 2013). Observations and models report different vertical resolution density, therefore we compute the mean and standard deviation of all 122 density profiles at ten 1 m thick vertical levels, beginning at Since the MERRA-2 bias in ISWR and ILWR is determined using only 9 AWSs located primarily in Dronning Maud Land (section 2.2), we were initially concerned with the bias' spatial representativeness. However, following bias correction, the statistically significant reduction in both RMSE and bias (p < 0.01) magnitude when compared to 10 m depth temperatures,  3.3 SNOWPACK density profiles sensitivity to atmospheric forcing uncertainty.
Although we have performed an evaluation of MERRA-2 SEB, and introduced a bias correction for ISWR and ILWR, un-understand the effect of these uncertainties on SNOWPACK simulated density, we perform an ensemble of SNOWPACK sim- "Enhanced Compaction", and 3) the default alpine SNOWPACK snow physics, which includes only a modest wind-enhanced 295 surface compaction using the strain rate (Bartelt and Lehning, 2002;Lehning et al., 2002a, b), denoted "Default".
From our analysis, it is immediately clear that the default SNOWPACK snow physics significantly underestimate observed surface density (> 150 kg m −3 ), suggesting that the wind effect described in the new snow density parameterization (Eq. 1) and the increased strain rate in this variant are insufficient for drifting-snow dominated polar environments. However, in the case of "Default" snow physics, this underestimation of density decreases with depth, which can be explained by the 300 reduced viscosity of relatively lower density snow, and therefore more rapid densification due to overburden stress. For the "Enhanced Compaction" simulation, we see similar behavior, however modeled surface density errors are reduced compared to the "Default" simulation. Compared to "Redeposit", the "Enhanced Compaction" simulation significantly underestimates near-surface density. This demonstrates that the enhanced wind compaction mechanism via the strain rate is not sufficient to capture observations. Instead, describing compaction of surface snow due to wind using wind-driven erosion and subsequent 305 redeposition at a higher density ("Redeposit"), which resembles the experimentally determined mechanism (Sommer et al., 2017), adequately captures observed near-surface density.

Surface snow density
A comparison between both SNOWPACK, IMAU-FDM, and GSFC-FDM simulated surface snow density, defined as the top meter, with 79 unique observations is shown in Figure 7. Observed densities range from 272 -507 kg m −3 , with a mean of 362 310 kg m −3 and standard deviation of 39 kg m −3 . Since it is known that meteorological conditions including annual accumulation and temperature influence Antarctic snow and firn density (Herron and Langway, 1980), we tested this hypothesis. We find a modest relationship between MERRA-2 1980 -2017 mean annual SMB and observed surface snow density (p < 0.01, R 2 = 0.23) as well as a significant, but weak correlation between elevation and observed surface snow density (p < 0.01, R 2 = 0.16), likely due to the relationship between elevation and surface climate. Note that we find no significant correlation between 315 MERRA-2 1980 -2017 mean wind speed and observed surface density (p = 0.14, R 2 = 0.03) but do in fact find a significant, albeit weak relationship between 1980 -2017 maximum wind speed and surface density (p < 0.01, R 2 = 0.09). Because drifting snow compaction is known to partially control snow density on daily to hourly timescales (Sommer et al., 2018), the lack of a significant relationship between mean annual wind speed and surface snow density combined with the significant relationship between maximum wind speed and surface density indicates the importance of resolving drifting snow compaction with high 320 temporal resolution (daily to hourly) meteorological forcing as opposed to annual means or climatology. According to our analysis, SNOWPACK better captures the range in observed surface density when compared to IMAU-330 FDM, as evidenced by a lower magnitude bias, higher R 2 , and a linear fit slope closer to unity, but the typical error is larger, as expressed by a slightly larger RMSE (Table 1). Likewise, when compared to GSFC-FDM, SNOWPACK exhibits a lower surface density bias magnitude and has a linear slope closer to unity, while it shows a larger RMSE and lower R 2 . Additionally, we calculate a p-value of 0.06 for a two-sided t-test in which the null hypothesis states that SNOWPACK and IMAU-FDM have identical average biases. This information combined with SNOWPACK's larger RMSE compared to IMAU-FDM and 335 GSFC-FDM, and smaller R 2 compared to GSFC-FDM, leads us to conclude that neither model is superior, and instead that all three models perform comparably with regard to surface density. Note that all three models perform particularly poorly when performance at high observed surface densities is still unknown, but we can come up with at least three potential explanations.

340
First, we speculate that high surface densities are found in areas with locally steep topography. However, we find no significant relationship between surface slope derived from a 1 km horizontal resolution digital elevation model (Helm et al., 2014) and modeled surface density bias (not shown), so this mechanism is an unlikely candidate to explain the bias. Second, this model degradation at high observed surface densities could potentially be explained by an underestimation of new snow density and/or the frequency/intensity of simulated drifting snow compaction due to the presence of local meteorological phenomena 345 not captured by MERRA-2 or RACMO2. Lastly, we cannot rule out the possibility of larger errors in the observational data for densities above 400 kg m −3 . For example, some surface density observations, particularly those from deep firn and ice cores, may not be representative of an undisturbed snow surface and could therefore report an artificially high surface density.
pieces, therefore confounding density measurements.

Near-surface snow and firn densification
In a comparison with 122 observed density profiles, SNOWPACK exhibits a lower bias compared to IMAU-FDM for the entire near-surface, and a lower bias compared to GSFC FDM from 0 -7 m depth (Fig. 8). correctly predicts observed variability in mean density, as shown by similar standard deviations between SNOWPACK and observations ( Fig. 8, panel A). Alternatively, the IMAU-FDM and GSFC-FDM both underestimate observed mean density variability at the surface, but converge towards the observed variability with increasing depth. Note that 95% confidence intervals on the SNOWPACK mean bias contain zero in the top 7 m (Fig. 8), thus we conclude that SNOWPACK is consistent with observations at these depths. However, the IMAU-FDM and GSFC-FDM mean bias 95% confidence intervals never 365 contain zero and are therefore not consistent with observations. Despite these contrasting conclusions, SNOWPACK's mean bias confidence intervals occasionally overlap with that of the IMAU-FDM and GSFC-FDM (e.g. 0 -1 m and 8 -10 m) in the case of IMAU-FDM) and therefore should not be considered uniformly statistically indistinguishable.
To test for the potential effect of compensating biases on our analysis of near-surface densification, we partition the 122 ob- we find that SNOWPACK, when compared to IMAU-FDM, shows lower absolute density biases throughout the entire nearsurface firn column, from the surface down to 10 m depth. Compared to GSFC-FDM, SNOWPACK exhibits a smaller absolute near-surface density bias only from 0 -8 m and 9 -10 m at low accumulation sites. For high SMB sites, SNOWPACK mean 375 biases range from -55.9 kg m −3 at 8 -9 m to -3.9 kg m −3 at 0 -1 m, while for low SMB sites, mean biases are generally Error bars represent 95% confidence intervals on bias means. reduced in magnitude and range from -14.3 kg m −3 at 9 -10 m to 8.3 kg m −3 at 5 -6 m. Given SNOWPACK's mean bias reduction, when compared to IMAU-FDM, and comparable peformance relative to GSFC-FDM, we have demonstrated our physics-based modeling approach is capable of reliably capturing Antarctic near-surface snow and firn density. Consistent with our analysis of all profiles in Figure 8, SNOWPACK, as measured by the mean bias 95% confidence intervals, is gen-380 erally consistent with observations and statistically distinguishable from IMAU-FDM and GSFC-FDM at low accumulation sites. However, at high accumulation sites, SNOWPACK exhibits significantaly worse performance as indicated by a strongly negative bias, and confidence intervalls that do not contain zero, particularly below 6 m depth.  application, amd could yield an advantage over calibrated firn models, where the application for non-calibrated sites is not guaranteed to yield good performance.

SNOWPACK and semi-empirical models hierarchical complexity
Although all three models are presented on similar footing, it is important to contextualize their performance by noting the different original purposes for developing SNOWPACK, IMAU-FDM, and GSFC-FDM, as well as well as their different level of 405 complexity in representing physical processes. The IMAU-FDM and GSFC-FDM are relatively simplified semi-empirical firn models designed to represent spatial and temporal evolution of firn density and surface height. On the other hand, SNOWPACK is a higher order complexity physics-based land-surface snow model originally intended for simulating snow microstructural properties relevant for avalanche formation in seasonally snow covered terrain. IMAU-FDM and GSFC-FDM are usually used to study the entire ice sheet firn column, which can be greater than 100 m thick, while SNOWPACK is typically used to simu-410 late seasonal snowpacks a few meters thick. Despite the much thicker simulated firn column in IMAU-FDM and GSFC-FDM (40 -120 m) when compared to our imposed restriction for SNOWPACK to the near-surface (i.e. depths 0 -10 m), SNOW-PACK is considerably more computationally expensive to run. It is therefore important to acknowledge the inherent trade off between process representation and cost, particularly for simulations which require long spinups or involve long time frames, e.g. paleoclimate applications or future climate scenarios. days) and GSFC-FDM (5 days), can resolve processes relevant for surface snow density which act on timescales of less than 30 and 5 days, respectively (e.g. snowfall and drifting snow compaction).

SNOWPACK simulated surface density variability
SNOWPACK predicts a wide range of simulated surface densities (mean of top meter), ranging from 230 -684 kg m −3 with a mean of 370 kg m −3 and standard deviation of 64 kg m −3 across 186 simulations (Fig. 11a). Furthermore, SNOWPACK 435 exhibits large seasonal surface density variability (mean magnitude = 53 kg m −3 ), when measured as the maximum difference between summer and winter surface density found in the period 1980 -2017 (Fig. 11b). Seasonal surface density variability is primarily driven by seasonal cycles in temperature, accumulation rate, and compaction rate. However, superimposed onto these processes at even shorter timescales (e.g. hourly to weekly) are individual accumulation events and snow microstructural evolution (Fig. 11c). Because satellite altimetry relies on instantaneous measurements of surface height, rather than temporal 440 means, firn models used to reliably convert changes in volume into mass must capture high frequency variability in firn density and surface height.

Conclusions
Accurate snow and firn density models are required to reliably determine ice sheet mass balance using satellite altimetry.
However, firn densification models currently used in altimetry studies do not resolve observed temporal increases in surface 445 snow density under the influence of wind. In this study, we demonstrate improved simulation of Antarctic near-surface (depths ≤ 10 m) snow and firn density upon implementation of a new drifting snow compaction routine into SNOWPACK, a detailed, physics-based land surface snow model. In particular, we show that when compared to two other semi-empirical firn densification models (IMAU-FDM and GSFC-FDM), SNOWPACK exhibits a lower mean snow and firn density bias throughout most of the near-surface at 122 observed density profiles across the Antarctic ice sheet. Despite this improvement, SNOWPACK 450 generally underestimates density, with a mean surface density bias of -8.2 kg m −3 and mean near-surface biases ranging from -0.6 kg m −3 at 2 -3 m depth to -25.3 kg m −3 at 8 -9 m depth. Meanwhile, IMAU-FDM exhibits a mean surface density bias of -20.4 kg m −3 and mean near-surface biases between -20.4 kg m −3 at depths 0 -1 m to -40.1 kg m −3 at depths 8 -9 m while GSFC-FDM shows a mean surface density bias of 20.4 kg m −3 and mean near-surface biases from 8.5 kg m −3 at depths 8 -9 m to 23.7 kg m −3 at 2 -3 m. SNOWPACK exhibits a lower density bias than IMAU-FDM throughout the entire 455 near-surface, however the GSFC-FDM shows a lower magnitude bias than SNOWPACK below 7 m depth and throughout the entire near-surface at high accumulation sites. Note that our analysis is limited to the top 10 m in order to focus on the most dynamic and variable part of the Antarctic firn layer. Because SNOWPACK is a physics-based model, extensive model tuning in order to fit observations is not required. By analyzing the simulations, excluding the sites used for calibration of the semi-eprical models, we found SNOWPACK to also have the lowest absolute bias. IMAU-FDM showed a lower bias when ex-