Performance of MAR (v3.11) in simulating the drifting-snow climate and surface mass balance of Adélie Land, East Antarctica

Abstract

Drifting snow, or the wind-driven transport of snow particles originating from clouds and the surface below and above 2 m above ground and their concurrent sublimation, is a poorly documented process on the Antarctic ice sheet, which is inherently lacking in most climate models. Since drifting snow mostly results from erosion of surface particles, a comprehensive evaluation of this process in climate models requires a concurrent assessment of simulated drifting-snow transport and the surface mass balance (SMB). In this paper a new version of the drifting-snow scheme currently embedded in the regional climate model MAR (v3.11) is extensively described. Several important modifications relative to previous version have been implemented and include notably a parameterization for drifting-snow compaction of the uppermost snowpack layer, differentiated snow density at deposition between precipitation and drifting snow, and a rewrite of the threshold friction velocity above which snow erosion initiates. Model results at high resolution (10 km) over Adélie Land, East Antarctica, for the period 2004–2018 are presented and evaluated against available near-surface meteorological observations at half-hourly resolution and annual SMB estimates. The evaluation demonstrates that MAR resolves the local drifting-snow frequency and transport up to the scale of the drifting-snow event and captures the resulting observed climate and SMB variability, suggesting that this model version can be used for continent-wide applications.

1 Introduction

A significant portion of the surface area of Antarctica is affected by wind-driven ablation or accumulation. The net snow accumulation at the ice sheet surface, i.e. the surface mass balance (SMB), is the resultant of mass gains (precipitation and riming), mass losses (water runoff and surface sublimation), and wind-driven snow transport which can either result in mass gain (deposition) or loss (erosion). The snow mass sublimated during transport by the wind is lost by the ice sheet surface to the atmosphere when particles originate from the surface. Although wind-driven snow sublimation has most often been described as an independent term in the SMB equation e.g., its contribution to surface mass loss is inherently included in the spatially integrated erosion–deposition balance. Drifting and blowing snow are usually conventionally distinguished as the wind-driven transport of snow particles respectively below and above a height of 2 m above ground. In this study both processes are combined into the single denomination of drifting snow for convenience, and erosion, deposition, horizontal and vertical transport of wind-driven cloud (i.e. that have not yet reach the surface) and eroded (raised from the surface) snow particles and their concurrent sublimation are all referred to as drifting-snow processes.

The net erosion–deposition balance in areas subject to drifting snow is mainly governed by the interactions between the complex ice surface topography, near-surface flow and surface snowpack state. Katabatic winds flowing over the surface of the Antarctic ice sheet accelerate down steep surface slopes, causing erosion when the wind shear stress is high enough to dislodge particles from the surface. Wind redistribution of snow occurs more generally on every spatial scale at which topographic features generate acceleration or deceleration of the near-surface flow, enhancing or reducing local SMB gradients up to sub-kilometre scales and occasionally forming extensive areas of near-zero to negative SMB in windy and dry interior regions of the ice sheet . From a modelling perspective, this means that drifting-snow processes vary as a function of the horizontal resolution . This is of particular importance for resolving the spatial variability in SMB at the ice sheet margins where drifting-snow processes can be major components of the local SMB , resulting in the transport and sublimation of large volumes of snow and/or their export off the continent boundaries .

As drifting snow triggers interactions between the atmosphere and the ice sheet surface, a comprehensive evaluation of drifting-snow processes using snow-transport models requires consistency between model results and observations for both mass transport and net accumulation rates simultaneously. The much larger availability of SMB observations compared to drifting-snow measurements in Antarctica and the wider applications of modelled SMB products have led model development and evaluation exercises to focus primarily on the representation of the SMB, with a secondary or most often non-existent consideration for drifting-snow processes . SMB and drifting-snow transport, however, are not independent variables. Arbitrary adjustments of model parameters favouring one can be made at the expense of the other e.g.. Similarly, neglecting or underestimating drifting-snow processes induces a smoothing of the modelled SMB gradients across areas of complex topography and can ultimately lead to overestimation of the snow mass input in regionally integrated SMB calculations . It is worth mentioning that the need to explicitly describe drifting-snow processes also increases with the general tendency towards higher horizontal resolution of the atmospheric models used to study the Antarctic SMB.

Measurements of drifting-snow mass fluxes are particularly interesting for evaluating snow-transport models since they constitute the integrated result of all the feedback and dynamical mechanisms (i.e. precipitation, local erosion, horizontal advection from upwind areas and sublimation) that contribute to the presence, amount and residence time of snow particles in the air. The general scarcity in drifting-snow measurements in polar regions is however a constraint to the development of parameterization schemes for large-scale applications and currently hinders quantitative evaluations of contrasting, continent-wide model estimates of drifting-snow mass and sublimation fluxes in Antarctica see. Numerical challenges associated with modelling drifting snow at the regional scale also arise from the numerous interactions of drifting-snow particles with the atmosphere and the snow surface organized in a complex system of positive and negative feedback mechanisms. The difficulty involved in capturing the resulting strong non-linearity of drifting-snow processes depends on the representation and number of feedbacks accounted for and is mirrored through a high sensitivity of model results to parameter choices and significant discrepancies between simulated and observed snow mass fluxes .

The polar-oriented regional climate model MAR includes a drifting-snow scheme initially developed to improve the representation of the Antarctic SMB . However, the drifting-snow scheme has only been used so far to study separately wind-driven ablation and individual drifting-snow events in coastal East Antarctica. Unpublished preliminary experiments with former physical parameterizations of drifting snow in MAR did not lead to a continent-wide agreement between model simulations and both drifting-snow and SMB observations. As a result, drifting snow has been kept disabled in recent decade-long investigations of the SMB of the Greenland e.g. and Antarctic e.g. ice sheets with MAR, despite the potentially missing aspects related to the important feedback mechanisms induced by drifting snow.

In this paper a modified version of the original drifting-snow scheme implemented in MAR is assessed through a concurrent evaluation of the drifting-snow climate and SMB reproduced by the model against a multi-year database of drifting-snow mass fluxes and SMB estimates. The evaluation focuses on the marginal slopes of Adélie Land, a katabatic wind region of East Antarctica which experiences drifting snow frequently and where the SMB exhibits a high spatial variability related to drifting-snow processes . The coupled atmospheric and snowpack components of MAR are presented in Sect. 2. The drifting-snow scheme is described in Sect. 3 together with the new developments and main changes relative to the original version. Section 4 provides information on the study area, the available data and the evaluation strategy. The modelled near-surface climate, drifting-snow frequency, mass transport and SMB are evaluated in Sect. 5. The sensitivity to the model version and input parameters of the drifting-snow scheme are discussed in Sect. 6 before concluding the paper.

2 Model descriptions

2.1 Atmospheric model

MAR is a hydrostatic atmospheric model originally developed to simulate the climate over high-latitude regions. The atmospheric dynamics in MAR are described in . Cloud microphysical processes and resulting precipitation are simulated by solving conservation equations for specific humidity, cloud droplets and ice crystals, raindrops, and snow particles . The radiative transfer through the atmosphere is adapted from , and cloud radiative properties are computed from the concentration of cloud droplets and cloud ice crystals. Turbulence is resolved in the surface layer following the Monin–Obukhov similarity theory and above the surface layer using a local closure scheme adapted to the stable boundary layer .

In this study MAR version 3.11 in its Antarctic setup is used with the updates described in , simply referred to as MAR hereafter. The simulations are performed on a grid of 80 $\times$ 80 cells at 10 km horizontal resolution to reduce computational cost and facilitate development and sensitivity experiments. The time step is set to 60 s, for a computational cost of 72 CPU hours per year of simulation in the chosen configuration. The topography is obtained through aggregation of the 1 km Bedmap2 surface elevation dataset . The model is driven at its lateral boundaries (pressure, wind speed, temperature, specific humidity), at the top of the troposphere (temperature, wind speed) and at the ocean surface (sea ice concentration, sea surface temperature) by 6-hourly ERA5 reanalysis fields . The atmosphere is described on a stretched grid with 24 vertical terrain-following levels, of which 8 and 5 are respectively located in the lowest 100 and 20 m of the atmosphere with a lowest level at 2 m. The model is relaxed towards the forcing solutions of wind speed and temperature from the top of the troposphere (i.e. above 10 km) to the uppermost atmospheric level following .

2.2 Snowpack model

The atmospheric part of MAR is coupled with the one-dimensional multi-layer surface model SISVAT (Soil Ice Snow Vegetation Atmosphere Transfer), which handles energy and mass transfer between the surface and the atmospheric boundary layer . SISVAT includes a representation of snow and ice layers and comprises subroutines for snow metamorphism, surface albedo, meltwater percolation, retention and refreezing. In the present study, 30 snow/ice layers are used to describe the snowpack with a fixed total resolved snowpack of 20 m in thickness. An aggregation scheme automatically manages the stratification of the snowpack due to precipitation, erosion–deposition of snow, mechanical compaction, and thermal and melting/refreezing metamorphism, enabling a dynamical evolution of the physical characteristics (temperature, density, water content, grain shape and size) of the different layers over time. If precipitation or deposition occurs when the snowpack already comprises the maximum number of layers, the formation of a new layer at the surface is achieved through aggregation of internal sub-surface layers. More generally, aggregation of adjacent layers is permitted according to the similarity of their physical properties. Thick layers can also be split to refine the discretization of the snowpack when the number of layers is lower than 10. Maximum layer thicknesses of the four uppermost layers of the snowpack are also prescribed (0.02, 0.05, 0.1 and 0.3 m) to ensure a fine discretization adapted to the description of sub-surface processes such as heat exchange with the surface and diffusion within the snowpack. Mass and heat are conserved along the snowpack stratification procedure. The snowpack was uniformly initialized with snow grain shape parameters of fresh snow (see Sect. 3.3 for definition) and a density of 500 kg m${}^{-\mathrm{3}}$ assuming a null liquid water content. The initial surface snowpack temperature is set to the reanalysis near-surface air temperature and then discretized along a predefined layer thickness profile as a function of distance to the surface to determine the temperature of internal snowpack layers. The model was then run from 1994 so that the snowpack had reached equilibrium with the climate preceding the period of interest (2004–2018) after a spin-up time of 10 years.

