## Abstract

The purpose of this paper is to quantify the influence of turbulence in collision statistics by separately studying the impacts of computational domain sizes, eddy dissipation rates (EDRs), and droplet sizes and eventually to develop an accurate parameterization of collision kernels. Direct numerical simulations (DNS) were performed with a relatively wide range of EDRs and Taylor microscale Reynolds numbers . EDR measures the turbulence intensity levels. DNS model studies have simulated homogeneous turbulence in a small domain in the cloud’s adiabatic core. Clouds clearly have much larger scales than current DNS can simulate. For this reason, it is emphasized that obtained from current DNS is fundamentally only a measure of the computational domain size for a given EDR and cannot completely describe the physical properties of cloud turbulence. Results show that the collision statistics are independent of the domain sizes and hence of the computational for droplet sizes no bigger than 25 *μ*m as long as the droplet separation distance, which is on the order of the Kolmogorov scale in real clouds, is resolved. Instead, they are found to be highly correlated with EDRs and droplet sizes, and this correlation is used to formulate an improved parameterization scheme. The new scheme well represents the turbulent geometric collision kernel with a relative uncertainty of 14%. A comparison between different parameterizations is made, and the formulas proposed here are shown to improve the fit to the collision statistics.

## 1. Introduction

High-resolution simulations and in situ measurements (such as observations by aircrafts) of cloud properties continue to advance our understanding of cloud and precipitation processes. However, difficulties persist in resolving certain puzzles in warm rain. Classical condensation theory stipulates that the growth rate of a cloud droplet radius is inversely proportional to the radius. As a result, a rain drop cannot be formed by condensation in a realistic time scale. Gravitational collision can speed up the growth process, but it is inefficient until a droplet attains a size of about 30 *μ*m (Pruppacher and Klett 1997). Classical theory, therefore, is unable to account for the observed rapid onset of warm rain (Szumowski et al. 1997; Knight et al. 2002). Over the past 2–3 decades, a substantial amount of work has been conducted on studying the effect of turbulence on the generation of warm rain [see the reviews by Devenish et al. (2012) and Grabowski and Wang (2013)]. Studies from laboratories (e.g., Fessler et al. 1994) and direct numerical simulations (DNS) (e.g., Squires and Eaton 1991; Ayala et al. 2008a) have demonstrated that air turbulence can accelerate droplet collisional growth by increasing the relative velocity and enhancing coagulation effects. The large-eddy simulation (LES) studies (e.g., Seifert et al. 2010; Wyszogrodzki et al. 2013; Franklin 2014; Grabowski et al. 2015) showed that turbulence can accelerate the onset of warm rain in shallow convective clouds as well as increase its intensity. In this paper, our focus is on the effects of turbulence on droplet growth by collision–coalescence. These effects are particularly important in the early stage of droplet growth and may be crucial for warm rain initiation.

It is known that there are three major effects of turbulence in enhancing the chance of droplet collision. The first one is the turbulent transport effect, which increases the droplet radial relative velocity (RRV) via the local shear and air acceleration. The second effect, termed the clustering effect, redistributes the cloud droplets such that they cluster in regions of low vorticity and high shear because of their inertia (Maxey 1987; Wang and Maxey 1993; Reade and Collins 2000; Vaillancourt et al. 2002). Chen et al. (2006) proposed that clustering can also occur in regions of low Lagrangian acceleration. The third effect is related to droplet–droplet aerodynamic interactions, which affect the collision efficiency. Because of the complexity of the problem and the lack of accurate and consistent representation of the disturbance flow in a turbulent background (Rosa et al. 2011), there are only a few DNS studies concerning this last effect (Wang et al. 2005; Ayala et al. 2007, 2014). Most of these treatments are still under development and remain computationally expensive.

