**Research article**
09 Nov 2018

**Research article** | 09 Nov 2018

# On the suitability of the Thorpe–Mason model for calculating sublimation of saltating snow

Varun Sharma Francesco Comola and Michael Lehning

^{1},

^{1},

^{1,2}

**Varun Sharma et al.**Varun Sharma Francesco Comola and Michael Lehning

^{1},

^{1},

^{1,2}

^{1}School of Architecture, Civil and Environmental Engineering, Swiss Federal Institute of Technology, Lausanne, Switzerland^{2}WSL Institute for Snow and Avalanche Research SLF, Davos, Switzerland

^{1}School of Architecture, Civil and Environmental Engineering, Swiss Federal Institute of Technology, Lausanne, Switzerland^{2}WSL Institute for Snow and Avalanche Research SLF, Davos, Switzerland

**Correspondence**: Varun Sharma (varun.sharma@epfl.ch)

**Correspondence**: Varun Sharma (varun.sharma@epfl.ch)

Received: 08 Feb 2018 – Discussion started: 15 Mar 2018 – Revised: 03 Oct 2018 – Accepted: 12 Oct 2018 – Published: 09 Nov 2018

The Thorpe and Mason (TM) model for calculating the mass lost from a sublimating snow grain is the basis of all existing small- and large-scale estimates of drifting snow sublimation and the associated snow mass balance of polar and alpine regions. We revisit this model to test its validity for calculating sublimation from saltating snow grains. It is shown that numerical solutions of the unsteady mass and heat balance equations of an individual snow grain reconcile well with the steady-state solution of the TM model, albeit after a transient regime. Using large-eddy simulations (LESs), it is found that the residence time of a typical saltating particle is shorter than the period of the transient regime, implying that using the steady-state solution might be erroneous. For scenarios with equal initial air and particle temperatures of 263.15 K, these errors range from 26 % for low-wind, low-saturation-rate conditions to 38 % for high-wind, high-saturation-rate conditions. With a small temperature difference of 1 K between the air and the snow particles, the errors due to the TM model are already as high as 100 % with errors increasing for larger temperature differences.

- Article
(1587 KB) -
Supplement
(74109 KB) - BibTeX
- EndNote

Sublimation of drifting and blowing snow has been recognized as an important component of the mass budget of polar and alpine regions (Lenaerts et al., 2012; Liston and Sturm, 2004; Vionnet et al., 2014; van den Broeke et al., 2006). Field observations and modelling efforts focused on Antarctica have highlighted the fact that precipitation and sublimation losses are the dominant terms of the mass budget in the katabatic flow region as well as the coastal plains (van den Broeke et al., 2006). Even though precipitation is challenging to measure accurately, methods to measure it exist, for example, using radar (Grazioli et al., 2017) or snow depth change (Vögeli et al., 2016). In comparison, sublimation losses are even harder to measure and can only be calculated implicitly using measurements of wind speed, temperature and humidity. Thus, in regions where sublimation loss is a dominant term of the mass balance, it is also a major source of error. This error ultimately results in errors in the mass accumulation of ice on Antarctica, which is a crucial quantity for understanding sea-level rise and climate change (Lenaerts et al., 2012; Rignot et al., 2011; Rémy and Frezzotti, 2006).

Aeolian transport of snow can be classified into three modes, namely creeping, saltation and suspension. Creeping consists of heavy particles rolling and sliding along the surface of the snowpack either due to form drag or bombardment due to impacting particles. Saltation consists of particles being transported along the surface via short, ballistic trajectories with heights mostly less than 10 cm and involves mechanisms of aerodynamic entrainment along with rebound and splashing of ice grains (Comola and Lehning, 2017; Doorschot and Lehning, 2002). Suspension, on the other hand, refers to transport of small ice grains at higher elevations and over large distances without contact with the surface. Current calculations of sublimation losses are largely restricted to losses from ice grains in suspension. This is true for both field studies (Mann et al., 2000), where sublimation losses are calculated using measurements, usually at the height of 𝒪(1 m), and in mesoscale modelling studies (Déry and Yau, 2002; Groot Zwaaftink et al., 2011; Vionnet et al., 2014; Xiao et al., 2000), where the computational grids and time steps are too large to resolve flow dynamics at saltation length and timescales. Mass loss in the saltation layer is hard to measure and is neglected based on the justification that the saltation layer is saturated. However, recent studies using high-resolution steady-flow, Reynolds-averaged Navier–Stokes (RANS)-type simulations (Dai and Huang, 2014) claim that sublimation losses in the saltation layer are not negligible, particularly for wind speeds close to the threshold velocities for aeolian transport, wherein a majority of aeolian snow transport occurs via saltation rather than suspension.

The coupled heat and mass balance equations of a single ice particle immersed in turbulent flow are

where *m*_{p}, *T*_{p} and *d*_{p} are the mass, temperature and diameter of
the particle respectively that vary with time, *c*_{i} is the specific heat
capacity of ice, *L*_{s} is the latent heat of sublimation, 𝔎 is
the thermal conductivity of moist air and 𝔇 is the mass
diffusivity of water vapour in air. Transfer of heat and mass is driven by
differences of temperature and vapour density between the particle surface
(*T*_{p}, *ρ*_{w,p}) and the surrounding fluid (*T*_{a,∞},
*ρ*_{w,∞}). The vapour density at the surface of the ice particle is
considered to be the saturation vapour density for the particle temperature.
The transfer mechanisms are enhanced by turbulence, the effect of which is
parameterized by the Nusselt (𝒩*u*) and Sherwood (𝒮*h*)
numbers. 𝒩*u* and 𝒮*h* are related to the
relative speed $\left(\left|{\mathit{u}}_{\mathrm{rel}}\right|\right)$ between the air and the
particle via the particle Reynolds number (*Re*_{p}) as

where *ν*_{air} is the kinematic viscosity of air, and *Pr* and *Sc* are the
Prandtl and Schmidt numbers.

