Dynamic characteristics of snowfall particles in atmospheric turbulent boundary layer and its effect on dust wet deposition

Wet deposition by snowfall refers to the scavenging of atmospheric dust by snow particles. Existing models only consider vertical scavenging under still-air conditions, neglecting turbulence-induced complex vertical and horizontal motions of snowfall particles in the actual atmosphere boundary layer, which leads to inaccurate estimation of wet deposition flux. Currently, precise quantitative analysis of dust collection mechanism during snow particle setting remains lacking under turbulence. We employ the Euler-Lagrange numerical method, simplifying snow particles as spherical particles, to simulate and analyze snow particles dynamic characteristics and dust collection in turbulent boundary layer. The study introduces a dimensionless parameter ${\mathit{\alpha }}_{\mathrm{d}}=V_{\mathrm{t}}/\mathit{\kappa }{u}_{*}=$ 0.2 (V_t is the terminal settling velocity of snow particles, and κ= 0.4 is the von Kármán constant) to characterize the dynamic behavior of snow particles in turbulence. This parameter reflects the relative strength between gravitational setting and turbulent diffusion. Results show that when α_d> 0.2, the vertical relative motion dominates (with stronger dominances as α_d increases); when α_d< 0.2, horizontal relative motion becomes dominant. This shift in dynamic characteristics significantly enhances total dust collection capacity of snow particles and causes the dominant collection mechanism from vertical to horizontal: for α_d≥ 1, vertical collection for over 75 % of the total; when horizontal dominance, horizontal collection contributes over 50 %. The study demonstrates that neglecting horizontal collection underestimates wet deposition flux. Thus, we establish a quantitative wet deposition model, providing a theoretical basis for predicting atmospheric dust wet deposition and artificial dust removal.

1 Introduction

Atmospheric dust particles influence climate through mechanisms such as absorption or reflection of radiation and modification of cloud properties (Fuzzi et al., 2015; Rosenfeld et al., 2014), and they also lead to respiratory and cardiovascular diseases (Dockery and Stone, 2007; Tang et al., 2017). Variations in dust concentration in the downstream regions are primarily governed by the removal rate via atmospheric wet deposition (Seinfeld and Pandis, 2006), which occurs through in-cloud and below-cloud scavenging (also referred to as the “washout” process). The contribution of below-cloud scavenging to total wet deposition range from 50 % to 60 % (Ge et al., 2021) and is considered a primary sink for atmospheric aerosols (Textor et al., 2006). The efficiency of below-cloud scavenging is closely associated with the dynamic behaviours of precipitation particles (rain drops or snowfall particles): snowfall particles capture airborne dust through dynamic processes such as Brownian diffusion and inertial impaction. These scavenging mechanisms during precipitation events lead to high-flux dust deposition (Henzing et al., 2006).

The dynamic behaviour of precipitation particles in the turbulent boundary layer forms the basis for studying dust wet deposition, with their motion governed by the interplay between particle properties and surrounding turbulence. Dorgan and Loth (2004) demonstrated that in a horizontal turbulent boundary layer, as particle size increases, the dominant mechanism of particle motion shifts from turbulent diffusion to gravitational settling. This transition fundamentally arises from a reversal in the relative dominance of turbulent fluctuation forces and gravity acting on the particles. Studies have found that for small heavy particles fall through turbulence, their mean settling velocity can be significantly altered compared to quiescent fluid conditions (Wang and Maxey, 1993; Balachandar and Eaton, 2010). The fall can be hindered, e.g. if weakly inertial particles are trapped in vortices (Tooby et al., 1977), or if fast-falling particles are slowed down by nonlinear drag (Mei et al., 1991), or loiter in upward regions of the flow (Good et al., 2012). Field observations of snowfall events with significant variations in turbulence intensity by Li et al. (2021) revealed that turbulence is a critical factor influencing the falling velocity and spatial distribution of snowfall particles. Consequently, the regulatory effect of turbulence on the snowfall particle settling process significantly influences the efficiency of dust particles scavenging.

The scavenging coefficient Λ, representing the rate of change in aerosol mass concentration due to below-cloud scavenging (Seinfeld and Pandis, 2016):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}N\left(d_{\mathrm{p}}\right)}{\mathrm{d}t}}=-\mathrm{\Lambda }\left(d_{\mathrm{p}}\right)N\left(d_{\mathrm{p}}\right),\text{(2)}& {\displaystyle }\mathrm{\Lambda }\left(d_{\mathrm{p}}\right)={\int }_{\mathrm{0}}^{\mathrm{\infty }}A(V_{\mathrm{D}}-v_{\mathrm{d}})E(d_{\mathrm{p}},D_{\mathrm{p}})N\left(D_{\mathrm{p}}\right)\mathrm{d}{D}_{\mathrm{p}},\end{array}$$

where N(d_p) is the mass concentration of aerosols with particle diameter d_p, and N(D_p)dD_p is the number of precipitation particles with diameter between D_p and D_p+dD_p in a unit volume of air (m⁻³). V_D and v_d are the terminal velocities (m s⁻¹) of precipitation particles and dust particles, respectively, E(d_p, D_p) is the collection efficiency (dimensionless) between dust particles of size d_p and a precipitation particle of size D_p, and A is the effective cross-sectional area of precipitation particles projected normal to the fall direction (m²). The settling velocity and collection efficiency of precipitation particles, among other key parameters, are closely related to their dynamic characteristics during fall. Previous experimental and theoretical results have shown that E(d_p, D_p) is the results of net action of various forces influencing the relative motion of aerosols and precipitation particles. Langmuir (1948) first proposed the theoretical framework for “collision efficiency” through a dynamic model. Subsequently, Greenfield (1957) elucidated the effects of Brownian diffusion, interception, and inertial collision as three key mechanisms influencing E(d_p, D_p). Building on this, Slinn and Hales (1971) integrated theoretical and observational data, developing the widely applied Slinn formula. However, classical models assume a quiescent fluid condition, neglecting the effects of turbulence on particle transport. As research on multiphase flow dynamics has advanced, the regulatory role of turbulence in the motion of precipitation particles has been recognized to further influence the collection efficiency E(d_p, D_p) of dust aerosols. In recent years, studies on the mechanisms affecting E(d_p, D_p) have begun to incorporate the roles of physical processes such as atmospheric turbulence and the wake effects of raindrops. Hua et al. (2017) found that, in the context of atmospheric turbulence, the randomness of particle trajectories increases due to pulsating airflow, thereby enhancing the collection efficiency E(d_p, D_p) of stationary raindrops for aerosol particles, providing a new definition of aerosols capture by raindrop under the turbulent effect.

Compared to raindrops, research on the aerosol scavenging efficiency of snowfall particles remains relatively scarce. The wide variety of snow-particle shapes, sizes, and density, which results in significantly different terminal falling velocity, cross-sectional area, and surrounding flow field structure, leading to greater uncertainty in the scavenging coefficient (Zhang et al., 2013). Currently, relevant research can primarily be divided into two type. The first type uses field-collected snow sample data to establish relationships between precipitation, snowfall particle size distribution, and ground-level aerosol mass concentration, thereby evaluating the scavenging efficiency of snowfall on aerosol. The second type focuses on the collision scavenging mechanisms between a single snowflakes (or snow crystals) falling at terminal velocity V_D and aerosol particles, investigated through laboratory simulations or direct observations (Pruppacher et al., 1998). It should be noted that while these two approaches are widely represented in current literature, they do not exhaust all possible research directions. Studies have found that the terminal velocity of snowfall particles is a critical dynamic parameter determining scavenging efficiency. To quantitatively describe the evolution of this parameter, early research on snowfall particle terminal velocity established an empirical power-law formula based on the maximum dimension of snowfall particles through experimental measurements (Langleben, 1954). Subsequent studies progressively refined theoretical models (Mitchell and Heymsfield, 2005; Pruppacher and Klett, 1997), with the formula proposed by Mitchell and Heymsfield (2005) becoming the benchmark for current research due to its general applicability. The distinct physical reasons for the influence of V_D and the cross-sectional area (A) on the scavenging coefficient are straightforward: a higher fall velocity or a larger cross-sectional area of the collector (raindrop, snow, or ice) results in a faster collection process.