To assess the effects of turbulence, there is a need to define the collision kernel, the relative velocity between droplets, and the degree of clustering. In a stagnant flow, the droplet collision rate is customarily described by the gravitational collision kernel: , where *R* is the sum of the radii of droplets involved in collisions (, and are the terminal velocities, and is the collision efficiency. Since we are only concerned about the geometric collision kernel (collision without considering droplet–droplet interactions), will not be pursued here.

Saffman and Turner (1956) introduced two expressions for the kinematic collision kernel by using cylindrical geometry and spherical geometry, respectively. The kernel denotes the surface of a droplet’s swept-out volume within which collisions are expected to occur [see Wang et al. (1998) for details]. Sundaram and Collins (1997) proposed a turbulent collision kernel as the product of Saffman and Turner’s (1956) cylindrical form and the radial distribution function (RDF) *g*(*R*):

Here, is the relative velocity between two droplets when they collide; *g*(*R*) measures the clustering effect of turbulence, with values larger than unity if clustering is present. In stagnant air, is just , and *g*(*R*) = 1 when droplets are uniformly distributed. Wang et al. (1998) extended this turbulent collision kernel to a spherical form and demonstrated that this formulation is more appropriate for the turbulent case since the cylindrical form overpredicts the collision kernel by about 20% or more:

Here, is the RRV, defined as , with being the separation vector.

It is widely accepted that the overall turbulent enhancement of collision statistics (i.e., the RRV, the RDF, and the collision kernel) depends on the turbulent eddy dissipation rate (EDR) *ε* and the droplet Stokes number (St). The latter, defined as the ratio of the particle response time to the Kolmogorov time scale , quantifies how fast a droplet reacts to the turbulent flow.

There are also a significant number of DNS studies demonstrating that the collision statistics are sensitive to the Taylor microscale Reynolds number . However, one should be very careful when interpreting their results, as there are two types of we use in the cloud community: one is real, and the other is computational. By definition, depends not only on the Taylor microscale *λ* but also on the turbulent fluctuating velocity at the largest scales. For the real , the largest scale is the overall size of the entire cloud; this , usually with an order of approximately 10^{4} (Pinsky et al. 2006), represents a physical property of the cloud turbulence. Clearly, clouds have much larger scales than current DNS models can simulate. We, and previous authors, have merely simulated homogeneous turbulence in a small computational box, presumably located somewhere in the adiabatic core of the cloud. For this reason, obtained from these models refers only to the computational domain size *L*; this computational , determined by the size of the domain, currently only attains values around 600 or even less. Franklin et al. (2005) found collision statistics increased with both and EDR from their DNS. Based on these simulations, Franklin et al. (2007) developed an empirical parameterization for the collision statistics that scales with . EDR terms were excluded from the parameterization by replacing them with using an empirical power law fitted to their data [(12) in Franklin et al. (2007)]. Since this power law is only valid in their particular model configuration and cannot be applied to the general case, the replacement may not be justifiable. In their model configuration, and EDR increased simultaneously, and thus the separate effects given by EDR and were unknown. Specifically, their *L* was fixed to for all simulations, and kinetic energy was pumped into each simulation at a different rate to create different turbulent intensities. Accordingly, spatial resolution was modulated in such a way that the turbulent dissipation range (scales below the Kolmogorov length scale) would be well resolved. However since EDR is negatively correlated to *η*, the scale separation between *η* and *L* depends on the turbulent intensities. In this case, is positively correlated with EDR such that the influence of the two cannot be separated. It follows that the cause of increments to the statistics was unclear. Therefore, a controlled study that can separate *ε* from is called for.

Ayala et al. (2008a) found that the effect of is secondary after reaching a certain value, and the EDR effect is dominant. Rosa et al. (2012, 2013) also found that collision statistics scaled with but eventually reached an asymptotic value after exceeded some threshold that depends on the size of the droplets involved in the collision. Collision statistics of small droplets () have smaller threshold values than large ones (), since the latter have larger Stokes numbers and interact with a wider range of turbulent eddies. However, Rosa et al. (2012) also commented that this dependency was a result of “inadequate flow scale separation” in low-resolution simulations, because it only appeared at low . For this reason, caution should be exercised when interpreting the dependency, as the integral scale of the turbulence that determines the real Reynolds number of the cloud cannot be resolved by current DNS. In addition, previous DNS studies were either limited to a narrow range of domain sizes or to only a few dissipation rates (see Table 1). There is a need to extend these results to a wider range of parameters.

In the past few decades, there were several attempts to parameterize the collision kernel, but all show limitations [see summary by Ayala et al. (2008b)]. They either excluded gravitational effects or clustering effects, or the results were only applicable to very large or small droplets or to a small range of . As discussed earlier, the scaling in Franklin et al.’s (2007) empirical parameterization is inconclusive because the effect of EDR cannot be isolated. In addition, their parameterization is only applied to the cross-sized collision case, as same-sized collisions (collisions between droplets with the same size) were not included in their model. Ayala et al. (2008b) derived analytical formulas for RRV and RDF that also depend on . They extended the range beyond the DNS limitation by taking into consideration the general representation of fluid velocity correlations in their formula. However, the integral length scale that determines the cloud-scale Reynolds number cannot be resolved by DNS; their parameterization might therefore be inaccurate. The parameterizations from Wang et al. (2000) and Zhou et al. (2001) also depend on but were developed based on low-resolution runs that might not provide sufficient scale separation. In Wang et al. (2000), ranges from 24 to 75, and Zhou et al. (2001) only studied the case of = 45 and 58. Wang et al. (1998) proposed a scaling law for the collision kernel and stated that it is independent of when is sufficiently large, as in the case of real clouds.

This paper represents our continued exploration of cloud droplet collisions under turbulence begun by Franklin et al. (2005, 2007). Our purposes are 1) to better understand the influence of on the collision statistics by separating its effect from that of the EDR and to clarify the physical meaning of in DNS; 2) to quantify the influence of turbulence intensity on cloud droplet collisions by simulating turbulent flow fields over a broad range of EDRs; and 3) to seek a better parameterization for collision statistics with respect to different turbulent conditions. The paper is organized as follows. Section 2 describes the DNS model and the methods used; section 3 provides the results, showing that the collision statistics are independent of while monotonically increasing with EDR. Based on these results, we propose a new parameterization of the collision statistics in section 4, and section 5 forms the summary and conclusions.

## 2. Model and methodology

Our model mainly follows the framework of Franklin et al. (2005) and Vaillancourt et al. (2001). A brief description is given below. Interested readers may refer to the original papers for details. Improvements and comparisons between our model configuration and that of Franklin et al. (2005) will be given in section 2b.

### a. Flow, particles, and numerics

The turbulent flow field was generated by solving the incompressible Navier–Stokes equations using a triply periodic pseudospectral technique (Orszag 1969). The governing equations are