Thorpe and Mason (1966) solved the above coupled Eqs. (1) and (2) by (a) neglecting the thermal inertia of the ice particle, thus effectively stating that all the heat necessary for sublimation is supplied by the air and (b) considering the temperature difference between the particle and surrounding air to be small, thereby allowing for Taylor series expansion of the Clausius–Clapeyron equation and neglecting higher-order terms, resulting in their formulation for the mass loss term as

where *ρ*_{s}(*T*_{a,∞}) is the saturation vapour density
of air surrounding the particle, saturation-rate ${\mathit{\sigma}}_{*}={\mathit{\rho}}_{\mathrm{w},\mathrm{\infty}}/{\mathit{\rho}}_{\mathrm{s}}\left({T}_{\mathrm{a},\mathrm{\infty}}\right)$, *M* is the
molecular weight of water and *R* is the universal gas constant. The above
formulation has been used extensively to analyse data collected in the field
(Mann et al., 2000), wind tunnel experiments (Wever et al., 2009) and numerical
simulations of drifting and blowing snow (Déry and Yau, 2002; Groot Zwaaftink et al., 2011; Vionnet et al., 2014).
In the modelling studies, the mass loss term is computed using
Eq. (4) and is added, with proper normalization, to the
advection–diffusion equation of specific humidity, while the latent heat of
sublimation multiplied by the mass loss term is added to the corresponding
equation for temperature (Groot Zwaaftink et al., 2011).

Two observations motivated us to investigate the suitability of the TM model
for sublimation of saltating snow particles. Firstly, the TM model assumes
that all the energy required for sublimation is supplied by the air. This
assumption was tested by Dover (1993), who compared the potential rates
of cooling of particles with that of the surrounding air due to sublimation.
Using scale analysis, Dover (1993) formulated the quantity $\mathit{\xi}=\mathrm{6}{\mathit{\rho}}_{\mathrm{air}}{c}_{\mathrm{p},\mathrm{air}}/\mathit{\pi}{\mathit{\rho}}_{\mathrm{i}}{c}_{\mathrm{i}}{\stackrel{\mathrm{\u203e}}{{d}_{\mathrm{p}}}}^{\mathrm{3}}N$, where
$\stackrel{\mathrm{\u203e}}{{d}_{\mathrm{p}}}$ is the mean particle diameter, *N* is the particle number
density and *ρ*_{i} is the density of ice, and showed that, for *ξ**>>*1, it
can be accurately considered that the heat necessary for sublimation comes
from the air. For standard values for an ice particle in suspension,
$\stackrel{\mathrm{\u203e}}{{d}_{\mathrm{p}}}=\mathrm{50}\phantom{\rule{0.125em}{0ex}}\mathrm{\mu}\mathrm{m}$ and *N*∼𝒪(10^{6}),
this condition is easily met (*ξ*∼𝒪(10^{3})).
However, if we input values typical for saltation, i.e. $\stackrel{\mathrm{\u203e}}{{d}_{\mathrm{p}}}=\mathrm{200}\phantom{\rule{0.125em}{0ex}}\mathrm{\mu}\mathrm{m}$ and *N*∼𝒪(10^{8}), *ξ*∼𝒪(1), and the condition is not met. Thus, for sublimation of
saltating particles, it is important to consider the thermal inertia of the
particles. A similar conclusion was reached in other modelling studies on
topics of heat and mass exchange between disperse particulate matter in
turbulent flow such as small water droplets in heat exchangers
(Russo et al., 2014) and sea sprays (Helgans and Richter, 2016).

Secondly, Eq. (4) computes mass loss as being directly proportional
to *σ*_{*} and neglects the temperature difference between the particle
and air. Equation (4) thus predicts a mass loss even in extremely
high-saturation-rate conditions, whereas immediate deposition of water vapour would
occur on a particle even slightly colder than the air. Indeed, some field
experiments have reported deposition as opposed to sublimation which was
expected on the basis of the measured undersaturation of the environment,
particularly near coastal polar regions (Sturm et al., 2002). A simple
everyday observation illustrates this fact clearly: there is immediate
deposition of vapour and formation of small droplets on the surface of a cold
bottle of beer even in room conditions with moderate humidity!

Motivated by the observations described above, in this article we describe four numerical experiments for which we compare differences between the fully numerical and the Thorpe and Mason (1966) solutions (referred to as NUM and TM approaches respectively). In experiments I and II, we numerically solve Eqs. (1) and (2) and compare the results with Eq. (4) for physically plausible values of a saltating ice particle. Results of these tests are presented in Sect. 2. High-resolution large-eddy simulations (LESs) of the atmospheric surface layer with saltating snow are performed for a range of environmental conditions to compute the differences between the NUM and TM approaches in realistic wind-driven saltating events. These results are presented in Sect. 3. A summary of the article is given in Sect. 4.

We consider an idealized scenario where a solitary spherical ice particle is held still in a turbulent airflow with constant mean speed, temperature and undersaturation. The evolution of the mass, diameter and temperature of the ice particle is calculated using both the NUM and TM models and an intercomparison is made. This scenario is similar to the wind-tunnel study performed by Thorpe and Mason (1966), who measured mass loss of solitary ice grains suspended on fine fibres. In this scenario, we consider that the heat and mass transfer between the ice particle and the air changes the mass and temperature of the particle only, while the mass and energy anomalies in the air are rapidly advected and mixed. This implies that the environmental conditions subjected to the ice particle remain constant. While it can be expected that the environmental conditions will vary along the trajectory of a ice particle undergoing saltation or suspension, it is nevertheless useful to perform this analysis, as it reveals important characteristics about the heat and mass evolution of a ice particle during sublimation and about the approximations used to derive the TM model.

