The validity of the kinetic collection equation revisited

. The kinetic collection equation (KCE) describes the evolution of the average droplet spectrum due to successive events of collision and coalescence. Fluctuations and non-zero correlations present in the stochastic coalescence process would imply that the size distributions may not be correctly modeled by the KCE. In this study we expand the known analytical studies of the coalescence equation with some numerical tools such as Monte Carlo simulations of the coalescence process. The validity time of the KCE was estimated by calculating the maximum of the ratio of the standard deviation for the largest droplet mass over all the realizations to the averaged value. A good correspondence between the analytical and the numerical approaches was found for all the kernels. The expected values from analytical solutions of the KCE, were compared with true expected values of the stochastic collection equation (SCE) estimated with Gillespie’s Monte Carlo algorithm and analytical solutions of the SCE, after and before the breakdown time. The possible implications for cloud physics are discussed, in particular the possibility of application of these results to kernels modiﬁed by turbulence and electrical processes.


Introduction
The kinetic collection equation (KCE) describes the temporal change of the mean number of particles of mass x i in a given volume of fluid through the process of coalescence and is written as Correspondence to: L. Alfonso (lesterson@yahoo.com) Here n(i, t) can be obtained for t>0 from a given initial spectrum n(i, 0). The coagulation kernel K (i, j ) contains the probability of coalescence of two drops of masses x i , x j . The KCE gives the time rate of change of the average number of i droplets as the difference of two terms, the first term describes the average rate of production of droplets of mass x i due to coalescence between pairs of drops whose masses add x i ,and the second term describes the average rate of depletion of x i droplets due to their coalescences with other droplets. As was pointed out by Gillespie (1975), the KCE is only an approximate time-evolution equation for n (i, t) because the numbers of droplets of different masses are statistically correlated, and the KCE equation contains no definite information concerning the size of the fluctuations from the average, which would be observed in independent realizations of the coalescence stochastic process. Furthermore, for certain collection kernels, the KCE gives nonphysical solutions in which the total mass of the system is not conserved (Drake, 1972;Aldous, 1997). For example, the solution of the KCE using a kernel proportional to the product of the masses of the colliding droplets, features unrealistic behavior such as failure to conserve mass, and divergence of the second moments.
The main goal of our work is to test a numerical criteria for the validity time of the KCE (Inaba et al., 1999), with analytical results obtained for the KCE with kernels for which analytical solutions existed. Because of that, we are not using a realistic collection kernel determined from either laboratory measurements or theoretical flow modeling. The idea Published by Copernicus Publications on behalf of the European Geosciences Union. is to test the numerical results with simple kernels in order to extend the results to real kernels in future works. Drake (1972) carefully analyzed the solutions of the KCE for polynomial kernels of the form A+B(x i +x j )+Cx i x j , and concluded that any polynomial containing an x i x j term is a poor approximation, based on the fact that the non-linear term leads to a time when the second moment of the raindrop distribution becomes infinite, and liquid water content in no longer conserved.
The reason for these behaviors has been previously explained and results from the fact that the KCE are valid only for infinite systems, e.g., systems with large number of particles in large volumes. However, as droplets grow by coalescence, the number of them inevitably decreases, and, as a result, the KCE becomes invalid to describe the process.
This problem is relevant to cloud physics, since the evolution of the large end of the spectrum is crucial in the description of precipitation development. The KCE assumes that the number of particles n (i, t) is a continuous variable. If the collection kernel increases steeply with the mass of the particles, then the collection rate at the high-mass end of the spectrum is significant. A single drop can acquire a mass much larger than the rest of the system and becomes separated from the smooth mass spectrum. In such a situation, the statistical fluctuations at the high-mass end of the spectrum must be taken into account.
A numerical approach to the collection process that takes into account statistical fluctuations is Monte Carlo (MC). Gillespie (1975) first developed an exact Monte Carlo framework for simulating the stochastic coalescence process. Within this framework, all assumptions included in the stochastic collection equation are avoided.
Another way to handle this problem is to study the probability P (n 1 , n 2 , ..., n N ; t) that the system has a drop spec-trumn=(n 1 , n 2 , ..., n N ) at time t. The evolution of the probability distribution P is described by the stochastic coalescence equation (Bayewitz et al., 1974;Lushnikov, 1978;Tanaka and Nakazawa, 1993;Inaba et al., 1999;and more recently Wang et al., 2006). This equation has the form: The first two terms in the right-hand side of Eq.
(2) show the probabilities of transition from other states into the statē n=(n 1 , n 2 , ..., n N ); the last two terms represent those from the staten=(n 1 , n 2 , ...n N ) to other states. The solution of Eq. (2) will produce the complete probabilistic picture of the process and can be used to calculate the true averages as suggested by Bayewitz et al. 1974) and Tanaka and Nakazawa (1994) as: The KCE results from taking the first moments for the particle number distribution (Eq. 3) and assuming that n i n j = n i n j . Under these assumptions Eq. (2) reduces to the KCE. Eq. (2) is very difficult to solve, even numerically since the number of states increases extremely fast with N 0. Analytical solutions were obtained for three cases: sum kernel B(x j +x j ) product C(x i ×x j ) and constant kernel. Bayewitz et al. (1974) obtained an evolution equation for the true mean of the total number of particles for the constant kernel case. Tanaka and Nakazawa (1993) compared the size distributions calculated from Eq. (2) for the three cases, with analytical solutions of the KCE and examined the conditions under which the kinetic collection equation is valid. The stochastic completeness of the KCE was also studied by Valioulis and List (1984).
Going further in this direction we will use the Monte Carlo approach in order to examine the conditions under which the kinetic collection equation is valid. Special attention will be paid to the time evolution and fluctuations of droplet concentration in the large end of the size distribution, which is crucial in precipitation development.
The main result of the present paper will be the test of the numerical criteria suggested by Inaba et al. (1999), to calculate the validity time for the KCE. This result is compared with analytical results previously obtained by Drake (1972) and Tanaka and Nakazawa (1994). We were lead to this conjecture on the basis of numerical simulations with the Monte Carlo algorithm presented in Sect. 2. In Sect. 3, approximating polynomials and analytical solution for the KCE and SCE are presented. Simulations and a comparison with analytical solutions are described in Sect. 4. Finally, in Sect. 5 we discuss the results and possible applications to more general kernels of importance for cloud physics.

