Hysteresis and Surface Shear Stresses During Snow-Particle Aeolian Transportation

Abstract

Surface shear stresses produced by wind and particle collision play a key role in aerodynamic entrainment and splash processes. The fluid shear stress at the surface during aeolian transport has been researched for decades; however, the equilibrium property reported in the literature, numerical simulations, and experiments is inconsistent. To discuss this discrepancy, this study investigates fluid and particle shear stresses at the surface during the aeolian transport of snow particles using a two-dimensional random-flight model of drifting snow. The simulations are performed for various friction velocities on a loose snow bed. By varying the wind conditions in stages, the transport hysteresis is confirmed, and the impact threshold is estimated from the particle transport rate ($0.206\,\hbox {m\,s}^{-1}$). The friction velocity at the surface during transport decreases marginally with an increase in wind speed caused by the impact threshold, revealing that our results do not contradict Owen’s second hypothesis. The total shear stress, which is calculated by summing the fluid and particle shear stresses, is vertically uniform in the equilibrium state; thus, the increase in the particle shear stress decreases the fluid shear stress at the surface. The equilibrium property of the fluid shear stress near the surface changes significantly with height (from a decreasing trend to an increasing trend) because the particle shear stress decreases rapidly in the height range of 1–10 mm. Our findings suggest that it is difficult to accurately measure the fluid shear stress in the surface vicinity using anemometers, and a new methodology is needed.

1 Introduction

Surface shear stresses are related to the start and stop of snow transport, and the accurate quantification of these stresses is essential in predicting the occurrence distribution of blowing snow and the snow redistribution that occurs due to the flow in plain and mountain regions. The aeolian transport of snow or sand particles is generally developed and maintained through four physical sub-processes: aerodynamic entrainment, wind-blown particle dynamics, splash processes (i.e., particle–surface collision), and wind modification (Bagnold 1941; Kok et al. 2012). The transport of snow is more complex than that of sand; the shape and size of snow-bed particles tend to change with the snow transformation. Moreover, the bond between particles is a result of snow sintering, and snow sublimation occurs during the transport. Strictly speaking, there are differences between the transport processes of snow and sand particles; however, we have described them interchangeably hereafter. For the start and stop of aeolian transport, two different types of critical thresholds exist: fluid and impact thresholds of the friction velocity (i.e., $u_*^\mathrm{f}$ and $u_*^\mathrm{i}$, respectively). The friction velocity $u_*$ is defined using the fluid shear stress $\tau$ and the fluid density $\rho _\mathrm{f}$ as $u_* = \sqrt{\tau /\rho _\mathrm{f}}$.

When the friction velocity acting on the surface exceeds the fluid threshold, the grains on the surface migrate downwind, and aeolian transport is initiated due to aerodynamic entrainment. The fluid threshold of sand and snow particles is proportional to the square root of their diameter (Bagnold 1941; Clifton et al. 2006), $u_*^\mathrm{f} \propto \sqrt{d}$, which can be derived by balancing the aerodynamic drag and gravity force acting on the deposited particle. The relationship between $u_*^\mathrm{f}$ and d holds for coarse particles with diameters greater than approximately 100 $\upmu$m, whereas the fluid threshold of finer particles incorporates lift and interparticle forces (Iversen and White 1982; Shao and Lu 2000; Kok et al. 2012). For snow particles, the fluid threshold is affected by the interparticle force that depends on the particle diameter and the bond between particles (Clifton et al. 2006). Once aeolian transport is initiated, it can be maintained at a friction velocity less than the fluid threshold. Airborne grains are deposited on the surface by decreasing the friction velocity (or wind speed), and aeolian transport ceases at the impact threshold. Recent numerical models incorporating turbulent fluctuations have reproduced the intermittent saltation around thresholds (Okaze et al. 2018; Pähtz et al. 2020; Rana et al. 2020, 2021; Zheng et al. 2020); thus, aeolian transported particles at the non-steady state exhibit more complicated transport properties. Bagnold (1941) estimated that the relationship between the fluid and impact thresholds is approximately $u_*^\mathrm{i}/u_*^\mathrm{f}\approx 0.85$, which implies that two stable states (transport and no transport of particles) exist at friction velocities between the impact and fluid thresholds. This phenomenon is defined as transport hysteresis.

During the aeolian transport of particles, the fluid shear stress (or friction velocity) at the surface constantly acts on the deposited particles and plays a key role in aerodynamic entrainment and splash processes on the surface. The equilibrium property of the fluid shear stress at the surface during aeolian transport has been researched for decades. Owen (1964) first proposed this theory—universally known as Owen’s second hypothesis—from an aerodynamics’ perspective. According to this theory, if the friction velocity at the surface, $u_*^\mathrm{s}$, exceeds the impact threshold, the enhanced splash entrainment of particles reduces the wind speed and consequently reduces the friction velocity at the surface. Additionally, when $u_*^\mathrm{s} < u_*^\mathrm{i}$, the friction velocity at the surface increases with wind strengthening because of the absence of airborne particles. Therefore, Owen’s second hypothesis suggests that the friction velocity at the surface is equal to the impact threshold at the equilibrium state of aeolian transport: $u_*^\mathrm{s} = u_*^\mathrm{i}$.

Numerous geophysical models have been applied using Owen’s second hypothesis (Raupach 1991; Marticorena and Bergametti 1995; Sauermann et al. 2001; Ho et al. 2011), owing to its suitability for analytical approaches. Concurrently, numerical simulations without Owen’s second hypothesis have also been conducted (Werner 1990; Nemoto and Nishimura 2004; Kok and Renno 2009; Kok et al. 2012; Niiya and Nishimura 2017), and several researchers have challenged Owen’s second hypothesis. The friction velocity in the saltation layer decreases once the particles are mobilized due to a lesser relative velocity, and the friction velocity in the surface vicinity is less than the impact threshold (i.e., $u_*^\mathrm{s} \le u_*^\mathrm{i}$).

In wind-tunnel experiments and field observations, the friction velocity has been calculated using several methods, such as the profile method and eddy-correlation method (Bagnold 1941; Gerety 1985; Spies et al. 1995; Nemoto and Nishumura 2001; Nishimura et al. 2014; Aksamit and Pomeroy 2018). The vertical profile of the friction velocity above the height of a few centimetres can be obtained by measuring the wind velocity using an ultrasonic anemometer. Nishimura et al. (2014) measured the wind velocity at a height of 3 cm above the snow surface in wind-tunnel experiments. However, the ultrasonic anemometer is unsuitable for accurately measuring the wind velocity in the vicinity of a surface, owing to the path length of the ultrasonic anemometer and the effect of the high particle concentration.

