Effects of turbulence on the geometric collision rate of sedimenting droplets. Part 2. Theory and parameterization

There is a typo in equation (76) of the paper which is corrected here. We also wish to take this opportunity to comment on the numerical equivalence of equation (80) and the alternative form mentioned on page 37. Please see the PDF for details.


Introduction
Understanding the geometric collision rate is a necessary step in determining the collision-coalescence growth of cloud droplets [1]. In part 1 [2], direct numerical simulations (DNS) were used as numerical experiments to better understand the interaction of sedimenting cloud droplets with three-dimensional, time-dependent turbulent motion at the dissipation-range scales. DNS simultaneously provided dynamic collision kernel and kinematic pair statistics, allowing direct validation of a kinematic description of the geometric collision kernel. DNS results revealed that both the radial relative velocity and the radial distribution function (RDF) could be moderately enhanced by air turbulence. However, the level of enhancement depends, in a complex manner, on the size of droplets (which in turn determines the response time and the terminal velocity) and the strength of air turbulence (i.e. the dissipation rate, Reynolds number, etc).
In this paper, our objective is to develop an analytical parameterization for the turbulent geometric collision kernel. There are several motivations for this. Firstly, DNS are expensive so the number of runs are necessarily limited by the available computational resources. At the same time, even for a relatively simple bidisperse system, the collision rate of cloud droplets depends on a large number of droplet parameters and air flow characteristic scales. Secondly, DNS has its limitation in resolving the range of length and time scales, which translates to the very limited range of flow Reynolds numbers it can handle. For example, DNS alone will not be able to resolve the effect of flow Reynolds number on the collision rate. Most importantly, to include the effects of air turbulence on the growth of cloud droplets by collision-coalescence, a theoretical parameterization of the enhanced collision kernel is desired. Carefully developed theoretical parameterization can also improve our understanding of the complex interaction of droplets with air turbulence, possibly over a much wider parameter range than is achievable in DNS. 3 A comprehensive literature survey in section 2 reveals that no general theory is available to predict the radial relative velocity for finite-inertia, sedimenting droplets in turbulent flow. Previous theoretical treatments often greatly simplify the problem by either completely neglecting the gravity or assuming the particle inertia is very small or very large. In the context of cloud droplets, the air flow dissipation rate, droplet response time and droplet terminal velocity could each vary by two to three orders of magnitude, most of the existing theoretical results can at best be viewed to be either only qualitatively correct or applicable to some limiting cases for cloud droplets.
In section 3, we will first develop a theory for radial relative velocity of sedimenting droplets in turbulent air, using a statistical approach assuming that the differential sedimentation governs the relative motion of droplets before geometric collision contact. We will show that the resulting theory correctly reproduces several limiting cases. It will be also compared to DNS results reported in part 1 as well as other previous theoretical formulations.
In order to predict the collision kernel, we also need a model for the RDF. This will be constructed in section 4 by extending the theory of Chun et al [3] and by fitting DNS data of the RDF. Finally the theory for radial relative velocity and the model for RDF are combined to predict the DNS dynamic collision kernel in section 5. This integrated parameterization of the geometric collision kernel has already been used in our recent study of warm rain initiation [4].

Theoretical background
As in part 1, we focus here on geometric collision of sedimenting droplets without considering the local droplet-droplet aerodynamic interactions. For a dilute, bidisperse system containing droplets of radii a 1 and a 2 , the average geometric collision rate per unit volume,Ṅ c , can be expressed asṄ where n 1 and n 2 are average number concentrations of the two size groups and 12 is known as the collision kernel. In part 1, it was confirmed by DNS that the collision kernel for sedimenting droplets in a turbulent air can be expressed in terms of the kinematic pair statistics as K 12 = 2π R 2 |w r |(r = R) g 12 where R = a 1 + a 2 , the radial relative velocity w r is defined in terms of the relative velocity w between two droplets with separation vector r as w r ≡ w · r/|r|, with r ≡ |r|. The angular brackets denote an average over all possible directions of r in addition to the usual ensemble averaging or averaging over space and time for a homogeneous and stationary system. The average radial relative velocity at contact |w r |(r = R) represents the average relative flux on the geometric collision sphere, a concept introduced initially by Saffman and Turner [5]. The other factor g 12 (r ) in equation (2) is known as the RDF measuring the effect of preferential concentration on the pair number density at separation distance r . It is the ratio of the actual pair concentration to that in a uniform suspension. For the base case of sedimenting droplets in stagnant air, the kinematic properties and the collision kernel are given as |w r | g = 1 2 |v p2 − v p1 |, g g 12 = 1.0, where v p1 and v p2 are the terminal velocities corresponding to the two different sizes. 4 DNS results [2] show that air turbulence can enhance the collision kernel by increasing both the relative velocity and the RDF, relative to their respective value for gravitational collision. The overall enhancement can be expressed as where |w r | / |w r | g represents the enhancement by air turbulence on the relative pair motion and will be referred to as the turbulent transport effect, g 12 represents the effect of local nonuniform pair concentration due to interaction of droplets with local flow structures, and this contribution will be referred to as the accumulation effect.
The task of predicting the turbulent collision kernel can then be divided into two parts: a theory for the turbulent transport effect and a formulation for the accumulation effect. Each of these will be considered in this paper. We will first review previous studies on these two aspects for both nonsedimenting and sedimenting particles, before proposing our own developments. Most of the theoretical studies discussed below were published outside the atmospheric sciences literature.

