Reynolds-number dependence of turbulence enhancement on collision growth

This study investigates the Reynolds-number dependence of turbulence enhancement on the collision growth of cloud droplets. The Onishi turbulent coagulation kernel proposed in Onishi et al. (2015) is updated by using the direct numerical simulation (DNS) results for the Taylormicroscale-based Reynolds number (Reλ) up to 1140. The DNS results for particles with a small Stokes number (St) show a consistent Reynolds-number dependence of the socalled clustering effect with the locality theory proposed by Onishi et al. (2015). It is confirmed that the present Onishi kernel is more robust for a wider St range and has better agreement with the Reynolds-number dependence shown by the DNS results. The present Onishi kernel is then compared with the Ayala–Wang kernel (Ayala et al., 2008a; Wang et al., 2008). At low and moderate Reynolds numbers, both kernels show similar values except for r2 ∼ r1, for which the Ayala– Wang kernel shows much larger values due to its large turbulence enhancement on collision efficiency. A large difference is observed for the Reynolds-number dependences between the two kernels. The Ayala–Wang kernel increases for the autoconversion region (r1, r2 < 40 μm) and for the accretion region (r1 < 40 and r2 > 40 μm; r1 > 40 and r2 < 40 μm) as Reλ increases. In contrast, the Onishi kernel decreases for the autoconversion region and increases for the rain–rain self-collection region (r1, r2 > 40 μm). Stochastic collision– coalescence equation (SCE) simulations are also conducted to investigate the turbulence enhancement on particle size evolutions. The SCE with the Ayala–Wang kernel (SCEAyala) and that with the present Onishi kernel (SCE-Onishi) are compared with results from the Lagrangian Cloud Simulator (LCS; Onishi et al., 2015), which tracks individual particle motions and size evolutions in homogeneous isotropic turbulence. The SCE-Ayala and SCE-Onishi kernels show consistent results with the LCS results for small Reλ. The two SCE simulations, however, show different Reynolds-number dependences, indicating possible large differences in atmospheric turbulent clouds with large Reλ.


Introduction
Several mechanisms have been proposed to explain the rapid growth of cloud droplets, which often result in fast rain initiation in the early stages of cloud development.Examples of these mechanisms include the turbulence-enhanced collision rate of cloud droplets (Falkovich and Pumir, 2007;Grabowski and Wang, 2013), turbulent entrainment (Blyth, 1993;Krueger et al., 1997), giant cloud condensation nuclei (Yin et al., 2000;Van Den Heever and Cotton, 2007), and turbulent dispersions of cloud droplets (Sidin et al., 2009).The first mechanism, which has received the most attention, has led to extensive research on particle collisions in turbulence (e.g., Sundaram and Collins, 1997;Wang et al., 2000;Saw et al., 2008;Onishi et al., 2009;Dallas and Vassilicos, 2011).
One direction taken by the research in this area is the simulation of collisional growth by solving the stochastic collision-coalescence equation (SCE).Such research relies on accurate collision-coalescence models, which consist of models for the collision kernel K c (r 1 , r 2 ) (where r i is the particle radius), the collision efficiency E c (r 1 , r 2 ), and the coalescence efficiency E coal (r 1 , r 2 ).To consider the influence of turbulence, several turbulent collision models have Published by Copernicus Publications on behalf of the European Geosciences Union.R. Onishi and A. Seifert: Re dependence of turbulence enhancement been proposed.Saffman and Turner (1956) analytically derived a collision kernel model for particles with no inertia or with a very small Stokes number (St = τ p /τ η , where τ p is the particle relaxation time and τ η is the Kolmogorov time), while Abrahamson (1975) derived a model for St 1.For moderate Stokes numbers, i.e., St ∼ 1, one difficulty is the preferential motion of inertial particles.Inertial particles preferentially cluster in regions of low vorticity and high strain if St 1 (Maxey, 1987), and they cluster in a way that mimics the clustering of zero-acceleration points by the sweep-stick mechanism if 1 St τ p /T I , where T I is the integral timescale of the turbulence (Coleman and Vassilicos, 2009).This matters because clustering increases the mean collision rate (Sundaram and Collins, 1997).To quantify the clustering due to the preferential concentration effect, a model is formulated for finite-inertial particles.However, the model requires several empirical parameters that should be determined from reference data, e.g., results from a direct numerical simulation (DNS).
One serious problem is that the Reynolds-number dependence of turbulent collisions has not yet been clarified.In fact, many authors ignore the Reynolds-number dependence and assume a constant collision kernel regardless of the Reynolds number (e.g., Saffman and Turner, 1956;Derevyanko et al., 2008;Zaichik and Alipchenkov, 2009) or assume a convergence to a constant collision kernel with increasing Reynolds number (e.g., Ayala et al., 2008a).Onishi et al. (2013) observed that the clustering effect, and consequently the collision kernel, decreases as the Taylor-microscale-based Reynolds number (Re λ ) increases for St = 0.4.Onishi and Vassilicos (2014) later clarified that the Reynolds-number dependence of the clustering effect for 1/3 St 1 is due to internal intermittency of the turbulence.Because a robust theoretical model for turbulent collision kernels is not yet available, we need empirical models for the investigation of turbulence enhancement on cloud development.As an example, the Ayala-Wang kernel (Ayala et al., 2008a;Wang et al., 2008) is a widely used turbulent kernel model.
Recently, Onishi et al. (2015) proposed an empirical kernel model based on DNS data for the wide range of 49 ≤ Re λ ≤ 530, where Re λ is the Taylor-microscale-based Reynolds number.Onishi et al. (2015) also conducted stochastic and direct collision simulations to investigate the turbulence enhancement on drop size evolution.They investigated the energy dissipation ( ) dependence for the range of 100 ≤ ≤ 1000 cm 2 s −3 and the Re λ dependence for the range of 66 ≤ Re λ ≤ 206.The results showed good agreement of the dependence between the stochastic simulations with the Ayala-Wang and Onishi kernels, but a significant discrepancy for the Re λ dependence between the two kernels.The discrepancy in Re λ dependence may become a critical issue for cloud simulations because Re λ is typically as large as O(10 3−4 ) in atmospheric turbulent clouds.However, Onishi et al. (2015) did not provide a detailed discussion on the difference of the Ayala-Wang and Onishi kernels in Re λ dependence.
This study, therefore, aims to compare the Ayala-Wang and Onishi kernels by focusing on their Re λ dependence.First, the Onishi kernel is updated by using the reference collision statistics obtained by the DNS for Re λ up to 1140.The Ayala-Wang and the present Onishi kernel values are compared in detail.The SCE simulations with the Ayala-Wang and Onishi kernels are also compared with each other and with the reference results from the Lagrangian Cloud Simulator (LCS; Onishi et al., 2015), which tracks individual particle motions and size evolutions in homogeneous isotropic turbulence.The collision growth simulation with the LCS is conducted for Re λ up to 333.
2 Turbulent coagulation kernel models