For the NUM approach, Eqs. (1) and (2) are solved in a coupled manner using a simple
first-order finite-differencing scheme for time stepping with a time step of
50 µs. For the TM approach, Eq. (4) is used with a
similar numerical set-up as for the NUM approach. In the TM approach, particle
temperature is not considered and the mass and energy transfer is determined
only by air temperature and saturation rate. The initial particle diameter
(*d*_{p,IC}) is 200 µm and the airflow
temperature is 263.15 K for both the NUM and TM approaches. We
use a constant air speed of 5 m s^{−1} resulting in
*Re*_{p}=80, 𝒩*u*=6.7 and 𝒮*h*=6.5 (using
Eq. 3). The values used here are typical of a saltating ice
particle (Kok and Renno, 2009; Thorpe and Mason, 1966; Vionnet et al., 2014).

In Experiment I, we study the heat and mass output from a sublimating ice
grain as a function of time. In the first case, Experiment I-A, we consider
the effect of three different values of airflow saturation rate
(${\mathit{\sigma}}_{*}=\mathrm{0.8}$, 0.9 and 0.95) on differences between NUM and TM
solutions. The NUM approach requires specification of the initial condition
for the particle temperature (*T*_{p,IC}). In Experiment I-A,
(*T*_{p,IC}) is taken to be the same as the airflow temperature
for the NUM approach, i.e. 263.15 K. Results for Experiment I-A are shown in
Fig. 1a–c, with subfigure (a) showing the mass output rate,
*F*_{M} and subfigure (b) showing the heat output rate, *F*_{Q}. Note that in
this figure and subsequent figures, $+(-)$ signifies mass and heat gained
(lost) by the air. Since we keep the temperature and undersaturation of the
air constant, the solutions of the TM approach are steady-state solutions
with constant heat and mass transfer rates as seen in Fig. 1a and b. On the
other hand, since the NUM approach solves the coupled equations that consider
the evolution of the particle temperature, the heat and mass transfer rates
evolve with time.

It can be seen that the NUM solutions initially evolve with time and
reconcile with the steady-state TM solutions after a transient regime of
about 0.3 s. Since the initial temperature of the particle is the same
as the air, there is no heat transfer between the air and the particle (see
the second term of the right-hand side of Eq. 1) initially. Thus, all heat transfer
rates are initially zero for the NUM case in Fig. 1b. The undersaturation of
the air forces mass transfer from the ice particle to the air and the energy for
the phase change comes from the internal energy of the ice particle. This
causes the particle temperature to drop (see Fig. 2 below). With the particle
now colder than the air, heat transfer from the air to the particle commences
and ultimately, the energy for sublimation comes entirely from the heat
extracted from the air. The initial dynamics of the heat and mass transfer
are completely neglected by the TM approach. In subfigure (c), the errors $\left(\mathrm{Err}\left(t\right)=\left({\int}_{\mathrm{0}}^{t}{F}_{\mathrm{NUM}}\mathrm{d}t/{\int}_{\mathrm{0}}^{t}{F}_{\mathrm{TM}}\mathrm{d}t-\mathrm{1}\right)\times \mathrm{100}\right)$ for mass, Err_{M} and heat, Err_{Q} are
shown. The errors reduce dramatically with time (for example, 15 % at
0.3 s) and interestingly do not depend on the saturation rate of the
airflow.

In the following case, Experiment I-B, similar simulations as in Experiment
I-A are performed, but with ${\mathit{\sigma}}_{*}=\mathrm{0.95}$, while the initial temperature
difference between the particle and the air is varied as ${T}_{\mathrm{p},\mathrm{IC}}-{T}_{\mathrm{Air}}=-\mathrm{2},-\mathrm{1},\mathrm{1},\mathrm{2}\phantom{\rule{0.125em}{0ex}}\mathrm{K}$. The results are shown in
Fig. 1d–f. It is interesting to note that, for each of the
four cases considered, the TM solution predicts sublimation of the particle
(consistent with ${\mathit{\sigma}}_{*}<\mathrm{1}$; see numerator of RHS of Eq. 3). On the
other hand, for cases with colder particles, the NUM solutions show that
there is initially deposition on the particle, along with larger values of
heat absorbed from the air. Correspondingly, in the cases with particles
being warmer than the air, the mass loss is much higher in the NUM solution
than that computed by the TM solution, while the heat gained by the particle
is also much higher. These higher differences are reflected in the Err_{M}
and Err_{Q} curves in subfigure (f), where errors are found to be an order
of magnitude higher that those in subfigure (c).

We define relaxation time (*τ*_{relaxation}) as the time
required for the NUM solution to reconcile with the TM solution. The
importance of this quantity lies in the fact that if the residence time of a
saltating ice grain in air is shorter than *τ*_{relaxation}, the TM
approach is likely to be erroneous and the NUM approach would be required. It
is intuitive that *τ*_{relaxation} increases with *d*_{p} on account of
increasing inertia and decreases with $\left|{\mathit{u}}_{\mathrm{rel}}\right|$ due to more
vigorous heat and mass transfer. Experiment I was repeated for values of
*d*_{p} and $\left|{\mathit{u}}_{\mathrm{rel}}\right|$ ranging between (50−1000 µm) and $\left(\mathrm{0}-\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}\right)$
respectively. The upper-bound of the wind-speed range is quite high and it is
extremely unlikely to find $\left|{\mathit{u}}_{\mathrm{rel}}\right|>\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ in
naturally occurring aeolian transport. Numerical results indeed confirm our
intuition and it is found that, for any given value of $\left|{\mathit{u}}_{\mathrm{rel}}\right|$,
*τ*_{relaxation} is found to be $\propto {d}_{\mathrm{p}}^{\mathit{\alpha}}$, where *α*(∼1.65). Furthermore, *τ*_{relaxation} decreases
monotonically with increasing $\left|{\mathit{u}}_{\mathrm{rel}}\right|$ for a given value of
*d*_{p}. For *d*_{p}=200 µm, the plausible values of
*τ*_{relaxation} are found to lie between 0.28 and 1.5 s (for $\left|{\mathit{u}}_{\mathrm{rel}}\right|$ = 10 and 0 m s^{−1} respectively). Interestingly,
*τ*_{relaxation} is not found to depend on either the saturation rate of
air or the difference between the initial particle and air temperature. Plots
of *τ*_{relaxation} are highly relevant to discussion in
Sect. 3 and presented there.

