Field investigation of 3D snow settling dynamics under weak atmospheric turbulence

Figures 2a and c illustrate the setup of our 3D PTV system. This system consists of four wire-synchronized cameras (Teledyne FLIR, FLIR Black Fly S U3-27S5C color unit with Sony IMX429 sensor: 1464 $\times$ 1936 pixels, 4.5 $\upmu\mathrm{m/px}$) strategically positioned around a light cone 5.5 m away, spanning a 90-degree angle range. This light cone is created by reflecting and expanding light from a searchlight using a curved mirror, similar to the setup used in our planar measurements. The cameras are tilted upward with 58-degree angles, leading them to image at a sample volume 10 m above ground. Each camera is then connected to its own data acquisition unit. These units are equipped with a board-level computer for issuing image capture commands to the cameras, a solid-state drive for storing both the system software and captured images, and a dedicated power supply. The cameras capture images with 2$\times$ decimation (732 $\times$ 968 pixels) to reach 200 Hz frame rate. As the standard checkerboard method cannot be applied to a field of view 10 m above ground after dark in the field, we use the wand calibration method described by Theriault et al. [43] for camera calibration. We use two colored light-emitting diodes (LEDs), set at a fixed distance apart on a carbon fiber rod, to act as the “wand” attached to an unmanned drone. This calibration process is conducted multiple times before and after each deployment to ensure the same field condition between camera calibration and snow particle imaging. We have developed custom-designed camera control software that synchronizes the image capturing process across all four cameras. The calibration of the cameras is conducted using the open-access software easyWand [43], which involves capturing images of the two colored LEDs as they move within the imaging volume. The software utilizes the trajectories of the two LEDs from all four cameras to conduct the calibration, resulting in a final reprojection error within 0.25 pixels.

For tracking snow particle trajectories, we utilize an open-source implementation of the shake-the-box (STB) method [44]. This version builds on the original STB method proposed by Schanz et al. [45], with enhancements specifically in the identification and removal of ghost particles. These improvements make it particularly suitable for our field data, which feature relatively high noise levels and a large field of view. This approach enables our system to capture snow particle trajectories within a considerable volume of approximately $4\times 4\times 6$ $\mathrm{m}^{3}$. The system achieves a spatial resolution of 6.3 mm/voxel and a temporal resolution of 200 Hz, allowing for detailed and precise tracking of snow particle movements. The tracked snow particle trajectories are then re-oriented as a group to have the average streamwise direction as the $x$ direction. Thus, the $y$ direction is defined as the spanwise direction, and the $z$ direction is defined as the vertical direction. From these trajectories, we obtain the Lagrangian velocity, $\bm{u}=(u_{x},u_{y},u_{z})$, the Lagrangian accelerations, $\bm{a}=(a_{x},a_{y},a_{z})$, and the resulting curvature, $\kappa=\|\bm{u}\times\bm{a}\|/\|\bm{u}\|^{3}$, where $\times$ represents the cross product. We use the second-order central difference method to calculate the Lagrangian velocity (first-order derivative) and acceleration (second-order derivative). The approximation introduces inherent errors, $O\left(\Delta t^{2}\right)$, which depends directly on the time step and is relatively small. However, the positioning errors of the snow particles can propagate and magnify in the velocity and acceleration calculation. As discussed by Schanz et al. [45] and Tan et al. [44], the iterative particle reconstruction, shake-the-box tracking, and trajectory filtering techniques significantly refine and reduce positioning errors. We quantify the root mean square of the difference between trajectory positions before and after filtering to be 0.3 pixels. This reduction potentially compensates for the positioning errors inherited from camera calibration, resulting in smaller errors in velocity and acceleration calculations. Consequently, the actual uncertainties in measuring velocity and acceleration are primarily influenced by the selection of filter length, which ranges from $45\pm 2$ frames (see Appendix A). This leads to an average acceleration uncertainty of 0.34 $\mathrm{m/s^{2}}$.

