Observations and simulations of new snow density in the drifting snow-dominated environment of Antarctica

Abstract

Owing to drifting snow processes, snow accumulation and surface density in polar environments are variable in space and time. We present new field data of manual measurements, repeat terrestrial laser scanning and snow micro-penetrometry from Dronning Maud Land, Antarctica, showing the density of new snow accumulations. We combine these data with published drifting snow mass flux observations, to evaluate the performance of the 1-D, detailed, physics-based snow cover model SNOWPACK in representing drifting snow and surface density. For two sites in East Antarctica with multiple years of data, we found a coefficient of determination for the simulated drifting snow of r² = 0.42 and r² = 0.50, respectively. The field observations show the existence of low-density snow accumulations during low wind conditions. Successive high wind speed events generally erode these low-density layers while producing spatially variable erosion/deposition patterns with typical length scales of a few metres. We found that a model setup that is able to represent low-density snow accumulating during low wind speed conditions, as well as subsequent snow erosion and redeposition at higher densities during drifting snow events was mostly able to describe the observed temporal variability of surface density in the field.

Introduction

Ice-sheet surface mass balance (SMB) comprises the sum of mass fluxes that act at the ice-sheet surface (Lenaerts and others, 2019). Mass accumulation is governed by precipitation, and deposition and condensation of water vapour. Locally, accumulation is impacted by drifting (lowermost 2 m of the atmosphere) and blowing snow (above 2 m), eroding the snow surface and creating new snow accumulations (Gallée and others, 2001; Déry and Yau, 2002; Libois and others, 2014). We use the term accumulation here for any mass added to the firn layer, irrespective of its source (drifting snow, precipitation, water vapour flux) and deposition for accumulations originating from drifting snow.

Snow transport by wind has a demonstrable impact on the spatial variability of Antarctic SMB (Lenaerts and others, 2012b), which has been found for the kilometre scale (Dattler and others, 2019), sub-kilometre scale (Kausch and others, 2020), as well as metre scale (Picard and others, 2019). These studies also demonstrated that total annual snow erosion and deposition by wind can locally exceed annual precipitation. For example, Picard and others (2019) showed that even the height of individual snow dunes formed during drifting and blowing snow can exceed annual accumulation over the East Antarctic plateau, such that patches of snow at the surface differ markedly in age. The occurrence of drifting and blowing snow also varies in time, governed by both atmospheric boundary layer processes (Huang and Wang, 2016; Paterna and others, 2016) as well as snow micro-structural properties (Schmidt, 1980; Lehning and others, 2000; Clifton and others, 2006). This adds complexity in assessing quantity and frequency of snow transport by wind.

The density of near-surface snow layers is, among other factors, impacted by drifting snow conditions (e.g. Brun and others, 1997; Sommer and others, 2018). With increased revisit frequencies, spatial resolution and vertical accuracy of satellite repeat altimetry, understanding the temporal and spatial variability of the density of newly formed snow depositions originating from drifting and blowing snow is crucial to translate changes in surface elevation to changes in mass (Shepherd and others, 2012; Montgomery and others, 2020; Zwally and others, 2021). It also has been repeatedly shown that determining the density of fresh snow accumulations is crucial to obtain a correct density profile with depth (Herron and Langway, 1980; Arthern and Wingham, 1998; Jay Zwally and Jun, 2002; Ligtenberg and others, 2011; Kuipers Munneke and others, 2015). For ice shelves, accurate density–depth profiles are important to determine the available pore space to store future surface meltwater (e.g. Kuipers Munneke and others, 2014; Trusel and others, 2015; Lenaerts and others, 2017). Drifting snow also frequently occurs in the Arctic tundra (e.g. Pomeroy and others, 1997), where the variations in density impact the thermal regime (e.g. Domine and others, 2012). Yet, field measurements of snow density are typically point measurements (Montgomery and others, 2018), insufficiently capturing the variability of density in space and time. For example, Weinhart and others (2020) and Wever and others (2021) show density variability at the surface of 60–100 kg m⁻³ over several tens of metres for the East Antarctic plateau, and sea ice, respectively. Over the scale of 100 km, Sugiyama and others (2012) and Weinhart and others (2020) report near surface density to vary up to 100 kg m⁻³ for this region. Sommer and others (2018) provide measurements from Antarctica of very low snow density at the surface of ~100 kg m⁻³, which existed for a few days after which drifting snow eroded this layer. Groot Zwaaftink and others (2013a) present measurements showing that intercepted solid precipitation on boards 1 m above the snow surface exhibited densities in a similar range (10–150 kg m⁻³), while the density of the uppermost 10 cm of the snow surface was found to be in the 300–400 kg m⁻³ range. Since the impact of drifting snow is strongest near the surface, those results were interpreted as the effect drifting snow has on the surface density, rather than wind speed alone.

New snow density parameterizations used in land surface schemes and snow and firn models such as Crocus (Vionnet and others, 2012), SNOWPACK (Lehning and others, 2002a; Schmucki and others, 2014), SnowTran-3D (Liston and others, 2007) and the Community Land Model, part of the Community Earth System Model (van Kampenhout and others, 2017), typically include meteorological parameters, such as temperature and relative humidity. This in recognition that atmospheric conditions influence snowfall hydrometeor size and shape, which in turn are found to control new snow density (e.g. Nakaya and others, 1958; Judson and Doesken, 2000; Colle and others, 2014; Ishizaka and others, 2016). Most parameterizations also include a term describing the influence of wind speed on new snow density, since with increasing wind speed, mechanical destruction of grains allow for deposition of denser snow (e.g. Vionnet and others, 2012; Schmucki and others, 2014; Comola and others, 2017; Morin and others, 2020). However, parameterizations based on fresh snow density measurements from alpine measurement sites were found to lead to poor performance under typical ice-sheet conditions (Groot Zwaaftink and others, 2013a; Steger and others, 2017). This likely is a result of the underrepresentation of high wind speed and low temperature regimes in those datasets. Therefore, parameterizations tailored to the polar regions have been developed (Liston and others, 2007; Groot Zwaaftink and others, 2013a).

Drifting and blowing snow resulting in new depositions frequently occurs in the absence of precipitation (Gossart and others, 2017; Palm and others, 2018), where the source of the drifting snow particles is erosion from the upwind snow surface. This implies that only considering the effect of wind on new snow density during periods of precipitation is insufficient to describe the full influence of wind on surface density. It is also important to note that recent wind tunnel experiments have demonstrated that densification due to wind only occurs under the presence of drifting snow (Sommer and others, 2017). More specifically, the freshly formed deposits exhibit higher density, whereas snow that is not mobilized does not undergo increased compaction from wind. This result has also been validated in the Antarctic field (Sommer and others, 2018). These findings suggest that the saltation process, and the mechanical destruction and subsequent sintering of drifting snow particles, is the governing process for compaction, rather than the previously proposed mechanism of densification of existing snow layers in the top of the snowpack (Craven and Allison, 1998).

