TRAIL part 2: A comprehensive assessment of ice particle fall speed parametrisations

Results are presented from a comprehensive experimental campaign studying the aerodynamics of 3D‐printed analogues of ice particles. Measurements of the drag coefficient of the analogues were acquired by using a novel experimental approach to digitally reconstruct the analogues' trajectory and orientation as they fall through a quiescent viscous liquid, using images acquired by a series of digital cameras. The data are used to evaluate commonly used parametrisations of ice particle fall speeds. We find that the accuracy of each of the parametrisations reduces at Reynolds numbers greater than approximately 100, as a result of effects associated with the onset of unsteady oscillating, fluttering or tumbling motions as the analogues fall. We propose a new fall speed parametrisation that predicts the measured drag coefficients at high Reynolds number with an accuracy consistent with that of the leading parametrisation at low Reynolds number.


INTRODUCTION
In part 1 of this study (McCorquodale and Westbrook, 2020) we reported an extension to the analogue method of investigating the aerodynamics of ice particles through the use of a novel algorithm (the Trajectory Reconstruction Algorithm implemented through Image anaLysis, TRAIL) that reconstructs the time-resolved trajectory and orientation of falling analogues. Using this system, the trajectories of 3D-printed analogues of ice particles are reconstructed across a sample volume of 200 × 200 × 200 mm to an accuracy of better than ±0.5 mm, whilst the orientation of the falling analogues within this volume are reconstructed to an accuracy of ±2.5 • . In part 2 of this study, we report new results from a comprehensive measurement campaign investigating the aerodynamics of 3D-printed analogues utilising the TRAIL methodology. Crucially, measurements are reported not only for analogues with simple geometries but also for analogues with complex irregular shapes. Here we focus on measurements of drag coefficients, which are used to test the accuracy of parametrisations of ice particle fall speeds.
Representations of ice particle fall speeds are used extensively in atmospheric physics, in applications includ-ing estimating microphysical properties from Doppler radar measurements (Matrosov, 2011) and from surface measurements of snowflake size and fall speed (von Lerber et al., 2017), predicting cirrus cloud properties in global climate models (Mitchell et al., 2008;2011;Jensen et al., 2018) and within numerical simulations of ice particle growth (Westbrook et al., 2004) that facilitate the calculation of scattering properties of hydrometers (Tyynelä et al., 2011). However, the representation of the particle fall speeds can represent a notable source of uncertainty within these applications. For example, numerical experiments by Jakob (2002) show that modest changes in representations of cirrus cloud ice crystal fall speeds can produce changes in the Earth's atmospheric radiation budget which are comparable to the changes associated with doubling the CO 2 .
Unfortunately, the accuracy with which parametrisations are able to predict the fall speed of natural ice particles remains unclear. This is due to the difficulties encountered in acquiring accurate measurements of the fall speed of natural ice particles; existing instrumentation (see e.g. Kruger and Krajewski, 2002;Schönhuber et al., 2007;Newman et al., 2009;Garrett et al., 2012) does not facilitate measurements of some key parameters. In particular, challenges are encountered in precisely determining the mass and geometric shape of falling ice particles (see part 1 for further details). In addition, whilst highly accurate and detailed data can be acquired in studies employing computation fluid dynamics (see e.g. Cheng et al., 2015;Nettesheim and Wang, 2018), only a limited range of simple, idealised geometries have been studied using this approach as a result of the significant computational resources required. Consequently, the absence of accurate and comprehensive data containing measurements of the mass, geometric shape and fall speed of natural ice particles has prevented rigorous evaluations of fall speed parametrisations. Westbrook and Sephton (2017) recently reported an assessment of ice particle fall speed parametrisations using measurements of 3D-printed analogues of ice particles falling in viscous liquids. The use of 3D-printed analogues yields distinct advantages when compared with natural ice particles; accurate measurements of fall speed can be collected under carefully controlled laboratory conditions using particles of known mass and geometric shape. Using this approach, Westbrook and Sephton (2017) reported that current parametrisations can systematically overestimate fall speeds by as much as 80%. However, the results of that study were subject to a number of limitations. In particular, only a small range of analogues were used in the experiments, the geometry of which was highly idealised. The range of analogues with irregular geometries was severely limited. In addition, estimates of the projected area, required to compute the drag coefficient used within fall speed parametrisations, were subject to large uncertainties because the projected area of each particle was calculated by estimating the observed orientation of the analogues in free fall by eye. Whilst this approach may yield reliable estimates of projected area for analogues with simple geometries (such as plates), this method is particularly unreliable for analogues with complex irregular geometries. Consequently, the accuracy with which parametrisations are able to predict the fall speed of ice particles with complex irregular geometries, such as polycrystals or aggregates, remains unknown, which is unfortunate since observations indicate that these particle types are the most prevalent in the atmosphere (e.g. Korolev et al., 2000;Field and Heymsfield, 2003).
In this study, we utilise the TRAIL methodology to evaluate ice particle fall speed parametrisations using analogues with a wide range of geometries, including complex irregular shapes representing polycrystals and aggregates. A brief summary of existing parametrisations of ice particle fall speeds is presented in Section 2. A summary of the measurement campaign conducted, including a description of the 3D-printed analogues used, is reported in Section 3. Measurements of drag coefficients, which comprise a new dataset that we believe to be more accurate and comprehensive than used in previous studies, are reported in Section 4. The measurements of drag coefficients are used to assess fall speed parametrisations in Section 5. A revised fall speed parametrisation with greater accuracy, particularly at high Reynolds number, is presented in Section 6. Conclusions are presented in Section 7.

OVERVIEW OF TERMINAL VELOCITY PARAMETRISATIONS
Parametrisations of the terminal velocity of ice particles are expressed through characterising the relationship between the drag coefficient of a particle C D and its Reynolds number Re, as a function of the particle shape. Here, the Reynolds number is given by Re ≡ u t D∕ , where u t denotes terminal velocity of an ice particle with maximum span D, falling in a fluid of kinematic viscosity . Likewise, the drag coefficient is given by where V is the volume of the particle, A is the projected area of the particle (normal to the fall direction), p and f are the densities of the particle and fluid, respectively, and g is the acceleration due to gravity. Commonly, parametrisations are re-cast as relationships between the Reynolds number Re and Best number Be ≡ C D Re 2 , since the Best number can be estimated directly from knowledge of a particle's shape and size and of the fluid in which it is falling.

Drag on a sphere
As a canonical case, the drag force acting on a spherical particle has been extensively studied, and comprehensive reviews are available, such as Clift et al. (1978). Here we focus on the study of Abraham (1970), which underpins current parametrisations of the drag coefficient of ice particles. Abraham (1970) proposed a theoretical model for the drag force acting on a spherical particle by equating the drag force experienced in a viscous fluid to that experienced by an "equivalent assemblage" of a particle and an attached boundary layer falling in an inviscid fluid. In this approach, the boundary layer comprises the region of flow in the viscous fluid in which the motion of the particle induces shear; the boundary layer thickness effectively denotes a point at which the effects of friction in the fluid are negligible. The effect of this boundary-layer model is to increase the projected area of the assemblage, relative to the original spherical particle, by a factor equal to (1 + 2 ∕D) 2 , where D denotes the diameter of the spherical particle (Abraham, 1970). By equating the forces acting in these equivalent systems, and assuming that the average depth of the boundary layer can be expressed as = 0.5D 0 Re −1∕2 , where 0 is a dimensionless coefficient, one obtains a relationship between the drag coefficient and Reynolds number of the form where C sph D denotes the drag coefficient of the spherical particles and C 0 is a dimensionless coefficient (see Abraham, 1970 for details). The corresponding Re-Be relationship (recall that Be ≡ C D Re 2 ) is The coefficients 0 and C 0 are evaluated using experimental observations. A comparison against experimental data indicates that Equations (2) and (3) are reliable for Reynolds numbers up to order 10 3 (Abraham, 1970). However, this model gradually becomes less reliable for Re ≳ 10 3 ; it predicts that drag coefficients continue to reduce as the Reynolds number increases, as shown in Figure 1, but experimental data indicate that the drag  Abraham (1970) and Clift and Gauvin (1970) (see legend) coefficient of spherical particles rises slightly then remains almost constant between Re ≈ 750 and Re ≈ 3 × 10 5 (see e.g. Clift et al., 1978). This discrepancy occurs as a result of significant changes in the structure of the flow around spherical particles as the Reynolds number increases. The structure of the flow around spherical particles has been extensively studied. For brevity, here we briefly discuss only key changes that occur as the Reynolds number increases; Clift et al. (1978) report a more detailed description of these changes. For Re ≪ 1, the flow past a spherical particle is approximately symmetric fore and aft. As the Reynolds number increases, the asymmetry of the flow becomes more pronounced. Flow separation occurs at Re ≈ 20, which gives rise to the formation of a vortex ring close to the rear stagnation point. As the Reynolds number increases, the separation point moves closer to the front of the sphere, so the attached wake becomes wider. At Re ≈ 130 the flow becomes unsteady and the vortex ring that makes up the wake begins to oscillate. Vortices are shed from the wake for Re ≳ 400. At this point, the wake consists of fairly steady flow close to the rear of the sphere and turbulent flow further behind the particle as a result of the vortex shedding, which can also cause fluctuations in the motion of the particle. As the Reynolds number increases, vortices are shed closer to the sphere, so a larger turbulent far-wake is formed. At Re ≈ 3 × 10 5 , a turbulent boundary layer forms towards the rear of the particle and the enhanced momentum transfer that occurs in the the turbulent boundary layer causes the separation point to move towards the rear of the particle. This abrupt change in the structure of the flow gives rise to a sharp reduction in drag coefficient -the so-called drag crisis (see e.g. Heymsfield and Wright, 2014).
From the above description, it is apparent that the reduction in accuracy of the Abraham (1970) parametrisation (i.e. Equation 2) coincides with the onset of unsteadiness in the wake of spherical particles. More complex curves are required to accurately model drag coefficients in this regime, such as that proposed by Clift and Gauvin (1970), which is accurate to within 6% up to Re ≈ 3 × 10 5 (Clift et al., 1978, p. 112). A comparison of the drag curves proposed by Abraham (1970) and Clift and Gauvin (1970) is shown in Figure 1.

