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

Drifting snow, or the wind-driven transport of snow particles and their concurrent sublimation, is a poorly documented process on the Antarctic ice sheet, 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 modifica5 tions relative to previous version have been implemented and include notably a parameterisation for drifting-snow compaction, differentiated snow density at deposition between precipitation and drifting snow, and a rewriting of the threshold friction velocity for snow erosion. Model results at high resolution (10 km) over Adelie 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. MAR resolves the local drifting-snow frequency and transport up the scale of the drifting-snow event and cap10 tures the resulting observed climate and SMB variability. This suggests that this model version can be used for continent-wide applications, and that the approach of drifting-snow physics as proposed in MAR can serve as a basis for implementation in earth system models.

regional scale also arise from the numerous interactions of drifting-snow particles with the atmosphere and the snow surface organised 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 (Gallée et al., 2013) and is mirrored through a high sensitivity of model results to parameter choices and significant discrepancies 60 between simulated and observed snow mass fluxes (Lenaerts et al., 2014;Amory et al., 2015;van Wessem et al., 2018).
The polar-oriented regional climate model MAR includes a drifting-snow scheme initially developed to improve the representation of the Antarctic SMB (Gallée et al., 2001). However, the drifting-snow scheme has only been used so far to study separately wind-driven ablation (Gallée, 1998;Gallée et al., 2001Gallée et al., , 2005 and individual drifting-snow events (Gallée et al., 2013;Amory et al., 2015) in coastal East Antarctica. Unpublished preliminary experiments with former physical parameterisations 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 decades long investigations of the SMB of the Greenland (e.g., Fettweis et al., 2017 and Antarctic (e.g., Kittel et al., 2018Kittel et al., , 2021Agosta et al., 2019) 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 70 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 Adelie Land, a katabatic wind region of East Antarctica which experiences drifting snow frequently (Amory, 2020a) and where the SMB exhibits a high spatial variability related to drifting-snow processes (Agosta et al., 2012). 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 75 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. MAR is a hydrostatic atmospheric model originally developed to simulate the climate over high-latitude regions. The atmospheric dynamics in MAR are described in Gallée and Schayes (1994). Cloud microphysical processes and resulting precipitation are simulated by solving conservation equations for specific humidity, cloud droplets and ice crystals, raindrops and snow particles (Gallée, 1995;Gallée et al., 2001). The radiative transfer through the atmosphere is adapted from Morcrette (2002), and cloud radiative properties are computed from the concentration of cloud droplets and cloud ice crystals. Turbulence is 85 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 (Gallée et al., 2015).
In this study MAR version 3.11 in its Antarctic setup (Agosta et al., 2019) is used with the updates described in (Kittel et al., 2021), simply referred to as MAR hereafter. The simulations are performed on a grid of 80 × 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 90 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 (Fretwell et al., 2013). 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 (Hersbach et al., 2020). The atmosphere is described on a stretched grid with 24 vertical terrain-following levels, of which 8 and 5 are respectively located 95 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 (van de Berg and Medley, 2016).

Snowpack model
The atmospheric part of MAR is coupled with the one-dimensional multi-layer surface model SISVAT (Soil-Ice-Snow-Vegetation-100 Atmosphere-Transfer) which handles energy and mass transfer between the surface and the atmospheric boundary layer (De Ridder and Gallée, 1998). SISVAT includes a representation of snow (Gallée and Duynkerke, 1997;Gallée et al., 2001) and ice (Lefebre et al., 2003) 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 105 to precipitation, erosion/deposition of snow, mechanical compaction, 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 110 to refine the discretisation of the snowpack when the number of layers is lower than 10. Maximum layer thicknesses of the 4 uppermost layers of the snowpack are also prescribed (0.02, 0.05, 0.1, and 0.3 m) to ensure a fine discretisation 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 initialised with snow grain shape parameters of fresh snow (see Sect. 3.3 for definition) and a density of 500 kg m −3 assuming a null liquid water content. The 115 initial surface snowpack temperature is set to the reanalysis near-surface air temperature and then discretised 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  3 Drifting-snow scheme 120 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 sketch ( Fig.   1) provides a general overview of the drifting-snow scheme.

Threshold friction velocity for initiation of snow erosion 125
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 −1 ) exceeds the threshold value for aerodynamic entrainment, i.e., the threshold friction velocity u * t (in m s −1 ) determined by the resistive gravitational and cohesive forces. Resistive forces depend on temperature (Schmidt, 1980) and metamorphism history (Gallée et al., 2001) 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 * t requires a detailed 130 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 parameterisations of u * t (e.g., Liston et al., 2007;Lenaerts et al., 2012a). The same approach is followed here; erosion in the model occurs when u * > u * t , where u * t is imposed by the uppermost snow layer density in which ρ s is the surface snow density (kg m −3 ), ρ i is the density of ice (920 kg m −3 ) and ρ 0 is the density of fresh snow (set to 300 kg m −3 ). The expression for u * t0 is retained from (Gallée et al., 2001) where i ER is an erodibility index describing the potential for snow erosion, d (dendricity) and s (sphericity) are the snow grain shape parameters, C D the drag coefficient for momentum, U the wind speed at the lowest prognostic level of the model (m s −1 ) and u * is classically obtained through integration of stability correction functions for momentum over the atmospheric 145 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 ER was defined as a function of surface snow characteristics (d, s and grain size). To reduce the number of sensitivity parameters, as in Lenaerts et al. (2012a) u * t is assumed independent on particle size and constant snow grain shape parameters are assigned (d = s = 0.5), implying an 150 erodibility index of 0.625 and a minimum u * t value of 0.3 m s −1 for ρ s of 300 kg m −3 . Additional criteria for snow erosion require that ρ s does not exceed ρ max = 450 kg m −3 and that no wet layer has formed at the top of the snowpack. The sensitivity to the formulation of u * t is quantified and discussed in Sect. 6.2.

Snow-transport modes
Particle motions in drifting snow are generally described through the two main transport modes consisting in saltation and 155 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 centimetres. 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 160 airborne snow mass. In the model the particle ratio in the saltation layer q salt (kg kg −1 ; mass of saltating snow particles per unit mass of atmosphere) is parameterised as a function of the excess of shear stress responsible for removal of snow particles from the surface following the expression of Pomeroy (1989) with e salt = 1/(3.25u * ) the saltation efficiency expressed as a dimensionless coefficient inversely proportional to the friction 165 velocity, g = 9.81 the gravitational acceleration (m s −2 ) and h salt =0.08436u 1.27 * the thickness of the saltation layer (m) as proposed by Pomeroy and Male (1992). Note that saltation is not explicitly resolved by the model and q salt , considered as constant throughout the saltation layer, only serves as a lower boundary condition for the suspension layer.
The formulation of q 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 (Aksamit and Pomeroy, 2018), numer-170 ical simulations suggest that steady-state saltation is achieved within an interval of a few seconds (Nemoto and Nishimura, 2004;Huang et al., 2016), well below the model time step of 60 s. With wind speeds at the surface typically reaching 10 m s −1 during drifting-snow occurrences (Fig. S3), 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.

175
Snow transport in turbulent suspension is computed by the tri-dimensional turbulence and advection schemes of MAR, enabling a discretisation 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 where q s * is the turbulent scale for snow particles (kg kg −1 ), ζ is the ratio of eddy diffusivities for suspended particles and momentum and q s is the snow particle ratio (kg kg −1 ) taken at the lowest model level. Because of fragmentation upon repeated 185 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 (Bintanja, 2000a). As MAR currently does not distinguish snow particles originating from the surface from those directly formed by its cloud microphysics, the factor ζ was introduced in the original formulation of Gallée et al. (2001) to enhance the upward turbulent particle transport and compensate for a likely overestimation of the settling velocity of drifting-snow particles (Gallée et al., 2005). Following 190 Bintanja (2000a), ζ is set to 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 distinction on 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 (Gallée et al., 2001). Sublimation of snow particles occurs along their residence into the atmosphere and is parameterised as a function of the snow particle ratio and undersaturation of air (Lin et al., 1983). In particular, the model assumes an exponential size distribution of suspended (cloud and eroded) snow particles (Gallée, 1995) n s = n 0 exp(−λ s D s ) 200 with n s the number of snow particles of diameter D s (10 −3 m) per unit volume, n 0 an empirical constant that corresponds to the intercept parameter of the size distribution (m −4 ), and λ s the (dimensionless) dispersion parameter where ρ p is the snow particle density (set to 100 kg m −3 ) and ρ a is the air density (kg m −3 ). Snow particles are considered as graupel-like snow of hexagonal type and the spectrally-averaged snow particle diameter D s is prescribed as a constant 205 following Locatelli and Hobbs (1974). 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 (Bintanja, 2001). Weakening of drifting snow in response to decreasing u * as turbulence declines is reflected through the 210 dependency of the surface snow turbulent flux u * q s * on the difference u 2 * − u 2 * t . Not distinguishing on the origin of particles despite differences in shape and size between cloud and eroded snow particles (Nishimura and Nemoto, 2005) 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 215 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 (Yamanouchi and Kawaguchi, 1984;Mahesh et al., 2003;Le Toumelin. et al., 2020). Representing the radiative contribution of drifting-snow clouds is achieved in the model 220 by including the snow particle ratio q s in the computation of cloud optical depth and emissivity (Gallée and Gorodetskaya, 2010).

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 225 of wind and saltation which break original crystal shapes and favour the formation of smaller, rounded snow grains (Sato et al., 2008), leading to enhanced sintering, more efficient mechanical packing and increased density (Vionnet et al., 2013;Sommer et al., 2018). Drifting-snow compaction, together with the exposure of denser snow or ice layers through erosion and/or sublimation, both naturally contribute 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 230 melting. Surface melting reduces or even inhibits the potential for erosion in summer by increasing water content, density and cohesion (Li and Pomeroy, 1997). 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 235 (60s) 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 240 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 (ρ 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 245 drifting snow (Huang et al., 2011). Drifting snow is deposited with a density ρ DR assumed to be that of the current surface layer, with the restriction that ρ DR does not exceed ρ max to account for maximum surface snow density values observed in Antarctica (Agosta et al., 2019) and enable for deposition of snow over more compacted snow and/or ice surfaces. The surface density ρ s is updated according to a relative contribution of both types of particles where f DR is the drifting-snow fraction varying between 0 and 1 and q 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 255 drifting-snow layers over Antarctica as retrieved from remote sensing techniques (Mahesh et al., 2003;Gossart et al., 2017;Palm et al., 2018). Snowfall conditions imply q s,z lim ∼ q s , low values of f DR and snow is deposited at the surface with a predominant contribution of ρ 0 . Conversely, drift conditions imply q s,z lim << q s , high values of f DR and ρ s tends towards The post-depositional increase in snow density through wind hardening is accounted for in the model by increasing the grid-average 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 parameterised according to a linear densification rate from the fresh snow value ρ 0 (assumed to be representative of snow that have been barely altered by post-depositional processes) to the prohibitive density value for snow erosion ρ max , i.e.
in which the characteristic time scale for drifting-snow compaction τ DR is set to 24 h. This value corresponds to the average duration of drifting-snow events reported in (Amory, 2020a) and is used here as the typical duration for exhaustion of erodible snow to be reached. The linear behavior of the densification rate follows the linear increase in surface snow density retrieved from measurements performed during a drifting-snow event in Adelie Land (Fig. S2). Further details on this experiment are provided in supplementary materials (Sect. S1). By fixing ρ 0 and parameterising u * t as an increasing function of ρ s (Eq. 1),

270
Eq. (11) 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 * t to endure through time and thus acts as a negative feedback for snow erosion.

Erosion
For each continental grid cell MAR calculates the actual snow mass eroded from the snowpack ER (kg m −2 ) during the current 275 time step according to the following chain of events ( Fig. 1): 1. The snow particle ratio q 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.
2. A potential maximum erosion ER max (kg m −2 ) is estimated from the surface turbulent flux of snow particles u * q s * 280 computed from step 7 at the previous time step: ER max =ρ a u * q s * dt.
3. Actual erosion ER (≥ 0) is calculated from removal of snow from the surface snowpack layer until ER max is reached or the layer has been entirely eroded.
5. The drift fraction is obtained from Eq. (10). Snow is deposited at the surface and surface density is adjusted according 285 to Eq. (9).
6. The threshold friction velocity u * t and the saltation particle ratio q salt are deduced from Eqs. (1) and (5).
7. The surface turbulent flux of snow particles u * q s * is computed from Eq. (6).
8. The contribution of erosion is added to q s , which is then transmitted to the turbulence, cloud microphysical and radiative schemes to compute advection and interactions of airborne snow with the atmosphere.

Surface roughness
Drifting snow is responsible for the development of surface microrelief, whose spatial arrangement combined to the orientation of the wind determines the roughness length for momentum z 0 . Because of wind-driven reshaping of the snow topography and the diversity in surface types, z 0 varies by several orders of magnitude with time and space across the Antarctic continent . The sensitivity of MAR to the parameterisation of surface roughness has been discussed in Amory roughness elements encountered in desert-like environments (Marticorena and Bergametti, 1995). 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 (Andreas and Claffey, 1995;Amory et al., 2016). Sensitivity experiments revealed that the drag partition scheme was responsible for a strong inhibiting effect on snow erosion beyond the observed magnitude of 310 the negative feedback mechanism, possibly as it does not account for the dynamical and erodible nature of snow microrelief (Amory et al., 2015). Consequently, it has been disabled in the current model version.

Field area, observation data and evaluation methods
The near-surface climate in coastal Adelie 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 315 the year (Amory, 2020a). 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 (Genthon et al., 2007;Favier et al., 2011). In a fairly narrow transition, net accumulation is observed a few kilometers inland despite significant drifting snow (Barral et al., 2014;Amory et al., 2017), with annual SMB values displaying a high, kilometre-scale variability as a result of wind redistribution (Agosta et al., 2012). The performance of MAR in reproducing the drifting-snow climate of Adelie 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.   1, Table 1). 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 325 water, and has been converted to be expressed with respect to ice in subfreezing conditions following Goff and Gratch (1945) 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.

320
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 in- 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., Mann et al., 2000;Nishimura and Nemoto, 2005). As the snow particle ratio q 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 µ OBS (kg m −2 s −1 ) 350 calculated through where µ i is the observed snow mass flux integrated over the exposed length h i of the corresponding 2G-FlowCapt™ sensor (Amory et al., 2015).
Modelled SMB is compared with observations obtained from annual measurements of 91 snow stakes distributed every ∼1.5 is a well-known feature of the study area supported by modelling efforts and observations (e.g., Parish and Wendler, 1991;Wendler et al., 1993).
Snow at the surface densifies with drifting snow (Eq. 11), causing the spatial distribution of ρ s to be linked with the variability in near-surface wind speed (Fig. 1b).  (Fig. 3a). 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 −1 and a rmse of 3.2 m s −1 (Fig. 3b). 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 380 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 −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 parameterisation of z 0 , locally still yielding too high values, while at D17 Fig. 4 illustrates that modelled z 0 values are closer to observations.   . 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.
The saturation vapour pressure of air is a strongly dependent function of temperature, and the relative humidity determines 395 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. 3c,d) is well represented at both stations (positive bias < 1°C and centered rmse < 2°C). As drifting-snow sublimation is a determining contributor to the atmospheric humidity  from several h to a few days at most (Amory, 2020a). To evaluate the model results at higher temporal resolution closer to characteristic time scales of drifting snow, half-hourly data are used to compute the probability of detection (P OD) and false alarm ratio (F AR) 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 (Rousseau, 1980) and already used to evaluate modelled driftingsnow occurrence in an alpine context (Vionnet et al., 2013), is also calculated. RI also takes into account the number of 430 occurrences for which the absence of drifting-snow is correctly simulated, d, such that 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 2), meaning the model ability to predict drifting-snow occurrences 435 at the monthly scale (Fig. 5b,c) arises from a reasonable reproduction of their actual timing and duration at the half-hourly resolution. MAR shows better results (higher P OD and RI) at D17 than at D47, but also simulates more unobserved occurrences (higher F AR) that compensate for missed occurrences in the calculation of monthly frequency values.

