Kinematic and dynamic collision statistics of cloud droplets from high-resolution simulations

We study the dynamic and kinematic collision statistics of cloud droplets for a range of flow Taylor microscale Reynolds numbers (up to 500), using a highly scalable hybrid direct numerical simulation approach. Accurate results of radial relative velocity (RRV) and radial distribution function (RDF) at contact have been obtained by taking advantage of their power-law scaling at short separation distances. Three specific but inter-related questions have been addressed in a systematic manner for geometric collisions of same-size droplets (of radius from 10 to 60 μm) in a typical cloud turbulence (dissipation rate at 400 cm2 s−3). Firstly, both deterministic and stochastic forcing schemes were employed to test the sensitivity of the simulation results on the large-scale driving mechanism. We found that, in general, the results are quantitatively similar, with the deterministic forcing giving a slightly larger RDF and collision kernel. This difference, however, is negligible for droplets of radius less than 30 μm. Secondly, we have shown that the dependence of pair statistics on the flow Reynolds number Rλ or larger scale fluid motion is of secondary importance, with a tendency for this effect to saturate at high enough Rλ leading to Rλ-independent results. Both DNS results and theoretical arguments show that the saturation happens at a smaller Rλ for smaller droplets. Finally, since most previous studies of turbulent collision of inertial particles concerned non-sedimenting particles, we have specifically addressed the role of gravity on collision statistics, by simultaneously simulating collision statistics with and without gravity. It is shown that the collision statistics is not affected by gravity when a < ac, where the critical droplet radius ac is found to be around 30 μm for the RRV, and around 20 μm for the RDF. For larger droplets, gravity alters the particle–eddy interaction time and significantly reduces the RRV. The effect of gravity on the RDF is rather complex: gravity reduces the RDF for intermediate-sized droplets but enhances the RDF for larger droplets. In addition, we have studied the scaling exponents of both RDF and RRV, and found that gravity modifies the RDF scaling exponents for both intermediate and large particles, in a manner very similar to the effect of gravity on the RDF at contact. Gravity is shown to cause the scaling exponents for RDF and RRV to level off for large droplets, in contrast to diminishing exponents for non-sedimenting particles.

kernel. This difference, however, is negligible for droplets of radius less than 30 µm. Secondly, we have shown that the dependence of pair statistics on the flow Reynolds number R λ or larger scale fluid motion is of secondary importance, with a tendency for this effect to saturate at high enough R λ leading to R λ -independent results. Both DNS results and theoretical arguments show that the saturation happens at a smaller R λ for smaller droplets. Finally, since most previous studies of turbulent collision of inertial particles concerned nonsedimenting particles, we have specifically addressed the role of gravity on collision statistics, by simultaneously simulating collision statistics with and without gravity. It is shown that the collision statistics is not affected by gravity when a < a c , where the critical droplet radius a c is found to be around 30 µm for the RRV, and around 20 µm for the RDF. For larger droplets, gravity alters the particle-eddy interaction time and significantly reduces the RRV. The effect of gravity on the RDF is rather complex: gravity reduces the RDF for intermediatesized droplets but enhances the RDF for larger droplets. In addition, we have studied the scaling exponents of both RDF and RRV, and found that gravity modifies the RDF scaling exponents for both intermediate and large particles, in a manner very similar to the effect of gravity on the RDF at contact. Gravity is shown to cause the scaling exponents for RDF and RRV to level off for large droplets, in contrast to diminishing exponents for non-sedimenting particles. demonstrated that small-scale turbulent motion can enhance the collision rate of droplets either by augmenting the relative velocity and collision efficiency or through inertia-induced droplet clustering; see the latest review articles on this subject by Grabowski and Wang [1] and Devenish et al [2]. A quantitative description of the effects of air turbulence remains challenging due to experimental difficulties in probing statistics at droplet scales and computational difficulties in simulating all relevant scales of the turbulent flow.
Direct numerical simulations (DNS) of multiphase flows have emerged as an important research tool for studying collision-coalescence of cloud droplets. Spectral-based DNS provides an accurate and quantitative representation of small-scale physical processes occurring in clouds. The rapid progress in both the numerical methods and availability of high-performance supercomputers allows DNS to address flows at a somewhat wider range of Reynolds numbers; thus its relevance to applications and its value for providing quantitative data continue to grow. Earlier DNS implementations of particle-laden turbulent flows were single-threaded applications designed for vector computers [3,4]. They are limited to low Taylor-microscale flow Reynolds numbers (typically ∼50 or less). As a result, only the dissipation range of the kinetic energy spectrum was adequately represented in the simulations. In addition, memory constraints did not permit tracking of a large number of particles; consequently collision statistics had significant numerical and physical uncertainties.
The emergence of multiprocessor computers with shared memory opens up new opportunities for developing multi-threaded codes. They allow higher grid resolutions to be handled within a reasonable wall-clock time. Standard OpenMP parallel programming libraries may be employed to develop multi-threaded codes. However, there are several limitations. Firstly, both the CPU memory and the number of processors are strictly dependent on machines and very limited. For example, each node in the IBM Power 575 cluster (4064 POWER6 processors running at 4.7 GHz) has 32 processors and a maximum of 64 GB memory. Secondly, the codes cannot scale well due to a bottleneck in CPU-to-memory connection. In terms of problem size, parallelization with OpenMP could efficiently handle problems up to 128 3 flow grid points with O(10 5 ) droplets, yielding a maximum Reynolds number of R λ ∼ 10 2 [5].
Since 2003, clock frequencies of CPUs had begun to stagnate or even decrease (due to energy constraints). The number of cores per CPU started to double every 2-3 years. Therefore, it was necessary to develop a new implementation suitable for computers with architecture based on distributed memory. The main difficulty in developing such codes is uniform distribution of tasks among processing units. In the simulation of droplet collision-coalescence in a turbulent flow, one must deal with two different representations in a single simulation: Eulerian grid for turbulent flow and Lagrangian motion of droplets (particle tracking). Our first DNS implementation on a supercomputer with distributed memory was based on one-dimensional (1D) domain decomposition (DD) and used an MPI (message passing interface) library for data communication. Scalability of the code (without droplet-droplet aerodynamic interaction) up to O(100) cores has been presented in [6,7]. This implementation was a significant step forward. It enabled utilization of a larger number of processors (up to the number of grid points in the decomposed direction), large memory size and improved cache utilization, leading to a higher overall computational efficiency.
The first objective of this study is to present a more-scalable MPI implementation, based on two-dimensional (2D) DD. The 2D DD implementation divides data in two of the three coordinate directions. This dramatically increases the number of processors one can deploy, at the expense of data communication through increased subdomain interfaces. This approach Figure 1. Sketches show two spatial DDs: the left panel is 1D decomposition with 8 subdomains, and the right one shows 2D decomposition with 16 subdomains. Computational fluid nodes and droplets in each subdomain are assigned to an individual processor. X , Y and Z are the three directions of the coordinate system and g points in the gravity direction. For 1D implementation, the domain is decomposed in the Y direction while for 2D implementation decomposition is performed in the Y and Z directions.
makes it possible to perform simulations at higher flow Reynolds numbers. Figure 1 illustrates the concepts of 1D and 2D DDs. The computational domain is decomposed along the directions perpendicular to gravity in order to minimize the number of droplets crossing the subdomain boundaries. MPI is used to communicate between the subdomains whenever the solution or data collection procedure requires data from neighboring processors. We will show scalability of the newly developed 2D DD implementation relative to the 1D DD implementation, and also validate the new implementation by comparing results with our previous results and results from other studies.
The second and major objective of the paper is to present some first results from the newly developed implementation addressing both kinematic and dynamics statistics related to turbulent collision of cloud droplets, for a wider range of flow Reynolds numbers than previously possible. In order to obtain more accurate data for the radial relative velocity (RRV) and radial distribution function (RDF) at contact, we also study the power-law scaling exponents of these pair statistics. Three specific but inter-related questions will be addressed in a systematic manner for geometric collisions of same-size droplets (of radius from 10 to 60 µm) in a typical cloud turbulence (dissipation rate at 400 cm 2 s −3 ). The first question concerns the sensitivity of the simulation results on the details of the large-scale driving mechanism. We will employ two different forcing schemes for all runs to study quantitatively any difference in the collision statistics. The second question concerns the effect of flow Reynolds number or equivalently the range of flow scales represented in hybrid direct numerical simulation (HDNS). Since most previous studies of the turbulent collision of inertial particles concerned non-sedimenting particles (i.e. gravity was not considered), our third question concerns the role of gravity on collision statistics. We address this by simultaneously simulating collision statistics with and without gravity. Together, we will not only clarify the validity and limitations of the simulation tool, but will also obtain a further understanding of all factors affecting the collision statistics.
The remainder of the paper is organized as follows. The numerical method and implementation details for scalable computation are described in section 2 for the background airflow turbulence and in section 3 for the dynamics of droplets, respectively. The scalability data for the newly developed 2D DD code are compared with the previous 1D DD code in Appendix B. Physical results concerning droplet pair statistics are presented in section 4, with a specific focus on large-scale forcing, gravity and the effects of flow Reynolds number. Key conclusions are summarized in section 5.

