Atmospheric Chemistry and Physics Discussions Interactive comment on “ The validity of the kinetic collection equation revisited ”

The paper is dedicated to the analysis of the properties of the kinetic collection equation (KCE) for several types of kernels allowing the analytical solution. The solution of the KCE is compared with a more precise solution of the stochastic collision equation. To estimate the validity of KCE, the stochastic collision equation is solved by the Monte Carlo method. The analysis of the behavior of the equation when particle concentration becomes very low is of special interest. The mass loss in the KCE in case of the formation of a few largest particles collecting all other particles is analyzed. Since the KCE is widely used in different branches of Physics, the results showing some unknown properties of the KCE are of interest, and I would like to recommend the paper for publication.


Introduction
Whether the formation of large droplets trigger the production of rain in warm cumulus clouds remains one of the open problems in cloud physics.Although several mechanisms have been proposed (Pruppacher and Klett, 1997), at present there is no complete explanation for the rapid growth of cloud droplets within the size range of diameters from 10 to 50 µm.Some existing hypotheses try to explain the formation of these large droplets by condensation of water vapor molecules onto droplet embryos (Khain et al., 2000).Many other studies include droplet coalescence as an important factor, mainly through two mechanisms: (i) the collision of large droplets growing on giant and ultra-giant nuclei, and (ii) the self-broadening of the droplet spectra by collisions between cloud droplets.Regarding this second mechanism, it has been emphasized by experimental (Vohl et al., 1999) and theoretical (Pinsky et al., 1999(Pinsky et al., , 2000) ) studies that there is a significant acceleration of droplet growth rate by collisions in a turbulent flow, with collision efficiencies that may reach values 10 times larger than in the pure gravity case.
In this study we focus on a model for the growth of cloud droplets by a fully stochastic turbulent collision-coalescence process and we will show that this model reveals a solgel transition in the system and the formation of runaway droplets.The term sol-gel transition (also known as gelation) is not very familiar in the context of cloud physics, and can be defined as a change from a system with a continuous droplet distribution to one with a continuous distribution plus a massive (runaway) droplet.The kinetic collection or coagulation equation (hereafter KCE) has long been used to model the time evolution of droplet size distributions due to collection events.The discrete version of this equation can be written as (Pruppacher and Klett, 1997) where N (i, t) is the average number of droplets with mass x i , and K(i, j ) is the coagulation kernel related to the probability of coalescence of two drops of masses x i and x j .In Eq. ( 1), the time rate of change of the average number of droplets with mass x i is determined as the difference between two terms: the first term describes the average rate of production of droplets of mass x i due to coalescence between pairs of droplets whose masses add up to mass x i , and the second term describes the average rate of depletion of droplets with mass x i due to their collisions and coalescence with other droplets.Although the term stochastic has been associated with the KCE for historical reasons, it is clearly deterministic and has no stochastic correlations or fluctuations included.
As a matter of fact, the average spectrum obtained from Eq. ( 1) and the ensemble averaged spectrum obtained over different realizations of the stochastic collection process are different.The solution to the KCE and the expected values calculated from the stochastic process are only equal if the covariances are omitted from the probabilistic model, as shown in Bayewitz et al. (1974) and Tanaka and Nakazawa (1993).It is only when this condition is fulfilled that the deterministic solution provided by Eq. (1) corresponds to the average value of n k over many realizations of the stochastic process.
The relevance of fluctuations compared to mean particle growth was discussed many years ago by Telford (1955), Robertson (1974) and Young (1975) and more recently by Kostinski and Shaw (2005).Telford (1955) introduced the probabilistic interpretation assuming that the concentration of droplets available for collection remains unchanged, and that the collecting droplets do not interact among themselves.Robertson (1974) basically follows the same approach as in Telford (1955), using a Monte Carlo procedure to calculate the collection process, also considered an idealized cloud of constant volume and assumed that only drop-droplet interactions are permitted.However, the method for selecting the time between coalescence events is not completely stochastic.Kostinski and Shaw (2005) adopted a version of Telford's approach to illustrate the influence of stochastic fluctuations, which can lead to a factor-of-10 acceleration in the growth of a few lucky drops.However, all these studies follow a quasistochastic approach (see the complete analysis in Gillespie, 1975a), which ignores fluctuations in the droplet mass spectrum.
The full stochastic approach was used by Bayewitz et al. (1974) and Gillespie (1975b), using a constant collection kernel.The stochastic completeness of the KCE was revisited more recently by Wang et al. (2006) for realistic collection kernels, resulting in a novel deduction of a master equation.They demonstrated that, for a system of finite liquid mass and narrow initial size distribution, both stochastic correlations and fluctuations are important.