where is the flow velocity, *P* is the pressure, , and is an external forcing to maintain the turbulence. The turbulent flow was assumed to be statistically homogeneous and isotropic, as is the characteristic of the adiabatic cores of cumulus clouds (Vaillancourt and Yau 2000). The adiabatic cores are regions free of entrainment and mixing with higher liquid water content than the rest of the cloud. Therefore, large droplets are expected to be found in these regions to initiate the collision and coalescence processes. To maintain statistical stationarity, the flow was forced in a low-wavenumber band following the approach of Sullivan et al. (1994) with , where is the wave vector. This forcing was improved in order to maintain a more stable average EDR, and details will be given in section 2b.

As for the cloud droplets, radii in the range from 5 to 25*μ*m were chosen, as these sizes are critical to forming larger droplets to set the stage for effective collisional growth. For sizes smaller than 30 *μ*m, Stokes’s drag law can be applied (Pruppacher and Klett 1997). Since we assume no droplet–droplet aerodynamic interactions, the disturbance velocity field caused by the droplets was ignored. Thus, the equation of droplet motion took on the following simple form:

where denotes the droplet velocity; is the flow velocity at the droplet center position; is the droplet response time, with *r* as the droplet radius; and and are the density of air and droplet, respectively. We assume that is constant with a value of . Therefore, the droplet terminal velocities are calculated as for small droplets.

To make our simulation close to reality, the initial droplet size distribution (DSD) of our polydispersive system was based on measurements made by Squires (1958) for trade wind cumulus clouds off the coast of Hawaii [refer to Fig. 5.9 in Rogers and Yau (1989)]. In the simulations, droplets were randomly placed throughout the domain after the turbulence had reached statistical stationarity. At each droplet position, the initial droplet velocity was given by the flow velocity obtained by linear interpolation (Yeung and Pope 1989). To ensure that the initial conditions do not influence the results, droplets were given enough time (on the order of for the largest droplet size) to respond to the flow before compiling the collision statistics.

We used the cell index method (Allen et al. 1987) to track droplet motions, following Franklin et al. (2005). The droplets were treated as ghost particles, meaning they were allowed to pass through each other after the collision and continue their original course. This algorithm has been validated in detail in Franklin et al. (2005) and is shown to be consistent with the Saffman and Turner (1956) formulation. Collisions were identified and counted by solving the quadratic equation [(9)] in Franklin et al. (2005).

Three statistics were examined to measure the turbulent effects on droplet collisions: the collision kernel, the RDF, and the RRV. The definition of RRV was given in section 1, and the RDF is defined as follows:

where is the domain size, and are the number of droplets from the two different size groups in the domain (for same-sized collision, becomes ), is the number of time steps, and is the number of pairs at contact found in the shell volume within . The shell volume had an outer radius slightly larger than *R* to allow enough sampling of droplet pairs at finite . It is safe to do so, because very little dependence of the RDF on the shell size is found as long as the outer radius is smaller than 1.1*R* (Wang et al. 1998).

In addition to its kinematic equation shown in (2), the collision kernel can also be computed by dynamically counting the number of collisions that take place during the simulation:

where is the number of collisions, and is the time step.

### b. Model improvements and validation

#### 1) Expansion of the domain size: MPI technique and validation

To run simulations with larger box sizes, we used the message passing interface (MPI) for parallelization instead of the OpenMP interface in the code of Franklin et al. (2005). The drawback of OpenMP is that it requires shared-memory machines with a small number of processors. On the other hand, the MPI technique can be used in both shared- and distributed-memory architectures and therefore allows a much larger number of processors. We used a one-dimensional decomposition technique to divide the computational domain equally in the *y* direction into vertical slabs such that each processor works on one slab. With this technique, we performed simulations with up to 589. Local dissipation rates of shallow cumulus clouds in the incipient or decaying region are usually moderate, with values around 10–100 cm^{2} s^{−3} (Seifert et al. 2010). Numerical studies show that the highest local dissipation rates can sometimes reach a few thousand at some regions of clouds (Franklin 2014; Benmoshe and Khain 2014). Therefore, the EDR used here spanned the range from 50 to 1500 cm^{2} s^{−3}.

To validate our MPI-based code, two tests have been performed with the following configuration. We fixed the box size to in the *x*, *y*, and *z* directions. We placed 25 000 droplets with a radius of 10 *μ*m and the same number of droplets of 20 *μ*m were placed randomly into the domain. The liquid water content was approximately , which is close to the typical value found in adiabatic cloud cores.

The first test was to validate the droplet tracking and collision detection scheme. To simplify the case, turbulence was turned off so that droplet motion was only driven by gravity along the *z* axis. The number of grid points was set to . Since the spatial distribution of droplets in the *x*–*y* plane does not change in time without turbulence, the initial condition for a single realization strongly affects our statistics. To avoid this problem, 20 realizations with different random initial droplet positions were used, and collision statistics were calculated from the ensemble average. We obtained an average RDF of and an average collision kernel of . These values are close to the respective theoretical values [ for a uniformly distributed case and according to (2)], with only of difference.