To describe the process of wind-driven compaction of already deposited snow, detailed physics-based models such as CROCUS and SNOWPACK include a wind compaction term when calculating compaction in near surface layers (Brun and others, 1997; Groot Zwaaftink and others, 2013a). With increasing wind speed, the densification rate in surface layers increases (Craven and Allison, 1998). For 1-D simulations, describing compaction by wind via snow settling may be adequate, however, when models are used to explicitly treat the spatial redistribution of snow by wind in wind-dominated environments (Lenaerts and van den Broeke, 2012a; Liston and others, 2018; Kausch and others, 2020; Amory and others, 2021), it may be critically important to allow for erosion of snow layers from the surface to distinguish drifting and no-drifting cases. We furthermore note that including wind speed in both the new snow parameterization and the compaction rate can be considered accounting for the same effect twice.

To better account for the role of drifting snow on surface density, Groot Zwaaftink and others (2013a) developed an event-driven deposition scheme, where snow is only deposited when the 100 h moving average wind speed is between 4 and 7 m s⁻¹. This approach, however, prohibits the model from depositing new snow with densities below 150 kg m⁻³, which is in contradiction with reports from the field. That motivated the most recent development by Keenan and others (2021) where snow is able to be eroded from the surface, and redeposited with higher density, when the model diagnoses drifting snow to occur. This approach is able to span the range from low fresh snow density when wind speeds are too low for drifting snow to occur (typically <5 m s⁻¹), up to high density typically found in wind-dominated environments. Although Keenan and others (2021) show the performance of the scheme with respect to near surface density and 10 m density profiles, we focus here particularly on the representation of drifting and blowing snow conditions, and compare detailed in situ measurements of snow density of freshly formed accumulations with simulated density during these conditions.

In the first part of Section ‘Data and methods’, we describe already published drifting snow field data we used, as well as newly obtained field data of snow accumulations and snow density. Section ‘Data and methods’ concludes with a description of the SNOWPACK model and the tested model setups, as well as the atmospheric forcing data. In Section ‘Results and discussion’, we discuss first the performance of SNOWPACK in diagnosing drifting snow, followed by observations of snow accumulation and a comparison of observed and simulated snow density for the field sites. Section ‘Sensitivity study’ investigates the sensitivity of SNOWPACK calculated drifting snow mass fluxes to poorly constrained parameters in the model.

Data and methods

Field sites

Observational data used here were collected at four field sites. The first two field sites are drift stations D17 (66.7° S, 139.9° E, 450 m above sea level (a.s.l.), 10 km from the coast) and D47 (67.4° S, 138.7° E, 1560 m a.s.l., 110 km from the coast) near the Dumont d'Urville base in Adélie Land, East Antarctica (Amory, 2020). Adélie Land is situated in a confluence zone of katabatic flow over steep slopes near the coast (see Fig. 1), making it among the windiest locations of Antarctica (Parish and Bromwich, 2007) with frequent drifting snow events (Amory, 2020).

Fig. 1. Locations of the study sites on the Antarctic ice sheet, including detailed maps showing the positions of the study sites on the REMA digital elevation model (Howat and others, 2019).

The two other field sites are located in Dronning Maud Land, East Antarctica. The first of these is located near the Princess Elisabeth Antarctica station (PEA) at 71.939° S, 23.315° E and ~1350 m a.s.l., where two drift stations were installed. The second field site is situated at Hammarryggen (HAM), an ice rise at 70.502° S, 21.874° E and ~360 m a.s.l., located near the Roi Baudouin Ice Shelf (Kausch and others, 2020), ~170 km away from PEA. As for stations D17 and D47, HAM is located in a confluence zone, exposing it to high wind speeds and frequent drifting and blowing snow (Kausch and others, 2020). In contrast, the PEA station is sheltered by the Utsteinen Nunatak, and, on a larger scale, by the Sør Rondane Mountains (see Fig. 1). Even though PEA experiences lower drifting and blowing snow occurrence than its surroundings, Gossart and others (2017) have demonstrated that blowing snow occurs here ~13% of the time.

Note that following the general definition of drifting snow being transported by the wind in the lowermost 2 m of the atmosphere and blowing snow above that layer, we refer to the stations measuring drifting snow as drift stations, since the measurements only cover the drifting snow layer.

Drifting snow measurements

The two sites D17 and D47 were instrumented with second-generation IAV Engineering acoustic FlowCapt^TM sensors (Trouvilliez and others, 2014). Each FlowCapt^TM sensor is 1 m long, and each site has two sensors stacked on top of each other, covering approximately the lowest 2 m of the drifting snow layer. A snow depth sensor mounted at the stations measured local snow depth changes, which can lead to varying burial of the drifting snow sensors over time. To convert the drifting snow mass fluxes f ₁ and f ₂ (kg m⁻² s⁻¹) from the two sensors, respectively, to the drifting snow mass transport F _2 m (kg m⁻¹) in the lowest 2 m of the atmosphere, we multiplied the mass fluxes by the exposed length of each sensor (h ₁ and h ₂, respectively) and the analysis time step Δt (s). The mass transport from both sensors was summed, applying the same correction for the burial of the FlowCapt^TM sensors as outlined in Amory (2020), which assumes a constant particle flux over the maximum possible measurement height of 2 m:

(1)$$F_{2\, {\rm m}} = \left(\,f_{{\rm 1}} h_{{\rm 1}} + f_{{\rm 2}} h_{{\rm 2}} \right){2\over h_{{\rm 1}} + h_{{\rm 2}}}\Delta t$$

As in Amory (2020), drifting snow mass fluxes below 10⁻³ kg m⁻² s⁻¹ were set to zero. The observed wind speed height also varied over time due to accumulation and height adjustments of the instruments. We constructed a time series of 10 m wind speed from the stations by assuming a logarithmic wind profile in a neutrally stable atmosphere with the roughness length parameterized as a function of air temperature as proposed in Eqn (7) in Amory and others (2017).

The long-term drifting snow measurements at sites D17 and D47 are published in Amory and others (2020). For site D47, measurements are available between 9 January 2010 and 27 December 2012. Site D17 was installed on 3 February 2010 and is still operational. Data included here end on 31 December 2018.

At PEA, two drifting snow masts were installed 350 m apart, a few kilometres west of the PEA station. These stations measured, among others, wind speeds at two levels (~1.1 and 3.6 m above the snow surface), air temperature and relative humidity (at ~1.1 m above the snow surface) and drifting snow, which was measured by Niigata Electric Snow Particle Counters (SPC-95) recording the number and size of saltating particles (Sato and others, 1993). The SPCs can rotate and vanes on the devices ensure that the sensors remain aligned with the wind direction. Sommer and others (2018) describe the data processing for these field sites in more detail. The measurement technique using SPCs provides a point measurement in the drifting snow layer, whereas the FlowCapt^TM sensors used at D17 and D47 provide an integrated value over the length of the exposed tube. For comparison with simulated saltation mass fluxes, we convert the observed drifting snow mass flux f _SPC (kg m⁻² s⁻¹) from the SPC to the drifting snow mass transport in the saltation layer F _salt (kg m⁻¹), by multiplying the mass flux with a typical depth of the saltation layer (h _salt) of 10 cm and the analysis time step:

(2)$$F_{{\rm salt}} = f_{{\rm SPC}} h_{{\rm salt}} \Delta t$$