Turbulent coagulation kernel
The geometric collision frequency per unit volume between particles with radius r 1 and those with radius r 2 , N c (r 1 , r 2 ), is expressed by the geometric collision kernel K c (r 1 , r 2 ) as where n p,i is the number density of particles with radius r i .
The coagulation kernel K coag can be expressed by the combination of the geometric collision kernel, collision efficiency E c , and coalescence efficiency E coal as The gravitational collision kernel describes the collisions due to the settling velocity difference in the form of where R 12 (= r 1 + r 2 ) is the collision radius and V ∞i is the gravitational particle settling velocity.Turbulence enlarges the geometric collision kernel, i.e., the turbulent geometric kernel K c,turb is larger than K c,grav .Turbulence also enhances the coagulation kernel through enlarging E c .The turbulence enhancement on the collision efficiency, η E , is defined as where [T] and [NoT] indicate the turbulent flow case and the stagnant (non-turbulent) flow case, respectively.It had been difficult to confidently discuss the collision efficiency in a turbulent flow until Ayala et al. (2007) developed a reliable superposition method, which iteratively solves the Stokes disturbance flows for a many-particle system.That superposition method is, however, computationally expensive due to its iteration procedure.Onishi et al. (2013) later developed a less costly method, named the binary-based superposition method (BiSM), which has been adopted in the LCS (Onishi et al., 2015).BiSM assumes that interactions via three or more particles are negligible.This dramatically reduces the computational cost but maintains reliability as long as the particle number concentration is small, as observed in atmospheric clouds.Sundaram and Collins (1997) showed, by means of a DNS, that the preferential concentration of inertial particles, the so-called clustering effect, increases the collision frequency.The clustering effect is expressed in the spherical formulation derived by Wang et al. (1998) as where . . .denotes an ensemble average, |w r (x = R 12 )| (|w r | hereafter) is the radial relative velocity at contact separation, and g 12 (x = R 12 ) (g 12 hereafter) is the radial distribution function at contact separation and represents the clustering effect.Droplet deformation and coalescence efficiency, which this study ignores, affect the collision growth of droplets with r > 100 µm, although such effects only become significant for droplets with r > 500 µm.It would, therefore, lead to some errors if extending the present results to such large droplets.et al. (2008a) provided a parameterization for the turbulent geometric collision kernel of finite-inertia sedimenting droplets by proposing an empirical model for g 12 in addition to a theoretical model for |w r | .

Ayala
By following the expression by Chun et al. (2005), the clustering effect for a monodisperse suspension of sedimenting droplets is expressed as where η is the Kolmogorov length.C 1 is a function of St, Re λ , and the non-dimensional parameter for gravity V ∞ /v η with the Kolmogorov velocity v η .This parameterization was extended for a bidisperse system in a manner similar to that in Chun et al. (2005): where r L = max (r 1 , r 2 ) and C 1 follows the same expression for the monodisperse case at St max = max(St 1 , St 2 ) = St(r L ), and r d is a length scale of the acceleration diffusion experienced by the particles.When two particles in a pair are two different sizes, any fluid acceleration or gravity will induce a relative velocity.This effect yields a diffusion-like process in the system and tends to smooth out inhomogeneities in the particle pair concentration.Thus, r d is larger for larger |St 1 − St 2 | for the bidisperse case and a monodisperse suspension form is recovered for the case r d r L .It should be noted for the discussion in Sect 4.4 that the g 12 model was designed to show maximum clustering at St ∼ 1 and a higher droplet clustering for larger Re λ (Ayala et al., 2008b).
In addition to the empirical g 12 model, Ayala et al. (2008a) developed a theory for |w r | that is applicable to inertial droplets sedimenting under gravity in a turbulent flow.The basic assumption was that the droplet relative trajectory is mostly determined by gravitational sedimentation.Following Dodin and Elperin (2002), they decomposed the radial relative velocity (between two particles falling under gravity in a homogeneous isotropic turbulent flow) 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 gravity axis.The random variable ξ(φ) is assumed to be normally distributed with a standard deviation σ (φ).
Using σ (φ = 90 • ) to approximate σ (φ), they obtained where σ is expressed in terms of τ p,i , V ∞i and flow parameters u rms (the rms of the velocity fluctuations) in terms of and Re λ .
2.3 Onishi model 2.3.1 Model for g 12 Onishi et al. (2015) proposed an original model for the clustering effect in monodisperse systems of non-sedimenting particles with Stokes' linear drag.
where A 1 and A 2 were empirically determined to be 110 and 0.38, respectively.The regime boundary St a is λ .A tanh smoothing function, z a , was employed to connect the two formulations in the equation as If we limited the discussion for the autoconversion regime, i.e., r < 40 µm, the range St ≤ 1 would be enough for the typical energy dissipation rate ≤ 1000 cm 2 s −3 observed in atmospheric turbulent clouds.However, as clearly shown in Fig. 1, St can be as large as 10 for r = 100 µm and = 1000 cm 2 s −3 .That is, in the discussion on the accretion process that describes the conversion from cloud to rain due to rain drops collecting cloud droplets, we need to deal with St > 1 as well.
Hence, this study modifies the parameterization in the original Onishi kernel to obtain better overall matching for a wider range of St.After trial and error, we finally obtained a modification of the form of C a as We confirmed that this form with a c = 0.046, b c = 0.36, a c2 = 0.094, and b c2 = 0.25 leads to an improvement, as shown later in Sect.4.2.The updated coefficients are summarized in Table 1.
To determine the clustering effect for bidisperse systems, the empirical formulation proposed by Zhou et al. (2001) is employed: where The gravitational settling affects the clustering effect for large St particles.The parameterization here does not consider the gravity effect.This would lead to some error in collision statistics.However, the error was not significant in this study and the present parameterization worked well for predicting the turbulence enhancement in size evolutions due to collisional growth as in Sect.4.5.Onishi et al. (2015) employed the model of Wang et al. (2000) for |w r | , which was based on the model by Kruis and Kusters (1997), as

