Statistical model for the orientation of non-spherical particles settling in turbulence

TL;DR

This paper presents a Gaussian statistical model to predict the orientation distribution of small anisotropic particles settling in turbulence, aiding understanding of cloud microphysics and related phenomena.

Contribution

The paper introduces a simple Gaussian model that accurately predicts particle orientation distributions considering inertia, advancing the understanding of particle behavior in turbulent flows.

Findings

Model accurately predicts orientation distribution

Applicable to ice-crystals in clouds

Enhances understanding of cloud microphysics

Abstract

The orientation of small anisotropic particles settling in a turbulent fluid determines some essential properties of the suspension. We show that the orientation distribution of small heavy spheroids settling through turbulence can be accurately predicted by a simple Gaussian statistical model that takes into account particle inertia and provides a quantitative understanding of the orientation distribution on the problem parameters when fluid inertia is negligible. Our results open the way to a parameterisation of the distribution of ice-crystals in clouds, and potentially leads to an improved understanding of radiation reflection, or particle aggregation through collisions in clouds.

Equations

$$m\ddot{\mbox{\boldmath$x$}}=\mbox{\boldmath$f$}+mg\hat{\mbox{\boldmath$g$}}\,,\quad\dot{\mbox{\boldmath$n$}}=\mbox{\boldmath$\omega$}\wedge\mbox{\boldmath$n$}\,,\quad\tfrac{{\rm d}}{{\rm d}t}\big{(}\mathbb{J}(\mbox{\boldmath$n$})\mbox{\boldmath$\omega$}\big{)}=\mbox{\boldmath$T$}\,.$$

$$\displaystyle\begin{bmatrix}\mbox{\boldmath$f$}\\ \mbox{\boldmath$T$}\end{bmatrix}=m\gamma\begin{bmatrix}\mathbb{M}^{(t)}&0&0\\ 0&\mathbb{M}^{(r1)}&\mathbb{M}^{(r2)}\end{bmatrix}\begin{bmatrix}\mbox{\boldmath$u$}-\mbox{\boldmath$v$}\\ \mbox{\boldmath$\Omega$}-\mbox{\boldmath$\omega$}\\ \mathbb{S}\end{bmatrix}.$$

$$\displaystyle\mbox{\boldmath$n$}_{t}$$

$$\displaystyle+\operatorname{Ku}\int_{0}^{t}{\rm d}t_{1}(1-{\rm e}^{C^{{(r1)}}_{\perp}(t_{1}-t)/\operatorname{St}})\big{[}\mathbb{O}(t_{1})\mbox{\boldmath$n$}_{0}+\tfrac{\lambda^{2}-1}{\lambda^{2}+1}\big{(}\mathbb{S}(t_{1})\mbox{\boldmath$n$}_{0}-(\mbox{\boldmath$n$}_{0}^{\sf T}\mathbb{S}(t_{1})\mbox{\boldmath$n$}_{0}\big{)}\mbox{\boldmath$n$}_{0})\big{]}$$

$$\displaystyle+\operatorname{Ku}^{2}\int_{0}^{t}{\rm d}t_{1}\int_{0}^{t_{1}}{\rm d}t_{2}\big{[}\mbox{\boldmath$c$}^{(OO)}_{ijkl}(\mbox{\boldmath$n$}_{0};t,t_{1},t_{2})O_{ij}(t_{1})O_{kl}(t_{2})+\mbox{similar terms with $OS$ and $SS$}\big{]}\,.$$

$$\displaystyle\mbox{\boldmath$v$}_{s}(\mbox{\boldmath$n$}_{0})$$

$$\displaystyle\langle(\mbox{\boldmath$n$}_{t}\cdot\hat{\mbox{\boldmath$g$}})^{p}\rangle_{\raisebox{-0.28453pt}{$\scriptstyle\mbox{\boldmath$n$}_{0}$}}\!\!\!=\!(\mbox{\boldmath$n$}_{0}\cdot\hat{\mbox{\boldmath$g$}})^{p}\!+\!p\operatorname{Ku}(\mbox{\boldmath$n$}_{0}\cdot\hat{\mbox{\boldmath$g$}})^{p-1}\langle\mbox{\boldmath$n$}_{t}^{(1)}\!\cdot\hat{\mbox{\boldmath$g$}}\rangle_{\raisebox{-0.28453pt}{$\scriptstyle\mbox{\boldmath$n$}_{0}$}}$$