Drag on an ice particle
Parametrisations of the terminal velocity of ice particles have used the model of Abraham (1970) to represent the underlying relationship between the drag coefficient and Reynolds number of falling particles. Consequently, the defining aspect of ice particle parametrisations is the approach taken to generalise this result for spherical particles to account for the more complex geometries associated with natural ice particles. Mitchell (1996) (hereafter M96) assumed that Equations (2) and (3) can be applied to falling ice particles by altering the empirical coefficients C 0 and 0 to universal values (M96 use C 0 = 0.6, 0 = 5.83). That is, M96 argued that a particle's fall speed explicitly depends on only three properties of the particle: its mass, projected area and maximum dimension. Other variables that directly characterise ice particle geometry, such as aspect ratio and porosity, were not incorporated into the parametrisation. In contrast, adjustments that amend the form of Equations (2) and (3) to account for different ice crystal geometries were proposed by Böhm (1989) (hereafter B89) and Heymsfield and Westbrook (2010) (hereafter HW10). In both cases, the authors introduce a parameter that quantifies the shape of the falling particle into Equation 2, which takes the modified form

Steady flow
where A r denotes the area ratio. (B89 report the use of empirical constants n = −3/4, C 0 = 0.6 and 0 = 5.83, whereas HW10 adopt n = −1/2, C 0 = 0.35 and 0 = 8.0.) The area ratio is defined as the ratio of the projected area of the particle (normal to the fall direction) to the area of a circle that circumscribes the particle projection. The corresponding adjustment to Equation 3 is given by substituting A −n r Be for Be. Consequently, in these parametrisations, a particle's fall speed also directly depends upon a quantity that characterises the ice particle's geometry (the area ratio A r ). An additional adjustment was also proposed by B89 according to a classification of particles into columnar or planar categories. For columnar particles, the equivalent circular cross-section diameter d * ≡ 2(A∕ ) 1∕2 replaces the maximum span D as the characteristic length scale in Re and A r (i.e. A r = 1). However, it is not clear how this additional adjustment can be implemented in practice; this adjustment requires a particle to be classified as either planar or columnar, a requirement which is not possible for many natural ice particles (e.g. irregular polycrystals and aggregates). Consequently, in this work, we neglect this additional adjustment, and all subsequent references to the B89 parametrisation refer specifically to the modified form of the C D -Re curves given by Equation 4. 1 In this study, we do not consider the parametrisation proposed by Böhm (1992), since this parametrisation relies on the idealisation of ice particles into porous prolate or oblate spheroids. In addition, Böhm (1992) introduce an aspect ratio as a control parameter that is not clearly defined and which is difficult to implement in practice for irregular ice particles. In any case, we note that the salient points regarding key limitations of existing ice particle fall speed parametrisations, which we identify and discuss in Sections 4 and 5, can also be applied to this parametrisation.

Unsteady flow
Amendments to these parametrisations have been proposed at high Reynolds number (i.e. Re ≳ 10 3 ). Recall that Equation 2 underpins these parametrisations, but this equation under-predicts the drag coefficient when the flow around falling particles is unsteady. This effect gives rise to overestimates of the true fall speed (Mitchell and Heymsfield, 2005). Amendments to parametrisations that increase the drag coefficient predicted under these conditions are sometimes referred to as "turbulent corrections", for the reasons explained below. Böhm (1992) report that fall speed parametrisations become less accurate due a transition to turbulent flow in the vicinity of irregular particles, which he believed to occur at 1,000 ≲ Re ≲ 2,000. Böhm (1992) proposed that the drag coefficient in these "turbulent" conditions, denoted by C Dt , can be computed using the formula 1 We believe there is a typographical error in Equation (9) of B89, which in Böhm's notation should read C D = 0.6(1 + 5.83Re −1∕2 ) 2 .
where C Dl denotes the drag coefficient under "laminar" conditions (i.e. predicted using parametrisations described above), Be 0 = 2.8 × 10 6 and = 1.6. In practice, this formula gives rise to drag coefficients that are unchanged when Re is small but which are enhanced by the factor when Re is large (i.e., Be ≫ Be 0 ). The choice of Be 0 is related to the physical argument above; a value of Be 0 = 2.8 × 10 6 is equivalent to Re ≈ 2,000. The origin of Böhm's assertion that the transition to turbulence occurs at 1,000 ≲ Re ≲ 2,000 can be traced to the study of Becker (1959), who reported that the wake of a freely falling particle becomes completely turbulent at a Reynolds number near 2,000. However, this statement must be interpreted with caution for several reasons; this conclusion was drawn from observations of the wake behind spheres and discs, whilst the definition of the characteristic length scale within the Reynolds number used by Becker (1959) is notably different from that used here and by Böhm (1992). (The characteristic length scale used by Becker (1959) is the diameter of a sphere of equal surface area to the particle.) In addition, work on spherical particles now confirms that the onset of turbulence in the wake is a more nuanced process than assumed by Böhm (1992) (see Section 2.1). Consequently, the Reynolds number at which the wake of an irregular ice particle becomes completely turbulent may be substantially different from the proposed range 1,000 ≲ Re ≲ 2,000. Moreover, the origin of the evidence for Böhm's choice of is unclear. It appears that this choice may have been informed using data on spherical particles; Equation 2 under-predicts the drag coefficient of spherical particles by approximately 60% at Re ∼ O(10 5 ) (consistent with the factor = 1.6). M96 adopts the same approach, but the factor = 1.6 used by Böhm (1992) is changed to = 1.3 to better conform with drag coefficient data for hail reported by Knight and Heymsfield (1983). M96 also argue that the physical mechanism that gives rise to the increase in drag coefficient is the onset of turbulence in the flow near Re ≈ 1, 000, in accordance with Böhm (1992), although this is not supported by experimental evidence.
In contrast, Mitchell and Heymsfield (2005) propose a turbulent correction that modifies Equation 3 by the subtraction of a power-law term −a 0 Be b 0 (where a 0 = 1.7 * 10 −3 and b 0 = 0.8) leading to smaller Re (and hence u t ) for a given Be and implicitly higher drag for a given Re. The modified drag coefficient curves are almost indistinguishable from their original forms when Re ≲ 100 but predict significantly higher drag coefficients at large Re. Unlike Böhm (1992) and M96, Mitchell and Heymsfield (2005) argue that this correction applies only to aggregate particles and arises through differences that occur in the boundary-layer structure. In particular, they argue that the increase in projected area associated with the assemblage of an aggregate particle and its boundary layer cannot be reliably approximated by the factor (1 + 2 D) 2 = (1 + 0 Re −1∕2 ) 2 . They attribute this to two effects: flow through the voids within an aggregate, and a breakdown of the boundary-layer assumption = 0.5 0 DRe −1∕2 . Mitchell and Heymsfield (2005) argue that this assumption is unreliable as it utilises the maximum span of the aggregate as its characteristic dimension, rather than a dimension that lies between the maximum span of the aggregate and the span of its component crystals. However, it is unclear why these mechanisms would be unique to aggregate particles when similar arguments could equally be applied to any ice particle with a geometry that exhibits branches or voids.