The existing scavenging coefficient expression is based on the assumption of a still-air condition and only considers the relative motion trajectories between snowfall particles and dust in the vertical direction. However, due to a 1–3 orders of magnitude difference in their particle sizes, their aerodynamic behaviours exhibit significant disparities: dust particles follow turbulent eddy motions, whereas larger snowfall particles, with greater inertia driving gravitational settling, exhibit limited responsiveness to turbulence. In the atmospheric turbulent boundary layer, this inertial disparity induces significant horizontal relative motion between particles, thereby influencing the scavenging efficiency of snowfall particles for dust. To address the limitations of traditional wet deposition models under quiescent condition in characterizing turbulence effects, this study employs numerical methods to dynamically reconstruct realistic turbulent fields, accurately simulating the motion trajectories of snowfall particles in turbulent environments. Based on Eulerian-Lagrangian coupled framework, we employ a gas-solid two-phase flow numerical simulations system that the OpenFOAM computational fluid dynamics platform with the Lagrangian particle tracking method. The system is used to simulate and analyse the motion behaviours of snowfall particles within the atmospheric turbulent boundary layer and their relationship with turbulence characteristics. Our work reveals the influence mechanisms of dust collection by snowfall particles and ultimately establishes a quantitative mathematical model for the snowfall scavenging process.

2 Numerical methods and validation

This study aims to quantitatively assess the influence of turbulence on the efficiency of snow particles in scavenging atmospheric dust. To achieve this, based on the assumptions and governing equations presented in Sect. 2.1, we first employ an OpenFOAM solver to simulate and generate a statistically stationary neutral atmospheric boundary layer turbulent wind field. Subsequently, a particle tracking module is developed based on the OpenFOAM Lagrangian particle tracking library. This module reads the pre-computed turbulent wind field as the background flow and solves the equation of motion for snow particles (Sect. 2.2) to track their trajectories and behaviour within the turbulent wind field. Finally, in Sect. 3, the impact of turbulence on the scavenging efficiency of snow particles is quantitatively analysed.

2.1 Flow field equations

In this study, atmospheric temperature variations are neglected, and the atmospheric boundary layer is assumed to be under neutral conditions. The three-dimensional incompressible Navier–Stokes (N–S) equations are employed to solve the turbulent flow field in the boundary layer, while the influence of energy exchange due to solar radiation on flow behaviour is disregarded. The mass conservation equation and momentum conservation equations implemented in OpenFOAM are expressed as:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial u_j}{\partial x_j}}=\mathrm{0},\text{(4)}& {\displaystyle \frac{\partial u_i}{\partial t}}+u_j{\displaystyle \frac{\partial u_i}{\partial x_j}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial p}{\partial x_i}}+\mathit{\nu }{\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_i}{\partial x_j\partial x_j}},\end{array}$$

where u_i represents the velocity component in the i-direction (with i = 1, 2, 3 corresponding to the x, y, z directions); x_i represents the spatial coordinate component (with i= 1, 2, 3 also corresponding to the x, y, z directions); and j is the summation index, indicating summation over all spatial directions (j= 1, 2, 3). p,ρ and ν represent pressure, air density, and kinematic viscosity, respectively.

The hybrid Reynolds-Averaged Navier-Stokes (RANS)/Large-Eddy Simulation (LES) strategy has been successfully applied to various complex flow scenarios, including the detailed resolution of high-Reynolds-number atmospheric boundary layer flows (Haupt et al., 2011), prediction of wind field characteristics around buildings (Liu and Niu, 2016), and investigation of near-field pollutant dispersion mechanisms (Lateb et al., 2014), with its reliability substantiated through extensive validation cases. This hybrid strategy employs RANS statistical modelling for small-scale eddies near the wall:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_j}{\partial x_j}}=\mathrm{0},\text{(6)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_i}{\partial t}}+{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_i{\stackrel{\mathrm{‾}}{u}}_j}{\partial x_j}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{p}}{\partial x_i}}+\mathit{\nu }{\displaystyle \frac{{\partial }^{\mathrm{2}}{\stackrel{\mathrm{‾}}{u}}_i}{\partial x_j\partial x_j}}-{\displaystyle \frac{\partial {\mathit{\tau }}_{ij}^{\mathrm{RANS}}}{\partial x_j}},\end{array}$$

here, the overbar “−” denotes the time-averaged quantity, and ${\mathit{\tau }}_{ij}^{\mathrm{RANS}}$ represents the Reynolds stress. In regions far from the wall, the model transitions to LES for resolving large-scale turbulence:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_j}{\partial x_j}}=\mathrm{0},\text{(8)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial x_i}}+{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{u}}_j}{\partial x_j}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{p}}{\partial x_i}}+\mathit{\nu }{\displaystyle \frac{{\partial }^{\mathrm{2}}{\stackrel{\mathrm{̃}}{u}}_i}{\partial x_j\partial x_j}}-{\displaystyle \frac{\partial {\mathit{\tau }}_{ij}^{\mathrm{LES}}}{\partial x_j}},\end{array}$$

where the superscript “ $\stackrel{\mathrm{̃}}{}$ ” denotes the resolvable-scale component of the physical quantity processed by the spatial filtering function. τ_ij represents the energy transfer between the filtered-out small-scale turbulence and the resolvable-scale turbulence, known as the subgrid-scale stress.

Due to the significant turbulence characteristics in the atmospheric boundary layer, accurate resolution of flow structures near the wall is required. We employ the Delayed Detached-Eddy Simulation (DDES) method based on the S-A model (Spalart et al., 2006). This method divides the simulation domain through the hybrid length scale l_DDES and establishes an adaptive turbulence resolution framework:

$$\begin{array}{}\text{(9)}& l_{\mathrm{DDES}}=d-f_{\mathrm{d}}max\left\{\mathrm{0},d-C_{\mathrm{DES}}\mathrm{\Delta }\right\},\end{array}$$

where d represents the Euclidean distance from a flow field grid point to the nearest solid wall, Δ= max (Δx, Δy, Δz) is the filter width defined as the local maximum grid spacing in the three direction, and the coefficient C_DES= 0.65. Spalart et al. (2006) introduced the shielding function f_d, improving the l_DDES length scale based on the DES method (Spalart, 2000). DDES addresses the critical issue of premature transition from RANS to LES in traditional hybrid methods.

2.2 Snowfall particle motion equations