$$\displaystyle\!+\!\tfrac{p}{2}\operatorname{Ku}^{2}(\mbox{\boldmath$n$}_{0}\cdot\hat{\mbox{\boldmath$g$}})^{p-2}\big{\langle}2(\mbox{\boldmath$n$}_{0}\cdot\hat{\mbox{\boldmath$g$}})(\mbox{\boldmath$n$}_{t}^{(2)}\!\cdot\hat{\mbox{\boldmath$g$}})\!+\!(p\!-\!1)(\mbox{\boldmath$n$}_{t}^{(1)}\!\cdot\hat{\mbox{\boldmath$g$}})^{2}\big{\rangle}_{\raisebox{-0.28453pt}{$\scriptstyle\mbox{\boldmath$n$}_{0}$}}$$

$$\displaystyle\langle(\mbox{\boldmath$n$}\cdot\hat{\mbox{\boldmath$g$}})^{2p}\rangle_{\infty}\!=\!\frac{1}{2p\!+\!1}\!+\!\sum_{i=1}^{\infty}\frac{\operatorname{G}^{2i}\sum_{j=1}^{i}p^{j}A^{(2i)}_{j}(\operatorname{St},\lambda)}{\prod_{k=1}^{i+1}(2p+2k-1)}\,.$$

$$P(n_g)$$

Full text

How non-spherical objects settle in a turbulent environment is a highly relevant question in several domains. An example is provided by very small ice crystals in clouds (size $\sim$100\text{,}\mathrm{\SIUnitSymbolMicro m}), which grow through aggregation to form precipitation size particles (size $\sim$1\text{\,}\mathrm{mm}) Pruppacher and Klett (1997); Cho et al. (1981); Chen and Lamb (1994); Hubbert et al. (2014). The settling of plankton in the ocean Ruiz et al. (2004); Cencini et al. (2013); Gustavsson et al. (2016) can induce patchiness of the population, therefore affecting mating, feeding and predation Guasto et al. (2012). In these problems, the orientational degrees of freedom clearly affect not only settling and collision properties, but also light reflection Yang et al. (2015). As a prerequisite to a description of these effects, this Letter provides an understanding of the orientation statistics of small spheroids settling in a turbulent environment based on a statistical model, under the assumption that fluid inertia can be neglected.

The interaction between turbulence and settling leads to intriguing phenomena, even in the simpler case of spherical particles. Maxey found that turbulence increases the settling speed of a single small particle Maxey (1987); Good et al. (2014). Substantial progress was recently achieved in understanding how two spherical particles settling together move relative to each other and collide  Gustavsson et al. (2014a); Bec et al. (2014); Ireland et al. (2016); Mathai et al. (2016); Parishani et al. (2015).

In a fluid at rest the orientation dynamics of slowly settling non-spherical particles is determined by weak torques resulting from fluid inertia Khayat and Cox (1989); Dabade et al. (2015); Candelier and Mehlig (2016); Roy et al. (2016). Turbulence affects the orientation of such particles through turbulent vorticity and strain. In the absence of settling this is well understood Jeffery (1922); Pumir and Wilkinson (2011); Parsa et al. (2012); Chevillard and Meneveau (2013); Gustavsson et al. (2014b); Byron et al. (2015); Zhao et al. (2015); Voth (2015); Voth and Soldati (2017); Gustavsson et al. (2016). Neglecting fluid inertia, the direct numerical simulations (DNS) of turbulence by Siewert et al. Siewert et al. (2014) demonstrated that settling induces a bias in the orientation distribution of the particles. The physical origin of this bias is not known, and it is not understood how the bias depends on the parameters of the problem: the turbulent Reynolds number, ${\rm Re}_{\lambda}$, the Stokes number (particle inertia), the gravitational acceleration, and the particle shape. Also, how significant are non-Gaussian, intermittent small-scale features of the turbulent flow Schumacher et al. (2014), such as intense vortex tubes Cho et al. (1981) in aligning the particles?

