TRAIL: A novel approach for studying the aerodynamics of ice particles

A novel experimental approach for studying the aerodynamics of ice particles is presented. The Trajectory Reconstruction Algorithm implemented through Image anaLysis (TRAIL) produces digital reconstructions of the trajectory and orientation of 3D‐printed analogues of ice particles that fall through a quiescent viscous liquid. Data extracted from this analysis can be used to test fall speed parametrisations and investigate the preferred orientation of analogues with complex irregular geometries. Experiments using thin circular discs are used to validate this approach. Measurements from a case study in which an analogue of an aggregate ice particle is used are also reported. Data extracted from the analysis clearly demonstrate important, but poorly understood, aspects of the aerodynamics of atmospheric ice particles. For example, the orientation that the analogue adopts in free fall is shown to depend not only on the geometric shape of the particle, but also on the Reynolds number at which it falls. In addition, measurements of the drag coefficient indicate that the accuracy of fall speed parametrisations reduces at high Reynolds number due to the onset of unsteady motions.


INTRODUCTION
The aerodynamics of ice particles is of fundamental interest as it influences almost every aspect of their microphysical evolution, including vapour growth, riming, aggregation, evaporation and melting. Consequently, the aerodynamics of ice particles is intimately linked with a number of problems within atmospheric physics. 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 . Assumptions related to the fall speed of ice particles are also required when estimating particle properties from Doppler radar measurements (Matrosov, 2011) and from surface measurements of snow particle size and fall speed (von Lerber et al., 2017). In addition, a detailed understanding of the orientations that ice particles adopt is required to correctly interpret polarimetric remote-sensing observations of clouds (Matrosov et al., 2017;Miao et al., 2003;Westbrook et al., 2010). However, the aerodynamics of ice particles remains poorly understood.
A range of approaches have previously been employed to study the aerodynamics of ice particles. Perhaps the most commonly used approach has been through surface observation of natural ice particles (see e.g. Zikmunda, 1972;Locatelli and Hobbs, 1974;Kajikawa, 1992). Over time, studies have progressed from manually recording the properties of individual ice particles to using automated in situ instrumentation, including the Particle Video Imager/Particle Imaging Package (Newman et al., 2009), 2D video disdrometer (Kruger and Krajewski, 2002;Schönhuber et al., 2007) and the Multi-Angle Snowflake Camera (Garrett et al., 2012). However, despite technological innovation that now facilitates the simultaneous acquisition of multiple images of ice particles from several different perspectives, alongside automated measurements of fall speed, existing instrumentation does not facilitate measurements of some key parameters. In particular, challenges are encountered in determining the mass and 3D geometric shape of falling ice particles, as well as in characterising the trajectory and orientation of the particles.
Another long-standing approach has been the use of ice particle "analogues" in laboratory experiments, where scaled-up models of ice particles are dropped in tanks of viscous liquids so that they fall at the same Reynolds numbers as ice particles falling in the atmosphere (see e.g. Jayaweera, 1972;List and Schemenauer, 1971). The advantage of this approach is that very accurate measurements can be acquired under precisely controlled conditions, whilst the shape and size of falling particles are also known. However, studies adopting this approach have been limited to analogues with simple, idealised geometries since traditional fabrication techniques precluded the fabrication of analogues with complex geometries. Laboratory experiments studying the aerodynamics of ice crystals grown in "cloud chambers" (e.g. Kajikawa, 1973;Bürgesser et al., 2016) can also provide reliable data on the aerodynamics of simple ice crystals. However, the range of particle sizes and shapes for which data is available is also limited; cloud chambers are unable to fully replicate the complex accretion processes that an ice particle experiences within the atmosphere, and consequently many natural ice crystal habits cannot be generated within cloud chambers, whilst maximum particle sizes are constrained by growth time limitations.
More recently, detailed data have been reported in studies employing computational fluid dynamics (see e.g. Cheng et al., 2015;Giovacchini, 2017;Nettesheim and Wang, 2018), in which the motion of a falling particle, and the surrounding fluid, are determined by solving the equations of motion of the system. A distinct advantage of this approach is the additional physical insight that can be derived from data describing the flow field around the falling particle; for example, these data can be used to investigate the enhanced growth rates experienced by precipitating ice particles (Ji and Wang, 1999;Wang and Ji, 2000). However, major computational resources are required for this approach, since the governing equations are solved on a computational grid that must be sufficiently resolved to precisely characterise the geometry of the ice particle and the surrounding flow field. Consequently, it is only feasible to sample a limited range of shapes in these numerical experiments; to date, studies adopting this approach have been limited to particles with simple, idealised geometries.
In summary, what is missing from previous work is an approach that combines the use of particles with complex irregular geometries with measurement techniques that enable key parameters that control the aerodynamics of the ice particles (such as mass and geometric shape) to be precisely determined. As such, there is a lack of detailed reliable data available on the aerodynamics of ice particles with irregular geometries, such as polycrystals or aggregates, 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). Westbrook and Sephton (2017) recently revived the analogue method for studying the aerodynamics of ice particles by employing modern 3D-printing techniques to fabricate analogues. This innovation overcomes the fundamental limitation inherent in previous analogue experiments by manufacturing ice particle analogues with complex geometries. In this study, we report a novel experimental approach that further extends the analogue method by digitally reconstructing the trajectory of falling analogues. This method, denoted the Trajectory Reconstruction Algorithm implemented through Image anaLysis (TRAIL), results in a time-resolved reconstruction of the location and orientation of ice particle analogues in free fall. Here, we demonstrate this approach using experiments in which scaled-up models of ice particles fall in quiescent viscous liquids, but we note that the technique could also be applied to particles falling in air. Crucially, TRAIL works not only for analogues with simple geometries, but also for analogues with complex irregular shapes. The detailed data that can be extracted from these digital reconstructions facilitate the measurement of drag coefficients, which are used in parametrisations of ice particle fall speeds (see Section 2.1), with superior accuracy when compared with previous studies. Moreover, the experimental approach also facilitates the study of other important aspects of the aerodynamics of ice particles, such as the preferred orientation of particles with complex irregular geometries.
The experimental approach underpinning TRAIL is described in Section 2. Validation of the experimental approach is reported in Section 3, whilst brief results from a case study investigating the aerodynamics of an irregular aggregate ice particle are reported in Section 4. Conclusions are presented in Section 5. Extensive results from the application of this technique to a large ensemble of analogues representing complex ice particles are reported separately in a companion paper.