Simulation of the background airflow
The basic ideas and algorithms for the hybrid DNS (HDNS) approach have been presented in [8]. Therefore, only a brief description of the approach is given here. The first step in HDNS is to develop homogeneous turbulent flow in a cubic domain of size (2π) 3 with periodic boundary conditions. The flow is simulated by solving the incompressible Navier-Stokes (N-S) equations: Here ω ≡ ∇ × U is the vorticity vector, P is the pressure, ρ is the fluid density and ν is the fluid kinematic viscosity (i.e. the air viscosity). To achieve stationary turbulence, the flow is driven by the forcing term f(x, t) which is non-zero only for a few low-wavenumber modes in Fourier space. Discretization of the domain is carried out by uniform division of the domain into N 3 grid points, where N is the number of grid points in each direction and takes the values of 2 n , n being a positive integer. The restrictions on N come from the use of the fast Fourier transform (FFT) algorithm. At these grid points, the turbulent fluid velocity is computed in spectral space yielding the physical fluid velocity through the inverse Fourier transform. Further details of the pseudo-spectral method are given in [3]. The use of the pseudo-spectral method for flow simulation has several advantages, such as the simplicity of implementing the large-scale forcing and imposing periodic boundary conditions, and high computational accuracy. On the other hand, use of the pseudo-spectral method involves significant computational complexity, namely, the intense (i.e. global) communication between processes due to parallel implementation of FFT. So far, only a few 3D FFT implementations for distributed memory computers have been documented in the literature. For our previous HDNS code based on 1D DD, we used the algorithm developed by Dmitruk et al [9]. Our new implementation utilizes a 2D DD of the 3D data field and a sequence of 1D FFT from the FFTW library (www.fftw.org) in each spatial direction. This new implementation represents the optimal balance between computation and communication on a scalable computer with O(100 000) processors, as explained in [10].
The key consideration of the implementation is a load-balanced efficient implementation of transpose operations, and the relevant details are presented in [10].
To initialize the velocity field we used a random phase algorithm with a prescribed Kolmogorov energy spectrum as E(k) ∼ |k| −5/3 [11], where k is the wave vector. Starting with such an initial setting and integrating the N-S equation with continuous injection of the kinetic energy at large scales, we obtain statistically stationary turbulent flows. Supplying the kinetic energy to the system at large scales is managed by a forcing scheme. Two different forcing schemes have been implemented. The first is a deterministic scheme (as in [12]) in which the energy levels of the two lowest wavenumber shells (0.5 < |k| < 1.5 and 1.5 < |k| < 2.5) are fixed to E(1) = 0.555 44 and E(2) = 0.159 843, respectively.
The second is the stochastic forcing scheme of Eswaran and Pope [13], in which random acceleration is added to each component of velocity. The specification of the acceleration forcing is based on six Uhlenbeck-Ornstein random processes [13]. There are two characteristic parameters defining the stochastic processes, namely, an acceleration variance σ 2 f and a forcing time scale t f . The time scale t f must be much smaller than the eddy turnover time (T e ), otherwise the stochastic forcing may not be time uncorrelated and the level of energy input will be reduced. Eswaran and Pope [13] showed that the average rate of energy input (which is also the average dissipation rate) could be expressed as where the number of modes being forced is N f = 80, the lowest wavenumber is k 0 = 1 and the fitting coefficient β was found to be 0.8 in Eswaran and Pope [13]. Since the term (σ 2 f t f N f k 2 0 ) −1/3 represents roughly the time scale of the large eddies, the group t f (σ 2 f t f N f k 2 0 ) 1/3 could be interpreted as the ratio of forcing time scale to the time scale of large-scale flow.
To keep the dissipation rate roughly the same for different grid resolutions, we set t f roughly equal to the Kolmogorov time scale using a reference dissipation 0 = 3600. This is done by setting Then equation (3) becomes where η 0 ≡ (ν 3 / 0 ) 1/4 . Clearly, as k 0 η 0 → 0, then → 0 . The time step size used for integration of the N-S equations is chosen to balance numerical stability, accuracy and computational efficiency. The spatial resolution of the simulation was monitored by k max η, which should be greater than unity for fine scales to be resolved. On the other hand, spatial resolution (k max η) should be kept close to unity in order to maximize turbulent Reynolds number R λ . The use of two different forcing methods for a given grid resolution results in somewhat different flow characteristics, allowing a sensitivity study of the collision statistics on the forcing scheme and large-scale flow. Ayala et al [5] postulated that the deterministic scheme may bring more coherence to turbulent eddies, which could lead to somewhat larger collision kernels.

