A Drifting and Blowing Snow Scheme in the Weather Research and Forecasting Model

Wind‐driven redistribution of surface snow is one of the key factors leading to heterogeneous accumulation of snow at small scales. Understanding the processes that lead to this heterogeneous accumulation is, therefore, of great importance to many glaciological and hydrological questions. High‐quality information on the wind field is necessary to realistically represent drifting snow. Here, we introduce a novel, intermediate‐complexity drifting and blowing snow module for the Weather Research and Forecasting (WRF) model that integrates seamlessly into the standard WRF infrastructure. The module also accounts for snow particle sublimation and considers the thermodynamic feedback on the atmospheric fields. In an idealized model environment the module was tested for physical consistency. Sensitivity experiments showed that drifting snow sublimation has on the one hand local effects on the erosion and deposition patterns and on the other hand non‐local effects on the larger‐scale atmosphere. The simple and computationally efficient implementation allows this module to be used for small‐scale and large‐scale simulations in polar and glaciated regions.


Introduction
The spatial and temporal variability of snow cover in complex mountain areas is an important component of the water budget of hydrological catchments and one of the key drivers for the mass balance of alpine glaciers.Larger-scale orographic enhancement of precipitation (Houze, 2012) as well as small-scale micro-physical processes and interactions between the near-surface flow field and hydrometeors (often referred to as preferential deposition) lead to heterogeneities in the initial accumulation of precipitation (e.g., Gerber et al., 2019;Lehning et al., 2008;Mott et al., 2014;Vionnet et al., 2017;Zängl, 2007Zängl, , 2008)).Redistribution of snow by the wind during and after the precipitation event further alters the small-scale variability in snow accumulation.Mott et al. (2011) demonstrated that variability in the snow accumulation at the ridge scale, in the order of tens of meters, is dominated by drifting snow processes and preferential deposition.Topographic features such as ridges and gaps shape the small-scale heterogeneity by locally accelerating the flow so that shear stresses and turbulence kinetic energy become strong enough to effectively remove snow from the surface (Mott et al., 2018).Turbulent motions transport the eroded snow away from the surface so that the particles can be advected over long distances with the background flow.During transport the particles are subject to gravity-driven subsidence and will accumulate in wind-protected areas.Depending on the atmospheric conditions, drifting snow particles often lose mass by sublimation, which in consequence modifies the vertical structure of temperature and moisture.
Several studies showed that redistribution of snow is essential to understand the spatial structure of local snow accumulation on timescales ranging from single events to entire seasons (Bernhardt et al., 2012;Musselman et al., 2015; J. W. Pomeroy & Gray, 1995;J. W. Pomeroy et al., 1998;Schirmer et al., 2011;Vionnet et al., 2021).Different studies implied that their calculations of glacier mass balance can be improved by accounting for drifting snow processes in the simulations (Abrahim et al., 2023;Mortezapour et al., 2020;Pradhananga & Pomeroy, 2022;Temme et al., 2023;Terleth et al., 2023).The importance of drifting snow sublimation was found to be highly variable depending on the environment and the flow conditions (e.g., Groot Zwaaftink et al., 2011;Groot Zwaaftink et al., 2013;Le Toumelin et al., 2021;Lenaerts & Van den Broeke, 2012;Lenaerts et al., 2010;MacDonald et al., 2010;Sauter et al., 2013;Strasser et al., 2008).Apart from alpine contexts, polar ice bodies have been shown to be greatly influenced by drifting snow on the one hand due to the redistribution and drifting snow sublimation (e.g., Amory et al., 2021;Lenaerts et al., 2010), but also by alternating the radiative budget (Le Toumelin et al., 2021).Nevertheless, both the measurement and modeling of drifting snow remain challenging due to the very small-scale nature and the complex interplay of many simultaneously acting processes.
In our study, however, we want to introduce a new, intermediate-complexity framework to simulate drifting snow.Our approach makes use of the already existing capabilities of the Weather Research and Forecasting (WRF) model including a standard land-surface model (Noah-MP) for the treatment of snow.Apart from that, the tracer module is employed for particle transport, hence all WRF-internal schemes and options for advection and turbulent diffusion can be used.With this intermediate-complexity design we are able to make use of high-quality wind fields and a reasonable representation of the snow pack while keeping the computational costs rather low, which will support our future goals of investigations on climatic timescales.This paper describes the physical and numerical fundamentals of the scheme and analyses the skill and physical consistency of this intermediate-complexity snow drift scheme.In our paper, we will use the term drifting snow for all wind-driven transport of snow without distinguishing additionally between drifting and blowing snow.
The remainder of this paper will be structured as follows: Section 2 gives an overview on the snow drift scheme and its implementation in WRF.Section 3 introduces the setup for the idealized test environment, of which the results are presented in Section 4. Finally, Section 4.5 gives an outlook for future real-world applications where we apply the module to a drifting snow event in the Alps and where we can verify the modeled results by comparing them to snow depth changes observed by Terrestrial Laser Scanning.