In Fig. 2, the evolution of particle diameter (*d*_{p}) and temperature (*T*_{p}) is presented with subfigures
(a) and (b) describing the evolution for simulations in
Experiment I-A, with (c) and (d) being the corresponding results from
Experiment I-B. In Experiment I-A, the particle diameters reduce linearly
with time for both the NUM and TM approaches with the more shrinking (or in
other words, sublimation) in the NUM solutions. More interesting is the
evolution of the particle temperature, where the particle undergoes
significant cooling due to sublimation and ultimately achieves a constant
temperature. For example, in the case for ${\mathit{\sigma}}_{*}=\mathrm{0.8}$, the particle
temperature is ultimately 0.85 K lower than the initial particle temperature
of 263.15 K. Note that, for the TM approach, particle temperature is of no
consequence and it is shown simply for reference.

Following the results of Experiment I, in Experiment II, we explore the parameter
space of $\left({\mathit{\sigma}}_{*},\phantom{\rule{0.25em}{0ex}}{T}_{\mathrm{p},\mathrm{IC}}-{T}_{\mathrm{Air}}\right)$ and compute the
total mass ($M={\int}_{\mathrm{0}}^{t}{F}_{\mathrm{M}}\mathrm{d}t$) and total heat ($Q={\int}_{\mathrm{0}}^{t}{F}_{\mathrm{Q}}\mathrm{d}t$) output by a sublimating ice grain for a finite time of *t*=0.5 s. Results shown in Fig. 3 subfigures (a) and (b) provide
a comparison of the total mass lost using the NUM and TM solutions
and the corresponding error is shown in subfigure (c). Similar
figures are presented for the total heat lost/gained by the air in subfigures
(d–f). The inclusion of the inertial terms essentially causes the contours to
be sloped for the NUM solution, while the TM solutions do not depend on
${T}_{\mathrm{p},\mathrm{IC}}-{T}_{\mathrm{Air}}$ as expected. The error between the NUM and TM solutions
are accentuated at high-saturation-rate regimes, with errors larger than 30 % for ${\mathit{\sigma}}_{*}>\mathrm{0.8}$.

In summary, experiments I and II highlight the fact that, during the
sublimation of an ice grain, there exists a finite, well-defined transient
regime before the NUM solutions match the steady-state TM solutions.
Furthermore, the NUM and TM solutions diverge rapidly with slight temperature
differences between the particle and the air and with increasing *σ*_{*}
(which is a cause of concern since in snow-covered environments, the air
usually is highly saturated). The results described above prompt an
interesting question: are the residence times of saltating ice particles
comparable to *τ*_{relaxation}? We use large-eddy simulations to answer
this question in the following section.

## 3.1 Experiments III and IV: simulation details

To further understand the implications of the differences between the NUM and
the TM approach, we performed LESs of the atmospheric surface layer with an
erodible snow surface as the lower boundary. We describe here only the main
details of the LESs that are relevant to our discussion with full model
description along with equations presented in Supplement Sect. S1. The
LES solves filtered Eulerian equations for momentum, temperature and specific
humidity on a computational domain of 6.4 m×6.4 m in the horizontal directions with vertical extent of the domain
being 6.4 m as well. The snow surface, which constitutes the
lower boundary of the computational domain, consists of spherical snow
particles with a log-normal size distribution with a mean particle diameter
of 200 µm and standard-deviation of 100 µm. The
coupling between the erodible snow-bed and the atmosphere is modelled through
statistical models for aerodynamic entrainment (Anderson and Haff, 1988),
splashing and rebounding of particle grains (Kok and Renno, 2009), which have been
updated recently by Comola and Lehning (2017) to include the effects of cohesion and
heterogeneous particle sizes. The use of these models essentially allows for
overcoming the immense computational cost of resolving individual
grain-to-grain interactions and allow us to consider the snow-surface as a
*bulk* quantity rather than a collection of millions of individual
snow particles. Once the ice grains are in the flow, their equations of
motion are solved in the Lagrangian frame of reference with only
gravitational and turbulent form drag forces included. Since the particle
velocities are known, $\left|{\mathit{u}}_{\mathrm{rel}}\right|$ is calculated explicitly and used
to compute *Re*_{p}, 𝒩*u* and 𝒮*h*. The horizontal
boundaries of the domain are periodic and the lower boundary condition (LBC)
for velocity uses flux parameterizations based on Monin–Obukhov similarity
theory, additionally corrected for flux partition between fluid and particles
between the wall and the first flow grid point
(Raupach, 1991; Shao and Li, 1999). The LBC for scalars (temperature and
specific humidity) are flux-free and thus the only source/sink of heat and
water vapour in the simulations is through the interaction of the flow with
the saltating particles.

All simulations are performed on a grid of $\mathrm{64}\times \mathrm{64}\times \mathrm{128}$ grid points with a uniform grid in the horizontal directions and a stretched grid in the vertical. A stationary turbulent flow is allowed to first develop, following which, the snow surface is allowed to be eroded by the air. All physical constants and parameters along with additional details of the numerical set-up are provided in Sect. S2. The LES in the configuration used in this study resembles the classic case of LES of a channel flow common in computational fluid dynamics research. The LES code has been developed in-house for last many years beginning with (Albertson and Parlange, 1999) and has been used and validated for various atmospheric boundary layer problems such as flows over heterogeneous surfaces (Bou-Zeid et al., 2004), hills (Diebold et al., 2013), diurnal cycles (Kumar et al., 2006), urban canopies (Giometto et al., 2016) and wind farms (Calaf et al., 2010; Sharma et al., 2017). The same code was used previously for modelling snow saltation by Groot Zwaaftink et al. (2014).

