Numerical analysis of the wake of complex-shaped snow particles at moderate Reynolds number

moderate Reynolds number Giorgia Tagliavini,1, 2 Mark McCorquodale,3, 4 Chris Westbrook,3 and Markus Holzner2, 5 1)Institute of Environmental Engineering, Department of Civil, Environmental and Geomatic Engineering, ETH Zürich, Zürich, Switzerland 2)Swiss Federal Institute for Forest, Snow and Landscape Research, WSL, Birmensdorf, Switzerland 3)Department of Meteorology, University of Reading, Reading, UK 4)Department of Civil Engineering, University of Nottingham, UK 5)Swiss Federal Institute for Aquatic Science and Technology, Eawag, Dübendorf, Switzerland


I. INTRODUCTION
An accurate prediction of snow precipitation is crucial to constrain climate models parametrization. The orientation of ice particles is critical to our understanding of polarimetric radar measurements (for weather prediction, monitoring of clouds and precipitation, and for investigating ice properties) [Lemke and Quante 1999;Geier and Arienti 2014], while their wake structure influences particles growth rates [Libbrecht 2009]. Snowflakes exhibit a wide range of sizes and shapes in nature [Kikuchi et al. 2013], from tiny snow crystal with a maximum dimension of few micrometers to large aggregates of snow crystals that reaches several centimeters in size. For large snow particles (D max > 100 µm), the C D − Re relation becomes non-linear, resulting from a complex interaction between the particle and the surrounding air (Re 1, where Re = u t D max /ν is the particle Reynolds number, D max is the particle maximum extension orthogonal to the flow direction [m], u t the snowflake terminal velocity [m/s] and ν the kinematic viscosity of air [m 2 /s]), and one cannot rely on Stokesian dynamics (valid for Re < 1) [Happel and Brenner 1983;Westbrook 2008;Zeugin et al. 2020]. Further complexity is added to snow particle falling motion because of their irregular shapes, which produce unsteady wake flow, and intricate falling trajectories, such as tumbling, oscillations and fluttering [Gunn and Marshall 1957;Nemes et al. 2017;McCorquodale and Westbrook 2021b;Zeugin et al. 2020].
To understand the complex free-falling behavior of snowflakes, the mutual influence of wake flow, particles shape, and drag needs to be investigated. Many numerical and experimental studies have investigated aerodynamic coefficients of complex-shaped particles [Magono and Nakamura 1965;Dietrich 1982;Mitchell and Heymsfield 2005;Westbrook and Sephton 2017]. In the atmospheric sciences community, experimental works have presented empirical relations to predict the drag coefficient and the terminal velocity of snow particles, both from field and laboratory measurement. Recently, McCorquodale and Westbrook 2021b performed laboratory measurement with 3D-printed snowflakes free-falling vertically in waterglycerol mixtures. They investigated different particle geometries (from plate-like crystals to aggregates) and their aerodynamic response, together with their falling attitudes. The study also discussed and improved the parametrization proposed by Heymsfield and Westbrook 2010 to predict snow particles terminal velocity, which includes parameters to quantify the particle shape (i.e., area ratio) and building on the work of Khvorostyanov andCurry 2005 andMitchell andHeymsfield 2005. The accuracy of these parametrized models are generally dependent on the measurements resolution and on the chosen parameters and empirical constants, often leading to large error in the prediction of the terminal velocity, as shown by McCorquodale and Westbrook 2021b. Moreover, for such models, the performance at high Reynolds numbers (Re 100) is generally poor due to the difficulty of introducing a suitable turbulence correction. On the numerical side, the research of Zastawny et al. 2012, Ouchene et al. 2016, and Sanjeevi et al. 2018 illustrated the aerodynamic coefficient prediction for non-spherical particles, such as disks, ellipsoids and fibers, starting from the work of Leith 1987. Tran-Cong et al. 2004, Binder et al. 2006, and Dioguardi et al. 2018 introduced new relations for much more complex geometries, based on sphericity or specific shape factors. The proposed relations, although well performing within the tested parameter range, are built on the specific models or experimental set-ups and are therefore difficult to generalize. Furthermore, these works focused on how the drag coefficient changes according to variations of the Reynolds number or particle orientations, without looking at particle's wake influence on drag. Various numerical models of free-falling non-spherical particles, such as disks and plate-like geometries [Auguste et al. 2013;Yang et al. 2014;Cheng et al. 2015;Nettesheim and Wang 2018], or columnar crystals and hexagonal plates [Hashino et al. 2016], were also employed to predict the drag coefficient and to investigate the particle falling attitude. Nonetheless, they focus mainly on investigating geometric shapes representative of "pristine" ice crystals, with falling trajectories limited to Re 1000 and a simpler wake than the one of irregular particles, such snowflakes.
Very few studies have been published so far on more irregularly shaped snow crystals or aggregates and their interaction with the environmental airflow. Continuing the research of Garrett et al. 2015 on the influence of turbulence on snow particle fall speed and orientations, the field observations performed by Li et al. 2021 showed that when snow vertical velocity is large compared to the characteristic turbulent velocity, strong accelerations are found among snow particles. Furthermore, preferential sweeping caused by turbulence enhances settling by bringing snowflakes in downward regions of the airflow. The recent numerical work of Zeugin et al. 2020 investigated realistic snow particles in the Stokes regime and proposed sphericity-based relations to estimate aerodynamic coefficients and the settling velocity of such particles up to Re < 10. Tagliavini et al. 2021 proposed a new method to evaluate the time-averaged drag coefficient of complex-shaped snowflakes based upon a Delayed-Detached Eddy Simulation (DDES) model of flow around a fixed particle and its inertia tensor for Reynolds number up to ≈ 4000. However, to the authors' knowledge, no studies so far have investigated in detail the wake flow of such complex geometries.
The wake that forms behind a particle influences both the particle drag and its falling motion. The most closely related studies have investigated the wake structure that forms behind simple, idealized particle geometries. In these studies, the wake that forms behind a particle influences both the particle drag and its falling motion. Uhlmann and Dusek 2014 numerically investigated a falling sphere in a viscous fluid at rest and identified different patterns of motion (from steady to chaotic) according to different flow regimes. They classified the sphere falling behavior by analyzing diverse structures in the wake. Liu et al. 2021 examined the settling behavior of dual spheres both experimentally and numerically. Their study reveals that particles terminal velocity and final separation distance is strongly influenced by the Reynolds number and their initial distance, whilst the repulsion between particles depends on the intensity and type of vortices formed by both particles while falling. Hashino et al. 2016 andToupoint et al. 2019 found a strong dependence of their motion on the aspect ratio of columnar crystals and disks, respectively, and discussed in detail the varying forces and torques in relation to the falling behavior. Recent experimental studies by Nedic et al. 2013 andNedic et al. 2015 investigated the drag coefficient of fractal, plate-like geometries and the influence of the fractal dimension on wake size and vortex shedding. The increase of the fractal dimension, i.e. longer particle perimeter, caused an increment in the drag coefficient (up to 7%) and in the shedding intensity. Cummins et al. 2018 looked at the separated vortex ring in the wake of a dandelion pappus to shed light on its unconditionally stable falling motion, showing that the pappus filaments produce separated vortex rings that cause the particle to fall steadily with higher drag coefficient when compared to a disk. Despite the use of non-spherical particles, limitations in the Reynolds number are once more present in these works and the wake properties of realistic snow particles across flow regimes are unknown.
Beyond fundamental interest, measurements of the wake are important as these can be used to accurately estimate drag force on a particle via analysis of the momentum flux [Durbin and Medic 2007]. The evaluation of the drag force through the integration of the wake momentum flux is commonly employed in experimental studies and computational models of bluff body flow, in cases when it is difficult to accurately estimate the drag force by direct measurements, and fits in the category of wake integration techniques [Shin et al. 2010;Terra et al. 2017]. This technique is sensitive to different factors, including grid density, type of grid, numerical smoothing, and the definition of the control volume [Giles and Cummings 1999;Zhu et al. 2007;Shirazi and Manshadi 2020]. Nevertheless, if these parameters are correctly set, it delivers a satisfactory prediction of the drag force and helps comprehending how different terms of the momentum deficit affect the total drag. To the best of our knowledge, only few works evaluated separately the momentum deficit terms to investigate how they individually contribute to the drag. In Terra et al. 2017, the momentum flux was decomposed, whilst experimentally measuring the drag of a sphere moving in stagnant air at Re = 10000. It was found that the velocity fluctuation term accounts for roughly the 10% of the momentum exchange in the wake, whereas the pressure term becomes negligible.
The objective of this paper is to gain a better understanding on the influence of wake flow and particle shape on drag and their potential impact on the falling motion of complex-shaped snow particles. With this aim, a DDES model, validated for the drag coefficient prediction with experimental data of 3D-printed falling snowflakes [McCorquodale and Westbrook 2021b;Tagliavini et al. 2021], is used. It solves for the airflow past a fixed, irregular snowflakes at both low (Re = 50, 75, 100) and moderate/high Reynolds numbers (Re = 500, 1000, 1500). These flow regimes are chosen such that they fall within the range of Reynolds numbers exhibited by falling snow particle in nature [Heymsfield and Westbrook 2010]. As in Tagliavini et al. 2021, the model is informed with laboratory data of free-falling 3D-printed snow particles (McCorquodale and Westbrook 2021b). The same experimental data are used in the discussion to relate the model results with empirical observations of falling behaviors. Firstly, the wake topology is investigated to identify peculiar flow features that may have an influence on the drag coefficient. Secondly, the momentum flux in the wake of each snowflake is analyzed and decomposed into its different terms (mean velocity, pressure, and fluctuating velocity) to assess their influence on the total drag at varying flow regimes. In the end, a novel approach that combines a DDES model with experiments is exploited to analyze the wake topology and the momentum flux decomposition to advance in the understanding of complexshaped particle drag coefficient and falling attitudes.