To complement our 3D PTV system, we also deployed a snow particle analyzer near the 3D PTV setup to assess the morphology and density of snow particles during each snow event (Figure 2b and d). All these measurements are crucial for accurately estimating the terminal velocity of snow particles in still air. As shown in Figure 2c, the snow particle analyzer employs a digital inline holography (DIH) system, which captures holograms of snow particles within a sample volume of $2.9\times 2.2\times 14.0$ $\mathrm{cm}^{3}$. This system achieves a spatial resolution of 14 $\upmu$m/pixel and a temporal resolution of 50 Hz. Through image analysis of the holograms, we obtain detailed information on particle size and shape, specifically the area equivalent diameter ($D_{\mathrm{eq}}$), major axis length ($D_{\mathrm{maj}}$), minor axis length ($D_{\mathrm{min}}$), area ($A_{e}$), etc. We also classify the shape of each particle into one of six types: aggregates, graupels, dendrites, plates, needles, and small particles. We define the characteristic particle size ($D_{p}$) as the area equivalent diameter for aggregates, graupels, and small particles, and as the major axis length for dendrites, plates, and needles. Additionally, a high-precision scale measures the weight of snow particles passing through the DIH sample volume, allowing us to estimate the average density of the particles. We also perform conditional sampling to achieve measurement of the density of individual snow particles.

For estimating the aerodynamic properties of snow particles, we follow the method proposed by Böhm [10]. This method involves calculating the Best number $X$ (also known as the Davies number), a dimensionless number that incorporates only the physical properties of snow particles and ambient air, and represents the equilibrium between gravity and drag forces. The Best number is defined as:

Noteworthily, unlike particle Reynolds number ($Re_{p}=\rho_{a}W_{0}D_{p}/\mu$) and drag coefficient ($C_{D}$), the definition of Best number eliminates the need to incorporate particle terminal velocity in the formulation, which is not readily available for complex snow particles. In equation 1, $\rho_{p}$ and $V_{p}$ are the density and volume of the snow particles, respectively, which are specific to the type of snow particle, as detailed in Li et al. [42]. Specifically, we approximate the complex snow particles as spheroids (graupels and small particles), combinations of small spheroids (aggregates), disks with thickness in correlation with their diameter (plates and dendrites), and thin cylinders (needles and columns). Such a method minimizes errors in volume estimation as compared to the typical spherical assumption used in the snow measurement community, resulting in the uncertainties of the volume estimation within 10$\%$ for snow particles with irregular shape (aggregates) and uncertainties of density within 20$\%$ for all demonstration cases in Li et al. [42]. $A_{e}$ is the effective snow particle imaged area, and $A$ is the circumscribed area of the enclosing circle or ellipse. Such an area ratio, $A/A_{e}$, serves as a simplified 2D measure of porosity and is instrumental in better predicting the drag of complex snow particles. As the snow morphological parameters are quantified by the snow particle analyzer while particles settle in various orientations, we assume that the ratio $A/A_{e}$ remains constant regardless of orientation. Finally, $\rho_{a}$ and $\mu$ are the density and viscosity of air, respectively.

Following the definition of the Best number, the drag coefficient of snow particles is modeled as a function of the particle Reynolds number, accounting for the unique morphology of snow particles. This approach indirectly incorporates the effect of snow particle density, which contributes to increasing the settling velocity. According to Stokes’ law for $Re_{p}\sim O(1)$, the correlation for the drag coefficient dependent on the particle Reynolds number is $C_{D}=24/Re_{p}$. However, the Stokes’ law becomes invalid as the Reynolds number increases, especially for complex snow particles. Researchers have made various attempts to model the drag coefficient of snow particles theoretically [10, 46, 47, 48]. As suggested by Böhm [10] and references therein, the drag coefficient of snow particles is modeled by considering the boundary layer surrounding the snow particles as a whole:

where $C_{0}=0.6$ is an inviscid drag coefficient and $\delta_{0}=5.83$ is a parameter controlling the evolution of the particle boundary layer, likely depending on the particle surface roughness, both empirically estimated. Equation 2 has the form of corrected Stokesian drag for a rigid sphere [49], but with different coefficients, modulating the transition from a linear drag at low $Re_{p}$ to a constant drag coefficient $C_{0}$ in the $Re_{p}$ independent regime. The effect of different snow morphologies is included in the $A/A_{e}$ term in equation 1, which is then used to predict the snow type specific drag coefficients, $C_{De}=\left(A/A_{e}\right)^{3/4}C_{D}$. The $C_{0}$ and $\delta_{0}$ parameters have been more recently updated, along with the dependency on the area ratio, by Heymsfield and Westbrook [50] and McCorquodale and Westbrook [51]. Additional corrections considering turbulent boundary layer, temperature, humidity, and accounting for different snow particle types, have been discussed in Khvorostyanov and Curry [46, 47] and Mitchell and Heymsfield [48].

As described in Böhm [10], we then obtain a semi-analytical and semi-empirical equation for the particle Reynolds number by combining equations 1 and 2:

The terminal velocity of the snow particles ($W_{0}$) in quiescent air is then calculated from the Reynolds number, $Re_{p}=\rho_{a}W_{0}D_{p}/\mu$. Once the terminal velocity is obtained, the aerodynamic particle response time is defined as $\tau_{p}=W_{0}/g$.

The analyses described have been meticulously applied to each snow particle type, leveraging the unique physical properties of individual particles, captured by the snow particle analyzer. Through a detailed examination of the collected holograms, we identify and classify each particle, subsequently analyzing their specific inertial properties, namely $D_{p}$, $\rho_{p}$, $A$, and $A_{e}$. By employing these properties within the Böhm model [10], we were able to estimate the aerodynamic properties of each particle. This rigorous method allows us to calculate the distribution and mean values of the terminal velocity ($W_{0}$) and drag coefficient for individual snow particles and specific snow types.

Over the course of the winter seasons from December 2021 to April 2023, we successfully carried out eight field deployments, encompassing a diverse range of environmental conditions. These deployments allowed us to study four major types of snow particles: aggregates, graupels, dendrites/plates, and needles/columns. We encountered wind speeds varying from a gentle 0.6 m/s to a more intense 8.4 m/s. Based on these wind speeds, we categorized the conditions into three turbulence levels: weak turbulence (wind speed less than 3 m/s with turbulent kinetic energy, TKE, below 0.3 $\mathrm{m}^{2}/\mathrm{s}^{2}$), moderate turbulence (wind speed between 3 and 6 m/s with TKE ranging from 0.3 to 2.0 $\mathrm{m}^{2}/\mathrm{s}^{2}$), and relatively strong turbulence (wind speed exceeding 6 m/s with TKE above 2.0 $\mathrm{m}^{2}/\mathrm{s}^{2}$). These turbulent properties were measured using the sonic anemometer positioned at a height of 10 m. Details of the estimation methods of these quantities can be found in Li et al. [40]. We use the second-order structure function of the streamwise velocity fluctuation to estimate the dissipation rate ($\varepsilon$). The Taylor microscale ($\lambda$) is then calculated as $\lambda=u^{\prime}\sqrt{15\nu/\varepsilon}$, where $u^{\prime}=\sqrt{\left(u_{x}^{\prime 2}+u_{y}^{\prime 2}+u_{z}^{\prime 2}\right)/3}$ is the representative scale of fluctuating velocity, and $\nu$ is the viscosity of air. Given the variety of snow particle types and wind speeds, our field data encompasses a total of 31 distinct conditions. To effectively separate the influences of snow morphology and atmospheric turbulence on snow settling velocity, a more systematic classification of the field snow and turbulence conditions is essential. We propose using the settling parameter $Sv_{L}=W_{0}/u^{\prime}$, which quantifies the relative impact of turbulence on snow gravitational settling [52, 53]. This parameter represents the ratio of the snow particle’s terminal velocity in still air ($W_{0}$) to the root-mean-square of the turbulent velocity fluctuations ($u^{\prime}$). A higher value of $Sv_{L}$ indicates that the influence of turbulence on the snow settling velocity is relatively minor.