Recently, Li and McKenna Neuman (2012) measured the wind velocity at a height of 3–270 mm during aeolian saltation with a laser Doppler anemometer. Their results show that the fluid shear stress decreases and approaches its impact threshold when the fluid moves towards the surface. Furthermore, Walter et al. (2014) measured the wind speed at a height of 5 mm during aeolian saltation using Irwin sensors (Irwin 1981). They reported that the friction velocity at the height of 5 mm decreased from the fluid threshold to the impact threshold when the wind speed above the saltation layer increased. It should be noted that the definition of impact threshold used in the experiments conducted by Walter et al. (2014) differs from the one in numerical simulations. The impact threshold was estimated as the minimum value of the friction velocity near the surface, which was determined in a short period after the onset of saltation, while the friction velocity above the saltation layer increased. Overall, these experimental results of aeolian saltation are evidence that the friction velocity in the first few millimetres above the surface exceeds the impact threshold (i.e., $u_* \ge u_*^\mathrm{i}$).

As mentioned above, the equilibrium property of the friction velocity (or fluid shear stress) in aeolian transport is inconsistently reported in the literature, numerical simulations, and experiments. This discrepancy, including the particle shear stress acting on the surface, has not yet been discussed in detail. In this study, we aim to investigate the equilibrium property of the fluid shear stress and its relationship with the particle shear stress in the surface vicinity using a random-flight model of drifting snow (Nemoto and Nishimura 2004; Niiya and Nishimura 2017). Sections 2 and 3 describe the governing equations of the random-flight model and the set-up of the numerical simulations, respectively. Section 4 shows the transport hysteresis and equilibrium properties of the surface friction velocity. Section 5 describes the mechanism of surface shear stresses based on the particle shear stress. Finally, Sect. 6 summarizes the study and presents the scope for future research.

2 Method

We simulated the aeolian transport of snow particles using a random-flight model of drifting snow (Nemoto and Nishimura 2004). The calculation process was reconstructed to describe the spatio–temporal transport structure with high accuracy (Niiya and Nishimura 2017). The model calculates the vertical profile of the wind speed and the two-dimensional trajectory of each particle in air. The aerodynamic entrainment and splash processes (i.e., particle–surface collision) were considered as sub-processes on the surface. Here, we present the governing equations of wind speed and particle trajectory. For further details, please refer to: Nemoto and Nishimura (2004) and Niiya and Nishimura (2017).

The vertical profile of the wind speed was reproduced using only the streamwise component of the time-averaged flow, ${\overline{u}}$. When assuming a horizontally uniform flow ($\partial /\partial x \approx 0$) and well-developed turbulent flow in the Reynolds-averaged Navier–Stokes equations, the advection and pressure terms can be removed. The governing equation of ${\overline{u}}$ is expressed as

$$\begin{aligned} \rho _\mathrm{f} \frac{\partial {\overline{u}}}{\partial t} = \frac{\partial \tau }{\partial z} - \sum ^{N} \frac{F_\mathrm{d}}{V_\mathrm{f}}, \end{aligned}$$

(1)

where $\rho _\mathrm{f}$ is the fluid density, t is the time, $\tau$ is the fluid shear stress, z is the height from the surface, and $V_\mathrm{f}$ is the volume of the fluid computational cell. Within each fluid computational cell, the number of particles and the streamwise component of the drag force acting on each particle are denoted N and $F_\mathrm{d}$, respectively. Thus, the second term on the right-hand side (r.h.s.) of Eq. 1 represents the momentum exchange between the wind and particles in air; however, the momentum exchange on the surface occurs due to the aerodynamic entrainment. In this study, the fluid shear stress is described using Prandtl’s mixing length theory and molecular viscosity

$$\begin{aligned} \tau = \rho _\mathrm{f} \left( \kappa z \frac{\partial {\overline{u}}}{\partial z} \right) ^2 + \mu \frac{\partial {\overline{u}}}{\partial z}, \end{aligned}$$

(2)

where $\kappa$ is the von Kármán constant (0.4) and $\mu$ is the dynamic viscosity. The first and second terms on the r.h.s. of Eq. 2 are the Reynolds and viscous stresses, $\tau _\mathrm{R}$ and $\tau _\mathrm{v}$, respectively. In particular, the viscous stress near the surface exceeded the Reynolds stress.

The vertical component of turbulent fluctuation, $w'$, is considered for each particle; however, the vertical component of the time-averaged flow is assumed to be zero: ${\overline{w}} = 0$. According to Wilson and Sawford (1996), the governing equation of $w'$ in the case of uniform-isotropic turbulence can be expressed using Euler statistics and the Kolmogorov similarity law

$$\begin{aligned} w'(t + \delta t) = \left( 1 - \frac{\delta t}{T_*^\mathrm{L}} \right) w'(t) + \sigma _w \sqrt{\frac{2 \delta t}{T_*^\mathrm{L}}} \eta (t), \end{aligned}$$

(3)

where $\delta t$ is the timestep of the numerical simulation, $T_*^\mathrm{L}$ is the Lagrangian time scale considering the particle inertia, $\sigma _w$ is the turbulence intensity, and $\eta (t)$ is a random number generated by the standard normal distribution. According to Hunt and Nalpanis (1985), $T_*^\mathrm{L}$ in the boundary layer of a neutral atmosphere is given as a functions of the friction velocity ($u_* = \sqrt{\tau /\rho _\mathrm{f}}$), the relative speed between the wind and particle ($V_\mathrm{R}$), and the height of the particle. In addition, $\sigma _w$ is given as a function of the friction velocity: ($\sigma _w = 1.3 u_*$). The initial condition of $w'$ in Eq. 3 is zero when the particle is subjected to aerodynamic entrainment and splash processes.

The trajectories of the airborne particles are calculated in the streamwise and vertical directions (i.e., the x–z plane). The particle is assumed to be an irrotational hard-spherical grain, and the only forces acting on the particle were gravity and drag forces. The governing equations of the particle velocity, $(v_x, v_z)$, are expressed using Newton’s equation of motion

$$\begin{aligned} m \frac{dv_x}{dt}&= \frac{\rho _\mathrm{f} S C_\mathrm{d}}{2} V_\mathrm{R} \left( {\overline{u}} - v_x \right) , \end{aligned}$$

(4a)