3 Drifting-snow scheme

This section describes the drifting-snow physics currently implemented in MAR. Details on the computation of the threshold friction velocity for snow erosion, snow-transport modes, interactions of drifting snow with the atmosphere and the surface, and then snow erosion and surface roughness are successively provided in the following subsections. A schematic diagram (Fig. ) provides a general overview of the drifting-snow scheme.

Figure 1

Schematic illustration of the drifting-snow scheme in MAR. Model variables are marked in bold black. The blue arrows denote mass and energy exchanges and drifting-snow processes are indicated in blue. The different computation steps listed in Sect. 3.5 are reported in red.

[Figure omitted. See PDF]

3.1 Threshold friction velocity for initiation of snow erosion

Erosion of snow is usually considered to initiate when the shear stress exerted by the flow at the surface (determined by the friction velocity $u_{*}$ in m s${}^{-\mathrm{1}}$) exceeds the threshold value for aerodynamic entrainment, i.e. the threshold friction velocity $u_{*\mathrm{t}}$ (in m s${}^{-\mathrm{1}}$) determined by the resistive gravitational and cohesive forces. Resistive forces depend on temperature and metamorphism history of the snow surface and involve various snow particle characteristics such as inter-particle cohesion, density, grain shape and size. It follows that an accurate prognostic of $u_{*\mathrm{t}}$ requires a detailed representation of these characteristics, which are particularly undocumented in Antarctica and for which measurements are generally very limited in the literature. As an alternative, density has often been proposed as a governing factor in parameterizations of $u_{*\mathrm{t}}$ e.g.. The same approach is followed here; erosion in the model occurs when $u_{*}>u_{*\mathrm{t}}$, where $u_{*\mathrm{t}}$ is imposed by the uppermost snow layer density:

1$$u_{*\mathrm{t}}=u_{*{\mathrm{t}}_{\mathrm{0}}}\exp \left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{0}}}}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{s}}}}\right),$$ in which ${\mathit{\rho }}_{\mathrm{s}}$ is the surface snow density (kg m${}^{-\mathrm{3}}$), ${\mathit{\rho }}_{\mathrm{i}}$ is the density of ice (920 kg m${}^{-\mathrm{3}}$) and ${\mathit{\rho }}_{\mathrm{0}}$ is the density of fresh snow (set to 300 kg m${}^{-\mathrm{3}}$). The expression for $u_{*{\mathrm{t}}_{\mathrm{0}}}$ is retained from : $$\begin{array}{}\text{2}& {\displaystyle }{u}_{*{\mathrm{t}}_{\mathrm{0}}}={\displaystyle \frac{\log \left(\mathrm{2.868}\right)-\log (\mathrm{1}+i_{\mathrm{ER}})}{\mathrm{0.085}}}{C}_{\mathrm{D}}^{\mathrm{0.5}},\text{3}& {\displaystyle }{i}_{\mathrm{ER}}=\mathrm{0.75}d-\mathrm{0.5}s+\mathrm{0.5},\text{4}& {\displaystyle }{C}_{\mathrm{D}}={\displaystyle \frac{u_{*}^{\mathrm{2}}}{U^{\mathrm{2}}}},\end{array}$$ where $i_{\mathrm{ER}}$ is an erodibility index describing the potential for snow erosion, $d$ (dendricity) and $s$ (sphericity) are the snow grain shape parameters, $C_{\mathrm{D}}$ is the drag coefficient for momentum, $U$ is the wind speed at the lowest prognostic level of the model (m s${}^{-\mathrm{1}}$), and $u_{*}$ is classically obtained through integration of stability correction functions for momentum over the atmospheric boundary layer. Dendricity represents the remaining initial geometry of fresh snow particles and varies from 0 to 1 with high values of $d$ describing fresh snow layers. Sphericity varies equally from 0 to 1 and defines the ratio of rounded to angular shapes in the snow layer. In previous studies involving drifting-snow applications with MAR, $i_{\mathrm{ER}}$ was defined as a function of surface snow characteristics ($d$, $s$ and grain size). To reduce the number of sensitivity parameters, as in , $u_{*\mathrm{t}}$ is assumed to be independent of particle size, and constant snow grain shape parameters are assigned ($d=s=\mathrm{0.5}$), implying an erodibility index of 0.625 and a minimum $u_{*\mathrm{t}}$ value of 0.3 m s${}^{-\mathrm{1}}$ for ${\mathit{\rho }}_{\mathrm{s}}$ of 300 kg m${}^{-\mathrm{3}}$. Additional criteria for snow erosion require that ${\mathit{\rho }}_{\mathrm{s}}$ does not exceed ${\mathit{\rho }}_{max}=\mathrm{450}$ kg m${}^{-\mathrm{3}}$ and that no wet layer has formed at the top of the snowpack. The sensitivity to the formulation of $u_{*\mathrm{t}}$ is quantified and discussed in Sect. 6.2.

3.2 Snow-transport modes

Particle motions in drifting snow are generally described through the two main transport modes consisting in saltation and turbulent suspension. Once the resistive forces have been overcome by the atmospheric drag force, erosion initiates through the saltation process, in which particles become mobile and periodically bounce on the surface within heights of the order of 10 cm. Turbulent suspension of snow occurs when snow particles obtain sufficient upward momentum to be entrained in the atmosphere by turbulent eddies from the top of the saltation layer without contact with the surface.

The drifting-snow scheme of MAR uses a set of semi-empirical formulations to predict the contribution of erosion to the airborne snow mass. In the model the particle ratio in the saltation layer $q_{\mathrm{salt}}$ (kg kg${}^{-\mathrm{1}}$; mass of saltating snow particles per unit mass of atmosphere) is parameterized as a function of the excess of shear stress responsible for removal of snow particles from the surface following the expression of :

5$$q_{\mathrm{salt}}={\displaystyle \frac{e_{\mathrm{salt}}}{g{h}_{\mathrm{salt}}}}\left(u_{*}^{\mathrm{2}}-u_{*\mathrm{t}}^{\mathrm{2}}\right),$$ with $e_{\mathrm{salt}}=\mathrm{1}/\left(\mathrm{3.25}{u}_{*}\right)$ the saltation efficiency expressed as a dimensionless coefficient inversely proportional to the friction velocity, $g$ $=$ 9.81 the gravitational acceleration (m s${}^{-\mathrm{2}}$) and $h_{\mathrm{salt}}=\mathrm{0.08436}{u}_{*}^{\mathrm{1.27}}$ the thickness of the saltation layer (m) as proposed by . Note that saltation is not explicitly resolved by the model, and $q_{\mathrm{salt}}$, considered constant throughout the saltation layer, only serves as a lower boundary condition for the suspension layer.

The formulation of $q_{\mathrm{salt}}$ is given for stationary conditions. Although the non-linear relationship between the flow and snow mass flux can induce fluctuations in particle concentration in non-stationary-conditions , numerical simulations suggest that steady-state saltation is achieved within an interval of a few seconds , which is well below the model time step of 60 s. With wind speeds at the surface typically reaching 10 m s${}^{-\mathrm{1}}$ during drifting-snow occurrences (Fig. S3 in the Supplement), this corresponds to characteristic lengths of a few tens of metres, i.e. 3 orders of magnitude lower than the horizontal resolution of 10 km used in this study, indicating that a formulation of a stationary particle ratio remains appropriate in this context.

Snow transport in turbulent suspension is computed by the three-dimensional turbulence and advection schemes of MAR, enabling a discretization of snow transport profiles along the vertical grid of the model. The suspension layer receives contributions from both cloud and eroded particles advected from overlying and upwind grid cells (top and lateral influx). Snow particles are suspended in the first model level through diffusion from the saltation layer (bottom influx). The mass actually removed from the surface then corresponds to the upward mass exchange between the saltation and the suspension layers and is expressed by the surface turbulent flux of snow particles $u_{*}{q}_{\mathrm{s}*}$ (m s${}^{-\mathrm{1}}$ kg kg${}^{-\mathrm{1}}$) assuming that it follows a bulk flux formulation: 6$$u_{*}{q}_{\mathrm{s}*}=C_{\mathrm{D}}U\mathit{\zeta }\left(q_{\mathrm{s}}-q_{\mathrm{salt}}\right),$$ where $q_{\mathrm{s}*}$ is the turbulent scale for snow particles (kg kg${}^{-\mathrm{1}}$), $\mathit{\zeta }$ is the ratio of eddy diffusivities for suspended particles and momentum, and $q_{\mathrm{s}}$ is the snow particle ratio (kg kg${}^{-\mathrm{1}}$) taken at the lowest model level. Because of fragmentation upon repeated collision between each other and with the snow surface, drifting-snow particles have smaller radius and thus have smaller settling velocities than snow particles that have not yet experienced contact with the surface . As MAR currently does not distinguish snow particles originating from the surface from those directly formed by its cloud microphysics, the factor $\mathit{\zeta }$ was introduced in the original formulation of to enhance the upward turbulent particle transport and compensate for a likely overestimation of the settling velocity of drifting-snow particles . Following , $\mathit{\zeta }$ is set to 3.