Model for |w r |
where ν is the kinematic viscosity and C w (St max ) = 1 + 0.6 exp −(St max − 1) 1.5 .The formulation of f KK was proposed by Kruis and Kusters (1997) as where θ i = τ p,i /T L with T L as the Lagrangian integral time and γ = 0.183u 2 rms /( ν) 1/2 .The Lagrangian integral time is parameterized as T L = 0.4T e , where T e (= u 2 rms / ) is the large-eddy turnover time (Kruis and Kusters, 1997;Zhou et al., 2001).In the equation, θ i shows the relative particle relaxation time to the particle-flow interaction time.Note that this |w r | parameterization is for non-sedimenting droplets.Onishi et al. (2009) concluded that gravitational sedimentation does not significantly influence turbulent collisions of cloud droplets.However, for this study, which extends the discussion to the small rain drop regime, the gravitational sedimentation cannot be ignored.Therefore, this study introduces a simple modification to make the model applicable to sedimenting droplets by considering the mechanism in which the gravitational settling shortens the interaction time of droplets with eddies (Onishi et al., 2009).Onishi et al. (2009) modeled the enlargement of the relative particle re- where f (θ ) is defined as the ratio of the particle velocity fluctuation to the flow velocity fluctuation, i.e., f (θ ) = v 2 p /u 2 rms , and s v = V p,∞ /u rms is a non-dimensional parameter quantifying the influence of sedimentation.By replacing θ i in Eq. ( 20) by θ i,sed , we obtain the radial relative velocity for droplets with gravitational sedimentation, |w r | turb,sed .
The above simple treatment is not yet complete.Ayala et al. (2008a) suggested the following two contributions of gravitational sedimentation on |w r | : (i) gravity reduces the interaction time of droplets with turbulent eddies, and therefore the variance of particle velocities is reduced, and (ii) gravity also decreases the correlation coefficient.The second contribution is missing in the present simple treatment.Nonetheless, since the present treatment leads to an improvement in the turbulent coagulation kernel, as shown in Sect.4.3, this study adopts this simple treatment and leaves more robust treatment to future work.
The turbulent collision kernel formulated from the above g 12 and |w r | turb,sed does not include the collision contribution due to the settling velocity difference.To include the contribution of the settling velocity difference, the following simple formulation was employed to obtain the total collision kernel.
Here, K c,turb denotes the turbulent collision kernel obtained by K c,turb = 2π R 2 12 |w r | turb,sed g 12 .

Turbulent enhancement on collision efficiency
Onishi et al. ( 2015) employed the collision efficiency values of Pinsky et al. (2001) (E c,PKS01 hereafter) and η E tabulated in Wang et al. (2008).These tabulated values spanned a relatively small range of particle sizes: the sizes of collector droplets (r 1 ) were 20, 30, and 50 µm and the size ratios (r 2 /r 1 ) were from 0.167 to 0.90.Later, Wang and Grabowski (2009) tabulated the preliminary values of the enhancement factor for a wider range of droplet sizes: r 1 = 20, 30, 40, and 50 µm and r 2 /r 1 from 0.0 to 1.0.Note that the data for r 2 /r 1 = 0.0 were simply set to the values for r 2 /r 1 = 0.0835.It should also be noted that Wang and Grabowski (2009) tabulated the enhancement factors against the Hall collision ef-ficiency (E c,Hall hereafter; Hall, 1980).Unfortunately, inconsistencies exist between the two collision efficiency models.We found differences that are sometimes much larger than 10 % of the mean between E c,PKS01 and E c,Hall , particularly for small and large r 2 /r 1 ratios, i.e., for r 2 /r 1 ∼ 0 and ∼ 1.These differences should be carefully compensated for in η E .Wang and Grabowski (2009) tabulated the enhancement on E c,Hall , η #Hall E .In fact, we observed an overestimation in turbulent enhancement on the autoconversion rate when we used η #Hall E for the SCE simulation with E c,PKS01 .For Table 2, we calculated Following Wang and Grabowski (2009), this study simply sets the values for r 1 ≤ 20 µm to those at r 1 = 20 µm, and similarly the values at r 1 = 60 µm to those at r 1 = 50 µm.The factor is set to unity for r 1 = 100 µm and larger.Also, following Seifert et al. (2010), for 100 ≤ ≤ 600 cm 2 s −3 , this study linearly interpolates/extrapolates between the values of η #PKS01 E at = 100 cm 2 s −3 and at = 400 cm 2 s −3 .For > 600 cm 2 s −3 the extrapolated values at = 600 cm 2 s −3 are used for η #PKS01 E .
3 Direct numerical simulations