ICE PARTICLE ANALOGUES
Experiments were conducted using the TRAIL methodology detailed in part 1, to which we refer the reader for a detailed description of the experimental procedure. In brief, the fall of 3D-printed analogues through a quiescent viscous liquid (water-glycerol mixtures) was recorded using a series of synchronised digital cameras. A custom algorithm (TRAIL) was used to digitally reconstruct the trajectory and orientation of the 3D-printed analogues from the images recorded, which enables the drag coefficient and Reynolds number to be estimated. This approach exploits dynamic similarity; when falling steadily, the aerodynamics of ice particles is a function of only the geometric shape of the particles and the Reynolds number at which they fall (List and Schemenauer, 1971). Consequently, when falling steadily, relationships between dimensionless quantities, such as C D , Re and Be, are the same for ice particle analogues falling in viscous liquids and natural ice particles falling in the atmosphere. Here, we report results from 844 experiments, comprising measurements of particles with 118 different geometries, each falling at a range of Reynolds numbers; for a given analogue geometry, the Reynolds number could be varied by changing the viscosity of the ambient fluid or changing the size of the 3D-printed analogue. A diverse range of ice particle analogues were used, as shown in Figure 2, typically ranging in size from 1 to 3 cm. (Full details of the geometry of the analogues used can be found in the supplementary information.) The analogues were designed such that their geometry has a microphysical basis, albeit idealised (see e.g. Kikuchi et al., 2013) and include representations of plates, plate polycrystals, columns, capped columns, bullet rosettes, mixed rosettes and aggregates. The geometries used were designed manually using a number of different approaches. Analogues representing ice particle F I G U R E 2 Ice particle geometries used in experiments habits with simple geometries were designed using an open-source computer-aided design package (OpenSCAD, www.openscad.org). A range of planar analogues were also designed using the numerical schemes proposed by Reiter (2005) and Demange et al. (2017) or taken from the numerical study of Nettesheim and Wang (2018). A range of analogues representing unrimed and rimed aggregate particles were obtained from numerical aggregation schemes (Westbrook et al., 2004;Leinonen and Szyrmer, 2015) or from micro-computed tomography (micro-CT) scans of natural ice particles (Giorgia Tagliavini, personal communication).
Some aspects of the geometry of the analogues were idealised to simplify the fabrication process. For example, for analogues representing unrimed aggregates, it was necessary to overlap the individual crystals of the aggregate very slightly at the point of contact in order to increase the structural integrity. Despite the use of steps to idealise the geometry of the analogues, we feel that these particles are a logical step to improving our understanding of the aerodynamics of both simple and complex natural ice particles. In particular, we stress that the analogues exhibit a much broader range of geometries than previous analogue experiments (for example List and Schemenauer, 1971;Kajikawa, 1971;Jayaweera, 1972;Westbrook and Sephton, 2017). Of particular note is the large range of analogues with irregular geometries representing polycrystals and aggregates; measurements from equivalent analogues have not previously been reported.
Following fabrication, care was taken to evaluate the accuracy with which the analogues were 3D-printed. In particular, the mass of each analogue was determined using precision scales and compared against a theoretical estimate calculated from the volume of the digital model of the analogue and the density of the printed analogues ( p = 1,174 ± 2 kg⋅m −3 ). Analogues were discarded if significant differences were observed. In these instances, the original digital models were further idealised to simplify the fabrication of the analogue and improve the accuracy of the printing process. However, modest differences (<5% on average) between the measured and theoretical mass remain as a result of the limitations inherent in the 3D-printing process. To mitigate any impact of these differences on the uncertainty in measurements of C D , Re and Be, the measured masses of the 3D-printed analogues are used in all subsequent calculations of C D and Re. However, since the small differences in mass correspond to differences in shape between each printed analogue and its digital counterpart, these printing errors may give rise to a modest uncertainty in measurements of A, D and, through error propagation, C D , Re and Be. That is, A and D were determined from digital reconstructions of the trajectory of the falling analogues, in which the geometric shape of each analogue is assumed to exactly match its digital model (see part 1 for details). Consequently, to better understand the nature of the printing errors, the dimensions of the 3D-printed analogues were measured and compared with dimensions obtained from the corresponding digital models. We note that, in general, printing errors principally arose through the inability of the fabrication technique to faithfully reproduce small-scale features of the analogues, typically manifesting as errors in the thickness of thin, plate-like elements or slender, column-like elements. The magnitude of the uncertainty associated with the printing errors on the measurements of A, D (and, through error propagation, C D , Re and Be) was estimated by comparing the measurements of A and D against a set of independent estimates computed directly from images acquired from a digital camera that recorded a projection of the analogues approximately normal to the fall direction (see part 1). A comparison of the two estimates indicates that the uncertainty in D associated with printing errors is negligible. In contrast, a comparison of the two estimates indicates that a conservative estimate of the uncertainty in measurements of A arising from print errors would be 5%. In part 1 we estimated that the TRAIL methodology results in measurements of C D and Re that are accurate to within approximately 5%. Through error propagation, by including uncertainty associated with printing errors when fabricating analogues, we estimate that measurements of Re and C D reported here are, on average, accurate to within approximately 5% and 7%, respectively.
Finally, we note that the data reported in subsequent sections are broken down using a division of particles into "plate-like", "column-like", "rosettes' and "aggregates" categories; the assignment of analogues into these categories is also shown in Figure 2. We stress that the sorting of analogues into these categories is not strictly based on mircophysical definitions; some plate polycrystals have been sorted into the plate-like category since they exhibit broadly similar aerodynamic properties.