Restoring the turbulence, the second test was to find a justifiable value of , where is the maximum resolvable wavenumber. To resolve the smallest scales of turbulence, while at the same time avoiding unnecessary computations on noncontributing scales far below the Kolmogorov length scale, we should minimize , while keeping it larger than unity. We performed 15 different simulations with different numbers of grid points varying from up to and with a fixed EDR of 1000 cm^{2} s^{−3}. In so doing, the turbulent flow parameters (such as *u*′, , and *η*) remained identical, and the grid size became the only factor influencing our collision statistics. With the grid size changing, spanned a range from 0.2 to 3.4, covering most values used by other researchers (Ayala et al. 2008a; Franklin et al. 2005; Wang and Maxey 1993; Yeung and Pope 1989). Figure 1 illustrates the dependence of the collision kernel and RDF on . All values were normalized by their respective averages over the data points with . As shown in the figure, the collision statistics varied within ±2% from the mean after . To be on the safe side, we settled on a value of 1.3 in all our subsequent simulations. Note that is determined solely by , since in our dealiased pseudospectral code, and *η* solely by the EDR, since . It follows that a one-to-one correspondence can be established between and *ε* once we fix the value of .

To test the scalability of the code, we performed seven groups of simulations with different numbers of processors (2, 4, 8, 16, 32, 64, and 128) on two different machines: the Silicon Graphics International Corporation (SGI) Altix 4700 supercomputer at the University of Montreal and the Mammouth parallèle II (MP2) machine at the University of Sherbrooke. Each group contained three members, and we took the average of the three as a reference wall time. All tests were performed at with , and we injected a total of 424 448 droplets into the flow, half of which have a radius of 10 *μ*m and the other half with 20 *μ*m. The code scales reasonably well on both machines. According to Amdahl’s law, the performance of a parallelized program is constrained by its sequential part when the number of processors increases (Amdahl 1967). If the code is fully parallelized, the true wall time will be half when the number of processors doubles. In other words, their relation can be precisely described by a power-law formula with the power equal to −1: . Here, is the wall time when the program is executed with a single processor. The wall times of our executions spanned from 60 min (with 128 processors) up to 2000 min (with 2 processors) on the SGI Altix, fitted into the curve of . The experiments on the MP2 gave a fitted curve of . Even though we did not attain perfect scalability on both machines, we saw a significant speedup of the simulations as the number of processors doubled.

#### 2) Improvements in computing statistics of flow parameters

In this paper, we used different approaches to obtain the EDR, the Taylor microscale *λ*, and the of the flow field. Franklin et al. (2005) estimated EDR from the energy spectrum, while we calculated it from the numerical integration of enstrophy (, where *Z* is the total enstrophy), which provides more accurate values. We calculated *λ* via in Franklin et al. (2005), where KE is the total kinetic energy in the system, and our *λ* was from . Here . In our study was computed via , while in Franklin et al. (2005). Of course, these quantities respect the same scaling in Kolmogorov phenomenology but not with the same proportionality constants.

To quantify the difference in the resulting flow statistics, we reran the four cases in Franklin et al. (2005). The results are listed in Table 2, along with the Franklin et al. (2005) results for comparison. The number of grid points is represented by *N*. Flow statistics show that *λ* in our study was around 3.0 times the value in Franklin et al. (2005) and was around 1.85 times. Since we did not change the actual properties of the flow but only the method of evaluating these statistics, similar collision statistics were expected under the same particle tracking and collision detection techniques. Indeed, the collision statistics in Table 2 show that our results were indeed close to those of Franklin et al. (2005), within a difference of .

#### 3) Separating the effects of and EDR

As explained in section 1, the conclusions regarding from Franklin et al. (2005) need further exploration, since the effects of EDR and cannot be separated. To solve this problem, we fix when simulating the same EDR, since EDR is determined by alone when is fixed [see section 2b(1)]. To increase , we increase the number of grid points *N* and hence the computational domain size *L*, thereby resolving larger scales of the flow while keeping EDR fixed. Since in our study is linearly proportional to *η* and scales with domain size *L*, it follows that . Therefore, when varying at fixed EDR, we changed *N* and fixed . Likewise, when varying EDR to study its influence, we tuned to obtain our fixed value and fixed *N*.

#### 4) Two forcing schemes

The total kinetic energy in the forcing band is fixed in Franklin et al. (2005). We refer to this forcing scheme as forcing 1. However, not knowing how much kinetic energy should be fixed to reach a desired average EDR, we have to manually adjust the initial energy in the forced wave band by trial and error, which is very inefficient. For this reason, we propose a modified version of the forcing (termed forcing 2), which injects equal amounts of energy into the forcing band at each time step. The energy supply at the forcing scales is balanced by dissipation at the viscous scales, when turbulence reaches statistical stationarity. It follows that the rate of energy injection can easily be determined by the EDR at the initial time step and held constant thereafter. In this way, we decrease the number of computations by specifying a steady average EDR directly.

To test the accuracy and efficiency of forcing 2, we reran the simulation of Franklin et al. (2005) at with this forcing and compared the results with those from forcing 1. The EDR for forcing 1 was 155 cm^{2} s^{−3} (see Table 2). We used this value to calculate the energy injection rate in the case of forcing 2. Since the turbulence is both statistically stationary and ergodic, ensemble averages can be represented by time averaging. Figure 2a demonstrates that forcing 2 reproduced a very similar ensemble kinetic energy spectrum as forcing 1 with a negligible difference of about 2% in flow properties. This result demonstrates that forcing 2 produces reliable results that are statistically identical to forcing 1 but arrives at the kinetic energy spectrum, obtained by specifying the EDR, much more readily. To examine further the efficiency and accuracy of forcing 2, we did another four simulations with . The computational box size was fixed, but resolution was different in each case. Figure 2b illustrates that the ensemble kinetic energy spectra from different resolutions lay on top of each other without any discernible shifts. This result lends support to the capability of forcing 2 to produce precise mean EDR covering a range of spatial resolutions, since DNS average EDRs vary within only . Although forcing 2 is more efficient in controlling the average EDR than forcing 1, most of our computations were obtained using forcing 1. To avoid repeating the costly calculations, our statistics will be analyzed from these data, bearing in mind that the EDRs show scatter within from the desired value, but there were only two cases with scatter beyond . Forcing 2 was used for the validation tests only.

