Wake characteristics of complex-shaped snow particles: Comparison of numerical simulations with fixed snowflakes to time-resolved particle tracking velocimetry experiments with free-falling analogs

The shape and the orientation of falling ice particles influence snow precipitation and are crucial for weather prediction and polarimetric radar measurements, respectively. Snow crystals display different sizes and shapes in nature (Kikuchi et al., 2013), and for snowflakes larger than few micrometers, Stokesian dynamics is not valid anymore (Brady and Bossis, 1988; Westbrook, 2008) [$Re≫$ 1, where $Re=(utDmax)/ν$⁠, with Re being the particle Reynolds number, u_t being the snowflake terminal velocity (m/s), D_max being the maximum span of the particle normal to the fall direction (m), and ν being the kinematic viscosity of air (m²/s)]. Hence, the interaction with the surrounding air becomes non-linear (large D_max implies large Re) (Abraham, 1970). Large snow particles (⁠$Dmax≫$ 100 μm) exhibit tangled falling trajectories due to the combined effects of irregular particle shape, size, and flow regime (i.e., Reynolds number).

A large variety of experimental techniques are nowadays available to investigate complex flows. PIV (particle image velocimetry) appeared more than three decades ago, and it became widely used as a quantitative flow measurement solution (Adrian, 1991; Fu et al., 2015; Chen and Chang, 2018). Its success can be explained by its ability to visualize and quantify the instantaneous velocity field within an entire plane of the flow domain (Brossard et al., 2015). When studying turbulent flow, there is the need of recording a three-dimensional velocity field and one can rely on volumetric techniques, such as tomographic PIV (Elsinga et al., 2006; Scarano, 2013; Qureshi et al., 2021), because of their capability to tackle the instantaneous 3D organization of the coherent structures in turbulence, in particular, for three-dimensional wake flow (Bearman, 1997; Chen and Chang, 2018; Terra et al., 2019). Lagrangian particle tracking methods are also available, alongside PIV (Wieneke, 2013). Among these techniques, time-resolved, three-dimensional PTV (4D-PTV) is becoming increasingly popular, especially after the development of the modern velocimetry algorithms such as shake-the-box by Schanz et al. (2016), which allows the use of much higher seeding densities and, thus, a measurement resolution comparable to that of PIV. Next to measuring the flow velocity field around falling snowflakes based on tracers, stereoscopic imaging can also be employed to measure the position and orientation of the snowflakes themselves. This was done by McCorquodale and Westbrook (2020b) to investigate the falling behavior and drag coefficient of freely falling, three-dimensional (3D)-printed snowflake analogs in a mixture of water and glycerol. This study also aimed to improve former parameterized models for snowflake terminal velocity (Mitchell and Heymsfield, 2005; Heymsfield and Westbrook, 2010), which are generally constrained by their specific experimental or numerical setup and confined to moderate Reynolds numbers (⁠$Re≲$1000).

Flow visualization methods are often combined with numerical approaches (DNS, LES, or RANS) to validate the latter, when studying turbulent flows. For instance, such a comparison has been used to test the applicability of LES to turbulent wakes of a rectangular cylinder (Kuroda et al., 2007), human-shaped manikin (Luo et al., 2018), and wind turbines wakes (Lignarolo et al., 2016; Wang et al., 2019). In these works, the experimental setup accurately reproduces the numerical one, and the object geometries are usually quite simple (spheres, cylinders, or airfoils), often at large Reynolds numbers (⁠$Re∼$ 4000). Very few works (Luo et al., 2018; Terra et al., 2019) that combine experiments and computational models are to be found on wakes of complex-shaped objects.

The falling behavior of non-spherical particles has been widely explored both experimentally and numerically for many years. Field et al. (1997) studied the behavior of falling disks and identified different falling motion (steady, periodic, and chaotic) according to dimensionless moment of inertia and Reynolds number. Other works on disks include the numerical investigation of Fabre et al. (2008) and Auguste et al. (2013), while and Vincent et al. (2016) and Esteban et al. (2019) experimentally investigated the influence of central holes in falling coins and the wake coherent structures of falling disks, hexagonal, and squared plates, respectively. The importance of investigating the wake characteristics of spherical and non-spherical particles lies in the influence of the latter on the drag and falling trajectory (Fernandes et al., 2007; Kim et al. 2018). Snow particles may exhibit fractal-like geometries (Stein et al., 2015; Dunnavan et al., 2019) that can affect their drag and wake. In this view, the studies by Nedic et al. (2013) and Nedic et al. (2015) looked into the impact of the fractal dimension of plate-like geometries on the drag coefficient, wake size, and vortex shedding. By increasing the fractal dimension, i.e., longer particle perimeter (while keeping the frontal area constant), the drag coefficient, and the shedding intensity increase, while the wake size decreases, together with the vortex shedding energy. Other works analyzed non-spherical particles such as plate-like crystals (Cheng et al., 2015) and columnar crystals (Hashino et al., 2016). Hashino et al. (2016) and Toupoint et al. (2019) found a strong dependence of the falling attitudes on the aspect ratio of the investigated geometry and described in detail the changing of forces and torques in relation to the particle falling motion. Furthermore, the wake that forms behind a snowflake impacts the snow crystal growth (or sublimation) (Wang and Ji, 1997; Ji and Wang, 1999) and its interactions with neighboring particles and the surrounding air (Nemes et al., 2017; Li et al., 2021). Notwithstanding this extensive literature, the high complexity of snow particles' geometries is only partially addressed in previous studies, and the flow regimes are generally moderate (⁠$Re≲$ 10³) due to limitations in the employed experimental techniques or in the computational capability of the chosen numerical approach. With regard to the first limitation (lack of studies tackling complex shapes), Zeugin et al. (2020) investigated realistic snow particles in the Stokes regime and proposed sphericity-based relations to estimate the particle aerodynamic coefficients and to set velocity up to Re ∼ 10. With regard to the second limitation (most studies tackle low Re), in our previous work, we proposed a novel approach based on a DDES model and the particle inertia tensor to accurately predict snow the particles drag coefficient and terminal velocity (Tagliavini et al., 2021a) at moderately high Re, and we investigated the influence of the particle shape on the wake flow and the drag coefficient of complex-shaped snowflakes at a fixed position (Tagliavini et al., 2021b).

