Elucidation of spatiotemporal structures from high-resolution blowing-snow observations

Systematic observations were conducted to investigate the spatiotemporal structures of blowing snow. Along a line perpendicular to the dominant wind direction on the lee side of a flat field, 15 snow particle counters (SPCs) and ultrasonic anemometers (USAs) were placed 1.5 m apart. Data were recorded at high frequencies of 100 kHz for SPCs and 1 kHz for USAs. The horizontal mass flux distributions, representing the spatiotemporal variability of blowing snow, exhibited non-uniformity in both time and space and manifested periodic changes akin to snow waves. Additionally, the presence of “snow snakes”, meandering near the snow surface, was observed. Quadrant analysis revealed predominant snow fluxes in quadrants Q1 ($u^{\prime }>\mathrm{0}$, $w^{\prime }>\mathrm{0}$) and Q4 ($u^{\prime }>\mathrm{0}$, $w^{\prime }<\mathrm{0}$). However, a more detailed parametric curve analysis indicated the existence of ejection events in Q2 ($u^{\prime }<\mathrm{0}$, $w^{\prime }>\mathrm{0}$) before snow waves and in front of snow snakes, shifting to Q1 and Q4 afterward, implying the consideration of both top-down and bottom-up mechanisms for burst sweep events.

1 Introduction

The transport of snow by wind, known as blowing or drifting snow, plays a significant role in engineering and climatological contexts. Near roads, blowing snow can cause snow drifts and reduce visibility, leading to traffic accidents. Wind-driven snow redistribution can result in localized snow accumulation and the formation of snow cornices, influencing snow avalanche danger. The impact of katabatic winds on the Antarctic ice sheet affects the mass and energy balance, becoming more critical when considering global climate change effects.

Numerous investigations have been carried out, both in the field and in wind tunnels, to explore the meteorological conditions that trigger blowing snow and to gain insights into its internal structures. However, the model of snow transport is divided into two main categories: Lagrangian models, which focus on the motion of individual snow particles in the saltation layer near the snow surface, and turbulence diffusion models, which deal with the suspension layer. Notably, Nemoto and Nishimura (2001) developed a random-flight model of blowing snow that incorporates turbulence, the inertia of snow particles, and particle size distribution. This model successfully captures the entire transition from the saltation layer to the suspension layer during a blizzard. However, it is limited to steady-state conditions and does not address temporal and spatial variations.

To gain a more comprehensive understanding, it is crucial to consider organizational structures such as turbulence sweeps and ejections, as well as intermittent dynamics and non-linear dependence on wind speeds, when discussing the onset and development of snow transport. Unfortunately, the existing time-averaged models based on local equilibrium theory and isotropic turbulence theory have limitations in determining coherent structures with regularity due to their restricted flow field representation.

Moreover, it has been observed that traffic accidents caused by reduced visibility on roads are strongly influenced by temporal and spatial variations rather than the intensity of blowing snow itself (Kajiya et al., 2001). Consequently, the demand for models capable of describing these temporal and spatial structures is widely recognized. Recently, efforts have been made to develop models that combine large-eddy simulation (LES) with Lagrangian snow transport models (e.g., Groot Zwaaftink et al., 2014; Okaze et al., 2018). However, these models lack comprehensive measurements that adequately address the high variability of these phenomena.

While a study by Aksamit and Pomeroy (2016, 2017, 2018) focused on temporal changes, it unfortunately lacked spatial information. As a result, the reliability of these models remains unconfirmed, and our understanding of the high-frequency fluctuations in the internal structures of blowing and drifting snow remains unsatisfactory. Against this backdrop, our study aims to systematically measure blowing and drifting snow to investigate their spatiotemporal structures. To achieve this goal, we have deployed 15 snow particle counters (SPCs) in designated test areas and conducted measurements using an equal number of ultrasonic anemometers, providing high-temporal-resolution data. While the dataset obtained from our series of observations is indeed substantial, this study primarily offers a brief overview of these observations and the results obtained.

2 Field observation

The field observations were conducted during the winter of 2018 in Tōbetsu, Hokkaido, Japan. Tōbetsu is located close to the coast of Ishikari Bay and experiences strong monsoon winds from the west during the winter. The measurement apparatus zone was situated in the windward area, where a flat snow field spreads over a distance of more than 200 m (Fig. 1).