The Monte Carlo algorithm
In this study we use the stochastic algorithm developed by Gillespie (1976) for chemical reactions. This algorithm was reformulated to simulate the kinetic behaviour of aggregating systems by Laurenzi and Diamond (1999), by defining species as a type of aggregate with a specific size and composition. In our case, species represent droplets of different sizes.
Within this framework, there is a unique index µ for each pair of droplets i, j that may collide. For a system with N species S 1 ,S 2, ...,S N µ ∈ N (N+1) 2 . The set (µ) defines the total collision space, and is equal to the total number of possible interactions. With this set the collision probability den-Atmos. Chem. Phys., 8, 969-982, 2008 www.atmos-chem-phys.net/8/969/2008/ sity function P (τ, µ) can be determined. This quantity is defined by: P (τ, µ) dτ =Probability that at time t the next collision in volume V will occur in the infinitesimal interval (t+τ, t+τ +dτ ) and will be a µ collision.
Gillespie derives this probability density function for a system of N species as Here µ ∈ N(N+1)

2
. The functions a µ are calculated according to a(i, j ) = V −1 K(i, j )n i n j dt = Probability that two unlike particles i and j with populations (number of particles) n i and n i will collide within the inminent time interval (5) Probability that, two particles of the same species i with population (number of particles) n i collide within the inminent time interval The collision probability density function is the basis of the Monte Carlo algorithm. For calculating the evolution of the system, two random numbers τ and µ must be generated. Equation (4) leads directly to the answers of the aforementioned questions. First, what is the probability distribution for times. Summing P (τ, µ) dτ over all µ (all possible collisions) results in The probability function for reactions can be obtained in a similar way, by integrating the probability density function (pdf) P (τ, µ) dτ over all τ from 0 to ∞ results in Equation (7) shows that the probability of a collision in time follows an exponential distribution. In order to obtain a random pair (τ , µ), according to the probability density function P (τ, µ) we first generate a random number r 1 distributed uniformly in the interval (0,1), then, the inversion method to obtain random numbers is applied. In the inversion method this random number is taken as the probability of a collision in the time period τ according to P 1 (τ ). This probability is obtained by integrating P 1 (τ ) from 0 to τ : Considering that 1−r 1 =r * 1 is also a uniformly distributed random number in the interval (0,1), then the time τ can be calculated from (Eq. 7) in the form: The collision number µ is calculated similarly. A random number r 2 uniformly distributed in the interval (0,1) is generated. Then the pdf P 2 (ν) (Eq. 8) must be integrated over ν until the sum of the µ probability exceeds the random number r 2. The inequality to obtain the collision index µ has the form (Gillespie, 1976) The former results lead to the Gillespie's direct algorithm: 1. Initialize (set initial numbers of species, set t=0, set stopping criteria).
3. Calculate µ according to the distribution P 2 (µ) = a µ α . 4. Change the numbers of species to reflect the execution of a collision.
5. If stopping criteria are not met, go to step 2.