In this work, we combine 4D-PTV measurements of freely falling, 3D-printed snowflake analogs and delayed-detached eddy simulation (DDES) that solves for the airflow around a fixed snow particle, previously validated for drag coefficient prediction (Tagliavini et al., 2021a). The measurements are compared with the numerical results of the same snowflake shape (kept fixed in the computational domain) at the same Reynolds number to assess the ability of the fixed-particle model to reproduce the velocity field of the same free-falling particle. In comparison, time- and space-averaged, non-dimensional flow quantities in the wake are evaluated. The aim is to, first, cross validate the two approaches for low Reynolds number cases, where close agreement of the wake characteristics is expected, and, second, to assess how strongly the unsteady falling motion perturbs the average wake patterns as compared to a fixed particle at the same Re. To do so, following the criteria employed in our previous works (McCorquodale and Westbrook, 2020b; Tagliavini et al., 2021b), two flow regimes are identified according to the type of falling motion displayed by the snow particles during laboratory tests: a flow regime at which the particles fall steadily (low Re) and another one at which the particles exhibit unsteady falling behavior (moderate/high Re). The combination of 4D-PTV data of free-falling particles and a DDES model will allow us to shed light on the capabilities and the limitations of both approaches and better understand differences between the wakes of fixed and freely falling particles.

In Sec. II, the experimental setup, the velocity field reconstruction, and the observations concerning the particle falling behavior are first described. Subsequently, the computational model is presented (Sec. II C), followed by the preparation of the datasets for comparison and analyzed quantities (Sec. II D). In Sec. III, the results from the comparison are provided and discussed. Initially, the steady falling cases are presented (Sec. III A), for which good agreement is found. Thereafter, the particle with unsteady falling attitudes are taken into account (Sec. III B). To facilitate the comparison, additional experimental data, whose reconstruction is done by including the particle rotation (Sec. II A 1), and the numerical data filtered on larger grids (2 and 4 mm) are also considered for unsteady cases. The conclusions from the comparison are drawn in Sec. IV.

This section includes the description of the experimental setup, the velocity field reconstruction, the laboratory observations, and the main features of the DDES model. Conform to our previous works (Tagliavini et al., 2021a; 2021b), the numerical model exploits the information about the Reynolds number and the final orientation of each snow particle from the experiments. The numerical solution is then compared with the measurements of the same particle geometry. The datasets preparation is also described in Sec. II D.