Approach of this study
We apply the pure stochastic model (Gillespie, 1975a) to perform a fully stochastic collision-coalescence calculation (see Appendix A) following Laurenzi et al. (2002), that inherently incorporates all stochastic correlations present in the collection process.This computational procedure is rigorously based on the probability distribution instead of using the kinetic collection equation.
As the collection process is stochastic in nature, it is, therefore, more accurately described by the master equation for the joint probability distribution P (n 1 , n 2 , ..., n k , ..., t) for the occupation numbers n = (n 1 , n 2 , ..., n k , ...) at time t.This equation can be written as (Bayewitz et al., 1974) K(i, j )(n i + 1)(n j + 1) P (..., n i + 1, ..., n j + 1, ..., n i+j − 1, ...; t) K(i, j )n i n j P ( n; t) (2) Note that the KCE can be recovered from Eq. ( 2) by considering the mean value of n k : and assuming, as in Bayewitz et al. (1974), that n i n j = n i n j .
The first moment of the distribution of N (i, t) corresponds to the total mass in the system (M 1 ); the second moment (M 2 ) for a number of droplet categories or sizes of the discrete distribution (N d ), is defined as Note that M 2 may become undefined if the initial number of droplets is small or if the kernel K(i, j ) increases sufficiently rapidly with x i and x j .In that case, the total mass in the system (M 1 ) starts to decrease.This is usually interpreted to mean that a macroscopic runaway droplet has formed (known as a gel) and the system exhibits a phase solgel transition (also called gelation).After this point in time, the average calculated from the stochastic process will differ from the average obtained from the KCE (Eq.1), and there is a transition from a system with a continuous droplet distribution to one with a continuous distribution plus a massive runaway droplet.Note that when gelation occurs, mass conservation is expected to break down.The gelation time, T gel , is defined as as the largest time such that the discrete model ( 1) has a solution with M 1 (t) ≡ M 1 (0) for t < T gel and M 1 (t) < M 1 (0) for t > T gel , where M 1 is the total mass of the system.
Since analytical expressions for the gelation time only exist for very simple kernels, Inaba et al. (1999) proposed that it could be estimated numerically by Monte Carlo simulations (an approach followed in Alfonso et al., 2008Alfonso et al., , 2010)).
The gelation time, T gel , is the point in time when the maximum of the ratio ρ (see Eq. 5) is reached: where M L1 is the ensemble mean of the mass of the largest droplet and σ (M L1 ) is the standard deviation for the largest droplet mass (σ ) over all the realizations.This standard deviation is calculated as In Eq. ( 6), M i L1 is the largest droplet mass for each realization and K is the number of realizations of the Monte Carlo algorithm.
The gelation time, estimated as described above, can be interpreted as the expected time of formation of the "lucky droplet" that becomes the embryo for runaway raindrops, and in this study, it is estimated from Monte Carlo simulations under turbulent conditions.

Results
It is well known that the KCE (Eq. 1) only has analytical solutions for a few selected kernels such as the product kernel K(i, j ) = Cx i x j .Note that the validity of the KCE will break down once gelation occurs, since the underlying assumption of a continuous distribution is no longer true.We demonstrate that the time when the parameter ρ (Eq.5) reaches its maximum value is a good estimate of the gelation time T gel .We first present our results using the fully stochastic model for a mono-dispersed initial droplet distribution (Sect.2.1), then compare the results from the fully stochastic simulation with the results from the KCE with turbulence for a bidispersed distribution, and (Sect.2.2) and (Sect.2.3) present the results for the non-turbulent case.Finally, in Sect.2.4 we calculate the time evolution of sol concentration, evaluate the performance of the stochastic algorithm and calculate the Monte Carlo errors.