MEASUREMENTS OF C D
Measurements of C D and Re are shown in Figure 3, in which the area ratio of the analogues is also shown. (A summary of the drag coefficient data is also provided in Tables S1 to S4 of the supplementary material.) Note that the measurements shown in Figure 3 exhibit significant variability in C D (Re) from shape to shape. In particular, at low Reynolds number (for Re < O(10 2 )), Figure 3 reveals a geometric effect whereby analogues with low area ratio exhibit higher drag coefficients than analogues with higher area ratio. This trend is apparent within each of the particle categories shown in Figure 2, but is most apparent for particles in the plate-like category on account of the greater range in area ratio of analogues in this category. Figure 3 indicates that this effect is the dominant source of scatter within the data for Re < O(10 2 ); once the influence of A r is considered, Figure 3 shows that the drag coefficient is only weakly linked to the particle category. This indicates that, at low Reynolds number (Re < O(10 2 )), the drag coefficients of ice particles may be accurately approximated using functions that characterise the shape of the analogue through a description of the projection of the particle in the fall-normal direction, such as the area ratio, irrespective of the full 3D shape of the particle. However, Figure 3 indicates that the observed relationship between C D and A r does not extend to Re ≳ O(10 2 ), as more complex geometric effects occur. These effects are most apparent in the results for analogues in the plate-like category shown in Figure 3a, which indicate that there is a sharp increase in the drag coefficient at Re ≈ O(10 2 ) for analogues with high area ratio (A r ≳ 0.6). In contrast, measurements indicate that the drag coefficient continues to reduce monotonically for Re ≳ O(10 2 ) for analogues with low area ratio (A r ≲ 0.3). Closer inspection of the data reveals that there is in fact a gradual transition between these extremes, whereby analogues with intermediate area ratio exhibit a less marked increase in drag coefficient at Re ≈ O(10 2 ). Consequently, as the Reynolds number increases, the overall dependence of the drag coefficient on the area ratio diminishes, as shown in Figure 3a. However, note that, for Re ≈ O(10 3 ), analogues with low area ratio do exhibit slightly lower drag coefficients than those with higher area ratio; we stress that this trend is the opposite to that observed for Re ≲ 10 2 .
Evidence for these trends can also be seen for analogues within the other categories (shown in Figure 3b-d), although increases in C D are in general smaller in these cases. The smaller increases in C D observed (relative to the low-Reynolds-number trend) can be attributed, in part, to the absence of analogues with very high area ratio (A r ≳ 0.7) in these categories. However, even when accounting for this, Figure 3 does indicate that the relative increases in C D within the column-like, rosette and aggregate categories are less significant than in the plate-like category, such that drag coefficients are lower on average at Re ∼ 10 3 . On the other hand, Figure 3c shows a series of measurements of C D for analogues in the rosette category that exhibit higher drag coefficients for Re ≳ 200, more consistent with the trend exhibited by analogues in the plate-like category (shown in Figure 3a). These measurements correspond to analogues that represent complex polycrystals consisting of a combination of column and categories (see Figure 2). Different symbols are used to show whether analogues fell in a steady ("+") or unsteady ("o") regime. In each case, the colour of the symbols denotes the area ratio of the analogues, as shown by the colour bar, which applies to all plots plane-type crystals (classification CP5 in Kikuchi et al., 2013). Similar, anomalously high drag coefficients are also shown in Figure 3d for Re ≳ 800, which correspond to measurements from some of the analogues representing rimed aggregates. We note that, in each of these cases, aspects of the geometry of the analogues are similar; they each exhibit large approximately flat surfaces which lie in approximately the same plane, akin to the geometry of analogues in the plate-like category.
Insight into the origin of the trends observed for Re ≳ 10 2 can be deduced from digital reconstructions of the trajectory of the analogues, which were determined from the experimental measurements (see part 1 for details of the reconstruction process). From the reconstructed trajectories, it was possible to determine whether particles fell steadily or unsteadily, wherein particles exhibited oscillating or fluttering motions. We note that many of the irregular particles (rosettes and aggregates) exhibited spiralling trajectories whereby they rotated around a vertical axis as they fell. We stress that spiralling trajectories are not necessarily an indication of the onset of unsteady motion or unsteady flow; rather, they are a result of a constant rotation of the particle due to the centre of mass and centre of drag not being co-located. Analogues that exhibited spiralling trajectories were judged to have fallen steadily when they exhibited a clear (stable) preferred orientation as they fell. (Full details on the observed kinematics, including the preferred orientation of the irregular analogues, are to be reported separately in a future paper.) The relative motion of the analogues (steady versus unsteady) is also shown in Figure 3, which clearly illustrates that the change in the shape of the C D -Re data at Re ∼ O(10 2 ) is coincident with the onset of unsteady motions as the particles fall. By means of explanation, we note that the onset of unsteady motions is a symptom of a change in structure of the wake of falling analogues. That is, the (unsteady) aerodynamic forces that give rise to the observed unsteady motions arise as a result of disturbances within the flow that occur through wake instability (Fernandes et al., 2007;Ern et al., 2012). Unlike Böhm (1992) and M96, we specifically refrain from referring to this process as the "onset of turbulence" in the wake. As indicated in Section 2.1 for spherical particles, the structure of the wake is thought to gradually change as the Reynolds number increases, whilst the onset of turbulent boundary layers actually gives rise to a sharp reduction in drag coefficient (Clift et al., 1978). Indeed, in the current experiments, we observed that, as the Reynolds number increases, some analogues exhibit several different phases in which a different type of motion was observed. Each type of motion is thought to arise as a result of a distinct change in wake structure. We note also that the anomalous measurements of drag coefficient in the rosette and aggregate categories can be linked to the type of unsteady motions observed. Recall that, in each of the anomalous cases, the analogues exhibit large approximately flat surfaces which lie in approximately the same plane; when in free fall at high Re this plane lies approximately normal to the fall direction and the analogues exhibit large-amplitude fluttering motions that are similar to those exhibited by analogues in the plate-like category (for comparable examples see Esteban et al., 2019). The onset of large-amplitude fluttering motions has previously been reported to be associated with an increase in drag coefficient (Ern et al., 2012). However, detailed measurements of the flow field around the falling particles are required to fully understand how the structure of the wake influences the drag coefficient and how this in turn is influenced by the large range of geometric shapes considered. Unfortunately, simulations and experimental measurements have been reported for simple, regular shapes (see for example Nettesheim and Wang, 2018;Esteban et al., 2019) but, as far as we are aware, detailed measurements of the flow field are not available for particles with irregular geometries similar to those used here.
Regardless, in the current context, a crucial aspect of this phenomena is the critical Reynolds number Re c at which we observe the onset of unsteady motions that give rise to a change in shape in the C D (Re) curve. Of particular interest here is the relationship between Re c and A r ; we have seen that much of the variability in C D (Re) from shape to shape is controlled by A r , so one may expect A r to be an important parameter in determining Re c . Analysis of digital reconstructions of the trajectory of the falling analogues enables us to estimate lower and upper bounds for the critical Reynolds number Re c , based on the fall regime observed in each experiment. We stress that the limited resolution (in the Re parameter space) of measurements for each analogue can result in a large range between the lower and upper bounds of the critical Reynolds number Re c . This uncertainty in the critical Reynolds number Re c , combined with the use of multiple analogues with very similar (or even the same) A r , means that the data cannot be clearly conveyed on plots showing Re c as a function of A r . Instead, here we report the measurements of Re c against an identification number, as shown in Figure 4, in which the area ratio of the analogues is also shown. Further information to enable the identification of each analogue geometry within Figure 4 is available in Tables S1 to S4 of the supplementary information. Figure 4 indicates that there is significant variation in Re c , which ranges between Re c ≈ 50 and Re c > 10 3 , dependent on the geometry. In general, the data indicate that analogues with the highest Re c values have low area ratios. However, Figure 4 shows that this relationship is not unique; we observe that the onset of unsteady motions is not only a function of Re and A r , but also depends on other details of the particle shape. For example, within the plate-like category, there are a number of analogues that exhibit A r ≈ 0.83, and across these particles the critical Reynolds number varies by a factor of approximately 5. This observation is consistent with results from a range of previous model studies. For example, previous studies investigating the aerodynamics of thin circular discs have reported that the fall regime is a function of Re, density ratio p ∕ f and aspect ratio t/D (Jayaweera, 1965;Field et al., 1997;Auguste et al., 2013). Vincent et al. (2016) further reported that the presence of a central hole stabilises the falling motion of thin discs. Moreover, in a model study of the aerodynamics of dandelion seeds, Cummins et al. (2018) reported that the critical Reynolds number Re c of porous discs sharply increases for A r ≲ 0.25. By means of explanation, it has been argued that the onset of the unsteady motions is related to the intensity of vorticity generated in the wake of a falling particle (Auguste et al., 2013). For example, a central hole stabilises the falling motion of thin discs as the hole results in a distinctly different wake structure (comprising a pair of counter-rotating vortex rings) from a solid disc, such that the flow-induced torque on the annular disc is reduced (Vincent et al., 2016). As such, deducing overall trends in Re c for particles of irregular shape is difficult, since the structure of the wake is a complex function of the 3D shape of a particle; distinctly different wake structures can even be observed for particles with very similar geometries, such as different types of thin porous discs (see for example Vincent et al., 2016;Cummins et al., 2018). Consequently, relationships between the critical Reynolds number at which we see the onset of unsteady motions and the geometric shape of ice particles cannot be reliably approximated in terms of simple shape factors, such as A r . This complex relationship contributes to the scatter in the drag coefficients shown in Figure 3 at high Reynolds number.
As noted in part 1, care should be taken when interpreting data in which analogues were observed to fall unsteadily, since the conditions for dynamic similarity, which we are exploiting here, are only strictly satisfied when the analogues fall steadily. When falling unsteadily, the aerodynamics of falling particles could be sensitive to the relative densities of the particle p and the ambient fluid f in which it is falling. Studies of spherical and cubic  Figure 2). Upward (downward) triangles are used to denote experiments for which only a lower (upper) bound of Re c could be determined; in these cases the actual value of Re c lies above (below) the symbol shown. Also shown is the area ratio of the ice particle analogues. Further information on the identification of each analogue is available in Tables S1 to S4 of the supplementary material particles have reported that the relationship between Re and C D is sensitive to the density ratio ≡ p ∕ f for Re ≳ O(10 2 ) (Jenny et al., 2004;Veldhuis et al., 2009;Horowitz and Williamson, 2010;Ern et al., 2012;Seyed-Ahmadi and Wachs, 2019). In these studies, for a spherical particle, the drag coefficient was found to be independent of for > c . In each study it was found that c < 0.5; by means of comparison, in the current work 1.03 ≤ ≤ 1.18. For cubic particles, for each value of , the drag coefficient was found to increase significantly at Re ≳ O(10 2 ), consistent with the results shown in Figure 3a, before attaining an approximately constant value at higher Reynolds number. However, more pronounced increases in C D were reported for cubes with low (Seyed-Ahmadi and Wachs, 2019); at Re ≈ 200, the drag coefficient for = 0.2 was found to be approximately 20% larger than when = 7. A similar trend has previously been observed for thin circular discs (Jayaweera, 1965). We stress that the increases in C D at Re ≳ O(10 2 ), relative to the trends identified at low Reynolds number, are large compared with the sensitivity of C D to previously reported in the literature. We therefore assume that C D (Re) measured in our experiments at high Re in unsteady conditions is indeed the same as C D (Re) for ice particles falling in air ( ≈ 1,000). Ongoing work, featuring experiments using micro-scale 3D-printed analogues falling in air and a comparison of our measurements with data acquired using computational fluid dynamics (Tagliavini et al.), will be used to assess the validity of this assumption.