To answer these questions we analyse a statistical model for the orientation of small heavy spheroids settling in homogeneous isotropic turbulence, for parameters relevant to cloud physics, and compare with results based on DNS of turbulence. Fig. 1 shows the predicted bias in the distribution of the vector $n$ pointing along the particle symmetry axis. The statistical-model predictions agree very well with the DNS results. This shows that that non-Gaussian turbulent fluctuations are not important. The statistical model explains the sensitive parameter dependence of the DNS results. This is important because it allows us to parameterise the bias, to quantitatively understand the physical properties of the system.

We analyse the model by an expansion in the ‘Kubo number’ ${\rm Ku}$, a dimensionless correlation time of the flow Gustavsson and Mehlig (2016). Padé-Borel resummation yields excellent agreement with numerical simulations at $\operatorname{Ku}=0.1$, and qualitative agreement with DNS of turbulence. At larger $\operatorname{Ku}$ the theory fails to converge, but the model still explains qualitatively the underlying mechanisms. Last, we discuss possible effects of fluid inertia.

Formulation of the problem. The equations of motion for translation and rotation of a particle reads

[TABLE]

Here $g$ is the gravitational acceleration (direction \hat{\mbox{\boldmathg}}), $x$ is the position of the particle, $n$ its symmetry vector, $m$ its mass, $\omega$ its angular velocity, and \mathbb{J}(\mbox{\boldmathn}) is its inertia tensor in the lab frame. In the point-particle approximation, force $f$ and torque $T$ on a spheroid are Kim and Karrila (1991); Marchioli et al. (2010); Gustavsson et al. (2014b):

[TABLE]

In Eq. (2), $v$ is the particle velocity, \mbox{\boldmathu}(\mbox{\boldmathx},t) is the turbulent velocity field, \mbox{\boldmath\Omega}\equiv\tfrac{1}{2}\mbox{\boldmath\nabla}\wedge\mbox{\boldmathu} is half the turbulent vorticity, $\mathbb{S}$ is the strain-rate matrix, the symmetric part of the matrix $\mathbb{A}$ of fluid-velocity gradients (its antisymmetric part is called $\mathbb{O}$), and $\mathbb{M}$ are translational and rotational resistance tensors: \mathbb{M}^{(t)}\equiv C^{(t)}_{\perp}\mathbb{I}+(C^{(t)}_{\parallel}\!-\!C^{(t)}_{\perp})\mbox{\boldmathn}\mbox{\boldmathn}^{\sf T}, \mathbb{M}^{(r1)}\equiv{K}^{{(r1)}}_{\perp}\mathbb{I}+({K}^{{(r1)}}_{\parallel}\!-\!{K}^{{(r1)}}_{\perp})\mbox{\boldmathn}\mbox{\boldmathn}^{\sf T}, and $\mathbb{M}^{(r2)}$ is a third-rank tensor. For a fore-aft symmetric particle, the equations of motion (1,2) are invariant under \mbox{\boldmathn}\to-\mbox{\boldmathn}, so that only the magnitude n_{g}\equiv|\mbox{\boldmathn}\cdot\hat{\mbox{\boldmathg}}| can play a role in the dynamics. The form of $\mathbb{M}^{(r2)}$ and of the $C$ and $K$-coefficients are known for spheroidal particles, see Supplemental Material (SM) SM and Ref. Fries et al. (2017). The parameter $\gamma\equiv 9\nu\rho_{f}/(2{a_{\parallel}a_{\perp}}\rho_{p})$ is Stokes constant, $\nu$ is the kinematic viscosity of the fluid, $\rho_{f}$ and $\rho_{p}$ are fluid and particle mass densities, $2a_{\parallel}$ is the length of the particle symmetry axis, and $2a_{\perp}$ is the particle diameter.