For the TM approach, Eq. (4) is used to compute the specific
humidity and (by multiplying with the latent heat of sublimation) heat
forcing due to each ice grain in the flow. On the other hand, for the NUM
approach, Eqs. (1) and (2) are
solved and only the turbulent transfer of heat between the air and the
particle (second term in RHS of Eq. 1) acts as a
heat forcing on the flow. An implication of the NUM approach is that the
particle temperature evolves during the ice grain's motion and this
necessitates providing an initial condition for the particle temperature
(*T*_{p,IC}).

The principle aims of Experiment III are to firstly quantify particle
residence times (PRTs) and their dependence on wind speeds and relative
humidities and secondly, compute the differences in the heat and mass output
between the NUM and the TM approaches during saltation of snow with complete
feedback between the air and the particles. PRT is defined as the total time
the particle is airborne and in motion, including multiple hops across the
surface. Note that the PRT is not computed for particles in suspension, i.e.
particles that stay aloft and never return to the surface. Towards this goal,
simulations are performed, each with a combination of initial surface stress,
${u}_{*}\in \left\{\mathrm{0.4},\mathrm{0.6},\mathrm{0.8}\right\}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and initial saturation rate, ${\mathit{\sigma}}_{*}\in \left\{\mathrm{0.3}\phantom{\rule{0.125em}{0ex}},\phantom{\rule{0.125em}{0ex}}\mathrm{0.6}\phantom{\rule{0.125em}{0ex}},\phantom{\rule{0.125em}{0ex}}\mathrm{0.9}\right\}$. These values are classified as low (L),
medium (M) and high (H) and correspond to wind speed at 1 m height above the
surface of 11, 16.3 and 21.8 m s^{−1}. Note that during fully
developed snow transport, the particles in the air impart drag on the flow
causing the flow to decelerate. The wind speeds at 1 m during fully developed
saltation are 8.77, 11.34 and 13 m s^{−1}. The simulations
are named as U*α*-R*β*, where $\left(\mathit{\alpha}\phantom{\rule{0.125em}{0ex}},\phantom{\rule{0.125em}{0ex}}\mathit{\beta}\right)\in \left\{\text{L, M, H}\right\}$. Each combination is
simulated independently for the NUM and TM approaches resulting in a total of
eighteen simulations. Experiment III is limited to simulating the
usual case where the initial air temperature (*T*_{Air,IC}) is the same as *T*_{p,IC}.

Experiment IV is aimed at exploring the implications of differences between
the two approaches in cases where *T*_{Air,IC} is significantly different
from *T*_{p,IC}. Such conditions can occur in nature during events such as
marine-air intrusions, katabatic winds, spring-season saltation events and
winter flows over sea-ice floes, where significant temperature differences
between the air and snow surface are likely. We repeat the low-wind case of
Experiment III with ${u}_{*}=\mathrm{0.4}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and choose the initial
saturation rate to be 0.95, motivated by results in Experiment II where
errors were found to increase with increasing saturation rate. Simulations
(named as UL-T(*γ*), where ${T}_{\mathrm{Air},\mathrm{IC}}-{T}_{\mathrm{p},\mathrm{IC}}=\mathit{\gamma}$) are
performed once again for each of two approaches with $\mathit{\gamma}\in \left\{\pm \mathrm{1},\pm \mathrm{2.5},\pm \mathrm{5}\phantom{\rule{0.125em}{0ex}}\mathrm{K}\right\}$ resulting in a total of 12 simulations. In all simulations
performed for experiments III and IV, ${T}_{\mathrm{p},\mathrm{IC}}=\mathrm{263.15}\phantom{\rule{0.125em}{0ex}}\mathrm{K}$. It is
important to note that the initial condition for particle temperature (*T*_{p,IC}) is fixed throughout the simulation period, which
essentially means that surface temperature is kept constant. This is
consistent with the imposed zero flux of heat at the surface. This imposition
will be justified a posteriori in the following section.

## 3.2 Results

In this section, results from the LESs performed for experiments III and IV are presented. Note that only the relevant results are presented, namely (a) particle residence times as a function of particle diameters and different forcing set-ups and (b) differences between the NUM and TM approaches for calculating average mass and heat transfer rates during saltation. Other results, for example, vertical profiles of mean and turbulent quantities, although interesting, are relegated to the Supplement as their analysis is out of scope of the current work. Additionally, two video illustrations (Supplement Movies M1 and M2; see Sect. S4) of an LES are provided to help visualize and frame the context of the simulations.

### 3.2.1 Particle residence times versus *τ*_{relaxation}

As mentioned in the concluding lines of the Sect. 2,
the principle quantity of interest is the PRT of saltating ice grains. In
Fig. 4a, the mean and median PRT of five different simulations
of Experiment III are shown as a function of the particle diameter.
Additionally, values of *τ*_{relaxation} computed in Experiment I for wind
speeds ranging from 0 to 10 m s^{−1} are also shown in the shaded
region. Recall that the shaded region represents all the plausible values of
*τ*_{relaxation} in naturally occurring aeolian transport. As examples,
*τ*_{relaxation} trends for three wind speeds, 0, 1 and 10 m s^{−1}
are shown and the power-law dependence can clearly be seen. It is found that
*τ*_{relaxation} is comparable to the PRT of saltating grains with
diameters between 125 and 225 µm. For 200 µm, the
mean PRT is found to be 0.6 s, while the median PRT is 0.2 s,
which is outside the range of admissible values of *τ*_{relaxation}. For
particles larger that 225 µm, the PRTs are an order of
magnitude smaller than plausible values of *τ*_{relaxation} and therefore
the TM model is likely to provide wrong values of mass loss. On the other
hand, lighter particles with diameters smaller than 100 µm have
much longer PRTs and the TM model is therefore valid. This proves that, while
the TM model is applicable for a majority of particles in suspension, it is
likely to cause errors for particles in saltation.