## 3. Results

To test the dependency and EDR dependency, our simulations employed five different eddy dissipation rates (*ε* = 50, 200, 500, 1000, and 1500 cm^{2} s^{−3}), and nine different box sizes (, and ). We managed to run most of the scenarios except those with both high resolution and low dissipation rates (see Fig. 3). The reason is that *N* is largest in the highest-resolution simulation, and is largest at the lowest EDR run [see section 2b(1)]. As a consequence, the slab volume and the number of droplets assigned to each processor become very large, causing an inordinate amount of communication between processors, yielding extremely slow computations. Nevertheless, the successful runs were sufficient for us to quantify the turbulence effects on droplet collisions as they span a broad range of EDR and .

### a. Collision statistics with different box sizes

The in our simulations ranged from 63 () to 589 (), which gave us a good indication of how collision statistics depend on or, equivalently, on the box size. Rosa et al. (2013) tested the convergence of collision statistics of 20–20-*μ*m droplet pairs and 30–30-*μ*m pairs with using two different forcing schemes: one was stochastic forcing, where random acceleration was added to the velocity component in a low-wavenumber band, and the other was deterministic forcing in which the energy level in the two smallest wavenumber bands was fixed. They concluded that the RDFs between large droplets () were harder to converge with than between small droplets () because they were under the influence of larger scales of motion. The growing trend at low in both forcings, as explained in their paper, was because more collision-related scales were resolved when increased until a saturation level was reached. In addition, there was another possible contribution. Although the average dissipation rate was set to 400 cm^{2} s^{−3} in all cases, some variations were observed. The curves of RDFs follow closely the fluctuations of average EDR in both forcing schemes. Under stochastic forcing, the EDR increased slightly with [see section 4.1 and Fig. 2 in Rosa et al. (2013)], while the corresponding RDFs were larger as well (see Fig. 15 in their paper). Under deterministic forcing, the EDR reached its peak around *R*_{λ} = 100–200, and a similar pattern was also found in the RDFs at . According to our DNS, the sensitivity of the RDF of 25–25-*μ*m droplets to EDR [i.e., ] was about 1.42 × 10^{−2} (cm^{2} s^{−3})^{−1} for EDR around 500 cm^{2} s^{−3} [i.e., increasing 1 cm^{2} s^{−3} in EDR will result in an increase of 1.42 × 10^{−2} in *g*(*R*)]. This implies that a change of 100 cm^{2} s^{−3} in EDR will result in a change of 1.42 in the value of RDF, which is on the same order of magnitude as that of the RDF fluctuations in their large regions. In section 3b, we will show in detail that collision statistics are very sensitive to the EDR.

We are in a position to compare the collision statistics with a wider range of under fixed EDRs. Figure 4 demonstrates the collision kernel between 10- and 20-*μ*m droplets with different , with each curve indicating different EDRs when using forcing 1. The curves fluctuate mainly because the mean EDR in each run scattered within 15% from the nominal value shown in the legend. The standard deviation of the statistics increases with , which can be induced by the intensified intermittency of the turbulence as domain size increases. Kolmogorov (1962) conjectured that the variance of local (instantaneous) dissipation rate normalized by its square mean value scales with Reynolds number as follows: , with the exponent confirmed by later experiments and atmospheric observations (e.g., Pope 2000; Siebert et al. 2010). Note that is also the vorticity kurtosis *K* and can be used to quantify the level of intermittency. We calculated *K* with different at the EDR of 1000 cm^{2} s^{−3} and found our exponent consistent with , with a relative error of 2.75%. Though intermittency experienced a power-law growth with the expansion of domain size, the collision kernel shown in Fig. 4 nevertheless is insensitive to , explicitly showing that the effect of intermittency is secondary, at least over the simulated range. This flatness with was also found in the collision kernel of all other droplet pairs, including both same-sized collisions and cross-sized collisions, as shown in Figs. 5a and 5b. The standard deviation was of a similar order of magnitude as in Fig. 4 and therefore is not shown. The fluctuations in the curves were caused by slight variations in the dissipation rate, as mentioned earlier. As can be seen from the gray dashed line in Fig. 5a(ii), the EDR varied around 500 cm^{2} s^{−3} over different simulations, creating a spike around . This, as a consequence, diminished the flatness of these collision statistics. Figure 5b demonstrates the same fact with another dissipation rate. As for RDF and RRV shown in Figs. 5c and 5d, respectively, the Reynolds number did not alter the collision statistics. Specifically, previous studies have shown that clustering happens in low-vorticity regions (e.g., Maxey 1987). Therefore, we found it interesting that the RDF is also independent of domain size, and we conclude that intermittency has a relatively small impact. However, this conclusion has to be interpreted cautiously. It is possible that the intermittency is still too weak, since the here may not be sufficiently large to show the high levels expected in real cumulus clouds. The independency implies that the large-scale motion has little effect on the droplets, and the collision-related scales have been fully resolved at even small (). Droplets of small Stokes number adjust quickly to the local fluid acceleration before colliding with surrounding droplets, and the local influence might greatly dominate over their acceleration history given by the large-scale motions. Thus, the effect, if any, of on droplet collision statistics is secondary compared to the effects from the local dissipation rate and Stokes number. As discussed in section 1, the dependency is, in fact, mainly the result of the model’s inadequate flow scale separation seen only at very low resolutions (i.e., the failure to resolve all the collision-related scales is because the DNS box size is too small). Cloud systems are characterized by an enormous span of scales, and their real Reynolds number is too large to be amenable to DNS. Thus, the calculated from DNS is determined by computational box size *L* rather than any real scale separation in clouds. The same conclusion can apply to , following the scaling . As a result, when interpreting the sensitivity of collision statistics to from DNS data, we should be aware of the other quantities that are changing with the model configuration. We believe that the scale of turbulence affecting droplet collisions is on the order of the mean separation distance of cloud droplets. Generally, cloud droplet number concentration is on the order of 10–1000 cm^{−3}, depending on the cloud type (e.g., Twomey 1959; Bennartz 2007; Leaitch et al. 1992). This gives an average droplet separation distance on the order of 0.1–1 mm, which is close to the order of the Kolmogorov length scale for the typical range of EDRs in clouds. As long as the model resolves scales down to the droplet separation distance, which is far below our smallest computational domain size, becomes irrelevant, and the turbulence can be characterized by EDR alone.