As discussed in Sommer and others (2018), the drifting snow measurements by the SPCs are expected to underestimate the mass flux, since they were positioned 13–24 cm above the snow surface, which is above the transition from the saltation to the suspension layer (Nemoto and Nishimura, 2004). The data from the drifting snow stations are published in Wever and others (2018). SPC data are available between 15 December 2016 and 22 December 2017 with some data gaps. A 10 m wind speed time series was constructed in a similar procedure as for the D17 and D47 sites, but assuming a roughness length of 7.5 × 10⁻³ m, as reported for this field site (Sommer and others, 2018), instead of using the parameterized roughness length.

As shown in Table S1 in the Supplementary material, we also include wind speed data recorded by the automated weather station (AWS) installed at PEA by the Institute for Marine and Atmospheric research Utrecht (IMAU), the Netherlands (Reijmer and Oerlemans, 2002; Jakobs and others, 2020), to extend the length of the time series. When referring to in situ wind speed, we use the IMAU AWS data for the period from 13 January 2010 until 15 January 2018. Starting 10 December 2016, wind speed data from the drifting snow masts are used when available. For the detailed study period between 10 December 2016 and 15 January 2017, only the data from the drifting snow masts are used, with nearest neighbour interpolation when data gaps are present.

Density measurements

The snow microstructure at HAM and PEA was surveyed with a SnowMicroPen (SMP) version 4, a high precision, constant-speed penetrometer (Schneebeli and Johnson, 1998; Proksch and others, 2015). The general relationship between the snow density and the SMP measurement is typically described by the median force F and distance between the structural elements L over a segment. We used the SnowMicroPyn software (WSL Institute for Snow and Avalanche Research SLF, 2018), version 1.0.1, with default settings for the ‘Proksch 2015’ algorithm (Proksch and others, 2015), which uses a segment length of 1.25 mm, for processing the SMP measurements. Since SMP version 4 is a newer version than the one used for the calibration developed by Proksch and others (2015) and was operated by us with different settings, we established a calibration based on manual density observations using a 100 cm³, 3 cm high box cutter (Proksch and others, 2016, see Section S1 in the Supplementary material for details). When acquiring a measurement, the instrument is placed on top of the snow surface using a frame, and the measurement tip travels vertically ~5–10 cm through the air before reaching the snow surface. The software's default algorithm for the detection of the snow surface was used to remove that first part of the measurement.

On 3 days in December 2016, fresh snowfall at PEA allowed for SMP measurements of light snow with densities below 100 kg m⁻³. In these cases, only the new snow layer, which was on top of a hard old snow layer, was surveyed by inserting the box cutter vertically in the new snow layer until it reached the old, noticeably much harder, snow surface and using a measurement of new snow depth at the location of the box cutter to calculate density. The manual snow pits and the SMP data, including the processing script, are published in Wever and others (2022).

Terrestrial laser scanning

The spatial variability of snow accumulations at PEA and HAM was determined with high precision using repeat terrestrial laser scanning (TLS). Sommer and others (2018) describe the scans obtained at PEA, and in Section S2 in the Supplementary material, we describe the scans obtained at HAM. At HAM, scans were acquired on 4 days and differences between scans are shown in Figure 2. The scans from 27 December, 4 January and 11 January are of high quality, whereas during post-processing, the scans from 2 January were found to have exhibited small tilt during the scan. This is visible by a transition from positive to negative snow depth change in the bottom right of Figures 2a, b. Figure 2c therefore additionally shows the snow depth change between 27 December and 4 January.

Fig. 2. Snow depth changes calculated from laser scans obtained on 4 days in the 2018–19 field season on the HAM ice rise, in an arbitrary local coordinate system. Areas that are coloured grey indicate areas that were not adequately covered by the scans. Black dots denote reflector positions for referencing laser scans on different days. The solid and dashed lines denote the SMP transects acquired on 4 January 2019 and 11 January 2019, respectively.

After the TLS scans on 4 and 11 January 2019, an undisturbed transect was chosen. Using a measurement tape, SMP measurements were taken at regular (~90 cm) spacing. To exactly locate the SMP transect acquired on 4 January in the scan, the first and last SMP measurements were marked with bamboo poles. These bamboo poles are visible in the scan from 11 January and provide exact coordinates of the first and last SMP measurement of the transect. The second SMP transect was acquired after the last scan, but here, locating the SMP transect inside the TLS field was done by measuring the distance to reference bamboo poles put in before scanning.

The laser scanning data from PEA and HAM are published in Wever and others (2018) and Wever (2022), respectively.

SNOWPACK model

For simulating the drifting snow and accumulation processes, we use the detailed, physics-based, multi-layer snow model SNOWPACK (Lehning and others, 2002a, b), and we compare three setups (see Table 1). Note that these three setups all contain modifications for polar regions and are therefore slightly different from the default settings for alpine snow cover simulations.

Table 1. SNOWPACK parameterizations for the albedo scheme and for parameters that differ between simulation setups

Note that these three setups all contain modifications for polar regions and are therefore slightly different from the default settings for alpine snow cover simulations.

In the first two setups, we use a new snow density scheme developed based on measurements in alpine terrain, but both setups differ in the treatment of wind compaction from drifting snow. In the first setup, which we refer to as ‘default’, the top layers (uppermost 7 cm) can undergo enhanced compaction when the actual wind speed u (m s⁻¹) exceeds a predefined threshold wind speed u ₀ (m s⁻¹) (Brun and others, 1997; Groot Zwaaftink and others, 2013a). Brun and others (1997) use a variable u ₀ based on surface snow properties, whereas Groot Zwaaftink and others (2013a) use a fixed threshold of 5 m s⁻¹. The enhanced compaction is achieved by modifying the default strain rate $\epsilon$, to define an enhanced wind compaction strain rate ($\epsilon _{{\rm enh}}$) as a function of depth below the snow surface d (m) when u > u ₀:

(3)$$\epsilon_{{\rm enh}} = \left(1 + A( d) \left(u - u_{{\rm 0}} \right)^n \right)\epsilon$$

with n = 3 and A(d) taken as:

(4)$$A( d) = 2.7 A_{0} \left(1-{d\over 1.25d_{{\rm max}}} \right)$$

with A ₀ and d _max set to 5 sⁿ m⁻ⁿ and 0.07 m, respectively.

The second setup achieves compaction of light near-surface snow layers by eroding them in the model and redepositing the associated mass with higher density. This scheme termed ‘redeposit’ has been proposed by Keenan and others (2021), who showed an improvement of the firn density in the uppermost 10 m of the firn over the Antarctic continent. It uses the alpine new snow density parameterization (Schmucki and others, 2014) for accumulations from precipitation and the new snow density parameterization developed for Antarctica (Groot Zwaaftink and others, 2013a) for deposited drifting snow mass. The drifting snow mass deposition is equal to the calculated eroded mass, derived from the mass transport in the saltation layer Q (kg m⁻¹ s⁻¹). Saltation initiates once the actual friction velocity exceeds the threshold friction velocity (Schmidt, 1982a), with the latter being calculated by SNOWPACK based on the snow properties of the uppermost layer (see Eqn (1) in Lehning and Fierz, 2008). Microstructural properties considered in the calculation of the threshold friction velocity are coordination number, which is a function of snow density, as well as grain size, bond size, sphericity and dendricity. Once the conditions for saltation are met, the saltation mass transport Q is calculated following Eqn (2) in Lehning and Fierz (2008). Note that Vionnet (2012) discusses that this equation is not implemented correctly in the SNOWPACK model, with the coefficients being expressed with inconsistent units, while also acknowledging that in spite of this error, the agreement with observations is still good (Vionnet, 2012; Vionnet and others, 2014). To estimate the locally eroded mass M _e (kg m⁻²) from the firn layer from the horizontal saltation mass transport Q, the latter is scaled by the fetch length L (m) and model time step Δt (s):