Results presented in Fig. 4a provide two additional insights.
Firstly, it is quite interesting to note that particles larger than 100 µm have the same mean PRT irrespective of low, medium or high
mean wind speeds. This means that the dynamics of the heavier particles are
unaffected by different mean wind speeds simulated in Experiment III, which
is consistent with the notion of self-organized saltation, which has recently
been shown by Paterna et al. (2016). For particles smaller than 100 µm, the mean PRTs increase with wind speed. Secondly, the
initial saturation rate does not seem to affect the PRT statistics for medium
and high-wind conditions and the UM- and UH- curves for different *R* values
overlap (this is the reason only five PRTs are shown in
Fig. 4a). In these cases turbulence is sufficient to rapidly
mix any temperature anomaly due to sublimation throughout the surface layer.
On the other hand, the mean PRTs of particles smaller than 75 µm, decrease with decreasing initial *σ*_{*}. A preliminary hypothesis
is that in low-wind conditions (UL), low initial saturation rate results in
more sublimation and cooling near the surface, resulting in suppression of
vertical motions. Even though this is an interesting result, further research
is needed to definitively link drifting snow sublimation to lower PRTs. We
leave further exploration of this phenomenon for a future study with some
preliminary analysis provided in Sect. S3.

The PRT distributions are found to be quasi-exponential with long tails, thus
resulting in large differences in mean and median PRTs shown in
Fig. 4a. These distributions are also strongly dependent on the
particle diameter. As an illustration, in Fig. 4b, cumulative
distributions of PRTs are shown for four particle diameters along with the
corresponding range of plausible *τ*_{relaxation} values. For the mean
particle diameter of 200 µm, we find that between 65 % and 85 %
of particles have PRTs shorter than *τ*_{relaxation}, whereas for the
75 µm particles, at most 30 % particles lie below the maximum
*τ*_{relaxation} threshold. This reinforces the fact that applying the
steady-state TM solution to sublimating ice grains in saltation could be
potentially erroneous.

### 3.2.2 Differences in total mass loss between NUM and TM models

We can directly assess the implications of differences in grain-scale
sublimation between the two approaches on total mass loss rates during
saltation at larger spatial scales as simulated using LESs in experiments III
and IV. In Fig. 5, we compare the total 15 min averaged rate of
mass loss computed in all cases in Experiment III (subfigure a) and
Experiment IV (subfigure c) using the NUM and the TM approaches with
corresponding errors shown in subfigures b and d respectively. Recalling the
adopted convention of $+(-)$ as gain (loss) of flow quantities, it can be seen
in Experiment III that sublimation increases with *u*_{*} and decreases with
*σ*_{*}. The errors, on the other hand, increase with increasing values of
both *u*_{*} and *σ**. The increase in error with *u*_{*} is mostly
due to the fact that an increase in *u*_{*} proportionally increases the
total mass entrained by air (see Supplement Fig. S1). The increase in
error with increasing *σ*_{*} is in accordance with analysis done in
Experiment II (see Fig. 3c, f) where it was shown that the
NUM and TM solutions diverge with increasing saturation rate. The least
error, 26 % is found for case UL-RL (i.e. ${u}_{*}=\mathrm{0.4}$, ${\mathit{\sigma}}_{*}=\mathrm{0.3}$),
while the largest error, 38 %, is found for UH-RH (${u}_{*}=\mathrm{0.8}$,
${\mathit{\sigma}}_{*}=\mathrm{0.9}$). Overall, for all the simulation combinations, the NUM
approach computes larger mass loss than the TM approach.

Experiment IV highlights the effect of temperature difference between
particle and air on sublimation. As shown in Fig. 5c, the mass
output is found to be negative (deposition) for the NUM solutions when the
air is warmer than the particles (i.e. cases UL-T(*γ*) with *γ*>0). This is contrary to the TM solutions, which indicate sublimation. In
cases with *γ*<0, the NUM approach shows a much higher sublimation
rate than the TM solutions. This occurs firstly due to higher vapour pressure
at the grain surface that results in enhanced vapour transport and secondly
because the warmer particles heat the surrounding air via sensible heat
exchange, causing the relative humidity to decrease. Errors increase
dramatically from an already high 100 % for UL-T(+1) to 800 % for UL-T(+5).
Simulations performed for medium- and high-wind cases in Experiment IV showed
even higher errors, which were similar to results in Experiment III and are shown here.

In this
article, we revisit the Thorpe and Mason (1966) model used to calculate
sublimation of drifting and blowing snow and check its validity for saltating
ice grains. We highlight the fact that solutions to unsteady heat and mass
transfer equations (NUM solutions) converge to the steady-state TM model
solutions after a *relaxation* time, denoted by *τ*_{relaxation}, which has a power-law dependence on the particle diameter and is inversely
proportional to the relative wind speed. Through extensive LESs of snow
saltation, we compute the statistics of the PRTs as a function of their
diameters and find them to be comparable to *τ*_{relaxation}. This helps
to explain the difference between mass output when using the NUM model to the TM
approach, also computed during the same LESs. The NUM approach computes
higher sublimation losses ranging from 26 % in low-wind, low-saturation-rate
conditions to 38 % in high-wind, high-saturation-rate conditions. Another set
of numerical experiments explore the role of temperature differences between
particle and air temperature in inducing differences between NUM and TM
solutions. We find the effect to be extremely dramatic with errors of 100 %
for a temperature difference of 1 K with increasing errors for larger
temperature perturbations. In general, the two solutions are found to diverge
rapidly as the saturation rate tends towards 1. The results showing
differences in mass output between the NUM and TM approaches in the LESs in
experiments III and IV, with complete feedback between particles and the air
are thus shown to be closely correlated to the results from extremely
idealized simulations of heat and mass transfer from a solitary ice grain in
experiments I and II.