3.3 Interactions with the atmosphere

Eroded snow is transmitted to the atmosphere by the surface scheme and added to the pre-existing airborne snow mass without distinguishing the source of particles. Thermodynamic interactions of airborne snow particles with the atmosphere are handled by the cloud microphysical scheme. The increase in air density due to the presence of snow particles is taken into account by modifying the formulation of virtual potential temperature . Sublimation of snow particles occurs along their residence into the atmosphere and is parameterized as a function of the snow particle ratio and undersaturation of air . In particular, the model assumes an exponential size distribution of suspended (cloud and eroded) snow particles :

7$$n_{\mathrm{s}}=n_{\mathrm{0}}\exp \left(-{\mathit{\lambda }}_{\mathrm{s}}{D}_{\mathrm{s}}\right),$$ with $n_{\mathrm{s}}$ the number of snow particles of diameter $D_{\mathrm{s}}$ (10${}^{-\mathrm{3}}$ m) per unit volume, $n_{\mathrm{0}}$ an empirical constant that corresponds to the intercept parameter of the size distribution (m${}^{-\mathrm{4}}$) and ${\mathit{\lambda }}_{\mathrm{s}}$ the (dimensionless) dispersion parameter: 8$${\mathit{\lambda }}_{\mathrm{s}}={\left({\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\mathrm{p}}{n}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{a}}{q}_{\mathrm{s}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{4}}},$$ where ${\mathit{\rho }}_{\mathrm{p}}$ is the snow particle density (set to 100 kg m${}^{-\mathrm{3}}$) and ${\mathit{\rho }}_{\mathrm{a}}$ is the air density (kg m${}^{-\mathrm{3}}$). Snow particles are considered as graupel-like snow of hexagonal type, and the spectrally averaged snow particle diameter $D_{\mathrm{s}}$ is prescribed as a constant following . The latent heat consumption and humidity release caused by atmospheric sublimation are directly accounted for in the energy and mass budget of each atmospheric layer in which sublimation occurs. This ensures that the model captures (i) the negative feedback of sublimation through the increase in relative humidity of the air, (ii) advective transport of humidity, and (iii) the sublimation-induced cooling increasing air density and inhibiting upward turbulent motions . Weakening of drifting snow in response to decreasing $u_{*}$ as turbulence declines is reflected through the dependency of the surface snow turbulent flux $u_{*}{q}_{\mathrm{s}*}$ on the difference $u_{*}^{\mathrm{2}}-u_{*\mathrm{t}}^{\mathrm{2}}$.

Not distinguishing the origin of particles despite differences in shape and size between cloud and eroded snow particles can affect the estimation of sublimation according to the predominance of one type of particles over the other in the actual airborne snow mass. Overestimation of atmospheric sublimation rates by the model within drifting-snow layers actually mainly consisting of eroded particles can thus be expected, which would be however partially counterbalanced by the enhanced negative feedback of sublimation and all the less pronounced as the relative contribution of cloud particles prevails.

Airborne snow particles interact with the radiative transfer through the atmosphere and affect the surface energy budget by modulating downwelling irradiance, similarly to optically thin, low-level clouds . The representation of the radiative contribution of drifting-snow clouds is achieved in the model by including the snow particle ratio $q_{\mathrm{s}}$ in the computation of cloud optical depth and emissivity .

3.4 Interactions with the surface

Alteration of surface characteristics through erosion–deposition of snow influences in turn the occurrence of drifting snow through various feedback mechanisms. Snow deposited at the surface during drifting snow is subject to the combined actions of wind and saltation, which break original crystal shapes and favour the formation of smaller, rounded snow grains , leading to enhanced sintering, more efficient mechanical packing and increased density . Drifting-snow compaction, together with the exposure of denser snow or ice layers through erosion and/or sublimation, naturally contributes to reduce the likelihood of additional drifting snow. The erosion–deposition process also influences the surface energy budget by modifying the surface albedo, which largely determines the energy available for melting. Surface melting reduces or even inhibits the potential for erosion in summer by increasing water content, density and cohesion . Capturing these effects is thus of particular significance to account for temporal variations in drifting-snow frequency over peripheral regions of the Antarctic ice sheet .

A different feature of the current drifting-snow scheme of MAR contrasting with earlier versions is that, instead of being simultaneously distributed over several upper snow layers, the influence of erosion and deposition at each model time step (60 s) is restricted to the uppermost snow layer only, under the consideration that only the surface snowpack layer can exchange momentum and mass with the atmosphere. For deposition, this reduces the computational cost by preventing rearrangements of several snow layers per time step. For erosion, this avoids numerical instabilities related to the likely removal of several snow layers deeper in the snowpack while the computation of the surface temperature and energy balance is based on the surface layer only. Snow layers with different characteristics may thus be deposited or exposed successively at the top of the snowpack during a drifting-snow event, thus influencing the simulated surface albedo.

The current version uses fixed values for the characteristics of deposited snow but implicitly accounts for differences between cloud and eroded particles. The characteristics of fresh snow (${\mathit{\rho }}_{\mathrm{0}}$, $d$ $=$ 1 and $s$ $=$ 0) differ from those of eroded particles which are assumed to have completely lost their initial shape through collision and sublimation during transport and to be fully rounded ($d$ $=$ 0 and $s$ $=$ 1), although numerical simulations suggest the coexistence of various particle shapes during fully developed drifting snow . Drifting snow is deposited with a density ${\mathit{\rho }}_{\mathrm{DR}}$ assumed to be that of the current surface layer, with the restriction that ${\mathit{\rho }}_{\mathrm{DR}}$ does not exceed ${\mathit{\rho }}_{max}$ to account for maximum surface snow density values observed in Antarctica and enable the deposition of snow over more compacted snow and/or ice surfaces. The surface density ${\mathit{\rho }}_{\mathrm{s}}$ is updated according to a relative contribution of both types of particles $$\begin{array}{}\text{9}& {\displaystyle }{\mathit{\rho }}_{\mathrm{s}}={\mathit{\rho }}_{\mathrm{0}}\left(\mathrm{1}-f_{\mathrm{DR}}\right)+{\mathit{\rho }}_{\mathrm{DR}}{f}_{\mathrm{DR}},\text{10}& {\displaystyle }{f}_{\mathrm{DR}}=\mathrm{1}-{\displaystyle \frac{q_{\mathrm{s},z_{\lim }}}{q_{\mathrm{s}}}},\end{array}$$ where $f_{\mathrm{DR}}$ is the drifting-snow fraction varying between 0 and 1, and $q_{\mathrm{s},z_{\lim }}$ is the snow particle ratio at the atmospheric level closest to $z_{\lim }$ (m), where the contribution of eroded particles to the mass ratio is assumed to be negligible compared to the contribution of snowfall. A value of 100 m above surface has been adopted for $z_{\lim }$ in accordance with the average depth of drifting-snow layers over Antarctica as retrieved from remote sensing techniques . Snowfall conditions imply $q_{\mathrm{s},z_{\lim }}\sim q_{\mathrm{s}}$ and low values of $f_{\mathrm{DR}}$, and snow is deposited at the surface with a predominant contribution of ${\mathit{\rho }}_{\mathrm{0}}$. Conversely, drift conditions imply $q_{\mathrm{s},z_{\lim }}\ll q_{\mathrm{s}}$ and high values of $f_{\mathrm{DR}}$, and ${\mathit{\rho }}_{\mathrm{s}}$ tends towards ${\mathit{\rho }}_{\mathrm{DR}}$.

The post-depositional increase in snow density through wind hardening is accounted for in the model by increasing the density of the uppermost snowpack layer in each grid cell exposed to drifting snow. The temporal evolution of surface density along the range of values for which snow remains erodible is parameterized according to a linear densification rate from the fresh snow value ${\mathit{\rho }}_{\mathrm{0}}$ (assumed to be representative of snow that has been barely altered by post-depositional processes) to the prohibitive density value for snow erosion ${\mathit{\rho }}_{max}$, i.e.

11$${\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{s}}}{\mathrm{d}t}}={\displaystyle \frac{{\mathit{\rho }}_{max}-{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{\tau }}_{\mathrm{DR}}}},$$ in which the characteristic timescale for drifting-snow compaction ${\mathit{\tau }}_{\mathrm{DR}}$ is set to 24 h. This value corresponds to the average duration of drifting-snow events reported in and is used here as the typical duration for exhaustion of erodible snow to be reached. The linear behaviour of the densification rate follows the linear increase in surface snow density retrieved from measurements performed during a drifting-snow event in Adélie Land (Fig. S2). Further details on this experiment are provided in the Supplement (Sect. S1). By fixing ${\mathit{\rho }}_{\mathrm{0}}$ and parameterizing $u_{*\mathrm{t}}$ as an increasing function of ${\mathit{\rho }}_{\mathrm{s}}$ (Eq. ), Eq. () does not necessarily enable a correspondence with actual snow surface densities but rather merely ensures a realistic time evolution of surface snow density. It also prevents large (positive) values of the difference $u_{*}-u_{*\mathrm{t}}$ to endure through time and thus acts as a negative feedback for snow erosion.

3.5 Erosion

For each continental grid cell MAR calculates the actual snow mass eroded from the snowpack (ER, kg m${}^{-\mathrm{2}}$) during the current time step according to the following chain of events (Fig. ).

• The snow particle ratio $q_{\mathrm{s}}$ in the first model level, which includes the contributions of snowfall, atmospheric sublimation and advection of snow as computed by the cloud microphysical and turbulence schemes, is transmitted to the surface scheme.