Computational methods
We now solve the three-dimensional continuity and Navier-Stokes equations for incompressible flows: The kinematic viscosity ν is set to 1.5 × 10 −5 m 2 s −3 , which is the value for atmospheric air at 1 atm and 298 K.The last term on the right-hand side represents the external forcing needed to achieve a statistically steady state.This study employs reduced-communication forcing (Onishi et al., 2011), which is suitable for massively parallel finite-difference models, to maintain the kinetic energy with |k| < 2.5, where k is a wavevector.Spatial derivatives are calculated using fourth-order central differences.The conservative scheme of Morinishi et al. (1998) is employed for the advection term, and the second-order Runge-Kutta scheme is employed for time integration.To solve the velocity-pressure coupling, we use the highly simplified marker and cell (HSMAC) scheme (Hirt and Cook, 1972), which iterates until the rms of the velocity divergence becomes smaller than δ/ , where is the grid spacing and δ is chosen to be 10 −3 .The governing equations are discretized by using a cubic domain of length 2π L 0 , where L 0 is the representative length.Periodic boundary conditions are applied in all three directions.The flow cube is discretized uniformly into N 3 grid points, resulting in = 2π L 0 N .Under the limit of a large ratio of the density of the particle material to that of the fluid (ρ p /ρ f 1), the governing equation for water droplets is given by where V is the particle velocity, U is the air velocity at the position of the droplet, u is the disturbance flow velocity due to the surrounding droplets, and τ p is the particle relaxation time defined as τ p = (2/9)(ρ p /ρ f )(r 2 /ν), in which r is the particle radius.F impulse denotes the impulsive force due to collisions and g is the gravity vector (= (−g, 0, 0), where g is the gravitational acceleration).The ratio of the density of the particle material to that of the fluid, ρ p /ρ f , is set to 8.43× 10 2 at 1 atm and 298 K, and f is the drag coefficient defined as the ratio between the nonlinear drag and the linear drag (Rowe and Henwood, 1961).It should be noted that Eq. ( 26), which adopts the point-particle assumption, is inaccurate for large St particles whose radii are not small enough compared to the Kolmogorov scale.
The second-order Runge-Kutta method is used for the time integration.The flow velocity at the droplet position U is linearly interpolated from the adjacent grid values.This simple linear interpolation is justified through comparisons with the cubic Hermitian, cubic Lagrangian, and fifth-order Lagrangian interpolations from Sundaram and Collins (1996).The disturbance flow u, which denotes the hydrodynamic interaction, is calculated by using the BiSM (Onishi et al., 2013).The particle mass and volume fractions are so dilute that the flow modulation is ignored.

Computation for turbulent collision statistics
After the background airflow has reached a statistically stationary state, monodispersed water droplets are introduced into the flow.After a period exceeding 3 times the eddyturnover time T 0 = L 0 /U 0 , collision detection is then started.Droplets are allowed to overlap (ghost-particle condition) and a collision is judged from the trajectories of a pair of droplets by assuming linear particle movement for the time interval t.
The detailed description of the procedures for calculating collision statistics can be found in Onishi et al. (2013), who conducted the DNS for Re λ up to 530.This study performed additional simulations to push the maximum Re λ forward, up to 1140.The computational settings for the present simulations are summarized in Table 3.