II. MATERIALS AND METHODS
In this section, the experimental set-up and the laboratory observations are briefly described. Then, the DDES model and its main features are introduced. Following our previous work [Tagliavini et al. 2021], the numerical model is informed, and its results then interpreted, by comparison with data from experiments [McCorquodale and Westbrook 2021b].

A. Experimental method
Laboratory experiments, exploiting the principle of dynamic similarity, were conducted on 3D-printed snowflake analogs dropped in a tank containing viscous fluid at rest (uniform mixtures of water and glycerol, with the volume fraction of glycerol between 0 and approximately 50%). The trajectory of falling analogs was recorded using a system of three synchronized cameras positioned to have orthogonal views of the tank. An algorithm (the Trajectory Reconstruction Algorithm implemented through Image anaLysis, TRAIL) allowed the reconstruction of the timeresolved trajectory and orientation of falling analogs over the recorded region (approximately 200 × 200 × 200 mm in size) [McCorquodale and Westbrook 2021a]. With this approach, it was possible to determine the precise orientation adopted by each particle in the free-fall and to represent these positions digitally with a mesh of triangular facets that correspond to the particle surface. Reynolds numbers and drag coefficients were also evaluated (within an accuracy of approximately 5% and 7%, respectively). For more details on the algorithm and on the accuracy of the orientations reconstruction, see McCorquodale and Westbrook 2021a.