• A potential maximum erosion ${\mathrm{ER}}_{max}$ (kg m${}^{-\mathrm{2}}$) is estimated from the surface turbulent flux of snow particles $u_{*}{q}_{\mathrm{s}*}$ computed from step 7 at the previous time step: ${\mathrm{ER}}_{max}={\mathit{\rho }}_{\mathrm{a}}{u}_{*}{q}_{\mathrm{s}*}\mathrm{d}t$.

• Actual erosion ER ($\ge$ 0) is calculated from removal of snow from the surface snowpack layer until ${\mathrm{ER}}_{max}$ is reached or the layer has been entirely eroded.

• Snow at the surface densifies following Eq. ().

• The drift fraction is obtained from Eq. (). Snow is deposited at the surface and surface density is adjusted according to Eq. ().

• The threshold friction velocity $u_{*\mathrm{t}}$ and the saltation particle ratio $q_{\mathrm{salt}}$ are deduced from Eqs. () and ().

• The surface turbulent flux of snow particles $u_{*}{q}_{\mathrm{s}*}$ is computed from Eq. ().

• The contribution of erosion is added to $q_{\mathrm{s}}$, which is then transmitted to the turbulence, cloud microphysical, and radiative schemes to compute advection and interactions of airborne snow with the atmosphere.

3.6 Surface roughness

Drifting snow is responsible for the development of surface microrelief, whose spatial arrangement combined with the orientation of the wind determines the roughness length for momentum $z_{\mathrm{0}}$. Because of wind-driven reshaping of the snow topography and the diversity in surface types, $z_{\mathrm{0}}$ varies by several orders of magnitude with time and space across the Antarctic continent . The sensitivity of MAR to the parameterization of surface roughness has been discussed in , who demonstrated that a dynamic representation of $z_{\mathrm{0}}$ is required to improve the modelling of wind speed and drifting-snow fluxes in Adélie Land. Due to inconsistencies between observed and modelled temporal variations in $z_{\mathrm{0}}$ values, the former parameterization of $z_{\mathrm{0}}$ was changed and is now computed as a function of the air temperature only . This parameterization was developed so that $z_{\mathrm{0}}$ fits the observed seasonal variations between high ($>$ 10${}^{-\mathrm{3}}$ m) summer and lower winter values in coastal Adélie Land, for air temperatures above $-$20 ${}^{\circ }$C. For lower temperatures, a constant $z_{\mathrm{0}}$ value of 2 $\times$ 10${}^{-\mathrm{4}}$ m is set in agreement with observations on the Antarctic Plateau .

Roughness features alter the spatial distribution of wind shear near the surface through pressure fluctuation gradients in their immediate vicinity (i.e. the form drag). This drag partitioning results in a loss of momentum from the near-surface flow by turbulent friction, which in turn reduces the energy budget available for erosion in the form of negative feedback. Previous versions of the drifting-snow scheme in MAR included a parameterization of drag partitioning developed for non-erodible roughness elements encountered in desert-like environments . While snow roughness features have been shown to effectively exert significant form drag inhibiting snow erosion, this mainly occurs for near-surface flows and roughness features of crosswise orientations and can essentially vanish through a rapid streamlining process of the microrelief under erosive conditions . Sensitivity experiments revealed that the drag partition scheme was responsible for a strong inhibiting effect on snow erosion beyond the observed magnitude of the negative feedback mechanism, possibly as it does not account for the dynamical and erodible nature of snow microrelief . Consequently, it has been disabled in the current model version.

4 Field area, observation data and evaluation methods

The near-surface climate in coastal Adélie Land is dominated by strong katabatic flows which drain cold air from the continental interior toward the steep coastal escarpment, enabling the regular incidence of well-developed drifting-snow events throughout the year . High erosion rates and export of drifting snow combine to melt and sublimation to produce local net ablation at the surface and resulting persistent blue-ice areas near the coast on the steepest part of the ice margin . In a fairly narrow transition, net accumulation is observed a few kilometres inland despite significant drifting snow , with annual SMB values displaying a high, kilometre-scale variability as a result of wind redistribution .

Table 1

Main characteristics of the two measurement locations used for the evaluation of MAR.

Station	Location	Elevation	Elevation	Observation period
		(m)	bias (m)
D47	67.4${}^{\circ }$ S, 138.7${}^{\circ }$	1560	$-$8	January 2010–December 2012
D17	66.7${}^{\circ }$ S, 139.9${}^{\circ }$	450	$-$66	February 2010–December 2018${}^{\ast }$

${}^{\ast }$ Station still operative.

The performance of MAR in reproducing the drifting-snow climate of Adélie Land is evaluated against 2 m wind speed and direction, air temperature, and air relative humidity observations collected at two locations 100 km apart, D47 and D17 (Fig. , Table ). Data at half-hourly intervals are available from automatic weather stations installed at both sites and operated by the Institut des Géosciences de l'Environnement (IGE) over the periods 2010–2012 for D47 and 2010–2018 for D17 . All datasets are reported as quality-controlled. At each station relative humidity is originally given with respect to water and has been converted to be expressed with respect to ice in subfreezing conditions following formulae and using the 2 m temperature record in the conversion. Climate variables extracted from the lowest model level (2 m) and the nearest grid cell to the observation location are used for comparison.

Half-hourly meteorological data also include drifting-snow mass fluxes measured almost continuously over their respective observation periods using acoustic second-generation FlowCapt™ sensors (hereafter referred to as 2G-FlowCapt™). The instrument consists of a 1 m long tube containing electroacoustic transducers that convert the acoustic vibration caused by the drifting-snow particles colliding with the tube into a snow mass flux integrated over the exposed length of the tube. At D47, two 2G-FlowCapt™ sensors were installed and superimposed vertically to sample the first 2 m above ground and detect the onset of drifting snow. Site D17 was initially equipped with only one 2G-FlowCapt™ set up close to the surface and completed with a second instrument in December 2012 to match the configuration of D47. As snow particles directly originating from clouds cannot be discriminated from saltating and/or suspended snow particles relocated from the ground, measured snow mass fluxes account for all forms of drifting snow along the sampling height. The IGE database then enables the evaluation of simulated drifting snow and the relative climate against time-averaged measurements of bulk flow and mass flux conditions, consistently with the steady-state drifting-snow physics implemented in MAR. An extensive description of the drifting-snow and meteorological equipment at D47 and D17 can be found in .

Thorough evaluation of drifting snow requires consistency between observed and modelled drifting-snow mass fluxes. Near-surface snow mass fluxes simulated by MAR can be inferred from $q_{\mathrm{s}}$ at the lowest vertical model level (2 m), which, similarly to the 2G-FlowCapt™, does not distinguish the origin of particles. By approximating the mean speed of suspended particles with the mean wind speed $U$ at the lowest model level, an average, horizontal (vertically integrated), near-surface drifting-snow mass flux ${\mathit{\mu }}_{\mathrm{MAR}}$ (kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) can be expressed as 12$${\mathit{\mu }}_{\mathrm{MAR}}=U{\mathit{\rho }}_{\mathrm{a}}{q}_{\mathrm{s}}.$$ With a lowest level at 2 m height, MAR does not capture the strong exponential decrease in snow mass flux with height existing close to the surface e.g.. As the snow particle ratio $q_{\mathrm{s}}$ has the same value throughout the model level, comparison between model and observations is performed by combining, when available, snow mass flux estimates at both measurement levels into an average, near-surface, drifting-snow mass flux ${\mathit{\mu }}_{\mathrm{OBS}}$ (kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) calculated through 13$${\mathit{\mu }}_{\mathrm{OBS}}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{1}}{h}_{\mathrm{1}}+{\mathit{\mu }}_{\mathrm{2}}{h}_{\mathrm{2}}}{h_{\mathrm{1}}+h_{\mathrm{2}}}},$$ where ${\mathit{\mu }}_{\mathrm{i}}$ is the observed snow mass flux integrated over the exposed length $h_{\mathrm{i}}$ of the corresponding 2G-FlowCapt™ sensor .

Modelled SMB is compared with observations obtained from annual measurements of 91 snow stakes distributed every $\sim$ 1.5 km along a 150 km long transect (see for description) that extends from the coast up to 1800 m a.s.l. and crosses the locations of D17 and D47. The SMB dataset includes observations collected over the period 2004–2018 and covers the strong SMB gradient existing between the coast and relatively drier and colder conditions inland. All annual observed values contained in each model grid cell (10 $\times$ 10 km${}^{\mathrm{2}}$) are averaged to produce a mean observed SMB value per grid cell.

5 Model evaluation

5.1 Near-surface climate

Accurate near-surface flows are required for a realistic representation of drifting-snow processes. The modelled annual mean (2004–2018) 2 m wind field (Fig. a) shows dominant southeasterly katabatic flows all over the integration domain. Modelled annual mean wind directions at D47 (149${}^{\circ }$) and D17 (162${}^{\circ }$) agree within less than 10${}^{\circ }$ with observations . Local flow acceleration causes the highest annual mean near-surface wind speeds to be simulated in confluent topography and over 150 km along the coast east of the D17 location. The occurrence of maximum wind jets along this coastal section of Adélie Land is a well-known feature of the study area supported by modelling efforts and observations e.g..

Figure 2

(a) Simulated annual mean (2004–2018) 2 m wind speed (colours) and direction (arrows) and (b) surface snow density (uppermost snowpack layer) on the model integration domain after removal of the relaxation zone (10 grid cells). The thin black curves show the ice sheet topography at 250 m contour intervals. The position of the model domain on the Antarctic ice sheet is indicated in the inset. The locations of measurement sites D47 and D17 are marked with red crosses.

[Figure omitted. See PDF]