(5)$$M_{{\rm e}} = Q {\Delta t\over L}$$

By default, L is set to 10 m. After calculating M _e, as many layers are eroded as necessary to satisfy M _e. In the sensitivity study, we test for a range of values for L (from 1 to 1 × 10⁵ m) and also test for several erosion limiting thresholds. It is important to recognize that in this approach, we assume that only the mass flux in the saltation layer (and not the suspension layer) contributes to local erosion and deposition. The saltation layer depth is typically ~10 cm (e.g. Doorschot and Lehning, 2002), while the drifting and blowing snow layer frequently extends beyond that. However, the saltation layer carries substantial mass. We consider the fetch length as a pragmatic tuning parameter that also encapsulates the fact that the suspension layer is ignored.

As a third setup, we use the scheme presented by Groot Zwaaftink and others (2013a), which they termed ‘event-driven’ deposition scheme. In that scheme, snow from precipitation is kept in a separate storage and is only added to the firn layer when the 100 h moving average wind speed is higher than 4 and lower than 7 m s⁻¹. The assumption is that for low wind speeds, any precipitation will be eroded during following high wind speed events, whereas for high wind speeds, snow particles are not sticking to the surface. In that study, a new parameterization for new snow density was developed based on field observations at Dome C, specifically describing high wind speed conditions. This parameterization is also used for depositing drifting snow mass in the second setup.

The SNOWPACK simulation time step was set to 15 min. For each site, simulations start at 1 January 1990 without snow and end at the last day of observational data, which is on 1 January 2013 for D47 and in 2018–19 for the other sites. For PEA, 5.7 m of snow builds up in the simulation. For all the other sites, the specified maximum simulated snow depth of 10 m is reached. To reduce computational efforts, SNOWPACK was configured to remove snow layers at the bottom of the column when the snow depth exceeds 10 m. We did not perform spin-ups, since we are only interested in the near surface layers (approximately uppermost 1 m of firn).

For PEA, HAM and D17, we aim to compare properties of newly added layers in the SNOWPACK model to observations of accumulations in the field. In the simulation, each layer gets assigned a deposition time stamp, which allows for calculating the deposited mass within a time period. By default, however, SNOWPACK merges layers for numerical efficiency when layers have similar properties, layers are thinned through sublimation or settling, or volumetric ice content decreases by melting. This makes tracking of deposited layers inaccurate and therefore, we switch off layer merging for the ~2 week study period. For the same reason, we also force SNOWPACK to add layers every time precipitation is reported in the forcing data, and we set the minimum layer depth for deposition to 1 mm.

Atmospheric forcing

We drive the simulations with hourly data of near-surface air temperature, relative humidity, wind speed, surface incoming shortwave and longwave radiation and precipitation, from the gridcell closest to the respective study sites from the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2, Gelaro and others, 2017). Among other atmospheric reanalyses, MERRA-2 has relatively low bias in accumulation rates in Antarctica (Gossart and others, 2019; Medley and Thomas, 2019). Gossart and others (2019) found that the accuracy of MERRA-2 wind speed is comparable to other models. Also, Keenan and others (2021) show that MERRA-2 does not tend to consistently over- or underestimate the highest wind speeds.

Note that these evaluations have been performed using AWSs, which generally are installed in locations that are representative for a larger area, consistent with the relatively coarse spatial resolution of MERRA-2. However, the area near PEA, for example, is characterized by strong topographic variability. Also ice rises exhibit small-scale variability in wind speed and accumulation when comparing windward and lee slopes, which could be captured in models with a higher resolution than MERRA-2 (Lenaerts and others, 2014). To further illustrate this, Figures 3b, d show that for PEA, MERRA-2 wind speed is higher than the observed wind speed and exhibits reduced temporal variability. This notion motivated us to also run simulations where we used in situ observed wind speed instead of MERRA-2 wind speed when available, while using MERRA-2 for the other variables and using MERRA-2 wind speed outside the observational period. Since observational data from Antarctica are relatively sparse, often reanalyses are used as meteorological forcing. By comparing both MERRA-2 with in situ data-driven simulations, we provide a first estimate of the uncertainties introduced by using forcing data from an atmospheric model.

Fig. 3. Meteorological data for the PEA field site from (a, b) 18 December 2016 to 1 January 2017 and (c, d) 1 January to 12 January 2017 for (a, c) bias-corrected cumulative precipitation, observed accumulated mass from density and depth measured using, respectively, a box cutter and ruler (crosses), SMP and ruler (circles) and SMP and TLS (squares), and air temperature from MERRA-2 (solid line) and measurements (dashed line). (b, d) MERRA-2 wind speed (solid) and in situ measured wind speed transformed to 10 m (dashed) and measured and modelled drifting snow amounts, scaled with the maximum in the study period, and only for times when measurements indicated drifting snow, or simulations predicted drifting snow (i.e. zero values not shown).

Even though we discussed earlier that MERRA-2 represents annual SMB adequately, MERRA-2 modelled SMB is the result of several processes, for the Antarctic ice Sheet most importantly precipitation and surface sublimation. This means that a correct SMB can also be achieved when both precipitation and surface sublimation are correspondingly over- or underestimated. Additionally, MERRA-2 does not consider drifting snow, and the associated substantial drifting snow sublimation that can occur (Schmidt, 1982b; Wever and others, 2009). Neglecting drifting snow sublimation is partly compensated for by a dryer atmospheric boundary layer and enhanced surface sublimation (Bintanja, 2001). Furthermore, the effect of drifting snow sublimation relative to annual precipitation is relatively poorly constrained and depends on the climatological conditions. Reported values range from ~2% for the Alps (Groot Zwaaftink and others, 2013b), 12% at Halley Station Antarctica (King and others, 2001), 22–34% in the Arctic (Liston and Sturm, 2002), to almost similar amounts as precipitation for the Arctic Ocean (Yang and others, 2010).

Additionally, local processes, convergence of drifting snow and advection from upwind areas can increase depositions compared to the surrounding area (e.g. Pomeroy and others, 1997; Essery and others, 1999). Lenaerts and van den Broeke (2012a) and Kausch and others (2020) demonstrate these effects for an ice rise east of HAM, where also orographic uplift over the ice rise increases accumulation. This occurs at smaller scales than captured by the MERRA-2 grid. Similarly, it is expected that the heterogeneous surroundings of PEA at scales smaller than the MERRA-2 grid cause local deviations in accumulation. From literature, it is clear that this could be overcome by using downscaling techniques for wind speed (e.g. Reynolds and others, 2021), combined with spatially explicit drifting snow models, such as SnowTran-3D (Liston and others, 2007) or Alpine3D with explicit drift (Lehning and others, 2008). To avoid the complexity and added computational costs of such models, we choose to bias-correct MERRA-2 precipitation when available field data allows, as we will now discuss.