Our DNS of turbulence use the code described in Voßkuhle et al. (2014) and in the SM SM . The Kolmogorov scales $u_{\rm K}$, $\eta_{\rm K}$, and $\tau_{\rm K}$ are determined by the dissipation rate $\varepsilon\equiv\nu\langle\operatorname{Tr}\mathbb{A}\mathbb{A}^{\sf T}\rangle$ (the average is along steady-state Lagrangian trajectories), and by \nu\approx$$1\text{\times}{10}^{-5}\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1} (air). The particle aspect ratio is $\lambda\equiv a_{\parallel}/a_{\perp}$. The simulations were done for spheroids of varying $\lambda$ and with {{\rm max}(a_{\parallel},a_{\perp})}=$$150\text{\,}\mathrm{\SIUnitSymbolMicro m}, much smaller than $\eta_{\rm K}$ for values of $\varepsilon$ pertaining to mixed-phase clouds (DNS: $\varepsilon\!\approx\!1,\,16,$ and $256\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{s}}^{-3}$). Particle inertia is measured by the Stokes number $\operatorname{St}_{\rm K}\equiv(\gamma\tau_{\rm K})^{-1}$. The mass-density ratio is ${{\rho_{p}/\rho_{f}}}\approx 1000$ (ice crystals in air), and the dimensionless gravity parameter is defined as ${\rm F}_{\rm K}\equiv g\tau_{\rm K}/u_{\rm K}$.

Statistical model. The model is appropriate for particles smaller than $\eta_{\rm K}$. We approximate the universal Schumacher et al. (2014) dissipative-range turbulent fluctuations by an incompressible, homogeneous, isotropic Gaussian random velocity field \mbox{\boldmathu}(\mbox{\boldmathx},t) with zero mean, correlation length $\ell$, correlation time $\tau$, and rms speed $u_{0}$ Gustavsson and Mehlig (2016) (details given in the SM SM ). In the persistent limit Gustavsson and Mehlig (2016), for $\operatorname{Ku}\equiv u_{0}\tau/\ell>1$, the model parameters $\operatorname{St}\equiv(\gamma\tau)^{-1}$ and ${\rm F}\equiv g\tau/u_{0}$ map to $\operatorname{St}_{\rm K}=\sqrt{5}\operatorname{Ku}\operatorname{St}$ and $F_{\rm K}=[F/(5\operatorname{Ku})]\ell/\eta_{\rm K}$. Here $\ell/\eta_{\rm K}$ is the ratio between the size of the dissipation range and the Kolmogorov length. In turbulence this ratio depends weakly on the Reynolds number ${\rm Re}_{\lambda}$ Calzavarini et al. (2009), $\ell=c\eta_{\rm K}{\rm Re}_{\lambda}^{1/2}$. For the data shown in Fig. 1 we have ${\rm Re}_{\lambda}=95$, and Fig. S1 in SM SM shows results for other values of ${\rm Re}_{\lambda}$. We find good agreement between the statistical-model results at large $\operatorname{Ku}$ and the DNS for $c\approx 1.3$. For $\operatorname{Ku}>1$, the model predictions depend on two parameter combinations only Gustavsson and Mehlig (2016), $\operatorname{Ku}\operatorname{St}$ and ${\rm F}/\operatorname{Ku}$. In terms of the DNS parameters this means that the orientation bias depends only on $\operatorname{St}_{\rm K}$ and $F_{\rm K}{\rm Re}_{\lambda}^{-1/2}$.

Perturbation theory. Eqs. (1,2) are solved by expansion in powers of $\operatorname{Ku}$ Gustavsson and Mehlig (2011, 2016). We outline the essential steps below, details are given in the SM SM . We use dimensionless variables: t^{\prime}\!\equiv\!t/\tau,\mbox{\boldmathr}^{\prime}\!\equiv\!\mbox{\boldmathr}/\ell,\mbox{\boldmathu}^{\prime}\!\equiv\!\mbox{\boldmathu}/u_{0}, and drop the primes. To calculate the steady-state distribution of n_{g}\equiv|\mbox{\boldmathn}\cdot\hat{\mbox{\boldmathg}}| we must evaluate the fluctuations of the fluid-velocity gradients along particle paths. This is achieved by an expansion in \delta\mbox{\boldmathx}_{t}\equiv\mbox{\boldmathx}_{t}-\mbox{\boldmathx}^{({\rm d})}_{t} around the deterministic solution \mbox{\boldmathx}^{({\rm d})}_{t} of Eqs. (1,2) for \mbox{\boldmathu}=0. This gives expansions in powers of $\operatorname{Ku}$ Gustavsson and Mehlig (2016):

[TABLE]