The LES results do come with a few important caveats. Firstly, the temperature and specific humidity fluxes at the surface are neglected. In other words, particles lying on the surface are considered to be dormant and do not exchange heat or mass with the air. A corollary to neglecting the scalar fluxes at the surface is that the initial condition for temperature of the particles entering the flow is fixed. This may be justified by considering that during drifting and blowing snow events, the friction velocity at the surface drops dramatically. This fact has been observed in both in experiments (Walter et al., 2014) and in our current LESs (see Fig. S2). This implies that direct turbulent exchange between the surface and air is curtailed and instead, the dominant exchange occurs between airborne particles and the air. In fully developed snow transport events, this is most likely to be true and only in intermittent snow-transport events will the surface fluxes be relatively important. This is nevertheless an important assertion that shall be more closely examined in future studies involving a full surface energy balance model, where the evolving temperature of the saltating ice grains prior to deposition is taken into account while calculating snow-surface temperatures.

Further work is required to make concrete improvements to modelling of sublimation of saltating snow, especially in large-scale models that do not explicitly resolve saltation dynamics. One potential approach is to modify the Monin–Obukhov-based lower-boundary conditions for heat and moisture to account for particle temperature during blowing snow events. An ancillary outcome of this study is the discovery that buoyancy can affect the dynamics of lighter snow particles (with diameters less than 75 µm) and decrease their residence times. Investigating this phenomenon requires a detailed analysis of turbulent structure within the saltation layer and is left for future publications.

In conclusion, analogous to the role played by saltating grains in efficient momentum transfer to the underlying granular bed, the NUM approach can be considered as an efficient transfer of heat and mass between the flow and the underlying snow surface, albeit with a closer physical relationship between the thermodynamics of the snow surface and that of the air. Thus, along with momentum balance of blowing snow particles, particle temperature and its thermal balance must also be taken into account.

The model outputs as well as the code for generating the data used in Sect. 2 can be requested from the authors.

The supplement related to this article is available online at: https://doi.org/10.5194/tc-12-3499-2018-supplement.

VS and ML formulated the research, VS and FC performed the simulations and analysis, and VS, FC, and ML wrote the paper.

The authors declare that they have no conflict of interest.

We acknowledge the support of the Swiss National Science Foundation (grant no.
200021-150146) and the Swiss Super Computing Center (CSCS) for providing
computational resources (project s633). We additionally thank Marco Giometto
for providing the original version of the LES code and Hendrik Huwald,
Tristan Brauchli, Annelen Kahl and Celine Labouesse for illuminating
discussions and important suggestions in improving the manuscript. No new
data were used in producing this manuscript.

Edited by: Florent Dominé

Reviewed by: two anonymous referees

Albertson, J. D. and Parlange, M. B.: Natural integration of scalar fluxes from complex terrain, Adv. Water Resour., 23, 239–252, https://doi.org/10.1016/S0309-1708(99)00011-1, 1999. a

Anderson, R. S. and Haff, P. K.: Simulation of Eolian Saltation, Science, 241, 820–823, 1988. a

Bou-Zeid, E., Meneveau, C., and Parlange, M. B.: Large-eddy simulation of neutral atmospheric boundary layer flow over heterogeneous surfaces: Blending height and effective surface roughness, Water Resour. Res., 40, https://doi.org/10.1029/2003WR002475, 2004. a

Calaf, M., Meneveau, C., and Meyers, J.: Large eddy simulation study of fully developed wind-turbine array boundary layers, Phys. Fluids, 22, 015110, https://doi.org/10.1063/1.3291077, 2010. a

Comola, F. and Lehning, M.: Energy- and momentum-conserving model of splash entrainment in sand and snow saltation, Geophys. Res. Lett., 44, 1601–1609, https://doi.org/10.1002/2016GL071822, 2017. a, b

Dai, X. and Huang, N.: Numerical simulation of drifting snow sublimation in the saltation layer, Sci. Rep.-UK, 4, 6611, https://doi.org/10.1038/srep06611, 2014. a

Déry, S. J. and Yau, M. K.: Large-scale mass balance effects of blowing snow and surface sublimation, J. Geophys. Res.-Atmos., 107, ACL 8-1–ACL 8-17, https://doi.org/10.1029/2001JD001251, 2002. a, b

Diebold, M., Higgins, C., Fang, J., Bechmann, A., and Parlange, M. B.: Flow over hills: a large-eddy simulation of the Bolund case, Bound.-Lay. Meteorol., 148, 177–194, 2013. a

Doorschot, J. J. J. and Lehning, M.: Equilibrium Saltation: Mass Fluxes, Aerodynamic Entrainment, and Dependence on Grain Properties, Bound.-Lay. Meteorol., 104, 111–130, https://doi.org/10.1023/A:1015516420286, 2002. a

Dover, S.: Numerical Modelling of Blowing Snow, PhD thesis, University of Leeds, 1993. a, b

Giometto, M., Christen, A., Meneveau, C., Fang, J., Krafczyk, M., and Parlange, M.: Spatial characteristics of roughness sublayer mean flow and turbulence over a realistic urban surface, Bound.-Lay. Meteorol., 160, 425–452, 2016. a

Grazioli, J., Madeleine, J.-B., Gallée, H., Forbes, R. M., Genthon, C., Krinner, G., and Berne, A.: Katabatic winds diminish precipitation contribution to the Antarctic ice mass balance, P. Natl. Acad. Sci. USA, 114, 10858–10863, https://doi.org/10.1073/pnas.1707633114, 2017. a

Groot Zwaaftink, C. D., Löwe, H., Mott, R., Bavay, M., and Lehning, M.: Drifting snow sublimation: A high-resolution 3-D model with temperature and moisture feedbacks, J. Geophys. Res.-Atmos., 116, D16107, https://doi.org/10.1029/2011JD015754, 2011. a, b, c

Groot Zwaaftink, C. D., Diebold, M., Horender, S., Overney, J., Lieberherr, G., Parlange, M. B., and Lehning, M.: Modelling Small-Scale Drifting Snow with a Lagrangian Stochastic Model Based on Large-Eddy Simulations, Bound.-Lay. Meteorol., 153, 117–139, https://doi.org/10.1007/s10546-014-9934-2, 2014. a