Radial relative velocity
The relative velocity between two colliding droplets is affected by the unsteady, threedimensional air velocity field and is usually larger than the differential terminal velocity in still air. Theoretical modeling of this enhanced relative motion has been attempted for over half a century. Most models were developed assuming no particle clustering (g 12 = 1.0), using either the cylindrical formulation or the spherical formulation [6,7]. For the cylindrical formulation the relative velocity was derived using where p(w) is the probability density of the vector relative velocity w when the pair is separated by a distance R. For the spherical formulation, the relative velocity was obtained in terms of the radial component where p(w r ) is the probability density of w r . Most studies assumed that p(w) and p(w r ) follow the Gaussian probability distribution. Furthermore, p(w) and p(w r ) were typically assumed to be independent of the direction of R. This is a reasonable approximation for p(w r ) in the absence of gravity, however, Wang et al [7] showed that this assumptions does not hold for p(w) and, otherwise, an overprediction may occur if the dependence of p(w) on the orientation of R is not correctly treated.
In the limiting case of zero-inertia particles, Saffman and Turner [5] derived an analytical expression for the radial relative velocity using a spherical formulation where is the average dissipation rate per unit mass and ν is the air kinematic viscosity. This simply represents the contribution from local flow shear at the dissipation scale. In the same 5 paper, using a cylindrical formulation, they simultaneously considered the relative motion of droplets due to the dissipation range fluid shear, fluid acceleration, and gravity and obtained the following expression: where τ p is the droplet inertial response time, u is the local air flow velocity, the first term is due to the acceleration mechanism, and the second is the gravitational term, and the third is due to the shear mechanism. It was assumed in the derivation of the above expression that τ p is much less than the flow Kolmogorov time τ k . Part 1 gives the typical values of τ k for cloud turbulence. Note that the shear term in equation (8) differs from the shear term in equation (7) due to an incorrect isotropy assumption for p(w), as shown in [7]. Furthermore, the above model does not properly account for the gravitational effect as the expression does not reduce to the expected result of |w| = |v p2 − v p1 | in the absence of turbulence. Hu and Mei [8] improved Saffman and Turner's equation based on an asymptotic analysis for small particle inertia in a turbulent flow. They derived an equation for |w r | which includes the shear mechanism, the acceleration mechanism, and a new coupling term which accounts for the combined effect of spatial variation of fluid acceleration and particle inertia where λ D is the longitudinal Taylor-type microscale of fluid acceleration. Wang et al [6] later extended Hu and Mei's result to suggest an expression including the gravity effect This is a more consistent expression than Saffman and Turner's as the gravity-only case is recovered when neglecting turbulent effects. If the gravity effect is neglected, equation (10) is identical to Hu and Mei's result. Using the spherical formulation, Dodin and Elperin [9] decomposed the relative velocity into turbulent and gravity-induced components and assumed that the turbulent component is normally distributed. They arrived at the following result: where The above result recovers the two special cases of inertialess particles without gravity and the gravity-only case without turbulence. They also showed that Wang et al's expression, equation (10), is in close agreement with their formulation.
In the limit of large particle inertia (τ p > T L with T L being the Lagrangian integral scale of air turbulence), particles move independent of one another, i.e. their individual motion is completely uncorrelated and similar to the random motion of molecules in the kinetic theory. Abrahamson [10] deduced the following result: where v (i) p denotes the particle rms fluctuation velocity which can be obtained from where u is the fluid rms fluctuation velocity. Abrahamson also generalized his theory to include gravity and obtained This expression is consistent with (15) when the terminal velocities are set to zero as the second term in the rhs of (17) becomes identical to the first term in (17). In the absence of air turbulence (i.e. v (i) = 0), this expression also recovers the expected result for the gravity-only case. The above studies primarily addressed the two limiting cases of particle inertia. For intermediate particle inertia (τ k < τ p < T L ), major theoretical difficulties arise in predicting the relative velocity. In this situation, it becomes necessary to consider the interaction of particles with the complete spectrum of turbulent eddies as well as to account for the correlation of the motion of neighboring particles.
Williams and Crane [11] proposed an expression of the radial relative velocity due to the acceleration mechanism, taking into account the correlation of particle velocities. Since analytical solutions were obtainable only for small and large particle response times, the relative velocity applicable to the whole range of particle inertia was interpolated with these limiting solutions. Williams and Crane considered the acceleration mechanism but neglected the shear mechanism. Their expression followed a spherical formulation and was written as with w 2 accel WC where θ i ≡ τ pi /T L . This interpolated result, however, fails to recover their own limiting case for weak-inertia particles. Yuu [12] took an approach similar to Williams and Crane's and derived a solution for intermediate particle inertia with the combined effects of the shear and acceleration collision mechanisms using a cylindrical formulation where . Yuu also included the effect of buoyancy through the added mass coefficient b = 3ρ/(2ρ w + ρ) which is essentially zero for our cloud droplet problem. Here ρ and ρ w are the air and water densities, respectively. The drawback of Yuu's model is that the velocities of equalinertia particles are perfectly correlated, which is a poor approximation for large particles [13]. Therefore, his model cannot be used to predict the collision rate of large identical particles. Furthermore, his result is not valid for very weak-inertia particles as the exponential form of the Lagrangian fluid correlation adopted in his derivation is not valid in the viscous subrange of turbulence [14].
Kruis and Kusters [14] derived expressions for the relative velocity for the limiting cases of weak-inertia and large-inertia particles using a cylindrical formulation. For the weak-inertia particle case, they followed Yuu's procedure but included the correct parabolic form of the correlation coefficient used by Williams and Crane [11]. Meanwhile for the large-inertia particle case, they extended Williams and Crane's procedure by including the added mass effect. By interpolating the limiting solutions, a general solution was obtained yielding where γ is defined as A drawback of their theory is that their expression for the shear mechanism vanishes when the buoyancy effect is neglected. Thus, in this sense, their theory is incomplete.

