The effect of turbulence on drifting snow sublimation

Sublimation of drifting snow, which is significant for the balances of mass and energy of the polar ice sheet, is a complex physical process with intercoupling between ice crystals, wind field, temperature, and moisture. Here a three-dimensional drifting snow sublimation model in a turbulent boundary layer is proposed. In contrast to most previous models, it takes into account turbulent diffusion of moisture from lower to higher elevations, allowing the air humidity near the surface to be undersaturated and thus sublimation to occur. From simulations with this model, we find that snow sublimation in the saltation layer near the surface dominates overall snow sublimation, despite an only marginal departure from humidity saturation ($<1\%$), because of a large particle concentration.


Introduction
[2] As one of the basic elements of Earth's environment, snow is widely distributed in high and cold areas, and its temporal and spatial evolution is crucial for the hydrological cycle [Zhou et al., 2014]. The mass and energy balance of ice shelves, which plays a key role in climate change [Cess and Yagai, 1991;Solomon et al., 2007], is largely affected by 'drifting snow', describing the transport of snow by wind, a phenomenon common in high and cold regions [Huang et al., 2016;Schmidt, 1972].
In fact, drifting snow may slow down global warming due to heat-trapping during snow crystal sublimation [Solomon et al., 2007], which occurs when the humidity of air is unsaturated [Schmidt, 1972]. As the sublimation rate of drifting snow particles is generally larger than that of ground snow due to the larger area of exposure [Dai and Huang, 2014;Huang et al., 2016], and since the accumulated time of drifting snow accounts for up to 1/3 of the duration of winter in the polar region [Mahesh et al., 2003;Mann et al., 2000], understanding the drifting snow phenomenon can play an important role in modeling mass and energy variations of the ice sheet.
However, because of moisture transportation, there are always slight deviations from the saturated state, which give rise to snow sublimation, and the impact of such deviations has not yet been quantified.
[4] Here we explore drifting snow sublimation in a natural turbulent atmospheric boundary layer and the impact of moisture transport by turbulent motions using the large eddy simulation (LES) model of the Advanced Regional Prediction System (ARPS) [Xue et al., 2001], where the drifting snow in the turbulent boundary layer is simulated using a model that accounts for the main particle transport processes [Lämmel et al., 2017]. The sublimation of mid-air snow particles is modeled by the sublimation rate formula of ice crystals [Thorpe and Mason, 1966], and its feedback effects on the air temperature and moisture is considered through introducing temperature and moisture source terms into the corresponding diffusion equations.