A schematic view of the experimental setup is shown in Fig. 1. The experiments are conducted in a transparent acrylic box, the cross section of which is a regular decagon with a circumscribing diameter of 400 mm. The tank is filled with uniform mixtures of water and glycerol, in which the volume fraction of glycerol is set between 0% and approximately 50% to a total liquid depth of approximately 1.44 m. The density and temperature of the mixtures are measured before each experiment; these measurements are used to estimate the viscosity of the mixture (Cheng, 2008; Volk and Kähler, 2018). The parameters of the water tank match with the ambient. Ambient humidity has no effect on the present design of the experiment. The temperature of the glycerin water mixture was noted for each set of experiments, and the values exhibited minimal variation. Across the period of several weeks in which experiments were conducted, the temperature of the glycerin water mixture was observed to lie between 14.3 °C and 17.5 °C. Since the 3D-printed ice-particle analogues are placed in the water for some time, equilibrium is obtained between the fluid and the particle analogues. During experiments, the analogs are submerged close to the top of the tank and placed inside a hemispherical cup. To initiate the experiment, the cup is tilted, releasing an analog with a random orientation. Care is taken to remove any air bubbles from the analogs before they are released. Qualitative observations indicate that following a short transient phase (typically over a vertical extent of less than 0.3 m), the subsequent fall of the analogs is independent of orientation upon release. That is, the analogs either (i) exhibit a clear preferred orientation in free fall that is insensitive to release conditions, at low Reynolds numbers, or (ii) exhibit a chaotic trajectory (moderate/high Re), the details of which are unreproducible even should the particle be released from a nominally identical orientation. To ensure quiescent conditions during experiments, the release of analogs is separated by a time period of at least 15 min. Measurements of instantaneous fluid velocities are obtained using four-dimensional particle tracking velocimetry (4D-PTV), employing the shake-the-box method (Schanz et al., 2016), applied to a volume approximately 10 × 10 × 10 cm³ in size approximately 1.25 m below the surface of the fluid. Small, tracer particles (fluorescent red/orange polyethlene microspheres with the diameter range 125–150 μm) are added to the water volume and thoroughly stirred. Seeding particles are selected to ensure they were approximately neutrally buoyant within the water–glycerol mixtures; tracer particles with densities of 1000, 1090, and 1200 kg/m³ are used in water–glycerol mixtures of approximate densities of 999, 1088, and 1144 kg/m³, respectively. A LED flash-light (LaVision LEDFlashlight 300 Blue), which emits blue light with a peak intensity at a wavelength of 451 nm, is used to illuminate the tracer particles located within the central region of interest (see Fig. 1). When excited by the blue light of the LED, the tracer particles fluoresce and emit light with a peak intensity at a wavelength of 606 nm. As analogs fall through the region of the interest, the motion of the illuminated tracer particles, as they are advected by the wake, is recorded using a series of high-speed digital cameras (VC-Phantom Miro M120) positioned to point horizontally into the tank's interior [see Figs. 1(b) and 1(c)]. The camera lenses (Zeiss Planar T^* 50 mm f/1.4 ZF.2), attached to the cameras via Scheimpflug adapters to apply perspective corrections, are fitted with a filter (Hoya MHC YA3 Orange) that permits the transmission of the red light emitted by the fluorescent tracer particles but blocks the blue light emitted by the LED flash-light. The cameras' capture is synchronized with the pulse of the LED-flashlight. The frame-rate of the flash-light and the cameras is varied according to the terminal velocity of the particle analogs; images are recorded at between 13 and 200 frames per second (at 1920 × 1200 pixel resolution). Tracer particle trajectories are determined using the LaVision's Shake-the-box algorithm. Tracer particle velocities are gridded to reconstruct velocity fields using a sub-volume size of 64 × 64 × 64 voxels (of order 0.6 × 0.6 × 0.6 cm³ in size), overlapped by 75% to achieve 16 voxel spacing between velocity vectors. This results in a physical spacing between velocity vectors of approximately 0.15 cm. A velocity vector is only returned in sub-volumes in which a particle track is present; spatial fitting, or interpolation across sub-volumes, is not applied. The velocity data are calculated and analyzed relative to the right-handed coordinate system (x; y; and z); here, z denotes the vertical axis, and (x; y) are the horizontal coordinates perpendicular to the fall-direction (see Fig. 1). The corresponding velocity components are denoted (u; v; and w). A custom algorithm (TRAIL) (McCorquodale and Westbrook, 2020a) is also used to digitally reconstruct the trajectory and orientation of the 3D-printed analogs from the images recorded by the digital cameras, which enables the falling analog to be co-located within the velocity fields measured through the Shake-the-box algorithm.

Since the particles employed in the simulations are fixed and the 4D-PTV measurements are performed for freely falling particles, it is necessary to transform the experimental data to facilitate the comparison. To do this, the coordinate system of each frame of the experimental data is adjusted to match the one used in the simulations (see Fig. 3), and the flow field is translated so that the center of mass of the particle is located at the origin. Next, the flow field is rotated such that the orientation of the reconstructed particle in the experiment matches the orientation in the simulation. Although it is known a priori how the particle is oriented relative to the fall direction (see Sec. II B), the orientation within the horizontal plane is not controlled in the experiment, and indeed complex particles are observed to slowly rotate over time. To match this to the fixed orientation in the simulation, different rotations around the fall direction are tested. The closest match between the 2D projections in the horizontal plane is selected, and then the complete flow field for that frame is rotated to match the simulation particle orientation. Finally, a sequence of such co-moving, co-rotating flow field measurements is time-averaged together. Following these steps, a second experimental dataset is created and added to the comparison to appraise the influence of two different types of measurements reconstruction approaches (with and without accounting for the particle rotation).

The snow particles analyzed in this study are shown in Fig. 2(a). Given that snowflakes exhibit a large variety of sizes and shapes in nature (Kikuchi et al., 2013), we took into consideration shapes that are examples of the most common classes of ice particles: monocrystals (D1), multi-habit crystals (CC), polycrystals (MR), and aggregates (AG). As a consequence of their diverse shapes, each snow particle displayed a peculiar falling behavior that is concisely described in this section within the Reynolds number range investigated in this study (60 $≲Re≲$1700).

