Stable and unstable fall motions of plate-like ice crystal analogues

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 3D-printed plate-like ice crystal analogues in a tank of water–glycerine mixture are observed with multi-view 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 back-and-forth 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 non-zero 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.

1 Introduction

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.

Figure 1 Axes of a crystal inclined by θ of 20° and pointing towards an azimuth, ϕ, of 45°, both relative to the laboratory frame of reference, given by a (x,z) horizontal plane and a vertical y axis. The crystal plane is represented by the a and a^′ axis, where the a axis is aligned with one of the crystal branches. The c axis is perpendicular to the crystal plane. The views provided are for an observer facing (a) the c–y plane (b), the a^′–y plane, and (c) the x–z plane.

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 dual-polarisation 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 radar-derived 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 ground-based 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 sub-sections.

1.2.1 Reynolds number

The Reynolds number is defined as

$$\begin{array}{}\text{(1)}& \mathit{Re}={\displaystyle \frac{V_{\mathrm{mean}}D}{\mathit{\upsilon }}},\end{array}$$

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 sub-types of behaviour, and more recently an in-between 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).

Figure 2 Phase diagram showing the stable (black) and unstable (purple) behaviour of falling particles as a function of I^* (dimensionless moment of inertia) and Re (Reynolds number). Data from TRAIL (this study and McCorquodale and Westbrook (2021a, b)) are in solid circles, and all other data points are for shapes relevant to ice crystals (Nettesheim and Wang, 2018; Kajikawa, 1992; Cheng et al., 2015). Solid lines and annotations are from Field et al. (1997) and dashed lines and rotated annotations are from Zhong et al. (2011), presenting the observed behaviour for circular discs. Acronyms in the legend refer to the shapes in Table 1. The hatched region is the expected range of I^* for dendritic planar ice crystals.

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 non-zero, 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

$$\begin{array}{}\text{(2)}& I_{\mathrm{disc}}^{*}={\displaystyle \frac{\mathit{\pi }}{\mathrm{64}}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{p}}}{{\mathit{\rho }}_{\mathrm{f}}}}{\displaystyle \frac{t}{D}},\end{array}$$

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 non-dimensional moment of inertia is

$$\begin{array}{}\text{(3)}& I^{*}={\displaystyle \frac{I_a}{{\mathit{\rho }}_{\mathrm{f}}{D}^{\mathrm{5}}}},\end{array}$$

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 a-axis 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⁻³ and the density of air is 1.2 kg m⁻³, 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 cross-sectional 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 Cloud-Aerosol 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 mixed-phase cloud layers. Similarly, Stillwell et al. (2019) found that horizontally oriented plates must occur in at least 25.6 % of all ice-only 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 zero-inclination 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, size-independent 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 (fern-like 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.

Multi-Angle 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 near-horizontal orientation but have non-zero modal values. Recent research into the MASC measurements by Fitch et al. (2021) has also reported preferential non-horizontal 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 non-horizontal, 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 pendulum-like 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 broad-branched 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 spin-up 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, multi-degree-of-freedom 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 – scaled-up 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 broad-branched models, small oscillations are observed at Re ≈ 200, although no oscillations are observed at this Reynolds number for stellar, dendritic, or stellar-with-plate 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 free-falling 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 plate-like 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.

2 Data

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.

Table 1 Particle shapes analysed in this study. Area ratio is calculated using the observed projected area divided by the circumscribing circle around the maximum diameter, as seen from beneath when fall motion is steady. Re is calculated using the observed mean velocity and maximum observed diameter as seen from beneath.

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 non-hexagonally 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 plate-like 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.

3 Results

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 sub-types – 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 near-zero inclination angle), with no measurable fluctuations in velocity and no horizontal movements.

Figure 3 Motion-type coverage of Reynolds number for each shape and aspect ratio. Stable, zigzag, transitional, and spiral are black crosses, pink circles, black diamonds, and blue triangles respectively. Particle shape labelling is defined in Table 1.

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⁻³–12 × 10⁻³ for stable motion and 0.28 × 10⁻³–11 × 10⁻³ 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³. 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 sub-types 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.

Figure 4 Case studies visualising the periodic motion sub-types. Side views of the particle motion (a, c, e, g) and linearly detrended centre of mass normalised by particle diameter (b, d, f, h) coloured by inclination angle θ (°) for the zigzag, zigzag–transitional, transitional, and spiral cases respectively.

Table 2 The observed behaviours of the presented case studies.

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 broad-branched 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 sub-type of periodic motion remains at a near-constant inclination angle and does not swing back and forth and instead precesses around its central point without touching its mean centre of fall.

The sub-types of periodic motion can be approximated by the sub-types 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

$$\begin{array}{}\text{(4)}& \mathit{\theta }=\left|{\mathit{\theta }}_{\mathrm{amp}}\sin \right(\mathit{\omega }t+{\mathit{\omega }}_{\mathrm{0}})+{\mathit{\theta }}_{\mathrm{tilt}}|,\end{array}$$

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 ω₀ is the phase shift.