The matrices $\mathbb{O}$ and $\mathbb{S}$ are evaluated along deterministic paths \mbox{\boldmathx}_{t}^{({\rm d})}=\mbox{\boldmathx}_{0}+\mbox{\boldmathv}_{s}(\mbox{\boldmathn}_{0})\,t (Fig. 2a) with settling velocity

[TABLE]

Also, $C^{(r1)}_{\perp}\equiv 5\lambda K^{(r1)}_{\perp}/(a_{\parallel}a_{\perp}(\lambda^{2}+1))$. Eq. (4) is the lowest-order solution of Eqs. (1,2). The terms in Eq. (Statistical model for the orientation of non-spherical particles settling in turbulence) that do not involve \delta\mbox{\boldmathx}_{t} depend only on the history of the fluid-velocity gradients along the paths \mbox{\boldmathx}_{t}^{({\rm d})} (‘history contribution’). The $c$-coefficients contain at most five powers of \mbox{\boldmathn}_{0}, and one must sum over all tensor products allowed by symmetry (Einstein convention). See SM SM .

The first integral shown in Eq. (Statistical model for the orientation of non-spherical particles settling in turbulence), by contrast, depends on \delta\mbox{\boldmathx}_{t}. It is therefore sensitive to how turbulence modifies the settling paths (‘preferential sampling’ Gustavsson and Mehlig (2016)).

We determine the steady-state moments \langle(\mbox{\boldmathn}_{t}\cdot\hat{\mbox{\boldmathg}})^{p}\rangle_{\infty} by first calculating the moments conditional on the initial orientation \mbox{\boldmathn}_{0}, using Eq. (Statistical model for the orientation of non-spherical particles settling in turbulence) and the relation

[TABLE]

where \mbox{\boldmathn}_{t}^{(i)} is the coefficient of $\operatorname{Ku}^{i}$ in Eq. (Statistical model for the orientation of non-spherical particles settling in turbulence). Eq. (5) is valid to order $\operatorname{Ku}^{2}$. We average over the fluid-velocity fluctuations as described in Ref. Gustavsson and Mehlig (2016). The moments are independent of the initial position \mbox{\boldmathx}_{0} due to homogeneity of the flow. We expect that effects of the initial velocity \mbox{\boldmathv}_{0} and angular velocity \mbox{\boldmath\omega}_{0} decay exponentially, so that they do not affect the steady state. We therefore set both to zero. Only the \mbox{\boldmathn}_{0}-dependence matters. In this way we obtain expressions for \langle(\mbox{\boldmathn}_{t}\cdot\hat{\mbox{\boldmathg}})^{p}\rangle_{\raisebox{-0.28453pt}{\scriptstyle\mbox{\boldmath$n$}_{0}}}, which involve secular terms that increase linearly with time as $t\to\infty$. But these terms must vanish since \mbox{\boldmathn}_{t} is a unit vector. This condition yields a recursion relation for the steady-state averages \langle(\mbox{\boldmathn}\cdot\mbox{\boldmath\hat{\mbox{\boldmath$g$}}})^{p}\rangle_{\infty}, independent of \mbox{\boldmathn}_{0}. This recursion is valid for arbitrary values of $\operatorname{G}\equiv\operatorname{Ku}\,{\rm F}\operatorname{St}/C^{(t)}_{\perp}$, and to order $\operatorname{Ku}^{0}$. Note that $\operatorname{G}$ can be large even if $\operatorname{Ku}$ is small. We solve the recursion by a series expansion in small $\operatorname{G}$:

[TABLE]

The coefficients $A^{(2i)}_{j}(\operatorname{St},\lambda)$ depend on the shape and inertia of the particle, but not on $\operatorname{G}$ or $p$. From Eq. (6) we obtain the Fourier transform of the probability distribution of n_{g}=|\mbox{\boldmathn}\cdot\hat{\mbox{\boldmathg}}|. Inverse Fourier transformation yields the distribution. To order $\operatorname{G}^{4}$ we find:

[TABLE]

The lowest-order term corresponds to a uniform distribution of \mbox{\boldmathn}_{t}. Let us examine the $\operatorname{G}^{2}$-term. It turns out that $A_{1}^{(2)}$ is negative for disks and positive for rods (see Fig. S2 in the SM SM ). This explains that the orientation of settling disks is biased: disks tend to fall edge on and rods settle tip first (as in Fig. 1).