The plate-like crystal (D1) falls steadily during the laboratory observations throughout the Re range. The largest projected area of the geometry is orthogonal to the falling direction (i.e., the maximum principal moment of inertia is aligned to the fall direction), as shown in Fig. 2(a). AG displays a spiraling trajectory even at low Reynolds numbers, although its falling behavior can be considered steady up to $Re≲$ 1000 in view of the small changes in the orientation (Fig. 2). It deviates from a steady motion as the Reynolds number increases and spirals around a vertical axis parallel to the falling direction at high Re (periodic motion). For what concerns CC, a steady falling attitude is present up to $Re≲$ 70; after that, the capped columns exhibits helical motions, whereby the long-axis of the particle is approximately aligned parallel to the fall direction. Initially, the helical motions show a clear periodicity, but as Re increases, the helical motions become chaotic. The same behavior is also seen in Kim et al. (2018) for rigidly linked disks. MR falls steadily up to $Re≈$ 250, then its motion becomes chaotic with huge variation in its orientation. More details on the falling behavior of each particle can be found in McCorquodale and Westbrook (2020b). To better illustrate the particles falling behavior, Fig. 2 reports the variations of the angle between the falling direction and the largest (α), the intermediate (γ), and the smallest (β) principal axis of inertia for each shape.

The computational model utilizes the Reynolds numbers and the observed orientations from experiments to define the inflow velocity and the particle position with respect to the flow direction, as described in Sec. II C. Figure 2(a) shows the final orientation of each snow particle. This orientation matches the one achieved by the particle after falling steadily for some time and after reaching its terminal velocity (Tagliavini et al., 2021a). This position is employed in the computational model also for cases in which particles are falling unsteadily (i.e., at moderate/high Reynolds numbers). This is motivated by the fact that the low Reynolds number orientation (i.e., steady falling) lies within the range of the recorded orientations even when the particles display unsteady falling attitude. Nevertheless, we showed that for certain particles, at high Re, the drag coefficient evaluated for these orientations differs significantly from the ones exhibited by the snowflake during experiments (Tagliavini et al., 2021a) because those particles frequently visited positions that are far from the final orientation. As a consequence, in our second work, we also explored “extreme orientations” [see Tagliavini et al. (2021b) for a detailed description]. In particular, in this study, we retain the extreme orientation of CC at $Re=$ 488, i.e., the position for which the angle between the long axis of the columnar crystal and the vertical fall direction is the greatest.

The computational model presented here is built up with Open source Field Operation And Manipulation (OpenFOAM 4.1), a C++ code based on the finite volume method (OpenFOAM, 2017). The model solves for the airflow domain, where a fixed, realistic snow particle is positioned (at its final orientation) at the origin of the domain coordinates (Fig. 3) and impinged by the flow. The fluid motion is attained by solving the transient Navier–Stokes equations

in which u is the flow velocity (m/s), p is the pressure (Pa), μ is the dynamic viscosity of the fluid (Pa $·$ s), ρ is the fluid density (kg/m³), and f are any external forces per unit mass (N/kg). Newton's second law of motion describes the one-way interaction between the air and the snowflake

with n being the unit vector pointing out of the particle surface A and τ being the viscous stresses (Pa). The computational domain and mesh are shown in Fig. 3 and are based on general guidelines for computational fluid dynamics (Ferziger and Perić, 2002). Both domain size and grid convergence studies are performed by means of Richardson's extrapolation (Roache, 1994). The volume of each snowflake is scaled to be equal to that of a sphere with a diameter of 1 cm, and the dimensions of the computational domain are defined as a function of the volume-equivalent sphere diameter (Tagliavini et al., 2021a). For each case, the Reynolds number matches the one of the laboratory experiments (Sec. II A) and is used to evaluate the uniform air speed set at the inlet boundary

where Re is the desired particle Reynolds number (from experiments), D_max is the maximum dimension of the snowflake orthogonal to the flow direction (m), $u∞$ is the inlet velocity (m/s), and ν is the kinematic viscosity of the fluid (m²/s). At the outlet, a zero static pressure is established. The snowflake is treated as a fixed wall, on whose surface a no-slip condition is imposed. A symmetry boundary condition is chosen for the lateral boundaries of the domain (Fig. 3). More details about the computational domain size, grid, and initial and boundary conditions can be found in Tagliavini et al. (2021a). In relation to the coarse-graining operations performed on the numerical data, we highlight that the grid size in the wake region [refinement region (b) in Fig. 3] and in the proximity of the snow particle [refinement region (c) in Fig. 3] is $≈$ 5 × 10^–4 m and $≈$ 3 × 10^–5 m, respectively.

The numerical model employs a DDES (Spalart et al., 2006) turbulence modeling approach (hybrid URANS-LES) with the Spalart–Allmaras turbulence closure model for $ν̃$ (modified eddy viscosity) for the URANS calculation (Spalart and Allmaras, 1992). The turbulence model was previously validated in Tagliavini et al. (2021a), where more details about the implementation, the solver, and the numerical schemes can be found.

The flow regimes (Re) at which the simulations are carried out are chosen such that they include either cases of steady (low Re) and unsteady (moderate/high Re) falling behavior. Namely, D1 is simulated at $Re=$ 206 (steady), AG at $Re=$ 62 (steady), 1691 (unsteady, periodic), MR at $Re=$ 1110 (unsteady, chaotic), and CC is evaluated at $Re=$ 488 (unsteady, chaotic). All the particles are simulated at their final orientation, except CC for which the additional extreme orientation is also considered (Sec. II B). For the summary of the different particle orientations and Re investigated, see Table I.