Principles underpinning the analogue method
Data collected from experiments in which scaled-up models of ice particles fall in quiescent viscous liquid can be used to understand the aerodynamics of real ice particles by exploiting 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, relationships between dimensionless quantities are the same for ice particle analogues falling in viscous liquids and natural ice particles falling in the atmosphere. This result can be exploited to study a range of problems related to the aerodynamics of ice particles. For example, 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 where u t denotes the terminal velocity of an ice particle with maximum dimension 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, these 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. That is, where is the dynamic viscosity of the fluid. Since C D , Be and Re are all dimensionless quantities, the relationships between these parameters can be evaluated using data acquired from analogue experiments through the principle of dynamic similarity. However, care should be taken when interpreting data collected from analogues falling in viscous fluids since the conditions for dynamic similarity are only strictly satisfied when the particles 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. For reference, we note that the onset of unsteady motions occurs between Re ∼ 10 2 and 10 3 ; typically discs and plate-like particles exhibit the onset of unsteady motions at a lower Reynolds number than particles with complex geometries such as aggregates (see companion paper). The influence of the density difference in this regime is poorly understood. However, since here we are principally concerned with describing the novel TRAIL methodology, which can be applied to both scaled-up models of ice particles falling in viscous liquids and life-size models of ice particles falling in air, we defer a detailed discussion of the influence of unsteady motions on the aerodynamics of ice particles to subsequent articles in which detailed results will be reported.

