Warm-air entrainment and advection during alpine blowing snow events

Blowing snow transport has considerable impact on the hydrological cycle in alpine regions both through the redistribution of the seasonal snowpack and through sublimation back into the atmosphere. Alpine energy and mass balances are typically modeled with time-averaged approximations of sensible and latent heat fluxes. This oversimplifies nonstationary turbulent mixing in complex terrain and may overlook important exchange processes for hydrometeorological prediction. To determine if specific turbulent motions are responsible for warm- and dry-air advection during blowing snow events, quadrant analysis and variable interval time averaging was used to investigate turbulent time series from the Fortress Mountain Snow Laboratory alpine study site in the Canadian Rockies, Alberta, Canada, during the winter of 2015–2016. By analyzing wind velocity and sonic temperature time series with concurrent blowing snow, such turbulent motions were found to supply substantial sensible heat to near-surface wind flows. These motions were responsible for temperature fluctuations of up to 1 ^∘C, a considerable change for energy balance estimation. A simple scaling relationship was derived that related the frequency of dominant downdraft and updraft events to their duration and local variance. This allows for the first parameterization of entrained or advected energy for time-averaged representations of blowing snow sublimation and suggests that advection can strongly reduce thermodynamic feedbacks between blowing snow sublimation and the near-surface atmosphere. The downdraft and updraft scaling relationship described herein provides a significant step towards a more physically based blowing snow sublimation model with more realistic mixing of atmospheric heat. Additionally, calculations of return frequencies and event durations provide a field-measurement context for recent findings of nonstationarity impacts on sublimation rates.

1 Introduction

At least 40 % of the world's population relies on the seasonal snowpack as a temporary reservoir of winter snowfall that then provides meltwater in spring and summer for downstream water use (Meehl et al., 2007). However, after snow has fallen, it is often subjected to sublimation while at rest or amplified in-transit sublimation during redistribution. Blowing snow redistribution can result in vast amounts of frozen water moving between basins or, in the case of sublimation, being removed entirely from the surface water budget in windswept regions. Blowing snow particles are highly susceptible to sublimation because of their high curvature, large surface-area-to-mass ratio, and high ventilation rates (Dyunin, 1959; Schmidt, 1982). While estimates may vary with climate, in the Canadian Rockies, blowing snow transport has been found to be responsible for sublimating 17 %–19 % of the yearly snowfall (MacDonald et al., 2010).

Snow sublimation is typically studied at large temporal and spatial scales within hydrometeorological modeling frameworks because of the complexity of the processes and the difficulty of particle transport tracking (e.g. Pomeroy et al., 1993; Pomeroy and Essery, 1999; Déry and Yau, 2002; Lenaerts et al., 2012; Groot Zwaaftink et al., 2013; Musselman et al., 2015). To accurately calculate all contributions to boundary layer energy balances, latent-heat-flux estimates rely on an accurate sublimation model and a precise understanding of how much energy is available for snow particle phase change. While latent and sensible heat exchanges between the turbulent atmosphere and snow particles can be represented by a system of coupled partial differential equations, they require forcing terms for dry-air entrainment and horizontal advection that are still poorly understood and have not been based on observed physical mechanisms if they are included at all (Bintanja, 2001). The decrease in temperature and increase in humidity in the atmosphere caused by snow sublimation may play a crucial limiting role in snow sublimation, but many blowing snow models struggle to capture the process of this feedback, which can result in unrealistic atmospheric conditions in near-surface boundary layers and subsequent errors in calculations of the blowing snow sublimation rate (Pomeroy and Li, 2000; Déry and Yau, 1999, 2001; Groot Zwaaftink et al., 2013; Musselman et al., 2015).