Figure 1 Observation field: Tōbetsu, Hokkaido, Japan (a). Measurement apparatus: 15 snow particle counters (SPCs) and ultrasonic anemometers (USAs) set on towers 1.5 m apart (b).

The 15 measurement towers were installed along a line perpendicular to the dominant wind direction, with each tower spaced 1.5 m apart. The towers were equipped with snow particle counters (SPCs; SPC-7 and SPC-95, Niigata Denki Co. Ltd) and ultrasonic anemometers (USAs; Young 81000) to measure both the spatiotemporal structures of blowing snow and the airflow (Fig. 1). The height of the SPCs was initially 39 cm, and the height of the USAs was 107 cm, but both heights decreased with increasing snow accumulation.

The SPC used in this study (Niigata Denki Co. Ltd) is an optical device (Nishimura and Nemoto, 2005) that measures the diameter and the number of drifting snow particles by detecting their shadows on a photodiode (assuming that drifting particles are spherical in shape). Electric pulse signals resulting from snow particles passing through the sampling volume (2 mm × 25 mm × 0.5 mm) are sent to a transducer and an analyzing data logging system (PC). In this way, the SPC is able to detect particles in the range of 40–500 µm. The analysis software divides the particles into 64 size classes and records the number of particles in each size class at 1 s intervals. The SPC is mounted on a self-steering wind vane, and hence the sampling region, which has a cross-sectional area of 2 mm × 25 mm (50 mm²), is maintained perpendicular to the horizontal wind vector. If the diameter of a snow particle is larger than that of the maximum-diameter class, the snow particle is considered to belong to the maximum-diameter class. Usually, SPCs are used to obtain snow particle size distribution and mass flux at 1 s intervals. However, in this study, the output signals from the SPC transducers were directly recorded with a high frequency of 100 kHz. This allowed the researchers to calculate not only the mass flux but also the particle size and speed with fine time resolutions (Nishimura et al., 2014). Interestingly, the SPC has also found application in sand transport research as a “sand particle counter” (Yamada et al., 2002; Mikami et al., 2005). All SPC sensors are accurately calibrated in advance using spinning wires of various diameters, enabling us to effectively account for any sensitivity differences. Detailed procedures are provided in Sato et al. (1993). Further, over the snow surface, we can generally avoid the effect of the dust and fine-soil contamination over the optical parts, unlike the positions over sand and soil surfaces. Numerous long-term observations have been conducted in the Alps (e.g., Naaim-Bouvet et al., 2010, 2014; Guyomarc'h et al., 2019), the Arctic (e.g., Lenaerts et al., 2014; Frey et al., 2020), and Antarctica (e.g., Sigmund et al., 2021; Wever et al., 2023), all of which attest to the reliability of the SPC. Generally, we can leave the system even for the whole winter without cleaning or wiping the optical parts. In particular this campaign was carried out only for less than 2 months. Thus, we believe the contamination on top of the optics can be reasonably neglected.

The wind speeds measured with USAs were stored with a frequency of 1 kHz. To ensure accurate synchronization of time for both recordings, precise calibration and alignment were performed. As elucidated further in subsequent sections, the heights of the sensors from the surface fluctuated due to snow accumulation. Consequently, we took the necessary step of manually measuring the sensor heights of both the snow particle counters (SPCs) and ultrasonic anemometers (USAs) at every observation period. This meticulous approach ensured that any variations in sensor height were accurately accounted for, maintaining the integrity and reliability of our measurements throughout the experiment.

Other meteorological data, such as mean wind speeds and directions, and air temperature were obtained from an automatic weather station (AWS) located approximately 1 km to the northwest (windward direction). Additionally, data on precise precipitation were available from the Double Fence Intercomparison Reference (DFIR) set at the Ishikari Blowing-Snow Test Field of the Civil Engineering Research Institute for Cold Region (CERI), located at a distance of nearly 4 km to the west. The DFIR data were primarily used to determine whether the observed phenomena were purely drifting or blowing snow or if they were accompanied by snowfall.

Overall, the field observations were conducted in a well-equipped measurement zone, allowing the researchers to collect detailed data on blowing snow and atmospheric conditions, which are crucial for their study of spatiotemporal structures of blowing snow.

Figure 2 Meteorological conditions at the observation site in February 2018. Red arrows indicate the duration of the intense observation period.

3 Results and analysis

In February 2018, the 10 min average of the meteorological conditions was recorded using an AWS and DFIR, as depicted in Fig. 2. The red arrows indicate the duration of the intense observation period. Despite not having an exceptionally high wind speed at 2 m above the surface, the westerly wind was sufficient for initiating drifting/blowing snow. Throughout the period, the air temperature remained mostly below 0 °C, causing the snow depth to increase from 80 to 110 cm in just 1 month. Consequently, the sensor height of the snow particle counter (SPC) decreased from 39 to 1 cm due to the snow accumulation.

During the intense observation period, the Richardson number (Ri), calculated at the AWS site, using vertical temperature and wind speed profiles (Sutton, 1953), indicated neutral conditions ($\mathit{Ri}<\left|\mathrm{0.01}\right|$).

Figure 3 Time series of wind speed, snow mass flux, particle speed, and particle diameter at the spanwise direction Y=9 m. (a) On 16 February, the sensor heights for the USA and SPC were 91 and 23 cm; (b) on 24 February, they were 76 and 8 cm.

Figure 3 illustrates the recordings of wind speed from the ultrasonic anemometer (USA) and snow mass flux, particle speed, and particle diameter obtained with the SPC at the spanwise direction of Y=9 m for 60 s. On 16 February, the sensor heights for the USA and SPC were 91 and 23 cm, respectively, while on 24 February, they were 76 and 8 cm. Although the sensor height of the USA is around 70 cm higher than that of the SPC, the change in snow mass flux aligned well with the wind speed on both days. In particular, the response on 24 February was more pronounced, suggesting that the snow flux in the saltation layer is more sensitive to wind speed fluctuations than the one in the higher suspension layer. Figure 3 reveals that some particles move rapidly, even surpassing the wind speeds. These particles likely originate from higher positions, where they are accelerated by stronger winds. However, in general, particle speeds are lower than the wind speed, as reported in Nishimura et al. (2014), who found that mean particle speeds are consistently 1 to 2 m s⁻¹ less than wind speeds at heights of 0.015 to 1 m. Although the sensor height of the SPC was lower than that of the USA, rough estimates based on the logarithmic wind profile support this trend. While particles smaller than 100 µm in diameter contributed more on 16 February, they displayed a wide distribution within the range of 40–500 µm on 24 February. Decreasing the height from the snow surface revealed larger particle numbers with a wide size distribution.

Figure 4 Horizontal mass flux distributions for 120 s, reflecting the spatiotemporal variability of blowing snow, on 12, 16, and 24 February. The measured heights from the snow surface were 31, 23, and 8 cm, respectively. Advected distance is also shown on the top x axis. F: snow mass flux.

Figure 5 Autocorrelations of horizontal mass flux shown in Fig. 4. Each line indicates the output from the individual SPC.

Figure 4 displays horizontal mass flux distributions, reflecting the spatiotemporal variability of blowing snow, measured with 15 SPCs for 120 s on 12, 16, and 24 February. The measured heights from the snow surface were 31, 23, and 8 cm, respectively. At the AWS point, wind direction was kept steadily as west and wind speeds at 2 m high on the 3 separate days were nearly the same: from 4.5 to 5 m s⁻¹. Further, precipitation observed with DFIR was quite small and can be neglected. The three panels in Fig. 4 depict maps of snow transport intensity, created by aligning time series data from the snow particle counter (SPC) transverse array in the spanwise y direction and applying Taylor's (1938) frozen hypothesis to substitute time for the streamwise direction. In the figure, the “advected distance” is presented on the x axis, estimated using the average wind speed provided by the ultrasonic anemometer (USA) and correcting for the height differential between the USA and the snow particle counter (SPC) sensor heights, assuming the logarithmic wind profile.

Despite steady westerly wind direction and similar wind speeds at 2 m height on the 3 separate days, the flux distributions were uniform with neither time nor space and exhibited periodic time – in other words, lateral to the streamwise direction – variation. The autocorrelations of horizontal mass flux in Fig. 4 are depicted in Fig. 5. Each line in the figures represents the output from an individual SPC. While the distributions on 24 February appear somewhat subdued, distinctive peaks are evident on 12 and 16 February, occurring approximately every 10, 35, 70, and 90 s. These correspond to distances of approximately 30, 115, 230, and 300 m, respectively.

Figure 6 “Snow waves” appeared on the field's upstream at night on 20 February (a), as did “snow snakes” (b). In both cases the wind blew from back to front. As shown in the figures, the snow waves are organized in a lateral or spanwise orientation and the snow snakes are quasi-parallel to the streamwise direction.

Figure 6, which shows pictures of the field's upstream at night on 20 February, highlights the periodical snow waves approaching the sensors from back to front. These waves correspond to the periodical change in snow flux observed in Fig. 4.

Although the formation mechanism of the snow wave has not been clarified yet, Kobayashi (1980) suggested that the snow waves strongly reflect the turbulent-flow structure of the wind, and its formation mechanism is similar to the wavelike cloud in the atmosphere and the wind wave in the ocean.

Figure 7 Horizontal mass flux distributions for 10 s at heights of 23, 8, and 1 cm from the snow surface. Advected distance is also shown on the top x axis.

Figure 8 The 2D autocorrelation of horizontal mass flux on 24 February is shown on the bottom of Fig. 7.

For 16 and 24 February, Fig. 7 displays the blowing-snow flux distributions for 10 s at heights of 23, 8, and 1 cm from the snow surface. Due to the substantial increase in snow depth, nearly 10 cm within a single day on 24 February, as depicted in Fig. 2, we have introduced two scenarios with SPC heights set at 8 and 1 cm. Contrary to Fig. 4, as the height decreased from the surface, long and narrow structures prominently emerged. These structures seem to correspond to the phenomena known as “snow snakes”, which meander laterally, merging and bifurcating as they move downwind as shown in Fig. 6, similar to observed phenomena in drifting sand movement. The height of the long structures along the snow surface varied, and taller parts occasionally protruded to higher positions, as seen in the cases of 23 and 8 cm in Fig. 6. Figure 8 illustrates the 2D autocorrelation of horizontal mass flux on 24 February depicted at the bottom of Fig. 7. The 1.5 m spacing between SPCs may not always be sufficient for capturing precise structures. However, the Delaunay triangulation (Cheng et al., 2012) used in Fig. 8 is a powerful tool for interpolation, mesh generation, and graphical applications. It is widely used in geographic information systems (GISs) to create terrain models, where triangulating elevation points construct a surface that accurately represents the terrain with minimal distortion. Consequently, the structures observed at 1 cm above the surface, with widths around 2 m, peak fluxes about 30 cm wide, and lateral spacing of approximately 5 m, as shown in Fig. 7, are quite plausible. Notably, the 2D autocorrelation of horizontal mass flux in Fig. 8 indicates a lateral spacing of about 5 m, which is more than 3 times the 1.5 m sensor spacing.

The propagation speed was roughly estimated to be 2 to 3 m s⁻¹ based on particle speed data from analyzing the SPC data (Fig. 3b).

Figure 9 Power spectra of mass flux and wind speeds at Y=10.5 m on 24 February.

Figure 9 displays the power spectra of mass flux at Y=10.5 m for 24 February in Fig. 4, along with the ones for wind speeds. The dominant frequencies and general trends of the power spectra for snow flux and wind speed appear very similar, implying that both vary in a correlated manner. Despite the sensor heights of SPCs being 68 cm lower than those of anemometers, these observations suggest that changes in snow flux reflect the structures of turbulence eddies near the snow surface.

Figure 10 Horizontal distributions of snow mass flux, wind speed, and ejection (Q4: red) and sweep (Q2: blue) events on 12 and 24 February 2018.

We applied quadrant analysis in this study and focused on turbulent sweeps and ejections to better understand the spatiotemporal structures of blowing snow. While previous discussions have centered on the amount of snow transport by wind, often using averaged wind speed over specific time periods, such as the 3rd to the 5th power of mean wind speed (u³ to u⁵) for 10 min (e.g., Dyunin, 1967; Mann et al., 2000; Nishimura and Hunt, 2000), in order to elucidate the spatiotemporal structures of blowing snow, it is necessary to set on our focus on the turbulent structures of the wind. To explore the correlation between snow transport and turbulent coherent structures in the boundary layer, quadrant analysis (Wallance, 2016) was applied. This classical analysis examines the fluctuating wind velocities of the streamwise and vertical components (u^′ and w^′) and their product, the Reynolds stress, to identify ejection in Q2 ($u^{\prime }<\mathrm{0}$, $w^{\prime }>\mathrm{0}$), where lower-momentum air is expelled away from the bed, and sweep in Q4 ($u^{\prime }>\mathrm{0}$, $w^{\prime }<\mathrm{0}$), where higher-momentum air is pushed towards the bed, events contributing to turbulent energy production. Figure 10 shows the horizontal distributions of snow mass flux, wind speed, and sweep (Q4) and ejection (Q2) events on 12 and 24 February 2018. It is evident that large snow mass flux, corresponding to the ridge of the snow wave (in Fig. 6), occurs under conditions of high wind speeds and sweep events (Q4; $u^{\prime }>\mathrm{0}$, $w^{\prime }<\mathrm{0}$) on both days.

Figure 11 Correlation between the horizontal snow flux at 1 cm and the absolute value of the kinetic shear stress $u^{\prime }{w}^{\prime }$ in each quadrant on 24 February 2018.

Figure 12 Relations between the snow flux, the fluctuating component of u^′ and w^′, and the product of u^′ and w^′ on 24 February 2018. All data and the observation period are the same as Fig. 11.

Table 1 Contributions of each quadrant on 12, 16, and 24 February at corresponding heights of mass flux from the snow surface: 31, 16, 8, and 1 cm.

Figure 11 illustrates the correlation between the horizontal snow flux at 8 cm and the absolute value of the kinetic shear stress $u^{\prime }{w}^{\prime }$ in each quadrant on 24 February 2018. The results indicate that snow fluxes were predominantly observed in quadrants Q1 ($u^{\prime }>\mathrm{0}$, $w^{\prime }>\mathrm{0}$) and Q4 ($u^{\prime }>\mathrm{0}$, $w^{\prime }<\mathrm{0}$), while they were far less prominent in quadrants Q2 and Q3. Table 1 displays the contributions of each quadrant, acquired on 12, 16, and 24 February, at corresponding heights from the snow surface: 31, 16, 8, and 1 cm; analysis was conducted with all the sensor data. Ejection events (Q2; $u^{\prime }<\mathrm{0}$, $w^{\prime }>\mathrm{0}$) and sweep events (Q4; $u^{\prime }>\mathrm{0}$, $w^{\prime }<\mathrm{0}$) contributed 12 % to 22 % and 33 % to 48 %, respectively. Q4 exhibited the largest contribution, while Q1 consistently followed in all events, regardless of the height from the surface. Moreover, the combined contributions of Q1 and Q4 account for over 64 %, while Q2 and Q3 contribute less than 36 %. The contributions from quadrants Q2 and Q3 ($u^{\prime }<\mathrm{0}$, $w^{\prime }<\mathrm{0}$) were notably smaller compared to the significant impact of Q1 and Q4. Figure 11 suggests a decrease in snow flux with increasing Reynolds stress, which may initially appear contradictory to conventional understanding. To further explore this relationship, Fig. 12 examines the correlation between snow flux and the fluctuating components of u^′ and w^′. Panel (a) demonstrates a clear increase in snow flux with u^′, while panel (b) illustrates a distribution concentrated around 0, resembling a normal distribution. Panel (c) reveals a negative trend between u^′ and w^′, indicating that w^′ decreases as u^′ increases. The relationship between the product of u^′ and w^′ and snow flux is depicted in panel (d). Although setting specific hole sizes in Fig. 11 may result in variations, Leenders et al. (2005) suggested that vertical fluctuations are inherently constrained by the distance from the bottom boundary and the overall scale of the structure, leading to a poor correlation between surface fluctuations and vertical fluid motions higher in the profile. Consequently, measurements of Reynolds shear stress around 1 m from the surface are typically poorly correlated to snow transport flux.

Figure 13 Distributions of the snow mass flux at 1 cm, as well as the wind speed, ejection, and sweep structures at 69 cm for 10 s on 24 February.

Figure 14 (a–e) Wind speeds and mass fluxes at Y=4.5, 10.5, and 15.0 m from 25 to 35 s and at Y=10.5 and 15.0 m from 60 to 70 s from Fig. 8. Furthermore, parametric curves of (u^′(t), w^′(t)) were plotted every 1 s with different colors for the period indicated by arrows in each figure: from 30 to 35 s for the former and from 67 to 70 s for the latter. The blue arrows indicate the onset of the snow flux increase.

Figure 15 Structures around the snow snake: wind speed and snow flux for 10 s at Y=15.0 m and the parametric curve from 1 to 5 s in Fig. 10. The blue arrow indicates the onset of the snow flux increase.

In contrast, Fig. 13 illustrates the distributions of the snow mass flux at 1 cm, as well as the wind speed, ejection, and sweep structures at 69 cm for 10 s on 24 February. The data indicate that the foremost region of the snow snake primarily aligns with the ejection phenomenon, trailed closely by the sweep motion. For a more detailed examination of the structures preceding and during the onset of transport escalation, specific 10 s sub-periods (marked with arrows in Fig. 10 and whole durations in Fig. 13) have been extracted and are presented in Figs. 14 and 15. Figure 14a–e present the wind speeds and mass fluxes at Y=4.5, 10.5, and 15.0 m from 25 to 35 s and at Y=10.5 and 15.0 m from 60 to 70 s in Fig. 10. Furthermore, using the approach introduced by Aksamit and Pomeroy (2017), parametric curves of (u^′(t), w^′(t)) were displayed, representing the duration of the time indicated by arrows in each figure: from 30 to 35 s for the former and from 67 to 70 s for the latter. These curves were plotted every 1 s with different colors.

For example, in Fig. 14a, no snow transport occurred between 30 and 31 s at Y=4.5 m, followed by an abrupt increase in snow flux. Afterward, the mass flux diminished to nearly zero from 31.5 to 32.5 s and then increased again from 33 to 34 s. The parametric curve of (u^′(t), w^′(t)), encompassing the entire first peak of the horizontal mass flux (pre-increase, increase, and decrease) between 30 and 33 s, was situated within quadrant Q2 (ejection region). However, the second increase occurred in quadrants Q1 and Q4, characterized by stronger horizontal wind speeds. Similarly, at Y=10.5 m in Fig. 14b, from 32 to 34 s, the parametric curve mainly remained in Q2, covering the duration before the increase and the flowing period leading to the maximum, and then moved to Q1. In Fig. 14c, which shows the case at Y=15.0 m, the snow flux started to increase at 33 s. The parametric curve was mostly in Q4 1 s before. However, at the onset of the increase at 33 s, it shifted to Q2 and later moved to Q1. Figure 14d illustrates the case at Y=10.5 m from 60 to 65 s, where the snow flux started to increase at 68 s and reached its peak at 69.1 s. The parametric curve from 67 to 70 s indicated that both periods, 1 s before the increase in snow transport and the path to the maximum, were in Q2 and subsequently moved to Q1 and Q4, where the horizontal wind speed was greater. The snow flux change at Y=15.0 m shown in Fig. 11e was almost similar to that at Y=10.5 m. Although the staying period of the parametric curve in Q2 was rather short compared to other cases, it was in Q2 during both the period just prior to the increase in mass flux and the front part and then shifted to Q1 and Q4. Overall, the parametric curve patterns indicated that the ejection part in Q2 exists prior to the snow waves, while Q1 and Q4, where horizontal wind speed is higher, follow.

In Fig. 15, structures around the snow snake are displayed, showing wind speed and snow flux for 10 s at Y=15.0 m and the parametric curve from 1 to 6 s. The front of the snow snake arrived at the sensor at around 2.5 s, and the snow flux increased rapidly. During 1 to 2 s, the parametric curve remained in Q1, but it shifted to Q2 just in front of the snake head and then moved to Q1 and Q4. It is worth noting that the curve briefly returned to Q2 before the second mass flux peak. During the flowing period, the curve predominantly stayed in Q1 and Q4, where u^′ is positive. In general, the parametric curve follows the ejection part in Q2 in front of the snake head and then shifts to Q1 and Q4, similar to the snow wave behavior in Fig. 14a–e.

4 Conclusions

We conducted systematic measurements of blowing snow, which allowed us to reveal the spatiotemporal structures of the phenomenon. In comparison, research on sand transport, particularly eolian streamers, which correspond to the snow waves and the snow snakes in this study, has been explored in studies such as Bauer et al. (1998, 2013) and Sterk et al. (1998). Bauer et al. (1998) used a hot-wire anemometer and load-cell sand traps, identifying three typical transport patterns: streamer families, nested streamers, and clouds with embedded streamers. However, these latter two patterns occurred only at high shear velocities. Comparing their classification with our study's findings in Figs. 4 and 7, the structures observed on 16 and 24 February likely correspond to streamer families and nested streamers. Yet, “clouds with an embedded streamer” were not found in the snow transport, likely because the wind speeds were relatively low and insufficient for transporting large amounts of snow particles. It is worth noting that the periodical changes in snow flux shown in Fig. 4 were not mentioned in Bauer et al. (1998, 2013) for the sand transport, possibly because their focus was on phenomena occurring within a short period of 5 s, while the snow waves recognized in Fig. 4, with frequencies less than 0.1 Hz, were observed over 120 s.

On the other hand, Sterk et al. (1998) utilized a Gill-type anemometer and acoustic sediment sensing instruments to investigate the relationship between sediment flux and Reynolds shear stress. They reported very low levels of statistical explanation, with correlation coefficients smaller than 0.14, which aligns with our findings, as depicted in Fig. 12d. Furthermore, Sterk et al. (1998) partitioned instantaneous sediment transport values and sorted them according to wind speed quadrant signatures. They observed that the largest mean saltation fluxes occurred during stress excursions located in quadrants Q1 and Q4, while excursions in quadrants Q2 and Q3 were unable to sustain appreciable saltation activity. This observation is consistent with the findings of Shonfeldt and von Lowis (2003), Leenders et al. (2005), and Wiggs and Weaver (2012), who employed sonic anemometers for sand transport observations. Our study corroborates these findings, as summarized in Table 1, where contributions from quadrants Q1 and Q4 were significantly larger than those from Q2 and Q3. This supports the notion that horizontal velocity fluctuations (u^′) play a pivotal role in eolian sediment transport compared to vertical fluctuations (w^′) or kinematic shear stress ($u^{\prime }{w}^{\prime }$).

Baas and Sherman (2005) applied the variable internal time-averaging (VITA) method to wind speed data, which identifies significant turbulent events using local variance within a moving time window (Blackwelder and Kaplan, 1976). They found that VITA events were generally poorly correlated with transport events and sweep-like motions with strongly positive u^′ signals dominated over ejection-like motions. Then, they interpreted the streamer patterns as having a probable origin in the large “eddy surface layer”, where turbulent eddies from higher regions of the atmospheric boundary layer travel downward, elongate and stretch out in the shear layer, and scrape across the surface. Later, Hunt and Morrison (2000) proposed a top-down model of boundary layer turbulence for very high Reynolds number flows. Aksamit and Pomeroy (2017) also applied quadrant analysis and VITA techniques to their snow transport measurements, revealing that sweeps were the dominant motion for initiating blowing snow and increasing the concentration and particle number flux near the surface. Active ejections occurred during periods of active blowing snow and after sweeps.

Our measurements also represent that the snow fluxes were observed largely in quadrants Q1 ($u^{\prime }>\mathrm{0}$, $w^{\prime }>\mathrm{0}$) and Q4 ($u^{\prime }>\mathrm{0}$, $w^{\prime }<\mathrm{0}$). However, we should note that the detailed analysis with the parametric curve revealed that the ejection part in Q2 exists both prior to the snow waves (in Fig. 14a to e) and in front of the snake head (in Fig. 15) and then shifts to Q1 and Q4, where horizontal wind speed becomes larger. Thus, in addition to the top-down mechanism proposed by Baas and Sherman (2005), bottom-up mechanisms of burst sweep events cannot be ruled out at this stage. Indeed, conducting detailed measurements with SPCs and USAs positioned at the same height from the surface would provide more accurate and directly comparable data. Additionally, incorporating the latest turbulence analysis techniques could further enhance our understanding and facilitate reaching conclusive findings. By leveraging advanced methodologies and instrumentation, such as high-resolution SPC arrays and state-of-the-art turbulence analysis algorithms, we can obtain more precise measurements and delve deeper into the complex interactions between blowing snow and turbulent-flow dynamics. These efforts would contribute significantly to advancing our knowledge in this field and ultimately lead to more comprehensive conclusions.

As stated in the Introduction, time-averaged blowing-snow models have limitations, necessitating the incorporation of relationships with organizational structures such as turbulence sweeps and ejections to discuss the onset and development of snowstorms more accurately. Models that can describe temporal and spatial structures are crucial. To address this, we have developed a large-eddy simulation coupled with Lagrangian snow transport model LLAST (Okaze et al., 2018). This model combines large-eddy simulation for turbulent flow and particle translational simulations based on Newton's equations of motion. The dataset obtained from the series of observations in this study is substantial, providing ample opportunities to derive concrete conclusions elucidating the spatiotemporal structures of blowing snow. However, it is important to note that our analysis is still quite limited at this stage, and only preliminary results are introduced here. Nevertheless, more detailed analyses, which will be presented in subsequent works, hold the promise of significantly contributing to the improvement and validation of the model.