For precipitation, our field data containing both depth and density allow to estimate accumulated mass (see Figs 4a and 3a, c). We found that precipitation amounts from MERRA-2 were significantly lower than the observed accumulated mass at the HAM and PEA field sites for the studied period. When accumulation amounts are over- or underestimated, the settling of the new snow can lead to a bias in density originating from the bias in precipitation amount. Based on the field observations, we increased the precipitation during the studied period by a factor 3 for the HAM field site. Similarly, also for the PEA field site, a factor 3 was required to achieve good correspondence for the precipitation events between 21 December and 24 December, and a factor 9 for the 30 December and 31 December precipitation event. Table S1 in the Supplementary material lists all the data sources and corrections for the meteorological forcing data.

Fig. 4. Meteorological data for the HAM field site from 27 December 2018 to 12 January 2019 for (a) bias-corrected cumulative precipitation, observed accumulated mass from density and depth measured using, respectively, SMP and TLS, and air temperature from MERRA-2. (b) MERRA-2 wind speed and modelled eroded and subsequently redeposited snow.

Figures 3 and 4 show the meteorological conditions during the study period for the PEA (18 December 2016 to 11 January 2017) and HAM (27 December 2018 to 11 January 2019) field sites, respectively. During these fieldwork periods, at both sites a total accumulation of ~40 kg m⁻² was reconstructed from TLS and SMP, as indicated by the blue squared markers in Figures 3a, c and 4a. As described above, the MERRA-2 time series for precipitation shown if Figures 3 and 4 have been adjusted to fit the observed accumulated mass. During both fieldwork periods, the accumulations built up over a few precipitation events.

Results and discussion

Drifting and blowing snow

We first verify the calculation of saltation mass transport by SNOWPACK compared to observed drifting snow mass transport, as this is a critical component of the redeposit scheme. Figure 5 compares the half-hourly sums of observed drifting snow mass transport with the half-hourly sums of simulated saltation mass transport using the redeposit scheme for the three sites. Figure S4 in the Supplementary material shows the same figure for daily sums. The simulated saltation mass transport over the full period is 0.63 and 0.86 times the observed drifting snow mass transport for D17 and D47, respectively. The underestimation of transported mass in the simulations is to be expected since the saltation layer depth is much smaller than the depth covered by the two FlowCapt^TM sensors. The coefficients of determination for all half-hourly time steps are r ² = 0.42 and r ² = 0.50, for D17 and D47, respectively, when using in situ wind speed (Figs 5a, b). These coefficients of determination are significant at the p < 0.001 level. Similar as reported by Amory (2020), the observed frequency of drifting snow is 67 and $82\percnt$ for D17 and D47, respectively. These frequencies are higher than the simulations reproduce (50 and $67\percnt$ for D17 and D47, respectively).

Fig. 5. Scatter plots showing half-hourly sums of observed drifting snow mass transport versus SNOWPACK-simulated saltation mass transport using the redeposit scheme for sites (a, d) D17, (b, e) D47 and (c, f) PEA. In (a, b, c), in situ observed wind speed was used for the simulations and in (d, e, f) MERRA-2 wind speed was used for the simulations. Different colours and markers are used for periods with and without precipitation (defined as <0.001 kg m⁻² h⁻¹) in MERRA-2. The solid black line denotes the 1:1 line and RMSE denotes the root mean square error. Note that zero values are also plotted.

The impact of using MERRA-2 wind or in situ wind to determine saltation mass transport in SNOWPACK is larger for D17 than for D47. For D17, the r ² decreases from 0.42 to 0.29, when using MERRA-2 wind instead of in situ wind, whereas for D47, the r ² is 0.50 irrespective of the wind speed source. The differences in r ² for D17 tested significant on the 95% confidence level using a two-sided test with the statistical tool cocor (Diedenhofen and Musch, 2015).

In Figure 5, we show periods without precipitation in MERRA-2 (defined as <0.001 kg m⁻² h⁻¹) in orange colour and different marker. The results show that the highest saltation mass transport in SNOWPACK and the highest observed drifting snow mass transport are under the presence of precipitation. Furthermore, the SNOWPACK calculated saltation mass transport has lower correspondence with observed drifting snow mass transport for periods without precipitation than for all periods. For D17 and D47, coefficients of determination for periods without precipitation are r ² = 0.22 and r ² = 0.23, respectively, both significant at the p < 0.001 level and statistically different from the value for periods with precipitation, based on cocor.