Numerical investigations of snow sublimation from a numerical modeling approach have recently provided new insights into non-steady-state aspects of sublimation (Dai and Huang, 2014; Li et al., 2017; Sharma et al., 2018) and the efficacy of the nearly universally used Thorpe and Mason (1966) model (see Schmidt, 1972) at high temporal and spatial resolutions (Sharma et al., 2018). Extending nonstationary sublimation models to alpine and other complex terrain environments could lead to reduced uncertainty in blowing snow sublimation models. Little research has been conducted to better understand the energy available for snow sublimation from entrainment or advection processes in natural atmospheric turbulence or the influence that resultant air temperature fluctuations may exert on sublimation rates. Recently, Grazioli et al. (2017) found that over long timescales, persistent katabatic winds in Antarctica provide a consistent supply of unsaturated air that can contribute to significant snow sublimation. In East Antarctica, they calculated up to 35 % of total yearly snowfall can be lost in this manner.

The objective of this research is to investigate turbulent structures down to subsecond timescales and identify their synchronization with near-surface temperature fluctuations. The study investigates the unsteady processes affecting blowing snow particle energy balances to better understand the form of advection and entrainment correction terms for sublimation calculations. To this end, a scaling relationship previously applied to near-neutral atmospheric-surface-layer data is tested to represent turbulent event frequency as a function of variable interval time averaging (VITA) parameters. Data used to validate this model consist of field measurements of three-dimensional wind velocities and sonic temperatures during blowing snow events at the blowing snow study site in Fortress Mountain Snow Laboratory (FMSL), Canadian Rockies (Fig. 1). These data are supplemented by observations of nearby temperature, relative humidity (RH), and wind speeds at three additional meteorological stations within FMSL. This provides a broader environmental context in which to understand potential thermodynamic feedback mechanisms beyond the blowing snow study site scale. The scaling relationship also gives a real-world context for recent numerical studies on the impacts of nonstationarity on blowing snow sublimation rates.

2 Methods

2.1 Field data

Fortress Mountain Snow Laboratory (FMSL) is located in the Kananaskis Valley in the Canadian Rockies of southwestern Alberta, Canada. FMSL is surrounded by very complex terrain, with multiple nearby 2900 m peaks having > 100 m vertical rock faces. The blowing snow study site is situated on a plateau at 2000 m at the base of a closed ski resort, providing ample upwind fetch with minimal obstruction from trees or buildings (Fig. 1 inset). Winter air temperatures at the FMSL blowing snow site typically range from −20 to +5 ^∘C, with frequent midwinter downslope Chinook (föhn) wind events. Snow depth at the blowing snow site remains fairly constant through the midwinter at approximately 1 m, with fresh snowfall frequently redistributed by wind events.

Figure 1 Overview of blowing snow study site and adjacent micrometeorological stations at Fortress Mountain Snow Laboratory in the Kananaskis Valley, Alberta, Canada. The prominent down-valley wind direction is noted on the map and on the site photo inset. Topographic map produced by the Canada Centre for Mapping and Earth Observation, Natural Resources Canada, © Her Majesty The Queen in Right of Canada.

Ultrasonic temperature and wind velocity time series were observed at the FMSL blowing snow study site using two Campbell Scientific CSAT3 sonic anemometers sampling at 50 Hz from November 2015 to March 2016. The anemometers were positioned on the same mast at heights above the snow surface varying over 0.1–0.4 and 1.4–1.7 m with snow surface accumulation or erosion. These anemometers will be referred to as the 20 and 150 cm anemometers throughout the remaining text. Extensive analysis of this dataset has already provided new insights into the turbulent mechanisms for blowing snow transport (Aksamit and Pomeroy, 2016, 2017, 2018). The turbulent structures scrutinized here have previously been coupled with particle tracking velocimetry (PTV) and high-speed video analysis of Aksamit and Pomeroy (2017, 2018) to better understand the wind–snow coupling. For each night (20 November, 4 December 2015 and 3 February, 3 March 2016), the time series spanning the entire duration of blowing snow video recording (from 18:00 local time to the end of video collection, approximately 23:59 local time) was divided into 15 min intervals and analyzed. One additional night of meteorological data was analyzed to compare energy transport mechanisms under much windier conditions, even though PTV analysis was not available. This additional night, 21 January 2016, had much stronger winds, gusting up to 15 m s⁻¹ because of the presence of a Chinook event. In this valley, these events have been previously associated with high blowing snow sublimation rates (MacDonald et al., 2018). The mean temperatures varied from −7 ^∘C during the previous 3 d to +3 ^∘C during the night of 21 January. This resulted in a much larger difference between air and snow surface temperatures and provided an interesting comparison of conditions that are critical for snow sublimation at short timescales (Sharma et al., 2018).