### b. Collision statistics under different eddy dissipation rates

Since we have demonstrated that collision statistics are not sensitive to , we averaged over all resolutions in our studies to examine the effects of EDR. We found that all collision and pair statistics increase monotonically with EDR, which was consistent with the finding of Ayala et al. (2008a). Among the three, RDF seemed to head toward an asymptotic value (Fig. 6b). The fluctuation of statistics among different resolutions grew slightly with increasing EDR, which can be attributed to the increasing intermittency of turbulence (e.g., Pope 2000).

To better evaluate how EDR impacts the different droplet pairs, we normalized the statistics by their respective values in a purely gravitational setting. Since the results are insensitive to , or *L*, only simulations at one resolution were examined. Figures 7a–c showed that the cross-sized collision statistics involving 20-*μ*m droplets increased with EDR. The strongest enhancement in collision kernel occurs in the case of 20–25-*μ*m droplet pairs, with a maximum enhancement of about 3.1 at the highest EDR. This increase was primarily attributed to the stronger clustering at high EDR, especially for the 20–25-*μ*m droplet pairs. Droplet pairs with similar size tended to cluster in the same regions of the flow because of similar droplet inertias and terminal velocities. This has been confirmed by previous studies (e.g., Ayala et al. 2008a; Franklin et al. 2005). An exception to this though, was pairs involving droplets with a radius *r* = 5 *μ*m, which always seem to have a small value of RDF (see Fig. 7b). This was probably because small droplets have small Stokes numbers, indicating that they adjust very quickly to changes in the flow and therefore act more like fluid tracers than inertial droplets. The normalized RRV grew nearly linearly for all sizes with the strongest enhancement for 20–25-*μ*m pairs as well (Fig. 7c). The results involving droplets of other sizes were similar and are thus not shown.

Same-sized collisions are impossible when gravity is the only influence on droplet motion because they have identical terminal velocities. However, when turbulence is present, collisions do occur. Before discussing the results, it should be noted that there was a variation of about 3% because of inadequate sampling, especially for droplets of *r* = 5 and 25 *μ*m. The 5-*μ*m droplets had less chance to collide with one another because of their minute size, and the 25-*μ*m droplets had the lowest concentrations among all the droplets, also leading to low collision rates. We can see from all panels of Fig. 7 that all three statistics increased with the size of droplet pairs, both for same-sized collisions and cross-sized collisions.

## 4. Can we formulate a better collision parameterization?

We arrived at the conclusion that collision statistics for droplets of radius smaller than 25 *μ*m are insensitive to , implying the EDR and the pair sizes are the determining factors. As a consequence, we are of the opinion that parameterizations scaling with are misleading. Here, we introduce a modified parameterization scheme for turbulent collision kernels based on DNS data and at the same time attempt to match the forms of theoretical expressions from Saffman and Turner (1956) and Wang et al. (1998).

### a. Turbulent collision kernel

Given the dependence of the collision kernel on the EDR and droplet size, and the sensitivity test of in section 3, we applied curve fitting to our data to obtain a formula for the geometric collision kernel in the following form:

where and are the Stokes number and response time of the colliding droplets, respectively. The first term on the right-hand side is the gravitational collision kernel, which is . The second term is the turbulence-enhanced component, which is the coupling of the clustering effect and the turbulent transport effect. The collision kernel, as demonstrated by DNS, is a function only of dissipation rate (as *η* and can be derived from it) and droplet size; and are therefore excluded from the formula since they are both determined by the arbitrary size of the computational domain.