We employ the Lagrangian particle tracking method to solve the motion of snowfall particles (Huang et al., 2024; Li et al., 2016, 2018). Snowfall particles are treated as a dilute phase with a volume fraction φ_V< 10⁻⁶, corresponding to a light snow event in terms of snow water equivalent rate. Under such light snow conditions, the disturbance of particle trajectories by turbulence becomes more significant, exerting a critical influence on the particle settlement distribution, considering only the one-way coupling of the fluid on the particles. As the particle diameter D_p is smaller than the Kolmogorov length scale η, particle collisions and rotational motion are neglected (Balachandar and Eaton, 2010). The shape of snowfall particles is simplified as spherical, with their motion governed by the following equation:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& m_{\mathrm{p}}{\displaystyle \frac{\mathrm{d}{\mathit{u}}_{\mathrm{p}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}\mathit{\pi }({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{a}})D_{\mathrm{p}}^{\mathrm{3}}g\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_{D_{\mathrm{p}}}{A}_{\mathrm{p}}\mathit{\rho }\left|{\mathit{u}}_{\mathrm{f}}-{\mathit{u}}_{\mathrm{p}}\right|({\mathit{u}}_{\mathrm{f}}-{\mathit{u}}_{\mathrm{p}}),\end{array}\end{array}$$

where u_p represents the velocity of snowfall particles; the snowfall particle density is ρ_p= 340 kg m⁻³, g= 9.81 m s⁻² is the gravitational acceleration, D_p is the snowfall particle diameter, ρ_a is the air density, A_p $=\mathit{\pi }{D}_{\mathrm{p}}^{\mathrm{2}}/\mathrm{4}$ is the projected cross-sectional area of the spherical particle, and u_f is the air velocity. The combined force of gravity and buoyancy acting on the snowfall particle is $F_{\mathrm{g}}+F_{\mathrm{b}}=\frac{\mathrm{1}}{\mathrm{6}}\mathit{\pi }({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{a}})D_{\mathrm{p}}^{\mathrm{3}}g$. $C_{D_{\mathrm{p}}}$ is the drag force coefficient for the snowfall particle; Re_p is the particle Reynolds number, with their respective expressions given by:

$$\begin{array}{}\text{(11)}& {\displaystyle }{C}_{D_{\mathrm{p}}}=\left\{\begin{array}{l}\frac{\mathrm{24}}{R{e}_{\mathrm{p}}}+\frac{\mathrm{6}}{\mathrm{1}+\sqrt{R{e}_{\mathrm{p}}}}+\mathrm{0.4},R{e}_{\mathrm{p}}\le \mathrm{1000}\\ \mathrm{0.424},R{e}_{\mathrm{p}}>\mathrm{1000}\end{array}\right.,\text{(12)}& {\displaystyle }R{e}_{\mathrm{p}}={\displaystyle \frac{\left|{\mathit{u}}_{\mathrm{f}}-{\mathit{u}}_{\mathrm{p}}\right|D_{\mathrm{p}}\mathit{\rho }}{\mathit{\mu }}}.\end{array}$$

The boundary condition for snowfall particles at the wall is set to trap, meaning that the particle trajectory calculation terminates when snowfall particles settle to the ground (wall). For all other boundaries, the condition is set to Escape: when a particle crosses such a boundary, its trajectory calculation is immediately terminated, and the particle is removed from the computational domain (neither reflected, trapped, nor rebounded), with no further tracking of its subsequent motion. It should be noted that the motion data generated by the particle prior to escape are still included in the statistical analysis.

2.3 Case setup

We conduct a systematic numerical simulation study on the motion of snowfall particles in the atmospheric boundary layer. The snowfall particle diameter was set to 10 discrete sizes ranging from 50 to 500 µm (specifically 50, 80, 100, 120, 150, 200, 300, 350, 400, and 500 µm), based on the typical sizes reported by Li et al. (2021). By specifying different initial inlet velocities (1, 5, 12, 15, and 20 m s⁻¹), we obtained the corresponding friction velocities: $u_{*}=$ 0.06, 0.31, 0.75, 0.93, and 1.18 m s⁻¹, with a focus on investigating the motion characteristics of snowfall particles under different u_* conditions. Considering that the scavenging coefficient is sensitive to particles dynamic processes within only a few hundred meters above the ground (Schumann, 1989), 2000 snowfall particles were released at 0.5 s intervals from a height of 100 m within a fully developed turbulent field. A time step of 0.05 s (CFL < 1) was set to ensure computational stability. Specific operating conditions and particle parameters are presented in Table 1.

Table 1 The particle Reynolds number (Re_p) of snowfall for various particle sizes under different friction velocity (u_*) conditions.

2.4 Turbulent inflow validation

To accurately simulate the turbulent wind field in the atmospheric boundary layer, it is necessary to establish inlet boundary conditions that conform to realistic turbulence characteristics. We construct a numerical wind tunnel based on the experimental model scale of Ishihara et al. (1999) (wind tunnel dimensions: 9 × 0.6 × 0.7 m³, with x, y and z representing the streamwise, vertical, and spanwise directions, respectively, as shown in Fig. 1), where the reference height above ground h= 0.04 m. A precursor simulation strategy based on the recycling method (Lund et al., 1998; Nozawa and Tamura, 2002; Vasaturo et al., 2018) is adopted, with a recycle station placed at the central streamwise cross-section of the computational domain (Fig. 1). Velocity data from this cross-section are collected in real-time and directly imposed at the inlet as the turbulent velocity boundary condition (inflow). The outlet is set as a free outflow (Outflow), the top and side boundaries are set as symmetry planes (symmetry), and the bottom is a no-slip wall (No-slip wall) (Zhou et al., 2022). A wall function (Wang and Moin, 2002) is employed to accurately capture the turbulent evolution near the rough wall. To maintain turbulence self-sustainability, a uniform grid resolution is used in the streamwise direction of the computational domain, with local refinement near the wall, resulting in a total of 2.65 million grid cells (as shown in Fig. 2).

Figure 1 Computational domain and boundary conditions.

Figure 2 Schematic of the overall grid division.

The wind profile generated by the recycling method is described by the logarithmic law:

$$\begin{array}{}\text{(13)}& \stackrel{\mathrm{‾}}{U}\left(z\right)=\left({\displaystyle \frac{u_{*}}{\mathit{\kappa }}}\right)\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right),\end{array}$$

where κ is the von Kármán constant with a value of 0.41, and z₀ is the aerodynamic roughness length. The root mean square of the streamwise fluctuating velocity is given by:

$$\begin{array}{}\text{(14)}& {\mathit{\sigma }}^{\mathrm{2}}={\sum }_{i=\mathrm{1}}^N(U_i-\stackrel{\mathrm{‾}}{U}{)}^{\mathrm{2}}/(N-\mathrm{1}),\end{array}$$

and turbulence intensity is defined as the ratio of the root mean square of velocity fluctuations (σ) to the mean wind speed:

$$\begin{array}{}\text{(15)}& I_u={\displaystyle \frac{{\mathit{\sigma }}_u}{\stackrel{\mathrm{‾}}{U}}}\times \mathrm{100}\mathit{\%},\end{array}$$

where I_u is the streamwise turbulence intensity (%) at a certain height, σ_u is the standard deviation of the streamwise velocity component (m s⁻¹), and $\stackrel{\mathrm{‾}}{U}$ is the time-averaged streamwise velocity (m s⁻¹).

The initial uniform inflow velocity is set to 5.2 m s⁻¹, with a time step of 0.001 s (satisfying CFL < 1). The turbulent field reaches a fully developed state after 10 s, and thus, the computational results from the 10–20 s time interval are selected for time-averaged statistical analysis. The fluctuating wind speed at the monitoring point exhibits typical random characteristics (Fig. 3). The simulated profiles of mean wind velocity and turbulence fluctuation based on these results are shown in Fig. 4. The simulation results of this study are compared with the numerical results of Zhou et al. (2022) and the experimental data of Ishihara et al. (1999).

Figure 3 Time history of fluctuating wind velocity at the monitoring point (x= 6 m, y= 0.6 m).

Figure 4 Characteristics of the turbulent boundary layer for the empty wind field obtained through numerical simulation, (a) normalized mean wind velocity, (b) normalized turbulent fluctuation, (c) power spectrum (streamwise).

The mean wind velocity is in good agreement with the simulated values of Zhou et al. (2022) and the experimental data, with a maximum relative error of approximately 4 %. The simulated standard deviation of the streamwise fluctuating velocity (${\mathit{\sigma }}_u/U_{\mathrm{ref}}$) is generally consistent with literature values and experimental data in the region $y/h>\mathrm{1.5}$, with localized deviations near the wall (the maximum relative error in turbulence fluctuations is approximately 15 %). The power spectrum shows good agreement with the simulation results of Zhou et al. (2022). Furthermore, the Reynolds number in this study (Re= 3.12 × 10⁶) is in exact agreement with the simulated value of Zhou et al. (2022), while the wind tunnel experimental value is 2.8 × 10⁶. Therefore, we adopt the recycling method based on the DDES turbulence model, which can accurately reproduce the characteristics of the atmospheric boundary layer turbulent wind field, validating the reliability and effectiveness of the method.

2.5 Grid independence verification

To reproduce the realistic atmospheric boundary layer, a numerical model of the empty wind field is established with dimensions: length (x) × height (y) × width (z) = 2000 × 400 × 600 m³. The boundary conditions are consistent with the previously validated case. To verify the reliability and accuracy of the numerical results, a grid independence study was conducted (Roache, 1993). The Grid Convergence Index (GCI) is defined as:

$$\begin{array}{}\text{(16)}& \mathrm{GCI}=F_{\mathrm{s}}{\displaystyle \frac{\left|\mathit{\epsilon }\right|}{r^m-\mathrm{1}}},\end{array}$$

where $\mathit{\epsilon }=(f_{\mathrm{1}}-f_{\mathrm{2}})/f_{\mathrm{1}}$ represents the relative error of grid convergence (f₁ and f₂ are the solutions from the fine grid and the previous coarse grid, respectively); the grid refinement ratio is $r_{k,k+\mathrm{1}}=(N_{k+\mathrm{1}}/N_k{)}^{\mathrm{1}/\mathrm{3}}$, where k is the grid level index used to distinguish computational grids with different densities, and N_k represents the total number of nodes at the kth grid level; the safety factor F_s= 1.25. We employ a second-order accurate scheme, thus m= 2.

We employ five grid schemes with different resolutions. Taking the calculated total relative travel distance (S_r) of snowfall particles with an initial velocity of 12 m s⁻¹ in air as an example, the GCI results are presented in Table 2. Using the criteria of ε≤ 1 % and GCI ≤ 5 % GCI, for grid numbers of 51.75, 68.48, and 79.17 million, the relative errors of S_r for snowfall particles with D_p= 500 µm are all less than 1 %, with corresponding GCI values of 4.09 % and 3.69 %, respectively. Considering both computational efficiency and resource costs, a grid number of 51.75 million is ultimately selected for the simulations.

Table 2 GCI Results for S_r of snowfall particles (D_p= 500 µm).

3 Numerical results

3.1 Statistical characteristics of the wind field

The wind field simulation results based on the aforementioned grid scheme (Fig. 5) demonstrate that the simulated atmospheric boundary layer turbulent wind field exhibits the irregularity and randomness characteristic of the real atmosphere. Figure 6 shows the corresponding mean wind velocity and turbulence intensity profiles. Comparisons with the international standard ASCE7-III (Zhou and Kareem, 2002) and the Chinese standard (Load Code for Design of Building Structures, GB 50009-2012) reveal that: the simulated mean wind velocity profile fall between the value of the Chinese standard and the international standard (Fig. 6a). As shown in Fig. 6b, the simulated turbulence intensity profile largely aligns with the lower bound of the range recommended by the ASCE 7-III standard (which permits variations within ±20 % of the baseline value, Kozmar, 2011), with a maximum relative error of 16 % in the near-ground region. Overall, the profile lies between the upper limits specified by the Chinese national standard and the ASCE 7-III standard. Thus, the simulation results effectively reproduce the wind field characteristics of a real atmospheric boundary layer.

Figure 5 Instantaneous streamwise velocity distribution in the x–y plane of the turbulent wind field at t= 3500 s.

Figure 6 Comparison of (a) the profile of the ratio of measured mean wind velocity to the mean wind velocity at the reference height of 10 m, and (b) the turbulence intensity profile with the Chinese standard and the international standard.

3.2 Analysis of snowfall particle dynamic characteristics

Analysis of the motion trajectory data over a 1000 s duration indicates that larger-diameter (D_p≥ 300 µm) snowfall particles, due to their greater mass and inertia, primarily exhibit stable vertical settling with a weak response to turbulent disturbances. In contrast, the motion of smaller-diameter snowfall particles is governed by both gravity and turbulent diffusion: gravity dominates vertical motion, while horizontal transport induced by turbulent eddies significantly broadens their settling range and increases uncertainty in settling positions (Fig. 7a). As the friction velocity increases, the motion trajectories of snowfall particles (particularly those with smaller diameters) are more strongly affected by turbulent disturbances, displaying more pronounced perturbations and distortions (Fig. 7b). We characterize the dynamic state of snowfall particles in the airflow by referencing the ratio α_d ($=V_{\mathrm{t}}/w$; dimensionless), where V_t is the terminal settling velocity of the snowfall particles and w ($w=\mathit{\kappa }{u}_{*}$, with κ= 0.4 as the von Kármán constant) is the vertical diffusion velocity of the fluid (Huang and Shi, 2017; Scott, 1995). This parameter reflects the competition between gravitational settling and airflow disturbances. The terminal settling velocity V_t can be calculated using the following equation (Carrier, 1953), which has been validated for snow particles within the typical size range (Huang and Shi, 2017):

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& V_{\mathrm{t}}=-{\displaystyle \frac{A}{D_{\mathrm{p}}}}+\sqrt{\left({\displaystyle \frac{A}{D_{\mathrm{p}}}}\right)+B{D}_{\mathrm{p}}},A=\mathrm{6.203}{\mathit{\nu }}_{\mathrm{a}},\\ & B={\displaystyle \frac{\mathrm{5.516}{\mathit{\rho }}_{\mathrm{p}}}{\mathrm{8}{\mathit{\rho }}_{\mathrm{a}}}}g,\end{array}\end{array}$$

where V_t is the terminal setting velocity of snow particles, D_p is the diameter of snow particle (m), ν_a is air viscosity coefficient (m² s⁻¹), ρ_p and ρ_a are the density of snow particles and air, respectively (kg m⁻³), and g is the acceleration due to gravity (m s⁻²).

Figure 7 Random motion trajectories of (a) snowfall particles with different diameters at $u_{*}=$ 0.75 m s⁻¹ and (b) snowfall particles (D_p= 100 µm) under different u_* conditions.

Figure 8 Temporal variations of snowfall particle velocity (the green solid line), wind field velocity at the snowfall particle location (the black solid line), and modulus of the difference between the wind speed at the particle location and the snow particle velocity (the blue solid line) under $u_{*}=$ 0.75 m s⁻¹.

Figure 8 shows the temporal evolution curves of u_f@p (where u_f@p is the fluid velocity at the particle's location), $|u_{\mathrm{f}\mathrm{@}\mathrm{p}}-u_{\mathrm{p}}|$ (modulus of the difference between the wind speed at the particle location and the snow particle velocity), and u_p for snowfall particles of two different diameters during their falling process. The speed of the snowfall particle relative to the wind field at its location is defined as the relative velocity of the snowfall particle (V_r,  m s⁻¹):

$$\begin{array}{}\text{(18)}& V_{\mathrm{r}}=\left|{\mathit{u}}_{\mathrm{f}\mathrm{@}\mathrm{p}}-{\mathit{u}}_{\mathrm{p}}\right|.\end{array}$$

The results shown in Fig. 8 indicate that the local flow fields experienced by snowfall particles of different diameters exhibit significant differences. Snowfall particles with D_p= 300 µm demonstrate a higher relative speed (V_r), suggesting a stronger ability to resist flow disturbances, more independent motion, and less susceptibility to turbulent influences. In contrast, smaller-diameter snowfall particles (D_p= 50 µm) exhibit significantly lower V_r, making them more susceptible to turbulence and displaying greater flow-following behaviour. This difference in V_r reflects the distinct dynamic effects of the flow field on snowfall particles of varying diameters. Further analysis reveals that the distribution of the dimensionless mean relative velocity ($V_{\mathrm{r}}/V_{\mathrm{t}}$; dimensionless) for snowfall particles follows a normal distribution, with the distribution shape remaining unaffected by changes in u_* and D_p (Fig. 9). The parameters u_* and D_p primarily influence the characteristic quantities (mean and variance) of the normal distribution.

Figure 9 Probability distribution of $V_{\mathrm{r}}/V_{\mathrm{t}}$ for (a) snowfall particles with Different diameter at $u_{*}=$ 0.75 m s⁻¹ and (b) snowfall particles with D_p= 100 µm under different friction velocities.

To quantitatively characterize the influence of different friction velocity (u_*) on the relative motion velocity of snowfall particles with varying diameter, Figs. 10 and 11 show the mean and standard deviation of the relative velocity as functions of u_* for different particle diameters. The model employs an exponential function of the form $e^{(-k{D}_{\mathrm{p}})}$, which decays with particle diameter, to capture the nonlinear modulation of particles inertia on turbulence response as particles size increases, thereby accurately describing the particles size distribution characteristics transitioning from “turbulence-following” to “inertia-dominated” behaviour.

As shown in Fig. 10a, the dimensionless mean relative velocity ($V_{\mathrm{r}}/V_{\mathrm{t}}$) increases exponentially with u_*, and the growth rate decreases significantly with increasing particles size. This trend indicates that small particles are highly sensitive to changes in turbulence, while large particles, dominated by inertial effects, exhibit a significantly weakened response to variations in u_*. To eliminate the influence of gravitational settling, the model first introduces a dimensionless relative velocity ($V_{\mathrm{r}}/V_{\mathrm{t}}$), placing particles of different sizes on a unified benchmark for dynamic comparison. On this basis, an exponential function $e^{\left(b{u}_{*}\right)}$ dependent on u_* is adopted to describe the variation in the mean relative velocity. The variation of the mean relative velocity of snowfall particles with friction velocity and particle diameter can be expressed by Eq. (19):

$$\begin{array}{}\text{(19)}& V_{\mathrm{r}}=V_{\mathrm{t}}\cdot a{e}^{b\cdot u_{*}},\end{array}$$

where parameters a and b are function of the snow particle diameter (D_p, mm), with their specific forms obtained through fitting:

$$\begin{array}{}\text{(20)}& a=(\mathrm{0.97}+\mathrm{19.85}{e}^{-\mathrm{57.6}{D}_{\mathrm{p}}}),b=\mathrm{2.27}{e}^{-\mathrm{11.02}{D}_{\mathrm{p}}},\end{array}$$

as illustrated in Fig. 10b, both parameters a and b decrease in a negative exponential manner with increasing snow particle diameter, further confirming the modulating effect of particle inertia on the turbulent response.

Figure 10 (a) Dimensionless mean relative velocity of snowfall particles with different D_p under various u_* conditions, and (b) the corresponding fitting parameters as a function of D_p (symbols represent simulated values, and solid lines represent fitted values).

Figure 11 (a) Standard deviation of the fluctuating relative velocity of snowfall particles with different D_p as a function of u_*, and (b) the corresponding fitting parameters as a function of D_p.

As shown in Fig. 11a, the standard deviation of the relative velocity fluctuations ($V_{\mathrm{r}}^{\prime }$) increases linearly with increasing friction velocity. In turbulent boundary layers, the root mean square of velocity fluctuations is typically proportional to the friction velocity, reflecting the fundamental principle that turbulent fluctuation energy originates from surface shear stress. Therefore, a linear function of the form $V_{\mathrm{r}}^{\prime }$ $=c_{\mathrm{0}}+k\left(D_{\mathrm{p}}\right)\cdot u_{*}$ is adopted to intuitively capture the response of particle velocity fluctuations to changes in turbulence intensity. The variation of $V_{\mathrm{r}}^{\prime }$ with friction velocity and snowfall particle diameter can be expressed by Eq. (21):

$$\begin{array}{}\text{(21)}& V_{\mathrm{r}}{}^{\prime }=\mathrm{0.003}+b_{\mathrm{1}}\cdot u_{*},\end{array}$$

where the intercept of the linear function is independent of snowfall particle diameter, and the parameter b₁ exhibits negative exponential growth with increasing particle diameter (D_p, mm).

$$\begin{array}{}\text{(22)}& b_{\mathrm{1}}=\mathrm{0.17}-\mathrm{0.05}{e}^{-\mathrm{3.25}{D}_{\mathrm{p}}}.\end{array}$$

The coefficient of determination (R²) for the mean relative velocity exceeds 0.90, while that for the fluctuating standard deviation of relative velocity surpasses 0.99. Consequently, Eq. (23) should be used for predicting snowfall particle relative velocity.

$$\begin{array}{}\text{(23)}& V_{\mathrm{r}}=\stackrel{\mathrm{‾}}{V_{\mathrm{r}}}+V_{\mathrm{r}}{}^{\prime }.\end{array}$$

Limited by the difficulty of fully simulating the trajectories of small-diameter particles, we introduce the product of the snowfall particle relative velocity V_r (for which a quantitative characterization method has been established) and the particle suspension time to quantify the relative travel distance of particles during their airborne phase. The mean particle suspension time is calculated based on the deposition velocity (V_d, m s⁻¹) theoretical framework of the two-layer model proposed by Zhang and Shao (2014):

$$\begin{array}{}\text{(24)}& V_{\mathrm{d}}={\left(r_{\mathrm{g}}+{\displaystyle \frac{r_{\mathrm{s}}-r_{\mathrm{g}}}{\exp (r_{\mathrm{a}}/r_{\mathrm{g}})}}\right)}^{-\mathrm{1}},\end{array}$$

where the gravitational resistance is denoted as $r_{\mathrm{g}}=\mathrm{1}/V_{\mathrm{t}}$, the aerodynamic resistance r_a:

$$\begin{array}{}\text{(25)}& r_{\mathrm{a}}=\left({\displaystyle \frac{B_{\mathrm{1}}\cdot S{c}_T}{\mathit{\kappa }{u}_{*}}}\right)\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right),\end{array}$$

and the surface collection resistance r_s:

$$\begin{array}{}\text{(26)}& r_{\mathrm{s}}={\left(w_{\mathrm{dm}}\left(S{c}^{-\mathrm{1}}+{\mathrm{10}}^{-\mathrm{3}/T_{\mathrm{p}}+}\right)+V_{\mathrm{t},\mathit{\delta }}\right)}^{-\mathrm{1}},\end{array}$$

where B₁ is an empirical constant, Sc and Sc_T are the Schmidt number and turbulent Schmidt number, respectively, w_dm is the conductance for momentum, $T_{\mathrm{p}}^{+}$ Is the dimensionless particle relaxation time, and V_t,δ is the terminal velocity of particles at the top of the laminar layer, typically approximated as V_t. Figure 12a shows a comparison between the theoretical and simulated values of snowfall particle deposition velocity V_d. At $u_{*}=$ 0.75 m s⁻¹, the simulated V_d (statistical mean of vertical velocity) is in close agreement with the theoretical values from Zhang and Shao (2014) (Fig. 12a), validating the reliability of the Lagrangian particle tracking method for snowfall particle transport simulations. Figure 12b shows that the theoretical V_d of snowfall particles increases with increasing particle diameter, but the dimensionless settling velocity $V_{\mathrm{d}}/u_{*}$ decreases as u_* increases, indicating that enhanced turbulence suppresses the settling efficiency of snowfall particles. The particle settling time is further calculated as:

$$\begin{array}{}\text{(27)}& T_{\mathrm{d}}={\displaystyle \frac{H_{\mathrm{0}}}{V_{\mathrm{d}}}},\end{array}$$

where H₀ (m) is the release height of the snowfall particles. As shown in Fig. 13, the dimensionless setting time ($T_{\mathrm{d}}\cdot u_{*}/H$), defined to eliminate scale effects (where H is the boundary layer height), for snowfall particles in this diameter range exhibits negative exponential decay with increasing D_p. This indicates that larger-diameter snowfall particles, characterized by a higher deposition velocity (V_d), undergo rapid vertical settling. In contrast, smaller particles, with a lower V_d and greater susceptibility to turbulent fluctuations, settle more slowly and display spatiotemporal instability. Meanwhile, an increase in u_* enhances the particle retention effect, making particles more likely to remain trapped for extended periods in low-kinetic-energy regions or eddy structures, significantly impacting the settling and transport processes.

Figure 12 (a) Comparison of simulated and theoretical snowfall particle deposition velocities at $u_{*}=$ 0.75 m s⁻¹, and (b) $V_{\mathrm{d}}/u_{*}$ as a function of particle diameter (Zhang and Shao, 2014).

Figure 13 Dimensionless settling time of snowfall particles as a function of D_p for different u_*.

To quantify the average motion distance of snowfall particles relative to the wind field during their settling time, namely the relative motion distance S_r:

$$\begin{array}{}\text{(28)}& \begin{array}{rl}{S}_{\mathrm{r}}& =V_{\mathrm{r}}\cdot T_{\mathrm{d}}=V_{\mathrm{t}}(\mathrm{0.97}+\mathrm{19.85}{e}^{-\mathrm{57.6}{D}_{\mathrm{p}}})\\ & e^{\mathrm{2.27}{e}^{-\mathrm{11.02}{D}_{\mathrm{p}}}\cdot u_{*}}\left({\displaystyle \frac{H_{\mathrm{0}}}{V_{\mathrm{d}}}}\right),\end{array}\end{array}$$

where V_t and V_d are the terminal settling velocity and settling velocity of snow particles (m s⁻¹), D_p is the snow particle diameter (mm), u_* is the friction velocity (m s⁻¹), and H₀ is the particle release height (m). The calculations are based on the release height H₀= 100 m and the definition of relative travel distance presented in this study (Eq. 28). The results indicate that the vertical relative travel distance (S_r1) of snowfall particles decreases monotonically with decreasing particle diameter. Further analysis incorporating the influence of friction velocity reveals that for D_p>200 µm, S_r1 increases with increasing friction velocity, suggesting that larger-diameter particles exhibit a greater vertical motion range in stronger airflow. In contrast, for D_p≤ 200 µm, S_r1 decreases with increasing friction velocity, indicating that smaller-diameter snowfall particles are easily entrained by turbulent eddies, displaying strong flow-following behaviour. Their lower relative velocity further affects S_r1, revealing the unique motion characteristics of smaller-diameter snowfall particles in the turbulent boundary layer. The horizontal relative travel distance (S_r2) decreases with increasing D_p (Fig. 14b): smaller particles are more readily transported over long distances by the airflow, resulting in larger S_r2, which increases significantly with u_*.

Figure 14 Snow particle (a) S_r1, and (b) S_r2 as a function of D_p for different u_*.

Further analysis of Fig. 15a shows that the ratio $S_{\mathrm{r}\mathrm{1}}/S_{\mathrm{r}\mathrm{2}}$ increases monotonically with the parameter α_d. The critical state ($S_{\mathrm{r}\mathrm{1}}/S_{\mathrm{r}\mathrm{2}}\approx$ 1) corresponds to an α_d value of 0.2, indicating that under turbulent conditions, snowfall particles with α_d = 0.2 achieve a maximum stable motion state under the combined effects of turbulence and gravitational settling – neither settling too rapidly due to excessive gravity nor remaining suspended too long due to strong turbulent entrainment. If α_d>0.2, vertical relative motion dominates, with the dominance becoming more pronounced as the parameter increases; if α_d< 0.2, horizontal relative motion predominates. Figure 15b demonstrates that the total relative travel distance (S_r) of snowfall particles decreases in a negative exponential manner with increasing D_p: the motion of larger-diameter particles is dominated by gravity, with S_r approaching H₀ and showing little sensitivity to changes in u_*; for smaller-diameter particles, S_r increases significantly with u_* due to enhanced turbulence, which extends their retention time and motion path in the wind field.

Figure 15 (a) Ratio of S_r1 to S_r2 for snowfall particle as a function of the particle parameter α_d, and (b) S_r of snowfall particle as a function of D_p for different u_*.

3.3 Analysis of the impact of snowfall particle dynamic characteristics on dust wet deposition

Based on Slinn's theory (Shao, 2008; Slinn, 1984), raindrops and dust particles exhibit relative motion in the atmosphere due to differences in their terminal velocities. Suppose a raindrop falls with a terminal velocity V_t and a dust particle falls with a terminal velocity v_t. Then, the relative speed at which the raindrop approaches the dust particle is (V_t−v_t). If the number density (number of particles per unit volume) of dust particle is n, then during the time interval Δt, the total number of dust particles that can be captured by the raindrop is $n(V_{\mathrm{t}}-v_{\mathrm{t}})\left(\mathrm{\Delta }t\right)\mathit{\pi }(R+r{)}^{\mathrm{2}}$. The number of particles actually captured by the raindrop is:

$$\begin{array}{}\text{(29)}& q_{\mathrm{0}}=e_{\mathrm{s}}n\mathit{\pi }(R+r{)}^{\mathrm{2}}(V_{\mathrm{t}}-v_{\mathrm{t}})\mathrm{\Delta }t,\end{array}$$

where the collection efficiency e_s is a function of both dust radius (r) and raindrop radius (R).

Existing studies commonly adopt the equivalent collision efficiency assumption (Pruppacher et al., 1998), implying that when a snowfall particle comes into contact with dust, the collision directly results in the capture of the dust. According to the wet deposition formula (Eq. 29), the total number of particles collected by a snowfall particle as it settles from the air to the ground, the collection amount (q) can be expressed as:

$$\begin{array}{}\text{(30)}& q=e_{\mathrm{s}}n\mathit{\pi }(R+r{)}^{\mathrm{2}}\sum _{i=\mathrm{0}}^{T_{\mathrm{d}}}(V_{\mathrm{t}}-v_{\mathrm{t}})\mathrm{\Delta }{t}_i,\end{array}$$

where ${\sum }_{i=\mathrm{0}}^{T_{\mathrm{d}}}(V_{\mathrm{D}}-v_{\mathrm{d}})\mathrm{\Delta }{t}_i$ represents the relative motion distance of snowfall particle relative to dust particles.

Figure 16 Schematic diagram of the dust collection process by a falling snowfall particle.

We assume that Aitken mode dusts (d_p= 20–100 nm) with a concentration of n= 1.5 × 10⁴ cm⁻³ is uniformly distributed within the boundary layer and fully follows the flow field motion, such that the velocity of dust particles equals the wind field velocity at their location, i.e., v_t=u_f@p. This concentration value is based on observational data from the Beijing region in China, where the number concentration of Aitken mode dusts is typically on the order of 10⁴ cm⁻³ under heavily polluted conditions (Liu et al., 2016). Therefore, the S_r is also defined as the distance of relative motion between snowfall particles and dust particles. The collection amount is the number of dust particles collected within the volume swept by the snowfall particles during their falling process, as illustrated schematically in Fig. 16a. Consequently, Eq. (30) can also be expressed as:

$$\begin{array}{}\text{(31)}& Q=e_{\mathrm{s}}n\mathit{\pi }(R+r{)}^{\mathrm{2}}{S}_{\mathrm{r}}=e_{\mathrm{s}}n\mathit{\pi }(R+r{)}^{\mathrm{2}}{V}_{\mathrm{r}}\cdot T_{\mathrm{d}},\end{array}$$

to investigate the influence mechanism of turbulence-induced relative motion between snow particles and dust on dust scavenging, we assume that a snow particle can completely collect all dust particles along its trajectory, i.e., the collection efficiency e_s= 1. Under this idealized assumption, the total number of dust particles contained within the spatial volume swept by a snow particle during its settling process is defined as the maximum potential collection amount (Q) by the snow particle, as illustrated in Fig. 16a. It should be noted that in the actual atmosphere, the actual collection efficiency of snow particles for dust is generally less than 1.0 due to aerodynamic effects (such as Brownian diffusion, interception, and inertial impaction). Therefore, the results calculated in this study represent the theoretical upper limit of the dust removal capacity. Correspondingly, its vertical and horizontal components represent the maximum potential vertical collection amount (Q₁) and the maximum potential horizontal collection amount (Q₂), respectively, which are used to evaluate the influence of turbulence on the dust removal capability in different directions.

Figure 17 Maximum potential collection amount of dust by snowfall particles in various directions as a function of D_p for different u_* (the vertical dashed line in the figure represents the boundary between vertically dominated and horizontally dominated snow particle collection mechanisms, and the corresponding particle size is defined as the critical diameter).

The maximum potential collection amount of dust by a single snowfall particle during its falling is calculated using the wet deposition formula (Eq. 31). Figure 17 reveals the influence of snowfall particle diameter on the maximum potential total collection amount (Q) and its directional components under different friction velocities: when $u_{*}\le$ 0.06 m s⁻¹, the motion of all D_p snowfall particles is dominated by vertical settling, resulting in Q being almost entirely contributed by the maximum potential vertical collection amount (Q₁). As u_* increases to 1.18 m s⁻¹, large-diameter snowfall particles still exhibit vertical collection dominance, with Q₁ approaching Q for these particles; whereas small-diameter snowfall particles (D_p <  120 µm) shift to being dominated by the maximum potential horizontal collection dominance (Q₂). Moreover, the critical diameter (vertical dashed line) increases from 50 to 120 µm with increasing friction velocities. These results demonstrate that increasing u_* drives the snowfall particle collection mechanism from vertical dominance to horizontal dominance: large-diameter snowfall particles retain vertical dominance due to gravitational settling, while small-diameter snowfall particles exhibit enhanced horizontal collection capability owing to intensified turbulent diffusion and prolonged settling time. As shown in Fig. 18, snowfall particles with α_d≥ 1, Q₁ accounts for over 75 % of the maximum potential total collection (with the corresponding Q₂ proportion being less than 25 %). In contrast, for snowfall particles with α_d< 0.2, the horizontal collection capability significantly strengthens with increasing u_*; the minimum proportion of Q₂ in the maximum potential total collection is approximately 50 % (at which point the maximum proportion of Q₁ is also around 50 %. Therefore, enhanced u_* markedly boosts the horizontal collection capability of small-diameter snowfall particles.

Figure 18 Proportions of the maximum potential vertical collection amount (Q₁) and maximum potential horizontal collection amount (Q₂) in the maximum potential total collection amount (Q) for snowfall particles as functions of α_d under different friction velocities.

Compared to gravitational settling, the instantaneous velocity fluctuations in turbulent wind fields significantly increase the complexity of interactions between snowfalls and the flow field. Turbulence causes snowfall particles to experience regions of varying velocities, resulting in acceleration, deceleration, or entrainment into vortices, thereby affecting their motion trajectories and dust collection capability. We define the growth rate of the maximum potential collection amount of snowfall particles as the collection enhancement efficiency relative to the gravitational settling, to quantify the impact of turbulence on the snowfall particles collection capability. The growth rate of the maximum potential collection amount of snowfall particles (r):

$$\begin{array}{}\text{(32)}& r={\displaystyle \frac{Q-Q_{\mathrm{0}}}{Q_{\mathrm{0}}}},\end{array}$$

where Q and Q₀ represent the maximum potential collection amount of dust by snow particles under friction velocity and gravitational settling conditions, respectively. As shown in Fig. 19, the results indicate that the growth rate of the maximum potential collection amount of snowfall particles exhibits a decreasing trend with increasing snowfall particle parameter α_d. When α_d≥ 4.5, r≈ 0; whereas when α_d≤ 0.5, the r can reach approximately 50 % (Fig. 19), indicating that the collection behaviour of large-diameter snow grains is primarily dominated by gravity, with limited influence from turbulence; under higher u_*, snowfall particles with α_d≤ 0.5 experience prolonged settling time due to turbulent effects, resulting in a significant enhancement of collection capability. Meanwhile, the growth rate of the maximum potential collection amount of snowfall particles increases with u_*. Further fitting analysis shows that r decreases as a negative exponential function with increasing dimensionless parameter α_d ($=V_{\mathrm{t}}/\mathit{\kappa }{u}_{*}$; dimensionless), which can be expressed as:

$$\begin{array}{}\text{(33)}& r=\mathrm{0.26}\cdot {\left({\mathit{\alpha }}_{\mathrm{d}}\right)}^{-\mathrm{1.20}}.\end{array}$$

The coefficient of determination (R²) for the growth rate of the maximum potential collection amount of snow particles exceeds 0.99, thereby establishing a mathematical model for the growth rate of the maximum potential collection amount of dust by snow particles. Given u_* and D_p, it can effectively predict the enhancement in snowfall collection capacity attributable to turbulent effects relative to gravitational settling.

Figure 19 Growth rate of the maximum potential collection amount of snowfall particles as a function of the parameter α_d (symbols represent simulated values, and solid lines represent fitted values).

Figure 20 Number concentration of snowfall particle as a function of D_p for different precipitation intensity.

Figure 21 Maximum potential total collection amount of the snow particle population (Q_N1) as a function of (a) D_p for different u_* at R= 8.5 mm h⁻¹, and (b) R for snowfall particles at $u_{*}=$ 0.75 m s⁻¹, respectively.

Snowfall particles settle in the form of a specific particles number size distribution during snowfall, which directly influences their scavenging effect on atmospheric dust. The snowfall particles size spectra is affected by ambient temperature, particle habit, precipitation intensity, and the stage of cloud and precipitation development (Harimaya et al., 2004; Woods et al., 2008). In practical applications, empirical formulas derived from raindrop size spectra are commonly used to approximate the size distribution of natural snow (Woods et al. 2008), among which the exponential distribution is widely applied in cloud microphysics to describe the snowfall particles size spectrum (Feng, 2009; Solomon et al., 2009), with its basic form expressed as:

$$\begin{array}{}\text{(34)}& N\left(D_{\mathrm{p}}\right)=N_{\mathrm{0}e}\exp (-{\mathit{\beta }}_e{D}_{\mathrm{p}}),\end{array}$$

where N_0e is the intercept parameter, representing the theoretical concentration when snowfall particle diameter approaches zero; β_e is the slope parameter, controlling the decay rate of the size distribution. Referring to the classification of daily precipitation intensity during the snow season by Suriano et al., 2023 (“Light”: ≤ 50.8 mm, “Moderate”: 50.8 $<R\le$ 152.4 mm, “Heavy”: 152.4 $<R\le$ 254.0 mm, “Extreme”: > 254.0 mm), we set different precipitation intensities R (2, 6, 8.5, 15, and 20 mm h⁻¹) to systematically analyse the wet deposition efficiency of snowfall particle populations on dust. As shown in Fig. 20, for the same particle size, the snowfall particle number concentration increases with increasing R, with the total number concentration spanning two orders of magnitude (Zhang et al., 2013), and exhibits exponential decay with increasing D_p. The variation of snowfall particle number concentration with D_p (µm) and R (mm h⁻¹) can be expressed as:

$$\begin{array}{}\text{(35)}& {\displaystyle }{N}_{\mathrm{1}}=(-\mathrm{0.8}+\mathrm{42}R)\times {\mathrm{0.986}}^{D_{\mathrm{p}}}\times {\mathrm{10}}^{\mathrm{6}},\text{(36)}& {\displaystyle }{Q}_{N_{\mathrm{1}}\left(D_{\mathrm{p}}\right)}=e_{\mathrm{s}}n\mathit{\pi }(R+r{)}^{\mathrm{2}}(V_{\mathrm{r}}\cdot T_{\mathrm{d}})N_{\mathrm{1}}.\end{array}$$

Based on the previously established results for the maximum potential collection amount of individual snow particles, the maximum potential total collection amount of a snow particle population can be calculated using Eq. (36). As shown in Fig. 21a, for snowfall particles with diameters less than 200 µm, the maximum potential total collection amount is highly sensitive to u_* and primarily dominated by turbulence: under high u_* conditions, turbulence prolongs the residence time of small-diameter snowfall particle populations in the atmosphere, thereby enhancing dust collision efficiency and significantly increasing their collection amount. In contrast, for snowfall particle populations larger than 200 µm, the maximum potential total collection amount is mainly dominated by inertia, with weaker influence from turbulence; the collection amount curves under different u_* tend to converge, indicating reduced dependence of large-diameter snowfall particle on turbulence. Within the $u_{*}=$ 0.06–1.18 m s⁻¹ range, the maximum potential total collection amount of snowfall particles in the 100–150 µm diameter interval reaches its peak, suggesting that this size range achieves optimal scavenging efficiency under the given precipitation intensity. Figure 21b shows that Q is linearly positively correlated with R, indicating that increasing precipitation intensity can significantly enhance the dust scavenging capability of snowfall particle populations. Among snowfalls of different D_p, the populations in the 100–150 µm range exhibits the largest slope, demonstrating that its scavenging efficiency is most sensitive to changes in precipitation intensity. It should be noted that the conclusions in Sect. 3.4 above are based on the assumption that the dust particle size falls within the Aitken nucleus mode of 20–100 nm. For the commonly observed atmospheric particle size range of 0.1–10 µm, the dust concentration corresponding to that specific particle size can be substituted into the wet deposition model established in this paper (Eq. 36) to derive the optimal scavenged snow particle size under such conditions.

4 Conclusion

In this study, we investigated the motion of snowfall particles in a turbulent boundary layer using the Eulerian-Lagrangian method, under the assumption that particles are spherical. The mechanism for dust collection by snowfall is investigated by analysing the snowfall particles behaviours and its dependence on turbulent characteristics. The main conclusions are as follows:

This study reveals the combined influence of particle diameter and friction velocity (u_*) on snow particle dynamics characteristics. The findings demonstrate that turbulence significantly prolongs the residence time of small snowfall particles in the atmosphere, thereby enhancing the randomness of their trajectories and consequently increasing the uncertainty in their final deposition. Furthermore, turbulence affects the particles' relative motion distance to the air: large particles, governed primarily by gravity, maintain stable vertical trajectories with the total relative motion distance close to the release height. In contrast, small particles experience a significantly expanded horizontal range due to turbulence, leading to a substantial increase in their relative motion distance. Overall, the total relative motion distance of snowfall particles exhibits a negative exponential decay with increasing particle diameter. As friction velocity increases, the dominance of vertical and horizontal motion of snowfall particles shifts: The ratio of the vertical relative motion distance (S_r1) to the horizontal relative motion distance (S_r2) of snowfall particles increases monotonically with the parameter α_d, where the critical state ($S_{\mathrm{r}\mathrm{1}}/S_{\mathrm{r}\mathrm{2}}\approx$ 1) corresponds to a snow particle α_d value of 0.2, indicating that snow particles with α_d= 0.2 can maintain long-term motion stability under turbulent conditions. If this parameter exceeds 0.2, vertical relative motion dominates, with the dominance becoming more pronounced as the parameter increases; if α_d≤ 0.2, the horizontal relative motion exerts the primary influence. The influence mechanism of snowfall particle dynamic characteristics in turbulent environments on dust wet deposition indicates that the scavenging efficiency of snowfall particles for dust is significantly different from gravitational settling: when the snow particle parameter α_d≥ 4.5, the collection growth rate (r) approaches zero, and the collection behaviour is entirely dominated by gravity; whereas when α_d≤ 0.5, the r can reach approximately 50%, with turbulence significantly enhancing the collection capacity of small snow particles. With increasing friction velocity, the dust collection capability of snowfall particles is markedly enhanced, and the collection mechanism shifts from predominantly vertical to mainly horizontal. Snowfall particles with α_d≥ 1 maintain a vertical collection advantage due to gravitational settling (accounting for over 75 % of the maximum potential total collection amount); while those with α_d≤ 0.2 shift to horizontal collection dominance under turbulence effects (accounting for at least 50 %). Within the u_* range of 0.06–1.18 m s⁻¹, the scavenging efficiency of 100–150 µm snowfall particle ensembles is optimal and exhibits a linear positive correlation with precipitation intensity. Among these, 100–150 µm snowfall particles are the most sensitive to variations in precipitation intensity.

A predictive model for snowfall particle motion (Eq. 23) is established, providing a quantitative theoretical basis for snowfall particle behaviour in wind-snow two-phase flow. A quantitative formula for wet deposition (Eq. 31) is also proposed, which can be utilized to quantify the enhancement effect of turbulence on the collection capability of snow particles. Consequently, the atmospheric dust wet deposition flux is accurately predicted, and significant application value is demonstrated in the fields of artificial dust removal and environmental assessment.

It should be noted that the actual shapes of snow particles are complex (such as dendritic, needle, plate, etc.), and differences in their shapes can lead to variations in terminal fall velocity, thereby affecting the threshold value of the dimensionless parameter α_d. Future work should conduct targeted studies on the α_d values corresponding to different snow particle shapes, in order to further expand the applicability of the model in real atmospheric environments.