M Z Sheikh

K Gustavsson

E Lévêque

B Mehlig

A Pumir

A Naso

Colliding ice crystals in turbulent clouds

Collisions, resulting in aggregation of ice crystals in clouds, is an important step in the formation of snow aggregates.

Here, we study the collision process by simulating spheroid-shaped particles settling in turbulent flows, and by determining the probability of collision. We focus on plate-like ice crystals (oblate ellipsoids), subject to gravity, to the Stokes force and torque generated by the surrounding fluid. We also take into account the contributions to the drag and torque due to fluid inertia, which are essential to understand the tendency of crystals to settle with their largest dimension oriented horizontally. We determine the collision rate between identical crystals, of diameter 300𝜇𝑚, with aspect ratios in the range 0.005 ≤ 𝛽 ≤ 0.05, and over a range of energy dissipation per unit mass, 𝜀, 1 cm 2 /s 3 ≤ 𝜀 ≤ 250 cm 2 /s 3 . For all values of 𝛽 studied, the collision rate increases with the turbulence intensity. The dependence on 𝛽 is more subtle. Increasing 𝛽 at low turbulence intensity (𝜀 ≲ 16 cm 2 /s 3 ) diminishes the collision rate, but increases it at higher 𝜀 ≈ 250 cm 2 /s 3 . The observed behaviors can be understood as resulting from three main physical effects. First, the velocity gradients in a turbulent flow tend to bring particles together. In addition, differential settling plays a role at small 𝜀 when the particles are thin enough (𝛽 small), whereas the prevalence of particle inertia at higher 𝜀 leads to a strong enhancement of the collision rate.

Introduction

As ice crystals collide in a turbulent cloud, they may stick to each other and coalesce to form larger aggregates [START_REF] Pruppacher | Microphysics of clouds and precipitation[END_REF]. This process is an important ingredient in the growth of ice crystals and the formation of snow aggregates [START_REF] Lo | The growth of snow in winter storms: an airborne observational study[END_REF][START_REF] Mitchell | Evolution of snow-size spectra in cyclonic storm. part i: Snow growth by vapor deposition and aggregation[END_REF]. Detailed field measurements have shown that the resulting aggregates have a characteristic exponential size distribution [START_REF] Field | Aggregation and scaling of ice crystal size distributions[END_REF], and are shaped approximately as prolate ellipsoids [START_REF] Jiang | Shapes and fall orientations of ice particle aggregates[END_REF]. The determination of the probability of collision of ice-crystals appears as a crucial element in the quantitative description of this process [START_REF] Khain | Physical processes in clouds and cloud modeling[END_REF][START_REF] Gavze | Gravitational collision of small non-spherical particles: swept volumes of prolate and oblate spheroids in calm air[END_REF].

Among the collision processes involving hydrometeors, those between ice crytals appear to be the least understood [START_REF] Wang | Physics and Dynamics of Clouds and Precipitation[END_REF]. Collisions between small droplets have been extensively studied theoretically and numerically over the past decades [START_REF] Shaw | Particle-turbulence interactions in atmospheric clouds[END_REF][START_REF] Franklin | Collision rates of cloud droplets in turbulent flow[END_REF][START_REF] Falkovich | Sling effect in collisions of water droplets in turbulent clouds[END_REF]Ayala et al. 2008b,a;[START_REF] Grabowski | Growth of cloud droplets in a turbulent environment[END_REF][START_REF] Voßkuhle | Prevalence of the sling effect for enhancing collision rates in turbulent suspensions[END_REF][START_REF] Saito | Turbulence and cloud droplets in cumulus clouds[END_REF]. In the case of ice crystals, the complexity due to the particle anisotropy and to the possibility of crystals to settle with different orientations leads to fundamental difficulties. The settling of ice crystals was studied experimentally from field observations [START_REF] Nakaya | Simultaneous observations of the mass, falling velocity and form of individual snow crystals[END_REF][START_REF] Matrosov | Inferring fall attitudes of pristine dendritic crystals from polarimetric radar data[END_REF] or tank measure-Corresponding author: Alain Pumir, alain.pumir@ens-lyon.fr