Computation for size evolutions due to collisional growth
To obtain reference data regarding droplet collisional growth, we tracked the growth of droplets that initially had the following exponential size distribution (e.g., Soong, 1974): where x m0 is the mass of a droplet with a radius of r m0 and n 0 is the initial number density.We carried out two cases: one with r m0 = 15 µm and n 0 = 1.42×10 8 m −3 and the other with r m0 = 10 µm and n 0 = 4.79×10 8 m −3 .The corresponding initial liquid water content was 2.0 g m −3 for both cases.
It was assumed that colliding particles immediately united without breakups and conserved mass and momentum.Table 4 summarizes the computational parameters for the flow calculation as well as the obtained flow statistics for the Table 3. Case configurations and typical turbulence statistics.Re = U 0 L 0 /ν, u is the rms of flow velocity fluctuation, k max (= N/2) is the maximum wavenumber, l η is the Kolmogorov scale, and Re λ is the Taylor-microscale-based Reynolds number.N p is the total number of particles.collision growth simulations.In cases T100, T, and T1000, the same grid configuration with the same Reynolds number was calculated, but the energy dissipation rates, which are in the typical range observed in turbulent atmospheric clouds, were 100, 400, and 1000 cm 3 s −2 , respectively.Cases T, TR127, TR206, and TR333 obtained flows with the same energy dissipation rate (400 cm 3 s −2 ) but with different Re λ values.Onishi et al. (2015) have already presented these cases, except for TR333 with r m0 = 15 µm.The present study additionally performed the case TR333 with r m0 = 15 µm to obtain a clear Reynolds-number dependence, as well as cases T, TR127, and TR206 with r m0 = 10 µm.and Vassilicos (2014) clarified that the Reynolds-number dependence of g 11 observed for 1/3 < St < 1 is due to internal intermittency of the three-dimensional turbulence.
To quantify the influence of intermittence on g 11 , we need to separate the local quantity from the global (average) quantity.Kolmogorov (1962) introduced the local energy dissipation as where superscript # denotes the local quantity.It was supposed that the probability density function (PDF) of l follows a log-normal distribution if l is much smaller than the flow integral scale.Assuming l ∼ η, we obtain where * = η .Parameters σ and µ appear in the first and second moments of * as * (= ) = exp µ + σ 2 /2 (30) respectively.
The intermittency is measured by the flatness factor F , defined as It is observed that F follows a power law relation with Re λ , for example F ∼ Re 3/8 λ (Pope, 2000).Given (33) Substitution of Eqs. ( 30) and (31) into Eq.( 33) yields R. Onishi and A. Seifert: Re dependence of turbulence enhancement Eq. ( 30) then yields That is, P LN ( * |µ, σ 2 ) can be rewritten as P LN ( * |Re λ ).We can define a local St, St * , as the PDF of which follows It should be emphasized that the shape of P LN (and consequently P ) depends on Re λ .If we assume a universal radial distribution function at contact separation against St * − g #univ 11 (St * )−, the global clustering effect can be obtained as It should be noted that g 11 depends on Re λ , whereas g #univ 11 does not (which is why it is called universal).For St * 1, the universal clustering effect would have the form g #univ 11 = A 1 St * 2 + 1 by following Eq.( 10).Substitution of this form into Eq.( 38) yields g 11 (St 1, Re λ ) = A 1 St 2 +1, regardless of the value of Re λ .This explains why the g 11 for St = 0.1 does not show a significant Reynolds-number dependence.For a moderate St * , we simply formulate the universal function by following Eqs.( 10) and ( 11) but without the smoothing operators, as follows: where A * 1 and A * 2 are empirical parameters and St * a is defined as (A * 2 /A * 1 ) 1/4 .Based on the DNS data for St = 0.1, 0.4, and 0.6 in the flow with Re λ = 130, we found that A * 1 = 110 and A * 2 = 0.073 work reasonably well.Although we have no justification for this universal function, it can provide g 11 for arbitrary St (< 1) and Re λ through Eq. ( 38).As we cannot analytically calculate the integration in Eq. ( 38), we have to numerically calculate it to obtain g 11 for a certain combination of St and Re λ .We calculated g 11 for St = 0.1, 0.4, and 0.6 with Re λ = 100, 200, 400, 1000, 4000, and 10 000.We then obtained the following empirical formulations by applying the least-squares method to the calculated results.
Figure 2 shows a comparison between g 11 values from the above equations and those from the DNS.The figure shows  40), ( 41), and ( 42), which were fitted to the sample values (+) with using the least-squares method.The error bars show ±one standard deviation obtained from more than three runs, with each run lasting for a time T 0 = L 0 /U 0 .Black solid circle (right y axis) denotes the RDF at r/η = 0.25 for St = 0.6, normalized by the value at Re λ = 88 reported in Ireland et al. (2016a).
that the empirical estimates can reproduce the Reynoldsnumber dependence of g 11 correctly.The figure includes the data for St = 0.6 recently reported in Ireland et al. (2016a).The two DNS datasets for St = 0.6 agree well.Ireland et al. (2016a) argue that the Reynolds-number dependence is weak for the Reynolds number range explored by DNS, but the argument should be carefully examined when extrapolating to atmospheric flows.Figure 2 shows indeed a weak decrease in g 11 for Stokes numbers between 0.3 and 1.

Modeling of clustering effect
Figure 3 shows a comparison between direct numerical simulation results and model predictions for g 11 .The dashed lines are the prediction by the Onishi model (Onishi et al., 2015), and the solid lines are the predictions by the present updated model.The DNS data for St ≤ 1 and for Re λ ≤ 530 were obtained from the table in Onishi et al. (2015).The data for St = 1.4,2, 4, and 8 were newly obtained.The results for Re λ = 874 and 1140 (these Reynolds numbers are the largest ever achieved for turbulent particle collision statistics) are included in the figure.The DNS data show a decreasing trend for St < 1 for the moderate Reynolds number range of 100 Re λ 1000.This decreasing trend with respect to Re λ is attributed to the flow intermittency (Onishi and Vassilicos, 2014)  The updated parameterization leads to improvement, particularly for the St ≥ 1 regime.For example, in the case of Re λ = 127, the rms values of the relative errors of the prediction with the original parameters for (i) St = 0.1, 0.2, 0.4, and 0.6 and for (ii) St = 1, 1.4, 2, 4, and 8 were (i) 0.081 and (ii) 0.239.The rms values with the present parameters were (i) 0.075 and (ii) 0.113.

Turbulent coagulation kernels for small
Reynolds-number flow Figure 4 shows a comparison between model predictions and DNS results of the coagulation kernel K coag (r 1 , r 2 ) for r 1 = 30 µm, Re λ = 127, and = 400 cm 2 s −3 .The kernel is normalized by the collision radius R and the local velocity gradient λ = ( /ν) 1/2 .The reference DNS considers the hydrodynamic interaction and the gravitational droplet sedimentations.We observe a large discrepancy for r 2 ∼ 30 µm (= r 1 ), where the turbulence enhancement on collision efficiency is difficult to define, because the collision efficiency for r 1 = r 2 cannot be defined for stagnant flow.Otherwise, the model predictions (Ayala-Wang model and Onishi model) agree well with the DNS results.As an example, we also observe a slight improvement in the Onishi model by including the sedimentation effect on |w r | (Sect.2.3.2) on the data for r 2 = 40 µm.
The Ayala-Wang model shows a local maximum around r 2 = r 1 .The DNS results also show a convex shape, but the value at r 2 = r 1 is much smaller than the prediction by the Ayala-Wang model.In contrast, the Onishi model does not show such a local maximum at r 2 = r 1 but does provide values much closer to DNS elsewhere.The convex shape is related to the diffusion effect denoted by r d in Eq. ( 7).Equation ( 16) for g 12 , employed in the Onishi model, was formulated for non-sedimenting droplets and this equation therefore leads to weaker acceleration-driven diffusion, i.e., smaller r d (Ayala et al., 2008a).This can explain why the Onishi model does not show the convex shape.
Figure 5 shows the ratio of the turbulent coagulation kernel to the Hall kernel for the turbulent flow with Re λ = 127 and = 400 cm 2 s −3 .The level of the ratio is basically similar for both the Ayala-Wang and Onishi models, and the ratio is nearly unity when the droplets are above 100 µm.