Turbulent transport and interaction of droplets
Droplet tracking starts from the moment when the flow becomes statistically stationary. Droplets are initially introduced into the flow at random locations. To achieve a truly random distribution of their location, a single core creates random locations for all droplets in the domain. Then, the coordinates of the droplets are broadcast to their host subdomains. This task is performed at the onset of the simulation and the cost is negligible compared to that of the full simulation.
Once the undisturbed fluid velocity U(Y (k) (t), t) is computed, droplets are advanced by solving their equation of motion, which for the kth droplet is written as where τ (k) p = 2ρ p (a (k) ) 2 /9µ is the Stokes inertial response time of the kth droplet, µ is the air dynamic viscosity, g is the gravitational acceleration and u (k) is the air perturbation velocity originating from the motion of surrounding droplets.
The location and velocity of each droplet are advanced by integrating the equation of motion with a fourth-order Adams-Moulton scheme (chapter 16.7 in [14]) for droplet velocity (equation (6)) and a fourth-order Adams-Bashforth scheme (chapter 16.7 in [14]) for droplet location (equation (7)). Two additional tasks have to be completed before equations (6) and (7) can be advanced in time. Firstly, interpolation of the fluid velocity from the regular grid to the location of the droplets must be performed. Secondly, Stokes disturbance velocities u (k) must be specified. These will be discussed next.

Velocity interpolation
Interpolation of the fluid velocity from the grid points to a droplet center requires substantial computational effort, since this needs to be done for a large number of droplets at every time step. For parallel implementations based on spatial DD, the interpolation algorithm requires intensive communication between processes in order to obtain data from a sufficiently large volume surrounding every droplet. Therefore, optimization of this task is essential to a highly scalable HDNS code. The size of the interpolation volume (and the number of grid points within) determines the cost of the interpolation. The size of the interpolation volume is specified by the order (thus the precision) of the interpolation method. Several different interpolation techniques have been developed in the past; see a brief review of the commonly employed methods in [6]. In the new code we use the six-point Lagrangian interpolation in each spatial direction. A detailed description of the algorithm is presented in [15].

Aerodynamic interaction between droplets
The disturbance flow caused by each droplet is assumed to be a localized Stokes flow. Ayala et al [8] showed that this assumption is appropriate for the cloud physics application for the following reasons. The droplet Reynolds number is typically small, so the disturbance fluid motion around a droplet follows Stokes flow. In all, 85% of the total viscous dissipation in this Stokesian-disturbance flow is contained in a spherical region of ten times the droplet radius. Therefore, the kinetic energy associated with the disturbance flow converts to viscous dissipation locally and quickly. Since the droplets are one to two orders of magnitude smaller than the Kolmogorov eddy, the disturbance flow is contained within a turbulent energydissipation eddy. If these Stokes disturbance velocities of all N p droplets in the system are superposed appropriately and the no-slip boundary condition on the surface of each droplet is considered, the disturbance velocity u (k) felt by the kth droplet is governed by a linear system [16,17] where the Stokes disturbance flow is given as In equation (8), the disturbance flow velocity felt by a target particle k is related to the disturbance flows of all neighboring particles (denoted by index m). In this paper, the terms particles and droplets are used interchangeably. Therefore, the disturbance velocity acting on a given droplet depends on the background turbulent flow, as well as the positions and velocities of all the other droplets in the system. It is important to note that the periodicity of the domain holds for particles and disturbance velocity u (k) as well as for the background flow.
Given particle radius a (k) , velocity V (k) , position Y (k) and the background flow field U(Y (k) , t) at a specified time t, (8) is a system of equations of 3N p unknowns for the three components of the disturbance velocities u (k) of the N p droplets. The disturbance velocity at the location of each droplet is coupled with the disturbance velocities of all the other droplets and the three spatial components of the disturbance velocity cannot be separated to form three independent smaller linear systems. Therefore, system (8) must be solved as a whole to yield the disturbance velocities u (k) . In our simulations, we solve the system using an efficient iterative scheme known as the generalized minimum residual algorithm GMRES [14] with constant coefficients in equation (8) computed in advance.
The solution of (8) could be significantly accelerated by truncating the aerodynamic disturbance flow of a droplet (with radius a) to some optimal distance H trunc × a (k) . H trunc is a non-dimensional value, which may be determined by analysis of the statistical data from previous numerical experiments. As far as collision statistics are concerned, Ayala et al [8] showed that one could truncate the summation in (8) and restrict the radius of influence of the disturbance flow. Their experiments showed that the computed collision efficiency is insensitive to H trunc if H trunc > 35. In the new implementation the truncation radius is set to H trunc = 50 which is more conservative than H trunc = 35.
To limit the data communication to the nearest-neighbor subdomains, we assume that the truncation sphere of any individual droplet at most covers a volume belonging to immediate neighboring subdomains. Since the droplet radii are in the range of 10-60 µm and are much smaller than the Kolmogorov length of the undisturbed flow, this assumption does not pose a severe restriction on the maximum number of cores. This assumption amounts to the upper bound on the number of subdomains in each decomposed direction, as where N subd1 is the number of subdomains in a single decomposed direction (namely, P y or P z ), a max is the maximum radius of the droplets in the system, L box is the domain size in the decomposed direction and [.] denotes the integer part of a real number.