Observations from experimental data
Snowflakes exhibit many shapes and sizes in nature [Kikuchi et al. 2013]. The snow particle shapes incorporated in this work are shown in Fig. 1 and are chosen such that they span across some key types of natural ice particles, namely, monocrystals (D1 -D1007), multi-habit crystals (CC -CC20Hex4), polycrystals (MR -MR172) and aggregates (AG -AgSt100) [McCorquodale and Westbrook 2021b]. The names in brackets refer to the extensive data reported in the supplemental material of McCorquodale and Westbrook 2021b. They do not refer to any physical property of the particle and are simply used to facilitate crossreferencing with data from previous works. Because of snowflake diverse geometries, a distinct falling be-havior is observed for each snow particle and will be briefly described.
During experiments, described in Section II A, the dendrite plate crystal D1 fell steadily for the tested Reynolds number range (10 Re 1500), and adopted an orientation with its largest projected area orthogonal to the falling direction (i.e., with the maximum principal moment of inertia aligned parallel to the fall direction), see Fig. 1. A steady falling behavior was also seen for AG at Re 1000, for MR at Re 150, and for CC only at Re 70. AG displayed a spiraling trajectory, rotating around a (vertical) axis aligned with the fall direction, as Re increases. CC also fell with a periodic trajectory in the range 70 Re 400 [McCorquodale and Westbrook 2021a], similar to a pair of rigidly connected disks [Kuroda et al. 2007]. For Re 400, CC and MR exhibited chaotic motion. The falling behavior of the different geometries is summarized in Fig. 2. Here, the variation (e.g., ∆α = α max − α min ) of the angles between the principal axes and the fall direction for each geometry is illustrated as a function of Re (α: angle between the fall direction and the largest principal axis, β: angle between the fall direction and the intermediate principal axis, and γ: angle between the fall direction and the smallest principal axis). Large variations (approximately greater than 5 • ) indicate unsteady falling behavior. More details can be found in McCorquodale and Westbrook 2021b.
The orientations shown in Fig. 1 represent the final orientations employed in the numerical model. These orientations are the ones observed during experiments when the particle was falling steadily (at low Reynolds number). They can also be considered as a good approximation for cases in which particles were observed to fall unsteadily during the experiments (i.e., at high Reynolds number), since the low Reynolds number orientation lies within the range of the orientations exhibited by the particles while falling unsteadily. Tagliavini et al. 2021 also showed that, at high Re, for some particles that undergo chaotic falling behavior, the drag coefficient in this assumed final orientation can strongly differ to the actual value. For these cases, they proposed a suitable average of two orientations with the largest and smallest projected area (according to the inertia tensor eigenvalues) that provided a good approximation of the chaotic behavior. Since particles CC and MR displayed chaotic motion already at moderate Reynolds numbers (Re 400), for this study, additional orientations observed during experiments were reconstructed for use within our numerical model to investigate more in depth the range of variability of the wake flow in cases of unsteady falling behavior. In particular, one additional orientation is chosen for CC at Re = 500, 1000, 1500, corresponding to the position for which the angle be-tween the long axis of the particle (i.e., aligned with the column) and the vertical fall direction is greatest. For MR (Re = 1500), three additional orientations are considered, which correspond to instantaneous orientations observed at different time instants as particles fell across the measurement region (namely time-steps t = 13, 94, 129). These orientations are limiting cases in terms of the angles between the principal axes of the particles and the vertical fall direction and therefore will be called " extreme" orientations throughout the paper (see Fig. 1).

B. Numerical model
The computational model presented in this study is based upon Delayed-Detached Eddy Simulation (DDES -hybrid LES-RANS approach) [Spalart et al. 2006], solves for the airflow around a realistic snow particle, and has been already validated for the estimation of the drag coefficient [Tagliavini et al. 2021]. The particle within the computational domain is scaled with the same volume-equivalent sphere of d eq = 0.01 m. In this way, the computational domain size remains constant throughout the runs. The DDES model employs the Spalart-Allmaras turbulence closure model forν (modified eddy viscosity) for the RANS calculation (described in Spalart and Allmaras 1992). Numerical simulations have been performed for different flow regimes from low (Re = 50, 75, 100) to moderate/high Reynolds numbers (Re = 500, 1000, 1500). The uniform inflow velocity u ∞ [m/s] is calculated as: where Re is the desired Reynolds number, ν is the fluid kinematic viscosity [m 2 /s] and D max is the maximum dimension of the snowflake (orthogonal to the flow direction) [m]. The characteristic length D max in Eq. (1) is chosen as such because it is an easily measurable geometrical feature for ice particles. Thus, our choice is motivated by this practical reason, rather than a physical one. For additional information regarding the size of the computational model, its grid features and its initial and boundary conditions, see

