the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Stable and unstable fall motions of platelike ice crystal analogues
Jennifer R. Stout
Christopher D. Westbrook
Thorwald H. M. Stein
Mark W. McCorquodale
The orientation of ice crystals affects their microphysical behaviour, growth, and precipitation. Orientation also affects interaction with electromagnetic radiation, and through this it influences remote sensing signals, in situ observations, and optical effects. Fall behaviours of a variety of 3Dprinted platelike ice crystal analogues in a tank of water–glycerine mixture are observed with multiview cameras and digitally reconstructed to simulate the falling of ice crystals in the atmosphere.
Four main falling regimes were observed: stable, zigzag, transitional, and spiralling. Stable motion is characterised by no resolvable fluctuations in velocity or orientation, with the maximum dimension oriented horizontally. The zigzagging regime is characterised by a backandforth swing in a constant vertical plane, corresponding to a time series of inclination angle approximated by a rectified sine wave. In the spiralling regime, analogues consistently incline at an angle between 7 and 28°, depending on particle shape. Transitional behaviour exhibits motion in between spiral and zigzag, similar to that of a falling spherical pendulum.
The inclination angles that unstable planar ice crystals make with the horizontal plane are found to have a nonzero mode. This observed behaviour does not fit the commonly used Gaussian model of inclination angle. The typical Reynolds number when oscillations start is strongly dependent on shape: solid hexagonal plates begin to oscillate at Re =237, whereas several dendritic shapes remain stable throughout all experiments, even at Re > 1000. These results should be considered within remote sensing applications wherein the orientation characteristics of ice crystals are used to retrieve their properties.
 Article
(6625 KB)  Fulltext XML

Supplement
(29577 KB)  BibTeX
 EndNote