In our analysis, we utilized the settling parameter, wind speed, and turbulent kinetic energy (TKE) as key criteria to categorize our 3D PTV and snow particle analyzer datasets. This approach led us to identify four distinct groups, which we labeled as ‘weak turbulence’ cases, with relatively smaller Taylor Reynolds number ($Re_{\lambda}$) and higher settling parameters ($Sv_{L}$) as shown in Figure 3. We assess the turbulence and micro-meteorological conditions for each group detailed in Table 1, employing estimation methods as outlined in the studies by Nemes et al. [6] and Li et al. [40]. Each group is dominated by one specific type of snow particle, which constitutes more than half of the snow population in the dataset. These types are aggregates, graupels, dendrites, and needles, as detailed in Table 2 and illustrated in the size and shape distributions in Figure 4. We leverage the capabilities of the snow particle analyzer, as detailed in Section II.1, to estimate these physical properties of snow particles. During a selected one-hour period characterized by dominant snow particle types, our analysis encompasses 200,000 holograms for each type of snow particle. This comprehensive dataset yields detailed information about approximately 28,000 aggregates, 13,000 graupels, 30,000 dendrites, and 21,000 needles. Complementing the snow particle analyzer measurements, our 3D PTV datasets include a total of 500 seconds of images for each dominant snow type, which are broken down into 50-second segments throughout the one-hour period selected. This rich dataset facilitates the identification of millions of snow particle trajectories, specifically around 322,000 for aggregates, 285,000 for graupels, 1,037,000 for dendrites, and 182,000 for needles, providing orders of magnitude more data than our previous studies. Specifically, our 3D PTV system measures the complete 3D velocity and particle acceleration components. The additional spanwise dimension of the data, compared to planar measurements, enables a thorough analysis of snow particle kinematics, including trajectory curvature and meandering. Furthermore, the integration of the 3D PTV system with the snow particle analyzer allows us to correlate the specific morphology of snow particles (e.g., size, shape, and type) with their settling behavior.

For the detected snow particles, their size and shape are measured through image analysis described in Section II.1. Given that these particles may present various orientations relative to the imaging plane, relying solely on the projected area (or equivalent diameter) falls short of providing a precise representation of the characteristic size of each particle, especially the non-spherical ones. We thus define the particle size as the equivalent diameter for aggregates and graupels, the major axis for dendrites (diameter) and needles (length), as detailed in Li et al. [42]. Upon closer examination, we observed notable differences among these types. Graupels and needles, for instance, tend to have a more uniform size distribution, with a smaller average size and standard deviation compared to aggregates and dendrites. Aggregates and dendrites, on the other hand, are generally larger, and their datasets include a mix of other particle types, resulting in a broader size distribution. We also analyzed the aspect ratio (i.e., the ratio between the minor and major axes lengths, $D_{\mathrm{min}}/D_{\mathrm{maj}}$) of these snow particles, defined as the ratio of their minor to major axis lengths, as measured by the snow particle analyzer. Graupels predominantly exhibit aspect ratios greater than 0.8, indicating their near-spherical shape. In contrast, the aspect ratios for the other types vary significantly from one, suggesting more anisotropic shapes. In this respect, note that the 2D holograms do not allow to accurately capture the averaged thickness of plate-like crystals due to the random particle orientation, unless further analysis is performed on selected particle images as in Li et al. [42]. Furthermore, we measured the average density ($\overline{\rho_{p}}$), together with the average particle size ($\overline{D_{p}}$) and aspect ratio ($\overline{D_{\mathrm{min}}/D_{\mathrm{maj}}}$), of the four datasets using the snow particle analyzer. Needles, being solid crystals with minimal riming, have the highest average density of 360 $\mathrm{kg}/\mathrm{m}^{3}$. Dendrites follow with an average density of 280 $\mathrm{kg}/\mathrm{m}^{3}$, as it is influenced by the gaps between branches, which contribute to the overall porosity of the particles. Graupels have an average density of around 220 $\mathrm{kg}/\mathrm{m}^{3}$, aligning with our previous measurements. Moreover, aggregates exhibit the lowest density of around 90 $\mathrm{kg}/\mathrm{m}^{3}$, as expected, attributable to their larger size and higher porosity.