8
Based on DNS data, Wang et al [15] and Zhou et al [16] inferred that no existing model can accurately predict the particle relative velocity. Therefore, they modified the model by Kruis and Kusters [14] using the spherical formulation and empirical fitting to DNS results, yielding The parameter γ was also modified to γ = αγ KK . Note that they modeled the shear mechanism using Saffman and Turner's weak-inertia expression. Their model was shown to be in good agreement with their DNS data for |w r | at several moderate flow Reynolds numbers. Zaichik et al [13] developed two statistical models to predict the turbulence-induced collision rate for monodisperse systems (i.e. a 1 = a 2 ). These models are valid over the entire range of Stokes numbers from the zero-inertia limit to the high-inertia limit. The first analytical model is based on the assumption that the joint probability density function of the fluid and particle velocities follows a correlated Gaussian distribution and the result is as follows: f u is the coefficient measuring particle response to the fluid velocity fluctuations, f (R) is the spatial correlation coefficient of the fluid velocities at two points separated by the collision radius R, and γ is defined by equation (23), with an improved model for the acceleration variance as a 01 = 11, a 02 = 205, a 0∞ = 7.
Their model was shown to recover the limiting cases of Saffman and Turner [5] and Abrahamson [10].
The parameterization, equation (26), was derived from the DNS data of fluid Lagrangian acceleration at low flow Reynolds numbers and the experimental data at moderate Reynolds numbers (R λ ∼ 2000), see figure 1 of Zaichik et al [13]. Hill [17] proposed asymptotic scaling formulae for the Lagrangian fluid acceleration as follows: Du Dt where S is the skewness of longitudinal velocity gradient and its value for R λ < 20 was taken from figure 8 of Herring and Kerr [18]. In figure 1, we compare Zaichik et al's parameterization and Hill's parameterization for the Lagrangian fluid acceleration. A smooth interpolation (a polynomial fit matching both the value and the gradient on the linear-log plot at R λ = 20 and 400) is used for Hill's model in the range of 20 < R λ < 400 (the region marked by two short vertical lines). The two models are very similar for low to moderate Reynolds numbers. At very high flow Reynolds numbers, the two models show different dependence on R λ . At the typical cloud flow Reynolds number of R λ = 20 000, the relative difference between the two models is about 30%. Zaichik et al [13] also presented a second model that stems from a kinetic equation for the particle-pair probability density function (PDF) of the relative velocity distribution. It predicts the statistical behavior of a pair of particles, namely, their mean and turbulent relative velocities and the RDF. The differential model provides a better prediction of |w r | when compared to the DNS data of Wang et al [15], whereas their first model, equation (26), overpredicts the radial relative velocity in general due to an underprediction of the two-particle velocity correlation, because they modeled the fluid velocity correlation seen by finite-inertia particles simply as Table 1. A summary of previous theoretical studies of two-particle relative velocity in a turbulent flow.
In a recent study, Pinsky et al [19] discussed the effect of nonGaussian Lagrangian fluid acceleration on the relative motion of droplets. Their theory is in spirit similar to Wang et al [6], namely, a leading-order formulation in terms of the droplet inertial response time. They showed that the nonGaussian distribution can lead to typically 10% reduction of the radial relative velocity. Unfortunately, their results were not given in a closed-form expression. Table 1 summarizes the previous efforts in modeling the two-particle relative velocity in a turbulent flow. We conclude that there is a lack of effort in modeling the pair relative-velocity for sedimenting particles. The available formulations for sedimenting particles address only the limiting cases of very large and very small particle response times. For the cloud physics problem, the droplet response time could be of the order of the Kolmogorov timescale (τ k ), while at the same time the gravitational effect is crucial. Thus, for the problem of cloud droplets, there is no theory available to predict the pair relative velocity. In section 3, we shall develop a theoretical model for |w r | that accounts for both particle inertia and gravitational effects.

RDF
In local regions of the flow where the air streamlines are severely curved (e.g. regions of high vorticity or high strain rate), sedimenting droplets, as a result of their finite inertia, can be somewhat nonuniformly distributed. This is shown in figure 8 of the part 1 paper.
Here, we review previous attempts to model this effect in terms of the (RDF) (or g(r ) for monodisperse systems and g 12 (r ) for bidisperse systems). The previous efforts were mostly based on observations from DNS. Reade and Collins [21] determined the dependence of monodisperse RDF g(r ) on St for R λ = 54.5 in the following form: where and Wang et al [15] proposed an equation that relates g(r = R) to St = τ p /τ k for several R λ in their DNS (45 < R λ < 75) but with R = η, where τ k and η are the flow Kolmogorov time and length. Their model is as follows: where Zhou et al [16] proposed a formulation for bidisperse systems by extending the formulation of Wang et al [15]. They related the bidisperse RDF to monodisperse RDF as with an empirically obtained correlation coefficient ρ n The monodisperse RDFs, g 11 (R) and g 22 (R), are given by equation (31).
The factor ρ n 12 in equation (35) signifies the rapid decorrelation of the concentration fields of different size droplets, leading to rapid drop of the magnitude of the bidisperse RDF relative to the monodisperse RDF [16]. Such a rapid reduction of g 12 was also noted by Bec et al [20] who studied clustering and collisions of heavy particles of different sizes in a random flow.
The above equations were developed for a limited range of R λ . Collins and Keswani [22] investigated the scaling of preferential concentration with the Reynolds number. They obtained g(r ) for a range of R λ from 65 to 152 and found that the RDF appears to be approaching a plateau with increasing R λ for any interparticle distance r much less than η. This appears to contradict the model by Wang et al [15] who show a linear growth with Reynolds number for 45 < R λ < 75. These could be due to the very different nominal particle sizes used in these studies. It is worth pointing out that both Wang et al [15] and Collins and Keswani [22] used ghost particles in their simulations. The differential model of Zaichik et al [13] appears to be capable of predicting DNS results of the RDF reported in Reade and Collins [21] and Wang et al [15]. But, again, the differential model does not provide a closed-form solution for g(r ).
Chun et al [3], following a similar procedure developed by Balkovsky et al [23] and Zaichik and Alipchenkov [24], derived expressions for g(r ) and g 12 (r ) in the limit of small particle inertia and r < η. They obtained the following expression for g 12 (r ) where c * 0 is a constant that was not explicitly derived by the authors but they claim that c * 0 depends upon the manner in which the flow around particle 1 transitions in space from local shear-driven small-scale to large-scale turbulence. r c is a length scale of the turbulent acceleration diffusion experienced by the particles. S 2 1 and R 2 1 are the average second invariant of the rate of strain and rate of rotation tensors seen by particle 1. a o is the normalized acceleration variance in isotropic turbulence (a o ≡ (Du/Dt) 2 /(v k /τ k ) 2 ), τ A is a characteristic correlation time for acceleration, and B nl is a nonlocal diffusion flux coefficient which was found to be 0.0926.
This expression also covers the monodisperse case when r c = 0, then which is similar to the power law equation proposed by Reade and Collins [21]. Although this model provides a better understanding of the physics of particle clustering, no explicit expressions were developed for c * 0 , S 2 1 , R 2 1 and τ A , thus the model remains somewhat incomplete.
In summary, the RDF has been studied either empirically or through PDF modeling approaches that tend not to yield a closed-form expression, due to the complexities of the clustering phenomenon. Unfortunately, the gravitational effect has not been included in any of the models discussed above. Gravity may further decrease particle clustering as it reduces the interaction time of the particle with the turbulent eddies. This was indeed observed in DNS by Vaillancourt et al [25], Wang et al [2,26] and Franklin et al [27] for sedimenting droplets. The only theoretical studies of RDF for sedimenting droplets in a turbulent flow is the work of Falkovich et al [28]- [30]. By extending the formulation for monodisperse RDF, Falkovich et al proposed the following expression for bidisperse RDF: where v k is the flow Kolmogorov velocity. The scaling exponent α should in general be a function of droplet response times, gravity and flow Reynolds number. The inclusion of gravity in the above expression was based on a heuristic argument. It is not clear which droplet size should be used to compute the τ p in the last term. Furthermore, since no explicit, closed-form expression for the scaling exponent α for a bidisperse, sedimenting system was provided, it is not clear how to compute g 12 using this theory. A comparison of Falkovich et al's theory with our DNS data was shown in figure 9(b) of the part 1 paper [2], using a constant α based on the computed monodisperse RDF.
In section 4, we will develop empirical models for g and g 12 when the gravity effect is considered.