Generation of drifting snow
[5] The fluid governing equations of the LES model are described in the supplementary material. In the turbulent boundary layer, drifting snow occurs when the local shear stress τ is larger than the critical value t τ . The rate of aerodynamic entrainment can be expressed as [Anderson and Haff, 1991 [Doorschot and Lehning, 2002], and Clifton et al., 2006], in which g is the gravitational acceleration, and p ρ and p d are, respectively, the density and diameter of snow particle. Note that the mass density of saltating snow particles is very close to that of ice cubes because most branches of snowflakes are knocked out in their impacts with the surface [Comola et al., 2019;Groot Zwaaftink et al., 2011]. The density of air obeys the state equation of moist air, d R is the gas constant of dry air, and T and q are, respectively, the temperature and the water vapor mixing ratio of air. Furthermore, where v R is the gas constant of water vapor. Note that aerodynamic entrainment is dominated by splash entrainment in the steady state [Kok et al., 2012], thus, the exact values of the parameters in equation (1) do not have a significant effect on the final results.
[6] Particle trajectories are calculated using the Lagrangian particle tracking method, and the wind speed at the particle location is determined by the linear interpolation of surrounding grid points. Treating the snow grain as a sphere and subjecting it only to the gravity and drag force, the motion equation of a snow particle can be written as [Anderson and Haff, 1988]: where p u is the particle velocity, t the time, ν the kinematic viscosity of air, r V the relative speed between snow particle and local wind, and Re / p p r d V ν = is the particle Reynolds number. Note that we use the center point of snow particle to calculate the grid location. The drag coefficient D C is calculated as [Bagnold, 1941] [7] To describe grain-bed interactions, we use the grain scale-dependent splash function of Lämmel et al. [2017] because it has been derived mostly from first physical principles and because, in contrast to other available splash functions [Comola and Lehning, 2017;Kok and Renno, 2009], it takes into account the energy partition occurring in subsequent binary collisions between bed particles [Crassous et al., 2007;Ho et al., 2012]. The rebound probability of an impacting particle is related to its impact velocity in v and incident angle in θ , which can be expressed as [Lämmel et al., 2017]: where ˆi n d and ˆp d are the diameters of the impacting and bed particles, respectively, which are both normalized by ( are parameters related to the microscopic normal and tangential restitution coefficient n ε and t ε , respectively, of the material and ( ) is the modified mass ratio of the collision partners.
[8] According to Higa et al. [1998], the normal restitution coefficient of ice typically has a constant value qe e in the quasi-elastic regime and a tendency to decrease in the inelastic regime, which can be described by a piecewise function: , n v is the normal relative velocity of the collision partners, and c v is the critical impact velocity for the transition from the quasi-elastic region to the inelastic regime, which can be expressed as: where 0 v and 0 d are constant parameters (all parameter values are specified in Table   1).
[9] The total restitution coefficient of a rebounding particle does not vary significantly, and thus is set to be its mean value [Lämmel et al., 2017]: A rebound angle re θ is picked from the following probability distribution function [Lämmel et al., 2017]: [10] The impacting particle also generates a number of low energy ejection particles e N [Lämmel et al., 2017]: is the kinetic energy of impact particle and m e p E m g d = is the minimum transferred energy for a bed particle to be counted as ejecta, in m and e m are the mass of impact and ejection particle, respectively, and is the average energy of ejection particles. In a manner analogous to Comola and Lehning [2017], equation (10) has been modified from Lämmel et al. [2017] to take into account the mean cohesion energy c f between bed particles. Furthermore, e µ and e σ are mean value and variance of the energy of ejected particles, respectively, which can be expressed as: [11] Then, the energy of ejection particles can be obtained from the log-normal distribution as following [Lämmel et al., 2017]: and the ejection angle follows the exponential distribution with a mean value of 50  [Kok and Renno, 2009].
[12] Because the turbulence frequency spectrum of atmospheric boundary layers and  Table 1.
[13] The wind field is initialized with the standard logarithmic wind profile: [14] The snow particle size distribution follows the gamma function [Nishimura and Nemoto, 2005;Schmidt, 1982]: where p α and p β are the shape and scale parameter, respectively. During the calculation, the diameter of each takeoff particle is picked from above function. The [15] The initial air temperature is uniformly 258.15 K. At the top boundary, the relative humidity is set constant (configuration 1: 60%; configuration 2: 70%). The humidity profile becomes steady once the sublimation rate of drifting snow particles matches the moisture flux at the upper boundary, as governed by the upper relative humidity boundary condition.

Results and discussion
[16] One feature of the LES model is that turbulence eddies larger than the grid scale can be resolved, which provides a means to explore the diffusion features of moisture and heat in the turbulent flow, and which is more sophisticated than modeling diffusion through the diffusion equation [Huang et al., 2016;Wever et al., 2009].
Drifting snow in the turbulent boundary layer is generated through entrainment of surface snow directly by the wind and/or the bombardment of snow particles that are already transported by the wind ('splash') [Clifton et al., 2006;Sugiura and Maeno, 2000]. Hence, when turning on the boundary layer flow, the number of snow particles increases until an equilibrium is reached, in which snow entrainment is balanced by snow deposition at the surface.
[17] For both configurations, the simulated steady state does not depend on the choice of the initial humidity profile. This steady state is reached when the production of water vapor balances turbulent diffusion. All simulation results are averaged over a period of 120 s during the steady state. Steady state profiles of transport rate and relative humidity obtained from our simulations for both configurations are shown in Figure 1(a-c) and compared with wind tunnel and field experiments from the literature [Bintanja, 2001b;Sugiura et al., 1998;Wever et al., 2009]. Both the associated equilibrium snow transport flux profile and humidity profile are consistent with observations. Furthermore, Figure 1(d) compares the overall simulated sublimation rate for various friction velocities * u with the measurements by Bintanja [2001b].
The data in Figure 1  The observed agreement between simulations and measurements gives us confidence that the numerical model captures the essential physics of drifting snow sublimation.
However, there is no certainty about the appropriateness of the model assumptions given that the experimental data used for the model evaluation cover only the region above the saltation layer.
[18] Numerous measurements and model predictions have shown that a nearly saturated layer exists near the surface [Bintanja, 2001a;Déry and Yau, 2002;Mann et al., 2000], thus, sublimation within this layer is believed to be negligible. To test this hypothesis, we identify the locations x z at which the relative humidity exceeds certain threshold values x (e.g., for all z < 99% z , the relative humidity is larger than 99%). Figure 2(a) shows that the thickness of two representative layers ( 99% z and 95% z ) well exceeds the mean saltation height s z (defined as the elevation below which 50% of the mass flux takes place [Kok et al., 2012]). However, despite an only marginal departure from humidity saturation, the vast majority of overall snow sublimation takes place within these layers, while sublimation underneath the mean saltation layer height accounts for more than half of the overall sublimation. The reason for this counter-intuitive behavior is that the particle concentration decays exponentially with height [Pähtz and Durán, 2018], which means that even a very slight departure from humidity saturation near the surface can have a strong effect compared with highly unsaturated layers on top where much less snow particles are transported. Note that we carried out test simulations without moisture transportation, in which case the humidity near the surface is fully saturated and snow sublimation thus does not take place there. Also, modest changes of the values of the model parameters, such as the mass density of snow, do not have an effect on the quality of the results.
[19] Because turbulent diffusion does not only drive moisture transport but also tends to keep transported snow particles suspended in air, we now estimate the contribution of suspended snow particle transport to overall transport in our drifting snow simulations (configuration 2). The standard way to do this is by examining the transport height of snow particles, with particles transported in lower elevations being classified as saltating and particles in larger elevations being classified as suspended [Bagnold, 1941]. However, also particles in lower elevations can be supported by turbulent diffusion, which is why we decided for a different estimation of particle particles are far from being suspended, whereas for large ratios, suspension plays a significant role. As shown in Figure 3, there is, indeed, a large peak at unity in the distribution of this ratio. Nonetheless, about 20% of all particles have hop times ratios that are 24.5 and larger, which means that suspended particle transport is considerable.
However, note that the limited vertical dimension of the simulation domain does not allow us to simulate deep suspension clouds that sometimes occur in nature (e.g., in Antarctica [Mahesh et al., 2003;Mann et al., 2000]).
[20] Recently Sharma et al. [2018] showed that the Thorpe-Mason model, which we use to model snow sublimation, underestimates the sublimation rate when the residence time of snow particles is too short assuming that air and ice have the same initial temperature. This underestimation affects large particles more strongly because of their small residence time, which are those particles that tend to be transported in lower layers. Hence, the contribution of near-surface layers to overall sublimation is expected to be even stronger than what we find from our model.

Conclusion
[21] In this study, drifting snow sublimation in the turbulent boundary layer has been investigated using a numerical model that considers the coupling effect between snow particles, wind field, and the effects of ice phase transition on air temperature and moisture. In particular, we carried out simulations to study the impact of turbulent diffusion on snow sublimation occurring near the surface where the relative humidity is very high (e.g., > 99%). We find that snow sublimation near the surface is the main contributor to overall snow sublimation and much larger than that of higher elevations where the humidity is highly unsaturated. Our results imply that snow sublimation near the surface, neglected in most previous studies, likely plays an important role for the mass and energy balances of snow cover in high and cold regions.     [Higa et al., 1998] t ε Microscopic tangential restitution coefficients -0.85 [Supulver et al., 1995]