This section presents an in-depth examination of the aerodynamic characteristics, including their terminal velocity, drag coefficient, and settling velocity, for each snow particle type. Table 3 consolidates key aerodynamic parameters derived from our analysis: the average settling velocity ($\overline{W_{s}}$) obtained through 3D particle tracking velocimetry (3D PTV), the average estimated still-air terminal velocity ($\overline{W_{0}}$) as outlined in Section II.1, the velocity fluctuation ($u^{\prime}$) and the Kolmogorov time scale ($\tau_{\eta}$) of the flow, the particle’s Stokes number ($St_{\eta}=\tau_{p}/\tau_{\eta}$), their settling parameter ($Sv_{L}=\overline{W_{0}}/u^{\prime}$), and the Froude number ($Fr_{\eta}=a_{\eta}/g$, where $a_{\eta}=u_{\eta}/\tau_{\eta}$ is the Kolmogorov scale acceleration). Needles exhibit the highest terminal velocity among all four types. With the same particle size, the cylindrical-shaped needles have the smallest projected area and the highest density, leading to larger terminal velocities. The Stokes number gauges the particle’s velocity response to sudden changes in flow, with values around one signifying a critical condition for turbulence-particle interactions. Settling parameters greater than one imply a weak influence of turbulence on the settling particles. The Froude number, a ratio of the characteristic flow acceleration ($a_{\eta}=u_{\eta}/\tau_{\eta}$) to gravitational acceleration, suggests that gravitational settling is more pronounced than the turbulence effect on the particles [54]. Comparatively, the settling velocity enhancements from the terminal velocities are moderate, ranging up to 32$\%$ for aggregates, 13$\%$ for dendrites, 4$\%$ for needles, and 3$\%$ for graupels. These findings indicate that the turbulence effects (e.g., preferential sweeping and loitering) on particle settling is generally weak under the examined conditions. Variations in settling enhancement across snow types may be largely attributable to differences in particle size, shape, and density.

Figure 5 presents a comparative analysis of the probability density functions (PDFs) for settling velocity ($W_{s}$) and estimated still-air terminal velocity ($W_{0}$) across various snow particle types. The estimate of $W_{0}$ is based on the Best number, $X=C_{D}Re_{p}^{2}$, which does not directly depend on the settling velocity of the snow particles, but rather on their physical properties and the ambient air. Following Böhm [10] approach, summarized by equations 1-3, we estimate the terminal velocity from measurable geometric and inertial properties by the snow particle analyzer. The PDFs for graupel, which are nearly spherical in shape, exhibit a close overlap between the settling and terminal velocities, indicating a minimal influence of turbulent eddies on their settling dynamics. Note that the ‘dent’ in the distribution of the terminal velocity of graupels in Figure 5b reflects their size distribution. On the contrary, for the other snow types—characterized by non-spherical geometries—the PDFs diverge despite the mean settling and terminal velocities for needles displaying only a 4$\%$ discrepancy. This variation suggests that the aerodynamic behavior of non-spherical particles is considerably affected by the randomization of their orientation due to flow disturbances and unsteady behavior. In quiescent conditions, particles falling stably tend to orient themselves to maximize the aerodynamic drag (i.e., preferential orientation), potentially due to the inertial forces of the surrounding media, presenting their maximal cross-sectional area perpendicular to the fall direction [22, 55]. However, in turbulent conditions, such a preferential orientation is not appreciable [55, 56], and the varying orientations result in a reduced effective cross-sectional area, potentially leading to an increased average settling velocity for non-spherical particles. Furthermore, while the settling velocity distributions for different snow types approximate a Gaussian profile, the estimated terminal velocities are rather skewed. This asymmetry arises from the inherent size distributions of the snow particles, which are typically modeled using a gamma distribution [57]. We also acknowledge the potential sampling differences between the 3D PTV measurements (likely under-representing the finest size fraction) and the snow particle analyzer data collection.