Three other FMSL stations near to the blowing snow measurement site provide complementary 15 min relative humidity (with respect to ice), air temperature, and wind speed measurements (Fig. 1). As relative humidity measurements were not available at the blowing snow study site during the 2015–2016 study season, these additional stations provided downwind test sites for evidence of the occurrence of large-scale thermodynamic feedbacks. The nearest complementary site is a sheltered forest station (Powerline) approximately 400 m away and 30 m higher in elevation (Smith et al., 2017). Additionally, there are two exposed sites, a ridgetop (Canadian Ridge) and lee side of a ridge (Canadian Ridge North) that are both approximately 600 m downwind and 200 m higher in elevation. The Powerline station receives much less wind than the exposed sites or blowing snow site and is much less susceptible to snow redistribution. Temperature and relative humidity from a previous study spanning January to March 2015 at the blowing snow study site showed a good correlation with the three nearby sites over a sample of 3700 15 min averages. R² values for air temperature between the blowing snow site and Canadian Ridge, Canadian Ridge North, and Powerline were 0.82, 0.83, and 0.97, respectively, and were 0.61, 0.62, and 0.80 for relative humidity, respectively. All correlation coefficients were statistically significant at 99.99 %. Meteorological variables at the blowing snow study site can be found in Table 1. The Monin–Obukhov stability parameter, ζ, was calculated following Monin (1970) and Stull (1988) such that

$$\begin{array}{rl}& \mathit{\zeta }={\displaystyle \frac{z}{L}},\\ & L=-{\displaystyle \frac{u_{\ast }^{\mathrm{3}}{\mathit{\theta }}_{\mathrm{0}}}{\mathit{\kappa }g\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}^{\prime }}}}\end{array},$$

