Dry-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 10 are typically modelled with time-averaged approximations of sensible and latent heat fluxes. This oversimplifies non-stationary turbulent mixing in complex terrain and may overlook important exchange processes for hydrometeorological prediction. To determine if warm and dry air advection during blowing snow events from atmospheric sweep and ejection motions can provide such exchange mechanisms, they were investigated at an alpine site in the Canadian Rockies and found to supply substantial sensible heat to blowing snow flows. These motions were responsible for temperature fluctuations of up to 1oC, a considerable 15 change for energy balance estimation. A simple scaling relation was derived that related the frequency of turbulent sweeps and ejections to the event magnitude. This allows 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 recurrence model modeled described provides a significant step towards a more physically-based blowing snow sublimation model. Additionally, calculations of return 20 frequencies and event durations provide a field-measurement context for recent findings of non-stationarity impacts on sublimation rates.

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]. 35 Represented as an energy balance residual, estimates of these turbulent fluxes rely on an accurate sublimation model, and a precise understanding of the energy 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 40 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 Li, 2000, Dery andYau, 1999;Groot Zwaaftink et al. 2013;Musselman et al, 2015].
Investigations of snow sublimation from a numerical modeling approach have recently provided new insights into non-steady 45 state aspects of sublimation Li et al., 2017;Sharma et al., 2018] and the efficacy of the nearlyuniversally used Thorpe and Mason [1966] model [see Schmidt, 1972] at high temporal and spatial resolution . 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. Extending non-stationary sublimation models to alpine and other complex terrain environments could lead 50 to reduced uncertainty in blowing snow sublimation models.
The objective of this research is to investigate turbulent structures down to sub-second timescales and identify their synchronization with near-surface temperature fluctuations. The study investigates the unsteady processes affecting blowing snow particle energy balances in order 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 55 is tested to represent turbulent event frequency as a function of Variable Interval Time Averaging (VITA) thresholds. Data used to validate this model consist of field measurements of three-dimensional wind velocities and sonic temperatures during blowing snow events at an alpine site in the Fortress Mountain Snow Laboratory (FMSL), Canadian Rockies ( Figure 1). These data are supplemented by observations of nearby temperature, relative humidity, and wind speeds at FMSL so as to analyze potential thermodynamic feedback mechanisms at a larger spatial scale. The scaling relationship of return frequency and event 60 magnitude also gives a real-world context for recent studies on the impacts of non-stationarity on blowing snow sublimation rates.

Fieldwork
Ultrasonic temperature and wind velocity time series were observed in FMSL using two Campbell Scientific CSAT3 sonic 65 anemometers sampling at 50 Hz from November 2015 to March 2016. The anemometers were on the same mast at heights above the snow surface varying over 0.1-0.4 m and 1.4-1.7 m with snow surface accumulation or erosion. Extensive analysis of this dataset has already provided new insights into the turbulent mechanisms for blowing snow transport [Aksamit and Pomeroy, 2016, 2018. The turbulent structures scrutinized here have previously been coupled with Particle Tracking Velocimetry (PTV) and high-speed video analysis of Pomeroy [2017, 2018] to better understand the wind-snow 70 coupling. For each night (20 Nov, 4 Dec, 2015 and3 Feb, 3 Mar 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) was divided into 15-minute 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, January 21, 2016 had much stronger winds because of the presence of a chinook (föhn) event, previously associated with high blowing snow 75 sublimation rates in the area [MacDonald et al., 2018]. The mean temperatures varied from -7 C during the previous three days to +3 C during the night of Jan 21. 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 .
Three other FMSL stations near to the blowing snow measurement site provide complementary 15-minute relative humidity, 80 air temperature and wind speed measurements ( Figure 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 protected forest (Powerline) station approximately 400 m away and 30 m higher in elevation [Smith et al., 2017]. Additionally, there are two exposed sites, include a ridgetop (Canadian Ridge) and lee side of ridge (Canadian Ridge North) that are both approximately 600 m downwind and 85 200 m higher in elevation. 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. " 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. Meteorological variables at the blowing snow study site can be found in Table 1.

Modified VITA Analysis 100
A modified VITA criterion [Morrison et al., 1989] that requires both VITA and quadrant analysis thresholds to be exceeded was applied to instantaneous Reynolds Stress signals to determine periods and types of turbulent motions. This approach harnesses the theory that deterministic coherent turbulent structures can drive blowing snow and are superimposed on incoherent inactive local turbulence [e.g. Sterk at al., 1998;Shih et al., 2017]. The mathematical basis of VITA analysis is the identification of periods of high variance in a turbulent time series, f(t). That is, we first identify periods when 105 where is a statistically or experimentally determined averaging time, is a user-defined threshold, overbar indicates a spatial or temporal average. To increase objectivity, and connect events to physically relevant turbulent structures, the modified VITA analysis used here also includes a quadrant hole analysis [Lu and Willmarth, 1973] criterion [Morrison et al., 1989] 110 which identifies the neighborhood around a VITA event where the Reynolds stress ( ) also exceeds a given threshold: Here, is the density of air, ' and ' are the fluctuating values of streamwise and vertical velocity, and is a user-defined threshold. 115 Conducting the modified VITA analysis over a range of thresholds (k U , k V ) and averaging times (T) provides a relatively objective gust identification scheme that classified significant events into sweeps (u' > 0, w' < 0) and ejections (u' < 0, w' > 0). The modified VITA algorithm categorized a turbulent event as a sweep or ejection if the parameterized curve ( ) = 〈 ] ( ), ] ( )〉 passes through only one of the two quadrants during the event. The same methods were used by Aksamit and Pomeroy [2017] to identify and couple turbulent gusts with blowing snow events. In this study, the concurrent sonic 120 temperature signal response was also measured and the fluctuation from the 15-minute 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 [Campbell Scientific, 2018]. Following Kailas and Narasimha [1994], events detected with larger thresholds are referred to as "stronger" or "more intense." 125 3 Results