Analytical solutions of the KCE and SCE using polynomial approximations.
The collection kernel for hydrodynamic interactions for the continuous case has the form: In (10), x and y are the masses of the colliding droplets, R(x) is the radius of the larger collector droplet, and r(y) is the radius of the smaller collected droplet, E(x, y) is the collection efficiency and is given by the product of the collision efficiency and coalescence efficiency. For general kernels of the form (12), the KCE has to be solved numerically.
Analytical solutions of the continuous KCE have been obtained by Golovin (1963), Scott (1968), Drake (1972) and  Drake and Wright (1972) for approximations of the hydrodynamic kernel given by the polynomials: Long (1974) calculated the coefficients for the polynomials (13) approximating the collection kernel (10) when the radius of the largest colliding drops is smaller than 50 µm. The results are displayed in Table 1. Other studies (e.g., Scott, 1968) used coefficients up to an order of magnitude larger.
Analytical size distributions of the KCE for the constant, sum and product kernels, are displayed in Table 2, and Table 3 shows the results for total concentration. For the stochastic collection equation (SCE), the true stochastic averages calculated from analytical solutions (Eq. 3) for the sum and product kernel are shown in Table 4 (Tanaka and Nakazawa, 1994).
We have tested the numerical code against the exact size distribution of the SCE reported in Tanaka and Nakazawa (1994) for the sum kernel case (K(x, y)=B(x+y), B=8.83×10 2 cm 3 g −1 s −1 ). The comparison was made for an initial concentration N 0 =100 cm −3 and an excellent agreement was founded (see Fig. 1). Long (1974) demonstrated that for small droplets (R≤50 µm) as terms of higher degree are included in a polynomial, the kernel K(x, y) in the continuous case is approximated adequately. Nevertheless, many authors (Drake, 1972;Pruppacher and Klett, 1997) have claimed that the terms xy give a "nonphysical" behaviour of the solutions, since solution does not conserve mass and there is a divergence of the second moment. The first (M 1 ) and second (M 2 ) moments (with respect to the droplet mass distribution) for the continuous case are defined by