where θ₀ is the potential temperature at the 20 cm anemometer and $u_{\ast }={\left({\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{4}}$. We use the sonic temperature in lieu of potential temperature as suggested by Stull (1988) as there were no atmospheric pressure measurements at the study site. Turbulence intensity was calculated as

$$I={\displaystyle \frac{\stackrel{\mathrm{‾}}{{u^{\prime }}^{\mathrm{2}}+{v^{\prime }}^{\mathrm{2}}+{w^{\prime }}^{\mathrm{2}}}}{{\left({\stackrel{\mathrm{‾}}{u}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{w}}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}}}.$$

Table 1 Meteorological Variables for five nights of observations. ^* Snow surface temperature taken from the nearby Powerline meteorological station.

2.2 Modified VITA analysis

Variable interval time averaging (VITA) is a method of time series analysis that identifies significant turbulent events as periods of high local variance. For a given time series f(t), the VITA method selects times such that

$$\begin{array}{}\text{(1)}& \widehat{f}(t,T)={\displaystyle \frac{\mathrm{1}}{T}}\underset{t-T/\mathrm{2}}{\overset{t+T/\mathrm{2}}{\int }}f(t{)}^{\mathrm{2}}\text{d}t-{\left[{\displaystyle \frac{\mathrm{1}}{T}}\underset{t-\frac{T}{\mathrm{2}}}{\overset{t+\frac{T}{\mathrm{2}}}{\int }}f\left(t\right)\text{d}t\right]}^{\mathrm{2}}>k_V\stackrel{\mathrm{‾}}{f^{\mathrm{2}}},\end{array}$$

where T is a statistically or experimentally determined averaging time, k_V is a user-defined VITA threshold, and the overbar indicates a spatial or temporal average. For the analysis conducted here, as in a previous blowing snow investigation at this site (Aksamit and Pomeroy, 2017), f(t) is taken to be the instantaneous Reynolds stress $\mathit{\tau }\left(t\right)=-{\mathit{\rho }}_{\mathrm{air}}{u}^{\prime }\left(t\right)w^{\prime }\left(t\right)$, where u and w are the instantaneous fluctuations around 15 min averages of streamwise and vertical velocities. Equation (1) needs two user-defined parameters, the temporal neighborhood T and the magnitude threshold k_v. To increase objectivity and connect our turbulent events to extensively studied and physical turbulent structures, a modified VITA analysis used here also includes a quadrant hole analysis criterion (Lu and Willmarth, 1973; Morrison et al., 1989). Subsequent to finding a time meeting the conditions defined by Eq. (1), we then identify the neighborhood surrounding that time where τ also exceeds a given threshold, k_Q, often called the “quadrant hole” value:

$$\begin{array}{}\text{(2)}& \left|\mathit{\tau }\left(t\right)\right|\ge k_Q{\mathit{\rho }}_{\mathrm{air}}\sqrt{\stackrel{\mathrm{‾}}{{u^{\prime }}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{{w^{\prime }}^{\mathrm{2}}}}.\end{array}$$

The modified VITA analysis was conducted over a range of thresholds ($\mathrm{0.01}\le k_V\le \mathrm{0.05}$, $\mathrm{0.05}\le k_Q\le \mathrm{4}$) and averaging times ($\mathrm{0.5}\le T\le \mathrm{40}$ s) found in the boundary layer and sediment transport literature (Morrison et al., 1989; Narasimha and Kailas, 1987, 1990; Bauer et al., 1998; Sterk et al., 1998; Wiggs and Weaver, 2012). This provides a relatively robust gust identification scheme that delimits significant turbulent events of varying duration and velocity magnitude. The subsequent analysis focuses on sweeps ($u^{\prime }>\mathrm{0}$, $w^{\prime }<\mathrm{0}$) and ejections ($u^{\prime }<\mathrm{0}$, $w^{\prime }>\mathrm{0}$) as they disproportionately contribute to the total surface Reynolds stress, are frequently used in models of turbulent boundary layer structures, and have been identified to play crucial roles in boundary layer heat flux and aeolian transport (e.g., Bauer et al., 1998; Adrian et al., 2000; Garai and Kleissel, 2011; Aksamit and Pomeroy, 2017). Please refer to Wallace (2016) for a recent review of the theory and experiments surrounding quadrant analysis and sweep–ejection cycling. The modified VITA algorithm categorized a turbulent event as a sweep or ejection if the parameterized curve $s\left(t\right)=\langle u^{\prime }\left(t\right),w^{\prime }\left(t\right)\rangle$ passes through only one of the two quadrants during the event (Q2 for ejections and Q4 for sweeps).

In this study, the concurrent sonic temperature signal response was also measured and the fluctuation from the 15 min mean air temperature was computed to identify the presence of relatively warmer or colder air during a particular event with respect to mean conditions. For the air temperatures during this study, CSAT3 anemometers have an error of less than ±0.002 ^∘C, which is considered negligible (Michael Clarke, Campbell Scientific, personal communication, 9 May 2018). Following Kailas and Narasimha (1994), events detected with larger thresholds are referred to as “stronger” or “more intense”.

Figure 2 Temperature and relative humidity 15 min average measurements at approximately 2 m above the snow surface during the five nights of investigation at three nearby micrometeorological stations. Flagged data have been removed from the time series and presented as gaps. Shaded gray areas are times when concurrent average wind speeds were above 3 m s⁻¹. Note the correlation between sites for both variables and the varying y axes between plots.

3 Results

3.1 Modified VITA results

During each blowing snow storm, there was no definitive evidence of humidity saturation or thermodynamic feedback at two of the three nearby weather stations (Fig. 2). Unfortunately, RH data were unusable at the Canadian Ridge site for four of the five nights. Increases in RH were typically coupled with decreases in air temperature and were transient in nature. The complex topography and enhanced turbulent mixing at FMSL may be responsible for this as indicated by the modified VITA analysis below. Indeed, though all three sites are situated in close proximity to each other, there is limited correlation for meteorological variables between all three, suggesting incredibly complex wind flow and energy fluxes in this alpine zone.

Figure 3 Bin-averaged temperature fluctuations of near-surface anemometer from the 15 min mean for events of specific return frequency and event duration for recordings over each blowing snow storm. Insets are plots of probability distribution functions of event duration, temperature deviation, and event frequency for each storm and type of event.

For each VITA-identified event, instantaneous temperature deviations from the 15 min mean were computed to represent the magnitude and sign of turbulent temperature mixing with respect to slower meteorological changes over the nights. Aksamit and Pomeroy (2017) noted that there is no objective choice of averaging time or event threshold for the modified VITA analysis. As such, sonic temperatures during active turbulent events were examined over a variety of thresholds to determine a range of behavior. The recurrence frequencies and average durations of sweep and ejection events for each threshold combination illustrated the average prevalence of sweep or ejection motions. Further sensitivity analysis of the impact of VITA parameters on wind–snow coupling has been conducted by Aksamit and Pomeroy (2017).

For each blowing snow storm, 3D point clouds of mean recurrence frequency, event duration, and event temperature deviation were calculated for the 20 cm sonic anemometer. Each point represents the values from one choice of averaging time and modified VITA thresholds for a 15 min observation period as discussed in Sect. 2.2. The 3D plots contain significant overlap, so for clarity, the mean temperature deviations were averaged over small ranges of event duration and frequency, as shown in Fig. 3. The insets in each subplot are three probability distributions computed from the original point clouds for each blowing snow storm: distributions of temperature deviation, event duration, and event frequency. Mean and skewness values are noted next to each distribution.

The analysis revealed for the four non-Chinook blowing snow storms (20 November 2015, 4 December 2015, 3 February 2016, 3 March 2016) that sweeps consistently brought warmer air to the near-surface anemometer. This can be seen as the colored temperature plots show average temperature deviations greater than zero for nearly all event duration and frequencies over each storm. Probability distributions show very few sweep events with negative temperature deviations, as well as a consistent positive mean and skewness. The Chinook storm on 21 January 2016 had a positive mean and skewness but exhibited short cold-air bursts as well. Mean temperatures for sweeps were warmer than ejections for all blowing snow storms.

Figure 4 The same measurements as shown in Fig. 3, now subtracting the upper-anemometer (140 cm above the surface) mean from the near-surface (20 cm) sweep and ejection temperatures for various blowing snow storms. Note the predominantly colder fluctuations.

Modified VITA analysis on the 140 cm anemometer wind and temperature time series showed similar results. Interestingly, analysis comparing 20 cm anemometer event temperatures to the 140 cm anemometer means revealed near-surface sweeps often occur with colder signatures than the nearby 140 cm anemometer means (Fig. 4). This is in contrast to what was found relative to surface temperatures. As these measurements were all made during the night and over a continuous snow cover with a slightly stable temperature profile, this indicates relatively warm upper-air mixing with cold near-surface air that resulted in a mixed temperature value between the two anemometer means. For example, blocks on the left side of Fig. 3i show a group of sweeps that were 1 ^∘C warmer than the mean temperature of the 20 cm anemometer (bright yellow), but in Fig. 4i, the same group of sweeps were 0.5 ^∘C colder than the 140 cm anemometer mean (dark blue). Color scales are equivalent in Figs. 3 and 4. This effect is further supported by the mean anemometer temperatures detailed in Table 1.

Figure 5 Example of variation in event frequency (N) over VITA thresholds (k_V) for different averaging times ($\mathrm{0.5}<T<\mathrm{12}$ s) at k_Q=1 for one 15 min study period on 4 December. Least-squares-fitted Eq. (3) curves are overlaid as dashed line for each collection of events: All modified VITA events, sweeps, and ejections as identified at the 20 cm (low) or 140 cm (high) anemometers.

Ejections present a less clear story but also appear very effective for surface layer mixing. Depending on the intensity of the events detected, ejections could be either warmer or colder than the 20 cm anemometer mean (Fig. 3). Over all blowing snow storms, ejection temperatures had a lower mean and were less positively skewed. This can be explained by their physical definition of moving air vertically away from the cold snow surface. During periods of greater atmospheric stability (20 November and 3 March) there was more variability in the temperature contributions from ejections. This may indicate stable layers of varying strength were able to form and cause less uniform mixing near the snow surface. When Monin–Obukhov coefficients were closer to zero (Table 1), indicating more neutral conditions, there was less variability in ejection temperatures, indicating a smaller range of temperatures during ejection-induced mixing. This can be seen by a comparison of Table 1 values and Figs. 3 and 4 probability insets. This mixing process is discussed in more detail in Sect. 4.

Over all nights, sweeps were of longer duration than ejections and had a higher frequency of occurrence. The probability curves in Fig. 3 show a second sweep return frequency peak of around 0.5 Hz for all nights. This is not present in the ejection frequency probabilities, which only have a single low frequency peak. These sweep and ejection motions have not been connected to a specific flow topology in these experiments (e.g., a hairpin bursting process) due to the complexity of the flow in this complex terrain. It very well may be the case that the sweep signatures are caused by both outer-layer and inner-layer motions as previously suggested by Aksamit and Pomeroy (2017). The ejections occur less often because of the rarity of large positive w values close to the snow surface and are thus present only during a less common generating mechanism.

3.2 Scaling relation

Though several differences in the datasets exist, the near-neutral and slightly stable conditions found during the blowing snow storms sampled suggest a Kailas and Narasimha (1994) scaling relationship may exist:

$$\begin{array}{}\text{(3)}& N=N_{\mathrm{0}}{e}^{-\mathit{\alpha }(k_V-\mathrm{1})}.\end{array}$$

Here, N is the recurrence frequency of a given modified VITA turbulence event type and N₀ and α are fitting parameters with N₀ known as the characteristic frequency. This scaling analysis focused on the case where k_Q=1 as this resulted in a good compromise between too many and too few events detected and is a standard value previously used for turbulent motion identification at this site (Aksamit and Pomeroy, 2017). Though the present modified VITA analysis involves an additional step in the identification algorithm as compared to the original work of Kailas and Narasimha (1994), a similar invariance (small standard deviation) in the log of the return frequency, log (N), was present over varying averaging times T for each VITA threshold k_V. This resulted in a good fit of Eq. (3) for the return frequencies of the total number of modified VITA events, as well as for sweeps and ejections individually. One example of this fitting for one 15 min period on 4 December is shown in Fig. 5. The squared l² norm of the residuals for each minimized least-squares fit are presented in the Supplement, as are the characteristic frequencies, N₀.

Table 2 Mean characteristic frequencies N₀ for turbulent events at both blowing snow site anemometers. The standard deviation of nightly means is shown in parentheses and indicates minimal changes between nights.

Total mean values (and standard deviations between nightly means) of N₀ are detailed in Table 2, for all turbulent events, only sweeps, and only ejections at both 20 and 140 cm anemometers. There was little variation in N₀ between nights of observation as seen in the relatively small standard deviation values. This suggests persistent flow features at this site from one night to the next that may be due to a persistent topographically induced flow feature or turbulence-generating mechanism at the study site. As could be expected from the analysis presented in Fig. 3, the characteristic return frequencies (N₀) of all turbulent events and for sweep events were greater than those for ejections. Of particular interest in this scaling relationship is a clear difference between N₀ for the 140 and 20 cm anemometer observations for both total events and solely sweep events. Over all nights, the characteristic frequency for total events was lower at the 140 cm anemometer, which corresponded with a drop in the number of sweeps, whereas the characteristic frequency of ejections was nearly identical at both heights. The threshold criteria in Eqs. (1) and (2) vary for measurement location and time, scaling by mean values calculated over each observation period at each anemometer. This implies that there were fewer relatively large sweeps away from the surface, and a possible shift in turbulent structure dynamics. Additionally, this supports the suggestion in Sect. 3.1 that the mechanisms generating sweeps and ejections may be different, with less common flow features resulting in the ejections.

4 Discussion

The same strong sweep events that have been previously found to be highly relevant for blowing snow initiation and transport at this site (Aksamit and Pomeroy, 2017) are also responsible for advecting warmer-than-average air to the near-surface layer. This is a critical insight for blowing snow sublimation modeling as the periods with greater than average blowing snow transport coincide with the presence of warmer-than-average air (sweeps).

Previous theoretical work has concluded that suppression of the sublimation of surface and blowing snow may occur if moisture fluxes near the surface are counterbalanced solely by diffusion (Bintanja, 2001). Mobbs and Dover (1993), Déry and Taylor (1996), Déry and Yau (1999, 2001), Groot Zwaaftink et al. (2013), and others have suggested that blowing snow sublimation could be a self-limiting process when thermodynamic feedbacks are included in a steady-state boundary layer model. However, these models did not account for warm- or dry-air entrainment or for the temporal correlation of transport bursts with warm-air entrainment. This missing forcing term may explain the lack of evidence of saturation in blowing snow field studies in the steppes of Russia, high plains of Wyoming (USA), prairies of Saskatchewan, alpine mountains of Alberta and arctic tundra of the Northwest Territories (Canada), and East Antarctica (e.g., Dyunin, 1959; Schmidt, 1982; Pomeroy and Li, 2000; Musselman et al., 2015; Grazioli et al., 2017). The evidence of frequent regeneration of warm air near the surface through advection or entrainment processes helps explain the discrepancy with diffusion-dependent models.

Figure 6 Fraction of time series occupied by sweep and ejection events of specific temperatures at the 20 cm anemometer. Refiguring of data in Fig. 3 with same color scale. Colors here also correspond to y-axis values.

Recent model simulations by Sharma et al. (2018) and Dai and Huang (2014) have shed light on the importance of temperature and wind speed fluctuations at the timescales of the sweep and ejection processes highlighted here. The comparison of the Sharma et al. (2018) large-eddy-simulation-driven sublimation model with the widely used steady-state model of Thorpe and Mason (1966) revealed that transient sublimation rates approached the steady-state model only after time periods ranging from 10⁻² to 10 s, depending on particle diameter and ventilation rates. At the velocities and particle sizes typical for the present study, their time to model relaxation was around 1 s. Furthermore, Dai and Huang (2014) found transient rates of sublimation in the saltation layer that reached a steady state only after 0.5–2 s. These modeled relaxation times are precisely in the range of turbulent warming and cooling events shown in Fig. 3. Figure 6 redisplays the data from Fig. 3 with the temporal fraction of modified VITA events calculated as the product of average event duration and frequency. For each night, one can find strong ejection events contributing air temperatures 1 ^∘C warmer than the mean for 15 % of the time and warm sweeps for up to 20 % of the time.

In addition to the timescale considerations and transient regimes in blowing snow sublimation calculations, Sharma et al. (2018) found temperature fluctuations of 1 ^∘C can affect instantaneous sublimation rates by as much as 100 %. Given that gusts causing temperature fluctuations of this order occur up to 35 % of the time, this advected energy warrants further investigation and inclusion in future models. Fortunately, a parameterization for mechanically explained advected energy may be possible through the simple exponential scaling relationships of Kailas and Narasimha (1994).

Over short timescales, there is a direct physical relationship between temperature profiles and temperature deviations during mixing events. This is physically intuitive if one considers the relative temperature change at a doorway when opening a door of a warm building to cold surroundings. Because of this dependence on instantaneous conditions during a mixing event, however, relationships between average temperature deviation magnitudes and long-term temperature gradients are not guaranteed. Comparing the nights of investigation, there is no monotonic relationship between increases in the average 140 and 20 cm sonic temperature differences and average sweep event temperature deviations. For example, on 3 March 2016 there was an average temperature difference of 0.9 ^∘C between anemometers, but the average downdraft (sweep) deviation was only +0.24 ^∘C. This is a smaller contribution than on 21 January 2016 where the air temperature difference was 0.5 ^∘C and the average sweep deviation was +0.28 ^∘C. This is almost certainly because long-time averages oversimplify the turbulent bursting process, and why eddy-covariance methods are suggested over bulk profile calculations of turbulent fluxes (Foken, 2008).

The present research has, however, suggested a simple similarity scaling of the return frequency of turbulent events of intensity k_V as identified by modified VITA analysis, through the exponential relationship of Kailas and Narasimha (1994). Such an empirical correction is compatible with the attached-eddy hypothesis (Townsend, 1976; Marusic and Monty, 2019) and other similarity scaling models of the turbulent boundary layer if the magnitude and frequency of bursts were to be defined to scale accordingly with an increase in the size of turbulent eddies away from the surface. This scaling is evident in a decrease in characteristic frequencies of turbulent events when moving from 20 to 140 cm measurements (Table 2, Supplement) and in a natural increase in modified VITA thresholds as the magnitude of turbulence measurements increases in Eqs. (1) and (2) for fixed k_V and k_Q.

This view of boundary layer mixing provides a simple platform with which to model and investigate a gust-driven regeneration function of warm–dry air at the near surface for blowing snow sublimation calculations. The inclusion of such a statistical recurrence model could provide an empirically defined quasi-periodic source of warm and dry air for blowing snow simulations. For example, this could be included in conservation of heat equations as a natural evolution of the constant entrainment and advection functions introduced by Bintanja (2001). In this way, it is possible to represent the mixing of distinct parcels of air of different temperatures through commonly studied turbulent structures. Such a recurrence model would be computationally efficient and a significant step towards a physically based blowing snow sublimation model.

Future high-temporal-resolution studies of air temperature and water vapor during sustained periods of above-snow-transport-threshold wind speeds would greatly benefit the research community. Short-timescale thermodynamic feedbacks to humidity from sublimation could come from similar high-frequency coupling analysis with closed path hygrometers or gas analyzers at multiple heights during blowing snow events. This would allow for a more complete understanding of the advection–thermodynamic feedback balance during blowing snow storms and advance the seminal profile studies of Schmidt (1982). As advection processes are local by nature (e.g., Harder et al., 2016), characteristic frequencies of turbulent events will vary with location and current atmospheric conditions. The small range of values of N₀ measured at this site during 5 months of this campaign suggests common flow phenomena will possibly dominate and aid in more universal applications of entrainment modeling, at least within specific seasons.

5 Conclusion

During an alpine blowing snow field campaign, analysis of turbulence time series and sonic temperatures indicates that exceptional warm-air entrainment and advection events can be associated with specific turbulent structures. Over five nights of investigation, sweeps brought relatively warm air to the snow surface, up to 1 ^∘C warmer than average near-surface temperatures. These parcels of air may also be relatively cold compared to temperatures measured only 1.2 m above, further adding to the complexity of the physics of blowing snow sublimation. Ejections also result in strong but less consistent temperature mixing. The current lack of understanding of advection or entrainment during snow transport may explain why the thermodynamic feedback parameterizations necessary in many blowing snow sublimation models are unphysical. In fact, field measurements of atmospheric conditions during these blowing snow events showed no evidence of significant sublimation feedbacks, let alone saturation of relative humidity. An enhanced influence of mechanical mixing in boundary layers with inhomogeneous temperature distributions, for example where there is topographically induced cold-air pooling or flow separation, may explain why sublimation rate observations and estimates can be high and can vary from study to study. The present research indicates that including a supply of warm and dry air from different near-surface regions of the flow is a physically accurate modeling assumption. A better representation of turbulent mixing in these regions is likely necessary for the improvement of sublimation rate estimates.

At present, further investigation of the connection of blowing snow sublimation to specific atmospheric structures would be beneficial. Specifically, vertical profiles of high-frequency temperature and humidity measurements are necessary to illuminate the impact of penetrating low-frequency gusts on warm- and dry-air regeneration at the surface during blowing snow sublimation in different environments. This analysis would require a closed-path-style water vapor measurement as snow particles could otherwise impact the signal quality. Such an experiment could provide high-resolution temperature and complementary water vapor measurements to more directly measure the influence of gusts on sublimation rates and begin to address discrepancies in sublimation found in different climates. Additionally, longitudinal studies of heat flux in near-surface layers would provide better insight into the connection between average boundary layer profiles and the presence of turbulent events of specific magnitude, frequency, and duration.