Radial relative velocity for sedimenting droplets
Here, we shall develop a theory for |w r | that is applicable to inertial droplets sedimenting under gravity in a turbulent flow. The basic assumption is that droplet relative trajectory is mostly determined by the gravitational sedimentation. The droplet inertia introduces a secondary correction to the relative trajectory. Therefore, the theory is applicable to cloud droplets where the terminal velocity, relative to the flow Kolmogorov velocity, is large. We will compare the theory with the DNS data in part 1 and other previous theories.

Theory
Following Dodin and Elperin [9], the radial relative velocity between two particles falling under gravity in a homogeneous isotropic turbulent flow at a specific angle of contact, θ , can be decomposed into a random part ξ caused by turbulent fluctuations and a deterministic part h due to gravity where the angle of contact, θ , is measured from the vertical axis z (see figure 2). The statistics of the radial relative velocity w r (θ ) are axisymmetric with respect to the z if the turbulent flow is isotropic and the gravity is aligned with the z-axis [9]. Assuming the Stokes drag acting on the droplets, h(θ ) can be written as The random variable ξ(θ ) has a mean equal to zero and a standard deviation denoted by σ (θ). For sedimenting droplets, ξ(θ) can depend on θ due to the coupling of turbulent fluctuations and gravity on the motion of droplets. The kinematic formulation of collision kernel requires the average magnitude of the radial relative velocity |w r | . Taking into consideration the axisymmetry and assuming a uniform distribution of pair density, an integral average over the surface of the collision sphere yields where |w r (θ )| is the mean absolute value of radial relative velocity at a given θ . Note Since the random variable (−ξ ) has the same distribution as ξ , the latter expression yields |w r (π − θ)| = |w r (θ )| as illustrated in figure 3. Therefore, equation (40) can be rewritten as where |w r (θ )| can be obtained through its probability distribution function here the mean value of |w r (θ )| does not depend on the sign of h, thus we choose to define h = |g||τ p1 − τ p2 |cos θ . Assuming that ξ is normally distributed, the probability distribution function P w r (θ ) can be written as Combining the above two equations, we have (44) in which the only unknown is the variance σ (θ) 2 defined as where v (i) r is the turbulent component of the ith droplet velocity along the line of centers. In terms of the Cartesian coordinates 15 note that due to the axisymmetry the y-component and the x-component of the turbulent relative velocity have identical statistics.
x ) 2 θ represent the variance of the z-component and the x-component, respectively, of the turbulent radial relative velocity at a given orientation angle θ .
The problem now reduces to determining the variance of the turbulent radial relative velocity. The terms (v (1) i − v (2) i ) 2 can be written as where (v (1) i ) 2 is the mean square of the ith velocity component (i = x, y, z) of the size-1 droplets (similar notation for size-2 droplets), while (v (1) i v (2) i ) is the cross-velocity covariance in the ith direction.
The turbulent component of the particle velocity can be found by integrating the particle equation of motion After separating out the mean motion due to gravity [31], we obtain v (k) Setting the reference time t, when the pair is at contact, to zero, we can express the particle velocity covariance and the mean-square particle velocity as v (1) i where is the two-time fluid velocity covariance seen by two particles located at Y (1) (τ 1 ) and Y (2) is the fluid velocity correlation seen by a single particle. Both covariances depend on the particle trajectories and therefore the gravitational settling, a manifestation of the coupling between turbulent fluctuations and gravity. Figure 4 illustrates the trajectories of two particles before colliding. It is important to note that due to the exponential terms in (51) and (52), the main contribution to variance and covariance of particle velocity comes from the time interval right before the contact, say −3τ pk < τ k < 0. Therefore, a good representation of the integrands is needed for that time interval. Now consider the two-particle fluid velocity correlation, during this period of time the droplet trajectory may be approximated by Y (2) (τ 2 ) ≈ X + V (2) (0)τ 2 and Y (1) (τ 1 ) ≈ X + R + V (1) (0)τ 1 , where X is the location of particle 2 at the moment of contact, R is the separation vector at contact, and V (1) (0) and V (2) (0) are the particle velocities at the moment of contact. As a leading-order approximation, we assume that the particle trajectories before collision (in the time interval −3τ pk < τ k < 0) follow straight lines with a velocity equal to the terminal velocity. The two-point, two-time fluid velocity correlation is approximated as