Drifting-snow transport
Qualitative evaluation of MAR against monthly cumulative near-surface drifting-snow mass fluxes at D47 and D17 (Fig. 6b,c) 440 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 7 compares observed and simulated snow transport during each observed drifting-snow event. As in Amory (2020a), a drifting-snow event is defined as 445 a period over which the observed snow mass flux is above the detection threshold of 10 −3 kg m −2 s −1 for a minimum duration of 4 h. Discrepancy between model and observations is larger for (shorter) events of lower magnitude (< 10 3 kg m −2 s −1 ) and reduces for longer, more important events that contribute predominantly to the local drifting-snow transport (Amory, 2020a). Table 3. Observed and simulated total drifting-snow transport (10 6 kg m −2 ) from 0 to 2 metres above surface at D47 and D17. A distinction is made between total snow transport computed from observed (OBS) and simulated (MAR obs ) near-surface drifting-snow mass fluxes cumulated over each drifting-snow event identified in the database and over each simulated drifting-snow event (MARsim).  Table 3 compares the total snow transport during observed drifting-snow events (OBS) to that estimated by the model during these observed events (MAR obs ) and during every simulated events (MAR sim ). Referring only to observed events, MAR 450 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. 5b), 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. 7a). 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 sim ∼ MAR obs > OBS). At D17, while 455 the simulated drifting-snow frequency (Fig. 5c) and total drifting-snow transport during observed events are in closer agreement with the observations, simulated but unobserved drifting snow is more significant (Table 3) and result in overestimation of total drifting-snow transport (MAR sim > MAR obs ∼ OBS).
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 460 cumulative near-surface transport during drifting-snow events are shown at both measurement sites for observed 2-m windspeed bins of 1 m s −1 in Fig. 8 for the most represented wind-speed categories (8 to 17 m s −1 ). Since wind speed and driftingsnow 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 8 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 465 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 (Trouvilliez et al., 2015), its ability to estimate drifting-snow mass fluxes remains to be assessed in Antarctic conditions (Amory, 2020a).
The mass transported in drifting snow has been shown to correlate with wind speed in a power-law fashion (e.g., Budd, 470 1966;Radok, 1977;Mann et al., 2000;Amory, 2020a). Figure 9, 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 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. 6a). The dispersion around a given wind speed value (

Surface mass balance 480
In the current version of MAR, cloud, eroded, and deposited snow particles are all included in the snow particle ratio q 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 which reflects the relative local proportion of eroded to deposited snow mass. Figure 10 compares the simulated and observed mean annual SMB and shows the elevation and i E/D along the stake transect whose location is marked with black dots on Fig. 11. MAR represents the general variability in SMB with a strong increase 490 over the first tens of kilometers from the coast and less variability further inland. The variability in i E/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 ∼67 kg kg m −2 y −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 E/D along the transect. The lowest, negative simulated 495 SMB value is not found at the coast but a few grid cells inland where i E/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 E/D displays maximum values (Fig. 11). 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 Adelie Land (Genthon et al., 2007;Favier et al., 2011).

500
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. 10). It also shows that the modelled SMB decrease a few kilometers 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 and knowing that different reanalysis (ERA5 vs ERA-Interim) were used as forcing, MAR displays enhanced capabilities in 515 simulating the timing, duration and magnitude of drifting snow (Fig. 12a) relative to the previous model version and setup (MARv2) used in Amory et al. (2015). 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 driftingsnow scheme. Improved drifting snow in MAR also leads to increased 2-m relative humidity as a result of enhanced atmospheric sublimation (Fig. 12b). Periods of near-saturated conditions, whose occurrence corresponds with that of drifting snow, are better

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 ρ 0 and characteristic time scale for drifting-snow compaction τ DR over a range of otherwise plausible values. Sensitivity to the dis-525 abling of drifting-snow compaction and to the parameterisation of the threshold friction velocity for snow erosion u * t and the saltating particle ratio q salt was also explored. Simulated drifting-snow frequency and transport is 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 summarised in Table 4.
The model exhibits a high sensitivity to ρ 0 . This parameter determines the density at which fresh snow is deposited but also 530 intervenes in the definition of u * t (Eq. 1) and drifting-snow compaction rate (Eq. 11). To conserve the physical consistency of the drifting-snow scheme, changes in ρ 0 must be applied to the whole set of parameterisations and necessarily lead to modifications that result from the combined sensitivities to each of the three parameterisations. This requires notably the adoption of values for ρ 0 that are probably above typical values for fresh snow, as u * t becomes a more restrictive criterion for erosion with lower values for ρ 0 (Fig. 13), and consequently that already partially integrate the effect of post-depositional (blue) for the case study in January 2011 at D47 discussed in Amory et al. (2015). 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. 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 5 kg m −2 y −1 during year 2010 at D47 and 2013 at D17. Values in brackets refer to the relative difference in % compared to observations. L07 and V12 denotes adaptation of the parameterisation of u * t according to Liston et al. (2007) andVionnet et al. (2012) respectively (see Appendix A). B00 denotes adaptation of the parameterisation of q salt according to Bintanja (2000b) depending on dominant local effects. Reducing (increasing) ρ 0 to 250 (350) kg m −3 at D17 leads to lower (higher) snow density values at deposition by snowfall but also increases (decreases) u * 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 ρ s is locally overcompensated by a corresponding strong increase in u * t . The model would become less and less sensitive to ρ 0 540 values below 250 kg m −3 because u * t would more rapidly reach high prohibitive values for snow erosion. The reverse situation is depicted at D47, where the model simulates lower ρ s (Fig. 2b). Lowering ρ 0 to 250 kg m −3 induces a local decrease in ρ s that improves drifting-snow frequency but produces drifting-snow transport beyond observed values (> 50 %) despite higher u * t . Overestimated and similar drifting-snow features are obtained at both locations with ρ 0 = 350 kg m −3 . Discussing the influence of ρ 0 beyond 350 kg m −3 would implicitly imply deposition of fresh snow in a well-advanced state of compaction 545 and is therefore not considered Changes in τ DR alter the timing and duration of drifting-snow events but appear to exert a moderate influence on modelled drifting snow. Halving and doubling τ DR to 12 h and 48 h respectively decreases and increases local drifting-snow frequency and transport with a relatively lower magnitude (< 30 %) than the tested changes in ρ 0 . The value of τ DR = 24 h is of local significance since it has been scaled on the observed median duration of drifting-snow events in Adelie Land. Different time 550 scales for drifting-snow compaction could be expected under different environmental conditions, particularly for different wind speed regimes. This would suggest that τ 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 τ DR to infinity and the density of deposited drifting snow ρ DR to the standard ρ 0 value of 300 kg m −3 produces almost permanent drifting snow and leads to the strongest overestimation (> 100 % at D47) 555 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 * t following the parameterisations of Liston et al. (2007) and Vionnet et al. (2012), hereafter referred to as L07 and V12. Appendix A gives numerical 560 details. The resemblance in the evolution of u * t with ρ s (Fig. 13) explains the similar model results obtained at both stations with L07 and V12. However, both parameterisations do not depend on the same set of variables. L07 only depends on ρ s and enables investigation of the influence of ρ 0 without inherently altering u * t nor τ DR . Overall, L07 is a stricter criterion for snow erosion than the control parameterisation to the extent that it cannot be offset by prescribing low ρ 0 values nor by increasing τ DR within the range of tested values (not shown). V12 corresponds to a weighting of the erodibility index i ER by ρ s and 565 remains dependent on ρ 0 . Like L07, V12 yields a general increase in u * t that becomes very restrictive and that a change in τ DR can hardly compensate for without increasing it substantially nor by increasing ρ 0 (not shown). The resulting inhibition at D17 is also quite close to that obtained with ρ 0 = 250 kg m −3 , where ρ s is high enough for u * t to drive the model response to changes in ρ 0 .
The model is virtually insensitive to changes in the formulation of the saltating particle ratio q salt . Equation (5) gives a 570 maximal value of q salt for intermediate values of the friction velocity u * . This is a potential limitation of this formulation that has been identified by Bintanja (2000a). Bintanja (2000b) tested another formulation, hereafter referred to as B00, in which the efficiency of saltation is assumed to be constant (e salt = 0.535), yielding a relation for q salt that increases monotonically with u * . Figure S4 illustrates that a significant sensitivity to the parameterization of q salt could be expected, especially for u * well above 0.6 m s −1 from which the control formulation and B00 start to diverge from each other. However, the comparison 575 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 salt . Table 4 shows that the best agreement between model and observations is obtained with the control setup. However, this 580 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 ρ 0 and the parameterisation of u * 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 585 consistency between erosion and accumulation rates. 3). The solid blue line highlights the current control parameterisation of u * t in MAR (ρ0 = 300 kg m −3 ).