Figure 5 Time series of inclination angles (°) for the periodic motion sub-types shown in Fig. 4 (black) and their fitted curves from Eq. (4) (blue). A histogram of angles is shown to the right of each panel, using 2.5° bins.

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 non-sinusoidal 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 single-frequency 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 sub-types 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 near-zero θ_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 non-zero θ_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).

Figure 6 Time series of azimuth angles, ϕ (°), for the periodic motion sub-types.

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 saw-tooth 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 square-wave 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 non-hexagonally symmetric (S, F, D, DP, and PB), the shortest branches are the axis the crystal pivots around.

Transitional cases have a non-constant $\mathrm{d}\mathit{\varphi }/\mathrm{d}t$, a combination of the saw-tooth 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:

$$\begin{array}{}\text{(5)}& U=\sqrt{u^{\mathrm{2}}+w^{\mathrm{2}}}.\end{array}$$

Sine waves can be fit to the vertical component of velocity, V, and the horizontal component of velocity, U, such that

$$\begin{array}{}\text{(6)}& {\displaystyle }V=V_{\mathrm{amp}}\sin ({\mathit{\omega }}_{V}t+{\mathit{\omega }}_{V\mathrm{0}})+V_{\mathrm{mean}},\text{(7)}& {\displaystyle }U=\left|U_{\mathrm{amp}}\sin \right({\mathit{\omega }}_{U}t+{\mathit{\omega }}_{U\mathrm{0}})+U_{\mathrm{offset}}|.\end{array}$$

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).

Figure 7 Time series of the horizontal and vertical velocity components for the periodic motion sub-types, normalised by the mean vertical velocity in each experiment.

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 non-zero 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.

Table 3 The experimental conditions for the presented case studies from Figs. 4–8.

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 quasi-steady mode with a non-zero near-constant inclination.

For the non-spiralling 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 non-spiralling 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 ground-based (usually close to zenith) or spaceborne (usually close to nadir). Variation in the projected area affects assumptions in backscatter cross-section 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 cross-sectional area of each analogue. The projected area as seen from below anti-correlates with θ in each time series, displaying an out-of-phase relationship: as θ increases, the aspect ratio of the analogues increases (becoming closer to 1), while the area ratio and projected area decrease.

Figure 8 Time series of reconstructed projected areas as seen from below for each case study, normalised by the observed projected area when θ=0 (solid blue line). θ (dashed black line) and cos θ (dashed light blue line) are provided for reference.

For all four cases, vertical velocity and projected area exhibit a 180° out-of-phase 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 cross-sectional 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 Motion-type 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 motion-type parameter, ϵ, is defined as

$$\begin{array}{}\text{(8)}& \mathit{ϵ}={\displaystyle \frac{{\mathit{\theta }}_{\mathrm{tilt}}}{{\mathit{\theta }}_{\mathrm{tilt}}+{\mathit{\theta }}_{\mathrm{amp}}}},\end{array}$$

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 non-zero θ_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 near-zero ϵ – 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 non-dimensionalise the frequency of oscillation of the experiments, the Strouhal number is calculated as

$$\begin{array}{}\text{(9)}& {\mathit{St}}_{\mathit{\theta }}={\displaystyle \frac{fD}{V_{\mathrm{mean}}}},\end{array}$$

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 ϕ:

$$\begin{array}{}\text{(10)}& {\mathit{St}}_{\mathit{\varphi }}={\displaystyle \frac{\mathrm{1}}{\mathrm{360}\mathit{°}}}\left|\stackrel{\mathrm{‾}}{{\displaystyle \frac{d\mathit{\varphi }}{dt}}}\right|{\displaystyle \frac{D}{V_{\mathrm{mean}}}},\end{array}$$

where $\left|\stackrel{\mathrm{‾}}{\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 non-rotational 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 non-circular disc shapes), the lowest onset of spiralling in particular is seen at Re =461 (for non-circular disc shapes), and the onset of spiralling for circular discs is seen at values as low as Re ≈ 300.

Figure 9 Inclination angle distributions of all unstable experiments, equally weighted by experiment, binned by Reynolds number and particle shape, with shapes in order of decreasing area ratio, excluding CD-P and D1. Experiments are split by spiral (ϵ>0.5; upper, orange) or zigzag (ϵ≤0.5; lower, blue). Quartiles (dashed) and mean values (solid) are inside each distribution. ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{amp}}$ ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{tilt}}$ are the mean values of θ_amp and θ_tilt for each set of experiments (separated by ϵ and particle shape). Reported values are the mean (±standard deviation of the mean) for each set of distributions.

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, Wang-BBP 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 high-amplitude 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).

Figure 10 As in Fig. 2, but with each experiment in TRAIL coloured by ϵ. Stable experiments are marked with black crosses. Solid lines are from Field et al. (1997) and dashed lines are from Zhong et al. (2013).