Reynolds-number dependence of kernel models
Figure 6 shows the ratio of the coagulation kernel for Re λ = 104 to that for Re λ = 10 3 .It should be noted that the E c and η E models employed in the Ayala-Wang and Onishi kernels do not consider the Reynolds dependence.Therefore, the figure actually shows the ratio of the geometric collision kernels, i.e., the ratio of |w r | g 12 .The Ayala-Wang kernel increases for the autoconversion region (r 1 , r 2 < 40 µm) and the accretion region (r 1 < 40 and r 2 > 40, and r 1 > 40 and r 2 < 40 µm).The Onishi kernel decreases for the corresponding autoconversion region, but increases for the rainrain self-collection region (r 1 , r 2 > 40 µm).
Figure 7 shows the ratio of g 12 for Re λ = 10 4 to that for Re λ = 10 3 .It should be noted that the form of  Eq. ( 22) violates the spherical form and we cannot rigorously define g 12,total and |w r | total that formulate K c,total = 2π R 2 12 |w r | total g 12,total .Here, we simply considered g 12 expressed by Eq. ( 12) as the g 12,total for the total kernel and obtained |w r | total = K c,total / 2π R 2 12 g 12 .As designed, the Ayala-Wang kernel shows the increase for increasing Re λ for both the autoconversion and the accretion regions.In contrast, the Onishi kernel shows a decrease for the autoconversion region, but a significant increase for the accretion region and the rain-rain self-collection region (i.e., r 1 , r 2 > 40 µm).This is due to the shift of the maximum clustering toward larger St with increasing Re λ .
Figure 8 shows the ratio of the radial relative velocity for Re λ = 10 4 to that for Re λ = 10 3 .The Ayala-Wang kernel shows little Reynolds-number dependence.In contrast, the Onishi kernel shows significant Reynolds-number dependence, which tends to be opposite to the Reynolds-number dependence of g 12 and thus weakens the Reynolds-number dependence of the collision kernel.
The Reynolds-number dependence of the clustering effect is larger than that of the radial relative velocity, and the contour shape of Fig. 6 is more similar to Fig. 7 than to Fig. 8 for both the Ayala-Wang and the Onishi kernels.That is, the Reynolds-number dependence of the two kernels can mostly be attributed to the g 12 parameterizations.
Note that the Fortran 90 code used to calculate the present Onishi kernel is provided in the Supplement.

Turbulence enhancement of autoconversion rate
We investigated the turbulence enhancement on the autoconversion rate, which is the conversion rate from the cloud category (r < 40 µm) to the rain category due to collisions between the small cloud droplets.The Ayala-Wang kernel and the present Onishi kernel were employed to calculate the coagulation growth of droplets modeled by the stochastic collision-coalescence equation (SCE): where m is the particle mass and n f is the number density function.The coagulation component of the spectral bin model in the Multi-Scale Simulator for the Geoenvironment (MSSG-Bin) cloud physics model (Onishi and Takahashi, 2012) was used to solve the SCE.The mass coordinate m was discretized as m k = 2 1/s m k−1 , where s was set to 16.The representative radius of the first bin was 2.7 µm and 528 classes were calculated, the largest class of which had a representative radius of 5.4 mm.The SCE solution is basically a mean-field approximation.In contrast, the LCS acts as a reference model as it includes all turbulence effects directly in its Lagrangian particle simulation.Due to the high computational cost, however, the LCS is restricted to moderate Reynolds number (here up to Re λ = 333).Following Seifert et al. (2010), Onishi et al. (2015) used a quantitative measure of the turbulence enhancement focusing on the timescale of the autoconversion process.The time required for a cloud to convert 10 % of its cloud mass into rain category drops is expressed as t 10 % , which can be used as a measure of the autoconversion timescale.Then, we can define the turbulence enhancement factor, E turb , as where the overbar indicates the mean value.
Figure 9a shows E turb as a function of for Re λ = 66 in the r m0 = 10 µm case.The LCS data show an almost linear increase with increasing .Both the SCE simulation with the Ayala-Wang kernel (SCE-Ayala hereafter) and that with the Onishi kernel (SCE-Onishi hereafter) show the same trend with the LCS data, although the SCE-Ayala slightly over-estimates the enhancement.The maximum relative difference between the SCE-Ayala and SCE-Onishi kernels was as small as 22 % at = 500 cm 2 s −3 .Both the SCE-Ayala and the SCE-Onishi kernels show a kink at = 600 cm 2 s −3 , where the turbulence enhancement on collision efficiency levels off. Figure 9b shows E turb as a function of Re λ for = 400 cm 2 s −3 in the case of r m0 = 10 µm.The SCE-Ayala and the SCE-Onishi kernels show different trends: the SCE-Ayala predicts an increasing enhancement with increasing Re λ , while the SCE-Onishi predicts almost constant or slightly decreasing enhancement.The difference between the two SCE predictions becomes larger for larger Re λ , with the LCS result closer to the SCE-Onishi prediction.The difference between the SCE-Ayala and the SCE-Onishi kernels can be explained by the Reynolds-number dependence of the two kernels, as discussed in Sect.4.4.This Reynolds-number dependence is relevant, because the SCE prediction becomes very different at large Re λ .For example, at Re λ = 2 × 10 4 , the SCE-Ayala prediction is 2.5 times larger than the SCE-Onishi prediction.The LCS results for Re λ ≤ 206 support the SCE-Onishi prediction.
Figure 10 shows E turb for the r m0 = 15 µm case, which was also discussed in Onishi et al. (2015).This study additionally performed the simulation for Re λ = 333 to investigate the Reynolds-number dependence more clearly.Basically, the results are similar to those in the previous fig-  ure.In Fig. 10, the SCE-Ayala and the SCE-Onishi kernels show closer results for Re λ = 66, and both SCE-Ayala and SCE-Onishi slightly overestimate the enhancement for > 400 cm 2 s −3 .The difference between the two predictions at Re λ = 2 × 10 4 is larger: the SCE-Ayala prediction is 3.0 times larger than the SCE-Onishi prediction.The LCS results for Re λ up to 333 clearly support the SCE-Onishi prediction.
In summary, both Figs. 9 and 10 show that the SCE-Ayala and the SCE-Onishi kernels produce consistent results for low Re λ with about a 20 % difference at most, but the two show very different values at large Re λ : the SCE-Ayala prediction becomes larger than the SCE-Onishi by a factor of up to 3 in cloud turbulence.This clearly suggests a strong demand for collision growth data with larger Re λ to construct a more robust turbulent kernel.