Conclusion
A dataset of drifting-snow and SMB observations collected in Adelie 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 0 as a varying function of air 590 temperature and disabling the drag partition scheme and its inhibiting effect on snow erosion, ii) rewriting the parameterisation of the threshold friction velocity for snow erosion u * t as a function of surface snow density only, iii) differentiating snow density at deposition between precipitation and drifting snow, iv) including a simple parameterisation for drifting-snow compaction of the uppermost firn 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-595 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.

600
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 remobilised from the ground from precipitating particles to explicitly account for different particle properties such as 605 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 610 atmospheric levels where drifting-snow transport and sublimation can supposedly be of substantial significance in terms of SMB (e.g., Amory and Kittel, 2019;Le Toumelin. et al., 2020) 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 (Markus et al., 2017) at an unprecedented level of spatial and temporal resolution would be a great candidate for such an exercise.

615
Code and data availability. The meteorological and drifting-snow data  can be downloaded at https://zenodo.org/record/3630497 (last access: 30 March 2021). The surface mass balance data are available at pp.ige-grenoble.fr/pageperso/faviervi/stakeline.php. This service is provided by the National Observation Service GLACIOCLIM (https://glacioclim.osug.fr/?lang=en) 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 webpage hosted by the OSUG website. However, we have also archived a specific version of the 620 dataset used in this study together with the MAR source code and outputs for replication of this study at https://zenodo.org/record/4314872 (Amory, 2020b Changes in the threshold friction velocity for snow erosion u * t was done following the parameterisations of Liston et al. (2007) and Vionnet et al. (2012). In Liston et al. (2007), u * t is parameterised a function of surface snow density ρ s only u * t =    0.005exp(0.013ρ s ) 300 < ρ s ≤ 450 0.1exp(0.003ρ s ) ρ s ≤ 300 In Vionnet et al. (2012) the definition of the erodibility index i ER is completed with a term depending on the surface snow density F(ρ s ) to account for high polar surface snow densities (Libois et al., 2014)