Corresponding author: Aurore Naso, aurore.naso@ec-lyon.fr ments of idealized objects settling in a fluid [START_REF] Willmarth | Steady and unsteady motions and waks of freely falling disks[END_REF]. The first numerical investigations of settling ice crystals assumed a fixed horizontal orientation (major axis perpendicular to the gravity direction) [START_REF] Pitter | A numerical study of viscous flow past a thin oblate spheroid at low and intermediate reynolds numbers[END_REF][START_REF] Wang | Numerical simulation of three-dimensional unsteady flow past ice crystals[END_REF]. [START_REF] Hashino | Numerical simulation of the flow fields around falling ice crystals with inclined orientation and the hydrodynamic torque[END_REF] imposed a steady flow around columns and plates with a fixed (but non-horizontal) orientation. More recently, it became possible to study numerically freely falling hexagonal ice plates [START_REF] Cheng | A numerical study on the attitudes and aerodynamics of freely falling hexagonal ice plates[END_REF] and planar crystals of more complex shapes [START_REF] Nettesheim | A numerical study on the aerodynamics of freely falling planar ice crystals[END_REF]. All these studies were carried out in laminar flows (fixed uniform flow around a fixed particle, or particle settling in a fluid at rest). It also became possible to simulate numerically the settling of ice crystals in turbulence (Siewert et al. 2014a,b;[START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF][START_REF] Naso | Collision rate of ice crystals with water droplets in turbulent flows[END_REF][START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF], and also to investigate numerically the collision of crystals in a turbulent flow (Siewert et al. 2014b;[START_REF] Gustavsson | Statistical model for the orientation of nonspherical particles settling in turbulence[END_REF][START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF][START_REF] Gavze | Gravitational collision of small non-spherical particles: swept volumes of prolate and oblate spheroids in calm air[END_REF]. Experimentally, crystal-crystal collisions were mainly studied by considering small ice crystals colliding with a fixed ice or artificial target [START_REF] Hosler | The aggregation of small ice crystals[END_REF][START_REF] Keith | The collection efficiency of a cylindrical target for ice crystals[END_REF]. More recent experiments focus on the aggregation efficiency [START_REF] Connolly | A laboratory investigation into the aggregation efficiency of small ice crystals[END_REF].

The goal of the present work is to provide an estimate of the collision kernel between settling ice crystals in turbulent flows. Notoriously, the turbulence intensity in clouds, as measured by the rate of dissipation of kinetic energy per unit mass, 𝜀, varies over a wide range: from 𝜀 ∼ 1 cm 2 /s 3 in weakly turbulent clouds up to ≳ 10 3 cm 2 /s 3 in cumulonumbi [START_REF] Grabowski | Comments on "preferential concentration of cloud droplets by turbulence: Effects on the early evolution of cumulus droplet spectra[END_REF].

2

We focus here on plate-like ice crystals. The convenient simplification, used here as well as in earlier works [START_REF] Pitter | A numerical study of viscous flow past a thin oblate spheroid at low and intermediate reynolds numbers[END_REF], consists in approximating the crystal as a thin oblate spheroid. The density of the crystals, 𝜌 𝑝 is much larger than that of the surrounding fluid, 𝜌 𝑓 . We took here the same values as in [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF], namely 𝜌 𝑝 = 0.9194 g/cm 3 and 𝜌 𝑓 = 1.413 × 10 -3 g/cm 3 . The description of anisotropic crystals requires an additional rotational degree of freedom, which affects the way the crystals sample the fluid. In particular, the settling velocity depends on the orientation of the crystal. This provides a mechanism that potentially enhances collision rate between crystals (Siewert et al. 2014b;[START_REF] Gustavsson | Statistical model for the orientation of nonspherical particles settling in turbulence[END_REF][START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF]. Models that correctly capture the orientation dynamics of the settling crystals are therefore essential to reliably describe the collision processes. Determining the orientation dynamics and the settling velocity of ice crystals is also a crucial ingredient for the development of accurate numerical methods of cloud microphysics simulation [START_REF] Shima | Predicting the morphology of ice particles in deep convection using the super-droplet method: development and evaluation of scale-sdm 0.2.5-2.2.0, -2.2.1, and -2.2.2[END_REF]). Here, we consider very small crystals, with a semi-major axis 𝑎 = 150𝜇m. This size is smaller than the size of the smallest eddies in the flows considered. The modeling of the force and torque acting on such crystals is crucial for the present study. Recent experimental work [START_REF] Lopez | Inertial effects on fibers settling in a vortical flow[END_REF] has shown the importance of fluid inertia in determining the torque acting on crystals. This aspect can be understood by estimating the various contributions to the torque acting on a crystal [START_REF] Sheikh | Importance of fluid inertia for the orientation of spheroids settling in a turbulent flow[END_REF], which leads to the conclusion that under a broad range of conditions, particles settle preferentially with their largest dimension oriented perpendicular to their direction of fall [START_REF] Kramel | Non-spherical particle dynamics in turbulence[END_REF][START_REF] Menon | Theoretical predictions of the orientation distribution of high-aspect-ratio, inertial particles settling in isotropic turbulence[END_REF][START_REF] Gustavsson | Theory for the effect of fluid inertia on the orientation of a small spheroid settling in turbulence[END_REF]. [START_REF] Gustavsson | Theory for the effect of fluid inertia on the orientation of a small spheroid settling in turbulence[END_REF] undertook a complete study to describe the statistics of crystals orientation in turbulent flow by taking into account fluid inertia in the expressions of the force and torque acting on a spheroid [START_REF] Cox | The steady motion of a particle of arbitrary shape at small Reynolds numbers[END_REF][START_REF] Khayat | Inertia effects on the motion of long slender bodies[END_REF][START_REF] Dabade | Effects of inertia and viscoelasticity on sedimenting anisotropic particles[END_REF]. This led to a full description of the statistics of orientations in the overdamped regime, i.e., when particle inertia can be neglected [START_REF] Gustavsson | Theory for the effect of fluid inertia on the orientation of a small spheroid settling in turbulence[END_REF], also derived by [START_REF] Anand | Orientation dynamics of sedimenting anisotropic particles in turbulence[END_REF]. A description of the statistics of orientation appropriate for a wide range of cloud conditions, consistent with measured orientation distributions, has been more recently derived [START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF], extending the earlier work of [START_REF] Klett | Orientation model for particles in turbulence[END_REF]. In this work, we perform direct numerical simulations of homogeneous isotropic turbulent flows, and follow the motion of spheroids, using the set of equations described in [START_REF] Dabade | Effects of inertia and viscoelasticity on sedimenting anisotropic particles[END_REF]; [START_REF] Sheikh | Importance of fluid inertia for the orientation of spheroids settling in a turbulent flow[END_REF], and recently used in a cloud microphysics context by [START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF]. We vary the energy dissipation range over the realistic range 1 cm 2 /s 3 ≤ 𝜀 ≤ 250 cm 2 /s 3 , as well as the aspect ratio of the crystals, keeping the diameter fixed and equal to 300 𝜇m, and we determine the collision kernels for these cases. We show that collisions may result from several possible effects. One mechanism, originally proposed by [START_REF] Saffman | On the collision of drops in turbulent clouds[END_REF], is due to the fact that the erratic turbulent motions may bring particles together. The possibility of a broad distribution of orientation favors collisions due to the different settling velocities of particles. Finally, at high turbulence intensity, the particle inertia becomes important, and makes it possible for particles not following the flow to collide. This effect, denoted as the 'sling effect' [START_REF] Falkovich | Particles and fields in fluid turbulence[END_REF][START_REF] Falkovich | Sling effect in collisions of water droplets in turbulent clouds[END_REF] or alternatively as a manifestation of 'caustics' [START_REF] Wilkinson | Caustics in turbulent aerosols[END_REF][START_REF] Wilkinson | Caustic activation of rain showers[END_REF], enhances the probability of collisions in the suspension.

The problem studied here is complex, with a large number of mechanisms and non-dimensional parameters at play. In this paper we focus on the effect of particle and fluid inertia on ice-crystal collisions, using a highly idealised model that does not take into account interactions, which may become important when the crystals are close: hydrodynamic interactions [START_REF] Shaqfeh | Orientational dispersion of fibers in extensional flows[END_REF], Coulomb interactions [START_REF] Magnusson | Collisions of micron-sized, charged water droplets in still air[END_REF], and the breakdown of hydrodynamics (the fact that the discrete molecular nature of the gas becomes important) when the gap between the particles is of the order of the mean free path of the ambient air [START_REF] Sundararajakumar | Non-continuum lubrication flows between particles colliding in a gas[END_REF]. The model describes how the crystals approach, but cannot predict collision outcomes. When counting collisions we use a simplified approach, the ghost-collision scheme [START_REF] Wang | On the collision rate of small particles in isotropic turbulence. i. zero inertia case[END_REF][START_REF] Gustavsson | Collisions of particles advected in random flows[END_REF][START_REF] Gustavsson | Statistical model for collisions and recollisions of inertial particles in mixing flows[END_REF], which simply counts how frequently crystals approach each other, disregarding any change in boundary conditions e.g. due to fragmentation or coalescence.

Our work is organized as follows. The simplifying assumptions made to treat the problem, as well as the numerical methodology are presented in section. 2. In section 3 we show results of our numerical simulations for the orientation, settling velocity and collisions of the crystals. The collision statistics are then discussed in section 4. Section 5 finally contains our conclusions.

Problem setup

The way we set up the problem here closely follows [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF]; [START_REF] Naso | Collision rate of ice crystals with water droplets in turbulent flows[END_REF].

Turbulent flow simulations.

Although turbulence in clouds is affected by various parameters, particularly by the flow stratification, we simplify the problem by restricting ourselves to homogeneous isotropic turbulence. Our work is based on direct numerical simulations of the Navier-Stokes equations:

𝜕 𝑡 u + (u • ∇)u = - ∇𝑝 𝜌 𝑓 + 𝜈∇ 2 u + f, (1)
∇ • u = 0, ( 2
)
where u is the velocity of the fluid, 𝑝 its pressure, 𝜌 𝑓 and 𝜈 are the mass density and kinematic viscosity, respectively.

We take the values of 𝜌 𝑓 = 1.413 × 10 -3 g/cm 3 , and 𝜈 = 0.1132 cm 2 /s. We use a triply periodic box, of size 𝐿 = 8× 𝜋cm ≈ 25 cm in each spatial direction. The number of grid points, 𝑁, or equivalently of Fourier modes in each direction is varied to maintain sufficient accuracy to allow us to reliably interpolate the velocity field and its derivatives at the particle positions, as necessary to solve the equations of motion for the particles. We insisted in our simulations that the largest wavenumber resolved, 𝑘 𝑚𝑎𝑥 , satisfies 𝜂𝑘 𝑚𝑎𝑥 ≳ 3, where 𝜂 = (𝜈 3 /𝜀) 1/4 is the size of the smallest eddies in the flow (𝜂 is known as the Kolmogorov scale). Further details about the code can be found in [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF]; [START_REF] Naso | Collision rate of ice crystals with water droplets in turbulent flows[END_REF]. We picked three values of 𝜀, 𝜀 ≈ 1 cm 2 /s 3 (flow I, runs 1-3), 𝜀 ≈ 16 cm 2 /s 3 (flow II, runs 4-7) and 𝜀 ≈ 250 cm 2 /s 3 (flow III, runs 8-10), further details are given in Table 1.

Ice crystals: elementary simplifications. We focus here on small planar ice-crystals, whose shape we approximate as oblate spheroids. More precisely, we consider an ellipsoid of revolution, with a semi-major axis 𝑎 and a semi-minor axis 𝑐, the aspect ratio 𝛽 = 𝑐/𝑎 being taken very small, 𝛽 = 0.005, 0.01, 0.02 and 0.05. We added the smallest value, 𝛽 = 0.005, compared to the set of values considered in [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF], to better illustrate the interplay between different physical effects, as discussed in Section 4.

As the volume fraction of particles in clouds is very small, it is justified to ignore the feedback of the particles on the flow, and to use the one-way coupling approximation [START_REF] Elgobashi | On predicting particle-laden turbulent flows[END_REF].

The motion of the crystals is obtained by solving Newton's equations in the inertial reference frame for the center of mass of the particle, x, and for its velocity, v, as well as for the orientation, n, defined as the unit vector parallel Table 1. Properties of the three turbulent flows used in this study. We simulated homogeneous isotropic flows in triply periodic boxes of size 𝐿 = 8 𝜋 cm ≈ 25 cm, with a fixed viscosity 𝜈 = 0.1132 cm 2 /s, injecting a fixed amount of energy 𝜀 in the system. We took three different values of 𝜀, resulting in three different flows, as listed in the to the axis of symmetry of the spheroid, as illustrated in Fig. 1, and for its angular velocity, ω:

𝑑x 𝑑𝑡 = v , 𝑚 𝑑v 𝑑𝑡 = f ℎ + 𝑚g, (3)
where f ℎ is the hydrodynamic force acting on the particle, and

𝑑 n 𝑑𝑡 = ω × n , 𝑚 𝑑 𝑑𝑡 I( n)ω = τ ℎ , (4)
where τ ℎ is the hydrodynamic torque acting on it. In Eqs.

(3) and (4), 𝑚 denotes the mass of the particles, and I( n) the moment of inertia tensor, whose expression is given in Appendix A. The modeling challenge consists in providing the correct expressions for the force f ℎ and for the torque τ ℎ exerted by the fluid on the particles.

Our approach rests on the assumption that the particle Reynolds number 𝑅𝑒 𝑝 = 𝑎|u-v|/𝜈, based on the difference between the flow velocity u and the particle velocity v, is small. Under this condition, it is often assumed that the hydrodynamic force and torque acting on the ice crystals can be derived in the Stokes limit, i.e. when all the nonlinear terms in the Navier-Stokes equations can be neglected [START_REF] Kim | Microhydrodynamics: principles and selected applications[END_REF][START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF]. In this approach, the effect of fluid inertia is therefore completely neglected (Siewert et al. 2014b;[START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF][START_REF] Gavze | Gravitational collision of small non-spherical particles: swept volumes of prolate and oblate spheroids in calm air[END_REF]. When the carrying flow is turbulent, such an assumption leads to the qualitatively incorrect conclusion that particles tend to settle with their major axis oriented parallel to gravity (Siewert et al. 2014a;[START_REF] Gustavsson | Statistical model for the orientation of nonspherical particles settling in turbulence[END_REF].

In fact, the Reynolds number happens to be of order unity in the realistic parameter regime considered here, see also [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF]. For these reasons, the first correction induced by fluid inertia must be accounted for in the expressions of f ℎ and τ ℎ , in a systematic expansion in powers of 𝑅𝑒 𝑝 . Whereas the correction to the force has been known for some time [START_REF] Brenner | The Oseen resistance of a particle of arbitrary shape[END_REF][START_REF] Khayat | Inertia effects on the motion of long slender bodies[END_REF], the expression for the torque in the small 𝑅𝑒 𝑝 limit has been derived more recently [START_REF] Dabade | Effects of inertia and viscoelasticity on sedimenting anisotropic particles[END_REF]. The explicit forms of the force f ℎ and of the torque τ ℎ including both the Stokes contribution and the first effect of fluid inertia are provided in Appendix A. They have been validated numerically by directly solving the Navier-Stokes equations around a spheroid [START_REF] Jiang | Inertial torque on a small spheroid in a stationary uniform flow[END_REF], and confronted to experimental measurements of cylindrical fibers settling in a viscous fluid [START_REF] Cabrera | Experimental validation of fluid inertia models for a cylinder settling in a quiescent flow[END_REF]. A systematic analysis of the resulting equations [START_REF] Sheikh | Importance of fluid inertia for the orientation of spheroids settling in a turbulent flow[END_REF][START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF]) (see also [START_REF] Kramel | Non-spherical particle dynamics in turbulence[END_REF]; [START_REF] Anand | Orientation dynamics of sedimenting anisotropic particles in turbulence[END_REF] shows that when the settling velocity 𝑈 𝑠 is large compared to the velocity of the smallest eddies in the flows, 𝑢 𝐾 = (𝜈𝜀) 1/4 , then the torque induced by fluid inertia prevails over the Stokes torque, which can lead to a preferential orientation of the spheroids with their largest dimension perpendicular to gravity. We stress that the novel aspect of the present work compared to the study of [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF] comes from the use of force and torque expressions including the effect of fluid inertia. This modification plays a very significant role in predicting the distribution of orientations of the crystals, as documented by [START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF] (recall that considering the Stokes torque only led to a preferential orientation of the spheroids with their major axis aligned with gravity). We will show in Sec. 3 and discuss in Sec. 5 the effect of this correction on the crystal-crystal collisions.

The equations of motion of the particles, Eqs. (3,4), involve several important dimensional quantities [START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF]). The first one is 𝜏 𝑝 , the characteristic time of the particle originating from the Stokes drag, and the second one is 𝑤 𝑠 , the typical settling velocity, which we take, following [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF], to be respectively equal to:

𝜏 𝑝 = 2 9 𝜌 𝑝 𝜌 𝑓 𝛽𝑎 2 𝜈 and 𝑤 𝑠 = 𝑔𝜏 𝑝 . (5)
It is convenient to make the time 𝜏 𝑝 and the velocity 𝑤 𝑠 dimensionless by dividing them respectively by 𝜏 𝐾 and by the velocity across the smallest eddies in the flow, 𝑢 𝐾 ≡ (𝜈𝜀) 1/4 . This leads to the definitions of the Stokes number St and of the settling number Sv, which characterizes the settling velocity [START_REF] Devenish | Droplet growth in warm turbulent clouds[END_REF][START_REF] Wang | Settling velocity and concentration distribution of heavy particles in homogeneous isotropic turbulence[END_REF][START_REF] Grabowski | Comments on "preferential concentration of cloud droplets by turbulence: Effects on the early evolution of cumulus droplet spectra[END_REF][START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF], as:

St = 4𝛽 9 𝜌 𝑝 𝜌 𝑓 𝑎 𝜂 2 and Sv = 4𝛽 9 𝜌 𝑝 𝜌 𝑓 𝑎 2 𝑔 (𝜈 5 𝜀) 1/4 , (6)
where we recall that 𝜂 = (𝜈 3 /𝜀) 1/4 is the Kolmogorov scale. Note that the quantities Sv and St, as defined by Eq. ( 6), differ by a factor 2 from the ones used in [START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF].

The approach used to solve the equations of motion for the crystals, Eq. (3,4), follows the method introduced by Siewert et al. (2014a) and by [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF]. Equation (3) has been integrated in time by using a secondorder Heun scheme (two-stage Runge-Kutta method). The fluid velocity and velocity gradients at the particle location are estimated from grid values by a tricubic interpolation. Equation (4) for the angular velocity is integrated in time in the particle-fixed coordinate system (center of mass and principal axes) where the moment of inertia tensor is constant and diagonal. The spatial rotation between the inertial and the particle-fixed frame is conveniently handled by using quaternions and tracking the orientation of each particle. Let us note that an additional inertial term arises in the (transformed) equation due to the non-Galilean nature of the particle-fixed frame. Equation ( 4) was also integrated in time with a second-order Heun scheme. Consistent with the periodic boundary conditions used, we reinjected particles leaving the box through one of the sides to the opposite side. In all the runs considered in this work, the characteristic time taken by a particle to settle through the box is larger than ≳ 2.5𝑇 𝐿 , where 𝑇 𝐿 is the large-eddy turnover time. Therefore, the flow has time to decorrelate between two passages of the particles.

Collision detection. Numerically, we determine the geometric collision rate [START_REF] Wang | Theoretical formulation of collision rate and collision efficiency of hydrodynamically interacting cloud droplets in turbulent atmosphere[END_REF]Ayala et al. 2008b;[START_REF] Grabowski | Growth of cloud droplets in a turbulent environment[END_REF]. This approximation neglects the hydrodynamic interactions, as well as any interaction between particles as they approach each other [START_REF] Pruppacher | Microphysics of clouds and precipitation[END_REF][START_REF] Batchelor | Hydrodynamic interaction of two small freely-moving spheres in a linear flow field[END_REF][START_REF] Chun | Clustering of aerosol particles in isotropic turbulence[END_REF]. Furthermore, we simply detect, from a large set of numerical trajectories of ellipsoids, the instants when two particles get into contact ("collide") with each other, while keeping the two colliding particles thereafter without any modification. The errors induced by this approximation, known as the "ghost-collision" approximation [START_REF] Wang | On the collision rate of small particles in isotropic turbulence. i. zero inertia case[END_REF][START_REF] Gustavsson | Collisions of particles advected in random flows[END_REF][START_REF] Gustavsson | Statistical model for collisions and recollisions of inertial particles in mixing flows[END_REF], do not exceed ∼ 15% in the case of spherical particles in a turbulent flow [START_REF] Voßkuhle | Multiple collisions in turbulent flows[END_REF].

To determine when a collision occurs, we implemented the algorithm developed by [START_REF] Choi | Continuous collision for ellipsoids[END_REF] to detect when two ellipsoids, moving in three dimensions, touch each other. The algorithm has been used successfully to detect collisions of settling ellipsoids in a turbulent flow (Siewert et al. 2014b;[START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF].

In practice, we simulated a million trajectories for each type of particle, for a time of ≥ 50𝑇 𝐿 for flow I, ≥ 30𝑇 𝐿 for flow II and ≥ 20𝑇 𝐿 for flow III. Note that we waited before recording any property for at least 50𝜏 𝜂 , so the particles could reach a statistically steady state. The data shown in the following sections correspond to ≳ 25, 000 collisions. We also checked the statistical significance of the numbers discussed below by dividing the total run into four or five identical sub-runs, and determining the number of collisions for each of these sub-runs. These comparisons 10 0 10 2 10 0 10 2 0.005 0.01 0.02 Fig. 2. Statistics of orientation: fluctuations of the angle of the settling spheroids with respect to the equilibrium position (𝜑 = 0 in still fluid). For each point, the statistical uncertainty has been evaluated by calculating the deviations between four subsets of the total data. The resulting deviations are ≲ 5% in most cases; they are ∼ 7% and ∼ 10% for the ( 𝜀 ≈ 1, 𝛽 = 0.005) and ( 𝜀 ≈ 250, 𝛽 = 0.02) cases, respectively. convinced us that the values of the collision kernel 𝐾 shown below are accurate to within a few percents.

Results

In this Section, we discuss the results of our model simulations. A summary of the runs, on which the following discussion is based, can be found in Table B1, see Appendix B.

a. Orientation and settling velocity of the crystals

The orientation of the settling crystals is a crucial aspect for the physical effects leading to collision between crystals, in addition to its importance to explain the reflection of light by crystals in clouds [START_REF] Klett | Orientation model for particles in turbulence[END_REF][START_REF] Bréon | Horizontally oriented plates in clouds[END_REF][START_REF] Yang | On the radiative properties of ice clouds: Light scattering, remote sensing and radiation parametrization[END_REF]. We recall that in the absence of any flow, oblate spheroids settle with their major axis oriented perpendicular to gravity. In the case of plate-like ice crystals, the unit vector parallel to the axis of symmetry of particle, n, is therefore parallel to ĝ, the vertical direction, parallel to the gravity g ≡ |g| ĝ. We pick the vector n to point downward, so that n • ĝ ≥ 0, and we denote the angle between n and ĝ as 𝜑, as illustrated in Fig. 1. Numerically, we recorded the fluctuations of 𝜑 of the simulated particles settling through the flow.

Figure 2 shows the rms of the angle fluctuations, ⟨𝜑 2 ⟩ 1/2 , for all our runs except run 7 (flow II and 𝛽 = 0.05). This figure shows the dependence of ⟨𝜑 2 ⟩ 1/2 as a function of the turbulence intensity 𝜀, for the three values of 𝛽: 𝛽 = 0.005 (rightwards pointing triangle), 𝛽 = 0.01 (leftwards pointing triangle) and 𝛽 = 0.02 (upwards pointing triangle). The behavior shown in Fig. 2 follows the predictions of our earlier 10 0 10 2 1 1.2 1.4

1.6 0.005 0.01 0.02 Fig. 3. Dependence of the averaged settling velocity of crystals as a function of 𝛽 and 𝜀. The three curves show the ratio between the averaged settling velocity 𝑈 𝑠 and the settling velocity of the same crystals in still fluid, 𝑈 0 , as a function of the energy dissipation. Each curve corresponds to a value of 𝛽 (𝛽 = 0.005: rightwards pointing triangles; 𝛽 = 0.01: leftwards pointing triangles, and 𝛽 = 0.02: upwards pointing triangles). The settling velocity increases with the intensity of turbulence. The statistical uncertainty has been evaluated using the same method as in Fig. 2. The deviations are always ≲ 3%, except for the ( 𝜀 ≈ 1, 𝛽 = 0.005) case, for which it is ∼ 7%.

studies [START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF]. Specifically, we find that the orientation distribution of the crystals is very narrowly peaked around 𝜑 = 0 at low values of 𝜀, and becomes broader as the turbulence intensity increases. Figure 2 also shows that the distribution of particle orientation is more peaked when 𝛽 is large, and becomes broader when 𝛽 decreases. In quantitative terms, we observe that the fluctuation of the particle orientation does not exceed 4 • for 𝜀 ≈ 1 cm 2 /s 3 and 𝛽 ≥ 0.01, and for 𝜀 ≈ 16 cm 2 /s 3 and 𝛽 ≥ 0.02, see Table B1.

Settling velocity. We begin by stressing that, as shown by [START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF], the model used here predicts settling velocities in still air in good agreement with those observed experimentally [START_REF] Pruppacher | Microphysics of clouds and precipitation[END_REF]. The first effect of the orientation fluctuations is to affect the settling velocity of the crystals. This is illustrated in Fig. 3, which shows the averaged settling velocity, 𝑈 𝑠 = -𝑣 𝑧 , where 𝑣 𝑧 is the vertical component of velocity (see Fig. 1). In Fig. 3, 𝑈 𝑠 is normalized by the velocity of particles settling in still fluid, 𝑈 0 . The ratios are very close to 1 at low values of 𝜀, and gradually increase with 𝜀. A first explanation for this is provided by the finding summarized above, that the orientation distribution broadens due to the increase of the turbulence intensity. In fact, the largest drag in the direction of gravity is obtained when 𝜑 = 0; so increasing 𝜑 leads to a smaller drag, and therefore to a larger settling velocity. A larger rms of the orientation fluctuation therefore increases the settling velocity 𝑈 𝑠 , compared to the reference value 𝑈 0 in the absence of any fluid motion. The effect of turbulence on the settling of particles has been widely investigated for spherical objects [START_REF] Gavze | Gravitational collision of small non-spherical particles: swept volumes of prolate and oblate spheroids in calm air[END_REF][START_REF] Ireland | The effect of Reynolds number on inertial particle dynamics in isotropic turbulence. Part 2. Simulations with gravitational effects[END_REF][START_REF] Rosa | Settling velocity of small inertial particles in homogeneous isotropic turbulence from high-resolution DNS[END_REF]): turbulence most often increases their settling velocity, as observed herein for spheroidal particles, and more rarely decreases it. Several physical mechanisms were proposed to explain these results: the preferential sweeping [START_REF] Maxey | The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields[END_REF], according to which heavy particles should sample preferentially downflow regions, is expected to increase the settling, whereas the loitering effect (falling particles spend more time in the regions with upward flow) [START_REF] Nielsen | Turbulence effects on the settling of suspended particles[END_REF] and vortex trapping should decrease the settling velocity.

The increase in the settling velocity when the orientation of the crystal differs from 𝜑 = 0 is further illustrated in Fig. 4, which shows the average of the settling velocity of crystals conditioned on the cosine of the angle 𝜑, or equivalently, on the 𝑧-component of the vector n aligned with the axis of symmetry of the spheroid, see Fig. 1. The orientation 𝜑 = 0 corresponds here to 𝑛 𝑧 = 1. The conditional averages shown in Fig. 4 correspond to the thinnest crystals, with the value of 𝛽 = 0.005, at the different values of 𝜀. For the three turbulence intensities, the numerical results reveal a clear increase of the settling velocity when 𝑛 𝑧 decreases towards 0. For the lowest turbulence intensity, the settling velocity at 𝑛 𝑧 ≈ 0, corresponding to a particle oriented with its major axis parallel to gravity, is roughly equal to 3/2 times the velocity of the particle settling with 𝑛 𝑧 = 1 (major axis perpendicular to gravity), which corresponds to the ratio between the drag of the particles settling with the two orientations in a quiescent fluid. As seen in a previous study, where the effect of fluid inertia on the dynamics of crystals was (incorrectly) neglected (Siewert et al. 2014b;[START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF], the systematic difference of the settling velocity on the crystal orientation plays an important role and tends to enhance collisions. Interestingly, we notice that the value of the settling velocity conditioned on 𝑛 𝑧 ≈ 1 is lower than the expected value in a quiescent fluid for the highest turbulence intensity. This result suggests that ice crystals with their major axis oriented perpendicular to gravity may be more sensitive to the loitering effect and vortex trapping, and/or less sensitive to preferential sweeping.

b. The collision kernel

In a dilute solution of volume 𝑉, containing 𝑁 particles, the average number of collisions, 𝑁 𝑐 , occurring over a duration 𝑇, is proportional to 𝑇, to the inverse of the volume, 𝑉 -1 , and to 𝑁 2 , and can therefore be written as [START_REF] Voßkuhle | Multiple collisions in turbulent flows[END_REF][START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF]: where 𝐾 is the collision kernel. The dimension of 𝐾 is a volume divided by a time.

𝑁 𝑐 = 1 2 𝐾 × 𝑁 2 𝑉 × 𝑇 (
Figure 5 shows the collision kernels 𝐾 determined numerically for our model system. Figure 5 𝛽 at 𝜀 ≈ 1 cm 2 /s 3 (+ symbols), 𝜀 ≈ 16 cm 2 /s 3 (diamond symbols) and 𝜀 ≈ 250 cm 2 /s 3 (square symbols). The statistical uncertainty has been evaluated using the same method as in Fig. 2. The deviations are always between 2 and 7%, except for the ( 𝜀 ≈ 250, 𝛽 = 0.02) case, for which it is ∼ 9%. creasing the turbulence intensity enhances the probability of collision. As a reference, we indicate as a dashed line the 𝜀 1/2 dependence, expected when collisions are due to the velocity gradients bringing particles together [START_REF] Saffman | On the collision of drops in turbulent clouds[END_REF][START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF]. Whereas the dependence of 𝐾 on 𝜀 is qualitatively close to this expectation when 𝜀 varies from ≈ 1 cm 2 /s 3 (flow I) to ≈ 16 cm 2 /s 3 (flow II), the evolution at higher values of 𝜀 shows a strong deviation from this behavior. The collision rate 𝐾 for the thickest particles, 𝛽 = 0.02, which was smaller for 𝜀 ≤ 16 cm 2 /s 3 , grows far more rapidly than 𝜀 1/2 . On the other hand, for the thinnest crystals, with 𝛽 = 0.005, the value of 𝐾, which was the largest at 𝜀 ≤ 16 cm 2 /s 3 , grows slower than 𝜀 1/2 between 𝜀 ≈ 16 cm 2 /s 3 and 𝜀 ≈ 256 cm 2 /s 3 . Remarkably, whereas the collision rate decreases at increasing 𝛽 for 𝜀 ≲ 16 cm 2 /s 3 , it increases at the largest value of 𝜀 simulated, as shown on Fig. 5(b). This complex dependence on the parameters of the problem can only be interpreted as the result of a competition between several different physical mechanisms in the dynamics of collisions. The following section is devoted to a closer analysis of these mechanisms.

Discussion

a. Collision mechanisms

The original Saffman-Turner theory [START_REF] Saffman | On the collision of drops in turbulent clouds[END_REF], valid for particles following exactly the flow, provides a convenient starting point to discuss the collision mechanisms at play in the present problem. We recall that the physics is based on expressing the collision kernel 𝐾 in a turbulent suspension of spheres, of radius 𝑎, as a flux of particles through a sphere centered at the position of a given particle in the suspension, of radius 2𝑎, the particles being brought together by the velocity gradient. Considering two colliding particles at positions x 1 and x 2 , with r defined as r = (x 2x 1 )/||x 2x 1 ||, and moving with velocity v 1 and v 2 , the relative velocity of the particles is then expected to scale as:

Δ𝑣 𝑟 ≡ (v 2 -v 1 ) • r ∼ 𝑎/𝜏 𝐾 . (8)
The quantity Δ𝑣 𝑟 (see Fig. 6(d)) is the relative velocity between the centers of mass of two colliding particles. A negative value of Δ𝑣 𝑟 means that these centers of mass are approaching each other, whereas a positive value reflects the fact that they are moving away from each other. For spherical particles, a collision is only possible when Δ𝑣 𝑟 < 0, but the rotational motion of spheroidal objects allows them to collide with Δ𝑣 𝑟 > 0.

We postulate here that the rotational degrees of freedom, i.e. the fluctuations in the orientation of the crystals, will not affect the estimate provided by Eq. ( 8) for spheres. In this spirit, Fig. 6(a-c) show the distribution of the relative velocity Δ𝑣 𝑟 between colliding crystals. To make Δ𝑣 𝑟 dimensionless, we divide it by the characteristic velocity induced by the velocity gradients when two spheroids come close to each other. We estimate this velocity as the ratio 𝑎/𝜏 𝐾 , where we recall that 𝑎 is the semi-major axis of the spheroid. The three different panels present our results at three different values of 𝜀, in increasing order of 𝜀 from (a) to (c).

Remarkably, at the lowest turbulence intensity, 𝜀 ≈ 1 cm 2 /s 3 , the PDFs of Δ𝑣 𝑟 almost perfectly superpose for 𝛽 = 0.01 (dashed black line) and 𝛽 = 0.02 (dashed-dotted red line). In fact, the center of the PDF of Δ𝑣 𝑟 at 𝛽 = 0.005 (continuous blue line) is also very similar to those at higher values of 𝛽. However, a difference can be seen by the appearance of a broad tail biased towards very negative values of Δ𝑣 𝑟 𝜏 𝐾 /𝑎, with a probability of ≈ 2 10 -3 times the maximum probability. To understand the origin of this broad tail for 𝛽 = 0.005, it is useful to recall that at this low turbulence intensity and for 𝛽 = 0.01 and 𝛽 = 0.02, the rms of the orientation fluctuations are extremely small (of the order of ∼ 3 • or less); they are slightly larger for 𝛽 = 0.005 (∼ 13 • ). This increased value of ⟨𝜑 2 ⟩ 1/2 implies a distribution of orientation that is not as sharply peaked around 𝜑 ≈ 0. As shown in Fig. 4, the settling velocity of crystals with a sufficiently small value of 𝑛 𝑧 ≡ | cos(𝜑)| (say for 𝑛 𝑧 ≲ 0.6) can exceed the averaged settling velocity by approximately 1 cm/s. This corresponds to Δ𝑣 𝑟 𝜏 𝐾 /𝑎 ≈ 22, which is close to the horizontal extent of the broad tail towards negative values of Δ𝑣 𝑟 . This horizontal extent of the blue curve in Fig. 6(a) is therefore close to the excess of settling velocity of the fastest ice crystals (those oriented with their major axis close to vertical) with respect to the averaged settling velocity. As a consequence, we interpret the broad negative tails of the PDFs as a signature of the differential settling effect: faster ice crystals, oriented predominantly "vertically", catch up with slower ones underneath, oriented "horizontally".

Further evidence for the role of differential settling can be obtained by comparing the collision kernel 𝐾 to its value 𝐾 0 in the same flow but turning off gravity (see the two last columns of Table B1). Whereas the values of 𝐾 and 𝐾 0 are very close when 𝛽 ≥ 0.01, the ratio between 𝐾 and 𝐾 0 becomes ≈ 3 for 𝛽 = 0.005, an effect that we attribute to the differential settling between crystals. We note that the values of 𝐾 0 are very close to the values of 𝐾 0 found in [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF], implying that the dynamics of orientation does not matter much to estimate the collision kernel if gravity was absent. On the contrary, in the realistic system, the much stronger alignment of the crystals with the direction of gravity in the present study leads to a strong reduction of differential settling, resulting in much lower values of 𝐾 compared to [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF].

The structure of the PDFs of Δ𝑣 𝑟 at the intermediate value of 𝜀 ≈ 16 cm 2 /s 3 , shown in Fig. 6(b), points to an increased importance of the differential settling between crystals, compared to the situation at 𝜀 ≈ 1 cm 2 /s 3 . Once again, the cores of the distributions for the three different values of 𝛽 superpose very well, which indicates that the vast majority of collisions between crystals at this intermediate turbulence intensity are due to particles being brought together by the action of the velocity gradients. In the cases 𝛽 = 0.005 (full blue line), and 𝛽 = 0.01 (dashed black line), one observes broader tails corresponding to very negative values of Δ𝑣 𝑟 . Compared to the lower turbulence intensity, 𝜀 ≈ 1 cm 2 /s 3 , the value of the probability where the broad tail develops at 𝛽 = 0.005 is much higher when 𝜀 ≈ 16 cm 2 /s 3 . This reflects the much broader distribution of crystal orientation, and in particular the larger probability of having particles with 𝑛 𝑧 < 0.6. A quantitative comparison shows that the broader tails observed for 𝛽 = 0.005 and 𝛽 = 0.01 in the distributions of Δ𝑣 𝑟 × 𝜏 𝐾 /𝑎 can be attributed to the fluctuations of orientation, with 𝑛 𝑧 ≲ 0.6. In comparison, for 𝛽 = 0.02, the distribution of Δ𝑣 𝑟 × 𝜏 𝐾 /𝑎 remains peaked around the center, without any visible deviation. The role of differential settling as a collision mechanism, clearly demonstrated for 𝛽 = 0.005 and 𝛽 = 0.01 by Fig. 6(b), can be further evidenced by comparing the collision kernel values in the presence and in the absence of gravity (Table B1): 𝐾 > 𝐾 0 for 𝛽 ≤ 0.01, whereas 𝐾 ≲ 𝐾 0 for 𝛽 ≥ 0.02.

Before discussing the last panel of Fig. 6, it is useful to notice that the particle Stokes numbers for 𝜀 ≲ 16 cm 2 /s 3 are small. The largest value, for 𝛽 = 0.02 and 𝜀 ≈ 16 cm 2 /s 3 is St ≈ 0.14. These values, which are multiplied by ≈ 4 at 𝜀 ≈ 250 cm 2 /s 3 , become 0.13 for 𝛽 = 0.005, 0.27 for 𝛽 = 0.01 and 0.54 for 𝛽 = 0.02. In this context, we recall that when the Stokes number is no longer small (∼ 0.3), the effect of particle inertia becomes prevalent, and as a consequence particles do not follow very closely the fluid velocity. This is the origin of the 'sling-effect' [START_REF] Falkovich | Acceleration of rain initiation by cloud turbulence[END_REF][START_REF] Falkovich | Sling effect in collisions of water droplets in turbulent clouds[END_REF]) (see also [START_REF] Wilkinson | Caustics in turbulent aerosols[END_REF]; [START_REF] Wilkinson | Caustic activation of rain showers[END_REF]), caused by particles colliding at a velocity Δ𝑣 𝑟 much larger than predicted by Eq. ( 8), therefore greatly enhancing the collision rates [START_REF] Gustavsson | Distribution of relative velocities in turbulent aerosols[END_REF][START_REF] Voßkuhle | Prevalence of the sling effect for enhancing collision rates in turbulent suspensions[END_REF][START_REF] Pumir | Collision aggregation due to turbulence[END_REF]. Consistent with the notion that the relative velocity between colliding particles can be greatly enhanced by the 'sling effect', the third panel of Fig. 6, corresponding to 𝜀 ≈ 250 cm 2 /s 3 , reveals in fact a broad widening of the tails of the distribution of Δ𝑣 𝑟 × 𝜏 𝐾 /𝑎 for 𝛽 = 0.01 and 𝛽 = 0.02.

At the smallest value 𝛽 = 0.005, the PDF of Δ𝑣 𝑟 × 𝜏 𝐾 /𝑎 appears as very similar to that obtained at lower turbulence intensity, in the absence of any marked effect of differential settling. In fact, the PDF of Δ𝑣 𝑟 𝜏 𝐾 /𝑎 superposes very well with those at 𝜀 ≈ 1 cm 2 /s 3 and 𝛽 ≥ 0.01, and 𝜀 ≈ 16 cm 2 /s 3 and 𝛽 = 0.02. The differential settling effect induces a very slight shoulder in the distribution at 𝛽 = 0.005 and 𝜀 ≈ 250 cm 2 /s 3 , which is to be expected since the contribution to Δ𝑣 𝑟 𝜏 𝐾 /𝑎 due to differential settling in this configuration is of order ∼ 3, and does not appear as major contribution to the distribution of relative velocity.

In comparison, the broadening of the distribution of Δ𝑣 𝑟 𝜏 𝐾 /𝑎 at 𝛽 = 0.01 and 𝛽 = 0.02 exceeds by far what is expected from differential settling. In fact, the values of 𝑈 𝑠 𝜏 𝐾 /𝑎 are of the order of 5 -6, much smaller than the extent of the PDFs shown in Fig. 6(c). This leads us to the conclusion that the mechanism leading to enhancement of the collision rate, much stronger than differential settling, is related to the generation of very large velocity differences, which we attribute to the above-mentioned effects of particle inertia [START_REF] Falkovich | Particles and fields in fluid turbulence[END_REF][START_REF] Wilkinson | Caustics in turbulent aerosols[END_REF].

Further evidence for the importance of such inertial effects can be obtained by studying the collision rate in the absence of gravity. The numerical results presented in Appendix C indicate that the collision rate for 𝑔 = 0, 𝐾 0 , scales as ∼ 𝜀 1/2 for all cases where the Stokes number St < 0.1. This scaling corresponds to the Saffman-Turner prediction, valid for weakly inertial particles. We recall that in the related problem of colliding spherical particles, inertial effects are found to be immaterial for St ≲ 0.1 [START_REF] Voßkuhle | Prevalence of the sling effect for enhancing collision rates in turbulent suspensions[END_REF]. Consistent with the results of [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF], obtained from a model that ignores the fluid inertia, the present results indicate that the same criterion applies in the case of spheroids. Namely, as shown in Table B1, inertia is found to be weak (St ≤ 0.1) in all simulations with 𝜀 ≈ 1 cm 2 /s 3 , and in those with 𝜀 ≤ 16 cm 2 /s 3 for 𝛽 ≤ 0.01. At 𝛽 = 0.02 and 𝜀 ≈ 16 cm 2 s 3 , deviations from the 𝜀 1/2 scaling are visible, see Fig. C1. The deviation from the ∼ 𝜀 1/2 is very clear at 𝜀 = 256 cm 2 /s 3 and 𝛽 ≥ 0.01, for which the Stokes numbers vary between St = 0.134 (𝛽 = 0.005) and 0.536 (𝛽 = 0.02). We observe that at 𝜀 ≈ 250 cm 2 /s 3 , the values of 𝐾 are close to the values of 𝐾 0 in the absence of gravity, providing further evidence that differential settling does not play a major role when the turbulence intensity is large. We stress that the discussion about the role of inertial effects in determining the collision rate in the absence of gravity is qualitatively completely consistent with the observations of Fig. 6.

It is worth stressing the role of fluid inertia in the problem. Comparing the results in Siewert et al. (2014b); [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF] with those of the present study clearly shows the role of fluid inertia. These effects, which are particularly significant in our study for 𝜀 ≲ 16 cm 2 /s 3 , can be traced back to the properties of preferential orientation. Taking into account the effect of fluid inertia on the orientation dynamics drastically changes the nature of the orientation bias. In particular, the distributions of orientation are extremely peaked at the smallest value of the energy dissipation, 𝜀 ≈ 1 cm 2 /s 3 . As a result, the strong dispersion of settling velocity observed in [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF] almost completely disappears when taking into account the effect of fluid inertia, except for the lightest crystals considered here (𝛽 = 0.005).

b. Representation of the different regimes in the parameter

space.

The approach proposed by [START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF] predicts a parametrization of the distribution of orientation, in particular of the variance of the angle 𝜑, as a function of two parameters, namely the Stokes number, St, and the settling number, Sv, defined by Eq. ( 6). In this subsection, we compare the prediction of that model with our own results, and discuss the implications for the various collision mechanisms introduced in Subsection 4a.

We begin by noticing that the parametrization in [START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF] actually involves the ratios St/𝐴 ⊥ and and Sv/𝐴 ⊥ , where 𝐴 ⊥ is one of the quantities defining the translational resistance tensor, defined in Appendix A. This coefficient 𝐴 ⊥ is essentially constant for the very small values of 𝛽 considered here. Quantitatively, the agreement between the fluctuations of orientation calculated in our simulations, and the predictions of the model used in [START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF], is generally quite good, see Table B1.

In the spirit of Figs. 3 and 5 of [START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF], Fig. 7 shows the variance of the orientation fluctuations, ⟨𝜑 2 ⟩, as a function of St and Sv. The figure also shows the points corresponding to the particles considered in the present work, color coded to stress the collisional mechanism involved for the corresponding runs, as discussed in Section 4a. The red points correspond to the purely Saffman-Turner mechanism. This collisional mechanism dominates when the angle fluctuations ⟨𝜑 2 ⟩ are very weak, for ⟨𝜑 2 ⟩ ≲ 10 -2 , and when St is low enough (St ≲ 0.1). The points for which the dominant collisional mechanism is provided by differential settling are marked in green. As the difference in the settling velocity is ultimately related to fluctuations in the orientation between crystals, this effect is only possible when the fluctuations of the angle 𝜑 are large enough. Judging from our numerical results, this occurs when ⟨𝜑 2 ⟩ ≳ 10 -2 (this value corresponds to an r.m.s fluctuation of slightly less than 6 • ). The case 𝛽 = 0.005 and 𝜀 ≈ 246 cm 2 /s 3 , shown by the full curve in Fig. 6(c), provides an example where the velocity differences due to differential settling were effectively immaterial, compared to ∼ 𝑎(𝜀/𝜈) 1/2 , the typical velocity differences due to turbulence. The critical value of Sv below which differential settling does not lead to appreciable velocity differences with respect to the typical velocity differences induced by turbulence can be estimated by comparing the difference in the settling velocity, and the fluctuations of the relative velocities of the colliding crystals, shown in Fig. 6. Judging from this Figure, we estimate the value of Sv where the effect of differential settling becomes immaterial to Sv ≲ 1.5. This transition region (like the two other ones) is indicated by a gray rectangle. For Sv < 1.5, we expect the Saffman-Turner mechanism to prevail. Last, the crystals for which the dominant collisional mechanism is due to inertial effects are indicated in blue, and all correspond to Stokes number St ≳ 0.1 -0.2. We have also added the point corresponding to 𝛽 = 0.05 and 𝜀 ≈ 16 cm 2 /s 3 (run 7 in Table B1), which appears as the point with the highest value of Sv. For the corresponding crystals, the Stokes number is St = 0.338, and the dominant collision mechanism is due to inertia. We notice that two of the crystals considered have a Stokes number St ≈ 0.13, right in the transition region when inertial effects become dominant. The two corresponding points have been marked by a magenta dot.

Conclusion

We have estimated numerically the rate of collisions between identical plate-like ice crystals in a turbulent flow. Our model involves a number of approximations, and is comparable with the study of [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF] with the major difference that we considered here the effect of fluid inertia. Earlier studies had demonstrated the crucial role of this physical effect in the dynamics of orientation of crystals settling in a turbulent flow, leading to the tendency for crystals to settle with their maximal dimension oriented perpendicular to gravity [START_REF] Lopez | Inertial effects on fibers settling in a vortical flow[END_REF][START_REF] Gustavsson | Theory for the effect of fluid inertia on the orientation of a small spheroid settling in turbulence[END_REF][START_REF] Sheikh | Importance of fluid inertia for the orientation of spheroids settling in a turbulent flow[END_REF][START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF], contrary to what was observed in Siewert et al. (2014a); [START_REF] Gustavsson | Statistical model for the orientation of nonspherical particles settling in turbulence[END_REF]; [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF].

Our numerical results concerning the collision statistics can be interpreted in terms of three main physical effects. The first one, as understood by [START_REF] Saffman | On the collision of drops in turbulent clouds[END_REF], is due to the effects of velocity gradients, which can bring particles together. At very large turbulence intensity, when the particle response time is larger than the characteristic times of the smallest eddies, the Kolmogorov time, particles do not exactly follow the flow, and can be "slung" out of eddies. The resulting sling effect tends to strongly increase the collision kernel of colliding particles. At low turbulence intensity, differential settling between particles of different orientations had been suggested to play a dominant role in studies that were not taking into account fluid inertia in the crystal dynamics (Siewert et al. 2014b;[START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF].

Remarkably, Fig. 7 suggests that the physical effects involved can be described by two parameters, St and Sv. This allows us to predict the main physical effects involved in the collision process in terms of these two parameters only. For this reason, we expect that the results of the present study can be expanded to describe plate-like crystals of various sizes (𝑎 and 𝛽 ≪ 1), over a broad distribution of turbulence intensity, provided the main hypotheses of the model are satisfied (size of the crystals 𝑎 smaller than the Kolmogorov scale, and small enough particle Reynolds number). On the other hand, our results leave open the question of the collision rate between objects with a very different shape, such as elongated crystals, describable in terms of prolate ellipsoids.

As demonstrated here, the three mechanisms play a role in determining the collision kernel when using the more realistic model of crystal dynamics. The main difference with earlier studies not taking into account the fluid inertia comes from the very different role of differential settling in the dynamics leading to collisions. At weak turbulence intensities, the collision rates observed in [START_REF] Jucha | Settling and collision between small ice crystals in turbulent flows[END_REF] were larger by approximately a factor of 10 compared to the collision rates observed here. This is due to the very peaked distribution of orientation dynamics at large values of Sv, which effectively reduces the dispersion of orientation, and therefore of settling velocity, as soon as the settling number is ≳ 5 -10. Using the definition of Sv, Eq. ( 6), and the values reported in Appendix B, we deduce that differential settling plays a role for 𝜀 ≳ 𝛽 4 × 1.6 10 8 cm 2 /s 3 .

The very strong bias in the orientation of crystals settling in turbulent flows is likely to also strongly affect the riming process, whereby small ice crystals collide in mixed-phase clouds with small droplets [START_REF] Naso | Collision rate of ice crystals with water droplets in turbulent flows[END_REF][START_REF] Jost | The effect of turbulence on the accretional growth of graupels[END_REF]. Further studies devoted to this problem will be reported separately.

As mentioned in the Introduction, our model is highly idealised. In particular, it disregards crystal-crystal interactions that may become important at small separations, and does not model collision outcomes. Nevertheless, we expect that the model describes at least qualitatively how particle and fluid inertia change the way in which the crystals approach (and collide). In the future we intend to test the relevance of the ghost collision approach and to formulate a more detailed collision model that can reliably predict collision outcomes. Developing such an approach, which is a necessary step to determine in quantitative terms the influence of the approximations used in the present work, is beyond the scope of this paper.

to the force and torque in a quiescent fluid are [START_REF] Brenner | The Oseen resistance of a particle of arbitrary shape[END_REF][START_REF] Cox | The steady motion of a particle of arbitrary shape at small Reynolds numbers[END_REF][START_REF] Khayat | Inertia effects on the motion of long slender bodies[END_REF][START_REF] Dabade | Effects of inertia and viscoelasticity on sedimenting anisotropic particles[END_REF]:

f (1) = -(6𝜋𝑎 ⊥ 𝜇) 3 16 𝑎𝑊 𝜈 3A-Id( Ŵ • A Ŵ ) AW , (A5a)
τ (1) = 𝐹 (𝛽)𝜇 𝑎 3 𝑊 2 𝜈

( n • Ŵ )( n× Ŵ ) . (A5b)
Here 𝑊 = |W | is the modulus of the slip velocity W = vu, Ŵ = W /𝑊 is its direction, Id is the identity matrix, and 𝐹 (𝛽) is a shape factor computed by [START_REF] Dabade | Effects of inertia and viscoelasticity on sedimenting anisotropic particles[END_REF]. For oblate spheroids with 𝛽 ≪ 1, the shape factor 𝐹 (𝛽) is positive, and tends to a finite limit of order ≈ 2.2 when 𝛽 → 0.

c. Moment of inertia tensor

The moment of inertia tensor of spheroids per unit mass, I, is given by:

𝐼 𝑖 𝑗 = 𝐼 ⊥ (

List of the runs

Table B1. Runs used to determine the collision properties between crystals. The flow properties are listed in Table 1. The number of crystals followed is 𝑁 𝑐 , 𝑇 𝑟𝑢𝑛 is the duration of the runs. The value of the averaged settling velocity is 𝑈 𝑠 . The table also lists the values of the Stokes number, 𝑆𝑡, and of the settling number, Sv, defined by Eq. ( 6). The rms of the fluctuations of orientation, ⟨ 𝜑 2 ⟩ 1/2 DNS , can be compared to the predictions from [START_REF] Gustavsson | Effect of particle inertia on the alignment of small ice crystals in turbulent clouds[END_REF], ⟨ 𝜑 2 ⟩ 1/2 stoch . We note that this model neglects the corrections of the crystal velocity due to fluid inertia, which explains much of the difference between the values reported here. The collision kernel 𝐾, defined by Eq. ( 7), and the collision kernel 𝐾 0 in the hypothetical case in the absence of gravity, are listed in the last two columns.

Runs

APPENDIX C

Collision rate in the absence of gravity

Studying collisions in the absence of gravity is a very simplified problem, not obviously appropriate in an atmospheric science context, as it completely misses the differential settling mechanism, which was shown to play an important role in Section 3b. On the other hand, the simplification is helpful in elucidating the physical mechanisms responsible for collisions. Figure C1 shows the collision kernel, as defined in a particle suspension by Eq. ( 7), in the absence of gravitational settling (the values are listed as 𝐾 0 in Table B1). The dashed line in the figure shows the 𝜀 1/2 behavior, corresponding to the Saffman-Turner prediction [START_REF] Saffman | On the collision of drops in turbulent clouds[END_REF]. The data indicates a 𝐾 0 ≈ 𝜀 1/2 behavior provided the Stokes number satisfies St ≲ 0.1, which is the case for all particles for 𝜀 ≈ 1 cm 2 /s 3 , and for particles with 𝛽 ≤ 0.01 at 𝜀 ≈ 16 cm 2 /s 3 . When St ≳ 0.1, clear deviations are visible. A similar picture was found in the case of colliding particles of spherical shapes [START_REF] Voßkuhle | Prevalence of the sling effect for enhancing collision rates in turbulent suspensions[END_REF]. . Dependence of the collision kernel, in the absence of gravitational settling, on the energy dissipation, for 𝛽 = 0.005 (right-pointing triangles), 0.01 (left-pointing triangles) and 0.02 (upper-pointing triangles). The dashed line corresponds to the dependence predicted in [START_REF] Saffman | On the collision of drops in turbulent clouds[END_REF], and provides a good description of the data when St ≲ 0.1.

Fig. 1 .

 1 Fig. 1. Representation of a planar ice crystal by an oblate spheroid. The vector n characterizing the orientation of the crystal is parallel to the axis of symmetry, and points downward. The angle 𝜑 is the angle between g and n.

Fig. 4 .

 4 Fig. 4. Averaged settling velocity conditioned on the orientation of the crystal, represented by 𝑛 𝑧 , the vertical component of n. The three curves correspond to 𝛽 = 0.005 and to the three values of 𝜀, as indicated by the legend. The horizontal gray dot-dashed line shows the settling velocity in a fluid at rest, with 𝑛 𝑧 = 1, which is the only steady and stable regime in this case.

Fig. 5 .

 5 Fig.5. Dependence of the collision kernel on: (a) the energy dissipation for 𝛽 = 0.005 (right-pointing triangles), 0.01 (left-pointing triangles) and 0.02 (upper-pointing triangles), and on (b): 𝛽 at 𝜀 ≈ 1 cm 2 /s 3 (+ symbols), 𝜀 ≈ 16 cm 2 /s 3 (diamond symbols) and 𝜀 ≈ 250 cm 2 /s 3 (square symbols). The statistical uncertainty has been evaluated using the same method as in Fig.2. The deviations are always between 2 and 7%, except for the ( 𝜀 ≈ 250, 𝛽 = 0.02) case, for which it is ∼ 9%.

Fig. 6 .

 6 Fig. 6. (a-c) Probability density function of the velocity difference between two colliding particles, defined by Eq. (8), and made dimensionless by dividing by 𝑎/𝜏 𝐾 . Increasing turbulence intensity from (a) to (c). The cores of the distributions superpose with each other. The broad tails reflect the effect of differential settling at low (𝜀 ≈ 1 cm 2 /s 3 , top) or at intermediate (𝜀 ≈ 16 cm 2 /s 3 , middle) turbulence intensity, or the sling effect, at the highest turbulence intensity (𝜀 ≈ 250 cm 2 /s 3 , bottom). (d) Illustration of the definition of the relative velocity Δ𝑣 𝑟 [Eq. (8)].

Fig. 7 .

 7 Fig. 7. Summary of the various regimes as a function of the two main parameters of the problem, namely the settling number, Sv (horizontal) and the Stokes number, St (vertical). The yellow dashed lines indicate the constant 𝛽-lines, and the orange dashed-dotted lines the constant 𝜀lines. The color code on the right corresponds to the fluctuations of the angle of orientation as a function of (Sv, St), as discussed in Gustavsson et al. (2021). The points used in the simulation are color coded according to the prevalent physical effect determining the collision: the dominant collision mechanisms is the Saffman-Turner mechanism (red points), the differential settling (green points), and the sling effect (blue points). The two magenta points correspond to transition regimes between several dominant collision mechanisms, which occur in the domains of parameters in transparent gray.

  Fig.C1. Dependence of the collision kernel, in the absence of gravitational settling, on the energy dissipation, for 𝛽 = 0.005 (right-pointing triangles), 0.01 (left-pointing triangles) and 0.02 (upper-pointing triangles). The dashed line corresponds to the dependence predicted in[START_REF] Saffman | On the collision of drops in turbulent clouds[END_REF], and provides a good description of the data when St ≲ 0.1.

  table. The Taylorscale Reynolds number 𝑅𝑒 𝜆 , the Kolmogorov time 𝜏 𝐾 , the integral time 𝑇 𝐿 , and the rms velocity fluctuation 𝑢 𝑟 𝑚𝑠 are listed for the three flows. The resolution for each run is 𝑁 3 .

	Flow	I	II	III
	𝜀 (cm 2 /s 3 )	0.976 15.62 246.4
	𝑅𝑒 𝜆	55.8	94.6	151.2
	𝜏 𝐾 (s)	0.341 0.085 0.021
	𝑇 𝐿 (s)	1.96	0.70	0.26
	𝑢 𝑟 𝑚𝑠 (cm/s)	2.18	5.72	14.4
	𝑁	384	768	1536

  𝛿 𝑖 𝑗 -𝑛 𝑖 𝑛 𝑗 ) + 𝐼 ∥ 𝑛 𝑖 𝑛 𝑗 , with 𝐼 ⊥ =

	1 + 𝛽 2 5	𝑎 2 ⊥ and 𝐼 ∥ =	2 5	𝑎 2 ⊥ .	(A6)
	APPENDIX B

Acknowledgments. BM thanks Fabien Candelier for discussions regarding the inertial torque. MZS, AP and AN acknowledge support from the IDEXLYON project (Contract ANR-16-IDEX-0005) under University of Lyon auspices. MZS's PhD thesis was supported by a grant from the Higher Education Commission Pakistan. KG and BM were supported by the grant Bottlenecks for particle growth in turbulent aerosols from the Knut and Alice Wallenberg Foundation, Dnr. KAW 2014.0048, and in part by VR grant no. 2017-3865 and Formas grant no. 2014-585. Computational resources were provided by PSMN, C3SE and SNIC.

Data availability statement. The dataset on which this paper is based is too large to be retained or publicly archived with available resources. Documentation and methods used to support this study are available on request to the authors.

APPENDIX A

Equations of motion

We provide here explicit expressions for the force, f ℎ , and for the torque, τ ℎ , acting on an oblate spheroid of semi-major axis 𝑎, semi-major axis 𝑐 and aspect ratio 𝛽 = 𝑐/𝑎 < 1.

As explained in the text, we express f ℎ and τ ℎ as the sum of the force and torque, f (0) ℎ and τ (0) ℎ , obtained in the Stokes limit, by neglecting all the nonlinear terms in the Navier-Stokes equations Eq. ( 1) and valid when the particle Reynolds number 𝑅𝑒 𝑝 is ≪ 1, plus a correction induced by taking into account the fluid inertia in a perturbative expansion in powers of 𝑅𝑒 𝑝 , f (1) ℎ and τ (1) ℎ .

a. Force and torque in the Stokes limit

The expression for the force f (0) ℎ reads:

where 𝜇 = 𝜌 𝑓 𝜈 is the dynamic viscosity of the fluid, A is the resistance tensor, and u is the fluid velocity at the particle location. For a spheroid, the expression of the resistance tensor is known exactly, see e.g. [START_REF] Kim | Microhydrodynamics: principles and selected applications[END_REF], and is recalled here for convenience. Namely, the translational resistance tensor A is equal to:

where we recall that n is a unit vector parallel to the axis of symmetry of the spheroid, and with coefficients

, where 𝛾 =

These expressions are consistent with those given in Tables 3.4 and 3.6 in [START_REF] Kim | Microhydrodynamics: principles and selected applications[END_REF]. In the case of interest here, the aspect ratio 𝛽 is ≪ 1, so 𝛾 ≈ 𝜋/(2𝛽), 𝐴 ⊥ ≈ 16/(9𝜋) and 𝐴 ∥ ≈ 8/(3𝜋). This implies that in a still fluid, crystals settling with their maximal dimension oriented vertically are 3/2 faster than crystals with their major axis oriented horizontally. Note that in a fluid at rest, only the horizontal configuration is stable.

The expression for the torque τ (0) ℎ , named after Jeffery [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF], is:

Here ω -𝛀 is the angular slip velocity, and 𝛀 = 1 2 ∇ × u is half the fluid vorticity at the particle position. It is related to the asymmetric part O of the matrix of fluid-velocity gradients by the relation Or = 𝛀 × r. The symmetric part of the matrix of fluid-velocity gradients is denoted by S. The resistance tensors C and H in Eq.(A3) have the following explicit expressions:

with

Here 𝜀 𝑖 𝑗𝑙 is the antisymmetric tensor, and repeated indices are summed over.

b. Force and torque due to fluid inertia

The expressions of the contributions due to fluid inertia are estimated by considering a particle falling through a quiescent fluid with a steady settling velocity. The slip generates fluid accelerations, and acts as a homogeneous background flow. To leading order in the particle Reynolds number, the resulting steady convective-inertia corrections