Understanding the motion of falling ice crystals is important to both the microphysical processes within clouds and their bulk characteristics, such as radiative and optical properties. However, their dynamics are not well understood; ice crystals have complex and irregular shapes and can exhibit fluttering, spiralling, and tumbling motions.
To quantify the orientation of analogues, the inclination angle, θ, is the angle made between the rotated ice crystal's c axis and the global vertical y axis (Fig. 1). Falling ice crystals, when stable, have a constant θ of 0° (List and Schemenauer, 1971). When unstable, it is commonly assumed crystals have Gaussian distributions of orientations, with a modal θ of 0°, and standard deviations varying between 10° (pristine ice crystals) and 40° (heavily aggregated snowflakes) (Melnikov and Straka, 2013; Ryzhkov et al., 2020). However, ice crystals exhibit a variety of unstable falling regimes, each corresponding to different distributions of orientations.
1.1 Importance of orientation of ice crystals
The orientations of falling crystals impact their projected area in the horizontal plane, their sedimentation rate, and the rate at which they can collide with other hydrometeors (Westbrook et al., 2010). Compared to ice crystals with purely vertical motion, ice crystals with horizontal motions in addition to the vertical will travel a farther distance, providing more opportunity to collide with other hydrometeors than ice crystals with vertical motion alone (Wang, 2021). This further impacts cloud macrophysical properties, such as radiative impacts and cloud lifetime.
Properties of ice crystal motion have important implications for radar and lidar observations: orientation directly influences signals sampled by dualpolarisation radar, as the orientation of crystals changes the differential reflectivity (Z_{DR}) (Bringi and Chandrasekar, 2001). Unstable motion causes fluctuations in the crystal velocity in the component of the crystal motion along the radar beam, broadening the Doppler spectrum width (Feist et al., 2019). Differences in assumptions of orientation can therefore impact the relationship between the derived ice crystal diameter and Doppler or polarimetric remote sensing observations, ultimately affecting radarderived precipitation rates (Matrosov, 2011; Schrom et al., 2023).
Horizontal orientation of ice crystals affects lidar observations, especially in the case of specular reflection, causing enhanced return for lidars pointing exactly at zenith or nadir (Sassen, 1977; Platt, 1977; Gibson et al., 1977; Hogan and Illingworth, 2003). The magnitude of the enhancement and its variation with elevation angle are strongly dependent on the chosen model for crystal orientation (Platt, 1977).
For clouds containing ice, crystal size, concentration, habit, and orientation all play a significant role in determining cloud radiative properties such as optical depth and albedo (Curry and Ebert, 1992; Ishimoto et al., 2012; Hirakata et al., 2014). Changes in these particle orientation assumptions can lead to high variation in the retrieval of cirrus properties from satellite observations. In certain cases, decreasing the assumed standard deviation of θ from 20 to 5° doubled the estimated optical depth (Masuda and Ishimoto, 2004). Horizontally oriented ice crystals have also been theorised to increase cloud shortwave albedo by up to 40 % (Takano and Liou, 1989).
When ice crystals are horizontally oriented, this gives them distinctive optical characteristics (Cho et al., 1981; Sassen, 1987). For instance, horizontal crystals can create a range of atmospheric optical phenomena such as sun dogs, light pillars, and Parry arcs, among others (Moilanen and Gritsevich, 2022). Additionally, spiralling ice crystals have been hypothesised to cause the rare “Bottlinger's rings” effect (Lynch et al., 1994; Tränkle and Riikonen, 1996).
Ice crystal orientation also impacts the apparent crystal properties (e.g. size, projected area, aspect ratio) inferred from analysis of 2D projections sampled by groundbased imagers such as PIP (Precipitation Imaging Package) (Jiang et al., 2017; von Lerber et al., 2017). To estimate the 3D parameters relevant for drag calculations from 2D projections of snowflakes, assumptions about particle orientation, shape, and motion must be made (Köbschall et al., 2023). Dunnavan and Jiang (2019) find that for highly eccentric particles (such as aggregates) that have large fluctuations in θ, very limited information can be inferred about a particle's 3D shape without specifying appropriate particle orientation distributions.
1.2 Phenomenology of circular discs
Analogies may be drawn between the aerodynamics of ice crystals and those of other idealised shapes, such as thin circular discs. There has been extensive experimental research on the aerodynamic behaviour of thin circular discs (e.g. Willmarth et al., 1964; Field et al., 1997; Ern et al., 2011; Zhong et al., 2011). Two dimensionless ratios have been proposed to characterise the motion of falling circular discs: the Reynolds number, Re, and the dimensionless moment of inertia, I^{*} (Willmarth et al., 1964; Field et al., 1997), discussed in the following subsections.
1.2.1 Reynolds number
The Reynolds number is defined as
where υ is the dynamic viscosity of the fluid, V_{mean} is the mean vertical velocity of the particle, and D is the maximum dimension of the particle.
Willmarth et al. (1964) identified that Re 100–200 is the critical point for the onset of unstable motions for circular discs, after which periodic behaviour begins. The value of the critical Reynolds number varies depending on particle shape. Field et al. (1997) report an experimental study of how metal circular discs fall through water and glycerol mixtures and how paper discs fall through air. Different falling regimes were observed depending upon the experimental parameters; discs could fall steadily, exhibit oscillating periodic motions, or tumble.
Periodic behaviour includes zigzag and spiralling subtypes of behaviour, and more recently an inbetween behaviour was identified as transitional through experimental investigations by Zhong et al. (2011). Zhong et al. (2011) find that for circular discs, at low Re and I^{*}, the most common behaviour is spiralling, whereas at high Re and I^{*}, zigzagging behaviour is most common, with transitional behaviour occurring at intermediate I^{*} and Re (Fig. 2).
At Reynolds numbers below the critical Reynolds number, flow around crystals is stable. For planar crystals in this regime, laboratory and field measurements have shown that the largest dimension becomes normal to the axis of gravity, and the plate crystals achieve a horizontal orientation, corresponding to a constant inclination angle of zero (Jayaweera, 1965; List and Schemenauer, 1971; Pruppacher and Klett, 1997 (see their Sect. 10), Kajikawa, 1992; Wang, 2021)
At a critical Reynolds number, the flow around the crystal becomes unstable, forming vortices as part of the boundary layer of fluid at the surface of the particle. When shedding of these vortices in the wake of crystals occurs, the distribution of pressure on the crystal changes, exerting forces that cause it to rotate (Zhong et al., 2013; Tagliavini et al., 2021a). These unstable motions are observed as oscillations in orientation, as well as in the vertical and horizontal velocities, such that they are nonzero, fluctuating, and have a distribution. There is a current lack of understanding about the orientation of ice crystals in unstable regimes, and one of the aims of this paper is to explore this.
1.2.2 Dimensionless moment of inertia, I^{*}
The dimensionless moment of inertia for a circular disc is defined as the ratio of the moment of inertia of a circular disc about its diameter and a quantity proportional to the moment of inertia of a rigid sphere of fluid of the same diameter, such that
where ρ_{p} is the density of the particle, ρ_{f} is the density of the fluid, t is the thickness of the disc, and D is its diameter (Willmarth et al., 1964). For more complex shapes, the more general nondimensional moment of inertia is
where D is the maximum dimension of the particle. The moment of inertia for rotation around the three principal axes of the crystals is calculated, where I_{a} is the smallest of these three moments, aligned in the aaxis direction of the crystal (Fig. 1) (Kajikawa, 1992). For reference, further information regarding the calculation of the moments of inertia can be found in Gregory (2006, p. 570).
1.2.3 Comparing ice crystals with discs: Re – I^{*} phase space and area ratio
Figure 2 presents a summary of the coverage of the data presented in McCorquodale and Westbrook (2021a, b), also used in this study, on the Re and I^{*} phase space and the key prior experiments on ice crystal shapes (Kajikawa, 1992; Cheng et al., 2015; Nettesheim and Wang, 2018) and circular discs (Field et al., 1997; Zhong et al., 2011) that are discussed in Sects. 1.3 and 1.4. Using a mass–diameter relationship from Nakaya and Terada (1935) for planar dendritic crystals and methods for estimating I^{*} from Kajikawa (1992), we find that a 10, 1, and 0.1 mm planar dendritic crystal, where the density of ice is 917 kg m^{−3} and the density of air is 1.2 kg m^{−3}, has an I^{*} of 0.02, 0.2, and 2.0 respectively. The hatched region of Fig. 2 displays this expected range of I^{*} for planar dendritic ice crystals – this matches well with the range of previous observations of ice crystals.
Our study focuses on planar crystals. These range from hexagonal plates (which present a solid obstacle to the flow at all Reynolds numbers) to stellar crystals and dendrites, which have much more open projections. One way to characterise this shape variability is by area ratio: the ratio of the maximum crosssectional area of the particle and the area of its circumscribing circle, which helps compare with circular discs. Fluid experiments on planar shapes report that the amplitude of oscillations in the descent velocity was maximum for circular discs and decreased with the area ratio, suggesting that unstable motions are inhibited by more complex shapes (Esteban et al., 2018).
1.3 Existing work on orientation of ice crystals
A variety of approaches have been developed to study the aerodynamics of ice crystals. Using the CloudAerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO), Zhou et al. (2013) simulated crystal distributions and orientations and found that horizontally oriented plates occurred in 60 % of optically thick ice and mixedphase cloud layers. Similarly, Stillwell et al. (2019) found that horizontally oriented plates must occur in at least 25.6 % of all iceonly column observations using polarisation lidar for their simulations to match the observations.
Common models of particle orientation distribution assume either uniform distribution, horizontal orientation with an inclination angle of zero, or a Gaussian distribution with a peak at zeroinclination angle (e.g. Borovoi and Kustova, 2009). However, models of particle orientation distribution for falling particles suggest modal inclination angles of approximately 10° (Klett, 1995). Melnikov and Straka (2013) attempt to retrieve the spread of fluttering angles from polarimetric radar data, assuming that the mean inclination angle is zero and that the distribution has a fixed, sizeindependent width, retrieving fluttering amplitudes on the order of 2–23°. However, remote sensing is an indirect measurement rather than a direct observation of the fall motion.
More direct measurements are possible, such as in situ observations of ice particles near the surface (e.g. Zikmunda, 1972; Locatelli and Hobbs, 1974; Kajikawa, 1992; Garrett et al., 2015; Fitch et al., 2021). Falling natural planar snow crystals placed into a tube were studied by Kajikawa (1992), who used a stereophotogrammetric method and found that hexagonal plate crystals exhibited stable and unstable motions, including a swing motion (zigzag) and a helical rotation motion (spiralling). The critical Reynolds number, above which crystals exhibited unstable motion, was found to vary depending on the specific crystal habit, as classified by Magono and Lee (1966). For crystals of classification P1a (hexagonal plates), the critical Reynolds number was found to be 47, while for P1f crystals (fernlike crystals), the critical Reynolds number was found to be 91. Nonetheless, the tendency of ice crystals to break, evaporate, and melt when handled led to high uncertainties in direct observations of ice crystals at the ground.
MultiAngle Snowflake Camera (MASC) observations by Garrett et al. (2015) found that the modes of the distribution of inclination angles were 20, 16, and 13° for graupel, rimed particles, and aggregates respectively, indicating that snow particles show a preference for nearhorizontal orientation but have nonzero modal values. Recent research into the MASC measurements by Fitch et al. (2021) has also reported preferential nonhorizontal inclinations for the orientation of snow particles, with a modal value of 12° observed for light wind speeds in shielded conditions.
These findings suggest that the assumption of Gaussian orientation distribution may not always be accurate and that the orientation of snow particles may exhibit preferential orientations that are nonhorizontal, even in quiescent environments. Lynch et al. (1994) proposed modelling the swinging motion of falling ice crystals similarly to that of a pendulum where the pivot of the pendulum falls vertically at constant velocity. This notion is supported by Esteban (2019), who found that oscillatory motions of discs and a variety of other planar shapes in both quiescent and turbulent fluids had pendulumlike motions, with turbulence simply adding noise to the oscillations.
We hypothesise that planar ice crystal analogues will behave similarly to this previous experimental work and test that hypothesis in this study. Falling ice crystals may be well approximated as falling pendulums, and there is a relationship between the distribution of angles and other fall motion aspects such as velocity fluctuations, perceived projected areas, and perceived aspect ratios.
Cheng et al. (2015) explored the behaviour of hexagonal plates using numerical simulation, with Re ranging from 46 to 974 and I^{*} ranging from 1.1 to 0.3. The plates are stable at Re = 46 and unstable at Re = 135. Smaller plates exhibit a zigzag motion, while larger plates exhibited spiralling, and none of the plates tumbled during the simulation, in contrast to the work by Field et al. (1997) and Zhong et al. (2011) on circular discs.
Nettesheim and Wang (2018) used numerical simulations to study the fall behaviour of branched crystals, showing unstable fall motions for sector plates at Re =384 and broadbranched plates at Re = 345. They also provided data on other experiments, with sine wave fits to the Euler angles of the particles during the experiments. However, the time series the angles are fit over include a spinup period between the initial “release” of the crystal and it settling into its preferred fall motion, which precludes a quantitative comparison with the results presented in our study.
As the unsteadiness of falling particles is a complex, nonlinear, multidegreeoffreedom phenomenon, numerical simulations impose significant computational cost and technical challenges. These simulations also rely on assumptions about turbulence, vortex shedding, and how these interact with falling particles, making it difficult to confidently simulate the wide range of conditions ice crystals experience.
Using analogues – scaledup models of natural crystals – presents a promising avenue for studying the fall behaviour of ice crystals in a laboratory environment. List and Schemenauer (1971) report measurements of machined analogues of snowflake particles falling in solutions of water and glycerine or salt water and exploit the dynamic similarity. This dynamic similarity only applies when falling steadily at terminal velocity, since the only dimensionless variables are Re and particle shape. When falling unsteadily, the ratio ${\mathit{\rho}}_{\mathrm{p}}/{\mathit{\rho}}_{\mathrm{f}}$, contained within I^{*}, is also significant. To compare the results with natural snowflakes falling in the atmosphere, the study considers five different designs of planar ice crystals and observes stable behaviour at Re < 100. For discs, hexagonal plates, and broadbranched models, small oscillations are observed at Re ≈ 200, although no oscillations are observed at this Reynolds number for stellar, dendritic, or stellarwithplate shapes. Köbschall et al. (2023) used analogues of aggregate snowflakes, finding that the area of complex snowflake analogues projected in the direction of flow is often maximised, and for many of their analogues, a rotation around the vertical axis was seen.
Building on previous work by Westbrook and Sephton (2017), McCorquodale and Westbrook (2021a) utilised modern 3D printing techniques to fabricate analogues for studying the aerodynamics of ice particles through the analogue method. In experimental studies, these analogues were analysed through a custom algorithm, producing digital reconstructions of the trajectory and orientation of the particle. From these experiments, analogues of aggregates are found to exhibit different preferential orientations depending on the Reynolds number for the same particle shape (McCorquodale and Westbrook, 2021c). Tagliavini et al. (2022) performed numerical simulations with dendritic crystals and compared results with freefalling analogues, using the particle tracking algorithms described in McCorquodale and Westbrook (2021a, b). They found that throughout the Re range in both numerical simulations and laboratory observations, the wake and motions of dendritic crystals were stable, even as high as Re = 1500, supporting the idea that the onset of unstable motions is sensitive to crystal geometry. This is a topic explored in the current article.
1.4 Investigating unresolved questions
It is evident that the representation of crystal orientation in many studies is not well constrained at present. There is evidence that unstable motions may be more complex than a simple zigzag motion, but the conditions under which this happens are not clear.
There are extremely limited data quantifying how the orientations of unstable crystals are distributed and what that distribution depends on as well as how frequent and large the velocity fluctuations (in both vertical and horizontal) are in response to the unstable wake of the falling crystal and how they are correlated with the variations in orientation. In this article we present new data to address these areas of uncertainty.
Building on previous work by Westbrook and Sephton (2017), McCorquodale and Westbrook (2021a) utilised modern 3D printing techniques to fabricate analogues for studying the aerodynamics of ice particles through the analogue method. To link the behaviour of real ice crystals to the theoretical behaviour observed by Esteban et al. (2019, 2018) in laboratory experiments, we further examine the experiments by McCorquodale and Westbrook (2021a), focusing on the fall behaviour of quiescent platelike particles, identify the angles at which ice crystal analogues fall, and test the potential relationship between the distribution of fall angles and other motion aspects.
The paper is organised as follows: in Sect. 2, we describe the experiment by McCorquodale and Westbrook and the data sets derived from it. In Sect. 3 we discuss the results, beginning with Sect. 3.1, discussing which particles fall steadily. Section 3.2 introduces and describes four case studies of periodic motion and how their orientations, velocities, and oscillation frequencies can be characterised. Section 3.3 discusses the broader trends and characteristics of the full data set, including how distributions of θ, oscillation frequencies, and motion type vary by shape and Reynolds number, as well as how velocity components vary with one another. Further discussion of these results, including a summary and conclusions, can be found in Sect. 4.
A diverse range of ice particle analogues were included in this study, ranging from hexagonal plates with an area ratio of 0.87 to open branched crystals with an area ratio as low as 0.23 (Table 1). The area ratio of the particles included is calculated using the observed projected area of the particle divided by the circumscribing circle at each time step during experiments when fall motion is stable. The mean calculated area ratio is then used to describe each result.
The ice crystal analogues were produced using a Form 2 3D printer (Formlabs), which achieves a high level of precision with a minimum layer thickness of 25 µm and a laser spot size of 140 µm. The maximum dimensions of particles ranged from 1 to 3 cm, with aspect ratios varying between 0.04 and 0.2 and area ratios varying between 0.2 and 1 (Table 1). Due to an artefact of how the numerical code from Reiter (2005) was used to create some of the crystal shapes, a few of the models (S, F, D, DP, and PB) were later realised to be nonhexagonally symmetric and instead have a horizontal aspect ratio (the diameter in the a axis to the diameter of the a^{′} axis) of 1, instead of 1:1.15 for a regular hexagon. We do not expect this to affect the broad behaviour of their fall motions, and indeed we observe zigzag, spiral, and transitional behaviour for these particles, but as noted later, this asymmetry may influence the details of the critical axis that zigzag motions are oriented around.
To replicate atmospheric conditions in the laboratory, the dynamical similarity experiment was conducted in a transparent acrylic tank with internal dimensions of $\mathrm{0.4}\times \mathrm{0.4}\times \mathrm{1.8}$ m. The tank was filled with uniform mixtures of water and glycerol, with the volume fraction of glycerol ranging from 0 % to approximately 50 %. By varying both the density and viscosity of the fluid, and the size of the analogues, it was possible to sample a wide range of Reynolds numbers for each shape (Table 1).
During the experiment, the ice particle analogues were allowed to free fall through the tank and were recorded using three orthogonal cameras. Each camera records the fall of the particle through a region approximately $\mathrm{0.2}\times \mathrm{0.2}\times \mathrm{0.2}$ m in size, 1.5 m below the surface of the fluid. By this point, the particles have reached their terminal velocities and their behaviour is insensitive to the initial release orientation.
The Trajectory Reconstruction Algorithm implemented through Image anaLysis (TRAIL) then produced digital reconstructions of the trajectory and orientation of the particle in free fall. The orientation of the particles was reconstructed using a set of Euler angles. More details on the fabrication of the analogues, experimental setup, and reconstruction algorithm can be found in McCorquodale and Westbrook (2021a).
These data, referred to as “TRAIL”, provide time series of the 3D positions and orientation of the falling analogues from which the 3D velocity vectors at each time step can be derived. The reconstructed orientations, described by the Euler angles, further enable the calculation of the inclination angle, θ, which is more widely used in atmospheric applications. A total of 354 experiments with platelike shapes were conducted, resulting in the range of values described in Table 1.
The instantaneous velocity at each time step is calculated by applying the central difference formula to the coordinate values, providing an estimate of the instantaneous velocity of the particle at each time point.
Motion observed in the laboratory was typically stable or periodic. Based on the variation of the particle inclination angle, θ (Fig. 1), the periodic behaviour can be divided into three subtypes – zigzag, spiral, and transitional behaviour – and will be analysed below.
3.1 Crystals which fall steadily
Particles in the TRAIL data set were diagnosed as exhibiting stable motion when the Euler angles that describe rotation about the a and a^{′} axes (see Fig. 1) fluctuate by less than ±2.5° across the measurement region; this threshold corresponds to the resolution of the 3D reconstruction. Stable particles fall horizontally with their a axis in the horizontal plane and the c axis oriented vertically (i.e. with a nearzero inclination angle), with no measurable fluctuations in velocity and no horizontal movements.
A total of 223 ice crystal analogues exhibited stable motion, while 131 exhibited unstable, periodic motion. Across all shapes, the Reynolds number alone cannot be used to predict stability: the Reynolds numbers observed ranged from 3 to 1615 for stable motion and from 197 to 2162 for unstable motion. The range of the dimensionless moment of inertia values was 0.14 × 10^{−3}–12 × 10^{−3} for stable motion and 0.28 × 10^{−3}–11 × 10^{−3} for unstable motion. With both variables exhibiting a considerable overlap in the presented behaviours, the onset of unstable motions for ice crystals cannot be considered the same as for circular discs, which become unsteady around Re = 100–200 (Field et al., 1997) and around Re = 200 for our results (Fig. 3).
Shape (approximated by area ratio) is found to have a large impact on instability. The coverage of stable and unstable behaviours for all ice crystal analogues in TRAIL is summarised in Fig. 3, and it can be seen that the onset of unstable motions can be at larger Re (by up to 1 order of magnitude) than the predicted onset of unsteadiness for circular discs. The spread of experiments and their motion types by Reynolds number, separated by shape, is presented in Fig. 3.
Some particles are stable for a much larger range of Reynolds numbers than others. A few shapes (D1 at all aspect ratios as well as D, DP, S, and F at aspect ratio 0.04) remained stable throughout all conditions, even at Re > 10^{3}. Previous studies report an increase in the drag coefficient when planar particles fall unsteadily (McCorquodale and Westbrook, 2021b). That is, the onset of unsteady motion is coupled with a change in wake structure (Zhong et al., 2011; Tagliavini et al., 2021b; Nettesheim and Wang, 2018), which in turn influences the drag coefficient. This change in C_{D} is more pronounced when the area ratio is high than when it is low (McCorquodale and Westbrook, 2021b), suggesting that unsteadiness is less vigorous in particles with low area ratios, such as dendrites.
3.2 Case studies of periodic motion
Four case studies are presented to illustrate the periodic motion subtypes seen in Fig. 4 and described in Table 2. These cases were picked by visual inspection as characteristic types of behaviour. In this section, we will quantitatively describe the four case studies and then objectively classify their motion based on inclination angle in Sect. 3.2.5.
Each case study experiment was conducted in pure water. The case study examples are of hexagonal plates except the spiralling case (Fig. 4g, h), which was a broadbranched plate, as none of the hexagonal plate studies exhibited pure spiralling behaviour with no wobble, but instead exhibited transitional spirals. Figure 4a, c, e, and g present side views of the particle cases, viewed from a laboratory frame of reference. Figure 4b, d, f, and h present the linearly detrended centre of mass of each particle at each time step, effectively subtracting the mean fall velocity, such that the particle is viewed from an observer falling at the same mean velocity as the particle.
3.2.1 Characteristics of periodic motion
The first of the periodic motion types seen is the zigzag case study: the particle swings back and forth in one plane, and as the particle swings away from its centre of fall, its inclination angle increases, akin to a planar pendulum motion. The zigzag–transitional case introduces an element of rotation around the vertical axis, such that the plane of swing slowly moves anticlockwise, and had the experimental run been longer, it may have rotated back to its original position. The transitional–spiral case is similar, but the rate of rotation around the vertical is faster, producing wider loops. Spiralling, the final subtype of periodic motion remains at a nearconstant inclination angle and does not swing back and forth and instead precesses around its central point without touching its mean centre of fall.
The subtypes of periodic motion can be approximated by the subtypes of spherical pendulums: zigzagging is similar to a planar pendulum, spiralling is comparable to a conical pendulum, and transitional motion captures the range of pendulum motion between the two extremes, with the horizontal displacement approximating a rhodonea curve (Helt, 2016). A conical pendulum characteristically traces out a circle in the horizontal plane, akin to spiralling cases, which also trace out a circle in the horizontal plane. Similarly, a planar pendulum serves as an analogy for the zigzagging motion, as they are both constrained to movement in a single plane, tracing out a line in the horizontal plane.
3.2.2 Time series of θ and ϕ
Series of inclination angles, θ, for the periodic motion types can be approximated as sinusoidal (Fig. 5). To distinguish between the regimes, rectified sine waves are fit to these inclination angle time series, using a fast Fourier transform as a first guess of the frequency of the sine wave and the SciPy curve fit function (Virtanen et al., 2020), such that
where θ is the inclination angle, t is the time in seconds, θ_{amp} is the amplitude of the wave, and θ_{tilt} is the angular displacement of the sine wave; the period of the sine wave is $\mathrm{2}\mathit{\pi}/\mathit{\omega}$ (in seconds) and ω_{0} is the phase shift.
All unstable motion presented in this study was observed to be periodic and is approximated through Eq. (4). We observed no complex tumbling, chaotic fluttering, or behaviour with significantly nonsinusoidal motion to it. However, we note that there are some experiments in the study in which additional (weaker) modes of oscillation, in addition to the primary frequency, seem to be present. For example, the transitional case in Fig. 5c has an amplitude that fluctuates slightly in time at a lower frequency than the primary mode of oscillation captured by the simple singlefrequency fit. We did not attempt to capture these finer details in our fitting procedure.
The root mean square error of the rectified sine wave fits to the four case studies is 1.1, 1.3, 2.7, and 0.9° for the zigzag, transitional–zigzag, transitional, and spiral cases respectively. The root mean square error of all fits to the data is 1° and is provided in the Supplement.
θ_{amp} and θ_{tilt} are found to summarise the motion types well, as they represent the variability and tilt of the particle respectively and constrain the pendulum model. They also allow the periodic subtypes to be distinguished quantitatively: a spiralling particle has a low θ_{amp} and a high θ_{tilt}, as it is consistently inclined and does not flutter.
Zigzagging behaviour is the opposite: a potentially high θ_{amp} and a nearzero θ_{tilt}, as it swings around a horizontal orientation, but flutters much more than a spiralling particle (Fig. 5). For instance, the zigzag example case (Fig. 5a) has a fitted θ_{amp} of 34° and a θ_{tilt} of 0°.
The transitional–zigzag case behaves similarly but never samples (close to) θ=0, and it is just starting to transition to having a slow rotational component. It should be noted that although the transitional–zigzag case has a higher θ_{amp} (42°) than the zigzag case (where θ_{amp} is 34°), this does not negate categorisation of behaviour for either case, as the trajectory of the transitional–zigzag case is close to zigzag behaviour but never samples exactly θ=0, and it is just starting to transition to having a slow rotational component.
The transitional case (Fig. 5c) has a smaller amplitude than both zigzag and zigzag–transitional cases; θ_{amp} is 9°, but θ_{tilt} is much higher, at 24°. The spiralling case (Fig. 5d) has an even smaller amplitude still, with θ_{amp} at 1°but θ_{tilt} at 20°. Spiralling behaviour can have nonzero θ_{amp}, although θ_{amp} is small. This small variation in θ is referred to as “wobble” for spiralling cases. It may be worth considering how much of this wobble is a physical phenomenon versus an artefact or experimental uncertainty; a wobble of 2.5° could easily originate from experimental uncertainty. Given the wobble in Fig. 5d appears to have a uniform period, we believe the wobble in this case is partly a physical phenomenon, but you can see the impact of experimental uncertainty in the traces within Fig. 5a–c at the limits of the inclination angle (e.g. for zigzag the angle often does not reach 0). This will partly be due to the finite time resolution of the measurements and partly due to the accuracy of the orientation reconstruction.
Figure 5 also displays the distributions of θ for each case study. The zigzag case has a distribution that is consistent with that expected for simple harmonic motion, where the most likely inclination is at the end of each swing where $\mathrm{d}\mathit{\theta}/\mathrm{d}t$ is smallest and least likely is an angle of zero. Spiral has an almost constant inclination angle, and hence a very narrow distribution of θ, centred on an angle significantly higher than zero. Transitional has a distribution that is in between the other two cases, and in common with spiral cases θ is always above zero. This is in significant disagreement with the common assumption that orientation is a Gaussian distribution where the most common orientation is horizontal (θ=0).
An azimuth angle, ϕ, represents where the c axis of the crystal when projected into plane view is pointing relative to the x axis (in the laboratory reference frame) (Fig. 1c). A spiralling particle has a linear increase (or decrease, in cases not shown) of azimuth angle, and $\mathrm{d}\mathit{\varphi}/\mathrm{d}t$ is constant (as seen in Fig. 6d). The sawtooth shape is produced by the angle being limited to ±180°.
Purely zigzag cases swing around one axis: when the particle goes from pointing one way to pointing another, ϕ changes by Δ180° seen by the squarewave shape (Fig. 6a). For pure zigzag cases, $\mathrm{d}\mathit{\varphi}/\mathrm{d}t$ is constant except for close to the instant where the particle becomes horizontal (θ=0) and ϕ becomes highly uncertain, as the c axis is momentarily pointing towards the vertical (and can be ignored). The axis that zigzagging cases pivot around tends to be the branches of the crystal, in the plane of the a and a^{′} axes. For shapes that are nonhexagonally symmetric (S, F, D, DP, and PB), the shortest branches are the axis the crystal pivots around.
Transitional cases have a nonconstant $\mathrm{d}\mathit{\varphi}/\mathrm{d}t$, a combination of the sawtooth and square waves seen in the spiral and zigzag cases. For the zigzag–transitional case (Fig. 6 b), ϕ increases during its time along each loop of the rhodonea curve and then jumps by a value close to 180° when the orientation of the c axis is close to vertical. The transitional case (Fig. 6c) has no visible jumps in ϕ, except for the aliasing at ±180°, as θ never becomes close to zero.
3.2.3 Velocity fluctuations
The amplitudes of the u and w components of velocity (in the x and z directions respectively) changed in the presence of any component of rotation around the vertical. Therefore, a combined horizontal component of velocity, U, was calculated as follows:
Sine waves can be fit to the vertical component of velocity, V, and the horizontal component of velocity, U, such that
These velocity components of the particles were found to follow a sinusoidal pattern consistent with the pendulum model (Fig. 7). U is fit with a rectified sine wave, as horizontal speed can become zero in zigzag cases but cannot be negative (Eq. 5).
For all experiments with periodic motion, θ, U, and V are all found to have sinusoidal patterns with the same period (i.e. ω_{V}=ω_{U}) but different offsets and amplitudes. U_{offset} is not always zero: it is nonzero for particles that drift as they descend or particles with a spiralling component.
A summary of the sine wave fit components can be found in Table 3. Horizontal velocity, U, was found to peak when the tilt was lowest (i.e. when the particle was flat or the inclination angle of the particle was closest to zero) just after vertical velocity, V, peaks.
The amplitude of V fluctuations relative to the mean vertical fall speed (${V}_{\mathrm{amp}}/{V}_{\mathrm{mean}}$) is 0.37, 0.47, 0.42, and 0.03 for the zigzag, zigzag–transitional, transitional, and spiral cases respectively. The amplitude of horizontal velocity fluctuations relative to the mean fall speed (${U}_{\mathrm{amp}}/{V}_{\mathrm{mean}}$) is similar to their vertical components (0.41, 0.55, 0.49, and 0.03 for each case respectively). In the case of the spiralling particle, U and V are held relatively constant compared to the other cases, effectively making the spiralling cases a quasisteady mode with a nonzero nearconstant inclination.
For the nonspiralling cases, large fluctuations suggest that the mean vertical velocity, V_{mean}, does not sufficiently characterise the velocity of a particle. The distributions shown in Fig. 7 for the nonspiralling cases are broad, suggesting that a broad spectrum of instantaneous velocities for a single type of oscillating particle should be considered when interpreting Doppler spectra.
3.2.4 Projected area fluctuations
One major application of this research is for dwelling radars and lidars, whether groundbased (usually close to zenith) or spaceborne (usually close to nadir). Variation in the projected area affects assumptions in backscatter crosssection and hence the retrieval of particle size and number. Fluctuations will also affect polarimetric measurements and retrievals using these, particularly when the particles are viewed from the side. Figure 8 presents the time series of reconstructed projected areas as seen from below for each case study, normalised by the planar crosssectional area of each analogue. The projected area as seen from below anticorrelates with θ in each time series, displaying an outofphase relationship: as θ increases, the aspect ratio of the analogues increases (becoming closer to 1), while the area ratio and projected area decrease.
For all four cases, vertical velocity and projected area exhibit a 180° outofphase relationship, where vertical velocity peaks just after the projected area reaches its minimum point. The projected area is more in phase with horizontal velocity, peaking just after U reaches its maxima. This supports the idea that vertical velocity increases when the projected area is minimised, as drag is minimised in the vertical direction, allowing the particle to accelerate.
In the case of an infinitely thin particle, the projected area as seen from below is equal to the crosssectional area of the particle multiplied by cos θ (and hence correlated with cos θ, as evident from Fig. 8). In the presence of particle thickness, the normalised projected area is expected to be greater than or equal to cos θ and therefore should never fall below 0.7, as θ never exceeds 45° for periodically oscillating cases. In fact, the observed projected area occasionally exceeds that of the projected area of the particle when horizontal (i.e. the ratio in Fig. 8b slightly exceeds 1): this can only be achieved through the influence of the finite thickness of the particle.
3.2.5 Motiontype parameter, ϵ
Since these four cases are not discrete classes of behaviour, we propose a parameter to characterise where each experiment lies on the continuum of motions between zigzag and spiral. To quantify the spectrum of periodic behaviour, a motiontype parameter, ϵ, is defined as
such that particles with ϵ=1 correspond to spiralling, as θ_{tilt}≫0 and θ_{amp}≈0 for spiralling cases. Zigzagging behaviour corresponds to ϵ=0, as zigzagging behaviour has high amplitudes, θ_{amp}≫0, and θ_{tilt}≈0. Transitional cases can have nonzero θ_{tilt} and θ_{amp}, and ϵ can therefore range between 0 and 1. The four case studies presented have ϵ= 0.01, 0.00, 0.73, and 0.94 for zigzagging, zigzag–transitional, transitional, and spiral respectively. The parameter ϵ characterises the distribution in θ and does not account for variation in ϕ. This explains why the transitional–zigzag motion case (Fig. 4b) has nearzero ϵ – the distribution shape of θ in that example is very close to that of a pure zigzag motion, even though there is also a weak azimuthal rotation superimposed which distinguishes it from the zigzag case. Although ϵ does not capture azimuthal variations, these details are often not practically significant when considering the statistics of a crystal population in a cloud; instead, it is the distribution of inclination angle θ which is of primary interest, and this motivates our definition of ϵ.
3.2.6 Oscillation frequencies
For bulk approximations (retrievals, microphysics schemes), it is useful to characterise the frequency of oscillatory behaviour. To nondimensionalise the frequency of oscillation of the experiments, the Strouhal number is calculated as
where f is the frequency of oscillation found by the sine wave fit to θ (Kajikawa, 1992). St_{θ}, the Strouhal number for θ, is therefore representative of the number of oscillations in θ of the particle in the time it takes for the particle to fall the vertical distance equal to its diameter.
St_{θ} (frequencies) of the four cases is 0.45, 0.34, 0.91, and 0.90 for zigzag, zigzag–transitional, transitional, and spiral respectively (Table 3).
A secondary Strouhal number can also be calculated using the rate of change of ϕ:
where $\left\stackrel{\mathrm{\u203e}}{\frac{\mathrm{d}\mathit{\varphi}}{\mathrm{d}t}}\right$ is the mean absolute rate of precession of the particle. St_{ϕ} is therefore the number of full turns around a vertical axis that a particle makes during the time it takes for the particle to fall the vertical distance equal to its diameter. St_{ϕ} for the zigzag and zigzag–transitional cases is 0.028 and 0.079 respectively. The zigzag case is effectively nonrotational around the vertical axis and therefore has the lowest St_{ϕ}. The transitional and spiral cases have substantial rotation around the vertical axis and therefore have St_{ϕ} of 0.57 and 0.46 respectively. The transitional case spirals faster and wobbles faster than the spiral case, but otherwise both St_{θ} and St_{ϕ} appear to increase with ϵ.
3.3 Characteristics of the full data set
3.3.1 Distributions of θ
Figure 9 displays the distributions of inclination angle across Re, excluding stable experiments, and separates the cases by particle shape and ϵ greater than or less than 0.5 to demonstrate the impact of ϵ on the distributions. The onset of unsteadiness is seen at Re ≈ 200 (Re=212 for noncircular disc shapes), the lowest onset of spiralling in particular is seen at Re =461 (for noncircular disc shapes), and the onset of spiralling for circular discs is seen at values as low as Re ≈ 300.
Different shapes show different distributions of inclination angles as well as exhibiting different θ_{tilt} and θ_{amp} (Fig. 9). For a particular shape, the distribution remains similar between adjacent Reynolds number bins; once a specific motion regime is reached for a particular shape, the same θ_{tilt} and θ_{amp} are maintained. For instance, WangBBP shapes spiral and have a mean θ of around 21° (along with a narrow distribution), while circular discs tend to have a much wider distribution, corresponding to highamplitude zigzag behaviour.
Many of the shapes that present spiralling behaviour at high Re first present zigzagging behaviour at intermediate Re (see also Figs. 3, 11, and 10).
3.3.2 Characterisation of motion
Sine waves were fit using Eq. 7 to all velocity components and to θ for all periodic experiments (Eq. 4), and the motiontype parameter ϵ was subsequently calculated (Fig. 12) (Eq. 8).
There is discrepancy between the literature on circular discs and our observations of ice crystal shapes when taking I^{*} into account. The motion parameter, ϵ, increases for increasing Reynolds number and dimensionless moment of inertia (Fig. 10), which is the opposite of the result reported by Zhong et al. (2011), who found that, for circular discs exhibiting periodic motions, spiralling occurred at lower Re and I^{*} and zigzagging occurred at higher Re and I^{*}. Our result instead agrees with Cheng et al. (2015), who found that hexagonal plates exhibit a zigzag motion at low Re, while larger plates at higher Re exhibited spiralling. Jayaweera (1965) also finds that for falling spheres, zigzagging occurs at lower Re than spiralling, which occurs at very high Re (Re > 10^{5}). Our results also deviate from the stable–periodic division line provided by Field et al. (1997), as the ice crystal shapes do not become unsteady until higher Reynolds numbers than circular discs, since shape has a strong impact on the conditions that the onset of unsteady motions occurs.
Whilst I^{*} is used successfully for studies of circular discs, shape must also be taken into account when considering the broad range of shapes that ice crystals exhibit. Differences in shape can be quantified by area ratio (the ratio of the maximum crosssectional area of the particle and the area of its circumscribing circle).
Across the experiments, ϵ increases for increasing Reynolds number and area ratio (Fig. 11. In agreement with Esteban et al. (2019) and Tagliavini et al. (2021a), stable fall behaviour is found to be much more likely for particles with lower area ratios; i.e. ice crystals with more dendritic or complex shapes are more likely to fall steadily when under the same Reynolds number.
Using linear support vector classification (Pedregosa et al., 2011), we identify a line of best fit that maximises the distance between stable and periodic behaviour, such that periodic behaviour occurs when
The critical point of this expression is displayed in Fig. 11.
Mean inclination angle was not found to distinguish well between periodic motion types (Fig. 12a). θ_{tilt} is found to be nearzero for zigzagging cases (where ϵ < 0.5), but θ_{amp} can reach up to 45° (Fig. 12b). For all experiments where ϵ > 0.5, θ_{tilt} is between 7 and 28° (Fig. 12c).
The distribution of ϵ was found to be bimodal (Fig. 12d), favouring either zigzag or spiralling behaviour, with transitional motion being less likely. Out of all 131 periodic platelike observed experiments, 65 experiments displayed ϵ < 0.2, 34 were found to have ϵ between 0.2 and 0.8, and 32 experiments had ϵ > 0.8. The potential cause of this is discussed in Sect. 3.3.4.
3.3.3 Velocity fluctuations
Across all periodic experiments, the amplitude of the vertical velocity, V_{amp}, was found to be approximately 85 % of the amplitude of the horizontal speed, U_{amp} (Fig. 13a), using leastsquares linear regression. In contrast to our findings, Kajikawa (1992) reported that the standard deviation of the horizontal velocity, U, was considerably larger (5 %–20 % of the fall velocity) than the standard deviation of the vertical velocity (< 3 % of the fall velocity) for dendriticshaped particles undergoing periodic oscillation. The reason for the difference between our findings and those of Kajikawa is unknown, and it is hard to understand why fluttering particles would have large horizontal velocity fluctuations but almost constant vertical velocity. More investigation of natural particles using modern observations, such as by Maahn et al. (2024), may help explore this in the future.
For zigzagging particles, as θ_{amp} increases, the amplitudes of both the vertical and horizontal speed components increase exponentially (Fig. 13b). For spiralling particles, as θ_{tilt} increases, the amplitude of vertical velocity increases slightly (i.e. there is more wobble). θ_{tilt} has no influence on the amplitude of vertical velocity for zigzagging particles, as it is nearzero (Fig. 13c). Other relationships between the variables mentioned in this study were also explored; however, no clear patterns or simple relationships were evident.
3.3.4 Strouhal numbers
Strouhal numbers were calculated for all experiments, as detailed in Eqs. (9) and (10). For experiments where ϵ≤0.2, St_{θ} has a mean of 0.29 and a standard deviation of 0.10, with no significant trend with Re (Fig. 14). Zigzagging particles (ϵ≤0.2) never have St_{θ} above 0.50. Despite having typically much smaller amplitudes than zigzagging particles, spiralling experiments (ϵ≥0.8) have a larger range of St_{θ}, with a potential for St_{θ} up to 1.50 and a mean St_{θ} of 0.5. St_{θ} for spiralling particles may be greater than for zigzagging; the wobbling of a spiralling particle is at a much smaller amplitude (θ_{amp}) than that of a zigzagging particle.
Similar to St_{θ} and previous literature on circular discs, there is no particular trend in St_{ϕ} with Reynolds number (Fig. 15). St_{ϕ} is close to zero for zigzagging particles (where ϵ≤0.2), as the rate of spiralling is very low (by definition), whereas for spiralling particles (ϵ≥0.8), St_{ϕ} can be as high as 0.7. No systematic relationship was found between either Strouhal number (St_{θ} or St_{ϕ}) and area ratio.
When spiralling analogues rotate faster, they tend to also wobble more frequently. For high ϵ (and correspondingly, small θ_{amp}), St_{θ} is approximately half of St_{ϕ}, corresponding to the classic observation that wobbling plates are found to wobble twice as fast as they rotate (Fig. 16) (Tuleja et al., 2007). When θ_{amp} is nonzero, the spiral motion that the centre of mass of the analogue makes is not a perfect circle: the smaller the angle of wobble, the closer the traces of the analogue are to circles. When the spin is not a perfect circle, St_{ϕ} no longer matches double the wobble rate, St_{θ}. For higher θ_{amp} (low ϵ) cases, St_{ϕ} and St_{θ} appear to both remain low and not depend on one another.
3.3.5 Mean inclination angle
For particles that are already spiralling, the higher θ_{amp} is, the more likely the behaviour is to be transitional (nonperfect spiralling) and the larger the swing the particle makes; therefore the wobble is less frequent and St_{ϕ} is lower. This may also be the cause of the bimodal distribution in ϵ (Fig. 12). Particles that have high St_{ϕ} spiral quickly relative to their vertical velocity. This fast rotation around the vertical means that any torque at 90° to the vertical axis of rotation (which causes wobble in θ) will have less effect because it is small relative to the gyroscopic torque. Therefore, the rotation of spiralling particles most likely inhibits any potential zigzagging motion. Across a set of experiments for a given particle, at low Reynolds numbers, mean inclination angle, $\stackrel{\mathrm{\u203e}}{\mathit{\theta}}$, is close to zero, as particles are stable, and at higher Reynolds numbers, particles become unstable and $\stackrel{\mathrm{\u203e}}{\mathit{\theta}}$ increases to a steady, nonzero value. To quantify and compare the onset of unstable motions for different particle properties, a logistic curve is fit to the data using a leastsquares method and the Trust Region Reflective algorithm from the SciPy curve fit function (Virtanen et al., 2020) for each particle shape (Fig. 17a) and area ratio (Fig. 17b) as follows:
where Re_{onset} is the value of the function's midpoint, ${\stackrel{\mathrm{\u203e}}{\mathit{\theta}}}_{\mathrm{unstable}}+{\stackrel{\mathrm{\u203e}}{\mathit{\theta}}}_{\mathrm{stable}}$ is the supremum of the values of the function, k is the steepness of the curve, and b is the minimum $\stackrel{\mathrm{\u203e}}{\mathit{\theta}}$ value of the function. Particles with higher area ratio typically have bigger oscillations with larger $\stackrel{\mathrm{\u203e}}{\mathit{\theta}}$. Dendrites with low area ratios and stellar crystals are stable at very high Re and typically have a smaller $\stackrel{\mathrm{\u203e}}{\mathit{\theta}}$ when unstable.
Ten different platelike snowflake shapes, in addition to circular discs, of up to three different aspect ratios each were allowed to free fall through a tank of water–glycerine mixture to simulate behaviours of real ice crystals in the atmosphere. The fall behaviour of these analogues was viewed by three orthogonal cameras, allowing for the digital reconstruction of their trajectories and orientations.
4.1 Summary
Four main falling regimes are observed: stable, zigzag, transitional, and spiralling. Stable motion has no measurable fluctuations, while other regimes involve periodic oscillations in both inclination angle and velocities. All unstable motions for the experimental series observed were of periodic behaviour: no tumbling behaviour is observed in this work. Stable analogues all had θ=0; i.e. their maximum dimension was in the horizontal plane. Zigzag motion involves swinging back and forth, while spiralling remains inclined at a constant, nonzero inclination angle.
Spiralling particles rotate steadily around the vertical axis at constant $\mathrm{d}\mathit{\varphi}/\mathrm{d}t$. That is, the rotation of spiralling planar particles does not result from a rotation around the c axis of the particle (see Fig. 1); rather, periodic rotations of equal amplitude, but phase difference of 90°, occur about the a axis and a^{′} axis such that θ is approximately constant. This rotation of the particle about the a and a^{′} axes causes the particle's centre of mass to trace a circular path. By contrast, zigzagging particles maintain a constant azimuth angle, ϕ, that has a square wave (such that the minima and maxima are spaced 180° apart). Transitional cases are a mix of the two behaviours: they swing back and forth but also rotate as they do so. Particle components of velocity (V and U for vertical and horizontal respectively) were also found to be sinusoidal with respect to time. Sine waves are fit to time series of θ and U and V, and the rate of spiralling, $\mathrm{d}\mathit{\varphi}/\mathrm{d}t$, is found through linear regression. The amplitude of the sine wave, θ_{amp}, was found to vary between 0° and 43.1°. Periodic motion is found to be analogous to the range of spherical pendulum behaviour, corresponding to simple harmonic motion. Time series of θ and velocities for periodic experiments are therefore sinusoidal, and distributions of θ have a nonzero mode. The results do not support the common assumption of Gaussian orientation distributions with a zeromodal angle during unstable motions: the distributions of θ have nonGaussian distributions and nonzero modes.
In the spiralling regime, components of velocity are held relatively constant compared to the other cases, effectively making the spiralling cases a quasisteady mode with a nonzero nearconstant inclination. When particles spiral, they are consistently inclined at an angle, observed to typically be between 7 and 28°. The central line of the sine wave fit, θ_{tilt}, is typically between 8 and 25° for spiralling behaviour for all particles, with a mean of 18.4 ° and a standard deviation of 6.8°.
Strouhal numbers (nondimensionalised frequencies) were found using the rate of spiralling, $\mathrm{d}\mathit{\varphi}/\mathrm{d}t$, and the frequency of the sine waves of θ, finding St_{ϕ} and St_{θ} respectively. These each represent the number of turns the particle makes and the number of wobbles the particle makes in the time taken for the particle to fall the vertical distance equal to its own diameter.
St_{ϕ} was found to be approximately half of St_{θ} for spiralling cases: particles were found to wobble twice as often as they made a full rotation around the vertical, consistent with previous work (Tuleja et al., 2007). For zigzagging cases, St_{θ} is found to be 0.29 ± 0.7, with no variation with Reynolds number.
The onset of unstable motions is found to be more likely for higher area ratios corresponding to less complex shapes (such as pristine hexagonal plates), in agreement with Esteban et al. (2018) and Tagliavini et al. (2021a). The shapes D1 (at all aspect ratios), DP, S, and F (at aspect ratio 0.04) remained stable throughout all experiments, even at Re > 10^{3}. The onset for circular discs was found to be as low as Re=197.
A motiontype parameter ϵ is calculated using θ_{amp} and θ_{tilt} to quantify the spectrum of behaviour from zigzag (ϵ=0) to spiral (ϵ=1). ϵ increases for increasing Reynolds number and dimensionless moment of inertia: spiralling is more likely when both parameters are higher. Particles were observed to exhibit zigzagging behaviour at lower Reynolds number than spiralling, and some particles (WangBBP) were found to exclusively spiral when unsteady. This contrasts with findings from Esteban et al. (2018) and Zhong et al. (2013), who expect spiralling to occur at lower Re and I^{*} than zigzagging.
4.2 Implication of results
The literature on ice crystal orientations often assumes that they can be modelled by a Gaussian distribution of the inclination angle with a mode at θ=0° and with a breadth that is independent of the particle size.
Our findings show that in fact we should expect the distribution of inclination to be nonGaussian, with a mode close to the maximum inclination that the particle experiences. For zigzag motions the distributions may be rather broad, spanning a few tens of degrees (Fig. 3a). For spiralling particles the distributions are very narrow (only a few degrees) but with a substantial systematic inclination (Fig. 3d).
It is also evident from our data that the distribution of inclination angle varies sharply around some critical Reynolds number. Small particles fall steadily with horizontal orientation. Large particles fall unsteadily with a substantial inclination on average. The data in Fig. 17 suggest that the transition between these two modes of fall is relatively sharp and is dependent on the shapes of the particles, with open shapes like stellar crystals and dendrites falling stably at higher Reynolds number (larger diameter) than hexagonal plates and broadbranched crystal forms.
Groundbased snowflake imagers have reported preferentially nonhorizontal orientations, and the laboratory results reported here may provide a means to understand that observation. In future work, we hope to make a more detailed comparison between field observations of fluttering crystals Maahn et al. (2024) and our expectations from the laboratory.
Incorporating this information into remote sensing retrievals and interpretation (e.g. polarimetric radar signatures) should help improve the accuracy and robustness of those analyses. New electromagnetic scattering databases provide increasing flexibility to integrate over arbitrary distributions of inclination angle (Brath et al., 2020). It is difficult as yet to make a simple prescription for what the distribution of θ should be as a function of crystal size and shape – as we have seen, there is significant variability in behaviour across crystals of different cross section and aspect ratio. However, Fig. 17 gives an indication of what a realistic mean inclination angle could be for various Reynolds numbers, while Fig. 9 provides more detail on the typical form of those distributions. We suggest choosing assumptions for crystal orientation distributions that are consistent with these data.
The observation that large fluctuations may occur in the fall velocity of unstable particles implies that a single mean speed cannot be used to approximate the velocity of a single fluttering crystal, and that a spread of velocities should be considered when interpreting such spectra. This is expected to appear as a broadening of the Doppler spectrum from a verticalpointing radar. An estimate of the magnitude of these fluctuations can be deduced from Fig. 13. Likewise, the fluctuating horizontal velocity of the crystals acts to broaden the Doppler spectrum for nearhorizontally scanning weather radars, which should be considered when inferring the distribution of turbulent air motions from such data, and it raises the intriguing prospect of retrieving crystal fluttering characteristics from horizontal Doppler data in conditions where turbulence and wind shear are weak or absent.
4.3 Limitations and future work
Strong turbulence, typically characterised by high kinetic eddy dissipation rates (e.g. $>\phantom{\rule{0.125em}{0ex}}\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}{\text{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\text{s}}^{\mathrm{3}}$), could affect crystal orientation when crystals are large (> 1 mm) (Klett, 1995; Garrett et al., 2015). In these cases, strong turbulence is found to widen the distribution of θ for both stable and unstable particles (Fitch et al., 2021). Our study only considers quiescent conditions, as we want to know under what conditions particle instability still occurs, even without the addition of turbulence.
Although strong turbulence is typical within convective clouds and at the ground, typical turbulenceinduced velocity perturbations across the faces of ice crystals within clouds are approximately 50 times smaller than the vertical velocity of the crystal, and other aerodynamic factors are involved (Westbrook et al., 2010). Studies of sun glints have shown that there are many cases in which turbulence does not dominate and conditions can be considered quiescent, such that ice particles have horizontal orientations (Marshak et al., 2017; Varnai et al., 2020). Turbulence does not typically dominate the fall behaviour of very small particles, as the scales of turbulence are not small enough to influence the orientation of the crystals, with only a slight wobble of up to 2° observed in some cases (Sassen, 1980; Klett, 1995).
In the literature on circular discs, Re is the key control on the onset of unsteadiness, while the parameter I^{*} is argued to modulate the form of the unsteady motion (recall Fig. 2). In our data, it is clear that I^{*} on its own is not the only relevant parameter or perhaps not even the leading control. Nevertheless, we acknowledge that I^{*} in our current experiment is significantly smaller than the case of ice crystals falling in air, largely due to the difference in density between the laboratory fluid (water) versus the atmosphere (air). To address this, we are currently undertaking a new set of experiments with much lighter analogues falling in air and will report these results in a future publication.
Many aspects of shape are not covered by area ratio, and other shape parameters could later be explored in addition to area ratio to capture the full variability of shape parameters. Future work therefore also includes exploration of the impact of shape, with the aim of understanding the influence of experimental conditions on unsteadiness more accurately than the results presented in this study.
Understanding the fall behaviour of ice crystals allows us to further understand the speed at which they grow, fall, and precipitate, allowing this behaviour to be modelled and parameterised more effectively. Further research exploring an even wider range of ice crystal shapes, sizes, and environmental conditions will help build on these findings and advance our overall understanding of ice crystal dynamics within the complex atmospheric system.
A1 Estimation of the magnitude of I^{*} in the atmosphere
We estimated the magnitude of I^{*} for real ice crystals falling in air by using a mass–diameter relationship from Nakaya and Terada (1935) for planar dendritic crystals:
where m is in milligrams and d is in millimetres. A predicted value for I^{*} can then be calculated following the method in Kajikawa (1992), such that
where I_{a} is the moment of inertia about the a axis of the crystal and M and D are the mass and diameter in kilograms and metres respectively. This can then be used to calculate I^{*} as
where we set the density of air to be 1.2 kg m^{−3}. We find that 10, 1, and 0.1 mm planar dendritic crystals have an I^{*} of 0.02, 0.2, and 2.0 respectively.
The mass–diameter relationship from Nakaya and Terada (1935) was chosen as an example to illustrate the order of magnitude of I^{*} that could be expected in the atmospheric case (we are not attempting to provide precise estimates of I^{*} for specific crystal shapes and dimensions). The approximate formula for I^{*} from Kajikawa (1992) was selected because it requires only the mass and diameter of the crystal as inputs, without detailed knowledge of the full mass distribution around the snowflake.
Data for the case studies and the characteristics of the full data set are available as a Supplement.
Videos from the four case studies presented in Sect. 3.2 are given in the Supplement. The supplement related to this article is available online at: https://doi.org/10.5194/acp24111332024supplement.
JRS contributed to conceptualisation, data curation, software, formal analysis of the work, and wrote the manuscript draft. CDW and THMS provided supervision, project administration, and helped write the original manuscript draft. MWM produced the data set and provided advice and data curation. All authors provided conceptualisation and review and editing of the manuscript.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This work was conducted with the support of a PhD studentship through the SCENARIO Doctoral Training Partnership, funded by NERC, project code F4114950. We would like to thank our reviewers and editor Ann Fridlind for their constructive suggestions, which improved the quality of this paper.
This work was conducted with the support of a PhD studentship through the SCENARIO Doctoral Training Partnership (grant no. NE/S007261/1).
This paper was edited by Ann Fridlind and reviewed by two anonymous referees.
Borovoi, A. and Kustova, N.: Display of ice crystal flutter in atmospheric light pillars, Geophys. Res. Lett., 36, L04804, https://doi.org/10.1029/2008GL036413, 2009. a
Brath, M., Ekelund, R., Eriksson, P., Lemke, O., and Buehler, S. A.: Microwave and submillimeter wave scattering of oriented ice particles, Atmos. Meas. Tech., 13, 2309–2333, https://doi.org/10.5194/amt1323092020, 2020. a
Bringi, V. N. and Chandrasekar, V.: Polarimetric Doppler Weather Radar, Cambridge University Press, ISBN 9780521623841, https://doi.org/10.1017/CBO9780511541094, 2001. a
Cheng, K. Y., Wang, P. K., and Hashino, T.: A numerical study on the attitudes and aerodynamics of freely falling hexagonal ice plates, J. Atmos. Sci., 72, 3685–3698, https://doi.org/10.1175/JASD150059.1, 2015. a, b, c, d
Cho, H. R., Iribarne, V., and Richards, W. G.: On the Orientation of Ice Crystals in a Cumulonimbus Cloud, American Meteorological Society, 38, 1111–1114, 1981. a
Curry, J. A. and Ebert, E. E.: Annual Cycle of Radiation Fluxes over the Arctic Ocean: Sensitivity to Cloud Optical Properties, J. Climate, 5, 1267–1280, https://doi.org/10.1175/15200442(1992)005<1267:ACORFO>2.0.CO;2, 1992. a
Dunnavan, E. L. and Jiang, Z.: A General Method for Estimating Bulk 2D Projections of Ice Particle Shape: Theory and Applications, J. Atmos. Sci., 76, 305–332, https://doi.org/10.1175/JASD180177.1, 2019. a
Ern, P., Risso, F., Fabre, D., and Magnaudet, J.: Wakeinduced oscillatory paths of bodies freely rising or falling in fluids, Annu. Rev. Fluid Mech., 44, 97–121, https://doi.org/10.1146/annurevfluid120710101250, 2011. a
Esteban, L. B.: Dynamics of NonSpherical Particles in Turbulence, in: Springer Theses, Springer Cham, PhD thesis, https://doi.org/10.1007/9783030281366, 2019. a
Esteban, L. B., Shrimpton, J., and Ganapathisubramani, B.: Edge effects on the fluttering characteristics of freely falling planar particles, Phys. Rev. Fluids, 3, 064302, https://doi.org/10.1103/PhysRevFluids.3.064302, 2018. a, b, c, d
Esteban, L. B., Shrimpton, J. S., and Ganapathisubramani, B.: Disks settling in turbulence, J. Fluid Mech., 883, A58, https://doi.org/10.1017/jfm.2019.922, 2019. a, b
Feist, M. M., Westbrook, C. D., Clark, P. A., Stein, T. H., Lean, H. W., and Stirling, A. J.: Statistics of convective cloud turbulence from a comprehensive turbulence retrieval method for radar observations, Q. J. Roy. Meteor. Soc., 145, 727–744, https://doi.org/10.1002/qj.3462, 2019. a
Field, S. B., Klaus, M., Moore, M. G., and Norl, F.: Chaotic dynamics of falling disks, Letters to Nature, 388, 252–254, 1997. a, b, c, d, e, f, g, h, i
Fitch, K. E., Hang, C., Talaei, A., and Garrett, T. J.: Arctic observations and numerical simulations of surface wind effects on MultiAngle Snowflake Camera measurements, Atmos. Meas. Tech., 14, 1127–1142, https://doi.org/10.5194/amt1411272021, 2021. a, b, c
Garrett, T. J., Yuter, S. E., Fallgatter, C., Shkurko, K., Rhodes, S. R., and Endries, J. L.: Orientations and aspect ratios of falling snow, Geophys. Res. Lett., 42, 4617–4622, https://doi.org/10.1002/2015GL064040, 2015. a, b, c
Gibson, A., Thomas, L., and Bhattacharyya, S.: Some characteristics of cirrus clouds deduced from laserradar observations at different elevation angles, J. Atmos. Terr. Phys., 39, 657–660, https://doi.org/10.1016/00219169(77)900794, 1977. a
Gregory, R. D.: Classical Mechanics, Cambridge University Press, 2006. a
Helt, G.: A Rose By Any Other Name, Tessellations Publishing, http://archive.bridgesmathart.org/2016/bridges2016445.html (last access: 24 January 2024), 2016. a
Hirakata, M., Okamoto, H., Hagihara, Y., Hayasaka, T., and Oki, R.: Comparison of Global and Seasonal Characteristics of Cloud Phase and Horizontal Ice Plates Derived from CALIPSO with MODIS and ECMWF, J. Atmos. Ocean. Tech., 31, 2114–2130, https://doi.org/10.1175/JTECHD1300245.1, 2014. a
Hogan, R. J. and Illingworth, A. J.: The effect of specular reflection on spaceborne lidar measurements of ice clouds, Tech. rep., Department of Meteorology, University of Reading, https://www.met.reading.ac.uk/~swrhgnrj/publications/specular.pdf (last access: 21 January 2024), 2003. a
Ishimoto, H., Masuda, K., Mano, Y., Orikasa, N., and Uchiyama, A.: Irregularly shaped ice aggregates in optical modeling of convectively generated ice clouds, J. Quant. Spectrosc. Ra., 113, 632–643, https://doi.org/10.1016/j.jqsrt.2012.01.017, 2012. a
Jayaweera, K.: The behaviour of small clusters of bodies falling in a viscous fluid, PhD thesis, Imperial College London, https://core.ac.uk/outputs/1590667/ (last access: 25 August 2022), 1965. a, b
Jiang, Z., Oue, M., Verlinde, J., Clothiaux, E. E., Aydin, K., Botta, G., and Lu, Y.: What can we conclude about the real aspect ratios of ice particle aggregates from twodimensional images?, J. Appl. Meteorol. Clim., 56, 725–734, https://doi.org/10.1175/JAMCD160248.1, 2017. a
Kajikawa, M.: Observations of the Falling Motion of PlateLike Snow Crystals Part I: The FreeFall Patterns and Velocity Variations of Unrimed Crystals, J. Meteorol. Soc. Jpn., 70, 1–9, https://doi.org/10.2151/jmsj1965.70.1_1, 1992. a, b, c, d, e, f, g, h, i, j, k
Klett, J. D.: Orientation Model for Particles in Turbulence, J. Atmos. Sci, 52, 2276–2285, 1995. a, b, c
Köbschall, K., Breitenbach, J., Roisman, I. V., Tropea, C., and Hussong, J.: Geometric descriptors for the prediction of snowflake drag, Exp. Fluids, 64, 4, https://doi.org/10.1007/s0034802203539x, 2023. a, b
List, R. and Schemenauer, R. S.: FreeFall Behavior of Planar Snow Crystals, Conical Graupel and Small Hail, J. Atmos. Sci., 28, 110115, https://doi.org/10.1175/15200469(1971)028<0110:FFBOPS>2.0.CO;2, 1971. a, b, c
Locatelli, J. D. and Hobbs, P. V.: Fall speeds and masses of solid precipitation particles, J. Geophys. Res., 79, 2185–2197, https://doi.org/10.1029/jc079i015p02185, 1974. a
Lynch, D. K., Gedzelman, S. D., and Fraser, A. B.: Subsuns, Bottlinger's rings, and elliptical halos, Appl. Optics, 33, 4580–4589, 1994. a, b
Maahn, M., Moisseev, D., Steinke, I., Maherndl, N., and Shupe, M. D.: Introducing the Video In Situ Snowfall Sensor (VISSS), Atmos. Meas. Tech., 17, 899–919, https://doi.org/10.5194/amt178992024, 2024. a, b
Magono, C. and Lee, C. W.: Meteorological Classification of Natural Snow Crystals, Journal of the Faculty of Science, 2, 321–335, 1966. a
Marshak, A., Várnai, T., and Kostinski, A.: Terrestrial glint seen from deep space: Oriented ice crystals detected from the Lagrangian point, Geophys. Res. Lett., 44, 5197–5202, https://doi.org/10.1002/2017GL073248, 2017. a
Masuda, K. and Ishimoto, H.: Influence of particle orientation on retrieving cirrus cloud properties by use of total and polarized reflectances from satellite measurements, J. Quant. Spectrosc. Ra., 85, 183–193, https://doi.org/10.1016/S00224073(03)002243, 2004. a
Matrosov, S. Y.: Feasibility of using radar differential Doppler velocity and dualfrequency ratio for sizing particles in thick ice clouds, J. Geophys. Res.Atmos., 116, D17202, https://doi.org/10.1029/2011JD015857, 2011. a
McCorquodale, M. W. and Westbrook, C. D.: TRAIL: A novel approach for studying the aerodynamics of ice particles, Q. J. Roy. Meteor. Soc., 147, 589–604, https://doi.org/10.1002/qj.3935, 2021a. a, b, c, d, e, f, g
McCorquodale, M. W. and Westbrook, C. D.: TRAIL part 2: A comprehensive assessment of ice particle fall speed parametrisations, Q. J. Roy. Meteor. Soc., 147, 605–626, https://doi.org/10.1002/qj.3936, 2021b. a, b, c, d, e
McCorquodale, M. W. and Westbrook, C. D.: The stable orientation of aggregate ice particles in freefall, 3rd International Summer Snowfall Workshop (ISSW3), 31 August–3 September 2021, virtual conference, https://www.met.reading.ac.uk/~sws04cdw/issw3_programme.html (last access: 31 September 2021), 2021c. a
Melnikov, V. and Straka, J. M.: Axis ratios and flutter angles of cloud ice particles: Retrievals from radar data, J. Atmos. Ocean. Tech., 30, 1691–1703, https://doi.org/10.1175/JTECHD1200212.1, 2013. a, b
Moilanen, J. and Gritsevich, M.: Light scattering by airborne ice crystals – An inventory of atmospheric halos, J. Quant. Spectrosc. Ra., 290, 108313, https://doi.org/10.1016/j.jqsrt.2022.108313, 2022. a
Nakaya, U. and Terada, T.: Simultaneous Observations of the Mass, Falling Velocity and Form of Individual Snow Crystals, Hokkaido Imperial University, 1, 191–200, http://hdl.handle.net/2115/34452 (last access: 2 June 2023), 1935. a, b, c
Nettesheim, J. J. and Wang, P. K.: A numerical study on the aerodynamics of freely falling planar Ice Crystals, J. Atmos. Sci., 75, 2849–2865, https://doi.org/10.1175/JASD180041.1, 2018. a, b, c, d
Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Duchesnay, E.: Scikitlearn: Machine Learning in {P}ython, J. Mach. Learn. Res., 12, 2825–2830, 2011. a
Platt, C. M. R.: Lidar Observation of a MixedPhase Altostratus Cloud, J. Appl. Meteorol., 16, 339–345, https://doi.org/10.1175/15200450(1977)016<0339:LOOAMP>2.0.CO;2, 1977. a, b
Pruppacher, H. and Klett, J.: Microphysics of Clouds and Precipitation, Springer Netherlands, ISBN 9780792342113, https://doi.org/10.1007/9780306481000, 1997. a
Reiter, C. A.: A local cellular model for snow crystal growth, Chaos, Solitons and Fractals, 23, 1111–1119, https://doi.org/10.1016/j.chaos.2004.06.071, 2005. a
Ryzhkov, A. V., Snyder, J., Carlin, J. T., Khain, A., and Pinsky, M.: What polarimetric weather radars offer to cloud modelers: Forward radar operators and microphysical/thermodynamic retrievals, Atmosphere, 11, 362, https://doi.org/10.3390/atmos11040362, 2020. a
Sassen, K.: Ice Crystal Habit Discrimination with the Optical Backscatter Depolarization Technique, J. Appl. Meteorol., 16, 425–431, https://doi.org/10.1175/15200450(1977)016<0425:ICHDWT>2.0.CO;2, 1977. a
Sassen, K.: Remote Sensing of Planar Ice Crystal Fall Attitudes, J. Meteorol. Soc. Jpn., 58, 422–429, 1980. a
Sassen, K.: Polarization and Brewster angle properties of light pillars, J. Opt. Soc. Am. A, 4, 570–580, https://doi.org/10.1364/JOSAA.4.000570, 1987. a
Schrom, R. S., Munchak, S. J., and Adams, I. S.: The Polarimetric Radar Scattering Properties of Oriented Aggregates, J. Atmos. Sci., 80, 1145–1161, https://doi.org/10.1175/JASD220149.1, 2023. a
Stillwell, R. A., Neely, R. R., Thayer, J. P., Walden, V. P., Shupe, M. D., and Miller, N. B.: Radiative Influence of Horizontally Oriented Ice Crystals over Summit, Greenland, J. Geophys. Res.Atmos., 124, 12141–12156, https://doi.org/10.1029/2018JD028963, 2019. a
Tagliavini, G., McCorquodale, M., Westbrook, C., Corso, P., Krol, Q., and Holzner, M.: Drag coefficient prediction of complexshaped snow particles falling in air beyond the Stokes regime, Int. J. Multiphas. Flow, 140, 103652, https://doi.org/10.1016/j.ijmultiphaseflow.2021.103652, 2021a. a, b, c
Tagliavini, G., McCorquodale, M., Westbrook, C., and Holzner, M.: Numerical analysis of the wake of complexshaped snow particles at moderate Reynolds number, Phys. Fluids, 33, 105103, https://doi.org/10.1063/5.0064902, 2021b. a
Tagliavini, G., Khan, M. H., McCorquodale, M., Westbrook, C., and Holzner, M.: Wake characteristics of complexshaped snow particles: Comparison of numerical simulations with fixed snowflakes to timeresolved particle tracking velocimetry experiments with freefalling analogs, Phys. Fluids, 34, 055112, https://doi.org/10.1063/5.0089759, 2022. a
Takano, Y. and Liou, K.N.: Solar Radiative Transfer in Cirrus Clouds. Part I: SingleScattering and Optical Properties of Hexagonal Ice Crystals, J. Atmos. Sci., 46, 3–19, https://doi.org/10.1175/15200469(1989)046<0003:SRTICC>2.0.CO;2, 1989. a
Tränkle, E. and Riikonen, M.: Elliptical halos, Bottlinger’s rings, and the iceplate snowstar transition, Appl. Optics, 35, 4871, https://doi.org/10.1364/AO.35.004871, 1996. a
Tuleja, S., Gazovic, B., Tomori, A., and Hanc, J.: Feynman’s wobbling plate, Am. J. Phys., 75, 240–244, https://doi.org/10.1119/1.2402156, 2007. a, b
Varnai, T., Kostinski, A. B., and Marshak, A.: Deep Space Observations of Sun Glints from Marine Ice Clouds, IEEE Geosci. Remote S., 17, 735–739, https://doi.org/10.1109/LGRS.2019.2930866, 2020. a
Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., İlhan Polat, Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., van Mulbregt, P., Vijaykumar, A., Bardelli, A. P., Rothberg, A., Hilboll, A., Kloeckner, A., Scopatz, A., Lee, A., Rokem, A., Woods, C. N., Fulton, C., Masson, C., Häggström, C., Fitzgerald, C., Nicholson, D. A., Hagen, D. R., Pasechnik, D. V., Olivetti, E., Martin, E., Wieser, E., Silva, F., Lenders, F., Wilhelm, F., Young, G., Price, G. A., Ingold, G.L., Allen, G. E., Lee, G. R., Audren, H., Probst, I., Dietrich, J. P., Silterra, J., Webber, J. T., Slavič, J., Nothman, J., Buchner, J., Kulick, J., Schönberger, J. L., de Miranda Cardoso, J. V., Reimer, J., Harrington, J., Rodríguez, J. L. C., NunezIglesias, J., Kuczynski, J., Tritz, K., Thoma, M., Newville, M., Kümmerer, M., Bolingbroke, M., Tartre, M., Pak, M., Smith, N. J., Nowaczyk, N., Shebanov, N., Pavlyk, O., Brodtkorb, P. A., Lee, P., McGibbon, R. T., Feldbauer, R., Lewis, S., Tygier, S., Sievert, S., Vigna, S., Peterson, S., More, S., Pudlik, T., Oshima, T., Pingel, T. J., Robitaille, T. P., Spura, T., Jones, T. R., Cera, T., Leslie, T., Zito, T., Krauss, T., Upadhyay, U., Halchenko, Y. O., and VázquezBaeza, Y.: SciPy 1.0: fundamental algorithms for scientific computing in Python, Nat. Methods, 17, 261–272, https://doi.org/10.1038/s4159201906862, 2020. a, b
von Lerber, A., Moisseev, D., Bliven, L. F., Petersen, W., Harri, A.M., and Chandrasekar, V.: Microphysical Properties of Snow and Their Link to Ze–S Relations during BAECC 2014, J. Appl. Meteorol. Clim., 56, 1561–1582, https://doi.org/10.1175/JAMCD160379.1, 2017. a
Wang, P. K.: Motions of Ice Hydrometeors in the Atmosphere Numerical Studies and Implications, Springer, 1 edn., ISBN 9789813344310, https://doi.org/10.1007/9789813344310, 2021. a, b
Westbrook, C. D. and Sephton, E. K.: Using 3Dprinted analogues to investigate the fall speeds and orientations of complex ice particles, Geophy. Res. Lett., 44, 7994–8001, https://doi.org/10.1002/2017GL074130, 2017. a, b
Westbrook, C. D., Illingworth, A. J., O'Connor, E. J., and Hogan, R. J.: Doppler lidar measurements of oriented planar ice crystals falling from supercooled and glaciated layer clouds, Q. J. Roy. Meteor. Soc., 136, 260–276, https://doi.org/10.1002/qj.528, 2010. a, b
Willmarth, W. W., Hawk, N. E., and Harvey, R. L.: Steady and Unsteady Motions and Wakes of Freely Falling Disks, Phys. Fluids, 7, 197–208, 1964. a, b, c, d
Zhong, H., Chen, S., and Lee, C.: Experimental study of freely falling thin disks: Transition from planar zigzag to spiral, Phys. Fluids, 23, 011702, https://doi.org/10.1063/1.3541844, 2011. a, b, c, d, e, f, g, h
Zhong, H., Lee, C., Su, Z., Chen, S., Zhou, M., and Wu, J.: Experimental investigation of freely falling thin disks. Part 1. the flow structures and Reynolds number effects on the zigzag motion, J. Fluid Mech., 716, 228–250, https://doi.org/10.1017/jfm.2012.543, 2013. a, b, c
Zhou, T., Huang, J., Huang, Z., Liu, J., Wang, W., and Lin, L.: The depolarization–attenuated backscatter relationship for dust plumes, Opt. Express, 21, 15195, https://doi.org/10.1364/OE.21.015195, 2013. a
Zikmunda, J.: Fall velocities of Spatial Crystals and Aggregates, J. Atmos. Sci., 29, 1511–1515, 1972. a