Helgans, B. and Richter, D. H.: Turbulent latent and sensible heat flux in the presence of evaporative droplets, Int. J. Multiphas. Flow, 78, 1–11, https://doi.org/10.1016/j.ijmultiphaseflow.2015.09.010, 2016. a

Kok, J. F. and Renno, N. O.: A comprehensive numerical model of steady state saltation (COMSALT), J. Geophys. Res.-Atmos., 114, D17204, https://doi.org/10.1029/2009JD011702, 2009. a, b

Kumar, V., Kleissl, J., Meneveau, C., and Parlange, M. B.: Large-eddy simulation of a diurnal cycle of the atmospheric boundary layer: Atmospheric stability and scaling issues, Water Resour. Res., 42, https://doi.org/10.1029/2005WR004651, 2006. a

Lenaerts, J. T. M., van den Broeke, M. R., Déry, S. J., van Meijgaard, E., van de Berg, W. J., Palm, S. P., and Sanz Rodrigo, J.: Modeling drifting snow in Antarctica with a regional climate model: 1. Methods and model evaluation, J. Geophys. Res.-Atmos., 117, D5, https://doi.org/10.1029/2011JD016145, 2012. a, b

Liston, G. E. and Sturm, M.: The role of winter sublimation in the Arctic moisture budget, Nord. Hydrol., 35, 325–334, https://doi.org/10.1029/2001JD001251, 2004. a

Mann, G. W., Anderson, P. S., and Mobbs, S. D.: Profile measurements of blowing snow at Halley, Antarctica, J. Geophys. Res.-Atmos., 105, 24491–24508, https://doi.org/10.1029/2000JD900247, 2000. a, b

Paterna, E., Crivelli, P., and Lehning, M.: Decoupling of mass flux and turbulent wind fluctuations in drifting snow, Geophys. Res. Lett., 43, 4441–4447, https://doi.org/10.1002/2016GL068171, 2016. a

Raupach, M.: Saltation layers, vegetation canopies and roughness lengths, Acta Mech., I, 83–96, 1991. a

Rémy, F. and Frezzotti, M.: Antarctica ice sheet mass balance, C. R. Geosci., 338, 1084–1097, https://doi.org/10.1016/j.crte.2006.05.009, 2006. a

Rignot, E., Velicogna, I., van den Broeke, M. R., Monaghan, A., and Lenaerts, J. T. M.: Acceleration of the contribution of the Greenland and Antarctic ice sheets to sea level rise, Geophys. Res. Lett., 38, L05503, https://doi.org/10.1029/2011GL046583, 2011. a

Russo, E., Kuerten, J. G. M., van der Geld, C. W. M., and Geurts, B. J.: Water droplet condensation and evaporation in turbulent channel flow, J. Fluid Mech., 749, 666–700, https://doi.org/10.1017/jfm.2014.239, 2014. a

Shao, Y. and Li, A.: Numerical Modelling of Saltation in the Atmospheric Surface Layer, Bound.-Lay. Meteorol., 91, 199–225, https://doi.org/10.1023/A:1001816013475, 1999. a

Sharma, V., Parlange, M., and Calaf, M.: Perturbations to the spatial and temporal characteristics of the diurnally-varying atmospheric boundary layer due to an extensive wind farm, Bound.-Lay. Meteorol., 162, 255–282, 2017. a

Sturm, M., Holmgren, J., and Perovich, D. K.: Winter snow cover on the sea ice of the Arctic Ocean at the Surface Heat Budget of the Arctic Ocean (SHEBA): Temporal evolution and spatial variability, J. Geophys. Res.-Oceans, 107, SHE 23-1–SHE 23-17, https://doi.org/10.1029/2000JC000400, 2002. a

Thorpe, A. D. and Mason, B. J.: The evaporation of ice spheres and ice crystals, Brit. J. Appl. Phys., 17, 541, https://doi.org/10.1088/0508-3443/17/4/316, 1966. a, b, c, d, e

van den Broeke, M., van de Berg, W. J., van Meijgaard, E., and Reijmer, C.: Identification of Antarctic ablation areas using a regional atmospheric climate model, J. Geophys. Res., 111, D18110, https://doi.org/10.1029/2006JD007127, 2006. a, b

Vionnet, V., Martin, E., Masson, V., Guyomarc'h, G., Naaim-Bouvet, F., Prokop, A., Durand, Y., and Lac, C.: Simulation of wind-induced snow transport and sublimation in alpine terrain using a fully coupled snowpack/atmosphere model, The Cryosphere, 8, 395–415, https://doi.org/10.5194/tc-8-395-2014, 2014. a, b, c, d

Vögeli, C., Lehning, M., Wever, N., and Bavay, M.: Scaling Precipitation Input to Spatially Distributed Hydrological Models by Measured Snow Distribution, Front. Earth Sci., 4, 108, https://doi.org/10.3389/feart.2016.00108, 2016. a

Walter, B., Horender, S., Voegeli, C., and Lehning, M.: Experimental assessment of Owen's second hypothesis on surface shear stress induced by a fluid during sediment saltation, Geophys. Res. Lett., 41, 6298–6305, https://doi.org/10.1002/2014GL061069, 2014. a

Wever, N., Lehning, M., Clifton, A., Rüedi, J. D., Nishimura, K., Nemoto, M., Yamaguchi, S., and Sato, A.: Verification of moisture budgets during drifting snow conditions in a cold wind tunnel, Water Resour. Res., 45, 1–14, https://doi.org/10.1029/2008WR007522, 2009. a

Xiao, J., Bintanja, R., Déry, S. J., Mann, G. W., and Taylor, P. A.: An Intercomparison Among Four Models Of Blowing Snow, Bound.-Lay. Meteorol., 97, 109–135, https://doi.org/10.1023/A:1002795531073, 2000. a