The fluid velocity correlations.
where the standard two-point fluid velocity correlation and the Lagrangian fluid velocity correlation are defined as where Y f (t) is a fluid Lagrangian trajectory. The approximation given by equation (53) becomes exact for two important limiting cases. Firstly, when the terminal velocities v p1 and v p2 are large, the trajectories are indeed vertical lines and as such the Eulerian two-point fluid correlation essentially determines the required Lagrangian fluid velocity correlation seen by the particle pair, regardless of the inertial response times of the particles. This similar concept was utilized in the study of turbulent particle dispersion in [33]. Secondly, when both the terminal velocities and the inertial response times are small, the particle pair essentially follows the motion of a fluid element, the fluid velocity correlation seen by the particle pair reduces to the Lagrangian fluid velocity correlation. In the context of cloud droplets, the inertial response time of the droplets is usually small, use of the Lagrangian time correlation is justified as a leading-order approximation. The spatial correlation also includes the local shear effect of turbulence on the relative motion of droplets [5]. In summary, the approximation (53) combines the effects of local shear, fluid Lagrangian time correlation and gravitational settling on the fluid velocity correlation relevant to the relative motion of droplets.
Similarly, the fluid velocity correlation seen by a single droplet is approximated as (56) For the fluid Lagrangian velocity correlation D L , we used the two-scale bi-exponential form [13,34] where T L is the Lagrangian integral timescale and z = τ T /T L . The Lagrangian integral timescale T L may be modeled as T L ≈ u 2 / . For 2z 2 > 1, the above expression may be written in terms of trigonometric functions as The Lagrangian Taylor microscale time τ T is defined as where a o ≡ (Du/Dt) 2 /(v k /τ k ) 2 is the nondimensional variance of fluid acceleration, and is modeled according to equation (26). The expression for D L satisfies the conditions when τ τ T and (τ T /T L ) → 0 (for high R λ flows). On the other hand, for small time delay (τ/τ T 1), the Lagrangian correlation becomes which is consistent with the definition of τ T in equation (59). For isotropic and homogeneous turbulence, the Eulerian two-point velocity correlation R i j can be described through the longitudinal and transverse two-point velocity correlations ( f (r ) and g(r )) as follows: The corresponding two-scale bi-exponential form for f (r ) is here, L e ≡ ∞ 0 f (r ) dr is the Eulerian integral length scale and β = √ 2λ/L e . For 2β 2 > 1, the equation may be rearranged to obtain an expression similar to (58). This expression recovers the limiting cases as It is important to note that, for collision of cloud droplets considered in this paper, the inertialrange dynamics is not directly relevant, but rather the dissipation-range flow dynamics is most relevant. One way to see this is to compare the relaxation distance τ 2 p g with the Taylor microscale λ. The ratio τ 2 p g/λ is approximately equal to 0.52St S v R −0.5 λ . Taking a typical cloud Reynolds number of R λ = 10 000 and dissipation rate of 400 cm 2 s −3 , this ratio is between 0.000 15 and 0.15. Therefore, the inertial-range scaling is not relevant. The general expression, equation (63), has the correct dissipation-range scaling as shown by equation (65). It was also designed to produce the correct integral length scale. The transverse correlation function g(r ) is related to f (r ) as [35] 3.1.2. The variance of particle relative-velocity fluctuation. We shall now derive a theory for the variance (σ 2 ) of turbulent component of the particle relative velocity. Equations (51) and (52) can be written as These properties are needed in order to evaluate σ (θ) through equations (46) and (47). Let us first consider the special case that the particles collide with θ = 90 • as shown in figure 5. In this case, equation (46) reduces to where R 11 (v pk e z τ ) is evaluated by using g(v pk τ ). In order to simplify the integrand, the two-particle fluid velocity correlation R 11 (R + v p1 e z τ 1 − v p2 e z τ 2 ) is approximated by f (R)g(v p1 τ 1 − v p2 τ 2 ). This turns out to be an excellent approximation as shown in figure 6, since, for our application, R η and f (R) = 1. Therefore, equations (70) and (71) can be written as Substituting the expressions for D L (τ ), f (r ) and g(r ), and after lengthy integrations the particle velocity correlation (v (1) x Figure 6. Two-point fluid velocity correlation along the trajectories of two particles colliding side-by-side at θ = 90 • . z is the vertical distance between the two particles before collision. The distance is given in cm and R = 60 µm.
where the constants b 1 , b 2 , c 1 , c 2 , d 1 , d 2 , e 1 and e 2 are defined as where it is assumed that v p1 > v p2 . Similarly, the variance of the kth particle velocity is given by where (α, φ) is defined as A similar procedure can be carried out to obtain the variance of ξ for θ = 0 • and θ = 45 • using equation (46) (see [36]). Figure 7 shows σ (θ) for three θ values as a function of a 2 with a 1 = 30 µm using the derived expressions. The turbulent flow used for the figure has R λ = 72.41 and = 400 cm 2 s −3 . The flow parameters T L , τ T , L e and λ were obtained from DNS, however, they could also be estimated in terms of and R λ as we will explain later. σ for the case of nonsedimenting droplets is shown, in which case σ is independent of θ. It appears that gravity enhances the fluctuation level of w r when compared to that of the nonsedimenting case.
If we take the case of nonsedimenting particles as the reference, the addition of gravity introduces two effects. Firstly, gravity reduces the interaction time of droplets with turbulent eddies, and therefore the variance of particle velocities (i.e. the first two terms in equation (47)) will be reduced ( figure 8). Secondly, gravity also decreases the correlation coefficient as shown in figure 9. This reduces the last term in (47). The net effect here is an increase of fluctuation of w r . Figure 7 also shows that σ is slightly larger for larger θ . Since this θ -dependence is relatively weak and the differential terminal velocity dominates the relative velocity for our application (figure 7), we shall no longer explicitly consider the θ -dependence. The result for θ = 90 • will be applied to all θ in deriving a closed-form expression for |w r | .

The final result for |w r |
Using σ (θ = 90 • ) to approximate σ (θ), equation (44) can be readily integrated to obtain where σ has been expressed in terms of τ p1 , v p1 , τ p2 and v p2 and the flow parameters T L , τ T , L e and λ are defined in the previous section. The above expression is identical in form to that of Dodin and Elperin [9]. However, our approach leads to a more general theory for σ that is applicable to arbitrary τ p1 and τ p2 with gravity. Dodin and Elperin [9] assumed τ p1 τ k and τ p2 τ k in their derivation of σ .
We shall now examine the limiting cases with our theory. For weak-inertia particles, using the fluid Lagrangian velocity correlation for short time delay D L = 1 − (τ/τ T ) 2 , we obtain (81) Figure 8. Variance of the particle velocity as a function of non-dimensional particle terminal velocity. The particle Stokes number is 0.5 and it is falling in a turbulent flow with R λ = 72.41 and = 400 cm 2 s −3 .
The first two terms represent the shear effect and the acceleration effect, similar to the results of Saffman and Turner [5]. The third term is new and it represents a coupling between gravity and fluid shear, a term that is absent in the previous theory by Dodin and Elperin [9]. This new term augments the turbulent transport effect. For typical conditions in clouds (e.g. R λ = 2 × 10 4 and = 100 cm 2 s −3 ), the ratio of the acceleration term to the shear term can be written as For a 10-20 µm or 20-30 µm pair falling in a turbulent flow with R λ = 2 × 10 4 and = 100 cm 2 s −3 , the value of this ratio is 1.89 × 10 3 , which means that for bidisperse systems of droplets, the turbulent acceleration effect is three orders of magnitude larger than the turbulent shear effect. The value of a o was obtained by using equation (26). The ratio of the coupling term to the shear term is This ratio is independent of , but increases quickly with the size of droplets. For instance, the ratio for a 10-20 µm pair in the above turbulent flow is 140.49. For a 20-30 µm pair, it becomes 949.75. The new coupling term is shown to be also important and cannot be neglected. The second limiting case is for very large particle inertia, namely τ p1 T e ≡ L e /u , τ p2 T e . This case was studied by Abrahamson [10]. In this limit, it is expected that two particles move independently or the last term in (47) approaches zero. The particle relative velocity is then determined entirely by the particle velocity fluctuations as given by the first two terms in (47). Our theory can be applied to reveal the nature by which this limiting case is approached. Figure 10 shows that both the variance of particle fluctuating velocity and the covariance decrease with increasing τ p2 /T e , but the magnitude of the covariance is much smaller than the variance recovering the qualitative picture illustrated in Abrahamson's theory. It is important to point out that Abrahamson [10] did not consider the coupling between gravity and turbulence, so the variance of particle velocity was overestimated in his theory. The third limiting case concerns very large terminal velocities (i.e. very large S v ). Under this limit, both (v (1) v (2) ) and (v (k) ) 2 tend to zero (because (α, φ) and (α, φ) go to zero) and our theory yields the expected radial relative velocity for gravitational collision. Physically, when particles settle very fast through the turbulent flow, they do not have time to interact with the flow so that the particle motion is unaffected by the flow.
Finally, for the case of negligible gravitational settling, we compare our theory with the DNS data from Wang et al [15] and Zhou et al [16] in figure 11. Our theory was not intended for this limiting case, and underpredicts the DNS results, particularly for the monodisperse system. This results from an overestimation of the covariance in the theory.
Coming back to cloud droplets, figure 12 demonstrates that the variance of droplet velocity is well predicted by our theory. Cloud droplets correspond to small particle inertia and large Figure 10. Comparison of (v (k) ) 2 and (v (1) v (2) ) predicted by the present study and by Abrahamson's theory for particles falling in a turbulent flow of R λ = 72.41 and = 400 cm 2 s −3 . τ p1 = 2.5 T e , Sv 1 = 9.3, and Sv 2 = 1.25. The solid line represents the two-particle cross velocity correlation for our model, while the other lines represent the mean square velocities of the particles.
settling for which the approximation used in our theory is expected to work best. When the terminal velocity is set to zero, our theory for (v (k) ) 2 is identical to that of Zaichik et al's [13] first approach. Zaichik et al [13] also found that the particle-velocity variance for nonsedimenting particles from DNS of Sundaram and Collins [37] can be well predicted. The comparison in figure 12 shows that our extension to include the gravity effect is very successful for single-particle energetics.

Comparison with DNS results
In figure 13, we compare DNS results in part 1 with our theory and previous theoretical expressions for |w r | that simultaneously considered turbulence and gravity. Figure 13(a) shows that our theory agrees with DNS data to within 5% relative difference for bidisperse cases. On the other hand, the relative error in the predicted |w r | for monodisperse systems at = 400 cm 2 s −3 is about 48%. As noted earlier, the approximation for pair trajectories used in our theory is not designed for droplet pairs of similar size (i.e. very small differential sedimentation). Interestingly, even for this case, figure 13(b) shows that the theory provides a qualitatively correct prediction. It should be noted that the weak θ -dependence discussed in figure 7 can be included in our theory, but the case of weak sedimentation is not explicitly considered.
The theory by Abrahamson gravely overestimates the radial relative velocity due to the exclusion of the covariance term and also overestimation of the variance. The weak-inertia Figure 11. Comparison of the predicted and simulated radial relative velocities of nonsedimenting particles in a turbulent flow. The particle radius is equal to half the Kolmogorov length scale (a = η/2). (a) Bidisperse system (similar to figure 12 in Zhou et al [16]). (b) Monodisperse system (similar to figure 8 in Wang et al [15]).  theory by Wang et al [6] underestimates |w r | for bidisperse systems but it overestimates |w r | for monodisperse systems.
We also compare our theory with several other theories for sedimenting droplets in figure 14(a) for droplet pair of 15 and 10 µm in radii, at a Reynolds number of R λ = 20 000. This plot is similar to figure 8(a) in Pinsky et al [19] and Hill's fluid acceleration model is assumed. All these previous theories shown in figure 14 assumed that the two droplets interact with the same local turbulent eddy, and the lifetime of such interaction is not affected by sedimentation. In other words, both droplets follow precisely the same Lagrangian fluid element during the whole interval of the interaction. Whereas in our theory, the relative position of the pair is continuously changing due to differential sedimentation, thus the local fluid velocities seen by the two droplets are different and are not perfectly correlated. Therefore, our theory considers the coupling between the eddy-droplet pair interaction and differential sedimentation. As a result, our theory predicts a smaller enhancement of the relative motion by air turbulence. In figure 14(b), we plot the Stokes numbers and the ratios of droplet terminal velocity to the Kolmogorov velocity (S v = v p /v k ). Since S v is at least one order of magnitude larger than St and St is finite for larger droplets, our treatment here could be more appropriate than the leadingorder treatment based on small St numbers.
The dependences of the radial relative velocity on flow Reynolds number and fluidacceleration parameterization are shown in figure 15. Our theory shows a much weaker dependence on flow Reynolds number than the results of Pinsky et al [19], for reasons discussed above. Hill's model of fluid acceleration yields a slightly more sensitive dependence on flow Reynolds number than Zaichik et al's model.   [19] are also shown for comparison. The flow Reynolds number is assumed to be R λ = 20 000.

RDF: the accumulation effect
In this section, we examine the second effect of air turbulence, namely, how the RDF depends on gravity and the turbulent flow parameters.
Several studies have attempted to correlate g(r ) with particle and flow parameters for a monodisperse suspension of droplets in a turbulent flow. Typically, a power law expression is used. For nonsedimenting particles, Chun et al [3] developed an expression for g(r ) as follows: This expression is only valid for r smaller than the Kolmogorov length scale (η), which is relevant for cloud droplets. We shall extend this form to sedimenting droplets by assuming that C 0 and C 1 are functions of St, R λ and a nondimensional parameter for gravity, |g|/(v k /τ k ). Using our DNS results, we obtained an empirical expression for C 1 by curve-fitting our DNS data with C 0 = 1. The resulting form is as follows: where In figure 16(a), we compare the empirical formula for g(r ) with DNS results for monodisperse particles in a turbulent air flow. It is observed that the level of particle clustering depends on the turbulent dissipation rate and R λ . The Stokes number increases with both and droplet radius a (see figure 1 in part 1). The largest clustering is developed by droplets whose St is about unity as it has been found in nonsedimenting cases [15]. We also observe a higher droplet clustering for higher R λ . The proposed empirical model correctly recovers several limiting cases. (i) For weakinertia particles (i.e. St 1), the model gives g(r ) = 1.0 as C 1 goes to zero. (ii) If the effect of gravity is very strong in comparison with the turbulence effect (i.e. |g| (v k /τ k )), C 1 also goes to zero which brings g(r ) = 1.0. (iii) The model mimics the behavior independent of R λ for small St (e.g. [15]). (iv) For very large R λ , the model implies that C 1 does not depend on gravity.
When a pair is made of two different sizes, any fluid acceleration or gravity will induce a relative velocity. This effect gives rise to a diffusion-like process in the system that tends to smooth out inhomogeneities in the particle pair concentration. Such a diffusion process contributes negatively to the pair clustering process, i.e. RDF will be reduced.
Chun et al [3] derived an expression for bidisperse pair RDF g 12 (r ) = C 0 η 2 + r c 2 r 2 + r c 2 where C 0 and C 1 are the same constants as in equation (84) for the monodisperse case. r c is a parameter that incorporates the diffusion process (see equation (36)). If r r c , the acceleration-driven diffusion predominates so that g 12 (r ) is independent of r . On the other hand, if r r c , the monodisperse case is recovered, i.e. droplets in the pair do not react differently to the local acceleration field. Extending the empirical model for the monodisperse system, we propose the following expression: where r η is implied. C 1 is assumed to follow the same expression found for the monodisperse case at St = max(St 2 , St 1 ). In order to include the effect of gravity, we modified the expression for r c in [3] as where a o g is and F(a o g , R λ ) is obtained by curve fitting No attempt was made to justify the functional forms here, except that a reasonable fit to our DNS data was realized. In figure 16(b), we compare the model with DNS results. The effect of acceleration-driven diffusion is very strong for bidisperse pairs, making g 12 (R) close to one if the size ratio a 1 /a 2 is larger than 1.3, even for figure 16(b), where the flow dissipation rate used is high in the context of clouds. The accumulation effect is important only when a 1 /a 2 < 1.2. In figure 16(b), we also include the model proposed by Zhou et al [16] for bidisperse systems, showing that their model overpredicts the DNS data as gravity was not considered in that study.
In figure 17, we compare the proposed model with previous models by Wang et al [15] and Reade and Collins [21] for a monodisperse system. The turbulent flow has R λ = 54.5 and = 400 cm 2 s −3 . We chose this R λ value because the model by Reade and Collins was developed based on this particular R λ .
The model by Wang et al [15] always underpredicts g(r ) as it was developed with R/η = 1. Since the preferential concentration decreases with separation distance [21], the underprediction by Wang et al's model [15] is anticipated. However, the overall shape of the present model is similar to that of Wang et al's.
The model by Reade and Collins [21] appears to overpredict g(r ) in the range from 20 µm to about 40 µm droplet radius. Since Reade and Collins' model was obtained for nonsedimenting particles, the overprediction is expected as gravity tends to drive particles out of vortical regions and decrease the tendency for particle clustering. The model by Wang et al [15] was also obtained for the case of no gravity. If this model were corrected for r < η, it would probably show a similar trend to Reade and Collins [21].

Geometric collision kernel
We shall now combine the theory for the relative velocity and the empirical model for RDF to predict the geometric collision kernel using the kinematic formulation for 12 , equation (2). In figure 18, we compare the geometric collision kernel predicted by the integrated model with the dynamic collision kernel from DNS. In general, the model provides a reasonable quantitative prediction. When the droplet Stokes number becomes large and the gravity effect is weak, namely for high flow dissipation rate ( = 400 cm 2 s −3 ), the model underpredicts the kernel by as much as 30%.
A known limitation of DNS is that the flow Reynolds number is much less than that in real clouds. We can use the integrated model to extend the range of R λ and to illustrate the enhancements by turbulence in real clouds. In figure 19, we show the results from the model at different R λ as a function of . Equation (26) was used to represent the turbulent acceleration variance ( (Du/Dt) 2 ), following Zaichik et al [13].
The model suggests that the contribution of R λ to the increase of collision rate is as important as the contribution of for the bidisperse collision between droplets of 10 and 20 µm in radii. The dependence on R λ is much weaker for monodisperse self-collisions of 20 µm droplets. While our model underestimates of bidisperse collisions for large , turbulence is expected to increase the overall collision kernel by a factor ranging from 2 to 3, relative to the gravitational kernel, when the flow dissipation rate is above 300 cm 2 s −3 and R λ is of the order of 10 4 . For 20 µm-droplet self-collisions, the collision kernel depends sensitively on . For example, the monodisperse collision kernel at = 600 cm 2 s −3 is about nine times the collision kernel at = 100 cm 2 s −3 . Typical dissipation rates in cumulus clouds range from 100 to 650 cm 2 s −3 [38,39]. Most of the dependence on flow Reynolds number results from its effect on RDF for the case shown in figure 19(a). This is demonstrated in figure 20, where the relative enhancements on RDF and the radial relative velocity are compared.
The theory developed above assumed a linear Stokes drag acting on the droplets. In part 1, it was shown that a nonlinear drag (NLD) must be used for droplets larger than 30 µm in radius, as the droplet Reynolds number is of the order of one. Here, we demonstrate that the effect of NLD can be included in the theory.
Following Wang and Maxey [32], we can apply our model to treat the nonlinear drag by modifying the particle response time and the terminal velocity as where the nonlinear factor f (Re p ) models the departure from the Stokes drag law and is given as f (Re p ) = 1 + 0.15Re 0.687 p [40]. Using these effective inertial response time and terminal velocity, the integrated model shows good agreement with the DNS results based on the NLD ( figure 21). For the bidisperse cases (figure 21), the droplet clustering is not significantly affected by the drag nonlinearity. The radial relative velocity is clearly decreased due to the reduced mean droplet terminal velocities. As a result, the collision kernel is also reduced by the NLD.

Summary and concluding remarks
A new theory has been developed to describe the radial relative velocity between a pair of droplets sedimenting under gravity in a turbulent flow. The theory is based on a statistical approach of the single-droplet and two-droplet Lagrangian velocity correlations. The main Figure 20. The relative enhancements of the radial relative velocity and the RDF by turbulence, as a function of R λ for droplet pair of 20 and 10 µm in radii (also see figure 19(a)). The dissipation rate is = 400 cm 2 s −3 . The horizontal dotted line marks the value of one. assumption is that the separation distance of a droplet pair before collision contact is governed by the differential sedimentation. This is a reasonable assumption for cloud droplets as the sedimentation dominates the droplet relative motion, with the air turbulence introducing a moderate modification to the relative motion. The theory yielded a reasonable quantitative prediction of the radial relative velocity of bidisperse pairs of different sizes when compared to the DNS data from part 1 [2]. For pairs of droplets of nearly equal size, the assumption is no longer valid; but interestingly, the theory could still provide a good qualitative prediction (see figure 13). In the weak-inertia limit, the theory reveals a new term making a positive contribution to the radial relative velocity resulting from a coupling between sedimentation and air turbulence on the motion of finite-inertia droplets. For droplets larger than 30 µm in radius, a NLD can also be included in the theory in terms of an effective inertial response time and an effective terminal velocity. This theory appears to fill a gap in the open literature concerning the prediction of radial relative velocity of finite-inertia, sedimenting droplets in a turbulent flow.
The relative sedimentation of droplets causes decorrelation of local fluid velocities and local fluid accelerations seen by the two droplets of different sizes. The explicit consideration of this aspect in our theory yields a weaker enhancement of the radial relative velocity by turbulence, when compared to previous results of Saffman and Turner [5], Wang et al [6], and Pinsky et al [19]. In these previous theories, a perfectly correlated local fluid element was used to describe the motion of two droplets. Our theory, which is supported by DNS data, implies that the relatively rapid sedimentation of droplets could reduce the applicability of approaches based on leading-order approximation in terms of inertial response time of droplets alone. Since our theory tends to underestimate the radial relative velocity for pairs of negligible differential sedimentation, it is likely that the true level of enhancement in radial relative velocity by turbulence lies between our theory and those of previous results.
In a few recent studies, the relative motion between droplets due to turbulence has been sub-divided into two general contributions: the relative motion due to local flow shear/fluid acceleration and the relative motion due to interaction of two droplets with different eddies before coming into contact that has been memorized by moderate droplet inertia. The latter is referred to as the sling effect by Falkovich et al [28,30]. It has been speculated that this sling effect could be significant compared to the local effects for large Reynolds-number turbulent air flow in clouds, leading to the formation of caustics of the particle velocity field and consequently, significantly enhanced collision rates [43]- [45]. The nonlocal effect is partially accounted for in our leading-order theory, at least for the case of large droplet sedimentation. This is because the integrals in equations (51) and (52) are over all previous time history. It may not be fully accounted for due to the leading-order approximation of the relative trajectory based solely on the relative sedimentation velocity. Since the integrals are over all time history, the local and nonlocal effects are not deliberately separated out.
An empirical model is also developed to quantify the RDF, by extending the theory of Chun et al [3] and by fitting DNS data of the RDF. This, together with the theory for the radial relative velocity, provides a theoretical parameterization for the turbulent geometric collision kernel. Using this integrated model, we find that turbulence could triple the geometric collision kernel, relative to the stagnant air case, for a droplet pair of 10 and 20 µm sedimenting through a cumulus cloud at R λ = 2 × 10 4 and = 600 cm 2 s −3 . For the selfcollisions of 20 µm droplets, the collision kernel depends sensitively on the flow dissipation rate. This parameterization of the geometric collision kernel has been used in our recent study of warm rain initiation [4], showing that the air turbulence plays a significant role in the initial phase of growth by collision-coalescence, namely, growth by autoconversion of cloud droplets.
The theory for the radial relative velocity extends the DNS results to higher flow Reynolds numbers and dissipation rates. However, it should be viewed as a first step due to the assumption used in its derivation. The approach used to approximate the pair separation before collision is similar to the first-order iteration in Reeks' approach for single particle dispersion [31]. Reeks [31] was able to derive the dispersion statistics using the successive second-order trajectory approximation. It may also be possible to develop a theory for the radial relative velocity using the successive second-order pair trajectory, but it would be a lot more laborintensive than what was done in this paper.
The modeling of the RDF remains a challenge. The recent advances made by Zaichik et al [13,41,42] based on a PDF theory offer a possibility for predicting RDF for finiteinertia sedimenting droplets, although so far they have not included the effect of gravitational sedimentation in their formulation. The gravity breaks the spherical symmetry of relative motion and pair concentration, so it is not yet clear whether their PDF approach can be applied to finite-inertia, sedimenting droplets. Certainly, more DNS results of both RDF and radial relative velocity at higher flow Reynolds numbers would be desired to build up sufficient data to test future theoretical developments. Specifically, the scaling exponent α for a bidisperse, sedimenting system has not been studied in DNS, which should be a priority for DNS studies.
Finally, we summarize the integrated parameterization of the turbulent collision kernel below. All relationships are compiled together so our proposed parameterization can be easily implemented by others in the future to include the effect of air turbulence when modeling droplet size distribution. For cloud droplets which satisfy the condition a η, ρ w ρ and Sv > 1, the geometric collision kernel can be calculated as  In the above parameterization, the input parameters are: (i) for the droplets, the radii a 1 and a 2 , and the water density ρ w ; (ii) for the turbulent air flow, the density ρ, the viscosity ν, the turbulence dissipation rate , the Taylor-microscale Reynolds number R λ , and (iii) the gravitational acceleration |g|. All other derived variables used in the model are listed in table 2.