$$\begin{aligned} m \frac{dv_z}{dt}&= \frac{\rho _\mathrm{f} S C_\mathrm{d}}{2} V_\mathrm{R} \left( w' - v_z \right) - m g, \end{aligned}$$

(4b)

where m and S are the mass and cross-sectional area of the particle, respectively; $C_\mathrm{d}$ is the drag coefficient defined as a function of the particle Reynolds number (White 1974); and g is the acceleration due to gravity. The mass and cross-sectional area of the particle are calculated using the particle diameter and particle density (d and $\rho _\mathrm{p}$) as $m = \rho _\mathrm{p} \pi d^3/6$ and $S = \pi d^2/4$, respectively. Here, ${\overline{u}}$ in Eq. 4 is calculated at the fluid computational cell to include the particle. The r.h.s. of Eq. 4a is equal to $F_\mathrm{d}$ in Eq. 1.

Aerodynamic entrainment occurs when the friction velocity at the surface $u_*^\mathrm{s}$ exceeds the fluid threshold $u_*^\mathrm{f}$. The number of wind-entrained particles is calculated based on wind-tunnel experiments (Shao and Li 1999). The streamwise component of the initial particle velocity is set as half the value of the friction velocity at the surface as $v_x = u_*^\mathrm{s}/2$, whereas the vertical component is set as $v_z = \sqrt{2 g d}$, so that the particle reaches a height equal to the diameter of wind-entrained particle. Splash processes occurs when the airborne particle collides with the surface. The particle dynamics are classified into three types (Fig. 1): ‘Deposition’ and ‘Rebound’ refer to the stop and rebound of the incident particle, respectively, and ‘Splash’ denotes particle ejection from the surface (i.e., surface erosion). These particle dynamics were simulated using statistical functions measured in wind-tunnel experiments of drifting snow, conducted by Sugiura and Maeno (2000); however, the value of $u_*$ was low. For descriptions of both these sub-processes, please refer to Niiya and Nishimura (2017).

For the splash functions, we apply the experimental results for disintegrated particles of natural compact snow that were obtained by Sugiura and Maeno (2000). The fragmentation of snow particles due to the impact was negligible, and the snow bed was loose. The mean and standard deviation of disintegrated particle diameters of natural compact snow were 0.35 mm and 0.28 mm, respectively. By utilizing the incident angle and speed ($\theta _\mathrm{i}$ and $v_\mathrm{i}$) in the splash functions, three values were calculated: number of ejected particles ($n_\mathrm{e}$) including the rebounded particle, and horizontal and vertical restitution coefficients ($e_\mathrm{h}$ and $e_\mathrm{v}$, respectively). The distributions of $n_\mathrm{e}$, $e_\mathrm{h}$, and $e_\mathrm{v}$ are expressed by binomial, normal, and gamma distributions, respectively. In this study, splash processes on a loose snow bed is reproduced.

Fig. 1

Schematic of the particle shear stress, $\tau _\mathrm{p}$, generated by splash processes. The streamwise velocity component of the incident and ejected particles is denoted by $v_x$

3 Numerical Setting

The numerical simulations of the random-flight model were conducted in a quasi-two-dimensional system, and the domain was a cuboid with a length, width, and height of 0.1, 0.01, and 1 m, respectively. We simulated the fluid computational cells to calculate the vertical profile of the streamwise velocity component ${\overline{u}}$ in Eq. 1. The vertical boundary between cells and the vertical cell size ($z_i$ and $c_i$, respectively) were given as

$$\begin{aligned} z_i&= 10^{-5 + {i}/{10}}, \end{aligned}$$

(5a)

$$\begin{aligned} c_i&= z_i - z_{i-1}, \end{aligned}$$

(5b)

where i is the cell number ranging from 0 to 51. In Eq. 5, $z_0$ represents the roughness length (basically, the surface asperity), which was fixed during the simulation, because the bottom of the domain was assumed to be a flat and loose snow bed. Hence, the no-slip boundary condition at the bottom of the domain was used; the streamwise velocity component below the roughness length is given as zero: ${\overline{u}}= 0$ at $z \le z_0$. The wind profile within the domain was determined by the boundary condition at the top of the domain, where the change in friction velocity ($u_*^\mathrm{t}$) was calculated stepwise to simulate the wind-tunnel experiments by Walter et al. (2014). The streamwise boundary condition for the particles was set to be periodic, whereas the particle exceeding the top of the domain was removed from the simulation.

Typical values for the air and snow particles were used in the simulation to represent specific materials for the fluid and particles. The fluid density and viscosity coefficient were fixed as $\rho _\mathrm{f}= 1.2\,\hbox {kg}\,\hbox {m}^{-3}$ and $\mu = 1.82 \times 10^{-5}$ Pa s, respectively. The particle density was given as $\rho _\mathrm{p}= 910\,\hbox {kg}\,\hbox {m}^{-3}$. According to observations (Schmidt 1982, 1984; Nishimura and Nemoto 2005) and experiments (Gromke et al. 2014; Nishimura et al. 2014) on blowing snow, the particle size distribution is well fitted by the gamma distribution. Thus, in the simulation, we assumed that the particle diameter of the loose snow bed adhered to the gamma distribution: $d \in \varGamma (3, 0.1)$ [units: mm]. Here, the mean and peak ${\overline{d}}$ are 0.3 and 0.2 mm, respectively, and the simulated particle diameters range from 0.01 to 1 mm. According to Bagnold (1941), the roughness length is approximately the particle diameter divided by 30; thus, we used $z_0 = {\overline{d}}/30 = 10^{-5}$ m. The fluid threshold is significantly related to the particle diameter. However, by examining the experimental data of Clifton et al. (2006), we assumed that $u_*^\mathrm{f}= 0.24\,\hbox {m\,s}^{-1}$ corresponded to the peak of the particle diameter (0.2 mm).

To study the transport hysteresis and properties of the surface shear stress, we varied the friction velocity at the top ($u_*^\mathrm{t}$) with time. The initial condition of $u_*^\mathrm{t}$ was set to $0.15\,\hbox {m\,s}^{-1}$ below the fluid threshold, and the initial wind profile was considered as a constant stress layer. Then, $u_*^\mathrm{t}$ increased by $0.01\,\hbox {m\,s}^{-1}$ per $T = 200$ s, which is the time required to reach the equilibrium state, according to the simulation conducted by Niiya and Nishimura (2017). After reaching $u_*^\mathrm{t}= 0.30\,\hbox {m\,s}^{-1}$, we decreased $u_*^\mathrm{t}$ to $0.15\,\hbox {m\,s}^{-1}$ at a rate of 0.01 $\hbox {m\,s}^{-1}$ per T, as illustrated by the green line in Fig. 2. This setting was limited to the condition of low friction velocities because of the range associated with the empirical formulae for splash processes (Sugiura and Maeno 2000). In the wind and particle calculations, the timestep was set as $\delta t= 10 ^{-7}$ s to avoid numerical instability.

4 Results

The results are displayed in three subsections: Sect. 4.1 presents the time-series data in no-transport and transport states of particles, Sect. 4.2 presents the vertical profiles of wind speed and friction velocity at the equilibrium state, and Sect. 4.3 presents the relationship between the surface friction velocity and the particle transport rate.

4.1 Time Evolution of Transport

The single indicator of the transport state was the particle transport rate, Q, which is the integral of the particle streamwise mass flux $q_x$ from the surface to the top boundary

$$\begin{aligned} Q = \int _{0\,{\text {m}}}^{1\,\mathrm{m}} q_x dz. \end{aligned}$$

(6)

In the simulation, $q_x$ from Eq. 6 was calculated using the particle mass passing through the streamwise boundary ($x= 0$ and 0.1 m). The fluid shear stress acting on the surface was evaluated using the friction velocity at the roughness length ($u_*^\mathrm{s}$) because the wind speed was zero below the roughness length. Figure 2 shows the time evolutions of the friction velocities at the top boundary and roughness length ($u_*^\mathrm{t}$ and $u_*^\mathrm{s}$, respectively) as well as the particle transport rate (Q). As mentioned in Sect. 3, the change in $u_*^\mathrm{t}$ per $T= 200$ s appears to be a step function. Moreover, the initial value of $u_*^\mathrm{s}$ is equal to $u_*^\mathrm{t}$, because the wind profile at $t= 0$ was set as the constant stress layer. Although we removed the particles exceeding the top boundary, this was not implemented during the simulation. That is, the saltation was dominant in the particle dynamics.

In the early stage of the simulation ($t < 10T$ in Fig. 2a), the particles could not be entrained by the wind because the friction velocity at the roughness length was lower than the fluid threshold ($u_*^\mathrm{s} \le u_*^\mathrm{f}= 0.24\,\hbox {m\,s}^{-1}$). Therefore, the particle transport rate is zero: $Q= 0$. This simulation, wherein no airborne particles were present, resulted in the fluid momentum transfer from the top boundary to the surface only. The friction velocity at the roughness length became similar to that at the top boundary with time and converged before the increase in $u_*^\mathrm{t}$: $u_*^\mathrm{s} \approx u_*^\mathrm{t}$.

Fig. 2

Time evolutions of friction velocities at the top boundary and roughness length ($u_*^\mathrm{t}$ and $u_*^\mathrm{s}$, respectively) as well as particle transport rate (Q): a entire time, b start of transport, and c end of transport. The friction velocity at the top boundary changes at the time interval of $T = 200$ s. The fluid threshold is given in simulations as $u_*^\mathrm{f}= 0.24\,\hbox {m\,s}^{-1}$, whereas the impact threshold is estimated from Eq. 7 as $u_*^\mathrm{i} \approx 0.206\,\hbox {m\,s}^{-1}$

After $t= 10T$, the friction velocity at the roughness length exceeds the fluid threshold, and the subsequent aerodynamic entrainment acted as a trigger for the particle transport (Fig. 2b). The particle transport rate increased significantly and then decreased slowly, whereas the friction velocity at the roughness length decreased rapidly and then remained less than the fluid threshold: $u_*^\mathrm{s} < u_*^\mathrm{f}$. According to the previous study on the random-flight model (Niiya and Nishimura 2017), the rapid decrease in $u_*^\mathrm{s}$ at the onset of saltation was caused by both the increase in the number of saltation particles and the rise in the mean saltation height. The development of saltation reduced the wind speed and friction velocity near the surface. The time required for decreasing $u_*^\mathrm{s}$ after the onset of saltation is much shorter in simulations ($\approx$ 5 s) than that in wind-tunnel experiments ($\approx$ 50 s) conducted by Walter et al. (2014). One of the possible reasons for this difference is that the saltation layer was not well developed in their experiments because of the fetch of 7 m from the measurement point of the friction velocity. Since the periodic boundary for the particles was used in our simulations, the long fetch was reproduced.

In the middle stage of the simulation ($10T \le t \le 25T$ in Fig. 2a), Q is positively correlated to $u_*^\mathrm{t}$. However, the correlation between Q and $u_*^\mathrm{s}$ is not clear. Transport hysteresis was confirmed in $22T \le t \le 25T$, where particle transport continued at the top boundary condition below the fluid threshold. The stepwise change of $u_*^\mathrm{t}$ is inspired by wind-tunnel experiments by Walter et al. (2014), but transport hysteresis was not observed in their experiments; the saltation stopped at the same level of friction velocity as that observed at the onset of saltation. The saltation process was triggered by the aerodynamic entrainment on loose sand (Walter et al. 2014); thus, the saltation layer could not be formed at the friction velocity below the fluid threshold. In the wind tunnel, the particle-supplying device should measure the transport hysteresis.

At $t= 25 T$, the friction velocity at the top boundary was lower than that at the roughness length, thereby terminating the particle transport (Fig. 2c). In the final stage of the simulation ($25T \le t$ in Fig. 2a), $u_*^\mathrm{s}$ and $u_*^\mathrm{t}$ converge. This case is similar to that in the early stage ($t < 10T$); however, the fluid momentum was released from the top boundary.

4.2 Wind Profile at Equilibrium State

Figure 3a, b shows the instantaneous vertical profiles of the streamwise velocity component ${\overline{u}}$ at the equilibrium state in Fig. 2. The dashed and solid lines denote the no-transport and transport states of the particles, respectively. We analyzed the effect of viscous stress on the wind speed because the fluid shear stress was calculated using the Reynolds and viscous stresses, $\tau _\mathrm{R}$ and $\tau _\mathrm{v}$, respectively, in Eq. 2. Considering only the Reynolds stress, the wind speed exhibits the relationship ${\overline{u}} \propto \ln z$. A typical case is denoted by the black line in Fig. 3a. Considering the viscous stress, the shape of the vertical profile of ${\overline{u}}$ changed below $z\approx 10^{-3}$ m. Therefore, the wind speed decreased completely (indicated by dashed lines in Fig. 3a, b). We also confirmed that the wind speed in the transport state does not significantly vary in the surface vicinity (indicated by the solid lines in Fig. 3a, b).

Figure 3c shows ratios of wind speeds ${\overline{u}}/{\overline{u}}_{t = 16T}$ in the transport state, wherein the denominator is the instantaneous streamwise velocity component at $t= 16 T$. The ratio of wind speeds is less than unity above the height of 1 cm, as shown in Fig. 3a, b; however, ${\overline{u}}/{\overline{u}}_{t = 16T}$ is close to unity at z = 3–6 mm. This implies that vertical profiles of the wind speed intersect at one point, which is known as the focus point observed by Bagnold (1941). According to Kok et al. (2012), the occurrence of a focus point is related to the decrease in the surface fluid shear stress that occurs during saltation as the friction velocity increases.

Figure 3d shows streamwise velocity components of wind and particle, ${\overline{u}}$ and $v_x$, as the average of the values calculated every 10 s at each fluid computational cell at $u_*^\mathrm{t}= 0.30\,\hbox {m\,s}^{-1}$. Considering the difference in wind and particle speeds, particles can be accelerated by winds above the height of 0.4 mm, whereas particles decelerate near the surface. In particular, the average particle speed tends to zero at $z < 70\, \upmu \hbox {m}$, because smaller particles are affected by the wind. The saltation-layer height at $u_*^\mathrm{t}= 0.30\,\hbox {m\,s}^{-1}$ is estimated from the vertical profile of the average particle speeds as 8 cm (indicated by the horizontal lines in Fig. 3).

Figure 4a, b, c indicates the vertical profiles of the friction velocity $u_*$ at the equilibrium state in Fig. 2, and the data obtained are the average values obtained 10 s before the friction velocity at the top boundary ($u_*^\mathrm{t}$) changes. The dashed and solid lines denote the no-transport and transport states of the particles, respectively. Figure 4d shows the vertical profiles of the Reynolds and viscous stresses ($\tau _\mathrm{R}$ and $\tau _\mathrm{v}$, respectively) in the no-transport state. In the no-transport state, the friction velocity is uniform; therefore, a constant stress layer is formed as per the definition of $\tau = \rho _\mathrm{f} u_*^2$. Note that the contributions of $\tau _\mathrm{R}$ and $\tau _\mathrm{v}$ to $\tau$ vary depending on the height. The viscous stress is dominant in the vicinity of the surface but is smaller than the Reynolds stress in higher regions (Fig. 4d). In the transport state, the friction velocity is not uniform and decreases near the surface, as shown in Fig. 4a, b, c. The decrease in $u_*$ is caused by the momentum transfer from the wind to the airborne particles (Fig. 4c).

Fig. 3

Instantaneous vertical profiles of the streamwise velocity component, ${\overline{u}}$: a increase and b decrease phases of $u_*^\mathrm{t}$ in Fig. 2; c ratios of speeds based on ${\overline{u}}$ at $t= 16 T$ ($u_*^\mathrm{t}= 0.30\,\hbox {m\,s}^{-1}$), and d ${\overline{u}}$ and particle speed $v_x$ at $u_*^\mathrm{t}= 0.30\,\hbox {m\,s}^{-1}$. The dashed and solid lines denote no-transport and transport states of particles, respectively. In d, particle speeds are averaged over 10 s. The saltation layer height at $u_*^\mathrm{t}= 0.30\,\hbox {m\,s}^{-1}$ is denoted by the horizontal line and determined from the vertical profile of $v_x$ as 8 cm

Fig. 4

Averaged vertical profiles of friction velocity and fluid shear stress ($u_*$ and $\tau$, respectively): a increase and b decrease phases of $u_*^\mathrm{t}$ in Fig. 2; c comparison of friction velocities at the same $u_*^\mathrm{t}$; d Reynolds and viscous stresses ($\tau _\mathrm{R}$ and $\tau _\mathrm{v}$, respectively). In a–c, the dashed and solid lines denote the no-transport and transport states, respectively. In d, the data for the no-transport state are plotted

4.3 Particle Transport Rate and Surface Friction Velocity

Transport hysteresis was confirmed in the time evolution of the particle transport rate Q (Fig. 2). Here, we estimate the impact threshold $u_*^\mathrm{i}$ from the dependence of the particle transport rate on the friction velocity at the top boundary $u_*^\mathrm{t}$. Figure 5a shows the equilibrium property of the particle transport rate, where the point and error bar indicate the average and standard deviation of 10 s, respectively. In the increase phase of $u_*^\mathrm{t}$, as illustrated by the blue points in Fig. 5a, the particle transport rate exceeds zero when $u_*^\mathrm{t} > u_*^\mathrm{f}$ ($= 0.24\,\hbox {m\,s}^{-1}$). The particle transport rate in the decrease phase of $u_*^\mathrm{t}$, as illustrated by the red points in Fig. 5a, is approximately the same as that in the increase phase; however, it exceeds zero below the fluid threshold. From this trend of Q, the particle transport rate was expected to be zero in the region of $u_*^\mathrm{t}= 0.20{-}0.21\,\hbox {m\,s}^{-1}$.

The fitting curve for Q in Fig. 5a was obtained using the least-squares method and the data of $Q > 0$

$$\begin{aligned} Q_{u_*^\mathrm{t} \le 0.30 \text {m s}^{-1}} = 0.0882 (u_*^\mathrm{t} - 0.206)^{1.15}. \end{aligned}$$

(7)

Generally, the particle transport rate of the saltation is proportional to the square and the cube of $u_*$ (Kok et al. 2012), and the scaling of $Q \propto u_*^{2.35}$ was determined by simulations of the random-flight model including higher friction velocities (Niiya and Nishimura 2017). We only use the data of lower friction velocities; thus, the power index of $u_*$ in Eq. 7 has a weak nonlinearity of 1.15. Here, we assume that the impact threshold equals the friction velocity corresponding to $Q= 0$ in Eq. 7: $u_*^\mathrm{i} \approx 0.206\,\hbox {m\,s}^{-1}$. From the least-squares method, the standard error of $u_*^\mathrm{i}$ is $0.001\,\hbox {m\,s}^{-1}$, as illustrated by the background colour (red) in Fig. 5. According to Bagnold (1941), the ratio between the impact and fluid thresholds is approximately 0.85; thus, our result ($u_*^\mathrm{i}/u_*^\mathrm{f}= 0.858$) is consistent with that of previous studies.

We studied the friction velocity at the roughness length $u_*^\mathrm{s}$, which was compared to the estimated impact threshold $u_*^\mathrm{i}$. Figure 5b shows the equilibrium property of the surface friction velocity, in which the point indicates the average value within 10 s. In the no-transport state ($Q= 0$ in Fig. 5), the friction velocity at the roughness length is equal to that at the top boundary because a constant stress layer was formed, as shown in Fig. 4. In the transport state ($Q > 0$ in Fig. 5a), the friction velocity at the roughness length decreases slightly with an increase in $u_*^\mathrm{t}$ around the impact threshold.

To further investigate the change in $u_*^\mathrm{s}$ in the transport state, we analyzed the variability. Figure 6 shows the box plots of $u_*^\mathrm{s}$, where the boxes indicate the lower and upper quartiles, the lines within the boxes indicate the median, and the error bars denote the minimum and maximum data points; outliers are indicated by points on the graph. In Fig. 6, outliers are defined as the data outside a closed interval that ranges from the lower quartile minus 1.5 times the interquartile range (IQR) to the upper quartile plus 1.5 times the IQR (where IQR is the distance between the lower and upper quartiles). Medians of $u_*^\mathrm{s}$ in $u_*^\mathrm{t} \le u_*^\mathrm{f}$ remain within the standard error of $u_*^\mathrm{i}$, as illustrated by the background colour (red) in Fig. 6, whereas those in $u_*^\mathrm{t} > u_*^\mathrm{f}$ decrease marginally and are located outside the standard error of $u_*^\mathrm{i}$. However, this decrease rate of $u_*^\mathrm{s}$ is quite small in comparison to the numerical simulations of blown sand and snow (Werner 1990; Nemoto and Nishimura 2004; Kok and Renno 2009; Kok et al. 2012; Niiya and Nishimura 2017). Therefore, our results of low friction velocities do not contradict Owen’s second hypothesis.

Fig. 5

Hysteresis in particle transport: a particle transport rate Q and b friction velocity at the surface roughness length $u_*^\mathrm{s}$. The periods of $t \le 16T$ and $t > 16T$ represent the increase and decrease phases of $u_*^\mathrm{t}$ in Fig. 2. The dashed and dotted lines indicate the fluid and impact thresholds ($u_*^\mathrm{f}= 0.24\,\hbox {m\,s}^{-1}$ and $u_*^\mathrm{i} \approx 0.206\,\hbox {m\,s}^{-1}$, respectively). The red colour in the background indicates the standard error of $u_*^\mathrm{i}$ in Eq. 7. In a, the error bar indicates the standard deviation of Q

Fig. 6

Box plots of the friction velocity at the roughness length, $u_*^\mathrm{s}$, in the particle transport state in Fig. 5b. The red colour at the background indicates the standard error of $u_*^\mathrm{i}$ in Eq. 7

5 Discussion

In this study, splash processes of the random-flight model was simulated using splash functions measured on the loose snow bed at low friction velocities (Sugiura and Maeno 2000). According to Niiya and Nishimura (2017), splash functions proposed by Sugiura and Maeno (2000) underestimated the saltation height, because the mean of the vertical restitution coefficient was observed to be less than unity at an increasing incident speed. Although the independent variables of splash functions were the incident angle and speed, the particle diameter was not considered. Despite the above-mentioned limited ability of the model to simulate splash processes, we discuss the equilibrium properties of the particle transport near the surface. It is noted that other splash functions have been proposed to overcome the limitation (Ammi et al. 2009; Comola and Lehning 2017), and splash processes enhanced by the snow sintering have been simulated by the discrete-element method (Comola et al. 2019).

5.1 Contribution of Particle Shear Stress to Total Shear Stress

Momentum transfer from wind to the airborne particles occurs in the transport state (Fig. 4c). Therefore, the particle shear stress $\tau _\mathrm{p}$ is essential for understanding the change in friction velocity at the roughness length, $u_*^\mathrm{s}$, as shown in Fig. 6. According to Raupach (1991), the total shear stress in the steady state $\tau _\mathrm{total}$ can be divided into the fluid shear stress $\tau$ and particle shear stress: $\tau _\mathrm{total} = \tau + \tau _\mathrm{p}$. Here, $\tau$ is calculated using Eq. 2, whereas $\tau _\mathrm{p}$ is estimated using particle dynamics. The specific form of $\tau _\mathrm{p}$ is given by Shao (2008)

$$\begin{aligned} \tau _\mathrm{p} = \sum ^{n\downarrow }m v_x - \sum ^{n\uparrow }m v_x, \end{aligned}$$

(8)

where $n{\downarrow }$ and $n{\uparrow }$ are the number of descending and ascending particles passing a certain height per unit area and unit time, respectively, while m and $v_x$ are the mass and streamwise velocity component of the particle, respectively. In this study, we calculated the particle shear stress at the boundary of each fluid computational cell by using all particle trajectories observed 10 s before the friction velocity changed at the top boundary ($u_*^\mathrm{t}$). However, we could not correctly estimate the particle shear stress near the surface because of the collision of particles into the surface (i.e., splash processes).

Figure 7 shows the vertical profiles of the fluid and particle shear stresses ($\tau$ and $\tau _\mathrm{p}$, respectively) obtained by analyzing two $u_*^\mathrm{t}$ values (0.23 and $0.30\,\hbox {m\,s}^{-1}$). The vertical dashed line indicates the fluid shear stress at the top boundary, which was converted as $\tau _\mathrm{t} = \rho _\mathrm{f} (u_*^\mathrm{t})^2$. Meanwhile, the horizontal dotted line indicates the mean particle diameter in the snow bed (${\overline{d}}= 0.3$ mm). Similar to the friction velocity, the fluid shear stress decreases near the surface, as shown in Fig. 4; however, the particle shear stress increases near the surface (Fig. 7). Moreover, the total shear stress is uniform and equalled the fluid shear stress at the top boundary, $\tau _\mathrm{total} = \tau _\mathrm{t}$. A constant stress layer is formed in the transport state, except near the surface. At $z < {\overline{d}}$, the particle shear stress decreases towards the surface because the particle coordinates were not located in fluid computational cells of the nearby surface.

To evaluate the particle shear stress acting on the surface, we considered splash processes, as shown in Fig. 1. In particular, the first term on the r.h.s. of Eq. 8 was calculated using the incident particle, and the second term was computed from the ejected particle as rebound and splash. Figure 8 shows the dependence of the fluid and particle shear stresses at the surface ($\tau$ and $\tau _\mathrm{p}$, respectively) on the fluid shear stress at the top boundary ($\tau _\mathrm{t}$), for which the data of $t > 15T$ in Fig. 2 are used. The fluid and impact thresholds converted into shear stress ($\tau _\mathrm{f}$ and $\tau _\mathrm{i}$) are denoted by the dashed and dotted vertical lines in Fig. 8. The total shear stress at the surface equals the fluid shear stress at the top boundary, giving a constant stress layer from the surface to the top. In the transport state of $\tau _\mathrm{t} >\tau _\mathrm{i}$, the particle shear stress marginally exceeds $\tau _\mathrm{t} - \tau _\mathrm{i}$, as illustrated by the dashed-dotted line in Fig. 8. If it is held in $\tau _\mathrm{p} = \tau _\mathrm{t} - \tau _\mathrm{i}$, then the fluid shear stress remains constant. Therefore, the result of $\tau _\mathrm{p} \ge \tau _\mathrm{t} - \tau _\mathrm{i}$ leads to a friction velocity at a roughness length smaller than the impact threshold, as shown in Fig. 6.

Fig. 7

Vertical profiles of fluid and particle shear stresses, $\tau$ and $\tau _\mathrm{p}$, respectively, in the particle transport state. The horizontal dotted line indicates the mean particle diameter in the snow bed, ${\overline{d}}$, and the vertical dashed line indicates the fluid shear stress fixed at the top boundary, $\tau _\mathrm{t}$

Fig. 8

Dependence of fluid and particle shear stresses at the surface ($\tau$ and $\tau _\mathrm{p}$, respectively) on the fluid shear stress at the top boundary ($\tau _\mathrm{t}$). The dashed and dotted vertical lines indicate the fluid and impact thresholds converted into shear stress ($\tau _\mathrm{f}$ and $\tau _\mathrm{i}$, respectively)

5.2 Change in Splash Processes

The particle shear stress at the surface was determined using splash processes, as shown in Fig. 1, because the particle shear stress was calculated using Eq. 8. To reveal the mechanism governing the increase in particle shear stress at the surface, as shown in Fig. 8, we analyzed the variables used to calculate the particle shear stress, $\tau _\mathrm{p}$, which were the numbers of incident and ejected particles per unit area and unit time ($n{\downarrow }$ and $n{\uparrow }$, respectively), particle mass or diameter (m or d, respectively), and streamwise velocity component $v_x$. Figure 9 indicates the relationships between these variables ($n{\downarrow }$, $n{\uparrow }$, d, $v_x$) and the fluid shear stress at the top boundary, $\tau _\mathrm{t}$. Here, the horizontal axis, $\tau _\mathrm{t} - \tau _\mathrm{i}$, was chosen to evaluate the transport state of $\tau _\mathrm{p} > 0$ in Fig. 8. The point and error bars represent the mean and standard deviation values determined for 10 s for the incident or ejected particles in the equilibrium state, respectively.

In Fig. 9a, the number of incident particles per unit area and unit time is equal to that of ejected particles, because the measurement was conducted in the equilibrium state. Moreover, a positive linear relation was determined between $n{\downarrow }$ (and $n{\uparrow }$) and $\tau _\mathrm{t} - \tau _\mathrm{i}$, as illustrated by the dashed line in Fig. 9a: $n{\downarrow } = n{\uparrow } \approx 76.3 \times 10^7$ ($\tau _\mathrm{t} - \tau _\mathrm{i}$). Most of the ejected particles represent the rebound process; however, approximately 13.8% of these particles correspond to splash processes. As shown in Fig. 9b, the mean diameter of the ejected particles is equal to that of the incident particles, because the saltation is dominant on the particle dynamics and ejected particles that collided with the surface. Additionally, the mean diameter was 0.151 mm and independent of $\tau _\mathrm{t} - \tau _\mathrm{i}$, but the standard deviation was significant. The particle mass, which was calculated from the mean diameter of 0.151 mm, corresponds to m = 1.64 $\times$ 10$^{-9}$ kg. In Fig. 9c, the mean streamwise velocity component of the ejected particles is lower than that of the incident particles, owing to the momentum loss caused by the collision. Splash processes on the loose snow bed are incorporated using the experimental results of Sugiura and Maeno (2000), in which the horizontal momentum of the incident particle is largely transformed into the vertical momentum of the ejected particle. Although the horizontal restitution coefficient $e_\mathrm{h}$ of the rebound particles depended on incident angle and speed, the mean of $e_\mathrm{h}$ was approximately 0.5. The mean streamwise velocity components of the incident and ejected particles in Fig. 9c are approximately 0.957 $\hbox {m\,s}^{-1}$ and $0.441\,\hbox {m\,s}^{-1}$, respectively, and independent of $\tau _\mathrm{t} - \tau _\mathrm{i}$. Such constant particle speeds near the surface have been observed in wind-tunnel experiments (Rasmussen and Sørensen 2008; Creyssels et al. 2009; Ho et al. 2011) and a numerical model (Kok and Renno 2009).

Based on the above result, we determine the particle shear stress at the surface in a mean-field approximation as

$$\begin{aligned} n{\downarrow } = n{\uparrow }&\approx 76.3 \times 10^7 (\tau _\mathrm{t} - \tau _\mathrm{i}), \end{aligned}$$

(9a)

$$\begin{aligned} m{\downarrow } = m{\uparrow }&\approx {\text {const.}}~(1.64 \times 10^{-9}~\text {kg}), \end{aligned}$$

(9b)

$$\begin{aligned} v_x{\downarrow } - v_x{\uparrow }&\approx {\text {const.}}~(0.516\,{\text {ms}}^{-1}), \end{aligned}$$

(9c)

where ($m{\downarrow }$, $v_x{\downarrow }$) and ($m{\uparrow }$, $v_x{\uparrow }$) are the mean mass and mean streamwise velocity component of the incident and ejected particles on the surface, respectively. By substituting Eq. 9 into Eq. 8, we obtain the relational expression

$$\begin{aligned} \tau _\mathrm{p} = n{\downarrow } \,m{\downarrow } \left( v_x{\downarrow } - v_x{\uparrow } \right) \approx 0.646 (\tau _\mathrm{t} - \tau _\mathrm{i}). \end{aligned}$$

(10)

Equation 10 indicates that the particle shear stress at the surface is proportional to the difference between the fluid shear stress at the top boundary and the impact threshold converted into shear stress. The linear relationship between $\tau _\mathrm{p}$ and $\tau _\mathrm{t} - \tau _\mathrm{i}$ was also confirmed by the data of $\tau _\mathrm{p}$ as shown in Fig. 8; therefore, the change in the number of airborne particles is a major factor in the particle shear stress at the surface. Equation 10, however, results in the underestimation of $\tau _\mathrm{p}$, as illustrated by the black solid line in Fig. 8; the dependence of $v_x$ on m (or d) and the wide distributions of m and $v_x$ are thought to cause this underestimation.

Fig. 9

Changes in the variables used to calculate particle shear stress, $\tau _\mathrm{p}$: a number of incident and ejected particles per unit area and unit time ($n{\downarrow }$ and $n{\uparrow }$, respectively), b particle diameter d, and c streamwise velocity component $v_x$. The horizontal axis, $\tau _\mathrm{t} - \tau _\mathrm{i}$, is the difference between shear stresses at the top boundary and the impact threshold in Fig. 8. In a, the number of ejected particles is the sum of the rebound and splash in Fig. 1

5.3 Distributions of Diameter and Streamwise Particle Velocity

In splash processes, the mean diameter and mean streamwise velocity component of the particles are unchanged, independent of the wind condition, as shown in Fig. 9b, c. However, the standard deviations are significant. Hence, the distributions of the diameter and streamwise velocity component should be analyzed to understand their effect on the particle shear stress at the surface.

Fig. 10

Probability density distributions of the particle diameter and streamwise velocity component, d and $v_x$, respectively, for splash processes: a, c incident particles and b, d ejected particles. In a, b, the dashed line is the particle-diameter distribution of the loose snow bed and the gamma distribution, given in simulations

Figure 10 indicates the probability density distributions of the diameter d and streamwise velocity component $v_x$ of the incident and ejected particles at four different friction velocities at the top boundary: $u_*^\mathrm{t} = 0.21$, 0.24, 0.27, and $0.30\,\hbox {m\,s}^{-1}$. These distributions were obtained from the data at 10 s for the incident and the ejected particles in the equilibrium state, displaying a single-peak pattern. The distribution at $u_*^\mathrm{t}= 0.21\,\hbox {m\,s}^{-1}$, as illustrated by the red line in Fig. 10, is not smooth, because of the small number of particles. In Fig. 10a, the distribution of the diameter of incident particles ranged from $d= 0.01$ mm to $d= 0.5$ mm and appeared to be independent of the friction velocity at the top boundary. Additionally, the diameter distributions of the ejected particles are the same as those of the incident particles, as shown in Fig. 10b, as rebounded and splashed particles collided with the surface again. These diameter distributions of incident and ejected particles deviate from the size distribution of the loose snow bed, as illustrated by the dashed line in Fig. 10a, b. Although the diameter of splashed particle was selected from the bed-size distribution, the larger particle was not well accelerated by the wind, and the incident velocity was insufficient to rebound on the surface. Thus, smaller particles account for a large percentage of diameter distributions. In Fig. 10c, the distribution of the streamwise velocity component of the incident particles ranges from $v_x = -0.1\,\hbox {m\,s}^{-1}$ to $v_x = 2.5\,\hbox {m\,s}^{-1}$. Moreover, the probability density of the lower and higher $v_x$ increases with an increase in the friction velocity at the top boundary. In contrast, the velocity distributions of the ejected particles are not dependent on the friction velocity at the top boundary, as shown in Fig. 10d, and shift to a lower velocity than that of the incident particles. The dynamics of the ejected particles were simulated using statistical functions, which were obtained by Sugiura and Maeno (2000). Based on these simulations, it was determined that the ejected particle distribution is unlikely to be influenced by minor changes of incident angle and speed.

To indicate the minor changes in the diameter and streamwise velocity component of the incident particles, we calculated the difference in the probability density. Figure 11 indicates the differences in the probability density of the incident particles, as shown in Fig. 10a, c. Here, we use the diameter and streamwise velocity component distributions at $u_*^\mathrm{t}= 0.24\,\hbox {m\,s}^{-1}$ as the baseline, because the number of particles at $u_*^\mathrm{t}= 0.21\,\hbox {m\,s}^{-1}$ is relatively small. In Fig. 11a, the four vertical lines represent examples of the correlation between particle mass and diameter: $m= 10 ^{-10}$, $10^{-9}$, $10^{-8}$, and $10^{-7}$ kg. As illustrated by the blue line in Fig. 11a, the probability density of fine particles with $m= 10 ^{-10}$–10$^{-9}$ kg decreases with an increase in $u_*^\mathrm{t}$, whereas the density of coarse particles increases with $m= 10 ^{-9}$–10$^{-8}$ kg. Additionally, the probability density of the lower and higher $v_x$ increases with an increase in $u_*^\mathrm{t}$ (Fig. 11b), which is also confirmed in Fig. 10c. Therefore, the dependence of the particle mass and particle velocity on the friction velocity should be considered when analyzing and understanding particle shear stress.

Fig. 11

Differences in the probability-density distributions between two states with different friction velocities at the top boundary: a diameter d and b streamwise velocity component $v_x$ of incident particles. In a, the typical particle mass m corresponds to d and is indicated by vertical lines

5.4 Equilibrium Property of Friction Velocity at Different Heights

The simulation enabled the calculation of the fluid shear stress at the roughness length and the particle shear stress acting on the surface. However, accurately measuring these factors via both the wind-tunnel experiments as well as field observations of snow blowing is challenging. Therefore, the total shear stress on the surface during snow blowing was measured using several types of drag meters. Nemoto and Nishumura (2001) measured the total shear stress on the surface using the load cell in the wind-tunnel experiment; however, the particle shear stress was calculated from the loss of the horizontal momentum of snow saltation. Additionally, the fluid shear stress in the surface vicinity was estimated using small-size anemometers. Walter et al. (2012, 2014) investigated the variability of the fluid shear stress near the surface in wind-tunnel experiments using Irwin sensors, which were developed to measure pedestrian-level wind (Irwin 1981). According to Irwin (1981), wind speed below $z= 5$ mm cannot be accurately determined because of noise; thus, the minimum measurement height of the Irwin sensor is $z= 5$ mm. Considering this technical problem, we discuss the equilibrium property of the friction velocity in the range of height from the roughness length $z_0$ to $z= 5$ mm.

Figure 12 indicates the equilibrium property of the friction velocity at four different heights: $z = z_0 (= 10^{-5}$ m), 1, 2, and 5 mm. The results for $z = z_0$ are the same as those presented in Fig. 5b for $t > 16T$. Each point is the mean of 10 s calculated at each friction velocity at the top boundary $u_*^\mathrm{t}$ in the equilibrium state. When $u_*^\mathrm{t}$ exceeds the impact threshold $u_*^\mathrm{i}$, the particle transport is maintained by the wind. In such a case, the friction velocity at $z = z_0$ decreases with an increase in $u_*^\mathrm{t}$ and is smaller than the impact threshold, owing to the effect of particle shear stress, as shown in Fig. 8. Conversely, the friction velocity at $z = 1$ mm is equal to the impact threshold and independent of $u_*^\mathrm{t}$; however, the friction velocities at $z= 2$ and 5 mm increase more significantly with $u_*^\mathrm{t}$ than the impact threshold. Therefore, the friction-velocity trend shows a significant correlation with the measurement height. This result is consistent with the wind-tunnel experiment of Li and McKenna Neuman (2012). As the particle shear stress decreases with height, as shown in Fig. 7, the fluid shear stress (or friction velocity) changes significantly in the range of z = 10$^{-3}$–10$^{-2}$ m, as shown in Fig. 4. Thus, we determine that the equilibrium property of the friction velocity is affected by the height.

Fig. 12

Equilibrium property of the friction velocity $u_*$ at various heights z. The dashed and dotted lines indicate the fluid and impact thresholds, respectively

6 Summary and Conclusions

This study aimed to investigate the mechanism of fluid and particle shear stress acting on the surface during the aeolian transport of snow particles. For a numerical model, we used the random-flight model of drifting snow to simulate the vertical profile of the wind speed and two-dimensional trajectories of airborne particles (Nemoto and Nishimura 2004; Niiya and Nishimura 2017). The numerical simulations of the model were conducted in a quasi-two-dimensional system, wherein the stepwise change of the friction velocity at the top boundary was set to reproduce the transport hysteresis. In the simulation, particle transport was initiated when the friction velocity at the roughness length exceeded the fluid threshold, which was set as 0.24 m s$^{-1}$ by reference to the snow experimental data (Clifton et al. 2006). Sintering between the snow particles was not accounted for in the simulation, but the fluid threshold for the snow bed was generally uncertain in the snow field. The particle transport was maintained when the friction velocity was less than the fluid threshold. Moreover, transport hysteresis is confirmed in our simulation, and the impact threshold is numerically estimated from the change in the particle transport rate to 0.206 m s$^{-1}$. The friction velocity at the surface roughness length during particle transport is marginally less than the impact threshold; thus, our results of low friction velocities do not contradict Owen’s second hypothesis. The particle shear stress was calculated from the dynamics of the descending and ascending particles, to determine the change in the friction velocity (or fluid shear stress). We also calculated the particle shear stress acting on the surface by considering the loss of particle momentum during splash processes. In the equilibrium state of particle transport, a constant-stress layer forms, that is, the total shear stress summing fluid and particle shear stresses are vertically uniform. The increase in the particle shear stress at the surface leads to a decrease in the fluid shear stress (or friction velocity) at the roughness length. Our analysis of the particle dynamics in splash processes reveal that the number of incident and ejected particles is proportional to the difference between the fluid shear stress at the top boundary and the impact threshold converted into the shear stress. Meanwhile, the mean diameters and mean streamwise velocity components were approximately constant. Therefore, the number of airborne particles is the dominant factor determining the particle shear stress at the surface at a low friction velocity. Notably, the equilibrium property of the friction velocity in the surface vicinity changes significantly with height, first decreasing and then increasing with wind speed, because the particle shear stress decreases in the 1–10 mm height range.

The random-flight model used in our study reproduced the drifting snow at the uniform-horizontal airflow and is suitable for investigating the steady-state transport properties. Recently, numerical models considering the turbulent structure have been developed for the aeolian transport of snow and sand particles (Groot Zwaaftink et al. 2014; Okaze et al. 2018; Pähtz et al. 2020; Rana et al. 2020; Zheng et al. 2020; Rana et al. 2021). The spatio-temporal inhomogeneous structure and variability of aeolian transport have been investigated using these unsteady models. Although the random-flight model is relatively simple compared to these models, we believe that our results provide useful insights into near-surface transport properties of fluid and particle shear stresses. Nevertheless, various issues still need to be addressed. For example, this study could not investigate the surface shear stresses at the high friction velocity; hence, the coherent structures in turbulence such as the hairpin vortices were not considered. Thus, the transport hysteresis and surface shear stresses should be verified using an unsteady model. Moreover, the surface shear stress during particle transport should be measured using various methods in wind-tunnel experiments, because the fluid and particle shear stresses are related to aerodynamic entrainment and snow compaction, respectively, owing to splash processes. Furthermore, most models and experiments have been designed to study the dynamics of blowing snow on loose snow; however, the snow-surface conditions in nature, such as the hardness and the size distribution, vary with time. Therefore, the effect of the snow-surface conditions on the surface shear stress should be studied. Finally, future research on this topic will include measuring the surface shear stress on loose and hard snow in a wind-tunnel experiment of drifting snow.