Before evaluating the quantities employed in the comparison, the numerical velocity field is modified to a pseudolab frame-of-reference by subtracting the inflow velocity $u∞$ to the vertical component of the velocity. In this way, the experimental and numerical flow field are much easier to compare. Moreover, the coordinates of both datasets are normalized using the particle maximum dimension D_max (m) orthogonal to the flow direction. As a consequence, the non-dimensional coordinate in the flow direction $x*$ is evaluated using the following formula:

where x is the coordinate in the flow direction (m). We then calculate the time-averaged, non-dimensional velocity $u*¯$ defined as

in which $u¯$ is the time-averaged velocity vector (m/s) and u₀ is the inflow velocity equal to $u∞$ (m/s) in the computational case and the snow particle terminal velocity u_t in the experimental one (m/s), respectively. Using the normalized coordinates and the time-averaged, non-dimensional velocity $u*¯$⁠, the non-dimensional velocity gradient $grad(u*¯)$ is obtained, whose L₂ norm (magnitude) is also calculated.

The snow particle wake is identified by an iso-volume of the time-averaged, non-dimensional velocity magnitude (⁠$u*¯$⁠). The iso-volume is delimited by threshold on $u*¯$ set at 90% of its maximum value. Given that the field of view captured in the experiments extends only a few D_max from the particle, the comparison focuses on the quantities in the near wake (⁠$x*≤$ 3 D_max). Furthermore, the particle center of mass is fixed at the coordinate system origin in both flow fields. In this way, the geometry of the wake in the numerical data overlaps one of the laboratory fields [Figs. 4(f) and 4(g)]. Once the wake is outlined, the time-averaged, non-dimensional velocity magnitude (⁠$u*¯$⁠) and the L₂ norm of the time-averaged, non-dimensional velocity gradient $grad(u*¯)$ are spatially averaged over different surfaces orthogonal to the flow direction along the wake iso-volume (⁠$U*$ and $gradU*$⁠, respectively, for which the time-averaging symbol is omitted to simplify the notation). For comparison between the numerical data and the measurements accounting for the particle rotation to speed up the measurements reconstruction, only the vertical component of the time- and space-averaged, non-dimensional velocity is taken into account (⁠$Ux*$⁠, see Sec. III B 1). The distance between these surfaces is $Δx*=$ 0.25 D_max, and the starting point is at $x*=$ 0.5 D_max for AG and D1, while for MR and CC is at $x*=$ 0.75 D_max due to the larger extension of those particles along the x axis (streamwise direction). Figure 4 depicts the comparison between the wake coefficients C_w of the numerical and experimental data that work as an indicator of the geometry matching of the wakes. The wake coefficient is evaluated as

where A_w is the wake area, i.e., the area of each surface normal to the flow direction along the wake iso-volume (m²), and A_p is the projected area of the snowflake in the flow direction (m²) (Nedic et al., 2013). The difference in the geometry of the experimental and numerical wake for steady falling cases [Figs. 4(a) and 4(b)] stays below 3%. Larger deviations can be found for particles with unsteady falling behavior, where the experimental wake is sometimes tilted in the y or z direction, and this is due to the dissimilarities in the wake shape between the numerical and experimental fields (Table I).This results in a less accurate matching of the wake iso-volumes [see Table I and Figs. 4(f) and 4(g)].

The estimate of the velocity uncertainty in the experimental flow field reconstruction is of around ±15%. The difference in the time- and space-averaged, non-dimensional velocity magnitude $U*$ between the measurements and the numerics is estimated using the relative percentage error

with ϵ_i being the error at each surface i normal to the flow direction along the wake iso-volume. By averaging the absolute values of ϵ_i, the averaged absolute error $ϵ¯$ is obtained. The data analysis and the post-processing are performed with ParaView 5.9 (Ayachit, 2015).

The spatial resolution of the computational field is higher than the one of the measurements; therefore, filtering operations are also carried out to facilitate the comparison of the unsteady cases (moderate/high Re), where the differences between the datasets are expected to be larger (Table I) (Buzzicotti et al., 2018; Corso et al., 2021). With this aim, two different uniform grids are employed with 2 and 4 mm resolution [Figs. 5(b) and 5(c), respectively]. These grid sizes are selected to be close to the experimental resolution (Sec. II A). Among the different methods available in ParaView (Ayachit, 2015), a Gaussian function is chosen for the interpolation over the coarser grid, varying the sampling radius according to the minimum grid size. The filtering of the numerical data and their comparison with experimental results are presented in Sec. III B 2.