Mitchell (1996) parametrisation
Mitchell proposed a parametrisation of ice particle fall speed of the form C D = C 0 (1 + 0 Re −1∕2 ) 2 . A comparison Percentage error when using the M96 parametrisation to predict the drag coefficient of ice particle analogues, where C M96 D denotes the predicted drag coefficient; data have been separated according to the Reynolds number at which the analogues fall (see panel labels for details). In each panel, analogues of different habits are shown separately in accordance with the categories shown in Figure 2 (see legend in (a) which applies to all panels). Different symbols are also used to show whether analogues fell in a steady ("+") or unsteady ("o") regime of this parametrisation with measurements of C D and Re is shown in Figure 5a, whilst the percentage error in the drag coefficient predicted by the M96 parametrisation is shown in Figure 5b, where C M96 D denotes the predicted drag coefficient.
At low Reynolds number (Re < O(10 2 )), Figure 5b indicates that the percentage error in the drag coefficient predicted by the M96 parametrisation is strongly correlated with the area ratio of the analogues, as reported in previous work (Heymsfield and Westbrook, 2010;Westbrook and Sephton, 2017). The observed dependence on area ratio can yield significant errors when using the M96 parametrisation to predict fall speed; errors in C D can reach 60-80% for analogues with high and low area ratio (A r > 0.8 or A r < 0.3). Moreover, across the ensemble of measurements collected, the M96 parametrisation under-predicts the drag coefficient by approximately 20% on average, such that the parametrisation overestimates fall speed on average.
Effects associated with the onset of unsteady motions for Re ≳ 10 2 give rise to a number of key features in the C D (Re) data shown in Figure 5. Figure 5a shows that the parametrisation proposed by M96 continues to underestimate, on average, the drag coefficients and therefore overestimate fall speed. However, the more pronounced increase in C D for analogues in the plate-like category results in higher drag coefficients for these particles when compared with particles within the column-like, rosette and aggregate categories. Consequently, the percentage error in the M96 parametrisation exhibits a clear dependence on particle category, whereby larger percentage errors are observed for plate-like particles; on average, the M96 parametrisation under-predicts the drag coefficient of plate-like particles by 47% when Re ≳ 500. In contrast, on average, the M96 parametrisation under-predicts the drag coefficient of aggregate and rosette particles by just 10% when Re ≳ 500. Finally, note that Figure 5b shows that the error in the M96 parametrisation is only weakly dependent on the area ratio of the analogues when Re ≳ 500.

Böhm (1989) and Heymsfield and Westbrook (2010) parametrisations
Both B89 and HW10 proposed parametrisations of ice particle fall speeds of the form C D = A n r C 0 (1 + 0 Re −1∕2 ) 2 . These parametrisations assume that effects associated with the geometric shape of ice particles can be modelled through the use of the term A n r . If correct, rescaling the drag coefficients as C D A −n r would lead to a reduction in the scatter of the C D -Re data that was observed in Figure 5. Correspondingly, when evaluating the B89 and HW10 parametrisations, we consider two (related) aspects of the approach: (i) whether the re-scaling of C D reduces scatter within the data, and (ii) the accuracy of the fit of the curves C D = A n r C 0 (1 + 0 Re −1∕2 ) 2 used within the parametrisations to the re-scaled data. Measurements re-scaled as C D A 3∕4 r and C D A 1∕2 r in accordance with the B89 and HW10 parametrisations, respectively, are shown

1000< Re
A r F I G U R E 6 (a) Measurements of Re and re-scaled drag coefficient C D A 3∕4 r for ice particle analogues. Also shown is the parametrisation proposed by B89. (b) Percentage error when using the B89 parametrisation to predict the drag coefficient of ice particle analogues, where C B89 D denotes the predicted drag coefficient; data have been separated according to the Reynolds number at which the analogues fall (see panel labels for details). In each panel, analogues of different habits are shown separately in accordance with the categories shown in Figure 2 (see legend in (a), which applies to all panels). Different symbols are also used to show whether analogues fell in a steady ("+") or unsteady ("o") regime in Figures 6 and 7, in which the parametrisations are also shown.
A comparison of Figures 5a, 6a and 7a indicates that the proposed re-scaling of the drag coefficient successfully reduces the scatter in the C D -Re data at low Reynolds number (Re ≲ 200) but actually increases scatter at high Reynolds number (Re ≳ 200). The increase in scatter at high Re is more pronounced when the drag coefficient is re-scaled as C D A 3∕4 r (B89) than C D A 1∕2 r (HW10). These trends occur on account of the effects associated with the geometric shape of the particles identified in Section 4; That is, both parametrisations assume that effects associated with the geometric shape of the particles give rise to a series of approximately parallel C D -Re curves, whereby analogues with low area ratio exhibit higher drag coefficients than analogues with high area ratio (see Figure 13 for a visualisation). However, this model is only accurate for Re ≲ O(10 2 ), on account of the changes in shape of the drag coefficient curve that occur at higher Re, which we associated with the onset of unsteady motions. In particular, recall that, at Re ∼ 10 3 , we observe that analogues with low area ratio can exhibit slightly lower drag coefficients than analogues with higher area ratio; these instances give rise to the increase in scatter at high Re shown in Figures 6  and 7 when the drag coefficients are re-scaled as C D A −n r (n = −1/2, − 3/4).
Turning to the accuracy of the fit of the B89 parametrisation to the re-scaled C D -Re data, Figure 6a indicates that, on average, the parametrisation proposed by B89 significantly overestimates the drag coefficient of the analogues and correspondingly underestimates the terminal velocity. At low Reynolds number (1 ≤ Re ≤ 100), the B89 parametrisation overestimates the drag coefficient by 94% on average; the predicted drag coefficients are too large by a factor of almost 2. At high Reynolds number, the average fit of the parametrisation improves on account of the increase in C D associated with the onset of unsteady motions. However, the parametrisation still overestimates the drag coefficient by 55% on average for Re ≳ 1,000. In addition, for Re ≳ O(10 2 ), errors in the predicted drag coefficients are strongly correlated with the area ratio of the analogues, on account of the effects associated with the geometric shape of the particles that are prevalent when particles fall unsteadily (see Section 4), so errors associated with individual analogues far exceed average errors. For analogues with low area ratio (A r ≈ 0.2) we record errors exceeding 250%, equivalent to predicting drag coefficients that are too large by a factor of approximately 3.5. (For context, note that ice particle habits with very low area ratios include needles, stellar crystals, rosettes with thin branches, and aggregates of these sorts of particles.) Moreover, Figure 6b shows that errors at high Re (Re ≳ 100) are also correlated with particle category; again, this is due to effects associated with the geometric shape of the particles that are prevalent when particles fall unsteadily.
In contrast, Figure 7a indicates that, on average, the parametrisation proposed by HW10 provides a good fit to the data at low Reynolds numbers (1 ≤ Re ≤ 100). Over this range, the HW10 parametrisation overestimates the drag coefficient by 9% on average, whilst the root-mean-square 1000< R e A r F I G U R E 7 (a) Measurements of Re and re-scaled drag coefficient C D A 1∕2 r for ice particle analogues. Also shown is the parametrisation proposed by HW10. (b) Percentage error when using the HW10 parametrisation to predict the drag coefficient of ice particle analogues, where C HW10 D denotes the predicted drag coefficient; data have been separated according to the Reynolds number at which the analogues fall (see panel labels for details). In each panel, analogues of different habits are shown separately in accordance with the categories shown in Figure 2 (see legend in (a) which applies to all panels). Different symbols are also used to show whether analogues fell in a steady ("+") or unsteady ("o") regime (RMS) error associated with individual measurements of the drag coefficient is 27%. This compares favourably with errors in the M96 parametrisation, in which drag coefficients were underestimated by 21% on average, whilst the RMS error is 38%, and the B89 parametrisation (94% and 105%, respectively). However, the accuracy of the parametrisation reduces beyond Re ≈ 10 2 on account of the more complex effects of the geometric shape of the particles that occur when falling unsteadily. In particular, for Re ≳ O(10 2 ), errors in the predicted drag coefficients are strongly correlated with the area ratio of the analogues and with the particle categories, which occurs for the reasons explained above. For Re ≳ 10 2 , we observe that errors for irregular particles (rosettes and aggregates) are typically approximately 20% less than errors for planar particles. For Re > 10 3 , the HW10 parametrisation underestimates the drag coefficient by 26% on average, whilst the RMS error is 38%; for comparison, the corresponding errors for the M96 parametrisation are 27% and 34%, respectively.
In summary, the rescaling of the drag coefficients proposed by B89 and HW10 accounts, in part, for effects associated with the geometric shape of particles that are prevalent at low Reynolds number (1 ≤ Re ≤ 100). At these low Reynolds numbers, the HW10 parametrisation predicts drag coefficients with greater accuracy than the M96 and B89 parametrisations. However, the rescaling of C D does not account for the more complex effects associated with the geometric shape of ice particles that arise with the onset of unsteady motions at high Reynolds number (Re > 100). Consequently, at high Reynolds number, the accuracy with which the HW10 parametrisation predicts C D reduces to a level comparable to that of the M96 parametrisation, whilst significant errors arise in the B89 parametrisation.

