On the suitability of the Thorpe-Mason model for Calculating Sublimation of Saltating Snow

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 (LES), it is found that the 5 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 air and surface temperatures, these errors range from 26% for low-wind low-saturation conditions to 38% for high-wind high-saturation conditions. With a small temperature perturbation of 1 K between the air and the snow surface, the errors due to the TM model are already as high as 100% with errors increasing for larger temperature perturbations. 10


Introduction
Sublimation of drifting and blowing snow has been recognized as an important component of the mass budget of polar and alpine regions (Liston and Sturm, 2004;van den Broeke et al., 2006;Lenaerts et al., 2012;Vionnet et al., 2014). Field observations and modelling efforts focused on Antarctica have highlighted the fact that precipitation and sublimation losses are the 15 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 Aeolian transport of snow can be classified into two modes, namely, saltation and suspension. Saltation consists of particles being transported along the surface via short, ballistic trajectories with heights less than 10 cm and involves mechanisms of aerodynamic entrainment along with rebound and splashing of ice grains (Doorschot and Lehning, 2002;Comola and Lehning, 2017). 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.

5
This is true for both field studies, where sublimation losses are calculated using measurements, usually at the height of O(1 m), and in mesoscale modelling studies, where the computational grids and time-steps are too large to resolve flow dynamics at saltation length and time scales. 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 large-eddy simulations (Dai and Huang, 2014) show that sublimation losses in the saltation layer are not negligible, particularly for wind speeds close to the threshold 10 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, K is the thermal conductivity of moist air and D is the mass diffusivity of water vapor in air. Transfer of heat and mass is driven by differences of temperature and vapor density between the particle 15 surface (T p , ρ w,p ) and the surrounding fluid (T a,∞ , ρ w,∞ ) and enhanced by turbulence, the effect of which is parametrized by the Nusselt (Nu) and Sherwood (Sh) numbers respectively. Nu and Sh are related to the relative speed (|u rel |) between the air and the particle via the particle Reynolds number (Re p ) as Re p = d |u rel | ν air ; Nu = 1.79 + 0.606 Re 1/2 p P r 1/3 ; Sh = 1.79 + 0.606 Re 1/2 p Sc 1/3 , where ν air is the kinematic viscosity of air and P r and Sc are the Prandtl and Schmidt numbers respectively. Thorpe and Mason (1966) solved the above coupled Eqs. (1) and (2) by, (a) neglecting the thermal inertia of the ice particle, 20 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 vapor density of air surrounding the particle, saturation σ * = ρ w,∞ /ρ s (T a,∞ ), M is the molecular weight of water and R is the universal gas constant. The above formulation has been used extensively to analyse data 25 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., 2011Vionnet et al., 2014).In the modelling studies, the mass loss term is computed using Eq. (4) and is added, with proper normalisation, 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 5 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 ξ = 6ρ air c p,air πρ i c i d p 3 N , where d p is the mean particle diameter, N is the particle number density, ρ 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, d p = 50 µm and Secondly, Eq. (4) computes mass loss as being directly proportional to σ * and neglects the temperature difference between 15 the particle and air. Eq. (4) thus predicts a mass loss even in extremely high saturation conditions, whereas immediate deposition of water vapor 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 under-saturation 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 20 with moderate humidity! Motivated by the observations described above, in this letter, we describe four numerical experiments where we compare differences between the fully numerical and the Thorpe and Mason (1966) solutions (referred to as NUM and TM approaches respectively). In Experiment 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 25 (LES) of 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 letter is made in Sect. 4

Comparison between NUM and TM solutions: EXPERIMENT I and II
We consider an idealised scenario where a spherical ice particle is held still in a turbulent air flow with constant mean speed, 30 temperature and under-saturation. 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. 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 of the heat and mass evolution of a ice particle during sublimation and about the TM model.
For the NUM approach, Eqs. (1) and (2) are solved using a simple first-order finite-differencing scheme. The initial particle diameter is 200 µm and the air-flow temperature is 263.15 K. We use a constant air speed of 5 ms −1 resulting in Re p = 80, Nu = 6.7 and Sh = 6.5 (using Eq. (3)). The values used here are typical of a saltating ice particle. In Experiment I, we study 5 the heat and mass output from a sublimating ice grain as a function of time. In the first case, we consider the effect of three different values of air-flow saturation (σ * = 0.8, 0.9 and 0.95) on differences between NUM and TM solutions. In this case, the initial particle temperature is taken to be the same as the air-flow temperature for the NUM approach. Results are shown in for mass, Err M and heat, Err Q are shown. The errors reduce dramatically with time (for example, 15% at 0.3 seconds) and do not depend on the saturation of the air-flow. In the following case, similar simulations are performed with σ * = 0.95 while the initial temperature difference between the particle and the air is varied as T P − T Air = −2, −1, 1, 2 K. It is interesting to note that while the TM solution predicts sublimation of the 15 particle (consistent with σ > 0), 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 gain 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). 20 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 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 |u rel | due to more vigorous heat and mass transfer. Experiment I was repeated for values of d p and |u rel | ranging between (50 − 1000 µm) and 0 − 10 ms −1 respectively. The upper-bound of the 25 wind-speed range is quite high and it is extremely unlikely to find |u rel | > 10 ms −1 in naturally-occurring aeolian transport.
Numerical results indeed confirm our intuition and it is found that for any given value of |u rel |, τ relaxation is found to be ∝ d α p , where α (∼ 1.65). Furthermore, τ relaxation decreases monotonically with increasing |u rel | 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 seconds (for |u rel | = 10 and 0 ms −1 respectively). Plots of τ relaxation are highly relevant to discussion in Sect. 3 and presented there.