In this section, the results of the comparison between numerical and experimental data are presented and discussed. The analysis is first carried out for snow particles with steady falling attitudes (low Re) for which good agreement between the experiments and the numerics is found (Sec. III A). Subsequently, the unsteady falling particles at moderate/high Reynolds numbers are considered. For these cases, two different analysis are performed: the first one involving the numerical data and the measurements with and without accounting for the particle rotation in the reconstruction (Sec. III B 1), and the second one including the experiments (without particle rotation), the numerical data, and the numerical data filtered on two larger grids (Sec. III B 2). This is done because, for unsteady falling snow particles, the differences between the numerical and experimental velocity field are much larger due to the high rotation rates displayed by the snowflakes. There is, therefore, the necessity to bring the two datasets closer for an easier comparison. Throughout this comparative study, the time-averaged, non-dimensional velocity field, together with the wake shape, is qualitatively compared. Then, for a more quantitative appraisal, the same quantities are spatially averaged over different sections in the streamwise direction along the wake iso-volume (⁠$U*$ and $gradU*$⁠, Sec. II D). The $U*$ relative error [Eq. (7)] is also reported to quantify the differences between the datasets. For comparison in Sec. III B 1, only the vertical component of $U*$⁠, i.e., $Ux*$⁠, is considered in the analysis.

Figures 6(a), 6(b), and 7(a)–7(c) illustrate the comparison for AG at $Re=$ 62. At this flow regime, AG has been observed to fall steadily. The time-averaged, non-dimensional velocity field shows similar values in both wakes [Figs. 6(a) and 6(b)] with slightly lower velocity close to the particle for the experimental case. In the latter, the wake shape is moderately tilted (Fig. 4), although the error in C_w is low (⁠$≈$3%, see Table I). The comparison of both $U*$ and $gradU*$ highlights good agreement between the two datasets with an averaged absolute error of 3.10% for $U*$ [see Figs. 7(a)–7(c)].

D1 at $Re=$ 206 also displayed a steady falling behavior during laboratory observations (Fig. 2). For this particle, the shape of the wake in the numerical field matches well the one from the experiments [see Table I and Fig. 6(c) and 6(d)]. The $U*$ experimental field moderately differs from the numerical case in the vicinity of the particle [Figs. 6(c) and 6(d)]. This difference is most likely caused by a lower density of the tracer particles close to the particle combined with its optical occlusion (i.e., the particle is not transparent and tracers are barely visible in some camera views near the snowflake analog, see Sec. II A), but this does not affect the quality of the comparison [Figs. 7(d)–7(f)] that results in an averaged absolute error of 1.79% for $U*$⁠.

At low Re, the investigated snow particles presented a steady falling behavior. This means that all the three angles (α, β, and γ) between the particle principal axes and its falling direction have values lower than 5° (see Fig. 2), i.e., the snowflake falls, on average, without a significant variation from its final orientation, as reconstructed from the experiments (Sec. II B). The 4D-PTV captures well the free fall of snow particles but presents lower resolution, especially in the vicinity of the particle, as compared to the numerical model. On the other hand, the DDES model of a fixed snow particle is a simplification of the real phenomenon, i.e., the particle is kept fixed and, thus, lacks the two-way coupling between the particle and the airflow. However, for particles that fall steadily, the numerical model succeeds in accurately representing the velocity field of a freely falling snow particle, despite the simplification. Therefore, it can be used as a robust tool to investigate snow particles with steady falling attitude.

At moderate/high Reynolds numbers (⁠$Re≳$400), all the examined snowflakes, except from D1, fall unsteadily. Figures 2(b)–2(d) show higher values of the averaged angles variation for CC and MR, which display chaotic falling behavior, while the angle variation for AG (periodic falling motion) stays below 20°. The large variation in the orientation of CC and MR is also visible in Figs. 8(a) and 8(b) that present the instantaneous positions manifested by both particles during their free fall. Here, we see that the final orientation, as reconstructed from the experiments, is displayed by CC and MR only at a couple of instants along the entire trajectory, while for the majority of their free fall time, other orientations are exhibited. The situation is much different for AG [Fig. 8(c)]. Due to its periodic falling attitude, it moderately changes its orientation, while its vertical trajectory is spiraling with respect to the falling direction. With this in mind, we analyzed the wake characteristics in cases of unsteady falling behavior performing two different comparisons: one that involves the numerical data and two different experimental datasets (with and without the particle rotation included in the velocity field reconstruction, as described in Sec. II A 1) presented in Sec. III B 1, and a second one that concerns the experimental data (without particle rotation) against the numerical data on two additional coarser grids, as compared to the computational one described in Sec. III B 2. The measurement reconstruction including the particle rotation and the filtering operations on the numerical data is performed to facilitate the comparison of unsteady falling cases for which larger differences are expected.