Parallel implementation of droplet tracking
In our new implementation, communication of data associated with droplets is implemented by proper use of droplet indices and careful rearrangement of droplet data when droplets leave or enter a subdomain. In order to efficiently locate the droplets and their immediate neighboring droplets inside the truncation sphere, we make use of the cell-index method and the concept of linked lists [18]. This is numerically implemented by dividing the domain into small cells and keeping track of droplet indices inside each cell. This speeds up the process of finding neighbor droplets, without affecting the physical results. To be able to track a specific particle over the whole domain, we also introduce the droplet global index. This global index is unique to each droplet and will be passed to the host subdomain if the droplet changes subdomain.

Physical results
In this section, we report the DNS results based on the 2D DD implementation. The results include characteristics of the background turbulent flow field and kinematic and dynamic statistics related to collision-coalescence of non-interacting droplets. We will limit our discussions to turbulent geometric collisions between droplets of the same size. Droplets of radius from 10 to 60 µm will be considered. The size increment is 2.5 µm for droplets smaller than 30 and 5 µm for the larger droplets. The physical flow dissipation rate is assumed to be c = 400 cm 2 s −3 in our simulations with physical viscosity set to that of air. These are used to scale the DNS units to match the conditions of cloud droplets, as explained in Ayala et al [5]. Under this physical condition, the Stokes number varies from 0.063 to 2.28 and S v from 0.446 to 16.1.

Background turbulent flow
Tables 1 and 2 list the average values of key flow parameters obtained from the simulations performed with both stochastic and deterministic schemes. The grid resolution is varied from N = 32 to 1024. The quantities shown in tables 1 and 2 are the kinematic viscosity ν, the time step size δt, the rms fluctuating velocity u , the energy dissipation rate , the Taylor microscale Reynolds number R λ = u λ/ν, the spatial resolution parameter k max η, the Kolmogorov length η, the Kolmogorov time τ k , the integral length scale L s [19], the integral length scale of the longitudinal spatial velocity correlation L f , the transverse Taylor microscale λ, the eddy turnover time T e , skewness S and flatness F of the velocity gradient, and the CFL number. In these tables, the statistical standard deviation σ A of any given time-averaged quantity A has been calculated, as in Eswaran and Pope [13], by σ 2 A = 2σ 2 A T / T , where A represents the time average of A (typically a volume-averaged quantity at a given time instant), T is the integral time scale of A, T is the total time duration and σ A is the standard deviation of A. The integral time scale T was estimated by the time delay when the correlation coefficient takes a value of 0.5.
For the simulated flows using the stochastic forcing, equation (5) can be used to predict the average dissipation rate, as shown in figure 2. It appears that a modified fitting parameter β = 2.0 provides a much better prediction than β = 0.8. This also explains why the dissipation rate increases with grid resolution N . As N is increased, the Kolmogorov scale is reduced relative to the box size. This leads to decreasing k 0 η 0 ; thus gradually approaches 0 from below. Analytical prediction of the energy dissipation rate given by equation (5) is based on  assumptions that scales (roughly) with N f 0 , and that the integral length scale scales with 1/k 0 . The larger standard deviation for the 1024 3 run is due to an inadequate time interval used for averaging. Based on the data in tables 1 and 2 we observe that for a given grid resolution, the deterministic forcing yields a larger R λ . The highest Reynolds number obtained at 1024 3 is R λ ∼ 500. Figure 3 shows the dependence of R λ on the mesh resolution N. Additionally, results from other independent studies [20][21][22][23][24] have been added for comparison. Using a deterministic forcing method, Franklin et al [24] obtained lower R λ mainly because they used a larger k max η  (in the range of 1.8-2.7). In our simulations, k max η has been chosen to be close to unity (see tables 1 and 2). Ishihara et al [20] used a deterministic forcing scheme yielding two sets of flow Reynolds numbers corresponding to two grid resolutions, namely, k max η ≈ 1 and ≈ 2. As shown in figure 3, their results with k max η ≈ 1 are in good agreement with our results. As for Bec et al [23], they tuned the fluid viscosity value to ensure η ≈ x. Furthermore, small differences in the details of the forcing method and also dealiasing might cause some variations in R λ . In order to validate the new (2D DD) DNS code for flow simulation, 1D energy spectra of the simulated flows have been computed and compared in figure 4 with experimental data of grid-generated turbulence from two different wind tunnel experiments [25,26]. In the experiments, the measurement stations were located far from the turbulence-inducing grids, to ensure their flow is homogeneous and isotropic. We chose the 1D energy spectrum for comparison because this quantity can be directly measured in the wind tunnel. Moreover, the 1D energy spectrum can be more accurately computed than for example the three-dimensional spectra, since the modes are more evenly distributed [3]. We conclude that, in the inertial  [25,26] are plotted for comparison. and dissipation subranges, the experimental spectra agree very well with energy spectra from our numerical simulations. Moreover, there is no significant difference between simulations performed with different forcing methods.