Snow Drift Scheme
The implemented snow drift scheme is based on the work of Sauter et al. (2013).
We employ the widely used approach of dividing the process of drifting snow into a saltation layer and snow particles in suspension (Mott et al., 2018).The saltation layer is fully parameterized and acts as a lower boundary condition for the flux of snow into the atmosphere, while snow particles in suspension are transported by the mean wind and turbulent diffusion and are subject to gravity-driven subsidence as well as to sublimation.This scheme therefore treats drifting of snow as continuous phases that interact only with the background flow while neglecting particle interactions (Bintanja, 2000;Durand et al., 2005;Gauer, 2001;Liston & Sturm, 1998;Naaim et al., 1998;J. Pomeroy & Gray, 1990;Sauter et al., 2013;Schneiderbauer et al., 2008).According to the continuum equation for mass conservation, the local rate of snow concentration change ∂ϕ s /∂t [kg m 3 s 1 ] is given by . (1) Here, (II) describes the advection by the mean wind with x i the cartesian coordinates and u i the cartesian components of the wind vector.The term (III) describes the turbulent flux divergence, (IV) the downward flux of snow particles due to gravity and (V) the concentration loss due to drifting snow sublimation.Advection and turbulent diffusion are handled by WRF-internal modules for the transport of passive tracers, while we parameterize terms (IV) and (V) as described in the following.The fallout velocity at height z [m] is following Vionnet et al. (2014) with the two parameters and Here ν air [m 2 s 1 ] represents the viscosity of air, ρ ice and ρ air [kg m 3 ] the density of ice and air, respectively, and g [m s 2 ] the acceleration due to gravity.The height-dependent mean particle radius r(z) [m] is calculated analogous to Sauter et al. (2013) with where r 0 [m] is the particle radius at ground level following Gordon et al. (2010) with: r 0 = 0.5( 7.8 ⋅ 10 6 u ⋆ 0.036 and u ⋆ [m s 1 ] the friction velocity.Optionally, V(z) and r 0 can be set to constant values in the namelist settings.

Journal of Advances in Modeling Earth Systems
10.1029/2023MS004007 SAIGGER ET AL.
Term (V) of Equation 1 accounts for the mass loss of suspended snow by sublimation and is calculated based on the formulation of Thorpe and Mason (1966) that has widely been used in earlier studies (e.g., Bintanja, 2000;Déry & Yau, 1999;Gallée et al., 2001;Sauter et al., 2013;Sharma et al., 2023;Vionnet et al., 2014).Sublimation from particles in the saltation layer is neglected similar to earlier modeling approaches.The sublimation-loss rate of suspended snow is approximated by ψ s ϕ s [kg m 3 s 1 ], where ψ s [s 1 ] is the sublimation-loss rate coefficient describing the change of snow particle mass due to heat exchange and ventilation effects while neglecting radiation effects: Here σ is the water vapor deficit with respect to ice, L s the latent heat of sublimation (2.838 • 10 6 J kg 1 ), K the thermal conductivity of air (0.024 W m 1 K 1 ), R v the gas constant for water vapor (461.5 J kg 1 K 1 ), D the molecular diffusivity of water vapor in air (2.25 • 10 5 m 2 s 1 ), T air the air temperature, and e i the saturation vapor pressure with respect to ice.The Nusselt number and Sherwood number Nu and Sh describe the transfer of heat and mass at the particle surface with where Re(z) is the particle Reynolds number: The model therefore reflects the ventilation effect of the falling snow particles, ventilation due to turbulence is neglected.The scheme considers the effect of sublimation on the vertical profiles of temperature and humidity.
For that, the local air temperature is reduced by the heat required for sublimation and the local specific humidity is increased accordingly.This feedback self-limits the sublimation process since its intensity depends on the saturation deficit and the air temperature of the environment (Sauter et al., 2013).
As described above, the drifting snow is divided into a saltation layer and the suspended snow particles.The snow mass concentration within the saltation layer is primarily gained by the aerodynamic entrainment of snow particles from the underlying snowpack.Once the surface shear stress exceeds the inertia and the cohesive bonds of the snow particles, transport takes place.Therefore, the vertical erosional mass flux from the surface into the saltation layer q e [kg m 2 s 1 ] is assumed to be proportional to the excess surface shear stress (Schmidt, 1980) with Here, ρ air [kg m 3 ] is the air density and u th [m s 1 ] the friction threshold velocity.The heuristic parameter e salt [ ] controls the efficiency of the erosion process.The threshold friction velocity u th is proportional to the snow density ρ snow [kg m 3 ] following Walter et al. (2004): Similar to Sauter et al. (2013) the drag exerted by the particles on the flow field is accounted for by limiting the saltation-layer snow concentration ϕ salt to a maximum value ϕ salt,max (J.Pomeroy & Gray, 1990): Additionally, a drag-corrected friction velocity u ⋆c following Naaim et al. (1998) is introduced that reduces the effective friction velocity to u th in the case of a fully saturated saltation layer: Replacing u ⋆ by u ⋆,c in Equation 10 results in the drag-corrected erosion flux q e,c q e,c = e salt ρ air Following J. Pomeroy and Gray (1990) the height of the saltation layer h salt is described as With the snow mass added to the saltation layer by q e,c and h salt the snow concentration in the saltation layer is updated.The saltation layer then serves as a reservoir for the flux of snow into the atmosphere f salt [kg m 2 s 1 ].Similar to Gallée et al. (2001) f salt is described as a bulk flux: where U [m s 1 ] and ϕ s,1 [kg m 3 ] are the horizontal wind speed and the snow particle concentration at the lowest model level of WRF, respectively.C D is the bulk transfer coefficient for momentum as calculated within the surface layer scheme.
As the friction velocity drops below the threshold, deposition takes place.Similar to Beyers et al. (2004) the deposition flux q d [kg m 2 s 1 ] corresponds to the downward flux at the lowest model level z 1 and the modified shear stress ratio

Implementation of the Snow Drift Scheme in the WRF Model
The implementation of the snow drift scheme consists of three basic modules, which include a module for erosion (module_snower.F), sublimation (module_snowsubl.F) and deposition (module_snowset.F) of snow particles as well as an auxiliary module for the calculation of the particle radius and the fall velocity (mod-ule_snowvel.F).The four modules are called by a high-level driver (module_snowdriver.F).The scheme is executed each time step between the land surface model (LSM) and the planetary boundary layer (PBL) scheme.The transport of snow particles is implemented as a passive tracer.This allows the user to choose between all available WRF-internal schemes and options for advection and turbulent diffusion.For real-case applications, the scheme requires the use of the LSM model Noah-MP (Niu et al., 2011).Noah-MP employs a snow model with up to three layers, that allows to represent multiple processes within the snow pack.Similar to other intermediate-complexity models, Noah-MP has been shown to substantially improve snow pack simulations both on a daily and seasonal scale compared to simpler bulk approaches (Arduini et al., 2019;Niu et al., 2011) while avoiding the relatively high computational costs of high-complexity models (Sharma et al., 2023).The snow drift scheme can be activated via the WRF namelist.The modules for erosion, sublimation and deposition can be activated separately.Basic parameters for calculating particle size, fallout velocity and the saltation coefficient can also be set (Table 1).Section 4.4 summarizes the sensitivity of simulations through setting different parameter values.
Before integrating the snow particle transport equation, the scheme passes through the following steps: 1. Calculation of the erosion flux for each snow-covered grid cell as a function of the maximum snow particle concentration in the saltation layer.The erosion flux decreases the snow on the ground and leads to the increase of the snow particle concentration in the saltation layer.(module_snower.F) 2. Then the mass flow between the saltation layer and the lowest level of the suspension layer is calculated.
(module_snower.F) 3. Now the particle radius and the fallout velocity are determined.These quantities are then passed on to the sublimation module.(module_snowvel.F) 4. The snow drift sublimation is now calculated in each grid cell.The energy and mass fluxes due to sublimation then modify the specific humidity and air temperature within the grid cells (if feedback has been activated in the namelist).(module_snowsubl.F) 5.In the last step, the deposition flux is determined, which is calculated on the basis of the particle concentration and fallout velocity.The deposition flux in the lowest model layer leads to the accumulation of snow particles and thus increases the snow on the ground.(module_snowset.F)

Model Setup
The WRF model has been described in numerous previous papers, as it is "the world's most-used atmospheric model" (Powers et al., 2017).In brief, it is a fully non-hydrostatic numerical model and includes a large suite of physics options (e.g., radiation, boundary layers, shallow convection) as well as the possibility to conduct idealized or realistic runs.To show the capability of the scheme, we perform idealized simulations with the WRF version 4.1.2(Section 2).The next paragraphs describe the settings for the control run (CTRL) and the sensitivity experiments in an idealized setup.An outlook toward real-case applications is presented in Section 4.5.
The model setup is inspired by the classical studies which perform idealized atmospheric simulations for downslope windstorms (Nappo, 2013) or orographic precipitation (e.g., Jiang, 2003;Kirshbaum & Durran, 2004;Mölg et al., 2009).Accordingly, we assume a Gaussian-shaped hill with a height h of where h m is the height of the mountain peak and a is the mountain half-width, for which we chose 500 m and 4 km, respectively.The domain size is 500 grid cells in the x-direction and 300 cells in the y-direction, with a grid spacing of 200 m in both horizontal directions; this eventually produces a domain 100 km (W-E) by 60 km (N-S) wide (Figure 1a).The domain is larger in W-E direction as the meteorological wind input is a zonal flow (see below) and, thus, the flow has ample space to develop in front and after the obstacle.In any case is the domain size well beyond 10 • a in the horizontal and even larger than 20 • a in the W-E dimension (e.g., Jiang, 2003).In the  vertical, 73 terrain-following levels reach up to an elevation of 10 km above the surface, and the smallest layer depth (highest resolution) of 10 m at the surface increases to a maximum of 200 m about 4 km above the surface using grid stretching (e.g., Wang & Kirshbaum, 2015).Rayleigh damping is activated for the upper 5 km of the domain to prevent gravity waves from artificial reflection.
The simulations are initialized horizontally homogeneously with profiles of potential temperature, specific humidity and the zonal and meridional wind component (Figures 1b-1d).The boundary conditions employed are periodic in y-direction and open in x-direction.The initial profiles are defined by a dry buoyancy (Brunt-Väisälä) frequency (N d ) of 0.015 s 1 , and a velocity U of 7.5 m s 1 of the zonal wind (meridional wind is zero) which results in a non-dimensional mountain height ĥ = N d h m U = 1 (Lin, 2007).In order to avoid artifacts at the domain boundaries a logarithmic wind profile is employed for the lowest 100 m of the initial profile.The entire initial atmospheric profile has a relative humidty of 80%.In consideration of representing an alpine environment at an altitude of about 2,000 m, the surface potential temperature is set to 290 K and the surface pressure to 800 hPa at the model initialization.With regard to the snow drift process, CTRL is initialized with a snow density of 100 kg m 3 over the entire domain, representing freshly fallen snow, and a snow depth of 1 m which ensures that at no point in the domain the entire snow pack gets removed.
Each model run is integrated 6 hr forward with a fixed time step of 0.9 s.After 3 hr of spin-up time a quasistationary state is reached, therefore we use only the output from hours 3 to 6 for the analysis.Due to the idealized nature of the runs and in order to attribute snow depth changes to snow drift as clearly as possible, most WRF-included physics options are deactivated.This is standard in idealized simulations, and the question of which options to keep depends on the specific processes that shall be examined by the model runs.In our case, the MM5 surface layer scheme (Fairall et al., 2003) is employed.We use the MYNN level 2.5 planetary boundary layer scheme (Nakanishi & Niino, 2006), however, we keep in mind that with our horizontal grid spacing of Δx = 200 m our simulations are within the "gray zone" of turbulence (Wyngaard, 2004).The land surface model is deactivated, which means that after initialization only erosion and deposition can alternate the snow on the ground.Other processes like melt, compaction or metamorphism therefore are not taking place in this case.In the snow drift scheme, CTRL has the sublimation of drifting particles and the associated feedback on the temperature and humidity fields activated.The saltation coefficient is set to 5 • 10 4 (Sauter et al., 2013) and the particle radius at the surface results from Equation 5, while the particle fall speed is calculated with Equation 2. In order to get a deeper understanding for the behavior of the model, several sensitivity experiments were conducted as will be discussed in Sections 4.3 and 4.4.

Mesoscale Flow Conditions
In this first subsection, we evaluate whether the modeled mesoscale flow features agree with the known physical principles of earlier studies (overviews e.g. in Durran, 1990;Jackson et al., 2013;Lin, 2007;Nappo, 2013;Smith, 1979) and whether these flow features favor drifting snow.After the spin-up, a steady-state wind field establishes with an accelerated near-surface flow in the lee of the summit, and gravity waves propagating both vertically and horizontally (Figure 2).The main accelerated flow separates from the surface about 2 km downstream of the summit.A slight blocking upstream of the hill and a dominant zone of reversed flow downstream of the hill are present, as well as local wind speed maxima downstream of the northern and southern flanks of the hill.With the acceleration of the flow in the summit region the shear stress on the surface exceeds the threshold value of 0.23 m s 1 , which is u th for ρ snow = 100 kg m 3 based on Equation 11.Thus, we can expect snow to be eroded in this region.
The flow features described here are in line with what we can expect based on the results of the aforementioned studies.The simulations are initialized with vertically constant profiles of wind and static stability, which combine to a non-dimensional mountain height of ĥ = 1.The horizontal wave number of the terrain k (k ≈ 1 a , Lin, 2007) is smaller than the Scorer parameter l ≈ N d U , thus allowing for vertically propagating gravity waves.At ĥ = 1 the flow is at the transition of establishing non-linear features such as steepening of waves, down-slope wind storms, flow blocking and reversed flows (Lin, 2007).With the accelerated flow over the summit, and the strong flow at the flanks and downstream of the flanks, the wind field shows features of both a flow-over and a flow-around regime, marking the transition between the two regimes around ĥ = 1 (Jackson et al., 2013).
In summary, the inclusion of the snow drift module to our simulation reproduces the main expected flow features and provides an environment that allows for drifting snow.

Snow Drift Patterns
The atmospheric fields govern the spatio-temporal patterns of snow redistribution.Here we will revisit this dependency, as it provides fundamental insights into the inherent behavior of the snow drift scheme.The sensitivity of the snow drift parameters will then be explained in more detail in Section 4.4.We start with the discussion on the erosion flux and then further address both the deposition and sublimation fluxes, followed by a short comparison to reported real-world drifting snow concentrations.
Based on Equation 10 we expect the erosion flux to be proportional to the excess of surface shear stress.This expectation is confirmed in Figure 3a, where erosion only takes place in regions where the friction velocity exceeds 0.23 m s 1 (dashed contour line, 0.23 m s 1 is the threshold friction velocity for a snow density of 100 kg m 3 based on Equation 11).Apart from that, the erosion follows the structure of the near-surface wind speed as seen in Figure 2 with the maximum erosion (about 2 kg m 3 in 3 hr) in the region of highest wind speeds directly in the lee of the summit.In the same way, the snow concentration in the saltation layer is dominated by the surface shear stress (Figure 3b) with maximum concentrations of about 5 • 10 3 kg m 3 in the area of highest wind speeds.In contrast to that, more processes play a role when considering the snow concentration in suspension.Snow particles are eroded from the surface and, thus, have their maximum concentration in the lowest layer.However, turbulent motions also mix the particles within a layer of approximately 100 m (Figure 3d).Apart from that, snow particles in suspension are affected by sublimation, subsidence and advection with the mean flow (see  1).Once steady-state is reached in the simulation, these processes balance each other, resulting in the suspended snow concentration as shown in Figure 3c.Here, the snow concentration does not have its maximum in the region of strongest erosion, but gradually increases in flow direction reaching its maximum of about 1.5 • 10 3 kg m 3 at the end of the erosion zone.While the snow particles are confined to a shallow layer in the erosion zone where the vertical transport is dominated by turbulence, they are advected further away from the surface in the region of flow separation (Figure 3d).
Net deposition can only occur when the friction velocity drops below the friction velocity threshold, provided that snow particles are present in the atmosphere and the downward flux due to gravity exceeds the vertical turbulent mixing.In our case, these conditions are met directly downstream of the zone of flow separation, resulting in a narrow strip where snow is deposited (positive values in Figure 3a).Comparing the total amount of accumulated erosion to the deposition, only about 3% of the eroded mass is deposited, while the rest is removed by sublimation and suspended by the prevailing flow.
Depending on the atmospheric conditions, snow particles lose mass due to sublimation.In our implementation, several factors determine the amount of sublimation: the available snow concentration, the particle size, the particle ventilation, as well as the air temperature and relative humidity.The influence of the particle concentration and relative humidity can be seen in the spatial distribution of accumulated sublimation in Figure 4.The highest sublimation amounts occur in the region of maximum particle concentration in the lee of the summit.Presumably also the warming and thus drying of the air due to subsidence and entrainment of air from above will enhance the sublimation here.On the other hand, the high relative humidity in the summit region inhibits sublimation there compared to regions at the northern and southern flanks of the hill where despite the lower snow concentration more sublimation is taking place (Figures 3c and 4a).In the experiments for parameter sensitivity in Section 4.4 we will explore further influences of particle size and ventilation speed.
While keeping in mind the idealized nature of the simulation presented above, a brief discussion of our results in the context of previous field measurements seems appropriate.With a simulated snow concentration in the saltation layer on the order of 5 • 10 3 kg m 3 (Figure 3b) our results roughly agree with concentrations reported by J. Pomeroy and Male (1992) or Doorschot et al. (2004) on the order of 10 3 -10 2 kg m 3 with comparable near-surface wind velocities.Our maximum snow concentration in the lowest layer of suspension (about 1.5 • 10 3 kg m 3 , Figure 3c) are rather high compared to measurements of J. Pomeroy and Male (1992) or Vionnet et al. (2014), who report concentrations on the order of 10 5 -10 3 kg m 3 at heights on the order of 1 m above the ground.Apart from this comparison, however, one has to keep in mind that the lowest mass level of our WRF domain is located at about 5 m above the ground and is representing a layer of roughly 10 m depth.The vertical structure with strong gradients in drifting snow concentration within the lower couple of meters above the ground can therefore not be accurately depicted and horizontal snow transport here is rather overestimated.

Local and Non-Local Effect of Drifting Snow Sublimation
In order to get deeper insights into the model behavior, sensitivity experiments were conducted, where each time one model parameter or setting was alternated, while the other ones were held constant.These experiments are summarized in Table 2.In this section we concentrate on the specific role of sublimation and its feedback onto the atmosphere (NO_SUBL: switched off entire sublimation process and NO_SFB: switched off sublimation feedback).In the next section, we will present the remaining experiments.As seen before in CTRL only 3% of the eroded snow mass is deposited on the ground while the rest is suspended or sublimated.Locally, the sublimation from drifting snow particles lowers the snow concentration, increases the specific humidity and decreases the temperature of the ambient air (Figures 5a-5c).That leads to a damping of the sublimation process as this is strongly controlled by relative humidity (see Section 4.2).These effects become evident when comparing the results of the CTRL simulation to the ones of NO_SUBL and NO_SFB.Switching off sublimation leads to a drastic increase in the downstream transport of snow, as well as in the area and the amount of snow deposition (Figure 6a, Table 2: about 11 times more deposition).Apart from that, snow also gets transported out of the model domain in NO_SUBL (not shown).Compared to CTRL a slight decrease of about 4.7% in snow erosion can be seen in NO_SUBL (Table 2).This decrease in erosion can be explained by two factors: on the one hand, removing snow by sublimation increases the gradient in Equation 16allowing for a stronger erosion.On the other hand, as will be seen later, the sublimational cooling leads to higher wind speeds close to the ground, which as well results in a stronger erosion in the CTRL simulation.When switching off the sublimation feedback, the aforementioned damping effect no longer exists.Thus, more snow gets sublimated, leading to smaller snow concentrations (Figure 5a) and almost no deposition of snow in NO_SFB (Figure 6b).
Although the sublimation process only takes place within a shallow layer in the summit region of the hill where snow particles are present, the effects of it can be observed over a much larger region.The direct cooling and moistening effect is visible in the entire downstream region where the air affected by sublimation is transported to (Figures 7a and 7b, below 500 m).Apart from that, the gravity wave field is affected, leading to an amplification of the gravity waves, intensified up-and downdrafts as well as increased horizontal winds (Figures 7a, 7c,  and 7d).
Despite the very different approaches to model the atmospheric fields, our observed local effects agree well with the results of Groot Zwaaftink et al. (2011).In both cases sublimation reduces the deposited snow mass, while the cooling and moistening dampens the sublimation process.Their approach, however, did not allow to represent the non-local response of the wind field to the local modifications.In conclusion, our results underline the findings of  Groot Zwaaftink et al. (2011), that the calculation of drifting snow and the sublimation feedback should be directly included in the atmospheric model in order to represent these complex interactions.

Parameter Sensitivity
We start the analysis of our sensitivity experiments with a doubling of the snow density (RHO_200).For erosion to take place, the surface shear stress has to overcome the inertia and cohesive bonds of the snow particles.Thus, a higher snow density requires a higher threshold friction velocity (see Equation 11) and results in a lower erosion flux (see Equation 10).When increasing the snow density to 200 kg m 3 this results in an overall reduction of about 35% in both the eroded and deposited snow mass compared to CTRL (Figure 8a, Table 2).Apart from the  total amounts, also the erosion area is reduced to the regions of strongest winds, where friction velocity exceeds the now required threshold.
The efficiency of the erosion process is represented by the saltation parameter e salt (Equation 10), which is set by default to 5 • 10 4 following Sauter et al. (2013).Doubling the efficiency to e salt = 10 3 (ESALT_2) results in an approximate doubling of the accumulated eroded snow mass (Table 2).However, a non-linear response is apparent in the snow deposition, which is higher by a factor of about eight compared to CTRL.This non-linearity is a result of the self-limiting damping effect in the sublimation.As seen in Section 4.3 sublimation increases the relative humidity of the ambient air, which in return dampens subsequent sublimation.Consequently, in an almost saturated environment a further increase in snow concentration only results in a small sublimation response but a stronger downstream transport and an increased deposition of snow (Figure 8b).A further sensitivity run with e salt reduced by a factor of 0.5 (ESALT_05) confirms this behavior: the erosion is reduced by about 52% while the deposition is reduced by 92% due to the stronger sublimation.
The experiments with fixed ground-level particle radii (R_1E4 and R_5E5) and fall speeds (V_01 and V_02) show the complex interplay of different processes in the model.In general, the particle size and fall speed can affect the model in several different ways.First, the fall speed is a function of particle size which on the one hand influences the downward flux of snow in suspension and thus the vertical distribution of snow.On the other hand, the ventilation speed of the snow particles that partly controls the sublimation is assumed to be equal to the particle fall speed (Sauter et al., 2013).Furthermore, the sublimation-loss rate ψ s is proportional to r 2 (Sauter et al., 2013) indicating that smaller particles are more prone to sublimation.The complex interplay of these factors is apparent in the non-linear response in the experiments with fixed fall speeds and ground-level particle radii.The faster downward flux outweighs the effect of increased ventilation speed when comparing V_02 and V_01.Larger particles are falling out more quickly and are less prone to sublimation leading to higher deposition (see R_1E4, R_5E5).
The experiments also show a weakness of the default setup.Here, the ground-level particle radius is calculated as a function of friction velocity (Equation 6) and is then given to Equation 5 to calculate the height-dependent particle radius.In a situation, where the flow detaches from the ground like in our case, this results in a drastic decrease in particle size and fallout velocity for particles that are actually decoupled from the ground.As a consequence, particle deposition is drastically decreased and sublimation increased.This problem can not occur if the ground-level particle radius is fixed and is less pronounced if the fallout velocity is fixed.We therefore recommend future users to think about fixing the ground-level particle radius if the expected flow situation is very complex, for example, in mountainous environments.Possible values might be taken from the size spectra of for example, Aksamit and Pomeroy (2016), Gordon et al. (2010), Naaim-Bouvet et al. (2011), Schmidt (1981Schmidt ( , 1986) ) or Vionnet et al. (2014) with average sizes on the order of 10 4 m.In other, less complex flow situations for example, over large ice bodies or in flat terrain the added benefit of the friction-velocity dependent particle size might still outweigh the other effects.

Outlook: Toward Real-Case Applications
Asides from idealized numerical studies, the newly implemented snow drift module can be also used for real case simulations.This application bears large potential for snow redistribution process studies in complex topography, where observations are usually sparse.An example for an area of interest embedded in complex topography with an interest in snow redistribution process is the Hintereisferner (HEF) glacier, located in the Ötztal Alps, Austria.
In an accompanying study (Voordendag et al., 2023) we apply the snowdrift module to a case study of a drifting snow event on 8 February 2021 to assess model performance with a real-case large-eddy simulation setup (Δ x = 48 m) introduced by Goger et al. (2022).Snow was initialized at the same time as the atmospheric model with the standard settings of the land-surface model NOAH-MP (Niu et al., 2011) assuming a fresh snowpack, mostly related to constraints in computational resources.The event is characterized by fresh snow fall followed by strong winds (more than 10 m s 1 ), which provides excellent conditions for drifting snow.While the majority of the study results are outlined in Voordendag et al. (2023), we show an overview of the simulated patterns of winddriven snow redistribution in the area of interest (Figure 9).The results suggest that the modeled snow depth change due to snowdrift is strongly related to the wind field and the underlying topography, as expected for this high mountain environment.In our study, snow erosion is the major contributor to snowdrift and the process removes on average 3.9 kg m 3 of snow throughout the 24 hr simulation over all glaciated areas in the domain.
The overall magnitude and patterns of snow redistribution are reproduced by the model reasonably well, while of course, fine-scale structures, that are smaller than the resolution of the model are missed.A more detailed analysis of the results and a detailed model evaluation with observations (Terrestrial Laser Scanning acquisitions, meteorological stations, etc.) can be found in Voordendag et al. (2023).

Conclusions and Outlook
In this paper we introduced a novel module to calculate drifting snow within the framework of the Weather Research and Forecasting (WRF) model.The basic idea was to have a computationally efficient module of intermediate complexity, which utilizes the existing capabilities of WRF as much as possible.For that purpose, the module is coupled to the snow scheme of the Land-Surface Model Noah-MP and employes WRF-internal modules for advection and turbulent diffusion in order to calculate the transport of drifting snow particles.The module is included directly into the calculations of WRF and, thus, is able to make use of maximum-quality wind fields and to allow for interactions of drifting snow with the background atmosphere, for example, through drifting snow sublimation.Test simulations within an idealized setting showed the physical consistency of the module.We were able to demonstrate that drifting snow sublimation influences both erosion and deposition patterns in a physically consistent way.Apart from that, we could show that the local cooling and moistening effect of drifting snow sublimation can have larger-scale non-local effects by alternating the gravity wave structure of the flow field.This underlines the necessity to include the computation of drifting snow and its feedback on the background flow directly within the atmospheric model.Due to the idealized nature of our main experiments we can only partly independently validate our results.Magnitudes of simulated snow concentrations roughly agree with field observations, while the simulated sublimation and the feedback onto the atmosphere can not be validate against field data.However, the first real-case simulations of Voordendag et al. (2023) using our snow drift module show an overall good agreement with observations.In future development steps we plan to include a radiation term into the calculation of drifting snow sublimation and to introduce a more sophisticated particle size distribution.Moreover, we intend to include the snow drift scheme into a coupling of WRF with the snow and ice model COSIPY (Sauter et al., 2020), which is currently under development.With its efficient design and incorporation into a widely used atmospheric model the module shows good potential to be used by a larger community in future studies on drifting snow for example, in mountain areas or over large ice bodies.

Figure 1 .
Figure 1.Overview of the model domain and initial conditions.The model domain is depicted in (a) with black contour lines of model topography with an interval of 100 m.The locations of the vertical cross section CS is indicated as a green line, the vertical profile P as a black X.The inital wind direction is indicated as black arrows in (a).The inital profiles of potential temperature, specific humidity and zonal wind are depicted in (b-d).

Figure 2 .
Figure 2. Average flow field for model time 3-6 hr in CTRL at 10 m above the ground (agl) in (a) and along the vertical cross section CS (cf. Figure 1a) in (b).The ucomponent of the wind is depicted as contour colors, arrows indicate the u-and v-component in (a) and u-and w-component in (b).The dashed contour line in (a) indicates a friction velocity of 0.23 m s 1 , black contour lines in (b) depict potential temperature with an interval of 2 K.

Figure 3 .
Figure 3. Snow mass change due to snow drift in CTRL accumulated from model time 3-6 hr as color contours (a).Arrows in (a) depict the horizontal wind field at 10 m agl as in Figure 2, the dashed contour line indicates a mean friction velocity of 0.23 m s 1 .Color contours depict mean snow concentration in the saltation layer (b), in 10 m agl (c) and along the vertical cross section CS (d, cf. Figure 1a).

Figure 4 .
Figure 4. Accumulated drifting snow sublimation in CTRL: vertically integrated (a) and local along the vertical cross section CS (b, cf. Figure 1a).Arrows in (a) depict the horizontal wind field at 10 m agl as in Figure 2. Thin (thick) contour lines in (b) indicate mean relative humidity of 95 (99)%.

Figure 5 .
Figure 5. Mean vertical profiles of suspended snow concentration (a), potential temperature (b), specific humidity (c) and u-component of the wind (d) at point P (cf. Figure 1a) for simulations CTRL, NO_SUBL and NO_SFB, averaged between model time 3-6 hr.Note that the profiles for NO_SUBL and NO_SFB overlap in (b-d).

Figure 6 .
Figure 6.Snow mass change due to snow drift in NO_SUBL (a) and NO_SFB (b) accumulated from model time 3-6 hr.

Figure 7 .
Figure 7. Vertical cross sections (CS, cf. Figure 1a) of the difference in mean potential temperature (a), specific humidity (b), u-component (c), and w-component of the wind (d) between CRTL and NO_SUBL, averaged between model time 3 and 6 hr (color contours).Mean potential temperature in CTRL is depicted as dashed black contour lines with 2 K interval, mean snow concentration of 5 • 10 3 kg m 3 as a thick black contour line.

Figure 9 .
Figure 9. Real-case model output from the lowest model level of accumulated total snowdrift since model initialization (colors) and horizontal 10 m wind vectors (arrows) from 8 February 2021 at eight different times.The black contours show the model topography, while the blue lines show the glacier outlines as they are represented in the model domain, data from Voordendag et al. (2023).

Table 1
WRF Namelist Options Used in the CTRL Simulation

Table 2
Summary of Sensitivity Experiments