Figure 11 Phase diagram showing the stable (black crosses) and periodic (coloured by ϵ) behaviour of falling particles as a function of area ratio and Reynolds number. Equation (11) is shown as the black line.

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 motion-type parameter ϵ was subsequently calculated (Fig. 12) (Eq. 8).

Figure 12 Mean inclination angle, θ_mean (a), θ_tilt (b), and θ_amp (c) by Reynolds number, coloured by ϵ. Histogram of ϵ (d).

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⁵). 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 cross-sectional 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

$$\begin{array}{}\text{(11)}& {\log }_{\mathrm{10}}\left(\mathit{Re}\right)>{\displaystyle \frac{\mathrm{2.82}-\mathrm{Arearatio}}{\mathrm{0.87}}}.\end{array}$$

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 near-zero 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 plate-like 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 least-squares 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 dendritic-shaped 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.

Figure 13 Variation of components of motion for all periodic experiments. Amplitude of horizontal velocity relative to mean vertical velocity compared to amplitude of vertical velocity relative to mean vertical velocity (a), θ_amp (b), and θ_tilt (c), coloured by ϵ. Dashed lines are fit to all experiments (a), experiments where ϵ < 0.5 (b), and experiments where ϵ > 0.5 (c).

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 near-zero (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.

Figure 14 Scatter plot of Strouhal numbers (St_θ) for each unstable experiment versus Reynolds number, coloured by ϵ. Box plots of observed Strouhal numbers for particles with θ_amp > 2.5° alongside.

Figure 15 Scatter plot of azimuth Strouhal numbers (St_ϕ) for each unstable experiment versus Reynolds number, coloured by ϵ. Experiments where θ_amp < 2.5° are marked with crosses. Box plots of observed Strouhal numbers for particles with θ_amp > 2.5° alongside.

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.

Figure 16 Scatter plot of azimuth Strouhal numbers (St_ϕ) versus inclination Strouhal numbers (St_θ). The slope of the dashed line is 2.

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 non-zero, 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 (non-perfect 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{‾}}{\mathit{\theta }}$, is close to zero, as particles are stable, and at higher Reynolds numbers, particles become unstable and $\stackrel{\mathrm{‾}}{\mathit{\theta }}$ increases to a steady, non-zero value. To quantify and compare the onset of unstable motions for different particle properties, a logistic curve is fit to the data using a least-squares 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:

$$\begin{array}{}\text{(12)}& \stackrel{\mathrm{‾}}{\mathit{\theta }}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{unstable}}}{\mathrm{1}+e^{-k\left({\log }_{\mathrm{10}}\right(Re)-{\log }_{\mathrm{10}}(\mathit{Re}{)}_{\mathrm{onset}})}}}+{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{stable}},\end{array}$$

where Re_onset is the value of the function's midpoint, ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{unstable}}+{\stackrel{\mathrm{‾}}{\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{‾}}{\mathit{\theta }}$ value of the function. Particles with higher area ratio typically have bigger oscillations with larger $\stackrel{\mathrm{‾}}{\mathit{\theta }}$. Dendrites with low area ratios and stellar crystals are stable at very high Re and typically have a smaller $\stackrel{\mathrm{‾}}{\mathit{\theta }}$ when unstable.

Figure 17 Mean inclination angle, $\stackrel{\mathrm{‾}}{\mathit{\theta }}$, against Reynolds number, Re, for each experiment (scatter) grouped by particle shape (a) and area ratio (b) with overlaid fitted logistic functions.

4 Discussion and Conclusion

Ten different plate-like 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, non-zero 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 non-zero mode. The results do not support the common assumption of Gaussian orientation distributions with a zero-modal angle during unstable motions: the distributions of θ have non-Gaussian distributions and non-zero modes.

In the spiralling regime, components of velocity are held relatively constant compared to the other cases, effectively making the spiralling cases a quasi-steady mode with a non-zero near-constant 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 (non-dimensionalised 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³. The onset for circular discs was found to be as low as Re=197.

A motion-type 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 (Wang-BBP) 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 non-Gaussian, 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 broad-branched crystal forms.

Ground-based snowflake imagers have reported preferentially non-horizontal 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 vertical-pointing 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 near-horizontally 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. $>\mathrm{0.1}{\text{m}}^{\mathrm{2}}{\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 turbulence-induced 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.

Appendix A: Appendix

Table A1 Output variables from logistic fits as presented in Fig. 17 for each shape.

Table A2 Output variables from logistic fits as presented in Fig. 17 by binned area ratio.

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:

$$m=\mathrm{0.0038}{d}^{\mathrm{2}},$$

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

$$I_a={\displaystyle \frac{\mathrm{1}}{\mathrm{16}}}M{D}^{\mathrm{2}},$$

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

$$I^{*}={\displaystyle \frac{I_a}{{\mathit{\rho }}_{\mathrm{f}}{D}^{\mathrm{5}}}},$$

where we set the density of air to be 1.2 kg m⁻³. 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.