Fabrication of ice analogues
Analogues of ice particles have been manufactured using a stereolithography (SLA) fabrication technique (Hull, 1986;Huang et al., 2013;Bikas et al., 2016), which builds plastic models up from a series of thin layers. Each layer is produced by using a laser to selectively cure regions of a liquid photopolymer resin. Initial layers of the analogues are formed against a metallic plate that rests within a bath filled with the photopolymer resin. After each layer is cured, the plate is raised upwards and an adjoining layer of the analogue is cured beneath the existing structure. Complex 3D structures are fabricated through the use of additional support structures that are removed at the end of the fabrication process.
Here a Form 2 3D printer (https://formlabs.com/) is used, featuring a minimum layer thickness of 25 m and a laser spot size of 140 m. The geometries of ice analogues are designed using computer-aided design packages and are converted into a set of instructions for the 3D printer using open-source software (Preform, https://formlabs.com/software/). The high precision of this technique, and the use of support structures, facilitates the fabrication of analogues with more complex geometries than used in previous studies based on ice particle analogues (such as Kajikawa, 1971;Westbrook and Sephton, 2017).

Sedimentation apparatus
An overview of the experimental set-up used to demonstrate the technique is shown in Figure 1. Experiments are conducted in a transparent acrylic tank with internal dimensions of 0.4 × 0.4 ×, 1.8 m. The tank is filled with uniform mixtures of water and glycerol, in which the volume fraction of glycerol is set between 0% and approximately 50%, to a total liquid depth of approximately 1.75 m. A thermocouple is used to confirm the absence of temperature gradients within the tank, which could otherwise induce large-scale circulations or density stratification. The density and temperature of the mixtures are also measured before each experiment; these measurements are used to estimate the viscosity of the mixture (Cheng, 2008;Volk and Kähler, 2018). During experiments, ice particle analogues are submerged close to the top of the tank and released with a random orientation; care is taken to remove any air bubbles from the analogues before they are released. Qualitative observations indicate that, following a short transient phase (typically over a vertical extent of less than 0.3 m), the subsequent fall of the analogues is independent of their orientation upon release. That is, the analogues either (a) exhibit a clear preferred orientation in free fall that is insensitive to the release conditions, or (b) exhibit a chaotic trajectory, the details of which are unreproducible even should the particle be released from a nominally identical orientation. To ensure quiescent conditions during experiments, the release of analogues is separated by a time period of at least 15 min.
The trajectory of falling analogues is recorded using a system of three synchronised cameras (Jai Go-5000M, 12.8 × 10.2 mm CMOS sensor, 35 and 50 mm focal length lenses at f # ∼ 5.6) positioned to have orthogonal views of the tank (Figure 1). Each camera records the fall of the particle through a region approximately 0.2 × 0.2 × 0.2 m in size, which is located centrally on-plan within the tank, approximately 1.5 m below the surface of the fluid.
Light-emitting diode (LED) panels are used to provide backlit illumination for the cameras (Figure 1), so each camera records a projected image of the particle (at 2560 × 2048 pixel resolution). Each camera's capture is synchronised using a signal generator; the frame rate of the capture is varied (between 10 and 50 frames per second) according to the fall speed of the particle.
We note that Esteban et al. (2019) reported that the dynamics of falling particles is influenced by the bottom of a tank when they fall to within a distance from the tank base of approximately 2D; in the current approach, this distance is typically equal to approximately 4 cm above the tank base (by comparison, the bottom extent of the recorded region is approximately 10 cm above the tank base). Analysis of trajectories reconstructed using the TRAIL methodology shows no evidence that the dynamics of falling particles changes as particles fall from the top to the bottom of the sample volume. Consequently, we conclude that hydrodynamic interactions between the falling analogues and the tank base are negligible within the recorded region. However, sedimentation apparatus are also known to suffer from "wall effects" at low Reynolds numbers, whereby the motion of a particle in free fall is retarded by the influence of the tank side walls (Brown and Lawler, 2003;Fidleris and Whitmore, 1961). The relative influence of the wall effect on the motion of a particle decreases either as the Reynolds number increases or the size of the particle decreases relative to the distance between the tank walls. In the current apparatus, the wall effect is thought to reduce the terminal velocity of particles by less than 1% at the lowest Reynolds number used (Fidleris and Whitmore, 1961). Consequently, we conclude that the hydrodynamic interactions between the falling analogues and the tank side walls are also negligible.

Co-ordinate systems
Since irregular particles in free fall exhibit both linear and angular displacements, a system of three coordinate systems are used to describe the resulting motion: (a) a (fixed) global reference frame x = (x 1 , x 2 , x 3 ) with origin at the centre of the region recorded by the cameras (shown in Figure 1), (b) a co-moving reference frame ) with origin at the particle's centre of mass and axes aligned with that of the global reference frame (also shown in Figure 1) and (c) a co-rotational reference frame with origin also at the particle's centre of mass but axes aligned with the principal axes of the particle (see for example Morin, 2008, p. 383). Here, the vertical axis is denoted x 2 , which increases in the direction of gravity. The orientation of the co-rotational reference frame x ′′ relative to the fixed global reference frame x is F I G U R E 1 Sketches and images showing key components of the experimental set-up, including the positioning of the cameras and the region of interest that they record. (a) Side and (b) plan view of the set-up. Also shown are the coordinate directions (x 1 , x 2 , x 3 ): note x 2 is the vertical axis. (c, d) Overview of the set-up; green dye has been used to make the fluid in the tank visible in the photographs, but this would otherwise not be present calculated via the co-moving reference frame, x ′ . That is, Euler angles, , , are used to describe the orientation of the co-rotational reference frame x ′′ through a series of three rotations applied to the co-moving reference frame x ′ . The Euler angles can be defined in a number of different ways; here we adopt a convention which consists of a rotation of about the positive x 3 ′ -axis, followed by a rotation of about the displaced positive x 2 ′ -axis, and subsequently a rotation of about the displaced positive x 1 ′ -axis, as shown in Figure 2.

Calibration
The trajectories of falling analogues are reconstructed using a custom algorithm described in Section 2.6. Central to this algorithm is the use of a function X = f (x) that describes the projection operation used to map between the global coordinate reference frame x = (x 1 , x 2 , x 3 ) and the pixel coordinates X = (X 1 , X 2 ) used in the images captured by the digital cameras. In effect, this projection operation describes the magnification that relates how distances in the sample volume translate to numbers of pixels in the acquired digital images. However, the projection operation also accounts for image distortion, which arises from the imaging optics and from refraction by the tank and the fluid within it. To describe this projection operation, we adopt a modified form of the mapping function proposed by Soloff et al. (1997), which takes the form where a i denotes a vector of six constants (two per camera). This mapping function was originally proposed as F I G U R E 2 Visualisation of Euler angle convention used; an example of the orientation of an analogue subject to rotations of ( , , ) = (70, 20, 30) is shown in (a-d). The analogue is shown (a) with the co-moving and co-rotating axes in alignment, (b) after rotation of about the positive x ′′ 3 -axis ("-"), (c) after rotation of about the displaced positive x ′′ 2 -axis ("-") and (d) after rotation of about the displaced positive x ′′ 1 -axis ("-"). In each panel, the orientation of the co-rotating axes are shown both before ("-") and after ("--") the corresponding rotation a simple technique for camera calibration, as part of the particle image velocimetry technique, that facilitated the reconstruction of three-dimensional velocity vector data from a stereoscopic pair of digital cameras.
During the calibration procedure for each camera, the pixel coordinates X that correspond to known reference points within the global coordinate reference frame x are determined by imaging a thin grid positioned in a number of different planes within the tank, as shown in Figure 3. A two-step multi-variate non-linear regression process is then used to determine the coefficients of the mapping function f (x) (Raffel et al., 2007). First, an estimate of the first-order terms of Equation 4 is determined by applying multi-variate regression to a truncated form of the mapping function f (x) ≈ a 0 + a 1 x 1 + a 2 x 2 + a 3 x 3 . This truncated form of Equation 4 physically represents a projection operation in which a single fixed magnification is assumed (i.e. the image plane exists at a fixed distance from the camera) and in which image distortion is neglected. The estimates of a 0 , a 1 , a 2 , a 3 are then used as initial guesses in multi-variate regression applied to Equation 4; the final values of these leading-order terms differ from their original estimates by no more than 2%.
A limitation of this calibration technique is that, if insufficient data are available for calibration, the use of third-order terms can result in undesirable oscillations in the mapping function near the edge of the images (Raffel et al., 2007). Initial testing during camera calibration indicated that the influence of third-order terms in the mapping function was small, but not insignificant. As such, only co-varying third-order terms of the form x i x i x j (no summation on repeated indices) are retained within Equation 4 in order to avoid these oscillations. Terms of the form x 3 1 , x 3 2 , x 3 3 and x 1 x 2 x 3 , which were originally included in the mapping function proposed by Soloff et al. (1997), were discarded. Coefficients of these terms were typically found to be an order of magnitude smaller than the third-order co-varying terms, so it was assumed that these terms were a result of noise within the calibration procedure. Analysis of the residuals of the final regression fit indicates that, on average, Equation 4 is accurate to sub-pixel precision, which corresponds to a physical distance of order 0.1 mm in size.

2.6
Summary of algorithm used to reconstruct particle trajectories The time-resolved trajectory and orientation of falling analogues are reconstructed from images of the analogue particle in free fall using a technique similar to that detailed by Marcus et al. (2014). Crucially however, unlike previous work, the algorithm described here is suitable for use with analogues of arbitrary geometry. Details of this technique are provided below.
Recall that, during an experiment, the free fall of a particle is recorded using a system of three synchronised digital cameras. For each recorded image, the falling particle analogue is identified and isolated from the background by applying a threshold to the grayscale pixel values (see for example Otsu, 1979), resulting in binary images. In this process a template image of the background is subtracted from each particle image to aid the identification of the falling analogue.
An approximate estimate of the trajectory of the analogue is then determined using Equation 4. First, the co-ordinates of the geometric centre of the particle's projection within the binary images, denoted X GC , are determined. Here the geometric centre is defined as the mid-point of the minimum bounding rectangle (with edges aligned with the axes of X) that encloses the projection. Next, combinations of three components of X GC F I G U R E 3 Schematic showing the camera calibration procedure used to determine the mapping function between the global reference frame x and the pixel coordinate system X; to calibrate the camera shown, images are acquired of the calibration grid positioned in a series of vertical (x i -x k ) planes. During the calibration procedure for each camera, x i , x j and x k correspond to different axes of the global reference frame (x 1 , x 2 , x 3 ) according to the location and orientation of the camera within the experimental apparatus (see Figure 1) are used to invert Equation 4 and obtain the corresponding location in 3D space, which we denote x GC . Note, however, that this inversion process is overdetermined; co-ordinates of X GC are required from only two cameras to invert Equation 4. The use of three cameras in the current experiments facilitates independent estimates of the particle's location, which are computed by using different combinations of particle images. A comparison of these independent estimates therefore indicates whether the location of the particles has been reliably determined. This comparison indicates that reconstructions in this apparatus are accurate to better than ±0.5 mm; for context, we note that particles are typically of the order 20 mm in size and are recorded falling across a vertical distance of approximately 200 mm. The preliminary estimate of trajectory is subsequently refined as part of the orientation reconstruction process (detailed below), in order to account for the distance between the particle's centre of mass and the point that corresponds to the centre of the particle projection.
Reconstruction of the particle's time-resolved orientation in free fall requires the use of digital models describing the geometric shape of the falling particles, such as the files used within the 3D-printing process. The geometric shape of each particle is analysed in order to determine the centre of mass, moments of inertia and principal axes; these parameters define the co-moving x ′ and co-rotational x ′′ reference frames (see Section 2.4).
The orientation of the falling particle is determined by comparing images acquired by the digital cameras against synthetic images that each camera would record if the particle had fallen at some postulated orientation, as shown in Figure 4. That is, for each time step, an initial guess of the particle orientation, described by Euler angles , , (see Section 2.4), is determined. (For the first time step, the initial guess of orientation is judged by eye, whilst for all subsequent times, the initial guess is given by the particle orientation calculated in the previous time step.) Synthetic images that each camera would record if the particle had Visualisation, in the co-moving frame x ′ , of the orientation reconstruction process showing the digital reconstruction of the falling analogue (black shading), projections equivalent to those recorded by each camera (grey shading), projections showing the locations of the particle COM x ′ COM = (0, 0, 0) ( ) and the geometric centre of the particle projection x ′ GC ( ), and the minimum bounding rectangles used to compute the geometric centre of the particle projection ("--") fallen at this assumed orientation are then reconstructed by projecting the location of the digital model of the particle onto 2D images using Equation 4, as shown in Figure 4. For each time step, the initial estimate of orientation is refined by repeating this process with a series of perturbed orientations and determining when the difference between the experimental and synthetic images is minimised. In the current apparatus, the particle orientation described by Euler angles , , is refined to within a tolerance of ±2.5 • . We note that, for analogues that exhibit radial symmetry (Section 3), any Euler angle that describes a rotation about the symmetry axis is neglected during the orientation reconstruction.
However, refinements to the preliminary estimate of the location of the particle (x GC ) must be incorporated into the orientation reconstruction in order to accurately reconstruct synthetic images using Equation 4. That is, Equation 4 requires the location of the falling analogue to be expressed in the global reference frame x. Consequently, the co-ordinates of the orientated analogue must be converted from the co-moving reference frame x ′ to the global reference frame x using an estimate of the location of the centre of mass of the analogue in the global reference frame x COM , that is x = x ′ + x COM . However, using x GC as an estimate of x COM can result in significant errors in the orientation reconstruction. As such, a correction is applied to account for the distance between the location of the geometric centre of the particle's projection (x GC ) and its centre of mass (x COM ). When expressed in terms of the co-moving frame x ′ , this distance is simply given by the location of the geometric centre of the particle's projection x ′ GC ; recall that the co-moving frame has its origin at the particle centre of mass. Here x ′ GC is estimated by assuming that the location of the centre of the particle's projection recorded by the digital cameras corresponds to the central point within the minimum bounding cuboid (with edges aligned with the axes of x ′ ) that encloses the digital model of the orientated particle, as shown in Figure 4. Combining the measurement of x GC (from images recorded by the digital cameras) and the estimate of x ′ GC (from the digital model of the particle), the location of the particle COM can be estimated as x COM = x GC − x ′ GC .

Measurements acquired using TRAIL methodology
The output of the data generated by the algorithm described in Section 2.6 facilitates the study of key aspects of the aerodynamics of ice particles. For each experiment, it is possible to rigorously analyse the reconstructed trajectory and determine whether, for example, the analogues fall steadily with a stable orientation or exhibit unsteady fluttering motions as they fall. Measurements from a series of experiments using analogues of the same geometric shape but different sizes then permits investigation into how this behaviour is influenced by the Reynolds number Re of the analogue. Note however, that measurements at a range of Reynolds numbers Re can also be acquired when using only a single particle by conducting a series of experiments in fluids of different densities and viscosities; by using both approaches in our apparatus, we are able to obtain measurements at regular intervals over Reynolds numbers spanning at least three orders of magnitude. In addition, the reconstructions directly facilitate analysis of the orientations that ice particles adopt and how this may vary with Re, which enables us to evaluate assumptions commonly used in remote-sensing applications. Moreover, fall speed parametrisations can be evaluated using measurements of Re, C D and Be that are straightforwardly determined using data extracted from the analysis.
In evaluating the key parameters Re, C D and Be, we note that u t , D and A can be evaluated directly from the digital reconstruction of the particle's trajectory. For particles that fall steadily, D is taken to represent the maximum span of the particle's projection normal to the fall direction (i.e. the maximum span in two-dimensional space), whilst A represents the projected area normal to the fall direction. When particles fall unsteadily, D and A represent the mean values of the corresponding properties across the reconstructed trajectory. We note that the definitions of D and A used here are not necessarily the same as those used in other studies. For example, when using instrumentation that provides automated measurements of fall speed from surface observations of natural ice particles, D and A are estimated from the projections of the falling particle, which are commonly viewed from the side (see for example von Lerber et al., 2017); in this case, measurements of Re and C D are not directly comparable to those recorded here.
The approach described above results in data with superior accuracy when compared with previous studies, including previous analogue studies. For example, in the current apparatus, at low Reynolds number (i.e. when particles fall steadily), the standard error of measured values of C D and Re across repeated experiments is less than 1% of their mean values. Larger standard errors (but still less than 3% of the mean values) are observed when particles fall unsteadily, which we attribute to two effects. Firstly, since the spatial extent of the region over which we reconstruct the particle trajectory is small, the average values of u t , A and D over the 20-cm sample are not fully converged when particles fall unsteadily. Secondly, results from studies conducted with simple circular discs suggest that, when particles fall at a Reynolds number close to the transition between fall regimes (i.e. periodic fluttering motion versus chaotic motion), the resulting motion observed is sensitive to small perturbations in the flow, such that different fall regimes (and terminal velocities) can be observed in repeat experiments (Auguste et al., 2013).
Of course, other sources of uncertainty also exist. Notably, the drag coefficient is linearly dependent on the density difference p − f , which is small in our case (varying between approximately 30 and 175 kg m −3 ). Consequently, small uncertainties in the particle or fluid density can propagate into large errors in C D . During experiments in the current apparatus, the density of the fluid is determined to an accuracy of ±0.5 kg m −3 using hydrometers. However, the density of the ice crystal analogues (1,174 kg m −3 ) could only be determined to an accuracy of ±2 kg m −3 by using a series of "sink/float" tests. That is, we determined the density of the analogues by determining the point at which they transitioned from sinking to floating within solutions of glycerol and water of different densities.
However, the significance of these sources of uncertainty varies between experiments. In particular, experiments conducted at low Reynolds number (i.e. Re ≲ 100) are conducted in a fluid with f ≈ 1,140 kg m −3 (i.e. large volume fraction of glycerol), such that the uncertainty associated with the density difference p − f is approximately 5%. However, in this case, the uncertainty associated with estimates of u t , A and D is small (i.e. <1%) since particles are observed to fall in a stable regime. At higher Reynolds number, the uncertainty associated with estimates of u t , A and D increases as particles are observed to fall unsteadily, but this is offset by a reduction in the uncertainty associated with the density difference p − f , since experiments are conducted in a fluid of lower density (i.e. with a lower volume fraction of glycerol). In light of these effects, we estimate that measurements of C D and Re in the current apparatus should be accurate to within approximately 5%.
Finally, we stress that the results derived through the application of the techniques described herein can provide insight into the aerodynamics of ice particles across a wide range of particle sizes. For example, in the companion manuscript, we report measurements of C D for ice particle analogues representing plates, plate polycrystals, columns, capped columns, bullet rosettes, mixed rosettes and aggregates at Reynolds numbers between Re ∼ 5 and 5000. For context, we note that previous studies investigating the fall speed of natural ice particles indicate that plates of 300 microns and 2 mm in size exhibit Reynolds numbers of order 5 and 50, respectively, aggregates of 1 and 10 mm in size exhibit Reynolds numbers of order 50 and 1, 000 respectively and graupel of 1 and 5 mm in size exhibit Reynolds numbers of order 100 and 1, 000, respectively (Kajikawa, 1972;Zikmunda, 1972;Locatelli and Hobbs, 1974;Kajikawa, 1975a;1975b). Thus, the TRAIL methodology can be used to derive insight into the aerodynamics of natural ice particles that exist under a wide range of atmospheric conditions.

VALIDATION OF EXPERIMENTAL PROCEDURES
To validate the experimental procedure described in Section 2, we conducted preliminary experiments in which measurements were acquired for circular discs of three different aspect ratios t/D = 0.04, 0.1 and 0.2, where t denotes the disc thickness, at Reynolds numbers between Re ∼ 1 and 10 3 . In these experiments, discs of six different sizes (D = 3, 5, 7, 10, 15 and 20 mm) were used in three different mixtures of water and glycerol (with 0%, 30.6% and 50.2% glycerol by volume). The discs were fabricated using the 3D-printing technique described in Section 2.2. The discs had negligible surface roughness (a layer thickness of 25 m was used during fabrication), so it was assumed that the surface roughness had no effect on the observed motion. The trajectories of the falling discs were also analysed and classified according to whether the disc exhibited (a) a steady fall with a stable orientation, (b) periodic motions or (c) chaotic motions. Examples of reconstructed trajectories classified into these categories are shown in Figure 5. (We note that error bars have not been shown in Figure 5 as they would be too small for the reader to identify them clearly within the plots.) Within the periodic regime, discs have previously been observed to exhibit both spiralling motions and fluttering motions, wherein discs oscillate along a single axis (as shown in Figure 5). The type of periodic motion observed is a complex function of the Reynolds number Re, aspect ratio t/D and density ratio p ∕ f (Zhong et al., 2011); in the current experiments, discs falling in the periodic regime principally exhibited linear fluttering motions.
However, Figure 5d illustrates that small horizontal motions can be observed even when particles fall steadily.
Such unreproducible low-frequency motions are commonly reported in experimental studies of falling discs (see for example Fernandes et al., 2007) and are thought to arise due to minor imperfections in the geometry of the falling discs or small disturbances within the ambient fluid leading to a change in the observed path of a given particle (Ern et al., 2012). In line with previous studies, these small displacements are neglected within subsequent analysis.
Measurements of Re and C D for the circular discs are shown in Figure 6, which are compared against the experimental measurements of Willmarth et al. (1964), Jayaweera (1965) and Jayaweera and Cottis (1969). The classification of the trajectories of falling discs into the categories described above is also shown in Figure 6, plotted as a function of Re and the non-dimensional moment of inertia I * ≡ p t∕64 f D. (We note that error bars have not been shown in Figure 6 as they would be too small for the reader to identify them clearly within the plots.) Also shown are the boundaries between different fall regimes that were identified in the experimental studies of circular discs reported by Willmarth et al. (1964) and Field et al. (1997). Early studies investigating the aerodynamics of falling discs suggested that their fall regime could be characterised using just Re and I * (Willmarth et al., 1964). However, more recent studies have indicated that the fall regime depends on both the density ratio p ∕ f and aspect ratio t/d and that these two effects cannot be reliably combined using I * alone (see for example Auguste et al., 2013).
For Re ≲ 100, the experimental data shown in Figure 6a indicate that C D is independent of aspect ratio and that our data are in close agreement with previous studies. Figure 6a further indicates that C D begins to increase for Re Re F I G U R E 6 (a) Measurements of Re and C D for thin circular discs. Discs of different aspect ratios are shown separately (see legend).
Also shown are experimental results from Willmarth et al., (1964), Jayaweera (1965) and Jayaweera and Cottis (1969). (b) Type of motion observed for falling thin circular discs (see legend) as a function of Re and I * . The transition points between fall regimes previously reported in the literature (Jayaweera, 1965;Field et al., 1997) are also shown Re ≳ 100 and is no longer independent of aspect ratio; larger C D values are observed for discs with a smaller aspect ratio (i.e. thinner discs). However, the data remain in good agreement with the experimental data reported by Jayaweera (1965), who used discs with similar aspect ratios to those used here (t/D = 0.007 − 0.3). Figure 6b illustrates that the observed increase in drag coefficient is coincident with the onset of unsteady motions, as observed in previous studies (see e.g. Jayaweera, 1965), which has been attributed to a change in wake structure (Auguste et al., 2013). However, Figure 6b also indicates that, in the current experiments, particles at I * ∼ 5 × 10 −3 were observed to fall in a steady regime at higher Reynolds numbers than observed in previous studies. We attribute this observation to differences in the properties of the discs used; the density ratio p ∕ f ≈ 1.03 − 1.18 in the current experiments is approximately half that of the experiments of Jayaweera (1965), in which glass discs were used, and an order of magnitude lower than in the experiments of Field et al. (1997), in which metallic discs were used. Given the previously reported sensitivity of fall regime to density ratio p ∕ f (Jayaweera, 1965;Auguste et al., 2013), this effect is likely to give rise to the difference observed between the current data and previous studies.
In light of these observations, we conclude that the experimental approach described in Section 2 gives rise to digital reconstructions of particle trajectories that yield measurements of the drag coefficient that are in close agreement (i.e. to within experimental uncertainty) with data previously reported in the literature. In addition, data arising from the digital reconstructions also enable us to classify the type of falling motions observed, which have also been found to be consistent with previously reported trends.

CASE STUDY
To further demonstrate the experimental approach described in Section 2, we also report results derived from a case study in which measurements were acquired for analogues of an aggregate ice particle generated using the model of Westbrook et al. (2004), at Reynolds numbers between Re ∼ 10 and 10 4 . In these experiments, analogues of five different sizes (volumes of 95, 185, 320, 756 and 1,477 mm 3 ) were used in three different mixtures of water and glycerol (with 0%, 30.6% and 50.2% glycerol by volume). The analogue chosen is shown in Figures 2 and 4; a full description of the geometric shape of the analogue is provided in the supplementary material. We note that experiments at Re ≈ 67, 440, and 1,800 were each repeated under nominally identical conditions a total of five times to illustrate the reliability of the measurements acquired. Error bars illustrating the range in measurements across these ensembles of experiments are shown in Figures 7e, f and 8; we expect these uncertainties to be representative of the other data shown in Figures 7 and 8 without error bars. We also stress that the relationships between dimensionless quantities reported here, measured using the 3D-printed analogues falling in a viscous fluid, are the same as we would observe for a natural ice particle of the same geometric shape falling in air, by virtue of dynamic similarity.
Examples of reconstructed trajectories of the falling analogues are shown in Figure 7a-c, which illustrate that the analogues fell in a spiralling trajectory, whereby they rotated around a vertical axis. There are numerous observations of irregular ice particles exhibiting spiralling trajectories in the literature (see e.g. Kajikawa, 1982;1989); the novel distinction here is that the digital reconstructions facilitate an analysis of whether the spiralling trajectories (e, f) Normalised radius and normalised rotation rate of spiral trajectories as functions of Reynolds number for analogues of an aggregate particle. Error bars illustrating the range in measurements recorded across repeated experiments at Re ≈ 67, 440 and 1,800 are also shown result from steady or unsteady motions. This distinction is important as different spiralling motions (i.e. steady versus unsteady orientation) manifest with different wake structures, which in turn will influence ice particle growth rates. Analysis of the reconstructed Euler angles indicates that, at low Reynolds number, the orientation with which the analogues fell is steady (with respect to the rotation axis), whilst unsteady (periodic) small-amplitude fluttering motions do occur as the Reynolds number increases. The occurrence of the small-amplitude fluttering motions is first identified at Re ∼ 10 3 . This critical Reynolds number is an order of magnitude larger than the commonly accepted value of Re ∼ 10 2 (Pruppacher and Klett, 2010, p. 445). This could, in part, be explained by noting that the Reynolds number at which the onset of unsteady motions occurs is geometry dependent (see companion paper for details). However, this can also be explained by noting that some previous studies have considered all spiralling trajectories to represent "unsteady motions", rather than distinguishing between steady spiralling trajectories, which occur when the COM and geometric centre of the particle are not aligned (Kajikawa, 1982), and unsteady spiralling trajectories, which occur as a result of wake instability at high Reynolds number. This distinction is made here since the onset of unsteady motions as a result of wake instability influences the fall speed of an ice particle (see discussion below). Further analysis of the reconstructed trajectories reveals more complex Reynolds number effects; trajectories of the COM within the horizontal x 1 -x 3 plane, normalised by the particle diameter D, for each of the previous examples are shown in Figure 7d. (Here the unreproducible low-frequency motions discussed in Section 3 have been subtracted from the trajectories to facilitate a  Mitchell (1996) and Heymsfield and Westbrook (2010). Error bars illustrating the range in measurements recorded across repeated experiments at Re ≈ 67,440 and 1,800 are also shown in each plot direct comparison between experiments.) Figure 7d illustrates that, for these examples, the normalised radius of the spiral trajectory R/D reduces as the Reynolds number increases; this reflects the wider trend found for this analogue, which is shown in Figure 7e. The reduction in the radius of the spiral trajectory coincides with an increase in the rotation rate of the analogues, as shown in Figure 7(f).
(Here the period of rotation T has been normalised by the time period D/u t it takes for each analogue to fall a distance equal to its size D.) In other words, as the Reynolds number increases, relatively speaking, the analogues fall in a narrower and more rapidly rotating spiralling trajectory. Similar results were reported by Kajikawa (1982) and Kajikawa (1989), although a weaker Reynolds number dependence was evident on account of the previous measurements' comprising ensembles of ice particles of different geometries such that the effect of Reynolds number could not be isolated. The change in spiralling trajectory described above coincides with changes in the orientation that the analogues adopt as the Reynolds number increases. That is, the analogues exhibit a clear preferred orientation as they fall, but the data reveal that this orientation does in fact vary as a function of Reynolds number. Measurements of the angle between the (vertical) x 2 axis and the largest principal axis , intermediate principal axis and smallest principal axis are shown in Figure 8a-c. (Where particles were found to fall unsteadily, the average values of , and are shown in Figure 8.) For reference, we note that these angles correspond to the angle between the x 2 axis and the blue, red and green axis, respectively, shown in Figure 2. We note that scatter shown in Figure 8a-c can largely be attributed to experimental uncertainty; recall that the orientation of the falling analogues is only resolved to within a tolerance of ±2.5 • . Regardless, Figure 8a-c reveals a clear and consistent trend relating to the orientation of the analogues, whereby at low Reynolds number analogues were observed to fall with the largest principal axis approximately aligned with the vertical fall direction (and the intermediate and smallest principal axes approximately aligned normal to the fall direction), in line with commonly used assumptions (see e.g. Ekelund and Eriksson, 2019). However, the data indicate that, as the Reynolds number increases, the orientation that the analogues adopt changes such that the largest principal axis is aligned further from the vertical fall direction whilst the intermediate principal axis is also aligned further from the horizontal plane; the smallest principal axes remains approximately aligned with the horizontal fall-normal direction. To the best of our knowledge, this is the first time that a relationship between the stable orientation that an ice particle adopts and the Reynolds number at which it falls has been identified. We attribute the change in orientation of the analogues to changes in the wake structure; experiments to validate this idea by directly measuring the flow fields around the falling particles are under way.
Finally, measurements of drag coefficients are shown in Figure 8d, which illustrate similar trends to the circular discs reported in Section 3. At low Reynolds number (when falling steadily), the drag coefficient is observed to decrease monotonically as the Reynolds number increases. However, the onset of unsteady motions at Re ∼ 10 3 results in a distinct change in the C D -Re curve, whereby further increases in Reynolds number give rise to an increase in drag coefficient. A comparison of the measured drag coefficients with the commonly used fall speed parametrisations of Mitchell (1996) and Heymsfield and Westbrook (2010) is also shown in Figure 8d. Figure 8d shows that the Mitchell (1996) parametrisation consistently under-predicts the drag coefficient, by 22% on average. In contrast, the Heymsfield and Westbrook (2010) parametrisation predicts the drag coefficient with remarkable accuracy for Re ≲ 10 3 , but fails to predict the increase in drag coefficient that occurs at the onset of unsteady motions. We note that an increase in drag coefficient at high Reynolds number can be incorporated into fall speed parametrisations through the use of a "turbulent correction" (Böhm, 1992;Mitchell, 1996;Mitchell and Heymsfield, 2005). In this case, the increase in C D at Re ∼ 10 3 shown in Figure 8d is broadly consistent with trends predicted by the "turbulent corrections", however the accuracy with which these parametrisations are able to predict the fall speed across a range of different irregular ice particles remains unclear; see the companion paper for details.

CONCLUSION
We have extended the analogue method of determining ice particle drag coefficients by developing a novel experimental approach for reconstructing the time-resolved trajectory and orientation of falling analogues. Using this system, the trajectories of 3D-printed analogues, of order 20 mm in size, are reconstructed across a sample volume of 0.2 m × 0.2 m × 0.2 m to an accuracy of better than ±0.5 mm, whilst the orientations of the falling analogues within this volume are reconstructed to an accuracy of ±2.5 • . The detailed data that can be extracted from the digital reconstructions result in measurements of the drag coefficient with superior accuracy when compared with previous studies. In addition, the experimental approach also facilitates the study of other important aspects of the aerodynamics of ice particles, such as the preferred orientation of particles with complex irregular geometries. Moreover, although the experimental approach was devised to study ice particle analogues falling in a quiescent viscous liquid, the approach could be deployed to study the aerodynamics of ice particle analogues falling in a range of gases and liquids, regardless of whether the fluid exists in a quiescent, laminar or turbulent state. The experimental approach is validated by conducting experiments on circular discs; measurements of the drag coefficient and fall regime are in close quantitative agreement with data previously reported in the literature. The detailed insight that can be derived using the experimental approach is demonstrated by reporting new experimental data from a case study in which an analogue of an aggregate ice particle is used. Analysis of the reconstructed trajectories reveals that some aspects of the aerodynamics of ice particles are poorly represented in commonly used models. For example, the analogue was found to fall in a spiralling trajectory with a steady stable orientation for Re ≲ 10 3 and only exhibited small-amplitude unsteady motions for Re ≳ 10 3 . These observations are inconsistent with the commonly accepted notion that unsteady motions occur at Re ∼ 10 2 (Pruppacher and Klett, 2010, p. 445). In addition, whilst the analogue was observed to exhibit a clear preferred orientation during free fall, the orientation adopted is shown to depend crucially on the Reynolds number at which it falls. To the best of our knowledge, this is the first time that a relationship between the stable orientation that an ice particle adopts and the Reynolds number at which it falls has been identified. Furthermore, a comparison of measurements of the drag coefficient with commonly used fall speed parametrisations reveals that the Mitchell (1996) parametrisation consistently under-estimates the drag coefficient, by 22% on average, whereas the Heymsfield and Westbrook (2010) parametrisation predicts the drag coefficient with remarkable accuracy for Re ≲ 10 3 but fails to predict the increase in drag coefficient that occurs at the onset of unsteady motions (i.e. for Re ≳ 10 3 ). The ideas briefly presented in this case study are to be explored in much greater detail, using a much more comprehensive dataset, in the companion manuscript and subsequent publications.