Modified VITA Results
During each blowing snow storm, there was no definitive evidence of humidity saturation or thermodynamic feedback at any of the three nearby weather stations ( Figure 2). 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 130 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. For each VITA-identified event, instantaneous temperature deviations from the 15-minute mean were computed to represent the magnitude and sign of turbulent temperature mixing with respect to slower meteorological changes over the nights. Aksamit 140 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 behaviour. The recurrence frequencies and average durations of sweep and ejection events for each threshold https://doi.org/10.5194/tc-2020-46 Preprint. Discussion started: 3 March 2020 c Author(s) 2020. CC BY 4.0 License. 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]. 145 For each blowing snow storm, 3D point-clouds of mean recurrence frequency, event duration and event temperature deviation were calculated for the lower sonic anemometer. Each point represents the values from one choice of averaging time and modified VITA thresholds for a 15-minute observation period as discussed in Section 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 Figure 3. Inset in each subplot are three probability distributions computed from the original point-clouds for each 150 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 (Nov 20, Dec 4, Feb 3, Mar 3) that sweeps consistently brought warmer air to the near-surface anemometer. This can be seen as the coloured temperature plots show average temperature deviations greater than zero for nearly all event duration and frequencies over each storm. Probability distributions 160 show very few sweep events with negative temperature deviations, as well as a consistent positive mean and skewness. The chinook storm on January 21 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. Modified VITA analysis on the upper-anemometer wind and temperature time series showed similar results. Interestingly, the analysis comparing low-anemometer event temperatures to the upper-anemometer means revealed near-surface sweeps often occur with colder signatures than the nearby upper-anemometer means (Figure 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 snowcover with a 170 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 Figure 3i show a group of sweeps that were 1°C warmer than the mean temperature of the lower anemometer, but in Figure 4i, the same group of sweeps were 0.5 C colder than the upper anemometer mean. Color scales are equivalent in Figures 3 and 4. This effect is further supported by the mean anemometer temperatures detailed in Table 1. 175