Radial distribution function
RDF is a measure of the effect of preferential concentration of droplets on the collision rate. Numerous efforts using different approaches, i.e. analytical [27,28], numerical [22,[29][30][31][32][33] and experimental [33,34], have been made previously to compute or measure the RDF. Most of the studies demonstrate that for separation distances r smaller than the flow Kolmogorov scale η, the RDF has the following power-law dependence: where the exponent c g (St, S v ) depends on the particle Stokes number St = τ p /τ k and velocity ratio S v = v p /v k . v k and v p are the Kolmogorov velocity and particle terminal velocity (settling velocity of an isolated particle in stagnant air). The second parameter c 0 is called the power-law pre-factor. A convenient method to compute the RDF on the fly during simulations was proposed in [35,36], with the RDF at contact (r = R where R = a 1 + a 2 ) being computed directly from the definition where i indicates a particular droplet radius and N pairs is the total number of pairs detected with separation distance (r ) falling in a spherical shell of inner radius equal to R − δ and outer radius equal to R + δ. Here δ is a small fraction (1%) of R [35,36]. V s is the volume of the spherical shell, . N i is the total number of a i droplets used in the simulation and V B is the volume of the computational domain. g ii (r ; t) is further averaged over time to obtain g ii (r = R).
This method is fairly efficient and simple to implement but may be subject to large numerical uncertainties if N i is not large, which could occur if the volume concentration of particles is low. Therefore, in this study we instead computed the RDF at r = R using a more accurate two-step approach. In the first step, the RDF is computed in the usual manner as in [35] but for varying separation distances. In the second step we fit the discrete bin data to the powerlaw equation (11) so the value of the RDF at contact can be more accurately estimated. In order to perform the first step, we generalize equation (12) to any shell with the inner radius equal to r − δ and the outer radius equal to r + δ, where r ∈ [R, 10R] is discretized into 180 bins of equal width 2δ. The similar methodology and power-law scaling have been used to process data obtained in the simulations with gravity turned on and off, to quantify the effect of gravity on the RDF and scaling exponent. Numerical uncertainties were estimated following the method described in chapter 16.7 in [14].

Scaling exponent.
In figure 5, we plot the RDF as a function of separation distance r . All kinematic statistics were computed in the simulations with 5 million particles at 512 3 grid resolution. The data exhibit a power-law behavior for all the droplet sizes considered, with the power-law exponent (i.e. the slope on the log-log plot) varying with the droplet radius.
For monodisperse droplets, the RDF is directly related to the level of clustering [1]. The power-law exponent c g in equation (11) is a key parameter measuring the level of clustering of the droplets at scales smaller than the Kolmogorov length. Theoretical analysis of clustering of the inertial particles at small distances, limited to small Stokes numbers, has been developed by Chun et al [28] where they showed analytically that in the limit St 1, c g is proportional to St 2 . Other studies have extended the expression for c g to cover larger Stokes numbers; see [37] for example. However, these theoretical analyses often involved various assumptions and their predictive capability must be checked against DNS and experimental data.
In figures 6 and 7 we plot c g data from our simulations and compare with theoretical estimations. Overall, our results without gravity are consistent with theoretical predictions derived in [28,31,37] and other independent numerical studies reported in [22,29,33]. Several important observations can be made. First, for Stokes numbers greater than ∼0.3, the analytical approximation c g ∼ St 2 [28] is no longer valid. This discrepancy between DNS data and theoretical prediction reflects a number of approximations used in the analysis of Chun et al [28], who assumed that the inward relative drift caused by the particle inertia scales linearly with the particle separation distance, and a pairwise diffusivity owing to turbulent fluid motion scales quadratically with the particle separation distance. Evaluation of the drift flux is based on an approximation that the second invariants of the rate of strain and rotation ( S 2 p and R 2 p ) depend linearly on the Stokes number, an assumption that may not be satisfied for particles with large inertia. The linear theory does not take into consideration nonlinear effects due to the depletion of particles in regions of high rotation. Additionally, the non-local diffusion flux has been evaluated using the limit of St 1 for which the leading order term reduces to a constant. These factors may explain the differences between the theory of Chun et al [28] and DNS data in figure 6(a). For the cases without gravity, c g reaches a maximum value at St ∼ 0.6 and then decreases with the Stokes number. The empirical formula of Derevyanko et al [37] appears to overestimate c g for intermediate Stokes numbers, and underestimates c g for large Stokes numbers.
The effect of gravity on the scaling exponent has received little attention so far, so we compare directly our results with and without gravity in figures 6 and 7. There is a significant effect of gravity on the value of c g , for both intermediate and large droplet sizes. For small to intermediate droplet sizes, gravity reduces c g slightly. However, for large droplet sizes, gravity increases c g significantly. In fact, for cases with gravity c g does not decrease much after reaching the maximum value as St is increased, which is different from the results without gravity. This could be related to the preferential sweeping [3] which causes strong small eddies to contribute more to the clustering under gravity. Furthermore, there is a weaker dependence on the flow Reynolds number when compared to results of non-settling particles. This may be explained by the fact that gravity reduces the interaction time of particles with highly intermittent eddies, and also causes eddies of certain size and strength to make a more dominant contribution to the clustering. The only published results on c g with gravity from Falkovich and Pumir [31] are shown for comparison. The qualitative trend of their results is very similar, although their values are larger than ours by around 10-15%. The results together suggest that the effect of gravity on the scaling exponent c g for cloud droplets needs to be parameterized in the future.

Radial distribution function at contact.
We shall now present the monodisperse particle RDF at contact which is directly relevant to the collision rate. Figure 8 shows the results as a function of particle radius for the two forcing methods. Without gravity, the RDF increases very quickly with particle radius a when a < 30 µm. A peak value is reached when 30 µm < a < 40 µm. The RDF then drops quickly as a is increased for a > 40 µm. This corresponds well to the behavior of c g with droplet size when gravity is removed. The peak RDF could reach a value of around 30, implying a significant effect of clustering on the collision rate of nonsettling monodisperse particles. The RDF also depends on flow Reynolds numbers, but the level of dependence changes with particle size. We will return to this point in section 4.5. The results are very similar for the two forcing schemes.
Gravity significantly affects both the value and dependence of RDF on droplet radius. Sedimentation suppresses the peak value by reducing particle-eddy interaction time, and the peak value drops to around 20. The large-scale forcing mechanism appears to have little effect on the particle sizes where the RDF peaks. Overall, the results depend less on the flow Reynolds number when the deterministic forcing is applied. We suspect that this is caused by the differences in the level of coherence and size of turbulence vortices derived from the different forcing methods. The RDF does not drop significantly for large droplets, very similar to the behavior of c g with gravity. Therefore, gravity tends to sustain the effect of clustering on the collision rate for larger droplet sizes, at least for collisions between particles of the same size. Also there is a weaker dependence on flow Reynolds number when the gravity is considered, compared to the results with gravity. forcing schemes. Two black-dashed lines at (a) represent theoretical prediction from [38]. G marks the cases for sedimenting droplets and others for nonsedimenting particles.
In figure 9 we compare the RDF with theoretical prediction derived by Ayala et al [38] as well as numerical results from independent DNS studies in [33,39]. The prediction by Ayala et al [38] is acceptable for small droplets (a < 25 µm), but tends to underpredict the RDF Figure 9. Comparison of simulated RDF with results developed by Saw [33] and Woittiez et al [39]. The thin dashed line represents the analytical solution developed by Ayala et al [38]. Simulations with gravity are labeled as G. D and S indicate the forcing method, i.e. deterministic and stochastic, respectively.
for larger droplets and when the deterministic forcing is used. For sedimenting droplets, our results using the deterministic forcing are in good agreement with the RDF results of [39] despite the fact that the flow Reynolds numbers are not exactly matched (81.4 versus 84.9). We note that the forcing method used in [39] is a modified version of deterministic forcing used in the present study. For non-sedimenting droplets, our results and those of [39] are in good quantitative agreement. It is quite surprising though that our results with stochastic forcing are almost identical to their data for non-sedimenting droplets. Apparently, the numerical results reported in [33] underpredict the RDF, considering the higher R λ employed in their simulations (R λ = 143).

Radial relative velocity
The radial relative velocity w r is defined in terms of the relative velocity w between two droplets with the separation vector r as w r ≡ w · r/|r|, with r ≡ |r|. The methodology for computing the RRV between particles is similar to that described in the previous section for the RDF. Again, we fit the normalized RRV from discrete bins to a power-law relationship defined as Figure 10 shows normalized RRV as a function of the separation distance. Evidently, the power-law scaling works well for all droplet Figure 10. RRV of sedimenting droplets as a function of separation distance normalized by Kolmogorov velocity. Black lines represent the power-law fit to the DNS data. The simulations were performed at grid resolutions of 512 3 using the deterministic forcing scheme.

Scaling exponent in the dissipation range.
sizes. In the limit of small droplet radius, the power-law exponent approaches unity, since droplets with very low inertia behave as fluid elements and the relative velocity is linearly proportional to r since the local strain rate is essentially uniform when r < η. The scaling exponent c w is obtained from figure 10 and plotted in figure 11 as a function of Stokes number for different R λ and two different forcing schemes. Since r is smaller than η, this is the scaling exponent for the dissipation range length scales. The variation of c w reflects both inertial filtering of local flow and non-local effects such as the caustics [40,41] or sling effect [31].
The results show that gravity plays a very important role in determining c w . This is not obvious; as for particles sedimenting in a still fluid, the monodisperse particle RRV would be zero for all r . While gravity does not alter much the scaling exponent for St < 0.4 when compared to non-sedimenting cases, it drastically affects c w for large droplets, with the scaling exponent maintaining a plateau, in sharp contrast to the case of non-sedimenting droplets where the exponent decreases monotonically with the particle radius. The plateau is reached for smaller droplet sizes when the flow Reynolds number is smaller. We also note a clear dependence of c w on R λ . At a given St, the scaling exponent for sedimenting droplets decreases with increasing R λ . In contrast, the scaling exponent for non-sedimenting droplets increases with increasing R λ . For non-sedimenting droplets of a given St, the dominant eddy size contributing to RRV is in the dissipation subrange and increasing R λ expands the range of larger scales of fluid motion that could partially contribute to the RRV. On the other hand, for sedimenting droplets, the contribution to RRV by the dissipation-subrange eddies is reduced by gravity (see also figure 13) and this relative reduction is likely to increase with R λ . This qualitatively explains the different  dependence of the scaling exponent c w on R λ . Further analysis is needed to fully understand this complex effect of gravity.
In figure 12 we compare the power-law exponent c w from our simulations with the results from Bec et al [42] for non-sedimenting droplets. We note that Bec et al considered a suspension of small, heavy and dilute particles without the gravity. Our results are in very good agreement with their results. They also found that the scaling exponent increases with flow Reynolds number.

4.3.2.
Radial relative velocity at contact. The monodisperse RRV of nearly touching particles as a function of their radius are shown in figure 13. Theoretical prediction of the RRV (black dashed lines) was taken from the formulation of Ayala et al in [38]. An alternative analytical study of the RRV but for non-settling particles has been recently developed by Pan and Padoan [43].
Several observations can be made from figure 13. Firstly, the RRV is insensitive to the details of the forcing scheme and the results from the two different forcing schemes are quantitatively similar. Secondly, the RRV for non-sedimenting droplets increases monotonically and roughly exponentially with particle radius, due to increasing contributions from large-scale turbulent motion and increasing non-local contributions (e.g. caustics produced by the sling effect [31]) when the Stokes number is increased. In addition, for particles larger than 30 µm, the RRV clearly depends on the flow Reynolds number, with smaller RRV corresponding to larger R λ .  [38]. The thin blue line shows the RRV for fluid elements. G indicates simulations with gravity.
The RRV from sedimenting droplets is similar to that of non-sedimenting droplets when a < 30 µm. The possible influence of large-scale fluid motion may be gauged by two scale ratios based on the large-scale flow as [3] τ p T e = St For the physical dissipation rate we considered ( c = 400 cm 2 s −3 ), τ p /T e ≈ 0.022 and v p /u ≈ 0.75 for 30 µm droplets at R λ = 100. As far as large-scale fluid motion is concerned, droplets of a < 30 µm are essentially passive (τ p T e ). The settling velocity is less than u so it does not alter the interaction time of droplets with these large eddies. Therefore, the large-scale fluid motion as well as the large-scale forcing should have a negligible effect on the RRV of droplets. We may refer to the interaction of droplets of a < 30 µm with large eddies as gravityindependent interaction, and droplets of a > 30 µm as gravity-modulated interaction.
On the other hand, droplets larger than 30 µm could have v p > u , as these cross through a large-scale eddy quickly. In this case, gravity plays a dominant role in determining the droplet-eddy interaction time. Indeed, figure 13 shows that the RRV for sedimenting droplets is significantly smaller than that of non-sedimenting particles when a > 30 µm. The theoretical prediction of Ayala et al [38], also shown in the figure, captures qualitatively this reduction. The gravity may also reduce the level of non-local contributions, but not completely. These effects cause the RRV for sedimenting droplets to grow at a much slower rate for 20 µm < a < 40 µm. For even larger droplets, the gravity gradually diminishes the effect of turbulent motion due to shorter droplet-eddy interaction time, leading to slowly decreasing RRV with increasing droplet size. It is not clear what is the precise effect of gravity-enhanced clustering for larger droplets seen in figures 8 and 9 on the RRV.
Another interesting observation for sedimenting droplets concerns the dependence of the RRV on R λ . For larger droplets, RRV increases with the flow Reynolds number, which is different from the dependence on flow Reynolds number for non-sedimenting particles. This difference between sedimenting droplets and non-sedimenting particles can be qualitatively explained as follows. In the case of non-sedimenting particles, increasing R λ leads to a smaller τ p /T e (see equation (14)). Since the contribution to the relative motion from large-scale fluid motion decreases with decreasing τ p /T e , it follows that increasing R λ yields a decreasing RRV. On the other hand, for sedimenting droplets, once v p > u , the large eddy-droplet interaction time is no longer T e , but rather is replaced by L f /v p ∼ T e u /v p . A typical response of a particle with an initial velocity v 0 to an eddy of given velocity u e is Therefore, qualitatively, the contribution of large-eddy motion is gauged by the factor In this case, increasing R λ will increase the factor in equation (16), leading to a larger RRV. The main conclusion here is that gravity alters the qualitative dependence on R λ when compared to the case of non-sedimenting droplets. For any given droplet size, there is a tendency for the RRV to saturate, namely, the RRV eventually becomes insensitive to flow Reynolds number. For small cloud droplets (roughly a < 30 µm), this saturation is reached as τ p /T e → 0 according to equation (14). For larger cloud droplets (roughly a > 30 µm), this saturation is reached by the factor in equation (16) approaching one. Regardless, this saturation for smaller droplets (meaning smaller St and S v ) is reached at smaller flow Reynolds number, since the range of flow scales affecting the pair statistics is more limited for smaller droplets due to their smaller Stokes number. This can also be clearly seen from the qualitative arguments expressed by equations (14) and (16).
The above theoretical arguments have two additional implications. Firstly, the Reynolds number dependence discussed above is a secondary effect. This is because for the gravityindependent interaction, τ p /T e 1 so the large-scale contribution is small. For the gravitymodulated interaction, the factor in equation (16) is in general close to one for R > 100. Secondly, the dividing droplet size a ≈ 30 µm here will change with the physical flow dissipation rate since St ∼ 0.5 and S v ∼ −0.25 (see [44]). Assuming the transition from gravityindependent interaction to gravity-modulated interaction takes place at v p /u ∼ 1, we then expect a larger dissipation will delay the transition from inertia-dominated interaction to gravitydominated interaction to a larger droplet size. Larger dissipation increases τ p /T e as well as St S v (∼ 0.25 ), thus will delay the saturation of RRV with R , to a larger R λ .

Kinematic and dynamic collision kernels
The collision rate between droplets can be described in terms of a dynamics collision kernel D 11 , namely, the ratio of collision rate to particle pair concentration. This dynamic kernel can be obtained in the simulation by detecting directly all collision events at each time step and then averaging over time [5]. The relative statistical uncertainty of the dynamic collision kernel, D 11 / D 11 , can be evaluated at the post-processing stage as [8,45] where T is the total time duration used for averaging, V B is the volume of the computational domain, N pairs = N p (N p − 1)/2 and N p is the total number of particles in the system. Alternatively, the collision kernel can be evaluated kinematically as [1] K 11 = 2π R 2 |w r |(r = R) g 11 (r = R), using the data described in sections 4.2.2 and 4.3.2.
In figure 14 we plot both the kinematic and dynamic collision kernels computed in our simulations at different Reynolds numbers. We confirm that, in all cases, results for the kinematic kernel (dashed lines) are in excellent agreement with those of the dynamic kernel (solid lines). Since increasing flow Reynolds number increases both the RDF (figure 8) and the RRV (figure 13), it has the effect of increasing the collision kernel. The deterministic forcing yields a slightly larger collision kernel, but the difference between results employing two different forcing schemes is negligible for droplets smaller than 30 µm. For droplets larger than 30 µm, the relative difference does not exceed 20%-it should be noted that part of this difference is due to the larger flow Reynolds number using the deterministic forcing.
The key finding concerns the influence of gravity on the collision kernel. We conclude that gravity reduces the collision kernel by a factor of 2-3. This can be attributed mainly to the smaller relative velocity between sedimenting particles. A secondary effect of gravity is to modify the RDF for both intermediate and large droplets. However, the effect of gravity on the collision kernel becomes negligible for droplets less than 30 µm.
To be complete, we also tabulate the dynamic and kinematic kernel data in tables 3-7, which extends our previous results reported in [5]. These data could be used to guide the development of theoretical parameterization of the turbulent geometric collision kernel for cloud droplets, especially the effects of gravity and flow Reynolds number on the collision statistics.

Dependence on R λ
Finally, we discuss briefly the effect of flow Reynolds number on the kinematic collision statistics. In figure 15, the two kinematic quantities (RDF and RRV) for droplets of intermediate sizes (20 and 30 µm) are shown as a function of R λ . The RRV is almost independent of the flow Reynolds number. The RDF increases with R λ for R < 100, but then appears to saturate for R > 100. In the case of deterministic forcing, there is even some evidence of realizing a peak value when 100 < R λ < 200, followed by a slight decrease in the RDF value as R λ is further increased. This non-monotonic behavior is not observed when the stochastic forcing scheme is applied. Our results extend the range of R λ studied in [7]. For larger droplets, the dependence on R λ is illustrated in figure 8 for the RDF and in figure 13 for the RRV. We note that previously, this saturation with flow Reynolds number has been shown for non-sedimenting particles [22]. For the smaller droplets of a = 20 µm, there is no significant difference between results (both the RDF and RRV) using different forcing schemes. This can be understood by the fact that the relative motion and distribution of small droplets are not governed by the large scales where the large-scale forcing could have a direct impact, as pointed out in section 4.3.2. For the 30 µm droplets, the deterministic forcing yields a larger RDF as a result of more coherent vortical structure and larger realized R λ .

Summary and conclusions
In this paper, we reported on our efforts in developing high-resolution simulations of the turbulent collision of cloud droplets. In order to perform high-resolution simulations, we have developed a parallel implementation of our hybrid DNS, based on 2D DD. This new development enables us to conduct hybrid DNS with a flow field solved at grid resolutions up to 1024 3 while simultaneously tracking up to ∼10 8 aerodynamically interacting droplets, although we limit the discussions in this paper to geometric collision statistics of monodisperse particles only.  developed code, and to extend our understanding of kinematic and dynamics collision statistics of cloud droplets. Three specific but inter-related questions have been addressed in a systematic manner. The first question concerns the effect of the large-scale forcing scheme. This was motivated by the use of both deterministic and stochastic forcing methods in previous studies which led to sometimes different quantitative collision statistics [5]. We employed both forcing schemes in a single code to test the sensitivity of the simulation results on the large-scale driving mechanism. The often used working assumption in DNS of the turbulent collision of cloud droplets is that the relative motion of droplets, due to their small size (a η) is small inertia (St ∼ 1 or less), is governed primarily by small-scale eddies of turbulence which are adequately resolved in DNS. This assumption may be questionable for larger droplets in low Reynolds number DNS turbulence as larger droplets could respond to a range of flow scales. We found that, in general, the results using the two forcing schemes are quantitatively similar, with the deterministic forcing giving a slightly larger RDF. The collision kernel is up to 20% larger (figure 14) for the largest droplets (a = 60 µm) considered in this study and 15% for droplets with radii a = 50 µm-it should be noted that part of this difference is due to the larger flow Reynolds number associated with the deterministic forcing. The difference is negligible for a < 30 µm.
A closely related question concerns the effect of flow Reynolds number or equivalently the range of flow scales represented in DNS. The high-resolution simulations provided a range of flow Reynolds number R λ , making it possible to study the dependence of pair statistics of droplets on R λ . We have shown that the Reynolds number dependence due to droplet interaction with large-scale fluid motion is of secondary importance. The interaction will eventually saturate, leading to R λ -independent results at some large R λ . We have shown by both DNS results and by theoretical arguments that the saturation happens at a smaller R λ for smaller droplets. It must be noted that there is another secondary effect of flow Reynolds number concerning the small-scale intermittent distribution of local dissipation rate, which has been discussed in [1,46].
Since most previous studies of the turbulent collision of inertial particles concerned nonsedimenting particles (i.e. gravity was not considered), we have specifically addressed the role of gravity on collision statistics, by simultaneously simulating collision statistics with and without gravity. We have shown (figures 8, 13 and 14) and also argued theoretically that the collision statistics is not affected by gravity when a < a c . For larger droplets, gravity alters the  particle-eddy interaction time and significantly reduces the RRV. The critical droplet radius a c is found to be around 30 µm for the RRV, and around 20 µm for the RDF. This critical size is expected to increase with the flow dissipation rate. We have proposed a concept of gravity-independent interaction versus gravity-modulated interaction to address the secondary effect of large-scale motion. The effect of gravity on the RDF is rather complex: gravity reduces the RDF for intermediate-sized droplets but enhances the RDF for larger droplets. This could be a combined result of reduced interaction time due to gravity and inertia-induced preferential sweeping [3].
In addition, we have also studied the scaling exponents of both RDF and RRV. For nonsedimenting particles, our results of scaling exponents are in excellent agreement with published results from other studies. We found that gravity modifies the RDF scaling exponents for both intermediate-sized and large particles (figures 6 and 7), in a manner very similar to the effect of gravity on the RDF at contact (figure 8). Gravity is shown to cause the scaling exponents, c g and c w , to level off for large droplets, in contrast with diminishing exponents for non-sedimenting particles. Further work is needed to understand such results and to develop a parameterization of RDF that includes the gravity effect. We have tabulated in tables 3-7 all results so that they could be used to guide the development of theoretical formulations to better understand the effects of gravity and flow Reynolds number. Greek letters β fitting coefficient D 11 dynamic collision kernel (cm 3 s −1 )

Appendix B. Parallel performance
The scalability performance of the six major tasks in 1D and 2D DD implementations is shown in figures B.1(a) and B.1(b), respectively. The timing measurements have been carried out for a benchmark problem of 2 × 10 6 droplets of radii 20 and 40 µm at 512 3 grid resolution. The energy dissipation rate of the physical flow was set to c = 400 cm 2 s −3 . All simulations were performed on NCAR's Bluefire machine (IBM Power 575 cluster) with identical compiler and optimization settings. Measurements are averaged over 1000 time steps of simulation. For this particular problem, the maximum number of processes that could be used in 1D DD implementation is 64. This limitation results from a shortcoming of the algorithm for treating aerodynamic interaction, namely, the distance between interacting particles H trunc cannot extend beyond the adjacent subdomains. As the number of processes increases, the thickness of the subdomains becomes very thin. Consequently, not all particle pairs separated at distance H trunc < 50 are realized by the algorithm. There is a less strict limitation in the 2D DD code as the thickness of the subdomains (columns) decreases roughly as N / √ P, when compared to N /P in 1D DD. Here P = P y P z is the total number of processes used.
For the implementation based on the 1D decomposition, all simulation tasks have an almost ideal scaling with the number of cores for a small number of processors, but quickly start to saturate at about N proc = 64. This suggests that 1D decomposition will not be the method of choice for larger problem sizes. The saturation is due to a systematic increase in the communication time as a result of the increase in the number of processes. The ratio of the communication time to the computation time scales as the ratio of surface area of a subdomain to its volume, which for the 1D DD code is proportional to ∼2P/N , where N is the grid resolution. For N = 512 and P = 64, this ratio is 0.25. The data size to be communicated is proportional to the surface area, which is ∼N 2 .
On the other hand, the 2D decomposition does not show such signs of apparent scalability saturation, for up to 1024 cores tested. This suggests that the 2D implementation has better scalability. This can be understood by the ratio of the surface area of the subdomain to its volume for 2D decomposition, which is 2(P y + P z )/N or roughly 4 √ P/N . This ratio is a factor of 2/ √ P smaller than that of 1D decomposition. For N = 512 and P = 1024, this ratio is 0.25. Furthermore, the data size for each communication is proportional to the surface area, which is N 2 / √ P. Both the smaller surface to volume ratio and the smaller data size make the 2D decomposition much more scalable than the 1D decomposition.
We note that the 2D DD code could not be run on a very small number of cores due to memory constraints. In contrast to 1D decomposition, the flow simulation in 2D decomposition is the bottleneck of the overall HDNS simulation. This is primarily due to the high computation and communication costs of 3D FFT.