Historical studies have demonstrated that the terminal velocity of snow particles exhibits a size-dependent characteristic, since the early research by Nakaya and Terada [15], Heymsfield [16], and Locatelli and Hobbs [9] fitting empirical data to establish a particle-mass-based approach to the settling. Specifically, Locatelli and Hobbs [9] conducted a thorough investigation of various snow particle types, deriving power-law empirical formulas to represent the size-dependent terminal velocity, expressed as $W_{0}=a{D_{p}}^{b}$, where $a$ and $b$ are constants that differ based on the snow particle type, based on shape and density. For our analysis, we employed the formulas relevant to aggregates of unrimed radiating assemblages of dendrite ($W_{0}=0.8D_{p}^{0.16}$), conical graupel ($W_{0}=1.2D_{p}^{0.65}$), rimed dendrites ($W_{0}=0.62D_{p}^{0.33}$), and rimed columns ($W_{0}=1.1L^{0.56}$, where $L$ is the length). Such power-law equations can be empirically obtained by fitting the size distributions in Figure 4 with the settling velocity in Figure 5. We optimize the linear coefficient with the same exponent to impose the same mean and similar distribution of the settling velocity for each snow type, and thus obtained the empirical equations: $W_{s}=1.45D_{p}^{0.16}$ for aggregates, $W_{s}=1.2D_{p}^{0.65}$ for graupels, $W_{s}=0.92D_{p}^{0.33}$ for dendrites, $W_{s}=1.66L^{0.56}$ for needles. We thus obtain the same equation for $W_{0}$ and $W_{s}$ for graupels, suggesting close alignment in the mean values and distributions between the measured settling velocities and the estimated terminal velocities, the same as predicted using equations from Böhm [10]. This also confirms negligible effects by the specific atmospheric turbulence conditions monitored during the settling of graupels. In contrast, for other non-spherical types of snow, the linear coefficients for $W_{s}$ are higher than those of $W_{0}$. This discrepancy highlights the morphology effects that modulate the settling velocity of these non-spherical snow particles, with potential turbulence effects considering the varying particle orientation because of turbulence disturbances, and the production of the Stokes number and settling parameter reaching critical condition [$St_{\eta}Sv_{L}\sim 1$, as suggested in 52, 53].

