The impact of fluctuations and correlations in droplet growth by collision-coalescence revisited. Part I: Numerical calculation of post-gel droplet size distribution

Abstract. The impact of stochastic fluctuations in cloud droplet growth is a matter of broad interest, since stochastic effects are one of the possible explanations of how cloud droplets cross the sizegap and form the raindrop embryos that trigger warm rain development in cumulus clouds. Most theoretical studies in this topic rely on the use of the kinetic collection equation, or the Gillespie stochastic simulation algorithm. However, the kinetic collection equation is a deterministic equation 15


1. Introduction
Although rain has been observed to form in warm cumulus clouds within about twenty minutes, 35 calculations that represent condensation and coalescence accurately in such clouds have had difficulty producing rainfall in such a short time except via processes involving giant cloud condensation nuclei (with diameters larger than 2 μm). One of the possible origins of this discrepancy is the stochastic nature of the collision coalescence process that is not well reflected in current models that rely almost exclusively in the kinetic collection (or Smoluchowski) equation, 40 hereafter referred to as KCE (Pruppacher and Klett, 1997): where N(i,t) is the concentration of droplets in bin i and K(i,j) is the collection kernel for droplets in bins i and j. Additionally, Eq. (1) fails to represent the droplet size distribution at the time when raindrop embryos are formed (Alfonso et al., 2008), as there is a transition from a continuous 45 distribution to a continuous distribution plus a massive raindrop embryo (i.e, the "runaway droplet"). At that point, the infinite system exhibits a sol-gel transition (also called gelation and in which the "runaway droplet" is labeled "gel"), the KCE breaks down and the total mass of the system calculated according to the KCE is no longer conserved.
One way to avoid this problem is to adopt the stochastic finite volume description of the coalescence 50 process by using the stochastic simulation algorithm proposed by Gillespie (1975). The stochastic simulation algorithm (hereafter referred as SSA), correctly accounts for fluctuations and correlations, and has been used in cloud simulation studies with realistic collection kernels 4 (Valioulis and List, 1984). However, the SSA has difficulties in accurately reproducing the large end of the droplet size distribution. This is due to the huge number of realizations required to obtain 55 a smooth behaviour at the large end of the droplet size distribution (Alfonso, 2015

  
(2) 60 The master equation (2) is a gain-loss equation for the probability of each state () Pn . The sum of the first two terms is the gain due to transition from other states, and the sum of the last two terms is the loss due to transitions into other states. This formulation was introduced in the pioneer works of Marcus (1968) and Bayewitz et al. (1974), and was studied in detailed by Lushnikov (1978, 65 2004) and Tanaka and Nakazawa (1993). However, these studies only offer analytical results for a limited number of cases (with constant, sum and product kernels), for mono-disperse initial conditions. Furthermore, most of these studies are limited to non-gelling conditions and do not provide a coherent framework for the general case.
The exception are the methods developed by Lushnikov (2004) from the analytical solution of the 70 master equation, and more recently by Matsoukas (2015), the later based on arguments from statistical physics. These methods, although also limited to very special cases (product kernel and mono-disperse initial conditions), are capable of obtaining solutions in the post-gel regime. For 5 example, in Lushnikov (1978Lushnikov ( , 2004, the coalescence process takes place in a system with a finite volume that includes a finite number of particles. Within this approach any losses of mass are, by 75 definition, excluded. In the infinite system described by the KCE (Eq. 1), the coagulation process instantly transfers mass to the gel, while in the finite system the gel coalesces with smaller particles decreasing their concentration not instantly by rather in a finite time.
In order to study the droplet size distribution after the formation of raindrop embryos (sol-gel transition), for systems with kernels relevant to cloud physics and arbitrary initial conditions, we 80 must rely on numerical methods capable of solving the master equation (Eq. 2). We can address this problem through a detailed comparison of the droplet size distributions obtained from the stochastic description for a finite system with the master equation (Eq. 2), and the deterministic approach for an infinite system by using the KCE (Eq. 1), using the numerical algorithm reported in Alfonso (2015). By the time the gel forms, certain differences are to be expected between the 85 two approaches at the large end of the droplet size distribution.
This analysis of the sol-gel transition problem in the cloud physics context could provide an alternative explanation of the differences between modeled and observed droplet spectra in clouds.
Several mechanisms have been proposed in the past (entrainment, presence of giant nuclei, supersaturation fluctuations, effects of air turbulence in concentration fluctuations and collision 90 efficiencies, effects of film forming compounds on droplet growth), and a large amount of literature exists regarding the variety of mechanisms that may explain this disparity, but a conclusive answer is still absent. This study does not attempt to dispute any of the mechanisms already proposed, but to explore another mechanism that has not yet been widely considered in the mainstream literature. 6 The paper is organized as follow: Section 2 presents an overview of the numerical algorithm 95 (following Alfonso, 2015). Numerical results (for the product and hydrodynamic kernels, respectively) with a detailed analysis of the method for calculating the sol-gel transition time, and a comparison with averages calculated with the KCE are presented in Sections 3 and 4. Finally, Section 5 presents a discussion of the limits of applicability of the KCE, an example of correlations in the critical region and conclusions. For a system consisting of N monomers at t=0, the number of possible configurations increases exponentially and can be approximated from the equation (Hall, 1967): For example, R(50)=217590 and R(100)=190 569 232.
The procedure is illustrated for a system with 5 monomers in the initial state, only for the purpose of demonstrating the method. As the system in this case has only 6 possible states, it is much easier to explain the details of the calculations. The 6 possible configurations generated from the initial 120 state (5, 0, 0, 0, 0) are displayed in Fig. 1.
In a second step, the probabilities of all generated configurations are updated according to the first order finite difference scheme (Alfonso, 2015): ; P n t t  at t=t0+Δt will increase if the states from which transitions are allowed, have a non-zero probability at t = t0 (second and third 8 terms in the right-hand size of Eq. (5)), and will decrease due to collisions of particles from the same state at t = t0 (fourth term and fifth terms in the right-hand side of Eq. (5)) if 0 ( ; ) P n t is positive.
The finite difference equation for (1,0,0,1,0) P is presented to illustrate the method. From the 130 generation scheme displayed in Figure 1, note that the only allowed transitions to (1, 0, 0,1, 0) are from the states (1, 2, 0, 0, 0) and (2,0,1,0,0) . Consequently, at t=t0+Δt, The time evolution of the probability of each state was calculated for the product kernel 140 ( , ) ij K i j Cx x  , considering C= 5.49x10 10 cm 3 g -2 s -1 following Long (1974) and for the initial condition (5, 0, 0, 0, 0;0) 1 P  . Due to the small number of droplets in the initial configuration (only 5), the simulated volume was set equal to 10 -2 cm 3 , with an initial droplet radius of 17µm. The time ( After the calculation for each state is completed, the expected values for each droplet mass can be found from the relation (Alfonso, 2015): where the discrete probability mass function is calculated from the state probabilities following the expression:

Estimating the time of gel formation
Lushnikov (2004) demonstrated that right after the sol-gel transition, the particle mass distribution splits into two parts: the thermodynamically-populated one whose behavior is described by the kinetic collection equation, and a narrow peak with a mass very close to the gel mass. For the infinite system described by the kinetic collection equation (Eq. 1) with kernel K(i,j)=Cxixj, the critical time 160 is calculated when the second moment of the distribution diverges: leading to the critical time of the sol-gel transition: After τ =Tgel the second moment becomes undefined, and the total mass of the system starts to 165 decrease.
For a finite system, the standard deviation ( ) of the mass of the largest droplet is an important quantity in order to calculate the critical time of the gel formation (Botet, 2011). At the critical time for the infinite system  must diverge, since it is proportional to the second moment of the This was explored in previous studies (Inaba, 1999;Alfonso et al., 2008Alfonso et al., , 2010Alfonso et al., , 2013, where  was calculated for a finite system from Monte Carlo simulations in order to estimate the sol-gel 175 transition times for the corresponding deterministic model of an infinite system. We can perform an example calculation of  by using the species formulation of the SSA (Laurenzi and Diamond, 2002), in this case: for each realization at a given time, and Nr is the  Figure 3 for a finite 11 system with N=40 droplets of 17 μm in radius (droplet mass = 2.058×10 -8 g) in a volume of 1 cm 3 .
For the product kernel with C= 5.49x10 10 cm 3 g -2 s -1 , the maximum occurs at T=1065 s, which is close to the sol-gel transition time Tgel=1075 s. The evolution of a system with an initial mono-disperse droplet size distribution of N0=40 droplets of 17 μm in radius at t0 , in a volume of 1cm 3 , and a corresponding liquid water content (LWC) of 190 0.823 gm -3 , was calculated using the master equation (Eq. 2). The initial condition for this case is (40, 0,..., 0;0) 1 P  and the time step was set equal to Δt=0.1s. The complete phase space of a population with N=40 monomers contains 37338 states, and the master equation must be solved for the 37338 states. Then, for each state there is a finite difference equation (5), or an equation similar to Eq. (6), but with 40-dimensional state vectors. 195 A discrete 40 bin grid was defined for our model. The mass for bin 1 is taken to be the mass of a 17 μm in radius droplet, and the mass of bin n is n times the mass of bin 1. Then, if all these 40 droplets in the initial distribution were to coalesce into a single droplet, the final droplet radius would be 58.14 μm in diameter and would belong to the mass bin 40.
The results for the droplet mass distribution are displayed in Figure 4 at t=300, 1000, 1800 and 200 2200s. Note that the gel is clearly seen in the distributions at 1800 and 2200s but not at 1000s. 12 To proceed further, the previous results are compared against the analytical size distributions from the KCE (Eq.1) calculated for the product kernel with mono-disperse initial conditions before (t=300s) and after (t=1200s) gel formation (Laurenzi and Diamond, 2002): 205 In Eq. (12), N0=40 is the initial concentration and v0 is the initial volume of droplets. The index i represents the bin size and C= 5.49x10 10 cm 3 g -2 s -1 .
The comparison of the droplet mass concentration (  Figure 5 at t=300, 1200s. Note that the KCE fails to capture the 210 gel formation after the critical time, and the droplet mass concentration calculated using the kinetic approach is much lower at the large end of the distribution. This is due to the fact that the total mass calculated according to Eq. 12 decreases after the sol-gel transition time. This decrease can be clearly seen in the time evolution of the LWC from the kinetic approach using the relation: where () mi is the mass for bin size i. After t~1000s ( Figure 6) the total mass of the system calculated according to the KCE starts to decrease, while the total mass calculated from the stochastic approach is conserved at all times.

13
Within the master equation approach, the expected value of the mass of the largest droplet M1 is approximately given by (Tanaka and Nakazawa, 1994): where i n are is the expected number for each droplet size values calculated from Eq. (7), () mi is the mass for bin size i, and the bin number i1 is defined from the relation: The mass of the gel, 1 M , is evaluated for t=1200, 1800 and 2200s from Eqs (14) and (15).
where Nr =1000 for this simulation and M i 1 is the largest droplet for each realization. The results obtained from both the master equation and the SSA are displayed in Table 1, showing a very good correspondence between the two approaches.
After the sol-gel transition the mass of the gel can also be estimated by using the infinite system 235 approach from the relation (Wetherill, 1990):

Results for the hydrodynamic collection kernel
Collisions between droplets under pure gravity conditions are simulated with the hydrodynamic kernel, which has the following expression: 245 The hydrodynamic kernel takes into account the fact that droplets with different masses (xi and xj and corresponding radii, ri and rj) have different terminal velocities, which are functions of their masses. In Eq. 18, E(ri ,rj) are the collection efficiencies calculated according to Hall (1980).  Figure 7 shows that there is a maximum at t=1310s, which is considered a good estimate for the sol-gel 260 transition time for the infinite system. Thus, the distributions obtained from the stochastic (master equation) and the deterministic (kinetic collection equation) approaches must be compared before and after 1310s. After the critical time, the gel mass was calculated using Eqs. (14) and (15). The results are displayed in Table 2, showing, again, a good agreement between the calculations performed with the SSA and the master equation. 265 Although for the hydrodynamic kernel the critical time was larger than 20 minutes, we must emphasize that, in general, for concentrations larger than 30-40 cm -3 , smaller critical times must be obtained. For example, for kernels proportional to the product of the masses, Malyshkin and Goldman (2001)  to cloud physics, we have a similar situation (a decrease in the time of occurrence of the phase transition as the number of particles in the initial distribution increases). A more detailed discussion of this problem for realistic collection kernels can be found in Alfonso et al. (2010Alfonso et al. ( , 2013. 275 with the deterministic (kinetic) approach

Calculation of the post-gel droplet size distribution from the master equation and comparison
The evolution of a system with the initial bi-disperse droplet size distribution described in the previous section is calculated here using the master equation (Eq. 2) with the initial condition 16 (20,10,..., 0;0) 1 P  , and a time step Δt=0.1s. The results for the expected droplet mass distribution as a function of radius are displayed in Figure 8 at four different times (t = 500, 1500, 1800 and 280 2500s).
Before the sol-gel transition, the mass spectrum exponentially decreases with increasing drop radius for both the KCE and the master equation. After the sol-gel transition, there are two behaviors in the droplet mass distribution of the master equation: i) an exponential decay that resembles the KCE description, and ii) a peak in the gel fraction of the distribution, in which the mass is calculated 285 according to Eqs. (14) and (15). As can be observed in Fig. 9, there are substantial differences between the kinetic and the stochastic approaches especially in the large end of the distribution after the critical time, with much higher values of the droplet mass concentration for the stochastic case.

290
In their pioneering studies using the stochastic framework, Marcus (1968) and Bayewitz et al. (1974) solved the stochastic master equation (Eq. 2) for a constant collection kernel and a monodisperse initial droplet distribution. The later study revealed significant deviations from the KCE when there is a small number of droplets in the initial distribution ( 50 Toral N  ). In our work, we extended the results obtained by these authors by calculating the expected droplet size distributions 295 for small systems within the stochastic framework, but by using mass dependent and relevant to cloud physics collection kernels (e.g multiplicative and hydrodynamic). Our results confirm the findings of Bayewitz et al. (1974) that in systems of small populations the results of the kinetic deterministic equations approach may differ substantially from the stochastic means at the large end of the droplet size distribution. 300 The application of the KCE to coagulating systems also requires that the particles be well-mixed (Bayewitz et al., 1974;Sampson and Ramkrishna, 1985), implying that every pair of droplets is always available for coagulation with each other (Sampson and Ramkrishna, 1985).
Another important assumption is that the droplet population is sufficiently large so that the existence of a droplet with particular properties is not conditionally dependent on the existence or non-305 existence of any other droplet. In other words, no correlations are assumed in the system, so that i j i j n n n n  .
The assumption that the system is sufficiently large is linked to the fact that the KCE is a Additionally, the KCE can fail even if the number of droplets is large when a raindrop embryo forms. At that critical time, there is a transition from a continuous droplet distribution to a continuous distribution plus a raindrop embryo (or runaway droplet). This sol-gel transition is well 18 known in other fields (e.g astronomy), but has not been sufficiently explored in the context of cloud 320 physics, where the gel would correspond to the raindrop embryo. This approach is developed in this paper through a detailed comparison of expected values calculated from the stochastic framework with averages obtained from the KCE for realistic collection kernels, before and after the sol-gel transition time.
The marked differences between these two approaches at the sol-gel transition can be related with 325 the increase of correlations at the critical point, and that can happen even for a large number of particles in the initial distribution (Malyshkin and Goldman, 2001). In principle, this analysis could be performed by using the SSA, which is an alternative tool for the master equation formalism (Eq. 2). However, the number of realizations required to obtain a smooth behavior at the large end in order to compare with averages from the KCE, would be extremely large.
It is necessary to emphasize that our method (although it can be computationally expensive) works 355 for any type of kernels, whereas the analytical techniques developed by Lushnikov (2004), and Matsoukas (2015) work only for very special cases.
The failure of the KCE to capture the gel-formation could provide an explanation of the inability of explicit microphysics cloud models to explain the droplet spectral broadening observed in small, 20 warm clouds. Therefore, even for large simulation cells, the use of the KCE is justified only in the 360 absence of the sol-gel transition.
For the small volume approach described in this paper, a model that considers the interaction between small coalescence volumes through sedimentation or other physical mechanisms for realistic collection kernels is needed. For a constant collections kernel, this theory was outlined by Stepanov (1990, 1991) based on a scheme proposed by Nicolis and Prigogine 365 (1977) for chemical reactions. Within this theory, the whole system is subdivided into spatially homogeneous sub-volumes (coalescence cells) that interact through the diffusion process, and the coalescence events are permitted only between droplets from the same sub-volume. As a result, we obtain a set of master equations for each sub-volume. Although very complex, it could be a starting point in order to consider the interactions between small coalescence volumes through different 370 physical mechanisms.