180
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 lower-anemometer mean (Figure 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 (November 20 and March 3) there was more variability in the temperature contributions from ejections. This may indicate stable layers of 185 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 Section 4.
Over all nights, sweeps were of longer duration than ejections and had a higher frequency of occurrence. It is interesting to 190 note that the probability curves in Figure 3 show a second sweep return frequency peak around 0.5 Hz for all nights. This is not present in the ejection frequency probabilities, which only has 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 outerlayer and inner-layer motions as previously suggested by Aksamit and Pomeroy [2017]. The ejections occur less often because 195 of the rarity of large positive ' values close to the snow surface, and are thus present only during a less common generating mechanism.

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: 200 = h jk(l m jn) . (

3)
Here, is the recurrence frequency of a given modified VITA turbulence event type, h and are fitting parameters with h known as the characteristic frequency. This scaling analysis focused on the case where 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 205 (small standard deviation) in the log of the return frequency, ( ), was present over varying averaging times for each VITA threshold v . This resulted in a good fit of Eq. 1 for the return frequencies of the total number of modified VITA events, as well as for sweeps and ejections individually. An example of this fitting for one 15-minute period on February 3 is shown in Figure 5. The squared ℓ " -norm of the residuals for each minimized least-squares fit are presented in the document supplement, as are the characteristic frequencies, h . 210 h varied little between blowing snow storms with mean values of 0. 80, 0.54, 0.49, 0.35, 0.1 and 0.09 for h,xy , h,z{ , h,xy |} , h,z{ |} , h,xy~• ,and h,z{~• respectively. This suggests persistent flow features at this site that may be due to a common characteristic topographically influenced flow. As could be expected from the analysis presented in Figure 3, h values for total events ( h,xy , h,z{ ) and for sweeps ( h,xy |} , h,z{ |} ) were greater than those for ejections. Of particular interest in this scaling relationship is a clear difference between h for the upper and lower anemometer observations for both total events 215 and solely sweep events. Over all nights, the characteristic frequency for total events was lower at the upper anemometer, which corresponded with a drop in the number of sweeps, whereas the characteristic frequency of ejections was nearly identical at both heights. An example of this can be found for one fifteen-minute period on December 4 in Figure 5, with further supporting evidence in the scaling relationship data table in the document supplement.
The threshold criteria in Eq. (1) and (2) varies for measurement location and time, scaling by mean values calculated over each 220 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. As well, this supports the suggestion in Section 3.1 that the mechanisms generating sweeps and ejections may be different, with less common flow features resulting in the ejections.

Discussion 225
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 sublimation of surface and blowing snow may occur if moisture 230 fluxes near the surface are only counterbalanced by diffusion [Bintanja, 2001]. Dover and Mobbs [1993]; Dery andTaylor [1996], Groot Zwaaftink et al. [2013] and others have suggested that blowing snow sublimation could be a self-limiting process when thermodynamic feedbacks were included in a steady-state boundary layer. However, these models did not account for warm-or dry-air entrainment, nor 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 studies in the steppes of Russia, high plains of Wyoming 235 (USA), prairies of Saskatchewan, alpine mountains of Alberta and arctic tundra of the Northwest Territories (Canada) [e.g. Dyunin, 1959;Schmidt, 1982;Pomeroy and Li, 2000;Musselman et al., 2015]. The evidence of frequent regeneration of warm air near the surface through advection or entrainment processes helps explain the discrepancy with diffusion-dependent models.
Furthermore, the regeneration process may explain why sublimation studies can result in a wide range of sublimation or 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. 255 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]. The potential improvement to 260 sublimation calculations by inclusion of a physically-based energy advection like this scaling relation is a topic of further investigation.
Future high temporal resolution studies of air temperature and water vapour during sustained periods with 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 265 multiple heights during blowing snow events. This would allow a more complete understanding of the advectionthermodynamic feedback balance during blowing snow storms and advance the seminal profile studies of Schmidt [1982].

Conclusion
Field measurements of atmospheric conditions during blowing snow in an alpine environment show no evidence of saturation 270 of water vapour during blowing snow events, during which frequent sweep and ejection motions enhanced atmospheric mixing of warm and cold air in the near-surface region. Sonic temperature signals indicate that sweeps bring relatively warm air to the 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.5 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 275 entrainment during snow transport may explain why the thermodynamic feedback parameterizations necessary in many blowing snow sublimation models are unphysical. An enhanced influence of mechanical mixing in boundary layers with inhomogeneous temperature distributions, for example where there is cold-air pooling, 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 280 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, 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 285 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.
The present research has suggested a simple similarity scaling of the return frequency of turbulent events of intensity v as identified by modified VITA analysis, through the exponential relationship of Kailas and Narasimha [1994]. This framework 290 provides a simple platform with which to model or investigate a gust-driven regeneration function of warm-dry air in the nearsurface for blowing snow sublimation calculations. With the inclusion of such a statistical recurrence model, 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. 295 Author Contributions NA and JW designed the experiment and contributed to the evaluation of the results. NA performed the field experiment and analysis. Both authors contributed to the writing of the manuscript.