Turbulent corrections
An increase in drag coefficient at high Reynolds number, associated with the onset of unsteady or turbulent flow in the wake of a falling particle, is predicted by the "turbulent corrections" proposed by Böhm (1992), M96 and Mitchell and Heymsfield (2005) (see Section 2.2). The increases in drag coefficient predicted by the turbulent corrections are shown in Figure 8, where C Dl denotes the drag coefficient predicted using parametrisations described above and C Dt denotes the "corrected" drag coefficients. Figure 8 indicates that the revised drag coefficients predicted after the Böhm (1992) and M96 turbulent corrections C Dt are essentially indistinguishable from the original values C Dl for Re ≲ 5 × 10 2 and exceed the original values by approximately 60% and 30%, respectively, for Re ≳ 5 × 10 3 . Recall, however, that the M96, B89 and HW10 fall speed parametrisations yield estimates of C Dl that decrease monotonically as Re increases. Consequently, when the Böhm (1992) and M96 turbulent corrections are applied, C Dt exhibits a modest increase (∼20% and 5%, respectively) in magnitude at Re ∼ 10 3 , before decreasing again for Re ≳ 5 × 10 3 . In contrast, Figure 8 indicates that the revised drag coefficients resulting from the Mitchell and Heymsfield (2005) correction C Dt exceed their original values C Dl by approximately 5% for Re ≲ 10 1 and that the proportion by which C Dt is larger than C Dl rapidly and consistently increases for Re ≳ 10 2 . Consequently, when applied to estimates of C Dl from the M96, B89 and HW10 parametrisations, the rapid increase in C Dt /C Dl for Re ≳ 10 2 results in values of C Dt that increase monotonically for Re ≳ 2 × 10 3 .
A comparison with the experimental data indicates that the turbulent corrections do not significantly improve the accuracy of fall speed parametrisations at high Re; That is, the turbulent corrections do not accurately represent the deviation from "steady flow behaviour" (see Section 4), which gives rise to the observed increase in drag coefficient. For example, our data show that the onset of unsteady particle motions, which is accompanied by the development of a unsteady wake, strongly depends on the geometric shape of an ice particle and, for some ice particle shapes, can occur at Reynolds numbers as low as Re ∼ 10 2 (see Figure 4). In addition, after the onset of unsteady motions, the trends in the C D (Re) data also strongly depend on the geometric shape of the ice particle; Figure 3 shows that more pronounced increases in C D occur for particles with high area ratio A r . We stress that these trends are absent from existing turbulent corrections, and we attribute the absence of these details to the limited range of high-accuracy data available to Böhm (1992), M96 and Mitchell and Heymsfield (2005); in particular, recall that results from quasi-spherical particles appear to have been used when deriving the turbulent corrections (see Section 2.2). Unfortunately, within the ice microphysics research community, the understanding of the underlying physical phenomena that give rise to these effects remains limited; following the onset of unsteady motions, the magnitude of the drag coefficient crucially depends on the structure of the wake (see Section 4), which remains poorly understood for natural ice particles. In particular, whilst the onset of unsteady particle motions at Re ≈ O(10 2 -10 3 ) is thought to be accompanied by the unsteady shedding of vortices in the wake, giving rise to a turbulent far-wake, the nature of the near-wake and boundary layers remains unknown. Regardless, it is evident that more sophisticated corrections are required to reflect the complex geometric effects on C D that are associated with the onset of unsteady particle motions at high Reynolds number.

REVISING THE FUNCTIONAL FORM OF ICE PARTICLE TERMINAL VELOCITY PARAMETRISATIONS
In Section 5, we demonstrated the need for a new ice particle terminal velocity parametrisation with sufficient flexibility to account for the complex geometric effects identified here. To this end, we propose that a revision to the functional form currently used for these parametrisations (i.e. Equation 4) is required. Here, we restrict our focus to modifying existing forms of Equation 4 and the use of A r as a parameter to describe the geometric shape of ice particles. We stress that more accurate parametrisations may be possible by adopting more complex shape parameters (such as sphericity, see Chhabra et al., 1999), but our aim here is to derive a parametrisation that can be seamlessly integrated into existing applications and which does not require the use of complex instrumentation to precisely determine the three-dimensional shape of falling ice particles.

Steady flow
As indicated in Section 2.2, both B89 and HW10 characterise the relationship between the drag coefficient and Reynolds number through the equation C D = A n r C 0 (1 + 0 Re −1∕2 ) 2 , which indicates that ice particles with low area ratios exhibit higher drag coefficients than ice particles with higher area ratios. This trend is consistent with the geometric effects observed in the experimental data when analogues were observed to fall steadily. However, the empirical values of n, C 0 and 0 proposed by B89 and HW10 were deduced by evaluating best fits of Equation 4 to limited datasets; the large uncertainties associated with these data have thus resulted in unreliable fits and unreliable estimates of n, C 0 and 0 to predict the drag coefficients of ice particle analogues for Re ≤ 100, plotted against A r , where C MW20 D denotes the predicted drag coefficient; different definitions of Reynolds number are used in (a) Re ≡ u t D∕ and (b)Re d ≡ u t d∕ . Best fits are plotted as a function n, which are shown for −1 ≤ n ≤ 0 at intervals of 0.05 (see colour bar, which applies to both plots). In both (a) and (b), the optimal value of n (n ≈ −0.6 and −0.4, respectively) is shown separately using the notation "− −". In (a), best fits in which n = 0 (as per M96), n = −0.5 (as per HW10) and n = −0.75 (as per B89) are also shown separately using the notation "− ⋅ −" (see for example Figure 6). In this section, we re-evaluate fits of the equation C D = A n r C 0 (1 + 0 Re −1∕2 ) 2 using the experimental data reported in Section 4 and determine optimal values of n, C 0 and 0 . We restrict our attention to low Reynolds number (Re ≤ 100) data, in which all particles were observed to fall steadily; unsteady affects are considered in Section 6.2.
Using the data reported in Section 4, best fits of Equation 4 to the experimental data were evaluated for different values of n in the range −1 ≤ n ≤ 0. For each fit, the residuals were analysed and the percentage errors evaluated, in accordance with the data shown in Figures 5b, 6b and 7b. In addition, average (RMS) errors and (second-order-polynomial) best fits characterising the average relationship between percentage error and A r , as a function of n, were also determined; a summary of these fits is shown in Figure 9a. Fits in which n = 0 (as per M96), n = −0.75 (as per B89) and n = −0.5 (as per HW10) are highlighted in Figure 9a. A comparison with Figures 5b, 6b and 7b indicates that the percentage errors associated with optimised forms of Equation 4 exhibit a similar dependence on A r as the M96, B89 and HW10 parametrisations, although the percentage errors are smaller on average in Figure 9a on account of using values of C 0 and 0 selected to best fit the experimental data. Figure 9a illustrates that the proper selection of n largely eliminates the relationship between the percentage error in the best-fit parametrisations and A r , or, in other words, accounts, in part, for the influence of geometric shape on the drag coefficient at low Reynolds number. The average relationship between the percentage error and A r is minimised when n ≈ −0.6, which lies between the values of n proposed by B89 and HW10. (When n ≈ −0.6, the best fit to the experimental data yields C D = 0.406A −0.6 r (1 + 6.32Re −1∕2 ) 2 ). Accordingly, the average (RMS) error associated with the best fits to the drag coefficient data (for Re ≤ 100) is also minimised when n ≈ −0.6, at 23%; for comparison, the average (RMS) error for the HW10 parametrisation was 27%. Consequently, optimising the value of n to account for the influence of geometric shape on the drag coefficient yields a small but not insignificant increase in the accuracy of ice particle terminal velocity parametrisations.
However, at this point it is instructive to re-analyse assumptions implicit within this parametrisation. In particular, so far we have adopted the conventional definition of Reynolds number Re ≡ u t D∕ , where the characteristic length scale D is defined as the maximum span of the analogue (normal to the fall direction). To the best of our knowledge, a clear and rigorous physical justification for this definition of the characteristic length scale is not present in the literature. Instead, this convention appears to be in use for practical reasons, as a geometric property of an ice particle that is easy to measure. Next we consider whether alternative definitions of the characteristic length scale may facilitate the derivation of more accurate fall speed parametrisations; indeed, this approach was used by B89 (see Section 2.2). We note that, for related problems in other fields, the characteristic length scale is commonly defined as the diameter of a sphere of equivalent volume d ≡ (6V∕ ) 1∕3 (see for example Dietrich, 1982;Tran-Cong et al., 2004;Bagheri and Bonadonna, 2016), such that the Reynolds number is defined as Re d ≡ u t d∕ . The revised characteristic length scale d also gives rise to a revised Best number Be d (Be d = C D Re 2 d ) as This length scale is an intuitive choice for irregular particles that are quasi-spherical in shape, but is less intuitive for use with complex geometric shapes exhibited by atmospheric ice particles. Regardless, as we shall see, the use of this revised characteristic length scale facilitates a better collapse of the C D (Re d ) data and thus more accurate fall speed parametrisations.
We now characterise the relationship between the drag coefficient and Reynolds number using a revised form of Equation 4 (i.e. using Re d in place of Re) The corresponding Re d -Be d relationship is We note that the alternative definition of the characteristic length scale employed in these parametrisations can be seamlessly integrated into applications employing existing fall speed parametrisations, since d is evaluated directly from the particle volume (or mass), which is already used as an input in the estimation of the Best number.
The process used previously to evaluate fits of Equation 4 is now used to evaluate Equation 7; a summary of the fits characterising the average relationship between the percentage error and A r as a function of n is shown in Figure 9b. Figure 9b again illustrates that the proper selection of n reduces the dependence of the percentage error of the best-fit parametrisations on A r ; That is, when adopting the revised characterise length scale d, ice particles with low area ratio still exhibit higher drag coefficients than ice particles with higher area ratio. Consequently, qualitatively similar trends are observed in both Reynolds number Re is used). Consequently, adopting a revised characteristic length scale d in the definition of the Reynolds number Re d results in an improved accuracy. We henceforth adopt Equation 7 (with n = −0.4, C 0 = 0.498 and 0 = 3.71) as a revised functional form of ice particle fall speed parametrisation for use at low Reynolds number; this form is modified in Section 6.2 to account for more complex geometric effects associated with the onset of unsteady motions in free fall at higher Reynolds number.