Snow at the surface densifies with drifting snow (Eq. ), causing the spatial distribution of ${\mathit{\rho }}_{\mathrm{s}}$ to be linked with the variability in near-surface wind speed (Fig. b). The lowest annual mean density values amount to 360 kg m${}^{-\mathrm{3}}$ in areas of low wind speed where melting is absent and increase in areas of flow acceleration where drifting-snow compaction is the most active. The highest ${\mathit{\rho }}_{\mathrm{s}}$ values around 480 kg m${}^{-\mathrm{3}}$ are produced at the coast where both drifting-snow compaction and (summer) melting occur.

Figure 3

Density scatter plots of observed vs. simulated half-hourly (a, b) wind speed, (c, d) air temperature and (e, f) air relative humidity with respect to ice at 2 m height for stations (a, c, e) D47 and (b, d, f) D17. The coloured lines show the 1 : 1 line (dashed black) and the best linear fit (red).

[Figure omitted. See PDF]

Consistent evaluation of model results for drifting-snow applications requires a fine temporal resolution since drifting snow responds to high-frequency fluctuations in wind speed. At half-hourly resolution, MAR generally underestimates 2 m wind speeds at D47, particularly above 10 m s${}^{-\mathrm{1}}$, with a mean bias of $-$1.5 m s${}^{-\mathrm{1}}$ and a root-mean squared error (RMSE) of 2.2 m s${}^{-\mathrm{1}}$ (Fig. a). Strong winds are better captured at D17 at the expense of weak wind conditions and model exactitude, resulting in a positive mean bias of $+$1.3 m s${}^{-\mathrm{1}}$ and a RMSE of 3.2 m s${}^{-\mathrm{1}}$ (Fig. b). Note that these statistics, despite data quality control, might still be affected by some measurements likely subject to instrument malfunction (icing), especially at D17 due to the proximity of the ocean, when observed half-hourly wind speeds near zero are reported concurrently with considerably higher simulated values. Actual overestimation of low wind speed values would however not be expected to be detrimental since observations show that drifting snow is usually triggered by wind speeds above 5 m s${}^{-\mathrm{1}}$ at both locations, a feature well reproduced by the model (Fig. S3). Data dispersion then reduces with increasing wind speeds.

The general underestimation in near-surface wind speed at D47 could be caused by the temperature-dependent parameterization of $z_{\mathrm{0}}$, locally still yielding too high values, while at D17 Fig.  illustrates that modelled $z_{\mathrm{0}}$ values are closer to observations. Another explanation potentially involved in the underestimation of wind speed maxima is the model misrepresentation of large eddies in case of strong winds. The local turbulence scheme of MAR is adapted for stable atmospheric boundary layers in which small eddies develop and quickly dissipate. Due to the strong turbulent mixing induced by high wind speeds, the atmospheric boundary layer is mostly statically neutral in coastal Adélie Land , and frequent development of drifting-snow layers of several hundreds of metres in thickness in this area suggests the presence of large eddies of the height of the boundary layer. Local turbulence schemes commonly struggle to reproduce the well-mixed character of neutral atmospheric boundary layers more representative of coastal windy Antarctic regions and typically fail to represent downward entrainment of momentum by large eddies from higher atmospheric levels, leading to potentially erroneous prediction of near-surface wind gusts.

Figure 4

Observed (blue circles) and simulated (red squares) monthly median values of roughness length for momentum at D17 during year 2013. Observed values are taken from

[Figure omitted. See PDF]

.

The saturation vapour pressure of air is a strongly dependent function of temperature, and the relative humidity determines the potential for atmospheric sublimation. Adequate performance in modelling these fields is of prime importance for the representation of drifting-snow sublimation. Near-surface air temperature (Fig. c, d) is well represented at both stations (positive bias $<$ 1 ${}^{\circ }$C and centred RMSE $<$ 2 ${}^{\circ }$C). As drifting-snow sublimation is a determining contributor to the atmospheric humidity budget in the area , reasonable relative humidity statistics (Fig. e, f) despite a remaining significant dispersion along the whole range of values suggest a realistic reproduction of this process in the lowest model levels.

The vertical resolution could limit the general ability of the model to represent the atmospheric boundary layer near the coast. Refining the vertical discretization by doubling at the same time the total number of levels and the number of levels in the lowest 100 m however does not significantly improve model performance when evaluated against near-surface observations (See Sect. S2).

5.2 Drifting-snow occurrences

Since drifting snow preferentially occurs during strong wind events, general biases in the model representation of wind speed influence the representation of drifting-snow frequency (Fig. ). Monthly frequency values are computed as the fraction of half-hourly drifting-snow occurrences with a simulated or observed near-surface drifting-snow mass flux greater than 10${}^{-\mathrm{3}}$ kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ . Figure b and c illustrate that the model reproduces the observed spatio-temporal variability in monthly drifting-snow frequency, with a correlation coefficient $r$ of 0.49 and 0.65 respectively at D47 and D17. The underestimation of strong winds at D47 leads to underestimation of drifting-snow frequency (observed and modelled averages of 0.81 and 0.61), which is then better represented by the model at D17 (observed and modelled averages of 0.61 and 0.57). Nearly consistent underestimation of drifting-snow frequency at D47 could also be caused by a misrepresentation of surface snow properties and their temporal evolution. For instance, surface compaction could be locally too strong in the model (see Sect. 6.2), which would shorten the duration of events and inhibit further drifting snow until snow replenishment by snowfall, while better model agreement with the observed frequency at D17 suggests a value of ${\mathit{\tau }}_{\mathrm{DR}}$ likely better suited for this location. Note that modelled drifting snow is not the most frequent where the strongest wind speeds are found (Fig. a). Highest frequency values are simulated east of D17 and D47 over topographic crests upstream of the Mertz Glacier (67.5${}^{\circ }$ S, 144.8${}^{\circ }$ E) where orographic lifting produces enhanced precipitation rates which then contribute to favour the occurrence of drifting snow by (i) increasing the airborne snow mass independently from erosion and (ii) decreasing ${\mathit{\rho }}_{\mathrm{s}}$ and lowering the erosion threshold.

Figure 5

(a) Simulated annual mean (2004–2018) drifting-snow frequency. Monthly drifting-snow frequency from observed (circles, blue dashed curve) and simulated (squares, red solid curve) half-hourly drifting-snow mass fluxes at (b) D47 and (c) D17. Monthly frequency values are calculated as the fraction of half-hourly drifting-snow occurrences (near-surface drifting-snow mass fluxes greater than 10${}^{-\mathrm{3}}$ kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ in both the model and observations).

[Figure omitted. See PDF]

Figure  only compares monthly frequency values. Typical durations of drifting-snow events in Adélie Land, however, range from several hours to a few days at most . To evaluate the model results at higher temporal resolution closer to characteristic timescales of drifting snow, half-hourly data are used to compute the probability of detection (POD) and false alarm ratio (FAR): 14$$\mathrm{POD}=\mathrm{100}{\displaystyle \frac{a}{a+b}};\mathrm{FAR}=\mathrm{100}{\displaystyle \frac{c}{c+a}},$$ where $a$ is the number of half-hourly drifting-snow occurrences correctly simulated, $b$ the number of occurrences missed by the model and $c$ the number of occurrences simulated but not observed. The Rousseau index (RI), a measure of model predictability originally defined to assess rainfall forecasts and already used to evaluate modelled drifting-snow occurrence in an alpine context , is also calculated. RI also takes into account the number of occurrences for which the absence of drifting snow is correctly simulated, $d$, such that 15$$\mathrm{RI}=\mathrm{100}{\displaystyle \frac{ad-\frac{(b+c{)}^{\mathrm{2}}}{\mathrm{2}}}{\left(a+\frac{b+c}{\mathrm{2}}\right)\left(d+\frac{b+c}{\mathrm{2}}\right)}},$$ where RI varies between $-$100 and 100. A negative value means that the model is less successful than an estimation entirely based on climatology, 0 indicates no skill and 100 is obtained for a perfect simulation.

Positive RI values are obtained at both locations (Table ), meaning the model ability to predict drifting-snow occurrences at the monthly scale (Fig. b, c) arises from a reasonable reproduction of their actual timing and duration at the half-hourly resolution. MAR shows better results (higher POD and RI) at D17 than at D47 but also simulates more unobserved occurrences (higher FAR) that compensate for missed occurrences in the calculation of monthly frequency values.

Table 2

Statistical evaluation of the ability of MAR to simulate half-hourly drifting-snow occurrences at D47 (2010–2012) and D17 (2010–2018). See Sect. 5.2 for the definition of POD, FAR and RI.

Station	POD	FAR	RI
D47	64.5	13.4	9.1
D17	80.9	25.4	45.5

5.3 Drifting-snow transport

Qualitative evaluation of MAR against monthly cumulative near-surface drifting-snow mass fluxes at D47 and D17 (Fig. b, c) reveals that the model captures the general temporal evolution of drifting-snow transport at both locations ($r$ $=$ 0.64 and 0.89). As a result of interactions between modelled near-surface flow conditions and surface snow properties, drifting snow is subject to a high variability in space and time. This variability is simulated by the model with alternating underestimation and overestimation of drifting-snow transport depending on the time period and location. Figure  compares observed and simulated snow transport during each observed drifting-snow event. As in , a drifting-snow event is defined as a period over which the observed snow mass flux is above the detection threshold of 10${}^{-\mathrm{3}}$ kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ for a minimum duration of 4 h. Discrepancy between model and observations is larger for (shorter) events of lower magnitude ($<$ 10${}^{\mathrm{3}}$ kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) and reduces for longer, more important events that contribute predominantly to the local drifting-snow transport .