Numerical simulations
The DDES model is built with Open source Field Operation and Manipulation (OpenFOAM 4.1) C++ software based on the finite volume method [Open-FOAM 2017]. It solves transient Navier-Stokes equations for the fluid domain: in which u is the flow velocity [m/s], ρ is the fluid density in [kg/m 3 ], p the pressure in [Pa], µ is the dynamic viscosity of the fluid in [Pa· s] and f are any external forces per unit mass [N/kg]. The interaction between the airflow and the snowflake is expressed through Newton's second law of motion: where n is the normal and tangential unit vector on the particle surface A and τ represents the viscous stresses [Pa]. The transient solver, the discretization schemes, and the number of iterations are the same of Tagliavini et al. 2021.
The time-averaged drag force is evaluated from the integration of the momentum deficit in the wake. The wake volume is isolated using a threshold of 98% of the total time-averaged velocity deficit and the momentum flux is calculated at a distance of x ≈ 4D eq from the snow particle. These criteria are followed according to the guidelines proposed by Zhu et al. 2007 for complex-shaped bodies. Furthermore, they assure the minimum percentage error betweenthe time-averaged drag coefficient computed from the simulations and the drag coefficient evaluated from the momentum flux integration (C D,sim −C D,mom )/C D,sim · 100 (for more details, see Fig. 4). The formula for the momentum deficit is: in which F D is the time-averaged drag force [N],ū is the mean velocity [m/s], u is the fluctuating velocity [m/s], p 0 is the ambient pressure [Pa], andp is the time-averaged pressure [Pa]. The viscous stresses in the momentum flux, evaluated at S wake account for (orthogonal to the flow direction), account for 0.5% or less of the total drag force and are thus not shown. From the obtained drag force, the time-averaged drag coefficient is determined, together with the different drag coefficient components derived from each term of the momentum deficit: where C D,tot is the total drag coefficient, while C D,ū , C D,ū 2 , and C D,∆p are the contributions to the total drag of the first, second, and third term of the right hand side of Eq. (4), respectively, and A p is the particle projected area [m 2 ]. The contribution of each momentum balance term is then expressed as a percentage of the total drag coefficient (e.g., To appraise the influence of the particle shape on the wake flow and consequently on the drag coefficient, two geometrical features are considered in Section IV to characterize the snowflake geometries: the porosity and Corey's shape factor [Corey 2019]. Porosity is defined as the ratio between the empty area to the planar projected area of the enclosing disk: in which A R = A p /A disk is the area ratio of the particle and A disk the area of the enclosing disk [m 2 ] (McCorquodale and Westbrook 2021b). Porosity goes from a maximum value of 1 (void) to a minimum of 0 (solid disk). For the case of D1 and AG, porosity originates from two factors. First, it is due to voids or pores within individual dendrites of the geometry. Second, it is also due to the fact that the shape of the enclosing envelop of the geometry deviates from a circular shape. For the case of CC and MR there are no dendrite features and the porosity is entirely due to the second factor. Corey's shape factor is defined as: where S, L, and I are the shortest, longest, and intermediate particle dimension, respectively. Corey's shape factor is equal to 1 for a sphere (or a cube) and tends to 0 for a thin disk. It does not account for the particle roundness (cube and sphere have the same value), but delivers information on the particle flatness, which is important for the wake formation, since it influences the separation point of the shear layers [Corey 2019]. The values of and CSF for each snowflake are reported in Tab. I and Tab. II and note that they vary for the extreme orientations.

III. RESULTS
The objective of this work is to better understand the mutual influence of shape, wake flow, drag coefficient and falling motion of snowflakes at varying flow regimes (50 ≤ Re ≤ 1500). To this end, we first investigate qualitative aspects of the wake behind the snowflakes across flow regimes and subsequently we address the wake contributions to the drag coefficient by means of Eq (4).
In Fig. 4, the drag coefficient for each snow particle at different flow regime (Re) is reported. At low Reynolds numbers (Re = 50, 75, 100), the dendrite plate (D1) shows the highest drag coefficient, while MR has the lowest values. As Re increases, the difference between the drag of each particle becomes smaller up to Re = 1500, where C D ≈ 1 for all the particles. In the same plot, the drag coefficients of the extreme orientations of particle MR and CC are also reported with empty markers. For CC, the drag is slightly higher for the extreme orientations in comparison to the final one, in particular at Re = 500. MR displays the lowest drag coefficient for the extreme position at t = 94 (C D = 0.59), followed by t = 13, and t = 129, while the highest value is recorded for the final orientation (C D = 1.32).
The wake topology of the different snowflakes shapes is presented in Section III A, while the results from the analysis of the wake momentum flux are described in Section III B. Finally, Section III C concerns the results of wake topology and momentum flux decomposition of the extreme orientation of particle CC and MR. Fig. 5, 6, and 7 show the wake flow as line integral convolution (LIC) vector field visualization, which produces streamline patterns tangential to the velocity field [Loring et al. 2014]. The time-averaged field chosen for the color map is the non-dimensional velocity magnitude U * = U/U ∞ (where U is the velocity vectors [m/s]). The plane in which the wake of each particle has its maximum lateral extension is selected in the figures. For conciseness, only Re = 100, 500, 1500 are described in detail here, since they display distinct qualitative differences.

A. Wake topology of different snow particles at increasing flow regimes
At Re = 100 (Fig. 5), D1 has the widest recirculation region with the highest velocity deficit. A large, symmetric vortex ring appears separated from the particle and the high velocity deficit extends further downstream. Similar patterns can be found for Re = 50, 75 with the recirculation region stretching in the streamwise direction as Reynolds number increases. MR and AG create asymmetric wakes with recirculation areas attached to the particles. For AG, the vortex region is present only at one side of the wake. The asymmetry also reflects in the top part of the wake, which displays a higher velocity deficit compared to the bottom part. CC displays two vortical regions: one between the hexagonal plates and the other one right after the particle. At low flow regimes, all the particles present stable falling behavior, except CC (periodic).
With regards to the flow field at Re = 500 (Fig.  6), the vortex ring of D1 is smaller and closer to the particle compared to Re = 100, while maintaining its symmetry and increasing its velocity deficit. The latter is the highest among all particles also at this flow regime and expands much further downstream. AG has a higher velocity deficit compared to MR and both present an asymmetric wake, tilted at one side (top part). CC (chaotic motion) presents symmetric recirculation zones between the hexagonal plates, while the previous vortex rings see a decrease in size and and increase in the velocity. U * extends significantly further downstream also for CC, but much lower values if compared with D1.
At Re = 1500 (Fig. 7), all snow particles have small and asymmetric vortex rings mostly attached to the particle. For AG and D1, the velocity deficit in the wake is higher and extends much further downstream (D1 presents the lowest value of U * in the far wake) in comparison to CC and MR (Fig. 7). Both wakes of AG and MR are tilted at one side, while CC displays a small velocity deficit at wake sides in the far wake, while at the center it is close to 0.
To better illustrate the unsteady behavior of the wake flow and its features, Fig. 8 and Fig. 9 show the 3D (left) and 2D (right) instantaneous wake flow field. The 3D structures are extracted with the Q-criterion (Q > 1 for Re = 500 and Q > 50 for Re = 1500) and colored according to the non-dimensional, streamwise component of U * (u * x ), while a 2D side view of the wake is depicted using the non-dimensional vorticity (ω * = (D max ω)/U ∞ ).
At Re = 500, D1 (Fig. 8a) presents a steady wake with a central recirculation bubble. Vortex threads (with lower velocity deficit) extend in the near wake in correspondence with the branches of the dendrite. AG (Fig. 8b) displays an unsteady wake flow with moderate vortex shedding. The wake structures are mainly hairpin vortices that origin from the particle sides, are shed in pairs several D max downstream, and eventually deform, interact, and decay in the far wake. MR (Fig. 8c) has a stretched, singular hairpin votex in the near wake that breaks down in smaller vortices after few D max downstream, producing the vertically extended meandering in the far wake (Fig. 8c, left). Even though the particle is symmetric and aligned with flow direction, the wake of CC is now slightly oscillating (Fig. 8d, right). These small oscillations can be better seen in the 3D representation of Fig. 8d (left), which displays an asymmetric vortex pattern. Fig. 9 illustrates the wake flows at Re = 1500. All the particles exhibit unsteady wakes and the instantaneous fields capture the diverse wake dynamics for different geometries. D1 is the only shape that does not exhibit a distinct vortex street at this high Re. Its wake is compact with small hairpin vortices that originate from the snowflake branches and become more fragmented far away from the particle (Fig. 9a). AG presents the highest values of vorticity downstream, among the particles that display a meandering wake (AG, MR, and CC). Horseshoe vortices form at the side of the particle and move downstream, before dissipating and fragmenting in the far wake (Fig. 9b). MR displays two stretched vortices in the near wake that break down a few D max downstream in a vortex street (Fig. 9c). CC exhibits an elongated thread closer to the particle that breaks down into a large vortex street with smaller structures (Fig. 9d). Beside D1, we note that all the other snowflakes were observed to fall unsteadily in experiments at Re = 1500 (CC, MR: chaotic, AG: periodic). Fig. 10 depicts the non-dimensional drag force fluctuations and the wake turbulent intensity as a function of the Reynolds number for each particle. At Re ≤ 100, both drag force and velocity fluctuations (T I) stay below 0.5%. As Re increases, a sharp increment can be seen for the two quantities, with the highest values for MR, followed by CC and AG. D1 exhibits the lowest values for both std(F D )/F D and T I.

B. Evaluation of the drag coefficient from the integration of the momentum deficit
Next, the influence of different terms in Eq. (4) on the total drag is assessed. To do so, the time-averaged drag coefficient C D,tot is split into different components according to Eq. (6) -(8). Each contribution is expressed as a percentage over the total drag (e.g., C D,ū /C D,tot · 100) and listed in Tab. I. The comparison between different shapes is depicted in Fig. 11 as bar plots, to stress the relative contribution of the different momentum flux terms. For all the particles, the viscous term turns out to be negligible and thus it is not shown, while the mean velocity term provides the overall largest contribution in comparison with the other terms. D1 and AG have the highest values for the mean velocity term contribution (C D,ū ) (above 83% throughout the tested Re), whilst the pressure term (C D,∆p ) stays below 7%. D1 also exhibits the lowest values for both C D,∆p and C D,ū 2 (pressure and fluctuating velocity term) for the entire Re range. The mean velocity contribution accounts for 93.19% for MR at Re = 50 and drops down to 74.22% at Re = 1500. Above Re = 500, there is an increase in the contribution of C D,ū 2 (up to 18.45% for Re = 1500) and a slightly increase in the pressure term for MR. CC displays a trend similar to MR, although with slight higher values for the mean velocity and lower values of the fluctuating velocity contribution at high Reynolds numbers. For CC, C D,∆p stays below ≈ 6% for the simulated flow regimes.

C. Wake topology and momentum flux for extreme orientations of MR and CC
Previously, Tagliavini et al. 2021 showed that the fluttering of the particles exhibiting chaotic falling motions has an appreciable influence on the drag coefficient at moderate/high Re. Therefore, we investigate additional orientations for CC and MR, since they fell chaotically at Re 400 during the experiments. One extreme orientation is studied for CC at Re = 500, 1000, 1500, corresponding to the configuration where the angle between the long axis of the particle and the vertical fall direction is largest. Three additional, extreme positions are considered instead for MR at different time instants during free-fall (Re = 1500). These orientations are limiting cases in terms of angle variations (see Fig. 1 and Fig. 2). The results of these additional simulations are reported in Tab. II, together with their geometrical features (porosity and CSF). Fig. 12 depicts the comparison of the timeaveraged, non-dimensional velocity field U * between the final and the extreme orientation of CC at Re = 500 (Fig. 12a) and at Re = 1500 (Fig. 12b). The extreme orientation displays two large recirculation zones in correspondence with the particle hexagonal plates, causing a "bifurcation" of the velocity deficit which concentrates at these tips and extends further downstream (the furthest for Re = 500). At Re = 1500, the two vortex regions at the hexagonal extremities become more asymmetric. The dissimilarity between the final and the extreme orientation is better illustrated in Fig. 13. While at Re = 500 the final orientation manifests only small structures (fluctuations in a non-turbulent environment) several D max away from the particle (Fig. 8d), the extreme position forms a fully turbulent wake (Fig. 13a) with a pair of hairpin vortices rising from the particle extremities (plates) that become smaller in the far wake. The same type of structures is present at Re = 1500 (Fig.  13b), but more compressed and fragmented further away from the snowflake.
In Fig. 14 and Fig. 15, the time-averaged, nondimensional velocity field and the 3D wake structures of different orientation of MR are shown, respectively. For what concerns U * , the distinctness between the diverse positions is evident from the shape and size of the recirulation zones of each configuration. The largest one can be found for t = 129 (Fig. 14c), while t = 94 has the smallest one (Fig. 14b), but, in this case, a higher velocity deficit extends much further downstream. In general, the wake structures of the extreme orientations look more compact and smaller than the ones of the final orientation. All configurations present a pair of hairpin vortices that are advected, deformed, and dissipated downstream (Fig.  15). The configuration at t = 129 appears to have similar structures compared to the final orientation that breaks down at approximately the same distance downstream of the snowflake.
Finally, Fig. 16 reports the comparison between the final and extreme orientations of CC (a, b) and MR (c) in relation to the contribution of the terms of the wake momentum flux. CC shows a small decrease in the contribution of C D,ū and C D,∆p . For this geometry, the case at Re = 500 has the largest disparity regarding C D,ū 2 (7.80% (final orientation) and 11.79% (extreme orientation)), while the fluctuating velocity supply to the total drag stays approximately the same at higher flow regimes. In MR's extreme orientationsū andū 2 (and ∆p) generally contribute to a greater and lesser extent to the total drag coefficient, respectively, in comparison to the final orientation.

IV. DISCUSSION
The forces experienced by a falling object are closely related to the properties of the wake forming behind that object. Hence, there is an interdependence between the object shape, the forces and the wake flow. Many experimental and numerical studies have investigated the falling behavior and the wake of disks [Auguste et al. 2013;Vincent et al. 2016], cylinders [Toupoint et al. 2019], planar polygons [Esteban et al. 2019], and square cylinders [Dutta et al. 2003] to understand which geometrical feature causes the wake to oscillate and thus produces unsteady falling behavior [Vincent et al. 2016]. More complex shapes are taken into account by Nedic et al. 2015 (fractal plates) and Cummins et al. 2018 (plumed seed of a dandelion). In this view, the snowflake shape plays a paramount role in determining the particle drag, terminal velocity, falling behavior and wake properties. Simple geometries such flat plates [Castro 1971;Steiros and Hultmark 2018] or porous spheres [Emadzadeh and Chiew 2020] have been investigated and compared to their solid counterparts. From these studies, increasing porosity seems to cause an increase in the drag coefficient and to alter wake properties. Moreover, high porosity combined with sharp edges may produce changes in the separation point of the shed shear layers. Dutta et al. 2003 investigated squared cylinder at different inclinations and found that changes in the inclination delay the separation point, reducing the momentum deficit (smaller velocity deficit) and thus decreasing the drag coefficient. The flatness of an object is one of the geometrical features that influence the separation point and it is taken into account in this work with Corey's shape factor [Corey 2019]. The snowflake geometries investigated in this paper are irregular, 3D shapes, thus their wake are much more intricate than those of much simpler shapes, such as ellipsoids, disks, or even polygons. Therefore, more than one geometrical feature must be taken into account at the same time. In our study, as mentioned in Section II B 1, both porosity and flatness (CSF) are investigated.
D1 has the highest porosity ( = 0.62) and at Re 100 manifests the highest drag coefficient (Fig. 4). We highlighted that at low Re the wake of this snow particle contains a separated vortex ring which influences the stable fall of D1 and is reminiscent to the behavior of the dandelion's pappus investigated by Cummins et al. 2018. The branches of the dendrite may act in a similar way as the filament of the dandelion's pappus, thickening the boundary layer between the branches and reducing the pressure downstream, hence increasing the drag compared to particles without dendrite features. AG presents a much lower drag compared to D1. Besides the difference in dendrite features compared to D1, AG Corey's shape factor is much higher (less flat particle). This shifts the flow separation point further downstream, whilst the separated shear layers meet and roll up closer to the particle (Fig. 8), which is generally followed by a decrease in C D [Dutta et al. 2003]. Furthermore, the magnitude of the time-averaged, non-dimensional velocity is slightly higher (U * 0.3, thus smaller velocity deficit) and the lower porosity, compared to D1, brings the recirculation region closer to the particle [Steiros and Hultmark 2018]. An analogous condition can be found for MR, which has a much lower porosity and the highest CSF (CSF = 0.95, Tab. I).
The latter may influence the flow separation (further downstream) and causing a decrease in C D (Fig. 4). In this case, the symmetric vortex rings may impact the steady behavior of MR. CC was observed to fall with spiraling, periodic motion at this flow regime. Its shape has the lowest porosity ( = 0.15) and a high CSF (CSF = 0.72). This snowflake forms four recirculation regions attached to the particle, two of which are found between the hexagonal plates of the columnar crystal. Because of the unsteady free-fall of this particle, the final orientation used in this model is not well representative of the wake flow of a freely moving particle, especially at high Reynolds numbers. At low Re, all the snowflakes show a dominant contribution of the mean velocity term of the momentum flux to the total drag (C D,ū ≈ 90%), while C D,∆p and C D,ū 2 stay below 4% and 8%, respectively. Hence the different way the geometries affect the mean velocity velocity profile in the wake is the dominant factor for the resulting C D . D1, AG, and MR displayed a steady falling behavior during experiments, while CC exhibited a spiraling, periodic motion (Tab. I).
As Re increases (Re = 500), the separated vortex ring in the wake of D1 is now attached to the particle and significantly smaller, albeit with much higher U * (Fig. 6a). A drop in the magnitude of the time-averaged, non-dimensional velocity (high velocity deficit) is seen despite the still large contribution to the total drag of the mean velocity term (Eq. (4)). The low U * region extends much further downstream and is very stable (small contribution of C D,∆p and C D,ū 2 ), as seen in Fig. 11. Fig. 8a depicts the 3D wake of D1 and its vortex threads that originate from the dendrite branches, a wake feature that may contribute to the particle stable fall. Although AG and MR display similar wake features, among other things the presence of a vortex street (Fig. 8b), MR has an abrupt increase in the fluctuating velocity term for MR (Tab. I), while this term remains below 10% of the total drag coefficient for AG. The vortex rings forming behind AG are small and slightly detached from the particle, if compared with the ones behind MR. The separated vortices may promote the stable fall of AG, notwithstanding the high Re, and this wake characteristic is also consistent with AG being more porous than MR. An evidence of the more turbulent wake of MR can be seen in Fig.8c, in which the instantaneous flow field of the wake of AG (Fig. 8b) exhibits more regular patterns of the coherent structures (pairs of horseshoe vortices advected downstream). MR (Fig. 8c) exhibits hairpin vortices only in the near wake instead, since they suddenly dissipate into more chaotic structures a few D max away from the snowflake. This particle shows a chaotic free-fall at this flow regime, whereas AG still manifests a steady falling behavior. At Re = 500, for CC a third pair of non-symmetrical, lowvelocity vortices forms between the hexagonal plates of the capped column (Fig. 6d). In addition, the far wake presents a slightly oscillating motion (Fig. 8d). Like MR, CC fell with chaotic motion at this flow regime during experiments.
At Re = 1500, all tested shapes manifest a similar drag coefficient (C D ≈ 1, see Fig. 4), although having different falling behavior. As shown in Fig  11, a notable difference between the diverse shape is the increasing contribution of C D,ū 2 , which is higher than ≈ 10% for the particles that exhibited unsteady falling motion (AG, MR, and CC), even though MR and CC final orientation might not be representative for their chaotic motion at this flow regime [Tagliavini et al. 2021]. D1 maintains a steady falling behavior and the cause may be related to the absence a welldefined vortex street and to the large and packed coherent structures in the wake. The weak vortex street of D1 is associated with the high porosity of the particle [Castro 1971], for which the main contribution to the total drag is still represented by the mean velocity term of Eq. (4). AG fell with periodic motion at Re = 1500 and that may be inferred by the weak shedding and the moderate fluctuating velocity supply of Eq. (4) (slightly higher than the one of D1). The velocity deficit contribution for D1 is still the largest (see Tab. I), but the other particles have a higher contribution from the pressure and fluctuating velocity terms, which compensate the lower velocity deficit and match the drag of D1 at this high Re (Fig. 4). The lowest flatness (highest CSF) of MR might be the cause of the convoluted structures in the wake (Fig. 9c) which could relate to the largest angles variation (Fig. 2). All these factors suggest a strong chaotic motion of this particle. Similar features are also present in CC with slightly lower values. As for Re = 500, the orientation of MR and CC varies significantly during the free-fall and therefore additional, extreme orientations are examined thereafter.
A link between wake turbulence and drag force can also be inferred from Fig. 10. As the Reynolds number increases, std(F D )/F D and T I become larger, with a sharper increment for particles that presented chaotic falling behavior. Moreover, the upsurge of T I (fluctuation of u ) is consistent with the increase in the contribution of u term in the total drag coefficient (Tab. I) and angle variations (Fig. 2).

A. Comparison between final and extreme orientations at moderate and high Reynolds numbers
Regarding the extreme orientations of CC, the considered configuration (at different Re) marks a signif-icant change in the wake structure (Fig. 13) compared to the final orientation, especially at Re = 500. The decrease in the contribution of the mean velocity term and the increase of the pressure and fluctuating velocity contribution may relate to the particle chaotic falling motion and could not be captured by solely investigating the final orientation. The increase in C D can be caused by the higher porosity (Tab. II) of this extreme position (CSF does not change significantly). At higher flow regimes (Re > 500), beside a small increase in the drag coefficient (Fig. 4), there are no meaningful variations of the contribution of Eq. (4) terms in comparison with the final orientation.
The extreme orientations of MR, at different time instants (same Re), present a much lower drag coefficient if compared to the final orientation (Tab. II, Fig. 4). This is consistent with an increase inū and decrease inū 2 (and ∆p) contributions, respectively. Since these configurations are reconstructed at different instants during the particle free-fall, they do not deliver crucial information on the drag coefficient, but help to highlight how the wake structures vary (Fig.  15) while MR falls. In fact, the wake structures become more and more similar to the final orientation, which can be considered as an "averaged" orientation in case of particle with unsteady behavior like MR (see Fig. 15c). This also emphasizes the limitation of taking into account only the final orientation when studying particles with chaotic falling behavior.
We want to stress that the model concerns a fixed, complex-shaped particle impinged by airflow and, despite the adoption of DDES as an established and accurate technique to simulate bluff body flows [Celik 2003], the model presents some limitations in representing the wake flow of a freely falling particle. The main one regards the absence of two-way coupling between the fluid and the particle. For particles that are not falling chaotically the drag coefficient of a freely falling particle is well reproduced by our fixed particle simulations [Tagliavini et al. 2021], which suggests that wake features in the simulations are representative of the coupled system. For chaotic behaviors this is probably less the case, however, the wakes of fixed particle at their extreme orientations provide insights on the range of variability that is exhibited by wakes of particles falling chaotically. The accuracy assessment of the present model in reproducing the wake of a freely falling particle, by comparison with velocimetry experiments, will be the subject of future publications.

V. CONCLUSIONS
The objective of this paper was to gain a deeper understanding of the wake flow behind complex shaped snow particles, its impact on drag and possible implications for falling behavior. To this end a DDES model, validated with experimental data of 3D-printed falling snowflakes [McCorquodale and Westbrook 2021b;Tagliavini et al. 2021], was used. From that, the wake flow is analyzed and the wake contributions to the drag are quantified.
At low Re, particles with highest drag were the one with the highest porosity (D1), which also produced a separated vortex ring, in agreement with the work of Cummins et al. 2018. Lower CSF (Corey's shape factor) is consistent with a delay in the separation point, thus causing the drag coefficient to drop. For low Reynolds number cases the mean velocity term in Eq. (4) contributed to ≈ 90% of the total drag.
For moderate/high flow regimes, porosity seemed to have a role in suppressing (D1) or reducing (AG) the vortex shedding. This lowered the fluctuations in the velocity. As a consequence, particles with higher porosity still exhibited a large contribution of C D,ū to the total drag. An increase in the contribution of the fluctuating term was seen in correspondence with the appearance of vortex shedding. High vorticity and fluctuations in the velocity were related to high CSF, large angle variation (unsteady falling behavior), and large fluctuations in the drag force.
Additional numerical modeling utilizing a range of orientations observed during experiments provided a more realistic representation of the wake flow field for CC and MR at moderate/high Reynolds number, at which the particles are known to exhibit unsteady motions in free fall. They highlighted the limitations of the fixed particle model that considers only one orientation (final orientation) when investigating snowflakes that exhibit unsteady falling motion. The contribution of the different terms of Eq. (4) did not change significantly for these extreme orientations, despite some differences in C D which might be related to the changing in projected area due to the diverse particle orientations.
Our study can serve as a template to shed light on the influence of shape on the wake flow and on how the latter affects the drag coefficient and the falling behavior of complex-shaped snow particles over a broad range of geometries.

ACKNOWLEDGMENTS
The numerical simulations were performed using ETH -Euler cluster resources of Prof. Ro-    (C D,tot ) of each particle at different Reynolds numbers (Re), evaluated from the numerical model using the momentum flux integration (Eq. (4)). Each value is decomposed according to the momentum deficit terms to assess the contribution of each specific term, which is expressed as percentage. The falling behavior of the particle for different Re is also specified (S = stable, P = periodic, and C = chaotic).

Particle
Re Falling behavior C D,tot C D,ū

FIG. 4:
Time-averaged C D values as a function of the Reynolds number (semi-log scale) for each particle at its final orientation, as reconstructed from experiments. The empty markers represent the drag coefficients of the extreme orientations. For particle CC, the same extreme orientation is tested at Re = 500, 1000, 1500, while for MR different extreme positions at different instants during free-fall are simulated at the same Reynolds number (Re =1500). The error bars represent the error between the time-averaged drag coefficient computed from the simulations and the drag coefficient evaluated from the momentum flux integration ((C D,sim − C D,mom )/C D,sim · 100). The averaged percentage error is 6.82%, 7.90%, 9.69%, 8.85%, 6.49%, and 9.59% for D1, AG, CC, MR, CC ext , and M R ext , respectively. The 3D representation is obtained by using the Q-criterion with Q > 1 and colored according to the streamwise, non-dimensional velocity u * x . The 2D illustration shows the non-dimensional vorticity magnitude ω * = (D max ω)/U ∞ for (a) D1, (b) AG, (c) MR, and (d) CC. The selected plane is chosen according to the wake maximum width.   (8)) to the total drag coefficient of particle CC (at Re = 500 (a) and at Re = 1500 (b)) and MR (at different time-steps (c)) and their final and extreme orientations, respectively.