Results for the product collection kernel using a fully stochastic model
The calculations were performed for an initial mono-disperse distribution of 100 droplets of 14 µm in radius (droplet mass 1.15 × 10 −8 g), with C = 5.49 × 10 10 cm 3 g −2 s −1 (Alfonso et al., 2008) in a volume of one cubic centimeter.This initial concentration is typical of maritime cumulus clouds and corresponds to a liquid water content (LWC) of about 1.15 g m −3 .The time evolution of the parameter ρ was estimated from 1000 realizations (K = 1000) of the Gillespie (1976) Monte Carlo algorithm.The results of this simulation are displayed in Fig. 1 and we observe that the maximum of ρ (solid line) was obtained at τ = 1335 s.Independently, the gelation time can be obtained analytically (Drake and Wright, 1972) from and was found to be 1379 s, very close (less than 5 % difference) to the time when ρ reaches its maximum value.After this time, the largest droplet continues to grow by accretion of smaller droplets and the total mass M 1 predicted by the KCE starts to decrease (Wetherill, 1990), as seen in Fig. 1.Thus, the numerical method provides a reliable approximation of the gelation time.

Comparison with KCE results under turbulent conditions
In many models, it is usual to model collisions between droplets under idealized, pure gravity conditions with a collection kernel of the form This hydrodynamic kernel under pure gravity conditions (K g ) does not take into account the turbulence effects, and considers that droplets with different masses (x i and x j , and .Time evolution of the parameter ρ defined in Eq. 5 (dashed line and right axis) and the total mass (solid line and left axis) calculated from the numerical solution of the KCE.The results were obtained for the product kernel K(x,y)=Cxy, (C=5.49×10 10 cm 3 g -2 s -1 ).
Fig. 1.Time evolution of the parameter ρ defined in Eq. (5) (dashed line and right axis) and the total mass (solid line and left axis) calculated from the numerical solution of the KCE.The results were obtained for the product kernel K(x, y) = Cxy, (C = 5.49 × 10 10 cm 3 g −2 s −1 ).
corresponding radii r i and r j ) have different settling velocities.
In turbulent air, the hydrodynamic kernel K g (Eq.9) can be enhanced due to (a) an increase in relative velocity between droplets (transport effect) and (b) an increase in the collision efficiency (the droplet hydrodynamic interaction effect).These effects were considered by implementing the turbulence-induced collision enhancement factor P Turb (x i , x j ) as calculated in Pinsky et al. (2008) for a cumulonimbus, with dissipation rate ε = 0.1 m 2 s −3 and Reynolds number Re λ = 2 × 10 4 , and for cloud droplets with radii ≤ 21 µm.Consequently, the turbulent collection kernel has the form K Turb (x i , x j ) = P Turb (x i , x j )K g (x i , x j ). (10) The simulation of collisions in a turbulent cloud was performed considering a cloud volume of 1 cm 3 and an initial bi-dispersed droplet distribution: 150 droplets of 10 µm in radius, and another 150 droplets of 12.6 µm in radius, corresponding to a LWC of 1.9 g m −3 .Collision efficiencies E(r i , r j ) in Eq. ( 9) are calculated according to Hall (1980).The evolution of the total mass obtained by solving the KCE (Eq. 1) numerically under these conditions is shown in Fig. 2 (solid line).Note that the total mass (expressed in % of the initial total mass) is no longer conserved after 1000 s.The behavior of the parameter ρ (Eq.5) evaluated from 1000 realizations of the Monte Carlo algorithm is also shown in Fig. 2 (dashed line), indicating a maximum at 1055 s, again very close (about 5.5 % difference) to the time when the numerical solution of Eq. (1) breaks down.These results clearly indicate that the sol-gel transition and the formation of a runaway droplet took place around 1000 s and that the parameter ρ (Eq.5) can be used as an estimator of the gelation time when realistic turbulence collection kernels are used.