Figure 6

(a) Simulated annual mean (2004–2018) near-surface drifting-snow transport. Monthly near-surface drifting-snow transport from observed (circles, blue dashed curve) and simulated (squares, red solid curve) half-hourly drifting-snow mass fluxes at (b) D47 and (c) D17.

[Figure omitted. See PDF]

Figure 7

Scatter plots of observed vs. simulated near-surface drifting-snow transport for each drifting-snow event at (a) D47 (red circles) and (b) D17 (blue squares). As in , a drifting-snow event is defined as a period over which the observed snow mass flux is above the detection threshold of 10${}^{-\mathrm{3}}$ kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ for a minimum duration of 4 h.

[Figure omitted. See PDF]

Table  compares the total snow transport during observed drifting-snow events (OBS) to that estimated by the model during these observed events (MAR${}_{\mathrm{obs}}$) and during every simulated events (MAR${}_{\mathrm{sim}}$). Referring only to observed events, MAR respectively overestimates ($+$21.6 %) and underestimates ($-$6.1 %) drifting-snow transport at D47 and D17. Together with the consistent underestimation of drifting-snow frequency at D47 (Fig. b), this means that the model simulates the main events but underestimates the occurrence of events of lower magnitude and associated transport at this location (Fig. a). However, overestimation of drifting-snow transport by the model during observed events overcompensates for the missed events, with a slight contribution of false alarms to the simulated total drifting-snow transport (MAR${}_{\mathrm{sim}}$ $\sim$ MAR${}_{\mathrm{obs}}$ $>$ OBS). At D17, while the simulated drifting-snow frequency (Fig. c) and total drifting-snow transport during observed events are in closer agreement with the observations, simulated but unobserved drifting snow is more significant (Table ) and results in overestimation of total drifting-snow transport (MAR${}_{\mathrm{sim}}$ $>$ MAR${}_{\mathrm{obs}}$ $\sim$ OBS).

Table 3

Observed and simulated total drifting-snow transport (10${}^{\mathrm{6}}$ kg m${}^{-\mathrm{2}}$) from 0 to 2 m above surface at D47 and D17. A distinction is made between total snow transport computed from observed (OBS) and simulated (MAR${}_{\mathrm{obs}}$) near-surface drifting-snow mass fluxes cumulated over each drifting-snow event identified in the database and over each simulated drifting-snow event (MAR${}_{\mathrm{sim}}$).

Station	OBS	MAR${}_{\mathrm{obs}}$	MAR${}_{\mathrm{sim}}$
D47	1.51	1.84	1.88
D17	5.25	4.92	5.4

Analysing the relationship between modelled wind speed and drifting-snow transport demonstrates that the model performance varies according to the flow conditions leading to drifting snow. Median relative biases in mean 2 m wind speed and cumulative near-surface transport during drifting-snow events are shown at both measurement sites for observed 2 m wind speed bins of 1 m s${}^{-\mathrm{1}}$ in Fig.  for the most represented wind speed categories (8 to 17 m s${}^{-\mathrm{1}}$). Since wind speed and drifting-snow transport are monotonically related, biases in event-averaged wind speed can explain biases in cumulative drifting-snow transport when both are of the same sign. Figure  shows fluctuations in the sign of the bias in drifting-snow transport when biases in mean wind speed remain of constant sign. This suggests that model errors in terms of drifting-snow transport can also be attributed to the erosion, microphysical and/or turbulence schemes. Nevertheless, note that uncertainties in the observations can also affect the evaluation. In particular, while 2G-FlowCapt™ sensors can detect the occurrence of snow transport with a high level of confidence , their ability to estimate drifting-snow mass fluxes remains to be assessed in Antarctic conditions .

Figure 8

Relative bias in mean wind speed (solid curves) and cumulated snow transport (dashed curves) during drifting-snow events. Relative biases are expressed as the ratio of the difference between model and observation to the observation and shown as median values within observed wind speed bins of 1 m s${}^{-\mathrm{1}}$ at D47 (red curves) and D17 (blue curves).

[Figure omitted. See PDF]

The mass transported in drifting snow has been shown to correlate with wind speed in a power-law fashion e.g.. Figure , in which the modelled annual mean of horizontal near-surface drifting-snow transport is studied as a function of 2 m wind speed for all continental grid cells, shows that this relationship is well reproduced by the model. This enables us to assess the plausibility of the model results in the absence of observations in other locations of the integration domain and explains the spatial distribution of drifting-snow transport, with maximum values generally simulated in areas of highest 2 m wind speeds (Fig. a). The dispersion around a given wind speed value (Fig. ) results from variations in the erosion threshold . Departure from the power-law relation at low wind speeds is caused by a predominant contribution of snowfall to the airborne snow mass and corresponds to areas of near-surface flow deceleration, such as over the Mertz Glacier ice tongue or in topographic depressions.

Figure 9

Simulated mean annual (2004–2018) near-surface drifting-snow transport as a function of 2 m wind speed at all continental grid cells.

[Figure omitted. See PDF]

5.4 Surface mass balance

In the current version of MAR, cloud, eroded and deposited snow particles are all included in the snow particle ratio $q_{\mathrm{s}}$. Precipitation, erosion and deposition can occur simultaneously and repeatedly within the same grid cell. Simulated snowfall (SF) and erosion (ER) amounts then respectively account for the cumulative snow mass that is transferred to and removed from the surface. Each of the two components thus cannot be used individually to determine the integrative contribution of snowfall, erosion and deposition separately when the drifting-snow scheme is switched on. From that perspective an erosion–deposition index $i_{\mathrm{E}/\mathrm{D}}$ is defined as

16$$i_{\mathrm{E}/\mathrm{D}}={\displaystyle \frac{\mathrm{ER}}{\mathrm{SF}}},$$ which reflects the relative local proportion of eroded to deposited snow mass.

Figure  compares the simulated and observed mean annual SMB and shows the elevation and $i_{\mathrm{E}/\mathrm{D}}$ along the stake transect whose location is marked with black dots on Fig. . MAR represents the general variability in SMB with a strong increase over the first tens of kilometres from the coast and less variability further inland. The variability in $i_{\mathrm{E}/\mathrm{D}}$ is more pronounced where the terrain is steeper near the coast and exhibits more variability in topographic surface slope, suggesting that SMB variability is driven by drifting snow. The mean SMB bias is negative with $\sim$ 67 kg m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$ (20 %) but is mainly due to underestimation of the simulated SMB near the coast. This results from either overestimation of erosion or underestimation of precipitation, or a combination of both, leading to the highest values of $i_{\mathrm{E}/\mathrm{D}}$ along the transect. The lowest, negative simulated SMB value is not found at the coast but a few grid cells inland where $i_{\mathrm{E}/\mathrm{D}}$ is maximal, i.e. where erosion amounts to 93 % of the snow mass accumulated locally. Simulated surface sublimation (and negligibly runoff) also contributes but to a lesser extent to net ablation at the surface (not shown). More generally, ablation is simulated where $i_{\mathrm{E}/\mathrm{D}}$ displays maximum values (Fig. ). This is in accordance with the literature, which suggests that drifting-snow processes are the main driver behind the formation of blue-ice areas in Adélie Land .

Figure 10

(a) Simulated (squares, red solid curve) vs. observed (circles, blue dotted curve) mean annual SMB (2004–2018) along the stake transect. The location of the SMB observations is marked with black dots in Fig. . Observed annual SMB values are averaged on MAR grid cells with no interpolation or weighting. The spatial resolution is 10 km. Distance along the transect starts at the coast and is computed as the average distance of all annual observations contained in each grid cell. Uncertainty of observed SMB (blue shaded area) is approximated by the standard deviation of all annual observations contained in each grid cell and is shown together with the mean annual (2004–2018) raw SMB signal along the transect (thin black curve) as an illustration of subgrid variability. (b) Number of cumulated annual observations per grid cell. (c) Elevation. (d) Erosion–deposition index.

[Figure omitted. See PDF]

Figure 11

(a) Simulated mean annual (2004–2018) SMB and (b) corresponding erosion–deposition index.

[Figure omitted. See PDF]

Although modelled and observed SMB values are in relative good agreement, the raw observed SMB signal (i.e. not spatially averaged over model grid cells) highlights an important subgrid variability in SMB gradients that is still not resolved at 10 km resolution (Fig. ). It also shows that the modelled SMB decrease a few kilometres upstream of the coastline is also retrieved in the observations but occurs more locally. This suggests that an improved representation of the spatial variability in SMB in the area with the model would require an even higher resolution to better resolve the interactions between the atmosphere and surface topography.

Model limitations also relate to the size of the integration domain since the model is not driven at its lateral boundaries by drifting-snow fluxes. This might affect the results, notably the location of ablation areas close to or directly downwind of the model boundaries where drifting snow can only diverge in the absence of possible advection from upwind areas outside of the integration domain.

6 Sensitivity analysis

6.1 Sensitivity to the model version

Improvements in the simulation of drifting snow with the newest version of MAR are illustrated in Fig.  for the month of January 2011 at site D47, as studied in . Despite a coarser horizontal resolution (10 km vs. 5 km) and knowing that different reanalyses (ERA5 vs. ERA-Interim) were used as forcing, MAR displays enhanced capabilities in simulating the timing, duration and magnitude of drifting snow (Fig. a) relative to the previous model version and setup (MARv2) used in . Since the representation of 2 m wind speed only slightly improves with the current version (not shown), the better agreement with observations can be mainly attributed to a higher performance of the drifting-snow scheme. Improved drifting snow in MAR also leads to increased 2 m relative humidity as a result of enhanced atmospheric sublimation (Fig. b). Periods of near-saturated conditions, whose occurrence corresponds with that of drifting snow, are better captured in MAR than in the former model version.