We interpret this result as follows: with snow particles in the air from precipitation, drifting snow is easier to initiate, since there are already particles in the air to break loose additional snow particles from the surface. Furthermore, fresh accumulations have not yet sintered or settled and can be eroded at even low wind speeds. This apparently poses less of a challenge for the SNOWPACK model to calculate. Without precipitation providing new snow at the surface or snow particles in the air, the friction velocity required for initiating drifting snow depends on the evolution of the snow microstructure at the surface (Guyomarc'h and Mérindol, 1998; Vionnet and others, 2013). As time passes since the last precipitation or drifting snow event, the friction velocity required to initiate the next snow erosion from the surface increases (Vionnet and others, 2013). However, once snow particles are mobilized, it creates a chain reaction where the impact energy from drifting snow particles on the surface breaks cohesive bonds between surface snow particles, causing more particles to eject (e.g. Comola and others, 2019). These processes are inherently more challenging to simulate, since they rely on an accurate simulation of snow microstructure.

For PEA, the half-hourly sums of simulated saltation mass transport and observed mass transport are substantially smaller than for D17 and D47 (Fig. 5c). This is mostly due to the lower wind speeds at this site (Fig. S2c in the Supplementary material). Furthermore, when using MERRA-2 wind speed, the simulated mass fluxes are overestimated compared to the simulations using in situ wind speed (Fig. 5c versus Fig. 5f). Figure S2c in the Supplementary material shows that MERRA-2 wind speed is generally ~5 m s⁻¹ higher than observed. The observed drifting snow mass transport from PEA is substantially smaller than the simulated saltation mass transport, even when using in situ wind speed (Fig. 5c) since, as mentioned before, the SPC was located above the saltation layer where typically the largest mass fluxes occur.

Figure 6 compares the performance of the different model setups in reproducing drifting snow mass fluxes, for sites D17 and D47. The performance between the default and redeposit scheme is almost identical for periods with precipitation in MERRA-2, as confirmed by the statistical significance test for the correlation coefficients. Both schemes treat the density and microstructural properties of precipitation similarly, resulting in similar SNOWPACK calculated drifting snow mass fluxes. However, for cases without precipitation, defined as <0.001 kg m⁻² h⁻¹, the performance between both schemes is different. As will be shown later, the temporal evolution of near surface density varies between both schemes. The increasing discrepancy in surface layer properties over time leads to deviating drifting snow mass fluxes when precipitation, which provides for similar surface layer properties, is absent.

Fig. 6. Scatter plots showing half-hourly sums of observed drifting snow mass transport versus SNOWPACK-simulated saltation mass transport using the (a, d) default scheme, (b, e) redeposit scheme and (c, f) the event-driven scheme for sites (a, b, c) D17, (d, e, f) D47, using in situ observed wind speed. Different colours and markers are used for periods with and without precipitation (defined as <0.001 kg m⁻² h⁻¹) in MERRA-2. The solid black line denotes the 1:1 line and RMSE denotes the root mean square error. Note that zero values are also plotted.

The event-driven scheme shows lower agreement between observed and SNOWPACK calculated drifting snow mass fluxes, indicated by a lower correlation coefficients and higher RMSE (see Figs 6c, f). Note that these fluxes only concern the calculation based on snow properties of the surface layer, not the amount of precipitation kept in the separate storage when wind speeds are either below or above the prescribed wind speed range for deposition. We attribute these differences to the fact that low-density snow is never deposited in this scheme, and that there is a temporal discrepancy between actual precipitation simulated by MERRA-2, and the deposition of that mass on the surface in the event-driven scheme.

Observed accumulations at HAM

For the field site at HAM, snow accumulation and erosion patterns determined from TLS are shown in Figure 2. We find that these patterns are spatially variable, and that in some cases both erosion and accumulation features can be identified with length scales of 5–10 m (see e.g. Fig. 2d). Similar erosion/deposition patterns, albeit with a smaller typical length scale, have been shown to exist on the Antarctic plateau (Picard and others, 2019). The first precipitation event (accumulation pattern shown in Fig. 2a) started by the end of the day on 29 December, ending on 1 January, accompanied by relatively high near-surface wind speeds, exceeding 10 m s⁻¹ (see Fig. 4b). The accumulation pattern from the second accumulation event, which occurred during calm conditions on 3 January, is shown in Figure 2b. As described in Section S2 in the Supplementary material, the scan of 2 January is less accurate, giving us lower confidence in the existence of the net erosion area at the bottom left of Figure 2a. Figure 2c, which compares scans with adequate accuracy, shows that the first two precipitation and drifting snow episodes resulted in mostly accumulation in the scanned area. In contrast, the third accumulation event that started on 9 January, with wind speeds again exceeding 10 m s⁻¹, created both net accumulation areas as well as net erosion areas (Fig. 2d). The presence of erodible snow near the surface created similar conditions as reported for PEA (Sommer and others, 2018). The low-density snow was likely mostly removed and redeposited with higher density in barchan snow dunes.

The density profiles from two SMP transects inside the TLS field are shown in Figures 7 and 8. The exact locations of the transects are indicated by the solid and dashed lines, respectively, in Figure 2. The objective was to cover ~50 m distance, which should capture typical dimensions of surface features from metre scale (Picard and others, 2019) to 50 m scale (Sommer and others, 2018).

Fig. 7. Density from SMP along a transect at the HAM field site (solid line in Fig. 2), acquired on 4 January 2019. The snow surfaces of 2 January and 27 December are shown with a dashed and solid line, respectively. Note that after 32 m into the transect, it was decided to only measure the full depth of the SMP every second measurement, and restrict to ~40 cm depth otherwise, due to time constraints.

Fig. 8. Density from SMP along a transect at the HAM field site (dashed line in Fig. 2), acquired on 11 January 2019. The snow surfaces of 4 January and 27 December are shown with a dashed and solid line, respectively.

The first transect, measured on 4 January 2019, is covered by the TLS scans from 27 December and 2 January. The surfaces from these scans are also plotted in Figure 7. We find that the accumulation depth is spatially variable along the transect, ranging from ~15 cm to only a few centimetres. Furthermore, the accumulation occurring between 27 December and 2 January has a density of ~300–350 kg m⁻³, whereas the accumulation between 2 January and 4 January has much lower density (<200 kg m⁻³).

On 11 January, a similar transect was measured (see Fig. 8). The spatial sampling is lower than on 4 January, because of logistical time constraints. The position of the surface on 4 January and 27 December is depicted in that figure as well. Between the surface of 27 December and 4 January, we do not find the clear distinction anymore as in the transect acquired on 4 January. We argue that this illustrates that strong wind eroded the existing, low-density, snow layer and redeposited high-density snow in the form of dunes.

In the field, it was also observed that high air temperatures and sunny weather created discontinuous melt/freeze crusts during the days before 27 December (see Fig. 4a). Both transects show these layers as high-density layers near the surface from 27 December (solid line in Fig. 7a).

Observed and modelled accumulated snow density at HAM

The observations from the HAM field site of the density of accumulations during three periods during the 2018–19 field season: (1) 27 December to 2 January, (2) 2 January to 4 January and (3) 4 January to 11 January are shown in Figure 9. The violin plots show the distribution of all the SMP measurement points in between the two surfaces determined from TLS. For periods (1) and (2), data from the first transect, which were acquired after period (2) were used (Fig. 7), whereas for period (3), data from the second transect, which were acquired after period (3), were used (Fig. 8).

Fig. 9. Distributions of SMP surveyed new snow density (blue) and accumulation depth (red) of accumulations and simulation results for three simulation setups, for three time periods for the HAM field site. Distributions are shown as violin plots (Hintze and Nelson, 1998). The violin plot combines a box plot (shown in black, indicating the median by a white dot, the inter-quartile range by a black box and either the minimum or maximum value, or 1.5 times the inter-quartile range, whichever is closer to the median, by the black lines) with a symmetrically plotted rotated kernel density showing the full, smoothed, distribution. Note that the event-driven simulation does not reproduce an accumulation event between 2 January and 4 January.

The first episode exhibits the highest new snow density of all three periods, with a median value (± median absolute deviation (MAD)) of 319 ± 38 kg m⁻³. The accumulations during the following two episodes exhibited much lower density, with a median (±MAD) of 173 ± 59 and 221 ± 82 kg m⁻³, respectively. This corresponds well to the pattern of wind speed during this period, with the highest wind speed occurring in the period 27 December to 2 January. During the second period, wind speeds were very low, and drifting snow occurrence was correspondingly low, as observed in the field and suggested by the simulations shown in Figure 4b. The accumulated snow density from SNOWPACK simulations, also shown in Figure 9, follows a similar pattern over the three accumulation periods in the default and redeposit scheme. The event-driven mode does not reproduce the accumulation between 2 and 4 January, since wind speed was below the threshold for deposition. The redeposit scheme simulation exhibits densities of 274, 117 and 260 kg m⁻³ for the three deposition events, respectively. This is up to 56 kg m⁻³ less than observed. The low density snow observed during the fieldwork likely occurs infrequently, or is only short-lived (i.e. subsequent high wind speed events erode such layers). As shown in Figure S5, the density of the uppermost 5 cm of the firn is lower than 280 kg m⁻³ <15% of the time in the simulations using the redeposit scheme for HAM.

The red violin plots on the right-hand side in Figure 9 show the accumulation depths as determined from TLS and the simulated accumulation depth. Since the precipitation time series was bias-corrected to match observed accumulated mass, differences in simulated accumulation depths arise from differences in simulated density.

Both measured density and accumulation depth show large variability. For accumulation depth, this corresponds to the irregular snow depositions and dune formation at the surface. For density, we attribute the spread to inherent natural variability when accumulation patterns are highly variable (Fig. 2). However, inaccuracies in the registration of the TLS scans and location of SMP measurements inside the TLS field may also occasionally cause the old snow surface to be considered part of the fresh deposition when analysing the SMP measurements.

Observed and modelled accumulated snow density at PEA

For PEA, two marked precipitation events can be identified during the study period (Fig. 3a). The first accumulation period between 21 December and 23 December was accompanied by low wind speeds (generally <6 m s⁻¹), and air temperatures ~−5°C (Fig. 3a). This resulted in low fresh snow densities of <100 kg m⁻³. This is confirmed by both SMP as well as manual measurements using a density cutter and using the hard, old snow surface as reference (see Fig. 10). The median (±MAD) SMP measured density on 22, 23 and 27 December is 68 ± 17, 95 ± 14 and 112 ± 18 kg m⁻³, respectively. This general increasing trend in density during these 3 days, as also shown in Figure 10, is for a large part caused by settling of the fresh snow. Relatively high air temperatures during the period suggest that settling rates were likely significant (Fig. 3a).

Fig. 10. Distributions of SMP and density-cutter surveyed new snow density (blue and cyan, respectively) and accumulation depth (red and orange, respectively) of accumulations and simulation results for three simulation setups, for five time periods for the PEA field site. Distributions are shown as violin plots (Hintze and Nelson, 1998). See caption to Figure 9 for details.

On 28 December, some blowing snow was recorded (Fig. 3b), which promptly increased the density at the surface to 210 ± 61 kg m⁻³ in the SMP measurements. The second large accumulation event during the studied period at PEA on 30 and 31 December occurred during high wind speeds (exceeding 10 m s⁻¹), and significant amounts of drifting snow were observed (Fig. 3b). As already discussed by Sommer and others (2018), the high wind speed eroded virtually all of the low-density snow from the previous precipitation event, and combined with new mass input from precipitation, the resulting drifting snow formed high-density snow dunes in the shape of Barchan dunes. The SMP measurements show a density of 419 ± 18 kg m⁻³ after this event.

The temporal evolution is well reproduced by the simulation setups allowing for low-density snow, using alpine-based new snow density parameterizations. We find a similar increase in simulated density in both schemes using either an enhanced wind compaction of surface layers (default setup) or redeposition of eroded layers (redeposit scheme). For example, the simulated density on 22 December with the redeposit scheme is 91 kg m⁻³ (default setup: 92 kg m⁻³), which is higher than the observed density. On 23 and 27 December, the simulated density by the redeposit scheme increases to 107 and 136 kg m⁻³, respectively. This is a consequence of snow settling, and some minimal drifting snow calculated by the model on 24 December (see Fig. 3b). In the default setup, the density increases over a similar range, to 107 and 154 kg m⁻³, respectively. At the end of December, with increasing wind speed and consequently drifting snow, the simulated density with the redeposit scheme increases to 242 kg m⁻³ on 29 December and further to 357 kg m⁻³ on 31 December. This temporal trend in simulated density corresponds well with the observed trend, and the simulated density on 31 December agrees well with the observed density. In the default setup, the density increases to 211 kg m⁻³ on 29 December and to 308 kg m⁻³ on 31 December, which constitutes a smaller increase than in the redeposit scheme and results in a more than 100 kg m⁻³ lower density than observed. As shown in Figure S5, the density of the uppermost 5 cm of the firn exceeds 200 kg m⁻³ 13% of the time in the simulations using the redeposit scheme for PEA.

The event-driven scheme cannot reproduce the low-density snow. The high wind speeds between 18 and 21 December, just before the accumulation event started (see Fig. 3b), caused deposition with a density of ~230 kg m⁻³. After that until 27 December, no new deposits formed in the event-driven simulation. Once the drifting snow occurred ~28–31 December, accompanied by precipitation, the event-driven scheme reproduces densities for the new accumulations similar to the densities in the other simulations.

Similar as found for the HAM field site, the PEA accumulation depth data shown by the red violin plots on the right-hand side in Figure 10 exhibit spatial variability in accumulation depth, most pronounced for the period following drifting and blowing snow. The light snowfall in calm weather between 22 and 27 December created a new snow layer with more uniform depth. Also, the variability in observed density is low for this period.

Observed and modelled accumulated snow density at D17

The Supplement to Amory and others (2021) describes snow density measurements at D17 during snowfall and drifting snow conditions in January 2014. The presence of a melt-freeze crust from previous warm and rainy weather provided a clear delineation between the old snow surface and the new snow, similar to the older surfaces determined from repeat TLS presented for HAM and PEA. Figure 11 shows the observed density increasing steadily from ~220 kg m⁻³ on 28 January, 00.00 UTC to ~350 kg m⁻³ 18 h later during and due to drifting snow conditions.

Fig. 11. Observed (black dots) and simulated snow density using the default scheme and in situ wind data (cyan), the redeposit scheme and in situ wind data (blue) and redeposit scheme and MERRA-2 wind data (red) of fresh snow accumulations at D17 in January 2017. Note that the event-driven scheme does not produce accumulation during the observational time frame.

We analysed the SNOWPACK simulations for D17 of all layers deposited in the model domain after 27 January, 7.00 UTC. Figure 11 shows that the density in the default and redeposit scheme in SNOWPACK increases at a similar rate as observed when using in situ wind speed, albeit ~12 h before the observed density increase. The total increase in density with the default scheme is higher than with the redeposit scheme. Note that the event-driven scheme does not produce accumulation during the observational time frame, since the wind speed conditions for deposition, as set for this scheme, were not met. It is also clear that using MERRA-2 winds provides too high densities because these wind speeds are higher than locally observed during this event. This event illustrates the inherent complexities of reproducing individual deposition events in wind-dominated cold environments. For the windy environment at D17, the low-density accumulations are a rare occurrence, and only very short-lived (as in this example). As shown in Figure S5, the density of the uppermost 5 cm of the firn is lower than 300 kg m⁻³ <2% of the time in the simulations using the redeposit scheme for D17.

Sensitivity study

Several poorly constrained parameters control simulated erosion in the current version of SNOWPACK. To test the sensitivity of the model results to these variables, we performed a sensitivity study with the data from the D17 site, which provides the longest available time series of drifting snow. We run with MERRA-2 forcing only and test for the sensitivity to roughness length, fetch lengthand an erosion limit, which inhibits erosion of layers based on its properties (density, and threshold friction velocity).

For drifting snow, roughness length is a critical parameter (e.g. Amory and others, 2015, 2017). A wide range has been reported in literature, depending on whether the snow surface is flat (roughness lengths of ~10⁻⁵–10⁻⁴ m), or is shaped by dunes and sastrugi (roughness lengths of ~10⁻³ m). It may vary spatially and seasonally based on the formation of snow erosion and deposition structures (Amory and others, 2017). In turn, the roughness length impacts calculated friction velocities, thereby controlling both the initiation and the amount of erosion. The roughness length parameterization as a function of air temperature, as proposed in Eqn (7) in Amory and others (2017), aims to describe seasonal variations in surface roughness. Figures 12a, d show the sensitivity of simulated saltation mass transport as a function of roughness length for the total amount and r ², respectively. The parameterized roughness length is denoted by A. We find that generally the modelled saltation mass transport increases with increasing roughness length due to the increased friction velocity (Lehning and others, 2000, 2008). The simulated frequency of saltation increases from 42% for a roughness length of 1 × 10⁻⁵ m to 70% for a roughness length of 2.5 × 10⁻² m, after which a decrease is found (not shown). The r ² value, however, is insensitive for roughness lengths smaller than 2 × 10⁻² m. Above that value, the correlation between simulated saltation mass transport and observed drifting snow mass transport drops markedly. The parameterized roughness length performs reasonably well, but does not provide optimal performance. The surface density is fairly constant for low roughness lengths, but decreases for higher values (Fig. 12g). Since MERRA-2 wind is unaltered for the different simulations, increasing roughness length decreases the 3 m wind speed used in the new snow density parameterization.

Fig. 12. Sensitivity of the SNOWPACK model results to the choice of roughness length (a, d, g), fetch length (b, e, h) and a limit on erosion (c, f, i) for the sum of saltation mass transport (a, b, c), the correlation between simulated saltation mass transport and observed drifting snow mass transport (d, e, f) and simulated uppermost 0.1 m snow density (g, h, i). For roughness length, A denotes the parameterization as a function of temperature as proposed by Amory and others (2017). Erosion limit denotes a density threshold above which snow layers cannot erode, except for none, which denotes no limit, and u*, which denotes a limit based on threshold friction velocity. In (a–f), blue solid (orange open) markers indicate results for periods with (without) precipitation in MERRA-2. The default simulation setup used in this study is indicated by a squared marker.

Next, we evaluate the sensitivity to the choice of fetch length. Figures 12b, e show that fetch length has relatively little impact on the total amounts of simulated saltation (Fig. 12b) and the correlation between simulated saltation mass transport and observed drifting snow mass transport (Fig. 12e), particularly during precipitation. The calculation of saltation mass transport is independent of fetch length because the latter is only applied to convert horizontal mass transport into a vertical erosion flux from the firn layer. In contrast to the insensitivity on drifting snow, we find a strong decrease in surface density with increasing fetch length (Fig. 12h), because with increasing fetch length, less of the snow surface is eroded and redeposited to satisfy the saltation mass transport (see Eqn (5)), reducing the impact on density of the uppermost 0.1 m of the firn layer. We attribute the slight increase in r ² for days without precipitation with increasing fetch length (Fig. 12e) to the lower surface density, which provides for erodible snow for longer periods of time.

Finally, we test the use of an erosion threshold to determine if erosion of firn layers continues. The saltation mass transport is determined based on the properties of the uppermost firn layer in the model. The model then erodes as many layers as needed to satisfy the saltation mass. However, it can be argued that when deeper layers exhibit properties that hinder erosion, such as high density, or high bond strength, erosion should be halted. Amory and others (2021) apply a threshold of 450 kg m⁻³. When the snow layer to be eroded exceeds this density, erosion is halted. We test for a range of densities, and also test for halting erosion when the threshold friction velocity for a layer exceeds the actual friction velocity. Figures 12c, f show that the amount of mass in saltation and the correlation with observed drifting snow mass transport is a strong function of density threshold. With the default absence of an erosion threshold (denoted by none) the model provides the highest correlation with observations. Similarly high agreement is achieved with a high surface density threshold, simply because such high surface density occurs only infrequently and the threshold criterion would mostly remain unmet. Applying a criterion where the friction threshold velocity is checked for each layer before erosion also yields lower correlation with observations. It has been argued that the threshold friction velocity to initiate erosion is higher than the friction velocity needed to sustain erosion, since once erosion has been initiated, the momentum from drifting snow particles hitting the surface adds to the friction provided by the flow of air over the surface (e.g. Schmidt, 1980, 1982a; Trouvilliez and others, 2014; Comola and others, 2019). Similar conditions occur during precipitation. The precipitating particles in the air would also cause the snow surface to be eroded at lower wind speeds than without precipitation. While this is not explicitly taken into consideration in SNOWPACK, the model logic compensates for that: first precipitation is added to the firn layer with a relatively low density from the new snow density scheme. Only upon erosion and redeposition, the parameterization for drifting snow is used which generally leads to higher density. This allows for saltation to initiate more easily in SNOWPACK when precipitation is present. Also note that the surface density increases with increasing density to halt erosion, which logically follows since snow below the threshold is susceptible to be removed from the firn layer and redeposited with higher density.

Conclusions

We exploited detailed field observations of drifting snow mass fluxes and snow accumulations in terms of depth and density from drifting snow-dominated environments in East Antarctica for a comparison with model simulations using the 1-D, detailed, physics-based snow cover model SNOWPACK. Statistically significant correlation was found between previously published measured drifting snow fluxes from Adélie Land (see Amory and others, 2020) and SNOWPACK-simulated saltation mass transport. Although the contribution of divergence/convergence of drifting snow and transport in suspension to the total transport in the drifting snow layer cannot be accounted for in a 1-D modelling approach, this suggests that SNOWPACK can provide realistic estimates of the occurrence of saltation, particularly when the local wind climate is represented well in the meteorological forcing data. The saltation mass transport calculated by SNOWPACK was found to be 0.63–0.86 times the observed drifting snow mass transport. Episodes with drifting snow during precipitation are better captured by the simulations than drifting snow events without precipitation.

We determined the spatial variability of accumulation depth and density by combining snow accumulation from repeat TLS and snow micro-penetrometry, for two field sites in eastern Dronning Maud Land at the PEA station and the HAM ice rise. We also included observations from one of the field sites in Adélie Land (D17). We showed that there is a large temporal variability of fresh snow density. At all three field sites, precipitation events with low wind speeds were observed, which resulted in low snow density. These conditions can exist for hours and even days, before high wind events erode the low-density snow and high-density snow is deposited.

Simulations of the snow density of fresh accumulations using SNOWPACK showed that certain model setups are capable of describing the full range of near-surface density. We found that using a stronger compaction term for near-surface snow layers during high wind conditions yielded similar results as when using explicit snow erosion and redeposition of snow layers in the model. The latter setup has the advantage that it is a closer description of the actual processes that allow for surface compaction during drifting snow conditions (i.e. the presence of drifting snow). These results substantiate the earlier published improvement in simulated density profiles found for the uppermost 10 m of snow for a wide range of climatological conditions on the Antarctic ice sheet when using the redeposit scheme (Keenan and others, 2021).

However, it is also important to note that our validation of snow density is performed for the Austral summer period only, and for a very narrow range of climatological conditions, with all sites being in the coastal areas. Also, errors in simulated near-surface density during the case studies were found to be up to 62 kg m⁻³ for the redeposit scheme, and more than 100 kg m⁻³ for the other schemes, which is larger than the standard error in measured density (37 kg m⁻³ for the SMP). This error illustrates the general shortcomings in our understanding of the exact deposition mechanisms in windy environments. Furthermore, uncertainties caused by meteorological forcing (particularly wind and precipitation) for the SNOWPACK model can be substantial, given the biases we reported in the Supplementary material when in situ measurements were available. Ultimately, drifting and blowing snow is a multi-dimensional problem, with advection possible over long distances. Also the sublimation of airborne snow particles can lead to substantial mass loss, a process we did not consider in this study. A 1-D approach, as we applied here, is strongly limited in processes captured.

Our results indicate that the assumption of high near-surface density in windy environments can be violated on short timescales. For D17, which is located in the most windy environment from the sites shown in this study, modelled surface density in the uppermost 5 cm of the firn layer is <150 kg m⁻³ $0.36\percnt$ of the time. For HAM, with lower average wind speeds, this is $1.7\percnt$ of the time. This has implications for using repeat satellite altimetry, which provides a snapshot of accumulation in time, since the spatial and temporal variability in the density of the near surface layers introduces important uncertainties in translating depth changes to the local surface mass balance of the Antarctic ice sheet.