Comparison with KCE results under non-turbulent, pure gravity conditions
To emphasize the importance of the turbulence enhancement in the collection process, an additional simulation was performed for non-turbulent flow under the Earth's gravitational field with the same initial conditions as in Sect.2.2.The total mass from the numerical integration of the KCE after 2000 s was found to be equal to 99.88 % of the initial mass (see Fig. 3), illustrating mass conservation for this case.Furthermore, the parameter ρ (Eq.5) estimated from the stochastic model never reaches its maximum, confirming that the sol-gel transition does not take place under these conditions.Therefore, the results suggest that the collision-coalescence process under non-turbulent conditions does not show a solgel phase transition for time intervals relevant to the problem under discussion (warm rain initiation), with no generation of a runaway droplet.

Time evolution of sol concentration and performance of the Monte Carlo algorithm
In order to check the performance of the Monte Carlo algorithm for the turbulent collection kernel (Eq.10), the ensemble averages (at each time) for the sol concentration over 1000 realizations were compared with the numerical solution of the KCE.As can be observed in Fig. 4, the time evolution of the sol concentration predicted by the two models was almost identical until the sol-gel transition approaches and the breakdown of the KCE takes place.Although this is a good check of the Monte Carlo algorithm, a formal statistical test (Z-test) test was applied to check whether the solution obtained from the deterministic KCE and the averages over 1000 realizations of the Monte  Carlo method are equal.The null hypothesis would be H 0 : N = N KCE , where N KCE is the sol concentration calculated from the KCE, and N is the true stochastic average calculated using the Monte Carlo method (Eq.11) In Eq. ( 11), N i is the sol concentration for each Monte Carlo realization.As expected, at a 5 % significance level, the null hypothesis H 0 : N = N KCE is rejected only after approximately 900 s due to the breakdown of the deterministic KCE as the sol-gel transition approaches.Following Gillespie (1975b), the errors of the Monte Carlo algorithm were calculated in order to check the accuracy in the calculation of Eq. ( 11).They were estimated as the rms (root mean square) of the fluctuations that may be expected to occur about the average (Eq.11).Then, for the rms fluctuations we have where N is the ensemble average calculated according to Eq. ( 11).The fluctuations about the average are displayed in Fig. 5 as vertical bars together with the time evolution of the sol concentration obtained from the Monte Carlo simulations.If the condition is fulfilled, then the results found in separate realizations will be practically identical (Gillespie, 1975b) and the estimate N will be sufficiently accurate.The ratio (13) was calculated in the time interval [0, 1000] and the condition was always fulfilled (See Fig. 5), corroborating the accuracy of the calculations of N for the turbulent collection kernel.

Discussion and conclusions
One of the outstanding problems in cloud physics is to explain how raindrops can grow by condensation and collisioncoalescence in times as short as 20 min.In order to form a raindrop with a radius of 1mm in a warm cloud, a total of 10 5 droplets with radius of 10 µm must collide and coalesce.When droplets are small and of uniform size, collisions between them are not very efficient and collision events do not occur with the required rate to produce raindrops until some of the droplets grow to a radius of about 20 µm.
The fully stochastic turbulent process presented in this study generates a "lucky droplet" that grows faster than the rest of the droplet population.The appearance of a runaway droplet after gelation is a possible mechanism to explain the formation of raindrops.To further clarify this point, we calculate the time evolution of the mass of the largest droplet and second largest droplets as ensemble averages over all the realizations of the fully stochastic simulations.Figure 6 shows the results for the turbulent case, clearly indicating a significant gap between the mass of the largest and second largest droplets after 1000 s.In contrast, the difference in mass in the non-turbulent, pure gravity case shown in Fig. 7 remains much smaller and with no runaway behavior.
The gelation time (Eq.5) estimated from a set of random realizations constitutes an expected value for the time of the formation of the runaway droplet and has to be compared with the deterministic result obtained from the KCE for a particular system.However, the gelation time is a stochastic variable whose empirical distribution can be calculated from the different realizations of the Monte Carlo algorithm.The problem here is to find a criterion to assess whether a system (for each individual realization) is in runaway growth.For example, Malyshkin and Goldman (2001) declare runaway growth when the coalescence rate of the largest particle is 50 % of the total coalescence rate.Aldous (1997) estimated the runaway time from Monte Carlo simulations via the maximum (over time) of the size of the second largest droplet.A very suitable indicator (Ormel et al., 2010) is M 1 M 2 , i.e. the ratio between the mass of the largest to the second largest  M M (ratio between the largest and the second largest droplet) for one of the realizations of the stochastic process for the turbulent case.droplet in the system.When this ratio increases, the system is in runaway growth, otherwise it is not.We have estimated the time evolution of the ratio M 1 M 2 from each realization to determine the time when this ratio started to increase (Fig. 8 shows the time evolution of the ratio M 1 M 2 for one of the realizations).Thus, an empirical distribution for the gelation times was obtained by generating 1000 realizations of the Monte Carlo algorithm (Fig. 9).A similar approach was taken to generate the distribution for the runaway particles radii at the gelation time (Fig. 10).However, in our case the value of N varies, as the number of collisions required to form a runaway droplet is different for each realization of the Monte Carlo process.
The fully stochastic simulations under turbulent conditions performed here include a collision enhancement factor for collisions between droplets with radii ≤ 21 µm, so the role of turbulence in producing the runaway droplet is likely underestimated in the present study.Since the nucleation and condensation processes are not yet included in this model, future developments will attempt to include the combined effect of turbulent collection and condensation (McGraw and Liu, 2003) on droplet growth.
Fig.1.Time evolution of the parameter ρ defined in Eq. 5 (dashed line and right axis) and the total