Figure 12

(a) Near-surface drifting-snow mass flux and (b) 2 m relative humidity in MARv2 (green) and MARv3.11 (red) vs. observations (blue) for the case study in January 2011 at D47 discussed in . Note that MARv3.11 is run at 10 km resolution and forced by ERA5, while MARv2 is run at 5 km and forced by ERA-Interim.

[Figure omitted. See PDF]

6.2 Sensitivity to input parameters

The performance of MAR is dependent on the choices made for several input parameters of the drifting-snow scheme. Experiments were undertaken to assess the sensitivity of the drifting-snow scheme to changes in the fresh snow density ${\mathit{\rho }}_{\mathrm{0}}$ and characteristic timescale for drifting-snow compaction ${\mathit{\tau }}_{\mathrm{DR}}$ over a range of otherwise plausible values. Sensitivity to the disabling of drifting-snow compaction and to the parameterization of the threshold friction velocity for snow erosion $u_{*\mathrm{t}}$ and the saltating particle ratio $q_{\mathrm{salt}}$ was also explored. Simulated drifting-snow frequency and transport are evaluated for each experiment at site D47 for year 2010 and D17 for year 2013, restarting from the simulation obtained with the control setup. Results are summarized in Table .

Table 4

Sensitivity analysis to various model parameters of the annual drifting-snow frequency and accumulated drifting-snow transport up to 2 m above the surface in 10${}^{\mathrm{5}}$ kg m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$ during year 2010 at D47 and 2013 at D17. Values in brackets refer to the relative difference in percent compared to observations. L07 and V12 denote the adaptation of the parameterization of $u_{*\mathrm{t}}$ according to and respectively (see Appendix A). B00 denotes the adaptation of the parameterization of $q_{\mathrm{salt}}$ according to (see text).

Experiment	D47		D17
	Frequency	Snow transport	Frequency	Snow transport
2G-FlowCapt™	0.84	8.28	0.61	9.15
Control: ${\mathit{\rho }}_{\mathrm{0}}$ $=$ 300 kg m${}^{-\mathrm{3}}$, ${\mathit{\tau }}_{\mathrm{DR}}$ $=$ 24 h	0.63	8.81 ($+$6.4 %)	0.60	9.71 ($+$6.1 %)
${\mathit{\rho }}_{\mathrm{0}}$ $=$ 250 kg m${}^{-\mathrm{3}}$	0.84	13 ($+$57 %)	0.34	4.54 ($-$50.4 %)
${\mathit{\rho }}_{\mathrm{0}}$ $=$ 350 kg m${}^{-\mathrm{3}}$	0.91	15.4 ($+$86 %)	0.89	15.1 ($+$65 %)
${\mathit{\tau }}_{\mathrm{DR}}$ $=$ 12 h	0.59	7.90 ($-$4.6 %)	0.55	8.88 ($-$3 %)
${\mathit{\tau }}_{\mathrm{DR}}$ $=$ 48 h	0.72	10.81 ($+$30.6 %)	0.68	11.37 ($+$24.3 %)
${\mathit{\rho }}_{\mathrm{DR}}$ $=$ 300 kg m${}^{-\mathrm{3}}$, ${\mathit{\tau }}_{\mathrm{DR}}$ $=$ $\mathrm{\infty }$	0.98	17.52 ($+$111.6 %)	0.97	16.47 ($+$80 %)
$u_{*\mathrm{t}}$ $=$ L07	0.46	4.71 ($-$43.1 %)	0.32	3.94 ($-$56.9 %)
$u_{*\mathrm{t}}$ $=$ V12	0.46	4.50 ($-$45.6 %)	0.31	3.95 ($-$56.8 %)
$q_{\mathrm{salt}}$ $=$ B00	0.63	8.89 ($+$7.4 %)	0.60	9.82 ($+$7.3 %)

The model exhibits a high sensitivity to ${\mathit{\rho }}_{\mathrm{0}}$. This parameter determines the density at which fresh snow is deposited but also intervenes in the definition of $u_{*\mathrm{t}}$ (Eq. ) and drifting-snow compaction rate (Eq. ). To conserve the physical consistency of the drifting-snow scheme, changes in ${\mathit{\rho }}_{\mathrm{0}}$ must be applied to the whole set of parameterizations and necessarily lead to modifications that result from the combined sensitivities to each of the three parameterizations. This requires notably the adoption of values for ${\mathit{\rho }}_{\mathrm{0}}$ that are probably above typical values for fresh snow, as $u_{*\mathrm{t}}$ becomes a more restrictive criterion for erosion with lower values for ${\mathit{\rho }}_{\mathrm{0}}$ (Fig. ), and consequently that already partially integrate the effect of post-depositional processes. Contrasting results between D47 and D17 demonstrate that the model response can be spatially heterogeneous depending on dominant local effects. Reducing (increasing) ${\mathit{\rho }}_{\mathrm{0}}$ to 250 (350) kg m${}^{-\mathrm{3}}$ at D17 leads to lower (higher) snow density values at deposition by snowfall but also increases (decreases) $u_{*\mathrm{t}}$, resulting in a general dampening (enhancement) of simulated drifting snow around 50 % in both frequency and transport. This indicates that lower surface snow density ${\mathit{\rho }}_{\mathrm{s}}$ is locally overcompensated for by a corresponding strong increase in $u_{*\mathrm{t}}$. The model would become less and less sensitive to ${\mathit{\rho }}_{\mathrm{0}}$ values below 250 kg m${}^{-\mathrm{3}}$ because $u_{*\mathrm{t}}$ would more rapidly reach high prohibitive values for snow erosion. The reverse situation is depicted at D47, where the model simulates lower ${\mathit{\rho }}_{\mathrm{s}}$ (Fig. b). Lowering ${\mathit{\rho }}_{\mathrm{0}}$ to 250 kg m${}^{-\mathrm{3}}$ induces a local decrease in ${\mathit{\rho }}_{\mathrm{s}}$ that improves drifting-snow frequency but produces drifting-snow transport beyond observed values ($>$ 50 %) despite higher $u_{*\mathrm{t}}$. Overestimated and similar drifting-snow features are obtained at both locations with ${\mathit{\rho }}_{\mathrm{0}}$ $=$ 350 kg m${}^{-\mathrm{3}}$. Discussing the influence of ${\mathit{\rho }}_{\mathrm{0}}$ beyond 350 kg m${}^{-\mathrm{3}}$ would implicitly imply deposition of fresh snow in a well-advanced state of compaction and is therefore not considered.

Figure 13

Threshold friction velocity ($u_{*\mathrm{t}}$) as a function of surface snow density (${\mathit{\rho }}_{\mathrm{s}}$) for three sets of parameterizations and different fresh snow density values (${\mathit{\rho }}_{\mathrm{0}}$). A constant roughness length for momentum $z_{\mathrm{0}}$ of 10${}^{-\mathrm{3}}$ m is assumed here for the computation of $u_{*\mathrm{t}}$ (see Sect. 3). The solid blue line highlights the current control parameterization of $u_{*\mathrm{t}}$ in MAR (${\mathit{\rho }}_{\mathrm{0}}$ $=$ 300 kg m${}^{-\mathrm{3}}$).

[Figure omitted. See PDF]

Changes in ${\mathit{\tau }}_{\mathrm{DR}}$ alter the timing and duration of drifting-snow events but appear to exert a moderate influence on modelled drifting snow. Halving and doubling ${\mathit{\tau }}_{\mathrm{DR}}$ to 12 and 48 h respectively decreases and increases local drifting-snow frequency and transport with a relatively lower magnitude ($<$ 30 %) than the tested changes in ${\mathit{\rho }}_{\mathrm{0}}$. The value of ${\mathit{\tau }}_{\mathrm{DR}}$ $=$ 24 h is of local significance since it has been scaled on the observed median duration of drifting-snow events in Adélie Land. Different timescales for drifting-snow compaction could be expected under different environmental conditions, particularly for different wind speed regimes. This would suggest that ${\mathit{\tau }}_{\mathrm{DR}}$ should be varying regionally and, if needed, could be adapted to more practical values depending on the local climate of the region of interest.

Disabling drifting-snow compaction by setting ${\mathit{\tau }}_{\mathrm{DR}}$ to infinity and the density of deposited drifting snow ${\mathit{\rho }}_{\mathrm{DR}}$ to the standard ${\mathit{\rho }}_{\mathrm{0}}$ value of 300 kg m${}^{-\mathrm{3}}$ produces almost permanent drifting snow and leads to the strongest overestimation ($>$ 100 % at D47) of drifting-snow transport among the sensitivity experiments. These model results illustrate the necessity of taking the negative feedback of drifting-snow compaction into account to simulate drifting-snow frequency and mass fluxes satisfactorily at the two measurement sites and prevent removal of almost all the accumulated snow.

