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

Experimental and numerical approaches have their own advantages and limitations, in particular when dealing with complex phenomena such as snow particles falling at moderate Reynolds number ( 𝑅𝑒 ). Time-resolved, three-dimensional Particle Tracking Velocimetry (4D-PTV) experiments of free-falling, three-dimensional (3D)-printed snowflakes analogs shed light on the elaborate falling dynamics of irregular snow particles, but present a lower resolution (tracer seeding density) and a limited field of view (domain size) to fully capture the wake flow. Delayed-Detached Eddy Simulations of fixed snow particles do not realistically represent all the physics of a falling ice particle, especially for cases with unsteady falling attitudes, but accurately predict the drag coefficient and capture wake characteristics for steadily falling snowflakes. In this work, we compare both approaches on time-and space-averaged flow quantities in the snowflake wake. Firstly, we cross-validate the two approaches for low 𝑅𝑒 cases where close agreement of the wake features is expected and, secondly, we assess how strongly the unsteady falling motion perturbs the average wake pattern as compared to a fixed particle at higher 𝑅𝑒 . For steadily falling snowflakes, the fixed-particle model can properly represent the wake flow with errors within the experimental uncertainty ( ± 15%). At moderate/high 𝑅𝑒 (unsteady falling motion) larger differences are present. Applying a co-moving