The predicted collision kernel calculated from the formula is shown in Fig. 8a along with the results of our DNS. The parameterization reproduced the collision kernel correctly except for the case of same-sized collisions between 5-*μ*m droplets (the points in the lower-left corner of the figure). Those collision events are not the major concern of our parameterization since they contribute very little to the overall droplet growth because of their low chance of colliding compared to other droplet pairs. Otherwise, high accuracy in predicting the collision kernel is achieved by our parameterization. The relative error, which is calculated by normalizing the absolute difference between values from the parameterization and the DNS and taking the average, was used to measure the uncertainty of our parameterization scheme. The uncertainties came from a variety of sources. One was from poor sampling of same-sized collisions, especially for both the largest and smallest droplets, as mentioned in section 3b; another source of uncertainty is simply the scatter resulting from finite-size ensembles of realizations. The result shows that the uncertainty in predicting the cross-sized collision case was very small, with a value of 6%; however, for same-sized collisions, it was still significant (28%).

### b. Radial relative velocity

Wang et al. (1998) modified the analytical equation of RRV based on the studies of Saffman and Turner (1956) and Hu and Mei (1997) and proposed a more accurate formula:

where *Du*/*Dt* is the local fluid acceleration, is the longitudinal Taylor microscale of fluid acceleration defined by and is the longitudinal correlation function of the velocity derivative. We fit curves to match our DNS data on the basis of their formula by keeping the local fluid shear term and gravity term (first and last terms inside the square brackets) and modifying the fluid acceleration term and coupling term (second term and third terms inside the square brackets) in a way that does not depend on [i.e., without and without the term]. The local fluid acceleration and longitudinal Taylor microscale are replaced with expressions involving EDR, which are the only way to use *ε* with and while keeping the same dimensions. We obtained the fit as follows:

The first term inside the curly brackets measures the effect of local shear from turbulence, which is also the only remaining term when one considers only noninertial particles. The second term denotes the different inertial response of the droplet pair to the local fluid acceleration. The third term is the overall response of the droplet pair due to local fluid acceleration, which is critical for same-sized collisions; the last term is the gravity term. Figure 8b demonstrates that the predicted RRV using the parameterization asymptotically approaches the observed value. For cross-sized collisions, the gravity term has the largest contribution to RRV as a result of different terminal velocities; when colliding, droplet sizes get closer until they are equal, the gravitational effect diminishes, and the local acceleration term (the third term) dominates.

Our parameterization scheme for RRV performed with high accuracy for both the cross-sized collision case and the same-sized collision case. For cross-sized collisions, the average root-mean-square error (RMSE) of the normalized RRV (i.e., the RRV normalized by the value in the corresponding stagnant air case) was 0.03, which is much lower than Franklin et al.’s (2007) value of 0.28. The relative error was within 13%, which is slightly higher than that in Ayala et al. (2008b) (5%). For same-sized collisions, the relative error was 15%, which improves on Ayala et al.’s (2008b) parameterization, where the relative error was 48% for . For a closer comparison, our relative error for was 14.73%. Table 3 gives the relative errors of the collision statistics from our parameterization.

### c. Radial distribution function

To quantify the clustering effects in a turbulent cloud system, we also parameterized the RDF. We attempted to simply apply the analytical formula of the collision kernel (2). Therefore, , where and were obtained from the parameterizations derived in the previous section. However, the quotient of and led to amplified error. Even though the average relative error was small (0.19), the parameterized RDF had a group of values below unity, which is unreasonable. Therefore, this simple method was proved incapable to evaluate the clustering effects in a quantitative sense, but it can still help to make a qualitative assessment, which is easier to calculate than the following method we will introduce.

Ayala et al. (2008b) developed an empirical formula for *g*(*R*) that includes both same-sized collisions and cross-sized collisions [see (11) below]. However their parameterization is a function of . This caused a spread of *g*(*R*) under the same EDR but over a broad range of . The average error using this method was 0.18 and was largely caused by the spread. To remove the effect, we modified their formula by replacing the terms that contained with constants in order to fit our data. The overall structure was the same, but the calculation of terms of and was simplified. In their original paper, and depended on the . Here we used constants to replace those terms, and their formulas are given below:

Here, is the Kolmogorov velocity and St = max(St_{1}, St_{2}).

The parameterized matched well the simulations with an average relative error of 0.12 (see Fig. 8c). The red dots from the DNS simulations with values below 1 originated in 5–5-*μ*m collisions and were caused by poor sampling, as these droplets have little possibility to collide with each other, and their pair statistics were based on less than 3000 collisions.

### d. Comparisons with DNS data and other models

In this section, a brief comparison of our prediction model with other studies (Franklin et al. 2007; Saffman and Turner 1956; Wang et al. 1998; Zhou et al. 2001; Ayala et al. 2008b) is given. We compared the turbulence enhancement of the collision kernel between 10–20-*μ*m droplets and that of 10–15-*μ*m droplets with different eddy dissipation rates.

The DNS results of Ayala et al. (2008a) were obtained from their largest (=72.4) case, where we assumed sufficient scale separation has been achieved such that does not affect the collision statistics. Figure 9a shows that the difference-between our DNS results and their results was insignificant. The collision kernel from Franklin et al. (2005) grew faster than the other models by a factor of 1.65 as a result of the difference in calculating flow parameters discussed in section 2b(2), suggesting that their parameterization yields significant overestimation of the collision kernel (see blue dashed line in Fig. 9a). The studies from Wang et al. (1998) and Saffman and Turner (1956) did not consider the clustering effect and consequently underestimate the values (see dashed curves with light green and gray, respectively). The RDF by Zhou et al. (2001) was parameterized based on a nonsedimenting assumption. Thus, one would expect some overestimation (Xue et al. 2008) because droplets would stay in clouds longer without sedimentation. However, even after including the clustering effect given by Zhou et al. (2001), the curve of Wang et al. (1998) was still lower than the DNS data, even though it displayed a similar slope (see dark green curve in Fig. 9a). Pinsky et al. (2006) developed a model in which the turbulence was represented using a statistical formulation. Therefore, large Reynolds numbers were achieved in their study. Since their model did not include clustering effects, underestimation was expected in the collision kernel. We compared their curve of the collision kernel of 10–15-*μ*m-pair statistics with the DNS data (not shown in figure); larger underestimation was found at low intensities (stratiform clouds and cumulus clouds as categorized by their paper in the bottom of Fig. 7). By contrast, our parameterization scheme performed well at low EDRs and deviated slightly from the DNS data as the dissipation rate became higher (see Fig. 9), implying a good performance in predicting stratiform to cumulus clouds (with an EDR less than approximately ). Overall, most of the parameterization schemes underestimated the turbulent enhancement except for the curve of Ayala et al. (2008b) (purple dashed curve). Although our parameterization showed some uncertainties, its results still show a good fit to the DNS.

## 5. Summary and discussion

In this study we examined the major question of warm rain initiation: What role does turbulence play in droplet collisional growth? And can we obtain a more accurate parameterization? With the help of DNS, we quantified the turbulence effect in droplet growth by collision–coalescence based on 39 successful runs (including 5 different eddy dissipation rates and 9 different Reynolds numbers or computational box sizes) with 15 different combinations of droplet pairs in the appropriate range for real clouds to initiate collision–coalescence processes. We studied three collision statistics (RRV, RDF, and turbulent collision kernel) between droplets of radius less than 25 *μ*m using a model configuration different from Franklin et al. (2005), which allowed us to successfully separate the effect of from *ε*, or EDR. We found that all collision statistics were insensitive to over the range of box sizes we can afford to simulate and as long as the model resolved all the smallest scales of turbulence (i.e., the Kolmogorov length scale). This implies that the droplet motions and their collision statistics are mainly affected by the smallest fluid scales and that the large-scale motion has negligible influence. Moreover, the collision statistics were found to be strongly correlated with the EDR and the size of droplets. Therefore, we developed a parameterization that is independent of and and instead is dependent on turbulence intensity (*ε* or EDR), the fluid properties (, , etc.), and the droplet sizes.

We proposed a form of the collision kernel that includes both same-sized collisions and cross-sized collisions and thus is able to describe the continuous growth of droplets by collisions. The parameterized collision kernel was demonstrated to be highly consistent with DNS results. The formula predicted the turbulent collision kernel with an error of 14% on average, better than other proposed formulas. We also parameterized the RRV and RDF, which offers future model developers a measure of the turbulence effect on the relative motion and spatial distribution of droplets. Their average relative errors were 14% and 12%, respectively—again, an improvement over previous studies.

It is recognized that the droplet growth rate is also determined by other factors, such as the collision efficiency and condensational growth. To accurately parameterize the effect of turbulence on cloud droplet growth, the next step is to include these mechanisms.

## Acknowledgments

We thank the three anonymous reviewers for their valuable comments. Computations were mainly made on the supercomputer Mammouth parallèle II from Université de Sherbrooke, managed by Calcul Québec and Compute Canada. The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), NanoQuébec, RMGA, and the Fonds de recherche du Québec—Nature et technologies (FRQNT). Parts of the simulations were run on the supercomputer SGI Altix 4700 from Université de Montréal. This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC).

## REFERENCES

*Computer Simulation of Liquids.*Oxford University Press, 385 pp.

*Proc. AFIPS Spring Joint Computer Conf.*, Atlantic City, NJ, Association for Computing Machinery, 483–485.

*New J. Phys.*,

**10**, 075015, doi:.

*New J. Phys.*,

**10**, 075016, doi:.

*J. Geophys. Res.*,

**112**, D02201, doi:.

*Phys. Fluids*,

**6**, 3742–3749, doi:.

**62**, 2451–2466, doi:.

*Proc. 1997 ASME Fluids Engineering Division Summer Meeting*, Vancouver, BC, Canada, American Society of Mechanical Engineers, FEDSM97-3608.

_{DR}history of Florida cumulus

*J. Geophys. Res.*,

**97**, 2463–2474, doi:.

*Turbulent Flows.*Cambridge University Press, 802 pp.

*Microphysics of Clouds and Precipitation.*2nd ed. Kluwer Academic, 954 pp.

*Phys. Fluids*,

**12**, 2530–2540, doi:.

*A Short Course in Cloud Physics.*3rd ed. Butterworth-Heinemann, 290 pp.

*Parallel Processing and Applied Mathematics*, R. Wyrzykowski et al., Eds., Lecture Notes in Computer Science, Vol. 7204, Springer Berlin Heidelberg, 401–410, doi:.

*New J. Phys.*,

**15**, 045032, doi:.

*Phys. Fluids*,

**6**, 1612–1614, doi:.

*Geofis. Pura Appl.*,

**43**, 243–249, doi:.

*Phys. Fluids*,

**10**, 2647–2651, doi:.

*New J. Phys.*,

**10**, 075013, doi:.