Fig. 2 .
Fig. 2. Time evolution of total mass calculated from the numerical solution of the kinetic collection equation (KCE) for turbulent collision coalescence (solid line and left axis) and the parameter ρ (dashed line and right axis) estimated from the Monte Carlo algorithm.

Fig. 2 .
Fig. 2. Time evolution of total mass calculated from the numerical solution of the kinetic collection equation (KCE) for turbulent collision-coalescence (solid line and left axis) and the parameter ρ (dashed line and right axis) estimated from the Monte Carlo algorithm.

Fig. 3 .
Fig. 3.Same as Fig. 2 but for the pure gravity case.

Fig. 3 .
Fig. 3. Same as Fig. 2 but for the pure gravity case.

Fig. 4 .
Fig. 4. Time evolution of the sol concentration obtained from the numerical solution of the kinetic collection equation (solid line) and from the Monte Carlo algorithm (crosses).

Fig. 4 .Fig. 5 .
Fig. 4. Time evolution of the sol concentration obtained from the numerical solution of the kinetic collection equation (solid line) and from the Monte Carlo algorithm (crosses).

Fig. 5 .
Fig. 5. Time evolution of the sol concentration obtained from Monte Carlo algorithm (crosses and left axis) fluctuations about the average (vertical bars and left axis) and the ratio σ (N(t)) N(t) (solid line and right axis).

Fig. 6 .
Fig.6.Time evolution of the ensemble means over all the realizations for the largest (solid line) and second largest (dashed line) droplet masses (expressed in multiples of a 10µm droplet mass) for the turbulent case.

Fig. 6 .
Fig.6.Time evolution of the ensemble means over all the realizations for the largest (solid line) and second largest (dashed line) droplet masses (expressed in multiples of a 10 µm droplet mass) for the turbulent case.

Fig. 7 .
Fig. 7. Same as Fig. 4 but for the pure gravity case.

Fig. 8 .
Fig.8.Time evolution of M 1 M 2 (ratio between the largest and the second largest droplet) for one of the realizations of the stochastic process for the turbulent case.

Fig. 9 .
Fig. 9. Empirical distribution for the gelation times obtained for the turbulent case.The results were generated from 1000 realizations of the Monte Carlo algorithm.

Fig. 9 .Fig. 10 .
Fig. 9. Empirical distribution for the gelation times obtained for the turbulent case.The results were generated from 1000 realizations of the Monte Carlo algorithm.

Fig. 10 .
Fig. 10.Empirical distribution for the radius of the runaway droplets.The results were generated from 1000 realizations of the Monte Carlo algorithm.