To better model the terminal velocity, it is important to quantify the aerodynamic drag of snow particles for various morphological types. In Figure 6, we present the mean drag coefficients and mean Reynolds numbers, estimated using the average particle size and measured settling velocity. The error bars indicate the variability of these quantities, reflecting the distribution of snow particle sizes and settling velocities as represented by their standard deviations. The drag coefficient is calculated as $C_{De,\text{ mean }}=2\overline{\rho_{p}}\overline{V_{p}}g/\left(\rho_{a}\overline{W_{s}^{2}}\overline{A_{e,\mathrm{max}}}\right)$, where $\overline{\rho_{p}}$ is the average snow particle density, $\overline{V_{p}}$ is the average particle volume [different expressions for different snow particle types defined in 42], $\overline{W_{s}^{2}}$ is the mean square settling velocity; and $\overline{A_{e,\mathrm{max}}}$ is the average maximal projected area of the measured snow particles (e.g. a flat falling dendrite, see Appendix B). The snow particles have an average Reynolds number ($Re_{p,\mathrm{mean}}=\overline{W_{s}}\ \overline{D_{p}}/\nu$) on the order of 100, agreeing with typical field measurements [50]. The drag coefficients for aggregates, graupels, and needles agree well with the model predictions from Böhm [10, equation 2], $C_{De}=\left(\overline{A/A_{e}}\right)^{3/4}C_{0}\left(1+\delta_{0}/Re_{p}^{1/2}\right)^{2}$, as presented by the dotted lines. Note that the average area ratio, $\overline{A/A_{e}}$, is calculated from the snow particle holograms for each snow type, and it is necessary to rescale the generalized drag equation 2 to the specific snow morphologies [10]. As graupels show more sphere-like features, their drag coefficient leans towards that of spheres [corrected for high Reynolds number by 49]. Despite the smaller terminal velocity for the non-spherical particles considering the particle orientation, the drag coefficient is well-predicted by the Böhm [10] model for aggregates and needles. Potential contamination from other types ($\sim$ 20$\%$ after filtering out particles with close to an aspect ratio of one) within the datasets might lead to the mismatch between the PDFs of the terminal velocity and measured settling velocity for needles. The enhanced settling for aggregates could be a combined result of particle orientation, weak turbulence enhancement considering the critical condition of $St_{\eta}Sv_{L}\sim 1$, and contamination from other types that do not align with the statistically dominant group, as shown in Table 2. Moreover, the dendrites show, on average, a higher drag coefficient as compared to the other types, potentially due to their large frontal area and higher density. Such a discrepancy could also explain the underestimation of the terminal velocity of dendrites by the equations from Böhm [10], considering the higher, on average, settling velocity. We further compare our measurements with laboratory experiments by Tagliavini et al. [12], with squares representing the aggregates (AgCr77, Ag15P1, AgSt18), pentagrams representing the dendrite (D1007), and left-pointed triangles representing the columnar snow (CC20Hex2). Our measurements agree well with the laboratory experiments, with the dendrite type showing a larger drag coefficient due to its disk-like shape. Such observations provide insights for snow settling modeling, especially for the predominantly dendrite snow events, as they show a large deviation from the Böhm [10] model prediction.

Our comprehensive analysis elucidates the distinctive settling behaviors of snow particles with various morphologies, addressing the questions raised at the beginning of our results section. Concerning the influence of morphology on snow aerodynamic properties, we observe that the response times for all snow particles are broadly similar, averaging around 0.1 seconds, on the same order of the intercept in the acceleration autocorrelation function. However, needles exhibit a marginally increased response time attributed to their higher density, particularly when compared with the density of the surrounding air. This higher density contributes to the needles’ higher average terminal velocity in still air. Furthermore, the empirical models by Böhm [10] well predict the drag coefficients for aggregates, graupels, and needles, while dendrites emerge as anomalies, exhibiting drag coefficients exceeding model predictions, likely due to their large frontal area and oblate, disk-like, geometry. Notably, although dendrites and needles have similar aspect ratios, as measured by the snow particle analyzer, the dendrites are disk-like, oblate spheroids, while the needles are columnar, prolate spheroids.

The different aerodynamic properties and morphologies of these snow particles have strong effects on their settling kinematics. Specifically, dendrites display unique behavior compared to other types. Their non-linear drag and substantial frontal area result in the most prominent acceleration fluctuation magnitude, occurring at relatively low frequencies. While the acceleration probability density function (PDF) for dendrites closely resembles that of a fluid parcel in turbulence, the acceleration auto-correlation function indicates a slow response to the rapid fluctuation of the flow velocity. Conversely, needles exhibit minimal acceleration fluctuation magnitude, suggesting relatively large inertia and a tendency to avoid intense cross-flow drag and convoluted trajectories. Yet, the acceleration autocorrelation function indicates a moderate rate of change in the direction of acceleration. Such different behaviors from that of dendrites are possibly due to their streamlined shape aligning with fluid flow structures, as for fibers in turbulence. Overall, acceleration statistics appear to correlate with the shape factors of these particles considered as spheroids, with dendrites and needles representing the spectrum’s extremes and aggregates positioned in between.