Conversely, the strongest inhibition of drifting snow (by 50 %) is obtained through changes in $u_{*\mathrm{t}}$ following the parameterizations of and , hereafter referred to as L07 and V12. Appendix A gives numerical details. The resemblance in the evolution of $u_{*\mathrm{t}}$ with ${\mathit{\rho }}_{\mathrm{s}}$ (Fig. ) explains the similar model results obtained at both stations with L07 and V12. However, both parameterizations do not depend on the same set of variables. L07 only depends on ${\mathit{\rho }}_{\mathrm{s}}$ and enables investigation of the influence of ${\mathit{\rho }}_{\mathrm{0}}$ without inherently altering $u_{*\mathrm{t}}$ or ${\mathit{\tau }}_{\mathrm{DR}}$. Overall, L07 is a stricter criterion for snow erosion than the control parameterization to the extent that it cannot be offset by prescribing low ${\mathit{\rho }}_{\mathrm{0}}$ values or by increasing ${\mathit{\tau }}_{\mathrm{DR}}$ within the range of tested values (not shown). V12 corresponds to a weighting of the erodibility index $i_{\mathrm{ER}}$ by ${\mathit{\rho }}_{\mathrm{s}}$ and remains dependent on ${\mathit{\rho }}_{\mathrm{0}}$. Like L07, V12 yields a general increase in $u_{*\mathrm{t}}$ that becomes very restrictive and that a change in ${\mathit{\tau }}_{\mathrm{DR}}$ can hardly compensate for without increasing it substantially or without increasing ${\mathit{\rho }}_{\mathrm{0}}$ (not shown). The resulting inhibition at D17 is also quite close to that obtained with ${\mathit{\rho }}_{\mathrm{0}}$ $=$ 250 kg m${}^{-\mathrm{3}}$, where ${\mathit{\rho }}_{\mathrm{s}}$ is high enough for $u_{*\mathrm{t}}$ to drive the model response to changes in ${\mathit{\rho }}_{\mathrm{0}}$.

The model is virtually insensitive to changes in the formulation of the saltating particle ratio $q_{\mathrm{salt}}$. Equation () gives a maximal value of $q_{\mathrm{salt}}$ for intermediate values of the friction velocity $u_{*}$. This is a potential limitation of this formulation that has been identified by . tested another formulation, hereafter referred to as B00, in which the efficiency of saltation is assumed to be constant ($e_{\mathrm{salt}}$ $=$ 0.535), yielding a relation for $q_{\mathrm{salt}}$ that increases monotonically with $u_{*}$. Figure S4 illustrates that a significant sensitivity to the parameterization of $q_{\mathrm{salt}}$ could be expected, especially for $u_{*}$ well above 0.6 m s${}^{-\mathrm{1}}$, from which the control formulation and B00 start to diverge from each other. However, the comparison between the drifting-snow mass fluxes obtained with B00 and the control run shows very small differences ($<$ 2 %) in the cumulative drifting-snow mass transport at D47 and D17 and similar drifting-snow frequency values. This result is relative to the current model version and is caused by the limitation made on the erosion of only the surface layer per model time step (see Sect. 3.4), which is strong enough to reduce the sensitivity to the formulation of $q_{\mathrm{salt}}$.

Table  shows that the best agreement between model and observations is obtained with the control setup. However, this agreement also most likely results from error compensations which are difficult to quantify and identify in the model. Although model results remain in the same order of magnitude as observations in most sensitivity experiments, the simulated drifting-snow frequency and transport appear very sensitive to input parameters, in particular to the choice of ${\mathit{\rho }}_{\mathrm{0}}$ and the parameterization of $u_{*\mathrm{t}}$ with the control setup, which inherently also affect the SMB. It is likely that different combinations of these parameters could lead to equally satisfactory model results, providing that the SMB is also conjointly assessed to preserve consistency between erosion and accumulation rates.

7 Conclusion

A dataset of drifting-snow and SMB observations collected in Adélie Land has been used to evaluate the latest version of the regional climate model MAR (v3.11) equipped with a new drifting-snow scheme. Modifications and developments relative to the previous drifting-snow scheme mainly consisted in (i) computing the roughness length $z_{\mathrm{0}}$ as a varying function of air temperature and disabling the drag partition scheme and its inhibiting effect on snow erosion, (ii) rewriting the parameterization of the threshold friction velocity for snow erosion $u_{*\mathrm{t}}$ as a function of surface snow density only, (iii) differentiating snow density at deposition between precipitation and drifting snow, (iv) including a simple parameterization for drifting-snow compaction of the uppermost snowpack layer, and (v) restricting erosion within one model time step to the uppermost snowpack layer only. The model provides improved results compared to a previous version and demonstrates a good capability to reproduce drifting-snow occurrences and near-surface transport up to the scale of the drifting-snow event at two measurement sites 100 km apart over multiple-year periods, conjointly with the resulting variability in SMB along a 150 km long transect crossing the location of the measurement sites. These results constitute an encouraging step toward the use of MAR to study drifting snow and its influence on the Antarctic climate and SMB over a wide range of spatial and temporal scales from local-case studies to a continent-wide climatology.

Modelling near-surface drifting snow could be expected to require a particularly fine vertical resolution close to the surface to resolve the mass exchange with the saltation layer and the strongly decreasing snow mass fluxes in the first metres above ground. A lowest atmospheric level at 2 m however does not appear as a limiting factor compared to observed snow mass fluxes vertically integrated from 0 to 2 m above ground. Future developments should focus on the distinction between snow particles remobilized from the ground from precipitating particles to explicitly account for different particle properties such as snow grain shape and size or settling velocity and enable disentanglement of the contribution of snowfall from eroded particles to the local climate and drifting-snow features.

Model evaluation in this study is limited to the surface ($<$ 2 m). Snow-mass transport and atmospheric sublimation processes during drifting snow are for a large part supplied with particles originating from the surface in the lower atmosphere, so a good representation of drifting-snow features near the surface is of prime importance. The performance of the model at higher atmospheric levels where drifting-snow transport and sublimation can supposedly be of substantial significance in terms of SMB e.g. is however not known yet and remains to be assessed. Remotely sensed properties of drifting snow from space (layer depth, optical thickness) recently made available through the ICESat-2 mission at an unprecedented level of spatial and temporal resolution would be a great candidate for such an exercise.

Appendix A Additional parameterizations of the threshold friction velocity

Changes in the threshold friction velocity for snow erosion $u_{*\mathrm{t}}$ were done following the parameterizations of and . In , $u_{*\mathrm{t}}$ is parameterized a function of surface snow density ${\mathit{\rho }}_{\mathrm{s}}$ only A1$$u_{*\mathrm{t}}=\left\{\begin{array}{ll}\mathrm{0.005}\exp \left(\mathrm{0.013}{\mathit{\rho }}_{\mathrm{s}}\right)& \mathrm{300}<{\mathit{\rho }}_{\mathrm{s}}\le \mathrm{450}\\ \mathrm{0.1}\exp \left(\mathrm{0.003}{\mathit{\rho }}_{\mathrm{s}}\right)& {\mathit{\rho }}_{\mathrm{s}}\le \mathrm{300}.\end{array}\right.$$

In the definition of the erodibility index $i_{\mathrm{ER}}$ is completed with a term depending on the surface snow density $F\left({\mathit{\rho }}_{\mathrm{s}}\right)$ to account for high polar surface snow densities : $$\begin{array}{}\text{A2}& {\displaystyle }{i}_{\mathrm{ER}}=\mathrm{0.34}{i}_{\mathrm{ER}}+\mathrm{0.66}F\left({\mathit{\rho }}_{\mathrm{s}}\right),\text{A3}& {\displaystyle }F\left({\mathit{\rho }}_{\mathrm{s}}\right)=\mathrm{1.25}-\mathrm{0.0042}\left({\mathit{\rho }}_{\mathrm{s}}-\mathrm{50}\right).\end{array}$$

Both parameterizations are given as valid for ${\mathit{\rho }}_{\mathrm{s}}$ values up to 450 kg m${}^{-\mathrm{3}}$ in agreement with the maximum snow density that can be attained through drifting-snow compaction in MAR.

Code and data availability

The meteorological and drifting-snow data can be downloaded at 10.5281/zenodo.3630497 . The surface mass balance data are available at http://pp.ige-grenoble.fr/pageperso/faviervi/stakeline.php . This service is provided by the National Observation Service GLACIOCLIM (https://glacioclim.osug.fr/?lang=en, last access: 23 May 2021) piloted by the OSUG (Observatories of Earth Sciences and Astronomy of Grenoble, France). The National Observation Service GLACIOCLIM have requested that we make the data available through the dedicated web page hosted by the OSUG website. However, we have also archived the specific version of the dataset used in this study together with the MAR source code and outputs for replication of this study at 10.5281/zenodo.4314872 . The MAR source code is tagged as v3.11.3 on https://gitlab.com/Mar-Group/MARv3.7 . See http://www.mar.cnrs.fr for more information about downloading MAR (last access: 23 May 2021).

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-14-3487-2021-supplement.

Author contributions

ChA, CK and XF designed the study. ChA and VF collected data on the field. ChA implemented model developments and updates with the help of CK and XF. ChA ran the simulations, performed the analysis and wrote the first draft of the manuscript. ChA, CK, LLT, CeA, AD, VF and XF discussed and revised the manuscript.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

This work would not have been possible without the financial and logistical support of the French Polar Institute IPEV (programme CALVA-1013 and GLACIOCLIM-SAMBA-411). The authors thank all the on-site personnel in Dumont d'Urville and Cap Prud'homme for their precious help in the field, in particular Philippe Dordhain for electronic and technical support. Luc Piard and Christophe Genthon are also acknowledged for their investment in collecting data and maintaining the observation system in Adélie Land, as well as Hubert Gallée for fruitful discussions. Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.–FNRS) under grant no. 2.5020.11, and the Tier-1 supercomputer (Zenobe) of the Fédération Wallonie-Bruxelles infrastructure funded by the Walloon Region under grant agreement no. 1117545. This is PROTECT contribution no. 16.

Financial support

This research has been supported by the Fonds de la Recherche Scientifique de Belgique, the TROIS-AS project (ANR-15-CE01-0005-01) and PROTECT, a project which has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 869304.

Review statement

This paper was edited by Steven Phipps and reviewed by two anonymous referees.