Validity of the KCE for kernels of the form Cxy
where n(x, t) is the droplet mass distribution and x is the droplet mass. The evolution of M 2 with time for kernels containing terms of the form xy diverge as a consequence of the fact that the KCE are valid only for systems with large number of particles in large volumes. This assumption is adequate for most kinetic processes, but the neglect of small population corrections in the KCE causes unrealistic behaviors as the total population of particles become small. Laurenzi and Diamond (2003) studied the case with a xy kernel with a Monte Carlo method and demonstrated that M 2 shows a rapid but finite increase and a rigorous conservation of mass. Drake (1972) calculated the analytical solutions of the KCE for polynomials of the form f (x, y)=Cxy. In this case the second moment evolution is given by Note that when M 2 is undefined. Then for t→τ a single macroparticle remains and M 2 (τ )→∞. The time point t=τ when the deterministic KCE predicts a divergence of M 2 and a decrease of M 1 (first moment, liquid water content) is called the gel point.
We have calculated τ for an initial monodisperse distribution of 100 droplets of 14 µm in radius (droplet mass 1.1494×10 −8 g). The volume of the cloud was set equal to 1 cm 3 . Using the value of C from Table 1 (C=5.49×10 10 cm −3 g −2 s −1 ), then τ in (Eq. 16) is 1378.7 s.
For the same type of kernel (Cxy), by using analytical methods, Tanaka and Nakazawa (1994) (Long, 1974). Here x and y denote the masses of the colliding droplets.
Note: Parameters A, B and C are constants. N 0 is the initial concentration and v 0 is the initial volume of droplets. The index i represents the bin number. Here x i and x j denote the masses of the colliding droplets for bins i and j . KCE described well the coalescence process as long as the mass of the largest droplet was smaller than M 2/3 T , where M T is the total mass of the droplet population.
When the coalescence growth is described by the KCE, the mass spectrum is continuous, as shown by the near-equality of the mass of the largest and second largest droplets. However, as the mass of the largest droplet grows, the number of droplets inevitably decreases and the KCE becomes unable to describe the coalescence process. This larger droplet acquires a mass much larger than the rest of the particles, and becomes detached from the smooth spectrum (Lee, 2000).
We have calculated the ensemble mean of the largest and second largest droplets, M L1 and M L2. The ensemble mean is given by the expression: where N r is the number of realizations of the Monte Carlo process, and M i is the droplet mass in the i−realization of the stochastic algorithm. The average time evolution of the largest and second largest droplets for this particular case is shown in Fig. 2. Inaba et al. (1999), by using a statistical model, found that the stochastic property of the system changes around the stage when the largest droplet mass is in the order of M 2 / 3 T . By using a statistical code for modeling planetary accretion, they calculated the ratio of the standard deviation for the largest particle mass over all the realizations, to the averaged value evaluated from 1000 numerical simulations M L1,S =STD(M L1 )/M L1 . The standard deviation for the largest droplet mass is calculated for each time by using the expression: where M L1 is the ensemble mean of the mass of the largest droplet over all the realizations (given by Eq. 17), N r is the number of realizations of the Monte Carlo algorithm and M i L1 is the largest droplet for each realization. Inaba et al. (1999) found that M L1,S was maximum in the vicinity of M L1 /M Table 3. Analytical solutions of the kinetic collection equation for total concentration calculated with monodisperse initial conditions (Scott, 1968).
Note: Parameters A, B and C are constants. N 0 is the initial concentration and v 0 is the initial volume of droplets. using this magnitude in order to calculate the validity time of the KCE. In order to check when the largest droplet acquires a mass much larger than the rest of the droplets, and becomes detached from the continuous spectrum, we have calculated the behavior of M L1,S evaluated from 1000 realizations (N r =1000) of the Monte Carlo algorithm. The results are displayed in Fig. 3. The maximum of M L1,S was obtained for τ =1315 s, very similar to the analytical estimation from Eq. (16) (τ =1388). The ratio STD(M L1 )/M L1 seems to be a very reliable way for estimating the breakdown time of the KCE. Around this time (τ =1388), the growth of M L2 (second largest droplet) stops while M L1 (largest droplet) continues to grow rapidly (see Fig. 2). After that, M L2 decreases with time because large droplets first coalesce with M L1 and its mass approaches the total mass gradually.
The maximum of the statistic M L1,S was obtained when the largest droplet was about 20 times larger (in volume) than the initial 14 µm droplet. By evaluating this mass in the condition M L1 /M 2/3 T a value of 0.86 was obtained, which is very close to 1 (Fig. 4.), in agreement with the analytical findings of Tanaka and Nakazawa (1994).
As pointed out by Inaba et al. (1999), the time of the maximum depends on the functional form of the collisional kernel. For other type of kernels the maximum will be obtained for different exponents of the total mass of the system M T (in the vicinity of M L1 /M β T =1), here the parameter β has to be estimated.
To illustrate this procedure, the coefficient β in M L1 /M β T ≈1 was calculated for the kernel (where x i and x j are the masses of the colliding droplets). As in the previous cases, we have estimated the validity time for an initial mono-disperse distribution of 100 droplets of 14 µm in radius (droplet mass 1. obtained at 505 s. At this time, the ensemble mean of the largest droplet is 16.91 times larger than the initial 14 µm droplets. Then, the parameter β can be estimated as β=ln(M L1 ) ln(M T )=ln(16.91) ln(100)=0.6141, which is almost equal to 3/5. Then, the KCE is valid until the stage in which the mass of the largest droplet (M L1 ) becomes comparable or larger M 3 / 5 T . For this same kernel, Inaba et al. (1999) estimated analytically that the KCE remains valid before the largest particle (M L1 ) becomes ≈M 3 / 5 T .

Comparison of the solutions of the KCE and SCE for
Cxy kernels. Results for total concentration The degree of accuracy of the solution of the kinetic collection equation is measured by the square relative error, defined by where N are the true stochastic averages calculated from the MC and N the averages from the KCE. As can be observed in Fig. 5, the square relative error shows a sharp increase after M L1 /M 2/3 T =1 (τ ∼1300 s). That means that the expected values calculated according to the KCE will differ from the true averages calculated from the Monte Carlo algorithm. After that time, the KCE breaks down.
For earlier times (t<1300 s) and values M L1 /M 2/3 T ≪1, the Monte Carlo technique produces averages for total concentration that are almost equal to the solution of the KCE.

Comparison of the solutions of the KCE and SCE for
Cxy kernel. Results for the size distribution.
The problem of the size distributions was studied by Bayewitz et al. (1974) for constant kernel solutions of the kinetic and stochastic collection equations. They found that in systems of small population, or in a system partitioned into small compartments, the results of the KCE and SCE may differ substantially, particularly in the long-term tail of the distribution. The same situation was observed by Wang et al. (2006) while comparing the size distributions of the KCE and size distributions from the stochastically complete equation. According to Tanaka and Nakazawa (1994), for product kernel, the solutions of the KCE (n i ) and the SCE (<n i >) agree with each other if the condition x i /M 2/3 T ≪1 is fulfilled (here x i denotes the droplet mass in bin i). The corresponding condition for the sum kernel case is x i /M T ≪1. Here M T is the total mass of the droplet population.
The size distributions obtained from our MC calculations are presented in Figs. 6 and 7, for two different times: t=1000 s and t=1600. This times correspond to values of M L1 /M 2/3 T equal to 0.49 and 1.39 respectively. At earlier times, when M L1 /M 2/3 T =0.49, the KCE size distributions match quite well the SCE size distributions. In contrast, after 1600 s (Fig. 7), the size distributions differ substantially for bin numbers larger than 7.
According to Tanaka and Nakazawa (1994) for sufficiently small mass x i , the solution n i of the stochastic collection equation agrees with that of the KCE equation even in the late stage. Table 4. True stochastic averages calculated from analytical solutions of the stochastic collection equation with monodisperse initial conditions for the sum and product kernels (Tanaka and Nakazawa, 1994).
Note: Parameters Band C are constants. N 0 is the initial concentration and v 0 is the initial volume of droplets. C N 0 k is the binomial coefficient.
Functions f k (T ) can be found by solving successively the equation: ∂f k (t) ij δ i+j,kk C i e −ij t f i (t)f j (t) for monodisperse initial conditions. K(x,y)=Cxy t=1000 sec Analytical Solution KCE Monte Carlo Fig. 6. Size distributions obtained from the KCE and the stochastic approach at t=1000 s for the product kernel.
As Fig. 8 shows, for the product kernel, the KCE and the SCE solutions start to differ for i≥5. When t>1300 s, there is no agreement between the analytical and the Monte Carlo solutions for bin numbers as small as 5, a fact that explained the marked differences observed in Fig. 7. The disagreement increases as we move to the large end of the distribution and time advances.

Validity of the KCE for polynomials of the form B(x+y)
For K(x, y)=B(x+y), the analytical solution of the stochastic collection equation can be calculated easily (see Table 4), and there is no need to use the Monte Carlo integration. As seen in Fig. 10, both analytical solutions for the KCE and the SCE are in excellent agreement even for bin sizes as large as 10. Then, the results displayed in Fig. 10 are in agreement with the less restrictive condition for the sum kernel, that the KCE (n i ) and the SCE ( n i ) solutions agree with each other if the condition x i /M T ≪1 is fulfilled. According to Drake (1972), when K(x, y)=B(x+y), M 2 (t) will exponentially increase with time but still be finite at any time. For Tanaka and Nakazawa (1994), the KCE is valid until a drop with mass comparable with M T appears, i.e., almost until the limit of complete aggregation, when a single macroparticle remains.
We have analyzed this problem by calculating the statistics STD(M L1 )/M L1 . Surprisingly, there is a maximum at τ =1320 s (Fig. 11), indicating that the liquid water content is no longer conserved after 1320 s. We have calculated the evolution of the liquid water content by using the analytical solution of the KCE for monodisperse initial conditions dis-  played in Table 2 according to: At the same time, we have calculated M 1 by using (Eq. 22) and the true averages N(i, t) from the SCE (See Table 3). The results are shown in Fig. 12. After t∼1300 s the total mass calculated with the KCE, starts to decrease, while the total mass calculated with the true averages from the SCE is conserved all the time. This reflects the fact that the stochas-tic approach can predict the behavior of the coalescence process for all times.
The results from Fig. 11 contradict the generalized idea that the KCE with a sum kernel is valid for all times (Drake, 1972). Actually, after some time, the total mass is no longer conserved. After 2000 s we have 82% of the initial mass. The total mass for the product kernel is also plotted indicating that after 2000 s. the remaining mass is only 45% of the initial mass. The smaller reduction in total mass for the sum kernel explains the better agreement between the size distributions ( Fig. 10) for the sum kernel.

Other approximating polynomials containing a xy term
For polynomials of the form A+B(x+y)+Cxy where A=B 2 /C the KCE is valid until the time (Drake, 1972): where L is the initial liquid water content of the droplets (M 1 (0)) . In evaluating τ the values of A, B and C calculated by Long (1974) and displayed in Table 1. The liquid water content (first moment of the distribution) was equal to 1.149×10 −6 g cm −3 (we consider a cloud initially containing 100 droplets of 14 µm in diameter in 1 cm 3 ). The analytically predicted τ for this polynomial form of the kernel was 1134 s while the numerically evaluated value was 1260 s, and the ratio M L1 /M 2/3 T =1.03. In Tanaka and Nakazawa (1994) the condition that the KCE is valid until M L1 /M 2/3 T is smaller than unity was deduced for kernels of the form Cxy. Nevertheless, it seems to work quiet well in general for kernels containing an xy term.
To further study this trend, τ was also estimated numerically for polynomials of the form B(x + y) + Cxy. The coefficients that better approximate the kernel for small droplets are (A=0, B=4.16×10 2 cm 3 g −1 s −1 , C=2.24×10 10 cm 3 g −2 s −1 ). The analytical expression for τ is given by (Drake, 1972;Long, 1974): where L is the initial liquid water content (M 1 (0)) . For L=1.149×10 −6 g cm −3 , from expression (Eq. 24) was ob-tained τ =1508.5 s. From the MC calculations, the maximum of STD(M L1 )/M L1 was 0.504 reached at τ =1310 s and the ratio M L1 /M 2/3 T was found equal to 0.96 (see Fig. 9a and b). The above mentioned results support the fact that the criterium proposed by Tanaka and Nakazawa (1994) for the product kernel in general works well for polynomial kernels containing an xy term.
Simulations with monodisperse initial distributions with the mentioned above kernels, are not very realistic in the cloud physics context, but they are the only way to rigorously check the accuracy of the statistics proposed.

Discussion and conclusions
In this study we have represented the kernels for the continuous case by a series of polynomials and used a MC algorithm to obtain the solutions of the SCE. The solution of the KCE for polynomials containing an xy term predicts an infinite value of M 2 for t=τ . For kernels of the form B(x+y), there is a no conservation of the total mass, but a less pronounce divergence than for the product kernel.
A frequently-encountered process in many fields of science is the random coalescence of small bodies into larger ones, conserving total mass. In astrophysics, the nonconservation of mass after the breakdown of the KCE is usually interpreted to mean that a runaway planet has formed, also known as a gel because of applications in physical chemistry. Astrophysical examples include the coalescence of  planetesimals into planets and of stars into black holes (Lee, 1987), but the largest application area is Physical Chemistry (aerosols, polymerization, phase separation in mixtures).
Since this problem is important in other branches of physical sciences, it is useful to look at different approaches. For condensed matter physicists, the situation that arises for kernels containing an xy term is a phase transition, typically called gelation in the context of coalescence models. Then, when gelation occurs, the mass conservation is expected to break down in finite time i.e.: there exists a T g, called gelation time such that and M 1 (t)<M 1 (0) for t>T g The physical interpretation is that after gelation, some mass is lost under the form of a particle of infinite size, with mass M 1 (0)-M 1 (t), called gel part. The rest of the particles are then called the sol part. For astronomers, the non conservation of the first moment after T g , is usually interpreted to mean that a runaway particle (planet) has formed.
In reality, the KCE describes the continuous mass droplet spectrum (without the gel part or the largest droplet) all the time. When a single droplet is detached from this spectrum then we have a continuous spectrum plus a massive droplet. To further study this trend, we can analyze the time evolution of the liquid water content for a polynomial of the form Cxy for monodisperse initial conditions (see Fig. 12). After t∼1300, the liquid water content starts to decrease. As mentioned above, the neglect of small population corrections causes unrealistic behavior as the total concentration of particles becomes small.
When the gel is formed the largest droplet is detached from the continuous spectrum, the KCE describes only the continuous droplet spectrum (sol part). For example, for t=1700 s, the largest droplet mass is in average 36.52 times larger than the initial 14 µm droplet. Then its mass is equal to 4.197×10 −7 g. At this time, the total mass calculated from the KCE is 7.432×10 −7 g. On the other hand, the initial water content for our simulations was M 1 (0)=1.149410 −6 g cm −3 .
The physical interpretation is that after t=τ , some mass is lost under the form of a particle of big size, with mass M 1 (0)−M 1 (t), the gel part which is not represented by the KCE. The gel mass in this example is M 1 (0)−M 1 (1700)=1.1494×10 −6 g-7.432×10 −7 g =4.1×10 −7 g which is almost equal to the mass of the largest droplet calculated with the MC algorithm (4.197×10 −7 g).
The former analysis confirms the fact that for t>τ the KCE actually models the evolution of the continuous size of the spectrum. As the largest droplet continue to grow by accretion of smaller droplets, the mass of the continuous spectrum will decrease, together with the liquid water content predicted by the KCE. The values of the total water content for the continuous spectrum and the largest droplet (gel part) for several times are shown and the total water content calculated as the sum of the continuous spectrum total water content and the largest droplet mass are shown in Table 5.
Note that the missing mass (M 1 (0)−M 1 (t)) calculated from the KCE is equal (within a 90% accuracy) to the mass of the largest droplet detached from the continuous spectrum for t>τ , and estimated from the Monte Carlo algorithm according to expression (Eq. 17). Then, for t>τ , the mass conservation can be formulated in the form This expression reflects the fact that the "missing mass" actually is transferred to the largest droplet that becomes isolated for times larger than τ . The non conservation of the initial mass when the largest droplet becomes separated from the continuous spectrum n(i, t) predicted by the KCE, explain the differences between the KCE and SCE size distributions after t>τ (Figs. 6, 7 and 8), and the underestimation of the concentration for bin sizes larger than 5. The underestimation in this case is a consequence of non conservation of the liquid water content for the continuous spectrum, when mass is constantly transferred to the largest droplet (gel part). Wang et al. (2006) also observed marked differences between size distributions predicted by the KCE and the SCE (Fig. 7 of Wang et al., 2006). In fact the mass predicted by the KCE is smaller then the mass predicted by the SCE (True Stochastic Collection Equation) for bin numbers smaller than 80.
The numerical criteria STD(M L1 )/M L1 described in this work could be used for calculating τ in the general case, when there is no analytical solutions for the KCE or SCE. One interesting question that arises is the validity of the KCE when turbulence or electrical processes influence the collection process. In these situations, the Monte Carlo algorithm and the already analyzed statistics STD(M L1 )/M L1 for the largest droplet will be useful in the evaluation of the validity of calculations made with the KCE. The alternative, is to use the Monte Carlo algorithm instead of the deterministic tool (Eq. 1).
Another question is the possibility of existence of such  large drops, since the collisional and spontaneous breakup modes will tend to fragment them. In our particular situation the answer is positive, because the collisional and spontanoeus breakup mechanisms will act for larger sizes. For example, at τ =1315 s (calculated for a kernel of the form Cxy), the largest droplet (gel) has a radius of then 38 µm. At t=2000 s the radius of the largest droplet radius is 52 µm. When simulations are performed for larger volumes, the KCE remains valid before a droplet or a number of droplets grow to the mass comparable to the total mass of the system. To further clarify this question, a Monte Carlo simulation was run for a 100 times larger cloud volume (100 cm 3 ), and an initial monodisperse distribution of 10000 droplets of 14 µm in radius (droplet mass 1.1494×10 −8 g). In the simulation, the product kernel was used. The maximum of STD(M L1 )/M L1 was obtained at the same time that in the simulation with the cloud volume of 1 cm 3 , a fact that confirms that the validity time is the same for the two cases. Despite the differences in the initial number of droplets and cloud volumes, the largest droplets (runaway droplets) have a similar size (between 40-50 µm in radius). Then, for larger cloud volumes the KCE is no longer valid after the formation of relatively large droplets that grow faster than the rest.
From the theoretical point of view it will be interesting to check whether the phase transition approach adopted by condensed matter physicists and astronomers could work in a cloud physics context. Long (1974) demonstrated that K(x, y) increases as x 2 for small droplets (R<50 µm) and as x for larger ones. He concluded that for typical continental and maritime clouds, the evolution of the raindrop distribution is closely described if the kernel has the piecewise approximation: 9.44×10 9 x 2 +y 2 or 1.10×x 2 if R≤50 µm (25b) or by 5.78×10 3 (x+y) or 6.33×10 2 x if R>50 µm (26b) In (25) y (26) x and y are the masses of the colliding drops, and R is the radius of the larger droplet. By doing this, he avoided the inclusion of "non-physical" xy terms, that predict to rapid growth for the large drops in a cloud. In other words, by choosing this piecewise approximation, the breakdown of the KCE will be avoided for longer times, since the KCE solutions for kernels of the type B(x+y) are valid until the largest droplet has a mass the limit of complete aggregation (i.e., for all times).
According to (Eq. 12), the time interval for droplets to grow from 20 µm to 100 µm in radius will be in the order of an hour (Pruppacher and Klett, 1997). Nevertheless, for smaller droplets, Cxy approximates quite well the collection kernel. Then, one open question is the possibility of inclusion of xy terms in approximations of K(x, y) when small scale turbulence or other processes are present, in order to predict a faster growth of smaller droplets. From this point of view, precipitation formation could be interpreted as a solgel transition. Several mechanisms have been proposed in the past (entrainment, presence of giant nuclei, supersaturation fluctuations, and effects of air turbulence in concentration fluctuations and collision efficiencies). More recently, a new model of raindrop growth (McGraw and Liu, 2003;McGraw and Liu, 2004) shows how small drops can form and explains some of the differences between continental and maritime clouds. They attacked the problem by extending the theory of statistical crossing of a kinetic barrier to the processes of condensation and collection. But despite all efforts, a conclusive answer is still absent.