Unsteady flow
As indicated in Section 4, the onset of unsteady motions as an ice particle falls, which is coupled with the onset of unsteady flow in the particle wake, gives rise to a change in shape of the C D -Re curve, whereby the drag coefficient deviates from the trend observed at low Reynolds number. This change in drag coefficient can be incorporated into fall speed parametrisations by amending the form of Equation 7 to the form where C D,steady ≡ C 0 (1 + 0 Re −1∕2 d ) 2 and h(Re d ) is a cut-off function that acts to gradually alter C D A 0.4 r from the value C D, steady (i.e. that predicted using the function for steady flow derived in Section 6.1) towards some value C 1 as the Reynolds number increases. This functional form has similar characteristics to the curves in Figure 1; with proper selection of C 0 , 0 , h(Re d ) and C 1 , Equation 9 can produce a curve indistinguishable from these results for spheres for Reynolds numbers up to Re d ≈ O(10 5 ).
To facilitate the selection of h(Re d ) and C 1 , we assume that the C D -Re d curves for the ice particle analogues retain key features of the C D -Re d curve for spherical particles. Principally we assume that, for any given analogue, the drag coefficient tends to a (geometry-specific) constant value at sufficiently high Reynolds number (excluding drag crisis). This constant value is C 1 . We note that only some of the data reported in Section 4 exhibit this trend. However, we attribute this to the limited range of Reynolds number considered and assume that all analogues would have exhibited approximately constant drag coefficients in experiments at higher Reynolds number, as per the data of, for example, Becker (1959), Haider and Levenspiel (1989) and Seyed-Ahmadi and Wachs (2019). We expect C 1 to be influenced by the geometric shape of the particles. Specifically we assume that C 1 depends on A r . An appropriate form of C 1 can be deduced by analysing C D A 0.4 r data at the limit of high Reynolds number (in the limit of high Reynolds number, Equation 9 simplifies to C D A 0.4 r = C 1 ). We assume that insight into the form of this function can be deduced using drag coefficient data in the range 500 < Re d , as shown in Figure 10. Figure 10 reveals two distinct trends: one for analogues in the plate-like category and another for analogues in the column-like, rosettes and aggregates categories. In each case, the data are consistent with the trend C D A 0.4 r ≈ a 0 + a 1 A r , where a 0 and a 1 are constants, so we assume that C 1 can be approximated as C 1 = a 0 + a 1 A r . However, the two trends observed in Figure 10 indicate that distinct values of a 0 and a 1 are required for particles with plate-like geometries and particles with columnar or irregular (rosette and aggregate) geometries. We note the presence of four anomalous values of C D A 0.4 r in the rosette category at A r ≈ 0.5 and A r ≈ 0.57; these correspond to the anomalous measurements of C D discussed in Section 4.
However, it is not clear how best to represent the cut-off function h(Re d ). To offer flexibility in this approach, we adopt a cut-off function of the form e −(Re d ∕Re T ) m , where m denotes a constant that controls how sharply the transition in Equation 9 from the curve C D, steady to the constant C 1 occurs and Re T is a characteristic Reynolds number at which this transition occurs. We may expect Re T to be related to the critical Reynolds number Re c at which we observe the onset of unsteady motions. Recall however that the relationship between Re c and the geometric shape of ice particles cannot be reliably approximated (to first order) as a function of A r . Consequently, since we are neglecting the use of other, more complex shape parameters, we adopt a fixed value of Re T within the cut-off function h(Re d ) used in Equation 9.
The values of m, Re T and the parameters in C 1 = a 0 + a 1 A r can be deduced by using non-linear regression on the C D A 0.4 r -Re d data. In accordance with the results shown in Figure 10, regression is applied separately to data in the plate-like category and data in the rosette, aggregate and column-like categories. Analysis of the residuals of the regression fits indicates that the average percentage error is minimised when m = 1.6 (i.e. h(Re d ) = e −(Re d ∕Re T ) 1.6 ). The resulting fits yield Re T = 183 and C 1 = 0.40 + 1.66A r for data in the plate-like category and Re T = 142 and C 1 = 0.38 + 0.67A r for data in the column-like, rosette and aggregate categories. A comparison of the regression fits with the experimental data is shown in Figure 11. We note that the values of Re T used are of the same order of magnitude as the estimates of Re c shown in Figure 4. We also note that the function C 1 evaluated from these regression fits is similar to, but not identical to, the function that would result from applying linear regression to data reported in Figure 10.
The percentage errors that occur when using Equation 9 to predict the drag coefficients of the ice particle analogues are shown in Figure 12, where C MW20 D denotes the predicted drag coefficients. A comparison of the average errors in the drag coefficient predicted by Equation 9 and by HW10 is presented in Table 1. Figure 12 indicates that the percentage error when using Equation 9 to predict drag coefficients is not correlated with A r , or the particle categories, at either low or high Reynolds number. Crucially, this illustrates that Equation 9 reliably approximates key aspects of the underlying physical trends. However, Figure 12 does indicate that, in some cases, errors associated with individual measurements of C D remain large, as drag coefficients of individual particles can by be under-or over-predicted by factors as large as 2. On the other hand, Table 1 illustrates that  Equation 9 yields a significant improvement in accuracy with which the drag coefficient is predicted on average, relative to the HW10 parametrisation; Equation 9 overestimates the measured drag coefficients by less than 3% on average, whilst the RMS error is less than 20%, which compares favourably with the equivalent statistics for the HW10 parametrisation (7.8% and 31.9%, respectively). As expected, the most significant improvements in accuracy occur at high Reynolds number (Re ≳ 10 2 ), although modest improvements in accuracy are also observed for Re ≤ 100.
We note that, whilst the change in shape of the C D -Re curve at the onset of unsteady motions, which we have considered in this section, is thought to be related to wake instability (see Section 4), the exact physical mechanism that controls this remains unclear. Indeed, given the slightly different shape of fits evaluated for data in the plate-like category relative to data in the column-like, F I G U R E 11 Measurements of Re d and C D A 0.4 r for ice particle analogues in the (a) plate-like and (b) column-like, rosette and aggregate categories (see Figure 2). Different symbols are used to show whether analogues fell in a steady ("+") or unsteady ("o") regime. In each case, the colour of the symbols denotes the area ratio of the analogues, as shown by the colour bar, which applies to both plots. Also shown are Equations 7 ("-") and 9 ("--"), in which (a)  F I G U R E 12 Percentage error when using Equation 9 to predict the drag coefficient of ice particle analogues; data have been separated according to the Reynolds number at which the analogues fall (see panel labels for details). In each panel, analogues of different habits are shown separately in accordance with the categories shown in Figure 2 (see legend in (c), which applies to all panels). Different symbols are also used to show whether analogues fell in a steady ("+") or unsteady ("o") regime TA B L E 1 Summary of average percentage errors when using Equation 9 and the Heymsfield and Westbrook (2010) parametrisation to predict the drag coefficient of ice particle analogues. A negative average error indicates that the parametrisation under-predicts C D on average and correspondingly over-predicts terminal velocity ("--"), suitable for use with particles that fall steadily, and given by Equation 9 ("-"), which features a correction to account for particles that fall unsteadily; (a) plate-like ice particles (Re T = 183 and C 1 = 0.4 + 1.66A r ), (b) other ice particles (Re T = 142 and C 1 = 0.38 + 0.67A r ). Curves are shown as a function of A r over the range 0.2 ≤ A r ≤ 0.8 at intervals of 0.1