Periodicity influence
As noted in Woittiez et al. (2009) and discussed in Appendix A in Ireland et al. (2016b), the periodicity of the computational domain may lead to errors for the settling particles with large St. Ireland et al. (2016b) defined the critical St, St crit , as where Fr is the Froude number (= a η /g, where a η is the Kolmogorov-scale acceleration), L(= 2π L 0 in this study) is the domain size, l is the integral scale, and u η is the Kolmogorov-scale velocity.For St larger than St crit , the periodicity problem may arise.Figures 4, 9, and 10 are for settling particles.For those figures, we have calculated St crit to check the periodicity problem.(i) For Fig. 4, St crit = 3.7, which corresponds to r crit = 75 µm; r crit is the radius of particle with St = St crit .The two plots from DNS, which correspond to r 2 = 80 µm and 120 µm, exceed r crit .However, since the two plots are more or less similar with the gravitational (Hall) kernel values, the turbulent contribution would be small compared to the gravitational settling contribution.That is the error due to the periodicity would not significantly affect the results.(ii) For Figs. 9a and 10a, r crit are 50, 65, and 70 µm for = 100, 400, and 1000 cm 2 s −3 , respectively.For Figs. 9b and 10b r crit are 65, 75, 85, and 90 µm for Re λ = 66.1, 127, 206, and 333, respectively.The enhancement factor E turb , shown in Figs. 9 and 10, was evaluated by t 10 % , which is defined as the time required for a cloud to convert 10 % of its cloud mass into rain category drops.The threshold between cloud and rain categories was set at r = 40 µm.That is, 10 % of particles, in mass and volume, are larger than 40 µm in radius at t = t 10 % by definition.For example, according to the DNS results, 3 % of particles are larger than 50 µm and only 0.9 % of particles are larger than 60 µm at t = t 10 % .The percentage of particles that are larger than 50 µm in radius may have some impact on t = t 10 % and consequently E turb .In this sense, the plot for = 100 cm 2 s −3 in Figs.9a and 10a, whose r crit is 50 µm, may contain some error associated with the periodicity problem.However, since E turb for that plot is nearly unity, indicating small turbulence enhancement, the periodicity problem does not change the present findings.

Conclusions
This study investigated the Reynolds-number dependence of turbulence enhancement on the collision growth of cloud droplets.The Onishi turbulent coagulation kernel proposed in Onishi et al. (2015) was updated by using the present direct numerical simulation (DNS) results for the Taylormicroscale-based Reynolds number (Re λ ) up to 1140.The following three components were updated: (i) the radial distribution function at contact separation of a monodisperse suspension of droplets, i.e., the clustering effect, g 11 ; (ii) the Atmos.Chem.Phys., 16,[12441][12442][12443][12444][12445][12446][12447][12448][12449][12450][12451][12452][12453][12454][12455]2016 www.atmos-chem-phys.net/16/12441/2016/radial relative velocity at contact separation, |w r | ; and (iii) the turbulence enhancement on collision efficiency, η E .We confirmed that the updated g 11 parameterization agrees better with DNS results than the original parameterization for Re λ ∼ 100.We also confirmed that the updated parameterization has better agreement with the Reynolds-number dependence of g 11 for the estimated values of St = 0.4 and 0.6.The model of radial relative velocity was updated to include the effect of the gravitational sedimentation of droplets.The comparison with the DNS results confirmed that the updated model for |w r | is better than the original one.The Onishi coagulation kernel employed the turbulence enhancement on collision efficiency η E , tabulated in Wang et al. (2008).The updated kernel is intended to adjust to more recent η E values, tabulated in Wang and Grabowski (2009).It should be noted that the collision efficiency E c in Pinsky et al. (2001) (E c,PKS01 ), which the Onishi kernel employs, is different from the E c in Hall (1980) (E c,Hall ), particularly for r 2 /r 1 ∼ 0 or ∼ 1.We proposed a compensation such that η E (in Wang and Grabowski, 2009), which shows the turbulence enhancement against E c,Hall , is applicable to the kernel with E c,PKS01 .The proposed compensation is simply to multiply η E in Wang and Grabowski (2009) by E c,PKS01 /E c,Hall .
The present Onishi coagulation kernel was compared with the Ayala-Wang kernel (Ayala et al., 2008a;Wang et al., 2008) together with the DNS values for Re λ = 66 and the energy dissipation rate = 400 cm 2 s −3 .For K coag (r 1 = 30 µm, r 2 ), both kernels show similar values comparable to the DNS values except for r 2 ∼ r 1 .For the nearly monodisperse case, the Ayala-Wang kernel overestimates the kernel but provides a sharp convex shape, i.e., a clear local maximum at r 2 = 30 µm, that agrees with the DNS data qualitatively.The Onishi kernel does not show such a convex shape due to weaker acceleration-driven diffusion on the clustering effect g 12 , but the kernel values are in fairly good agreement with the DNS.The Reynolds-number dependence of the two kernels was also compared.It was shown that the Ayala-Wang kernel in-creases for the autoconversion region (r 1 , r 2 < 40 µm) and the accretion region (r 1 < 40 and r 2 > 40, and r 1 > 40 and r 2 < 40 µm).In contrast, the Onishi kernel decreases for the autoconversion region but increases for the rain-rain selfcollection region (r 1 , r 2 > 40 µm).These Reynolds-number dependences can be attributed to the Reynolds-number dependence of the clustering effect.
We also compared the stochastic collision-coalescence equation (SCE) simulations for both kernels; one with the Ayala-Wang kernel (SCE-Ayala) and the other with the present Onishi kernel (SCE-Onishi).Lagrangian Cloud Simulator (LCS; Onishi et al., 2015) simulations were also conducted to obtain reference data of the turbulent enhancement on collisional growth, in particular the enhancement on the autoconversion rate.The SCE-Ayala and SCE-Onishi kernels show consistent results for Re λ = 66 with about a 20 % difference at most, but the two SCE simulations show a different Reynolds-number dependence, resulting in large differences at large Re λ .It should be emphasized that the SCE-Ayala prediction can become larger than the SCE-Onishi by a factor of up to 3 in the typical large Re λ range observed in cloud turbulence.These simulations clearly suggest a strong demand for reference collision growth data with larger Re λ from DNS or laboratory measurement to construct a more robust kernel model.This is our goal in future studies.