Padé-Borel resummation. Now consider higher orders in the $\operatorname{G}$-expansion. The series (6) is asymptotically divergent and must be resummed. Fig. 2 demonstrates that Padé-Borel resummation Bender and Orszag (1978); Gustavsson and Mehlig (2016) of the series yields excellent results. Shown are results from a resummation of (6) to order $\operatorname{G}^{34}$ (thick solid lines). These results agree very well with numerical simulations of the statistical model for $\operatorname{Ku}=0.1$ and $\operatorname{St}=10$ (symbols). The resummed theory works up to ${\rm G}=10$, and in this range the bias increases with increasing ${\rm G}$. The resummed theory also predicts that the moments increase as $\operatorname{St}$ increases, for fixed $\operatorname{G}$. A more detailed analysis of the recursion leading to Eq. (6) reveals, however, that the limit ${\rm G}\to\infty$ is delicate. Perfect alignment requires $\lambda=\infty$ SM .

In summary, perturbation theory in $\operatorname{Ku}$ shows that turbulence gives rise to an orientation bias (Fig. 2), in excellent agreement with statistical-model simulations at $\operatorname{Ku}=0.1$ and in qualitative agreement with DNS (Fig. 1).

The calculations leading to Eq. (6) reveal that each moment \langle(\mbox{\boldmathn}_{t}\cdot\hat{\mbox{\boldmathg}})^{2p}\rangle_{\infty} is a sum of two contributions that stem from the ‘preferential sampling’ and ‘history’ terms in Eq. (Statistical model for the orientation of non-spherical particles settling in turbulence). For small $\operatorname{Ku}$ the history effect is dominant, the orientation bias is entirely determined by the history of fluid-velocity gradients along straight deterministic paths, Fig. 2a. Decomposing the leading-order contribution as $A_{1}^{(2)}\!=\!A_{1,\rm pref.}^{(2)}\!+\!A_{1,\rm hist.}^{(2)}$ we find that $|A_{1,\rm pref.}^{(2)}|\ll|A_{1,\rm hist.}^{(2)}|$ (Fig. S2b in the SM SM ). Fig. 3a leads to the same conclusion. It shows the distribution $P(n_{g})$ for $\operatorname{Ku}=0.1$. Also shown is $P(n_{g})$ computed for particles falling with constant velocity \mbox{\boldmathv}=\mbox{\boldmathv}_{\rm s}(\mbox{\boldmathn}_{0}). We choose the squared initial orientation \mbox{\boldmathn}_{0}\mbox{\boldmathn}_{0}^{\sf T} in (4) as the steady state average \langle\mbox{\boldmathn}_{0}\mbox{\boldmathn}_{0}^{\sf T}\rangle_{\infty}, evaluated using the small-$\operatorname{Ku}$ theory. This corresponds to keeping just the history contribution to $\langle n_{g}^{2p}\rangle_{\infty}$. We observe excellent agreement with the full statistical-model simulations. This shows that the history effect causes the orientation bias at small values of $\operatorname{Ku}$.