Finally, to answer the third question, it becomes apparent, from our analysis above, that the combination of turbulence and non-spherical particle morphologies can modulate the particle settling velocities even under weak atmospheric turbulence. For the cases investigated here, we hypothesize that dendrites exhibit an enhanced settling velocity (Figure 5c) that is due to an underestimation of the drag coefficient in the model by Böhm [10]; graupels settling velocity is well predicted even though spherical particles smaller than the Kolmogorov scale close to critical Stokes conditions were expected to exhibit settling velocity enhancement. The significant cross-flow velocities (considering the large settling parameter $Sv_{L}$) experienced by graupels may suggest that preferential sweeping was not the only mechanism in play during settling. Aggregates drag coefficient is well captured in the still air model, implying that the observed enhanced settling is likely due to combined effects of anisotropic particle orientation and the weak atmospheric turbulence. The observed conditions are marked by $St_{\eta}Sv_{L}\sim 1$ for which settling enhancement has been predicted and observed [52, 53]. Disentangling turbulence and morphology effects is challenging because turbulence-induced disturbances alter the preferential orientation (i.e., particles with their largest projected area facing the settling direction) of stably falling particles. The meandering motions of the non-spherical particles, whether fluttering or tumbling, likely affect their orientation, reducing their average projected area compared to a steady settling and thus enhancing settling velocity. However, direct measurement of particle orientation during settling is technically challenging and beyond our current capability.

Thus, we investigate the interconnection between the meandering motion and the vertical acceleration along the trajectories of snow particles. Recognizing that these fluctuations might not be perfectly synchronized (owing to the particles’ inertial response and the variability in drag force related to changes in projected area and settling velocity), we have considered a slight phase shift between the varying spanwise location $y(t)$ and vertical acceleration $a_{z}(t+\tau)$ along the trajectories. This adjustment aims to align the locations of maximum meandering with the smallest projected area, which typically corresponds to greater downward acceleration. To maintain the integrity of the correlation, we limit the phase shift to 0.15 times a quarter of the meandering period to avoid creating an inverse relationship between vertical acceleration and meandering motion. Note that this time lag ($\tau$) is on the order of $0.1\tau_{p}$, and close to the estimated Strouhal number for the disk-like particles. Thus, during a fraction of the anisotropic particle rotation and the corresponding translational response time, the particle is experiencing a reduction in drag area, leading to its acceleration downward. Figure 13a and b demonstrate that, following this phase shift, the greatest downward accelerations predominantly occur at the furthest extent of the spanwise meandering motion ($y^{\prime}_{a,\mathrm{max}}\sim 1$). Conversely, the least downward acceleration—or even upward acceleration—tends to happen near the central position ($y^{\prime}_{a,\mathrm{min}}\sim 0$), where anisotropic snow particles are likely to have their maximum projected area facing downwards. Subsequently, we calculate the correlation coefficient between vertical acceleration and spanwise position during the snow particles’ meandering motion. The results reveal substantial positive correlation coefficients for dendrites ($\overline{\sigma_{y,a}}=0.58$) and aggregates ($\overline{\sigma_{y,a}}=0.45$), which is consistent with their observed enhanced settling velocities. Needles display a moderate correlation coefficient ($\overline{\sigma_{y,a}}=0.33$), reflective of their anisotropic shape. Graupels, however, exhibit a low average correlation coefficient of $\overline{\sigma_{y,a}}=0.15$, as changes in particle orientation are not expected to significantly affect the drag force.