Data availability
Data for the present graphs are available from the corresponding author upon request.
The Supplement related to this article is available online at doi:10.5194/acp-16-12441-2016-supplement.

Figure 1 .
Figure1.Stokes number against the particle radius for various energy dissipation rates.

4
Results and discussion 4.1 Estimate for Reynolds-number dependence of clustering effect of small-St particles Onishi et al. (2013) observed that the clustering effect and consequently the collision kernel decreases as the Reynolds number increases for Re λ > 100 and St = 0.4.Later, Onishi

Figure 2 .
Figure 2. Radial distribution function (RDF) at the contact of monodisperse particles with St = 0.1, 0.4, and 0.6 against Re λ .The plotted open symbols are the reference DNS results.The lines are the results of Eqs.(40), (41), and (42), which were fitted to the sample values (+) with using the least-squares method.The error bars show ±one standard deviation obtained from more than three runs, with each run lasting for a time T 0 = L 0 /U 0 .Black solid circle (right y axis) denotes the RDF at r/η = 0.25 for St = 0.6, normalized by the value at Re λ = 88 reported inIreland et al. (2016a).

Figure 3 .
Figure3shows a comparison between direct numerical simulation results and model predictions for g 11 .The dashed lines are the prediction by the Onishi model(Onishi et al., 2015), and the solid lines are the predictions by the present updated model.The DNS data for St ≤ 1 and for Re λ ≤ 530 were obtained from the table inOnishi et al. (2015).The data for St = 1.4,2, 4, and 8 were newly obtained.The results for Re λ = 874 and 1140 (these Reynolds numbers are the largest ever achieved for turbulent particle collision statistics) are included in the figure.The DNS data show a decreasing trend for St < 1 for the moderate Reynolds number range of 100 Re λ 1000.This decreasing trend with respect to Re λ is attributed to the flow intermittency(Onishi and Vassilicos, 2014) as discussed in the previous subsection.The black solid line is the estimated g 11 for St = 0.4 and the black dashed line is for St = 0.6 (Eqs.41 and 42, respectively).The present Onishi model shows slightly better agreement with the DNS data in terms of the slopes in comparison with the

12450R.Figure 5 .
Figure 5. Ratio of the turbulent coagulation kernel to the Hall kernel in the turbulent flow with Re λ = 127 and = 400 cm 2 s −3 .(a) Ayala-Wang kernel and (b) the present Onishi kernel.

Figure 6 .
Figure 6.Ratio of the coagulation kernel for Re λ = 10 4 to that for Re λ = 10 3 .(a) Ayala-Wang kernel and (b) the present Onishi kernel.

Figure 7 .
Figure 7. Ratio of the clustering effect g 12 for Re λ = 10 4 to that for Re λ = 10 3 .

Figure 8 .
Figure 8. Ratio of the radial relative velocity at contact separation |w r | for Re λ = 10 4 to that for Re λ = 10 3 .

Figure 9 .
Figure 9. Turbulence enhancement factors for r c = 10 µm as a function of (a) the energy dissipation rate and (b) the Taylor-microscalebased Reynolds number Re λ .Re λ = 66 in (a) and = 400 cm 2 s −3 in (b).The error bars indicate the standard deviations.

Figure 10 .
Figure 10.Turbulence enhancement factors for r c = 15 as a function of (a) the energy dissipation rate and (b) the Taylor-microscale-based Reynolds number Re λ .Re λ = 66 in (a) and = 400 cm 2 s −3 in (b).The error bars indicate the standard deviations.

Table 1 .
Parameter values for g 11 model.

Table 4 .
Case configurations and typical turbulence statistics.Re = U 0 L 0 /ν, where U 0 is the representative velocity and L 0 is the representative length, u is the rms of the flow velocity fluctuation, k max (= N/2) is the maximum wavenumber, l η is the Kolmogorov scale, λ is the local shear rate, and Re λ is the Taylor-microscale-based Reynolds number.