Persistent limit. In the persistent limit we use numerical simulations with $\operatorname{Ku}=10$ to analyse the orientation bias in the same way as for small $\operatorname{Ku}$. The result is shown in Fig. 3b (parameters correspond to two curves in Fig. 1(a). We plot the full statistical-model distribution and results for particles with a constant velocity (4) that neglects preferential sampling. For the data in Fig. 3b, the average \langle\mbox{\boldmathn}_{0}\mbox{\boldmathn}_{0}^{\sf T}\rangle_{\infty} is computed using statistical-model simulations. We see that the history effect makes a substantial contribution to $P(n_{g})$. But since the distributions do not match, we infer that preferential sampling also contributes. This contribution is hatched in Fig. 3b.

Limit of large settling speeds. Fig. 3c shows the moments $\langle n_{g}^{2p}\rangle_{\infty}$ for $p=1,2,3$ in the persistent limit as functions of the DNS Stokes number $\operatorname{St}_{\rm K}$. Open symbols denote full statistical-model simulations, solid lines correspond to simulations based on straight deterministic paths. At intermediate Stokes numbers we see a clear difference between the two simulations, preferential sampling is important in this region.

As the Stokes number grows, however, the Figure demonstrates that preferential sampling ceases to play a role. In this limit the orientation bias is entirely caused by the history effect. The bias shown in Fig. 3c increases as $\operatorname{St}_{\rm K}$ increases. But as the perturbation theory indicates, the limit of large $\operatorname{G}$ is quite subtle. Statistical-model simulations for $\operatorname{Ku}=1$ show that the degree of alignment starts to decrease for very large $\operatorname{G}$.

Conclusions. We analysed a statistical model for the orientational dynamics of small heavy spheroids settling in turbulence. The predictions of the model agree well with our own numerical results based on DNS of homogeneous isotropic turbulence (Fig. 1). Our statistical-model analysis shows that there are two distinct competing mechanisms causing the orientation bias: preferential sampling and the history effect. The latter dominates for large settling speeds, but it makes substantial contributions also in other parameter regimes. Preferential sampling dominates only when the bias is negligibly small. When the bias is significant, the history effect explains at least about 50% of the bias observed in Fig. 1.

We have shown that the orientation alignment depends on combinations of dimensionless numbers: $\operatorname{St}_{\rm K}$ and $F_{\rm K}{\rm Re}_{\lambda}^{-1/2}$. Our analysis shows that it is the small-scale properties of the flow that determine the orientation alignment. The ${\rm Re}_{\lambda}$-dependence arises only because it determines the ratio between the smooth scale $\ell$ to $\eta_{\rm K}$. We note that $F_{\rm K}{\rm Re}_{\lambda}^{-1/2}$ equals the ratio of the settling velocity and the rms turbulent velocity fluctuations.

Our results pertain to small ice crystals settling in turbulent clouds, and allow us to model the sensitive dependence of the effect upon particle shape, size, and the turbulence intensity. This is important since turbulent dissipation rates vary widely in clouds. Our results predict strongly varying degrees of alignment. That the statistical model is in excellent agreement with the DNS opens a way to parameterise the orientation distribution of ice-crystals in clouds. This potentially leads to an improved understanding of the radiative properties of clouds, and of particle aggregation through collisions in clouds.

The present work is based on the point-particle approximation of heavy particles, which neglects the effect of fluid inertia. This requires the particle Reynolds number ${\rm Re}_{p}\equiv av_{c}/\nu$ to be small, where $a={\rm max(a_{\parallel},a_{\perp})}$. Estimating the slip velocity $v_{c}$ by the Stokes settling speed, we find that ${\rm Re}_{p}$ is of order unity for the data shown in Fig. 1, so the condition is marginally satisfied. The shear Reynolds number, ${\rm Re}_{s}$, must also be small. Since ${\rm Re}_{s}\equiv a^{2}\sqrt{\langle{\rm tr}\,\mathbb{S}^{2}\rangle}/\nu\sim(a/\eta_{\rm K})^{2}$ Candelier et al. (2016), this condition is satisfied for small particles.

Lopez et al. Lopez and Guazzelli (2017) analysed the orientational dynamics of rods settling in a vortical flow. For small ${\rm Re}_{p}$ they found a bi-modal distribution, with peaks at $n_{g}=0$ and $1$. They explain the peak at $n_{g}=0$ by the effect of fluid inertia. Our results may explain the peak at $n_{g}=1$. These results, although not for a turbulent flow, indicate that turbulent and fluid-inertia torques compete in general. How to model this competition is an open question. For small Stokes numbers one may formulate an ad-hoc model by simply adding turbulent and fluid-inertia torques, along the lines suggested in Ref. Lopez and Guazzelli (2017). But in general it remains a challenge to take into account effects due to fluid inertia from first principles, in a turbulent environment. Simulations resolving particle and fluid motion Homann and Bec (2010); Fornari et al. (2016) and experiments Traugott and Liberzon (2017); Variano (2016); Marcus et al. (2014); Kramel et al. (2016) for micron-sized particles in turbulence are needed to test the predictions, and to determine the orientational dynamics of larger particles where fluid inertia must matter Kramel et al. (2016). Finally, how to extend the ideas developed here to particles lighter than the fluid remains a challenging task.