First, we analyze CC at $Re=$ 488, shown in Figs. 9(a)–9(d), 10(a), and 10(b). For CC, both its final and its extreme orientations (see Sec. II B and Tagliavini et al., 2021b) are used in the numerical model. The extreme orientation corresponds to the position for which the angle between the long axis of the particle and the vertical fall direction is observed to be the greatest. This particular orientation is included, because it is closer to the orientations that the particle typically displays during its fall [Fig. 8(a)]. If we look at the wake non-dimensional velocity field [Fig. 9(a) for the final orientation model, (b) for the extreme orientation one, (c) for the experiments without rotation, and (d) the experiments with rotation in the reconstruction, left-side image] and the wake shape [Figs. 9(a)–9(d), right-side image] significant differences are present. In particular, in the experimental case without rotation, the wake is slightly tilted due to the unsteady falling motion of the particle. This is why the C_w comparison results in high discrepancies of 18.65% and 23.29%, as compared to the final and extreme orientation models, respectively [Table I, Fig. 4(d)]. For this reason, the experimental dataset reconstructed considering the particle rotation is included to improve the comparison by bringing the numerical data closer to the experiments. In this case, the differences in the matching of the wake geometries (C_w) are reduced, even if only for the extreme orientation model [see Table I and Fig. 4(c)], from an average error of 23.29% to 7.22%. For a more quantitative comparison, the values of $Ux*$ are shown in Fig. 10(a) and its relative error in Fig. 10(b). Implementing the rotation in the experimental field reconstruction allows us to properly take into account all the positions that the particle exhibits, including those similar to the final orientation. As a consequence, the two datasets show lower differences that drop from 37.79% (without rotation) to 20.71% (with rotation) for comparison with the final orientation model and from 12.00% (without rotation) to 11.70% (with rotation) for the extreme orientation one (see Table I). The reason why the dissimilarities for the extreme orientation model are lower, as compared to the final orientation model, can be explained by the helical motion of CC (see Sec. II B) with positions frequently similar to the extreme orientation along its trajectory [Fig. 8(a)]. The average error for the extreme orientation model stays below the experimental uncertainty (±15%).

The same comparison is carried out for MR at $Re=$ 1110 (chaotic falling behavior). For this case, the final orientation model is compared with the experiments reconstructed with and without particle rotation [Figs. 9(e)–9(g)]. Despite the qualitative differences in the non-dimensional vertical velocity [Figs. 9(e)–9(g), left], the geometry matching of the wake is improved if the particle rotation is taken into account measurements reconstruction [from 11.66% to 6.67%, see Table I and Fig. 4(d)]. This additional manipulation of the experimental data allows for a decrease in the error of $Ux*$ from 36.41% to 22.68%, as shown in Figs. 10(c) and 10(d) and in Table I.

For AG at $Re=$ 1691 [periodic falling motion, Figs. 9(h)–9(j)], the experiments with the particle rotation show better agreement with the final orientation model, and the improvement in the error of $Ux*$ is in the same order of magnitude as for MR (approximately 13%). The measurements with rotation [Fig. 9(j), left] present two clearly separated recirculation zones in the wake, which slightly resemble the numerical velocity field, despite the more marked separation. The wake shape in this case is also more similar to the numerical one [Figs. 9(h)–9(j), right], although more tilted in the negative z direction. C_w in the experimental case with rotation is also in good agreement with the numerical one (Table I), but the matching of $Ux*$ [Figs. 10(e) and 10(f)] still results in an averaged error larger than the experimental uncertainty of the velocity field reconstruction (±15%), as for MR.

The additional manipulation of the experimental data to include the particle rotation in the reconstruction of the velocity field (Sec. II A 1) brings the measurements values closer to the numerical ones. This is achieved because the orientation of the particle in the experimental field is forced to match the numerical one. Nevertheless, in the majority of the investigated cases, the differences are larger than the experimental uncertainty due to the unsteady falling behavior of the particle that significantly modifies the average wake structure as compared to the fixed particle model.

MR at $Re=$ 1110 displayed chaotic falling behavior at high flow regimes. For this second comparison, the numerical data with the particle at its final orientation are filtered on two larger grids (namely, 2 and 4 mm, as described in Sec. II D) to investigate the influence of different spatial resolutions on the comparison. While for the DDES model [Fig. 11(a)] and the two filtered numerical datasets [Fig. 11(b) for 2 and (c) for 4 mm, respectively], the differences in the non-dimensional velocity field and in the wake shape are barely visible (a part from a slight change in the wake shape at the bottom). The experimental field (without rotation) substantially differs from the previous one [see Table I and Fig. 11(d)]. Nevertheless, the wake matching error remains around 12% [Fig. 4(d)]. The gap between the datasets is, however, evident in the comparison of $U*$ and $gradU*$ [Figs. 12(a) and 12(c), respectively]. Interestingly, the $U*$ relative errors [Fig. 12(b)] show an improvement of $≈$13% for the 4 mm grid compared to the original numerical data.

The same approach is followed for the analysis of AG at $Re=$ 1691 [Figs. 11(e)–11(h)] for which two filtered numerical datasets are also included (using a 2 and a 4 mm grid, as for MR). From Figs. 11(e)–11(h) (left), we can see that the strong velocity deficit in the wake is progressively smeared out in the filtered numerical fields as the grid size is increased and reaches values comparable to the experimental field (without rotation). The wake matching for this comparison [Fig. 4(e), Table I] presents an error of roughly 11%, and the wake shape is flattened as the filtering grid is increased. The comparison of $U*$ and $gradU*$ [Figs. 12(d) and 12(f), respectively] shows that the filtering of the numerical data brings the values closer to the ones of the experimental flow field, as it is also illustrated in Fig. 12(e), where the filtering on a 4 mm grid lowers the differences for $U*$ of $≈$ 12%.