I. INTRODUCTION
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] ( ≫ 1, where  = (    )/, with  being the particle Reynolds number,   the snowflake terminal velocity (m/s),   being the maximum span of the particle normal to the fall direction (m), and  the kinematic viscosity of air (m 2 /s)).Hence, the interaction with the surrounding air becomes non-linear (large   implies large ) [Abraham 1970].Large snow particles (  ≫ 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 a 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 the one 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 snow flakes themselves.This was done by Mc-Corquodale and Westbrook 2020b to investigate the falling behavior and drag coefficient of freely falling, three-dimensional (3D)-printed snowflakes analogs in a mixture of water and glycerol.This study also aimed to improve former parametrized models for snowflake terminal velocity [Mitchell and Heymsfield 2005;Heymsfield and Westbrook 2010], which are generally constrained by their specific experimental or numerical set-up and confined to moderate Reynolds numbers ( ≲ 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 set-up accurately reproduces the numerical one and the object geometries are usually quite simple (spheres, cylinders, or airfoils), often at large Reynolds numbers ( ∼ 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 since 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. 2008and Auguste et al. 2013, while and Vincent et al. 2016and 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 andNedic et al. 2015 looked into the impact of 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 andToupoint 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 ( ≲ 10 3 ), due to limitations in the employed experimental techniques or in the computational capability of the chosen numerical approach.With regards 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 settling velocity up to Re∼ 10.With regards 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 particles drag coefficient and terminal velocity [Tagliavini et al. 2021a], at moderately high , 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 snowflakes analogs and a Delayed-Detached Eddy Simulation (DDES) that solves for the airflow around a fixed snow particle, previously validated for the 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.For the 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 .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 ), and another one at which the particles exhibit unsteady falling behavior (moderate/high ).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 Section II, the experimental set-up, the velocity field reconstruction, and the observations concerning the particle falling behavior are firstly described.Subsequently, the computational model is presented (Section II C), followed by the preparation of the data sets for the comparison and the analyzed quantities (Section II D).In Section III, the results from the comparison are provided and discussed.Initially, the steady falling cases are presented (Section III A), for which good agreement is found.Thereafter, the particle with unsteady falling attitudes are taken into account (Section III B).To facilitate the comparison, additional experimental data, whose reconstruction is done by including the particle rotation (Section 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 Section IV.

II. MATERIALS AND METHODS
This section includes the description of the experimental set-up, the velocity field reconstruction, the laboratory observations, and the main features of the DDES model.Conform to our previous works [Tagliavini et al. 2021a;Tagliavini et al. 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 data sets preparation is also described in Section II D.

A. Experimental set up
A schematic view of the experimental set-up is shown in Fig. 1.The experiments are conducted in a transparent acrylic box, the cross-section of which is a regular decagon with 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 to 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 ), 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 minutes.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 cm x 10 cm x 10 cm in size approximately 1.25 m below the surface of the fluid.Small, tracer particles (Fluorescent Red/Orange polyethlene Microspheres with diameter range 125 to 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 3 are used in water-glycerol mixtures of approximate densities of 999, 1088 and 1144 kg/m 3 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 Fig. 1 (b) and Fig. 1 (c)).The camera lenses (Zeiss Planar T* 50mm 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 LEDflashlight.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 x 1200 pixel resolution).Tracer particle trajectories are determined using LaVision's Shake-the-box algorithm.Tracer particle velocities are gridded to reconstruct velocity fields using a sub-volume size of 64 x 64 x 64 voxels (of order 0.6 x 0.6 x 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; 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; 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.

Reconstruction of the experimental flow field
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 Section 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 are time-averaged together.Following these steps, a second experimental data set 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).

B. Falling behavior of snow particles
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 ≲  ≲ 1700).
The plate-like crystal (D1) falls steadily during the laboratory observations throughout the entire  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  ≲ 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  (periodic motion).For what concerns CC, a steady falling attitude is present up to  ≲ 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 shows a clear periodicity, but as  increases the helical motions becomes chaotic.The same behavior is also seen by Kim et al. 2018 for rigidly linked disks.MR falls steadily up to  ≈ 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 Section II C. Fig. 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  = 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.

C. Computational model and simulations
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),  the pressure (Pa),  is the dynamic viscosity of the fluid (Pa• s),  is the fluid density (kg/m 3 ), 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 the unit vector pointing out of the particle surface A and  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 (Section II A) and is used to evaluate the uniform air speed set at the inlet boundary: where  is the desired particle Reynolds number (from experiments),   is the maximum dimension of the snowflake orthogonal to the flow direction (m),  ∞ is the inlet velocity (m/s), and  is the kinematic viscosity of the fluid (m 2 /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, 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 () at which the simulations are carried out are chosen such that they include either cases of steady (low ) and unsteady (moderate/high ) falling behavior.Namely, D1 is simulated at  = 206 (steady), AG at  = 62 (steady), 1691 (unsteady, periodic), MR at  = 1110 (unsteady, chaotic), and CC is evaluated at  = 488 (unsteady, chaotic).All the particles are simulated at their final orientation, except CC for which the additional extreme orientation is also considered (Section II B).For the summary of the different particle orientations and  investigated, see Tab.I.

D. Comparison between numerical and experimental data
Before evaluating the quantities employed in the comparison, the numerical velocity field is modified to a pseudo lab frame-of-reference by subtracting the inflow velocity  ∞ to the vertical component of the velocity.In this way, the experimental and numerical flow filed are much easier to compare.Moreover, the coordinates of both data sets are normalized using the particle maximum dimension   (m) orthogonal to the flow direction.As a consequence, the non-dimensional coordinate in the flow direction  * is evaluated using the following formula: where  in 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  0 is the inflow velocity equal to  ∞ (m/s) in the computational case and the snow particle terminal velocity   in the experimental one (m/s), respectively.Using the normalized coordinates and the time-averaged, non-dimensional velocity u * , the nondimensional velocity gradient  (u * ) is obtained, whose  2 norm (magnitude) is also calculated.
The snow particle wake is identified by an isovolume of the time-averaged, non-dimensional velocity magnitude ( * ).The iso-volume is delimited by threshold on  * set at 90% of its maximum value.Given that the field of view captured in the experiments extends only a few   from the particle, the comparison focuses on the quantities in the near wake ( * ≤ 3  ).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 the one of the laboratory field (Fig. 4 (f), (g)).Once the wake is outlined, the time-averaged, non-dimensional velocity magnitude ( * ) and the  2 norm of the time-averaged, nondimensional velocity gradient  ( * ) are spatially averaged over different surfaces orthogonal to the flow direction along the wake iso-volume ( * and  * , respectively, for which the time-averaging symbol is omitted to simplify the notation).For the 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, nondimensional velocity is taken into account ( *  , see Section III B 1).The distance between these surfaces is Δ * = 0.25  and the starting point is at  * = 0.5  for AG and D1, while for MR and CC is at  * = 0.75  due to the larger extension of those particles along the x axis (streamwise direction).Fig. 4 depicts the comparison between the wake coefficients   of the numerical and experimental data, that works as an indicator of the geometry matching of the wakes.The wake coefficient is evaluated as: where   is the wake area, i.e., the area of each surface normal to the flow direction along the wake isovolume (m 2 ), and   is the projected area of the snowflake in the flow direction (m 2 ) [Nedic et al. 2013].The difference in the geometry of the experimental and numerical wake, for steady falling cases (Fig. 4 (a) and (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 the experimental field (Tab.I).This results in a less accurate matching of the wake iso-volumes (see Tab.I and Fig. 4 (f) and (g)).
The estimate of the velocity uncertainty in the experimental flow field reconstruction is of around ±15%.The difference in the time-and spaceaveraged, non-dimensional velocity magnitude  * between the measurements and the numerics are estimated using the relative percentage error: with   the error at each surface  normal to the flow direction, along the wake iso-volume.By averaging the absolute values of   , the averaged absolute error  is obtained.The data analysis and the post-processing is 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 ), where the differences between the data sets are expected to be larger (Tab.I) [Buzzicotti et al. 2018;Corso et al. 2021].With this aim, two different uniform grids are employed with 2 mm and 4 mm resolution (Fig. 5 (b) and (c), respectively).These grid sizes are selected to be close to the experimental resolution (Section 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 Section III B 2.

III. RESULTS AND DISCUSSION
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 ), for which a good agreement between the experiments and the numerics is found (Section 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 (Section 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 (Section 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 data sets 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 ( * and  * , Section II D).The  * relative error (Eq.7) is also reported to quantify the differences between the data sets.For the comparison in Section III B 1, only the vertical component of  * , i.e.,  *  , is considered in the analysis.

A. Steady falling motion
Fig. 6 (a), (b) and Fig. 7 (a), (b), and (c) illustrates the comparison for AG at  = 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 (Fig. 6 (a) and (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   is low (≈ 3%, see Tab.I).The comparison of both  * and  * highlights a good agreement between the two data sets, with an averaged absolute error of 3.10% for  * (see Fig. 7 (a), (b), and (c)).
D1 at  = 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 Tab.I and Fig. 6 (c) and (d)).The  * experimental field moderately differs from the numerical case in the vicinity of the particle (Fig. 6 (c) and (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 Section II A), but this does not affect the quality of the comparison (Fig. 7 (d), (e), and (f)), that results in an averaged absolute error of 1.79% for  * .
At low , 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 (Section II B).The 4D-PTV captures well the free-fall of snow particles, but presents a 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.

B. Unsteady falling motion
At moderate/high Reynolds numbers ( ≳ 400), all the examined snowflakes, except from D1, fall un-steadily.Fig. 2 (b), (c), and (d) show higher values of the averaged angles variation for CC and MR, which display chaotic falling behavior, whilst 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 Fig. 8 (a) and (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 data sets (with and without the particle rotation included in the velocity field reconstruction, as described in Section II A 1) presented in Section 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 Section III B 2. The measurements reconstruction including the particle rotation and the filtering operations on the numerical data are performed to facilitate the comparison of unsteady falling cases, for which larger differences are expected.

Experimental data reconstructed accounting for the snow particle rotation
Firstly, we analyze CC at  = 488, shown in Fig. 9 (ab), (b), (c), and (d) and Fig. 10 (a) and (b).For CC, both its final and its extreme orientations (see Section 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 (Fig. 9 (a), (b), (c) and (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   comparison results in high discrepancies of 18.65% and 23.29%, as compared to the final and extreme orientation model, respectively (Tab.I, Fig. 4 (d)).For this reason, the experimental data set 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 (  ) are reduced, even if only for the extreme orientation model (see Tab.I and Fig. 4 (c)), from an average error of 23.29% to 7.22%.For a more quantitative comparison, the values of  *  are shown in Fig. 10 (a) and its relative error in Fig. 10 (b).Implementing the rotation in the experimental field reconstruction allows to properly take into account all the positions that the particle exhibits, including those similar to the final orientation.As a consequence, the two data sets show lower differences that drop from 37.79% (without rotation) to 20.71% (with rotation) for the comparison with the final orientation model and from 12.00% (without rotation) to 11.70% (with rotation) for the extreme orientation one (see Tab. 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 Section 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  = 1110 (chaotic falling behavior).For this case, the final orientation model is compared with the experiments reconstructed with and without particle rotation (Fig. 9 (e), (f), and (g), respectively).Despite the qualitative differences in the non-dimensional vertical velocity (Fig. 9 (e), (f), and (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 Tab.I and Fig. 4 (d)).This additional manipulation of the experimental data allows for a decrease in the error of  *  from 36.41% to 22.68%, as shown in Fig. 10 (c) and (d) and in Tab.I.
For AG at  = 1691 (periodic falling motion, Fig. 9 (h), (i), and (j)), the experiments with the particle rotation show a better agreement with the final orientation model, and the improvement in the error of  *  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 (Fig. 9 (h), (i), and (j), right), although more tilted in the negative z direction.  in the experimental case with rotation is also in good agreement with the numerical one (Tab.I), but the matching of  *  (Fig. 10 (e) and (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 (Section 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.

Numerical data filtered on coarser grids
MR at  = 1110 displayed chaotic falling behavior at high flow regimes.For this second comparison, the numerical data with the particle at its final orientation is filtered on two larger grids (namely, 2 mm and 4 mm, as described in Section II D), to investigate the influence of a different spatial resolution on the comparison.While for the DDES model (Fig. 11 (a)) and the two filtered numerical data sets (Fig. 11 (b) for 2 mm 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 ones (see Tab.I and Fig. 11 (d)).Nevertheless, the wake matching error remains around 12% (Fig. 4  (d)).The gap between the data sets is however evident in the comparison of  * and  * (Fig. 12 (a) and (c), respectively).Interestingly, the  * 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  = 1691 (Fig. 11 (e), (f), (g), and (h)), for which two filtered numerical data sets are also included (using a 2 mm and a 4 mm grid, as for MR).From Fig. 11 (e), (f), (g), and (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), Tab.I) presents an error of roughly 11%, and the wake shape is flattened as the filtering grid is increased.The comparison of  * and  * (Fig. 12 (d) and (f), respectively) shows that the filtering of the numer-ical data brings the values closer to the ones of the experimental flow filed, as it is also illustrated in Fig. 12 (e), where the filtering on a 4 mm grid lowers the differences for  * of ≈ 12%.
The comparison at moderate/high  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 , 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 .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 bring a significant decrease in the wake geometry matching error (  ), as compared to the filtering operation.Nevertheless, a 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.

IV. CONCLUSIONS
In this paper, velocity field data from a DDES model of a complex-shaped, fixed snowflake were compared with velocity measurements of freely falling, 3Dprinted 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 ( * and  * , non-dimensional velocity magnitude and non-dimensional velocity gradient magnitude, respectively).Before performing the comparison, the matching of the wake geometry of both data sets 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 data set and two reconstructed the 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 mm and 4 mm resolution, respectively.Both the approaches had the objective of bringing the two data sets (experiments and numerics) closer for a suitable comparison.
With respect to snow particles with a steady falling attitude (low Reynolds numbers,  ≲ 400), two cases were investigated: AG at  = 62 and D1 at  = 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 set-up.
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 set-up, 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 a 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.together with the averaged absolute error (Eq.7) of the wake coefficient and of the time-and space-averaged, non-dimensional velocity magnitude are listed.The error for   and  * are separated for the two different experimental data sets: with and without the particle rotation implemented in the reconstruction (Section II A 1). Furthermore, the results from the filtering operation (coarser grids) of the numerical data are also presented.In the comparison of the experimental data with the particle rotation included, the errors are evaluated with respect to the vertical component of the non-dimensional velocity  *  (Section II D).(b) and (e), respectively) is shown with dotted lines.Since the coordinate system is made non-dimensional using the particle maximum diameter   , the non-dimensional coordinate  * represents the distance from the particle in terms of its maximum extension normal to the flow direction.
FIG. 1: Sketches of the experimental set-up showing its key components, including the position of the cameras, relative to the LED illumination, and the region of interest that they record.(a) Side and (b) plan view of the set-up.Also shown are the coordinate directions (x, y, z), and (c) the real experimental set-up.

FIG. 4 :
FIG. 4:Wake coefficient   for each snow particle ((a-b) steady falling cases and (c-e) unsteady falling ones).  works as an indicator of the wake geometry matching.The averaged absolute error of   for each case is summarized in Tab.I.The matching of the wakes is shown in (f) and (g) for a steady falling particle (D1 at  = 206) and for an unsteady falling one (AG at  = 1691), respectively.The experimental wake is the solid one, while the numerical one is transparent.

FIG. 5 :
FIG. 5: Coarse grids employed in the filtering operation of the numerical data.(a) Original computational grid, (b) 2 mm grid, and (c) 4 mm grid for particle MR, as example.

FIG. 6 :
FIG. 6: Results of the comparison between experimental and numerical data for snow particles with steady falling behavior.The pictures illustrate the contour plots of the time-averaged, non-dimensional velocity magnitude  * and the wake iso-volumes for AG at  = 62 ((a) and (b)) and D1 at  = 206 ((c) and (d)).

FIG. 7 :
FIG. 7:Results of the comparison between experimental and numerical data for snow particles with steady falling behavior.The plots report the comparison of the time-and space-averaged, non-dimensional velocity magnitude  * , its relative error, and the time-and space-averaged, non-dimensional velocity gradient magnitude  * for AG at  = 62 ((a), (b), and (c)) and D1 at  = 206 ((d), (e), and (f)).The experimental uncertainty on the velocity field reconstruction (±15%) (plot (b) and (e), respectively) is shown with dotted lines.Since the coordinate system is made non-dimensional using the particle maximum diameter   , the non-dimensional coordinate  * represents the distance from the particle in terms of its maximum extension normal to the flow direction.

FIG. 8 :FIG. 12 :
FIG. 8: Falling orientations reconstruction for (a) CC at  = 488, (b) MR at  = 1110, and (c) Ag at  = 1691.The  * = /  axis is parallel to the falling direction in this case ( * is the falling direction in the comparison).(a) and (b)exhibit chaotic falling behavior, while (c) falls with a periodic attitude (see Tab.I and Fig.2).

TABLE I :
Summary of the cases studied in this work.The Reynolds number and the falling behavior of each particle,