30
Following results of Experiment I, in Experiment II, we explore the parameter space of (σ * , T P − T Air ) and compute the while the TM solutions do not depend on T P −T Air as expected. The error between the NUM and TM solutions are accentuated at high saturation regimes, with errors larger than 30 % for σ * > 0.8.
In summary, Experiments I and II highlight the fact that during the sublimation of an ice grain, there exists a finite, welldefined 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 5 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.  (Raupach, 1991;Shao and Li, 1999). The LBC for scalars are flux-free and thus the only source/sink of heat and water vapour 20 in the simulations is through the interaction of the flow with the saltating particles. The snow surface consists of 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. All simulations are performed on a grid 25 of 64 x 64 x 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 setup are provided in Supplementary Material S2.
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 30 the turbulent transfer of heat between the air and the particle (second term in R.H.S 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 (PRT) 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 realistic saltation of snow. PRT is defined as the total time the particle is air-borne and in motion, including multiple hops across the surface. Towards this goal, simulations are performed, each with a combination of initial surface stress, 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 10 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 * = 0.4) and choose the initial saturation to be 0.95, motivated by results in Experiment II where errors were found to increase with increasing saturation.
Simulations (named as UL-T(γ), where T Air,IC − T P,IC = γ) are performed once again for each of two approaches with γ ∈ {±1 K , ±2.5 K , ±5 K} resulting in a total of twelve simulations. In all simulations performed for Experiments III and 15 IV, T P,IC = 263.15 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.

20
As mentioned in the concluding lines of the Sect. 2, the quantity of interest is the PRT of saltating ice grains. In Fig. 2a, 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 ms −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 3 wind speeds, 0, 1 and 10 ms −1 are shown and the power-law dependence can 25 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 seconds while the median PRT is 0.2 seconds, 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 30 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. 2a provides 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 wind speeds. This means that the dynamics of the heavier particles are unaffected by larger-scale turbulence statistics, 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 increases with wind speed. Secondly, the initial saturation does not seem to effect 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. 2a). In these cases turbulence is sufficient to rapidly mix any temperature anomaly due to sublimation throughout the surface layer. On the other 5 hand, in low wind conditions (UL), low initial saturation results in more sublimation and cooling near the surface, resulting in suppression of vertical motions. This is reflected in the mean PRTs of particles smaller than 75 µm, which decrease with decreasing initial σ * .
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. 2a. These distributions are also strongly dependent on the particle diameter. As an illustration, in 10 Fig. 2b, 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% to 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.

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 field scale. In Fig. 3, 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 20 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 Supplementary Fig. S1). The increase in error with increasing σ * is in accordance with analysis done in Experiment II ( see Fig. 2(i,l) ) where it was shown that the NUM and TM solutions diverge with increasing saturation. The least error, 26% is found for case UL-RL (i.e., u * = 0.4, σ * = 0.3) while the largest error, 38% 25 is found for UH-RH (u * = 0.8, σ * = 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. 3c, 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 30 approach shows a much higher sublimation rate than the TM solutions. This occurs firstly due to higher vapor pressure at the grain surface that results in enhanced vapor 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, similar to results in Experiment III and are shown here).
The LES results do come with two caveats. Firstly, the scalar fluxes at the surface are neglected. This can 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 and in our current LESs (see Supplementary Fig. S2). This implies that direct turbulent 5 exchange between the surface and air is curtailed and instead, the dominant exchange occurs between air-borne particles and the air. This is nevertheless an important assertion that shall be more closely examined in future studies involving a full surface energy balance model. Secondly, we neglect collisions between particles in the saltation layer. While the effect of collisions has been shown to be important, incorporating their effects in the present work is not likely to affect our primary conclusion of the inadequacy of the TM approach for saltating snow.

Discussion and Conclusion
In this letter, 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 as τ relaxation that has a power-law dependence on the particle diameter and is inversely proportional to the relative wind speed. Through extensive 15 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 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 conditions to 38% in high-wind, high saturation 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 20 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 tends towards 1.
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. Future efforts are 25 geared towards including a high-fidelity snowpack thermodynamic model within the LES-saltation framework to further study sublimation of drifting and blowing snow in a more realistic setting and to compare simulations with field observations.