The comparison at moderate/high Re allows a deeper insight in the capabilities and limits of both approaches (experiments and computational model). With regard to the experimental data, they correctly capture the wake characteristics of a freely falling snowflake analogs but have a limited spatial and temporal resolution that lead to a lower accuracy (±15% error in the velocity) and temporal statistics, as compared to the numerical model. On the other hand, more fine scale details of the velocity field and larger temporal statistics are present in the numerical model, albeit the use of DDES (hybrid URANS-LES), which carries intrinsic modeling inaccuracies (Tagliavini et al., 2021b). The numerical model becomes problematic when dealing with unsteady falling particles, because it requires a fixed orientation to be specified. For the majority of the tested geometries at moderate/high Re, their orientation considerably varies, and the falling unsteadiness grows with increasing Reynolds numbers (McCorquodale and Westbrook, 2020b). In this view, the numerical model is not able to accurately reproduce the average wake of freely falling snowflakes at moderate/high Re. Including the particle rotation in the experimental field reconstruction and filtering the numerical data on coarser grids partially improves the comparison. Accounting for the particle rotation in the measurements reconstruction brings a significant decrease in the wake geometry matching error (C_w), as compared to the filtering operation. Nevertheless, good agreement is not achieved for unsteady falling particle, because their movement alters significantly the wake characteristics. Measurements with larger time and space windows will allow for more orientations similar to the ones used in the model to be visited by the particle during its free fall (smaller differences in the wake characteristics) and, therefore, may provide a better agreement with the numerical data. In this view, we assess that using extreme orientations as fixed orientation in the model may better represent the average wake of variable orientations visited during free fall and, together with additional investigations over larger space/time windows, may provide a further understanding to what extent the fixed-particle model still gives reasonable results for unsteady falling snowflakes.

In this paper, velocity field data from a DDES model of a complex-shaped, fixed snowflake were compared with velocity measurements of freely falling, 3D-printed snow particle analogs of the same shape and at the same flow regime. To do so, the data were made non-dimensional to facilitate the comparison, and time- and space-averaged quantities were considered (⁠$U*$ and $gradU*$⁠, non-dimensional velocity magnitude, and non-dimensional velocity gradient magnitude). Before performing the comparison, the matching of the wake geometry of both datasets was assured. Two different particle falling attitudes were identified (steady and unsteady falling motion), and the comparison was carried out separately for each falling attitude. To investigate the sensitivity of the comparison, for particles with unsteady falling behavior, two analyses were performed: one involving the numerical dataset and two reconstructed measurement fields (with and without accounting for the particle rotation in the reconstruction), and the other one with the experiments (without rotation) and the numerical data filtered on two coarser grids of 2 and 4 mm resolution, respectively. Both approaches had the objective of bringing the two datasets (experiments and numerics) closer for a suitable comparison.

With respect to snow particles with a steady falling attitude (low Reynolds numbers, $Re≲$ 400), two cases were investigated: AG at $Re=$ 62 and D1 at $Re=$ 206. Both particles fell without significant variations from their final orientation, as reconstructed from the experiments. The differences between the experimental and numerical data are much smaller than the experimental uncertainty (±15%). Therefore, the fixed-particle, DDES model succeeded in representing the wake dynamics of steady falling cases, despite the simplified numerical setup.

For unsteady falling snowflakes, although the introduction of the experimental data with particle rotation implemented or the filtering of the numerical data, the discrepancies between the measurements and the numerics remained above the experimental uncertainty even if not dramatically. This second comparison highlighted the limitations of both approaches: the limited time and space window of the experimental setup that, nevertheless, captures well the free fall of snowflake analogs, and the lack of two-way coupling between the particle and the airflow in the DDES model (fixed particle). Measurements with larger time and space windows and additional numerical data on extreme orientations may provide better agreement with the numerical data. Therefore, using extreme orientations as fixed orientation in the model may better represent the average wake of visited orientations during free fall and, together with experimental data over larger space/time windows, may provide a further understanding to what extent the fixed-particle model still gives reasonable results for unsteady falling snowflakes.

The numerical simulations were performed using ETH-Euler cluster resources of Professor Roman Stocker's Group from the Department of Civil, Environmental, and Geomatic Engineering at ETH Zürich.

Giorgia Tagliavini would like to thank Pascal Corso (ARTORG-University of Bern) for his inputs and suggestions on the filtering operation of the numerical data as well as on the comparison between experimental and numerical data.

The authors have no conflicts to disclose.

G.T. setup the computational model and carried out the data analysis of the simulations results, the comparison with experimental data, and wrote the final draft of the paper. M.H. supervised the project, together with C.W.. C.W. and M.M. conceived and performed the experiments. M.H.K. analyzed the experimental data and performed the data reconstruction. M.M and M.H.K. wrote Sec. II A and Sec. II A.I., respectively. All authors discussed the results and contributed with valuable inputs to the final manuscript.

The data that support the findings of this study are available from the corresponding author upon reasonable request.