Predicting ice particle fall speeds
To predict the fall speed of natural ice particles, Equation 9 must be expressed in terms of the Best number, Be d , which can be estimated directly from knowledge of a particle's shape and size and the fluid in which it is falling. By using the relation Be d ≡ C D Re 2 d to eliminate C D from equation 9, and re-arranging, one obtains valid for 1 ≲ Re d ≲ 10 3 and 0.2 ≲ A r ≲ 0.83. Thus, in principle, ice particle fall speeds can then be estimated in three steps: firstly estimating Be d using Equation 6, then estimating Re d by solving Equation 10 and finally estimating u t using the expression u t = Re d ∕d. For completeness, we note that the particle volume V that is used to estimate Be d (Equation 6) and d ≡ (6V∕ ) 1∕3 can be straightforwardly determined from the particle mass using the relation V = m∕ ice , where ice denotes the density of ice (i.e. 917 kg m −3 ). In addition, we note that, where direct measurements of particle mass m and projected area A are not available, estimates of m and A can be obtained by using existing m-D and A-D relations. However, the approach described above requires Re d to be evaluated numerically, since Equation 10 cannot be inverted analytically. To facilitate the practical implementation of this parametrisation, we report an alternative expression that can be used to estimate Re d . That is, since Equation 9 simplifies to C D A 0.4 r = C D,steady at very low Reynolds number and simplifies to C D A 0.4 r = C 1 at very high Reynolds number, we obtain the expressions Re d =  Figure 13 as a function of A r , alongside the relationships described by Equation 9, which includes the correction to account for particles that fall unsteadily. Two sets of curves are shown in Figure 13 in accordance with the results of Section 6.2, which use different values of Re T and C 1 . The appropriate values of Re T and C 1 are Re T = 142 and C 1 = 0.38 + 0.67A r for all ice particles, except where there is evidence that ice particles exhibit a plate-like geometry, for which Re T = 183 and C 1 = 0.40 + 1.66A r . For the avoidance of doubt, we consider the following particles from the classification system proposed by Kikuchi et al. (2013) to exhibit a plate-like geometry: P1-P8, CP1, CP2a/b, CP3, CP4a, R1b/c and R2b/c, which includes plane crystals (pristine and rimed), assemblages of plane-type crystals, capped columns/bullets and crossed plates. Figure 13 shows that each set of curves obtained using Equations 7 and 9 are essentially indistinguishable from each other when Be d ≲ 3 × 10 3 (i.e. Re d ≲ 50). Consequently, when Be d ≲ 3 × 10 3 , the classification of particles using the guidance detailed above can be disregarded and the original Equation 7 used in place of Equations 9 and 10. Figure 13b also shows that, for irregular ice particles, the drag coefficient is approximately constant and independent of A r at Be d ∼ 10 5 , taking a value of C D ≈ 1.
Finally, recall that the approach used in this study exploits dynamic similarity to apply the results obtained from experiments in which ice particle analogues fall in quiescent fluid to natural ice particles falling in the atmosphere. Consequently, the relationships described by Equations 7 and 9 are valid in quiescent (i.e. non-turbulent) regions of clouds. However, the accuracy of these relationships in turbulent regions of clouds remains unclear. In previous studies, atmospheric turbulence has been reported to both enhance and reduce the average settling velocity of ice particles (see for example Garrett and Yuter, 2014). The enhancement or reduction in average settling velocity of an ice particle as a result of atmospheric turbulence is thought to depend, in part, upon the relative magnitude of the terminal velocity of the ice particle and the root-mean-square velocity of the turbulent flow (Good et al., 2012). Unfortunately, the influence of particle geometry on this interaction remains poorly understood. Consequently, further work is required to understand the complex interaction between atmospheric turbulence and irregularly shaped ice particles before the accuracy of the relationships described by Equations 7 and 9 in turbulent clouds can be determined.

CONCLUSIONS
We have reported comprehensive measurements from an extensive experimental study investigating the drag coefficients of 3D-printed analogues of ice particles. Drag coefficients were determined through the use of the novel algorithm TRAIL (see part 1) that reconstructs the time-resolved trajectory and orientation of the analogues when falling in a quiescent viscous liquid. In total, we report results from 844 experiments, comprising measurements of particles with 118 different geometries each falling at a range of Reynolds numbers between Re ∼ O(10 0 ) and Re ∼ O(10 3 ). Using these measurements, we identify, for the first time, complex geometric effects on the drag coefficient that arise for Re ≳ O(10 2 ), which are associated with the onset of unsteady motions, such as oscillating or fluttering motions, during free fall. The onset of unsteady motions is thought to occur as a result of the formation of instabilities in the wake of the falling particles. However, the critical Reynolds number at which this occurs exhibits substantial variability; some analogues were observed to exhibit unsteady motions at Re < 100, whilst others were observed to fall steadily at Re > 1, 000. We find that the critical Reynolds number cannot be reliably approximated in terms of simple shape factors, such as A r . We attribute this result to the complex dynamics that give rise to wake instability, which are related to the intensity of the vorticity generated in the wake of the falling particles (Auguste et al., 2013).
Using the experimental data, we assess the accuracy of the fall speed parametrisations proposed by Mitchell (1996), Böhm (1989) and Heymsfield and Westbrook (2010) and the turbulent corrections proposed by Böhm (1992), Mitchell (1996) and Mitchell and Heymsfield (2005). We find that, at low Reynolds number (Re ≲ 100), the Heymsfield and Westbrook (2010) parametrisation is more accurate than either the Böhm (1989) or Mitchell (1996) parametrisation. In this case, the average error when using the parametrisation to predict the measured drag coefficients is approximately 27% for Re < 100, although in some cases the drag coefficients are underor over-predicted by factors as large as 2. Larger average errors are observed at higher Reynolds numbers on account of the more complex geometric effects that occur when particles fall unsteadily. We also find that the turbulent corrections proposed in the literature do not significantly improve the accuracy of the fall speed parametrisations at high Reynolds number. Principally, this is because our data indicate that the change in drag coefficient associated with the onset of unsteady motions depends crucially on the geometric shape of the ice particle.
We present a new ice particle fall speed parametrisation using our newly derived insight into geometric effects on the drag coefficient. At low Reynolds number, the parametrisation is similar to that proposed by Heymsfield and Westbrook (2010). However, the parametrisations differ significantly for Re ≳ 50, whereby the new parametrisation predicts drag coefficients that tend to a constant C D ≈ a 0 A −0.4 r + a 1 A 0.6 r at Re ∼ 10 3 . Different values of a 0 and a 1 are reported for analogues that exhibit plate-like geometries, when compared with analogues with irregular geometries such as aggregates. We note that these functions may not be suitable for use with hail or graupel on account of the different geometric shapes associated with these types of ice particles. This approach yields a significant improvement in accuracy with which the drag coefficient can be predicted at high Reynolds number, relative to the Heymsfield and Westbrook (2010) parametrisation, without requiring detailed knowledge of the geometric shape of the falling particles. The average error when using the parametrisation to predict the measured drag coefficients is less than 20%.