Stochastic coalescence in Lagrangian cloud microphysics

Stochasticity of the collisional growth of cloud droplets is studied using the super-droplet method (SDM) of Shima et al. (2009). Statistics are calculated from ensembles of simulations of collision–coalescence in a single well-mixed cell. The SDM is compared with direct numerical simulations and the master equation. It is argued that SDM simulations in which one computational droplet represents one real droplet are at the same level of precision as the master equation. Such simulations are used to study fluctuations in the autoconversion time, the sol–gel transition and the growth rate of lucky droplets, which is compared with a theoretical prediction. The size of the coalescence cell is found to strongly affect system behavior. In small cells, correlations in droplet sizes and droplet depletion slow down rain formation. In large cells, collisions between raindrops are more frequent and this can also slow down rain formation. The increase in the rate of collision between raindrops may be an artifact caused by assuming an overly large wellmixed volume. The highest ratio of rain water to cloud water is found in cells of intermediate sizes. Next, we use these precise simulations to determine the validity of more approximate methods: the Smoluchowski equation and the SDM with multiplicities greater than 1. In the latter, we determine how many computational droplets are necessary to correctly model the expected number and the standard deviation of the autoconversion time. The maximal size of a volume that is turbulently well mixed with respect to coalescence is estimated at Vmix= 1.5× 10−2 cm3. The Smoluchowski equation is not valid in such small volumes. It is argued that larger volumes can be considered approximately well mixed, but such approximation needs to be supported by a comparison with fine-grid simulations that resolve droplet motion.


Introduction
Coalescence of hydrometeors is commonly modeled using the Smoluchowski equation (Smoluchowski, 1916), often also called the stochastic coalescence equation.It is a meanfield equation that can be derived from the more fundamental stochastic description by neglecting correlations between numbers of droplets of different sizes (Gillespie, 1972;Bayewitz et al., 1974).These correlations are especially important in small volumes and neglecting them can lead to unphysical behavior.For example, when a single drop contains the majority of water in a coalescence cell (gelation), the Smoluchowski equation does not conserve mass for some coalescence kernels (Leyvraz, 2003).
Another limitation of the Smoluchowski equation is that it describes the evolution only of the expected number of droplets of a given size.It does not contain information about fluctuations around this number, which are suspected to be crucial for precipitation onset (Telford, 1955;Scott, 1967;Marcus, 1968).The rate of collisions between droplets depends on their sizes.Small droplets rarely collide with each other, because they are repelled by disturbance flow induced by their settling.Once a droplet reaches a threshold size, it becomes more efficient at collecting smaller droplets.The mean time for a droplet to reach the threshold size is long, but some lucky droplets could reach it much sooner through a series of unlikely collisions.Then they grow quickly, resulting in a quicker onset of precipitation.This effect cannot be described using the Smoluchowski equation.
Moreover, although the Smoluchowski equation can be written for the discrete number of droplets of a given size, it is more often used for droplet concentrations.This adds an additional assumption that the coalescence volume is large, somewhat in agreement with neglecting fluctuations in the Published by Copernicus Publications on behalf of the European Geosciences Union.13510 P. Dziekan and H. Pawlowska: Stochastic coalescence in Lagrangian cloud microphysics number of collisions and correlations in droplet numbers (Gillespie, 1972).
A number of methods alternative to the Smoluchowski equation exist.They are capable of addressing stochastic coalescence, but have some shortcomings that make their use in large-scale cloud simulations impossible.The most accurate one is direct numerical simulation (DNS).In DNS, trajectories of droplets are simulated explicitly and collisions occur when they come into contact.The downside of DNS is that it is computationally extremely demanding.Running a large ensemble of simulations from which statistics could be obtained would take a prohibitively long time.An alternative approach is to use a master equation (Gillespie, 1972).It describes the temporal evolution of the probability of observing a given number of particles of a given size.Collisions are allowed between all particles in some coalescence volume and are assumed to be Markovian; i.e., they only depend on the instantaneous state of the system and not on its history.This can only be justified if the volume is well mixed, i.e., if droplets are randomly redistributed within the volume after each collision.It is worth noting that DNS does not require such assumptions, so it reproduces correlations between positions and sizes of droplets.The master equation was analytically solved only for monodisperse initial conditions with simple coalescence kernels (Bayewitz et al., 1974;Tanaka and Nakazawa, 1993).A more general form of the Bayewitz equation is given in Wang et al. (2006), but cannot be solved for any realistic coalescence kernel.Solving the master equation numerically is extremely difficult due to a huge phase space to be considered.Recently, Alfonso (2015) developed a method to solve the master equation numerically, but was only able to apply the method to a system of up to 40 droplets (Alfonso and Raga, 2017).Alternatively, the stochastic simulation algorithm (SSA) (Gillespie, 1975;Seesselberg et al., 1996) can be used to model a single trajectory obeying the master equation, but obtaining large enough statistics would require very long computations.
Several Lagrangian methods have been developed to model cloud microphysics (Andrejczuk et al., 2008;Shima et al., 2009;Sölch and Kärcher, 2010;Riechelmann et al., 2012).Their common point is that they explicitly model microphysical processes on a small population of computational particles, each representing a large number of real particles.We will refer to these computational particles as super-droplets (SDs).The words "droplets" and "drops" are reserved for real hydrometeors.A thorough comparison of coalescence algorithms from Lagrangian methods was done by Unterstrasser et al. (2017).It led to the conclusion that the method of Shima et al. (2009) "yields the best results and is the only algorithm that can cope with all tested kernels".It was also found to be optimal in DNS tests (Li et al., 2017).In the light of these results, we choose to use the coalescence algorithm of Shima et al. (2009) in this work.Throughout the paper, by the name "super-droplet method" (SDM) we refer to this specific algorithm, and any conclusions regarding the SDM are valid only for the Shima SDM.The Shima et al. (2009) algorithm is not based on the Smoluchowski equation, but, similarly to the master equation, on the assumption that the volume is well mixed.The algorithm introduces some simplifications that may increase the scale of fluctuations in the number of collisions, as described in Sect. 2. These simplifications are not necessary in the limiting case of a single computational particle representing a single real particle, which we call "one-to-one" simulations.Then, the Shima et al. (2009) algorithm should be equivalent to the SSA; i.e., it should produce a single realization in agreement with the master equation.To show that this is true, we compare the Shima et al. (2009) algorithm with the master equation and the SSA in Sect.3. We also compare it with the more fundamental DNS approach in Sec. 4. Once the "oneto-one" approach is shown to be at the same level of precision as the master equation, we use it to study some physical processes that are related to the stochastic nature of coalescence.The way the sol-gel transition time changes with cell size is studied in Sects.3 and 6.We quantify how quickly the luckiest cloud droplets become raindrops.In addition, we use the "one-to-one" approach to validate more approximate methods.The Shima et al. (2009) algorithm with multiplicities greater than 1 is studied in Sect. 4. We determine how many computational particles are required to obtain the correct mean autoconversion time and correct fluctuations in the autoconversion time.Next, in Sect.5, we determine how large a cell has to be for the Smoluchowski equation to correctly represent the rate of rain formation.Throughout the paper we observe that the evolution of the droplet size spectrum strongly depends on the size of the coalescence cell.The size of a well-mixed air volume is estimated in Sect.7 and some implications for cloud simulations are discussed in Sect.8.

The super-droplet method
In this section we present how collision-coalescence is handled in the super-droplet method.Further information about the SDM can be found in Shima et al. (2009).Consider coalescence of water droplets in a well-mixed volume V .Other processes, like water condensation and evaporation, are not included.Thanks to the assumption that the volume is well mixed, all droplets within the volume can collide with each other, independently of their positions (Gillespie, 1972).Therefore droplet motion does not have to be explicitly modeled and droplet coalescence can be calculated in a stochastic manner, as is done in the master equation.Consider two randomly selected droplets i and j .The probability that they will collide during the time step t is P (r i , r j ) = K(r i , r j ) t/V , where r i and r j are their radii and K is the coalescence kernel.We use gravitational coalescence kernels, so the effect of turbulence on coalescence is not studied.
At the heart of the super-droplet method is the idea that many droplets with same properties within a well-mixed volume can be represented by a single computational entity, called the super-droplet (SD).As we are interested only in droplet coalescence within a single cell, it is sufficient if SDs are characterized by two parameters: radius r and multiplicity ξ , which is the number of real droplets that a SD represents.Only integer multiplicities are allowed.In the algorithm of Shima et al. (2009), two simplifications are made that may affect the amplitude of fluctuations in the number of collisions.The first simplification is that SDs collide in an "all-or-nothing" manner.If a collision happens, each real droplet represented by the SD with lower multiplicity collides with a single droplet represented by the SD with higher multiplicity.If the ith and j th SDs collide, their parameters are updated to where primes denote post-collisional values and we assume ξ j ≤ ξ i .Intuitively, one would expect the "all-or-nothing" procedure to lead to larger fluctuations than in a real system, because the number of collision trials is artificially reduced.The second simplification, that we will refer to as "linear sampling", is that instead of considering all N (N − 1)/2 collision pairs, only N/2 non-overlapping pairs are randomly selected.N is the number of SDs in the coalescence cell and x stands for the largest integer equal to or smaller than x.To keep the expected number of collisions equal to the real one, coalescence probabilities are scaled up.The probability of coalescence of two SDs i and j that belong to the same collision pair is P SD (r i , r j , ξ i , ξ j ) = max(ξ i , ξ j ) P (r i , r j ) (N (N − 1)/2)/ N/2 (Shima et al., 2009).Real droplets represented by the same SD cannot collide with each other, because they have the same sedimentation velocities.We will perform two types of simulations.In the "oneto-one" simulations, all SDs have multiplicity ξ = 1.That way the "all-or-nothing" simplification is removed.N 0 superdroplets are initialized by randomly drawing radii from the assumed initial droplet size distribution, where N 0 is the initial number of real droplets in a cell.Coalescence causes one of the SDs to be discarded.Unlike in the original method of Shima et al. (2009), timestep length is adapted at each time step to ensure that none of the collision pairs has a coalescence probability greater than 1.This approach is similar to the direct simulation Monte Carlo method used in diluted gas dynamics (Bird, 1994).In Sect. 3 we show that the "one-toone" method is in agreement with the master equation.
The second type of simulation, in which the number of SDs remains constant (with rare exceptions), is closer to the ideas of Shima et al. (2009).We will refer to it as the "constant SD" simulation.In this type of simulation, the number of SDs is prescribed as N SD and SDs have different multiplicities.Typically, N SD is much smaller than N 0 .We use a novel algorithm for initializing the radii and multiplicities of SDs.The aim of this algorithm is to avoid large differences in the initial droplet size distributions between realizations.Super-droplet radii are not completely randomly drawn from the assumed distribution as in the "one-to-one" simulations.Instead, the assumed distribution is divided into N SD bins and the radius of a single SD is randomly selected within each bin.The bins have equal sizes on a logarithmic scale.Consider an initial droplet size distribution n(ln(r)).We employ a notation in which we omit the division of radius by unit of length whenever the logarithm of a radius is taken; i.e., ln(r) stands for ln(r/µm).The concentration of droplets with radii in the range from r to r + dr is n| r,r+dr = n(ln(r))dln(r).The first step of the initialization is to find the largest and the smallest initial super-droplet radius, r max and r min .They are found iteratively, starting with r min = 10 −9 m and r max = 10 −3 m.We require that they satisfy the condition where r e is either r max or r min and l r = (ln(r max ) − ln(r min ))/N SD .In each iteration, if r min (r max ) does not satisfy Eq. ( 2), it is increased (decreased) by 1 %.Once r min and r max are found, N SD super-droplets are created.The radius of the ith SD is initialized by randomly selecting ln(r i ) in the range (ln(r min ) + (i − 1) l r , ln(r min ) + i l r ].The initial multiplicity of the ith SD is ξ i = n(ln(r i )) l r V + 0.5 .Please note that increasing N SD causes l r to decrease and this in turn gives relatively large values of r min and relatively small values of r max .It means that this initialization procedure does not represent tails of the distribution well, especially for large N SD .It also means that the "constant SD" initialization with N SD = N 0 is not equivalent to the "one-to-one" initialization.Since the large tail is important for coalescence, we add ∞ ln(r max ) n(ln(r))dln(r)V + 0.5 super-droplets with ξ = 1 to the cell.Their radii are selected by randomly drawing ln(r) from the distribution Cn(ln(r))H (ln(r) − ln(r max )), where C is a normalizing constant and H (x) is the Heaviside step function.This makes the actual number of SDs N higher than the prescribed value N SD , typically by ca. 1 %.We do not add SDs from the small tail of the distribution, because very small droplets are of little importance for rain formation.In this type of simulation, the timestep length is constant: t = 1 s.It is not adapted, as it is done in the "one-to-one" simulations, to make the simulation computationally more efficient.Using a constant timestep length can make the coalescence probability exceed unity.If it does, it is assumed that a pair of SDs collides more than once during the time step (Shima et al., 2009).Then, the procedure for calculating post-collisional parameters (Eq. 1) is applied γ = min (γ , ξ i /ξ j ) times, where γ ≥ 1 is the number of collisions between the ith and j th SD and www.atmos-chem-phys.net/17/13509/2017/Atmos.Chem.Phys., 17, 13509-13520, 2017 ξ i ≥ ξ j .Such handling of multiple collisions can cause the expected number of collisions to be lower than the real one if γ > ξ i /ξ j .Another inconsistency is that, rigorously, the probability of collision between SDs should change after each of the γ collisions.For these reasons timestep length should be carefully selected so that multiple collisions do not happen often.If two SDs with identical ξ collide, the multiplicity of one of them drops to zero.Then, the SD with ξ = 0 is used to split the SD with the largest ξ in the system into two.This is slightly different than in the Shima et al. (2009) algorithm, in which a SD with ξ = 0 is used to split the other SD that came out of the collision that caused the multiplicity to drop to zero.Super-droplets are discarded after collision only if all other SDs have ξ = 1.We use an implementation of the SDM from the libcloudph++ library (Arabas et al., 2015).
3 Comparison of the "one-to-one" SDM with the master equation The goal of this section is to show that the "one-to-one" SDM is at the same level of precision as the master equation.To this end, we calculate the average droplet size distribution and the standard deviation of mass of the largest droplet from an ensemble of "one-to-one" simulations.As a reference, we use the results from Alfonso and Raga (2017), who used the master equation approach to study the sol-gel transition.In a system of aggregating particles, the sol-gel transition (gelation) occurs when most of the total mass is located in a single agglomerate (Leyvraz, 2003).For some forms of the coalescence kernel, the Smoluchowski equation is known not to conserve mass after the transition.Alfonso and Raga (2017) present numerical solutions of the master equation for a small cloud volume undergoing the sol-gel transition, for which the Smoluchowski equation is not valid.We perform simulations for the same setup to test whether the "one-to-one" simulations are in agreement with the master equation approach.
Consider a 1 cm 3 volume containing 20 droplets with a radius of 17 µm and 10 droplets of radius 21.4 µm.A gravitational collision kernel is used with collision efficiencies from Hall (1980).Collision efficiencies are bilinearly interpolated in the radius ratio of radii space.Droplet terminal velocities are calculated using the formula from Beard (1976).
Figure 1 shows the average mass distribution obtained using the "one-to-one" simulations with and without linear sampling of collision pairs.The average is calculated from an ensemble of = 10 4 realizations for each case.In simulations without linear sampling, all N (N − 1)/2 collision pairs are considered and a constant time step t = 0.1 s is used.Both approaches give similar results, which shows that the linear sampling technique does not affect the average number of collisions.In addition, the "one-to-one" simulations are compared with the master equation solved using the method of Alfonso and Raga (2017).Both approaches are generally in agreement, with some differences at the large end of the distribution.These differences may be caused by the way the coalescence efficiency tables are interpolated.Another possible source of discrepancies is the numerical diffusion present in the finite-differences method of Alfonso (2015).To test whether the "one-to-one" method also gives correct fluctuations in the number of collisions, the relative standard deviation of mass of the largest droplet σ (m max )/ m max is plotted in Fig. 2.This value is of interest, because the sol-gel transition time coincides with the time at which σ (m max )/ m max reaches maximum (Leyvraz, 2003;Alfonso and Raga, 2017).In Fig. 2, "one-to-one" simulations, with and without linear sampling, are compared with the results of the master equation approach presented in Alfonso and Raga (2017).Please note that Alfonso and Raga (2017) obtained values of σ (m max )/ m max from an ensemble of SSA simulations rather than by solving the master equation, as was the case in Fig. 1.As in Fig. 1, we do not observe any negative effect of using the linear sampling technique, and the "one-to-one" simulations compare relatively well with the SSA.Possible sources of discrepancies are the same as in Fig. 1.Judging from Figs. 1 and 2, we conclude that the "one-to-one" approach is in agreement with the master equation approach.
It accounts for the correlations in the number of droplets per size bin and as such is more fundamental than the Smoluchowski equation approach.
The "one-to-one" SDM with linear sampling is computationally more efficient than solving the master equation directly, or using the SSA.It also puts no constraints on the initial distribution of droplets.Therefore we can use the SDM to predict gelation times for larger systems and more realis- "One-to-one" "One-to-one" no linear sampling SSA tic initial conditions.We use an initial droplet distribution that is exponential in mass n(m) = n 0 m exp(−m/m), where n(m)dm is the number of droplets in the mass range (m, m + dm) in unit volume, n 0 = 142 cm −3 , and m is the mass of a droplet with radius r = 15 µm.This is the same distribution as in Onishi et al. (2015).The total initial number of droplets in the system is N 0 = n 0 V .The results of the "oneto-one" simulations for N 0 up to 10 5 are shown in Fig. 3.For N 0 ≥ 10 2 , the relative standard deviation of the mass of the largest droplet decreases with increasing cell size.This can be understood if we look at a larger cell as an ensemble of 10 smaller cells.Comparing between independent realizations, the variability in the size of the single, largest droplet will be smaller if this droplet is selected from 10 cells in each realization than if it was selected from only a single cell per realization.Interestingly, for N 0 = 10 5 an inflection point appears around t = 500 s.It is not seen in smaller cells.This indicates that some new source of variability is introduced.We believe that it is associated with collisions between large raindrops.We will come back to this in Sect. 5. Intuitively, we would expect the time for most of the mass to accumulate in a single agglomerate to increase with increasing cell size.This turns out to be true for cells with N 0 > 10 3 .For cell sizes 10 2 < N 0 < 10 3 , gelation time is approximately the same, around 300 s.
4 Fluctuations in conversion to raindrops and validity of the "constant SD" SDM Fluctuations in time of conversion of cloud droplets to raindrops were studied using direct numerical simulations by Onishi et al. (2015).Following their notation, by t 10 % we denote the time after which 10 % of the mass of cloud droplets  is turned into droplets with r > 40 µm.Droplets of this size should then quickly grow through coalescence.The time t 10 % is used as a measure of the efficiency of rain production.We will compare the results of the "one-to-one" simulations with DNS and try to determine how many SDs are needed in the "constant SD" simulations to accurately represent coalescence.The initial droplet distribution and coalescence kernel are the same as in Sect.3.
"One-to-one" "Constant SD", "Constant SD" NLS, "Constant SD", "Constant SD", "Constant SD", "Constant SD", Figure 5. Relative standard deviation of t 10 % against cell size.SDM results are based on samples of at least 10 3 realizations.DNS results are taken from Onishi et al. (2015).Where not shown, error bars are smaller than plotted points.The value α = 6 was obtained by curve fitting to the "one-to-one" results.The abbreviation "NLS" in the legend stands for "no linear sampling".
In Fig. 4, values of mean t 10 % for different initial numbers of droplets are plotted against the number of SDs.The results of both the "one-to-one" (rightmost points in each series) and "constant SD" (rest of the points in the series) simulations are presented.For comparison, t 10 % obtained by solving the Smoluchowski equation using the Bott algorithm is plotted (Bott, 1998).In the Bott algorithm, we used t = 1 s and mass bin spacing m i+1 = 2 1/10 m i .The same parameters were used in any Bott simulation presented in this paper.Convergence tests were done for each case.The "one-to-one" results converge with increasing cell volume (i.e., increasing N 0 ) to a value higher than the Smoluchowski result.The difference is probably caused by the numerical diffusion of the Bott algorithm.In the "constant SD" simulations, the error caused by using SDs with ξ > 1 weakly depends on the cell size.Using 10 3 SDs gives t 10 % within 1 % of the "oneto-one" value.Using 10 2 SDs causes about 10 % difference.This shows that, in terms of computational cost, it is relatively cheap to obtain a good estimate of the average result of coalescence using the SDM.The SDM results are also compared with the results of DNS, in which air turbulence was not modeled but hydrodynamic interactions between droplets were accounted for.We choose this kind of DNS because it should be well described by the Hall kernel that is used in the SDM and in the Smoluchowski To analyze the amplification of fluctuations in the "constant SD" method, we plot the relative standard deviation of t 10 % in Fig. 5.For reference, results of DNS from Onishi et al. ( 2015) are shown.Results from our "one-to-one" simu-Figure 6. Points depict the minimal, limiting value of the relative standard deviation of t 10 % for a given number of super-droplets in "constant SD" simulations.For each value of N SD , the minimal value of σ (t 10 % )/ t 10 % is calculated as an average of the respective points to the right of the α/ √ N 0 curve in Fig. 5.The line depicts the fitted function β/ √ N SD with β = 2.
lations are in good agreement with the DNS.Small discrepancies are probably caused by the fact that the DNS included turbulence of various strength for different N 0 .Results of the "one-to-one" simulations were fitted with the function α √ 1/N 0 , resulting in α = 6. Figure 5 also presents fluctuations in the "constant SD" simulations for various N SD .This type of simulations gives correct amplitude of fluctuations only for relatively low values of the ratio N 0 /N SD .For constant N SD , as N 0 increases, the amplitude of fluctuations decreases correctly.Then, above some critical value of the N 0 /N SD ratio, fluctuations stop to decrease and remain constant independent of the cell size.This is a result of introducing unrealistic correlations between droplet sizes, which is a consequence of using low number of computational particles (Bayewitz et al., 1974).To show that the linear sampling technique does not contribute to this effect, we plot result of a "constant SD" simulation without linear sampling for N SD = 32, which is the same as for N SD = 32 with linear sampling.We show the limiting, minimal value of relative standard deviation of t 10 % in Fig. 6.It decreases as β √ 1/N SD , with β = 2.By comparing it with α = 6, we conclude that in order to obtain correct fluctuations in t 10 % using "constant SD" simulations, the number of SDs has to be N SD ≥ 1 9 N 0 .Using so many SDs is not feasible in large eddy simulations (LES), but is possible in smaller-scale simulations.Also, knowing α and β, we can estimate the magnitude of fluctuation amplification in the "constant SD" SDM.

Validity of the Smoluchowski equation
The Smoluchowski equation presents a mean-field description of the evolution of the size spectrum.It is exact only in "One-to-one", 0 = 10 2 "One-to-one", 0 = 10 3 "One-to-one", 0 = 10 4 "One-to-one", 0 = 10 5 Smoluchowski Figure 7. Rain content ratio θ for different cell sizes averaged over ensembles of = 10 3 simulations.Shaded regions show 1 standard deviation interval.
the thermodynamic limit (V → ∞).We will try to determine the minimal cell size for which the Smoluchowski equation can be used without introducing major errors.To do so, we analyze the evolution of θ , the ratio of rain water (r ≥ 40 µm) content to the total water content.Onishi et al. ( 2015) denote this value by τ .We do not adopt this notation to avoid confusion with the characteristic time.
We compare the results of the "one-to-one" simulations with solutions of the Smoluchowski equation for two cases -with fast and with slow rain development.In both cases collision efficiencies for large droplets are taken from Hall (1980), and for small droplets from Davis (1972).In simulations with fast development of rain, we use the same initial distribution as in Sects.3 and 4. As seen in Fig. 7, the Smoluchowski equation gives the correct mean rain development rate for cells with N 0 ≥ 10 4 .The Smoluchowski curve is slightly shifted left, probably due to the numerical diffusion of the Bott algorithm, as discussed in Sect. 4. In cells smaller than N 0 = 10 4 , rain develops slower than predicted by the Smoluchowski equation.Agreement of stochastic coalescence in large cells with the Smoluchowski equation for a similar initial distribution was shown using the SSA by Seesselberg et al. (1996).Onishi et al. (2015) present figures similar to Fig. 7 but obtained from DNS runs for N 0 = 7.24 × 10 4 (Fig. 1b therein).They show good agreement between DNS and the Smoluchowski equation with the kernel of Long (1974), at least up to t = 330 s.If the Hall kernel is used in the Smoluchowski equation, autoconversion is quicker than in the DNS, as discussed in Sect. 4.
Next, we validate the Smoluchowski equation in a setup with slow rain development.This time the initial droplet size distribution is below the size gap, i.e., the range of radii for which both collisional and condensational growths are slow.We use r = 9.3 µm and n 0 = 297 cm −3 as in Wang et al. (2006).In addition, we cut the distribution to 0 for r ≥ 20 µm."One-to-one", 0 = 10 3 "One-to-one", 0 = 10 4 "One-to-one", 0 = 10 5 "One-to-one", 0 = 10 6 "One-to-one", 0 = 10 7 Smoluchowski Figure 8.As in Fig. 7 but for an initial distribution with r = 9.3 µm, n 0 = 297 cm −3 and a cutoff at r = 20 µm.The ensemble size is This cutoff is used in the SDM modeling as well as when solving the Smoluchowski equation.That way we get rid of the occasionally very large SDs present at t = 0 in some realizations of the SDM.For these initial conditions, rain development takes much longer and fluctuations can play a bigger role.Results are presented in Fig. 8. Again, we see convergence of the "one-to-one" simulations to the Smoluchowski result, but in this case the cell has to be larger (N 0 ≥ 10 7 ) for the Smoluchowski equation to be valid.The way the "oneto-one" curves converge to the Smoluchowski curve is interesting.As in the first case, in smaller cells rain appears later than in larger cells.On the other hand, the rain formation rate (the slope of the curves in Fig. 8) in smaller cells starts to decrease at higher values of θ than in larger cells.
In consequence, smaller cells can produce a higher rain ratio than larger ones, although they started producing rain later (e.g., compare the curves for N 0 = 10 5 and N 0 = 10 7 for t > 4200 s).The decrease in the rain formation rate is associated with the decrease in the concentration of raindrops n r , plotted in Fig. 9.The number of raindrops decreases due to collisions between drops from this category.A single drop that is produced in such a collision is less efficient at scavenging cloud droplets than the two pre-collision drops.As a result, the growth rate of θ decreases.Using large well-mixed volumes may introduce additional, unrealistic rain-rain collisions.Consider two droplets within a large cell that independently grow to the rain category.They have to be separated enough so as not to deplete liquid water from each other's surroundings as they grow.If we assume that the cell is well mixed, they can collide immediately after becoming raindrops and generate an even larger drop.In reality, they could collide only after some time after becoming raindrops, because first they would need to overcome the initial separation.This means that using large well-mixed volumes, e.g., in the Smoluchowski equation, may result in an artificial de-Table 1.Average and standard deviation of time (in s) for lucky realizations to produce a single drop with r ≥ 40 µm.γ is fraction of the fastest realizations out of all realizations that were used to compute the respective t 40 γ and σ (t 40 ) γ .The sub-ensemble size γ is shown to give an idea about precision of estimation of the respective statistics.
In coalescence cells with N 0 ≤ 10 4 , we do not observe the decrease in the number of raindrops within 5000 s, probably because the sizes of the raindrops are similar.In larger cells, more raindrops with a broader distribution are formed.In consequence, they collide more often, which decreases their number and the rate of collection of cloud droplets.It is likely that the additional rain-rain collisions in large volumes are responsible for the additional inflection point around t = 500 s in the plot of the relative standard deviation of the largest droplet mass for N 0 = 10 5 (cf.Fig. 3).They could also lead to the deviation from the ∼ 1/ √ N 0 scaling shown in Fig. 5. Fluctuations in cells with N 0 = 10 7 are greater than predicted using this scaling.We also observe that for t ≤ 3000 s, the number of raindrops does not depend on cell size (cf.Fig. 9), but that the mass of rain water does (cf.Fig. 8).In larger cells raindrops acquire larger sizes through collisions with cloud droplets, but the rate of autoconversion of cloud droplets into raindrops is not affected much by cell size.

Lucky droplets
There is a well-established idea that some droplets undergo a series of unlikely collisions and grow much faster than an average droplet (Telford, 1955;Scott, 1967;Marcus, 1968;Robertson, 1974;Mason, 2010).These few lucky droplets are argued to be responsible for droplet spectra broadening and rain forming quicker than predicted by the Smoluchowski equation.Luck is supposed to be especially important during crossing of the size gap, when collisions happen rarely (Robertson, 1974;Kostinski and Shaw, 2005).A single droplet that would cross the size gap through lucky collisions could then initiate a cascade of collisions.Theoretical estimation of the "luck factor" was presented in Kostinski and Shaw (2005).We use the "one-to-one" simulations to test predictions from that paper.
We are interested in the time t 40 it takes until the first droplet grows to r = 40 µm.We perform ensembles of simulations for different cell sizes N 0 .The size of an ensemble is denoted by .The initial distribution is the same as in the second case in Sect. 5.The mean radius is r = 9.3 µm, well below the size gap.The liquid water content is 1 g cm −3 and the concentration is 297 cm −3 , so the smallest cell that has enough water to produce a droplet with r = 40 µm is N 0 ≈ 80. Therefore the smallest cell size we consider is N 0 = 10 2 .For each value of N 0 , we select sub-ensembles of the luckiest realizations, i.e., those with the smallest t 40 .We consider subensembles of size γ with log 10 (γ ) = −4, −3, −2, −1, 0. In each sub-ensemble, we calculate the mean t 40 γ and the standard deviation σ (t 40 ) γ , where the subscript γ denotes the size of the sub-ensemble from which the statistics were calculated.The results for different cell sizes are shown in Table 1.
There is a large variability in t 40 γ with cell size.This is caused by the fact that t 40 depends only on a single largest droplet.Larger cells contain more droplets, so probability of producing single large droplet increases with cell size.We notice that t 40 γ is approximately the same along the diagonals of Table 1.For example, a cell containing 10 6 droplets on average will produce first rain droplet in 30 min.If we divided it into 10 cells with 10 5 droplets each, the luckiest 1 out of 10 would also produce a droplet in 30 m on average.This shows that using large coalescence cells does not affect the formation of the first raindrops.The differences discussed in previous sections emerge later, when there are already some raindrops that can collide with each other.To show that the size of a coalescence cell does not affect the formation of the first raindrops, it is helpful to think of the test presented in this section as a simulation of a large system initially containing N tot = N 0 /γ droplets divided into coalescence cells of size N 0 .We are interested in the mean time t 40 for the first droplet out of all N tot droplets to grow to r = 40 µm.In Fig. 10 we plot t 40 against N tot for different sizes of coalescence cell.It is seen that t 40 does not depend on N 0 .Even using a coalescence cell with N 0 = 10 2 , in which there is barely enough water to produce a drop with r = 40 µm, does not change the results.
Next we calculate the "luck factor", i.e., how much faster the luckiest droplets grow to r = 40 µm compared to the average droplets.To calculate it we use the data for N 0 = 10 2 , because cells of this size can produce only a single droplet with r = 40 µm.Larger cells behave like an ensemble of cells with N 0 = 10 2 as far as t 40 is concerned, so calculating the "luck factor" using t 40 from larger cells would tell us how much faster the luckiest ensemble of droplets produces a rain droplet, compared to an average ensemble -a quantity that we are not interested in.We find t 40 1 / t 40 10 −3 ≈ 5 and t 40 1 / t 40 10 −5 ≈ 11.The value of t 40 10 −5 was estimated at 1366 s based on values along the diagonal for larger γ and larger N 0 .Kostinski and Shaw (2005) estimate that the luckiest 10 −3 fraction of droplets should cross the size gap around 6 times faster than average, and the luckiest 10 −5 around 9 times faster.These values are in good agreement with our observations.

Size of a well-mixed coalescence cell
In the previous sections we have seen that the size of the coalescence cell has a profound impact on the evolution of the system.In this section we estimate the size of a cell that can be assumed to be well mixed.All methods in which the probability of a collision of droplets depends only on the instantaneous state of the cell and not on its history rely on the assumption that the cell is well mixed.This includes the master equation, SSA, the SDM as well as the Smoluchowski equation.The assumption that a cell is well mixed is valid if τ mix τ coal , where τ coal and τ mix are the characteristic times for coalescence and cell homogenization, respectively (Lehmann et al., 2009;Gillespie et al., 2014).By well mixed we mean that droplets should be distributed homogeneously within the cell before every collision.Droplet coalescence generates inhomogeneities, i.e., correlations between droplet positions and sizes.
Rigorously, the characteristic time for coalescence is the mean time between coalescence events, as in diffusionlimited chemical systems (Gillespie et al., 2014).To estimate Figure 10.Mean time until a system of N tot droplets produces the first droplet with r = 40 µm.The system is divided into coalescence cells of size N 0 .The figure is based on the results of "one-to-one" simulations given in Table 1.its magnitude, consider a single large collector droplet falling through a field of smaller droplets.Using a geometric coalescence kernel with efficiency E, the mean time between collisions is τ coal = (Eπ(r l + r s ) 2 v r n s ) −1 , where r l and r s are radii of large and small droplets, v r is the relative velocity and n s is the concentration of small droplets.For r l = 100 µm, r s = 10 µm, v r = 70 cm s −1 , E = 1 and n s = 100 cm −3 , we get τ coal ≈ 0.4 s.
Droplets in the cell can be mixed through turbulence.Turbulence acts similarly to diffusion and its characteristic time for mixing is τ t mix = (V (2/3) /ε) (1/3) , where V is cell volume and ε is turbulent energy dissipation rate (Lehmann et al., 2009).Turbulent energy dissipation rate in clouds is in the range from 10 cm 2 s −3 for stratocumulus clouds to 10 3 cm 2 s −3 for cumulonimbus clouds (Malinowski et al., 2013;Grabowski and Wang, 2013).Let us assume that τ t mix τ coal is satisfied if τ t mix = 0.1τ coal .Even in the most turbulent clouds, this means that the coalescence cell has to be very small V ≈ 1.5 × 10 −2 cm 3 .On average, this volume would contain around one droplet, depending on concentration of droplets.For such small coalescence volumes, the Smoluchowski equation is not valid and the SDM would be very cumbersome, because extremely short time steps would be required.To use larger cells, we need to choose some less strict value of characteristic time of coalescence.Some larger cell sizes, which would be approximately well mixed, could be found phenomenologically through fine-grid simulations including droplet motion.One example of such reference simulations are DNS runs from Onishi et al. (2015) discussed in Sect µm, the Smoluchowski equation gives correct results.This suggests that cells with N 0 ≥ 10 4 can be used in this case.
Another process that can mix droplets is sedimentation.It is difficult to assess its timescale, because it strongly depends P. Dziekan and H. Pawlowska: Stochastic coalescence in Lagrangian cloud microphysics on droplet sizes.Droplets of similar sizes are not mixed by sedimentation, but it is efficient at mixing raindrops with cloud droplets.We can expect that it would prevent depletion of cloud droplets in the surroundings of a rain droplet that was observed for the smallest cells in Sects.3 and 6.Sedimentation acts only in one direction, so it could only allow us to use cells larger only in the vertical direction.

Conclusions
The super-droplet method can exactly represent stochastic coalescence in a well-mixed volume.It was compared with the master equation approach (see Sect. 3) and with direct numerical simulations (see Sect. 4).The precision of the SDM is controlled by the number of super-droplets used.Fluctuations in the autoconversion time are represented well if N SD ≥ N 0 /9.Using smaller N SD increases the standard deviation of the autoconversion time by a factor 1 3
The SDM was used to study stochastic coalescence for two initial droplet size distributions -with small (r = 9.3 µm) and with large (r = 15 µm) droplets.They result in slow and fast rain formation, respectively.Dependence of the system behavior on the size of the well-mixed coalescence cell was observed, especially in the small-droplet case.Cell size not only affects fluctuations in the observables, but also their expected values.If the coalescence cell is small, sizes of droplets are strongly correlated and depletion of cloud water plays an important role.In real clouds, these two effects are probably not manifested, because collector drop sedimentation acts against them.In relatively large cells, raindrops collide with each other more often than in small cells.This leads to a reduction in the rate of conversion of cloud water to rain water, because scavenging of cloud droplets becomes less efficient.In consequence, the highest rain content is produced in cells of intermediate sizes.Possibly, these additional rain-rain collisions can be justified by turbulent droplet motion and sedimentation, but they might also be an artifact caused by using an unrealistically large well-mixed volume.Fine-grid computer modeling with explicit droplet motion could be used to resolve this issue.If the additional collisions were found to be unrealistic, it would mean that cloud models that use large well-mixed cells, e.g., by using the Smoluchowski equation, produce too little rain.
The additional rain-rain collisions do not affect results if droplets are initially large.Then, collisions of cloud drops and raindrops and between cloud droplets are frequent, so the increase in the rate of collisions between raindrops is not important.The mean behavior of the system converges to the Smoluchowski equation results with increasing cell size.Good agreement with it is found for cells with N 0 ≥ 10 4 .
The picture is different if droplets are initially small.Conversion of cloud droplets into raindrops is slow, so the decrease in raindrop concentration due to the additional collisions is relatively more important.The Smoluchowski equation is found to be valid for N 0 ≥ 10 7 for the slow-coalescence case.One could expect that condensational growth will lead to initial conditions with high radii of droplets, for which the additional collisions are not important.Li et al. (2017) have shown that condensation can regulate differences between Eulerian and Lagrangian coalescence schemes.Discrepancies between these schemes that they observed in simulations with condensation and coalescence were smaller than in pure coalescence simulations.
Another aspect of the slow-coalescence scenario is that in it, some lucky droplets can grow much faster than average droplets.We found that a single luckiest droplet out of 1000 grows 5 times faster than average and the luckiest out of a 100 000, 11 times faster.These values are in good agreement with the analytical estimation of Kostinski and Shaw (2005).
We estimate a well-mixed (with respect to coalescence) volume in the most turbulent clouds to be only 1.5 × 10 −2 cm 3 .It is on the order of the volume occupied by a single droplet.Larger cells can be assumed to be only approximately well mixed.For example, in the fast-coalescence case, DNS modeling gives the same results as the Smoluchowski equation (Onishi et al., 2015).Box model simulations using a well-mixed volume with N 0 = 10 4 droplets also give the same results.Therefore it can be assumed that such a volume is approximately well mixed in the case of fast coalescence.In the slow-coalescence case, the well-mixed volume needs to be larger than in the fast-coalescence case for the Smoluchowski equation to be valid.The size of an approximately well-mixed cell for this case can be determined using DNS with initially small droplets.The volume of cells used in LES is typically 10 orders of magnitude larger than a well-mixed volume.The LES cells do not necessarily have to be well mixed.It is sufficient if they are homogeneous, i.e., they are an ensemble of identical, approximately wellmixed sub-cells.Some statistical moments for such ensembles were presented in this work.In general, it is not clear what the size of these sub-cells could be and whether the Smoluchowski equation is valid for them.grant 2012/06/M/ST10/00434. Numerical simulations were carried out at the Cyfronet AGH computer center, accessed through the PLGrid portal.We are grateful to Wojciech W. Grabowski for fruitful discussions.
Edited by: Ilona Riipinen Reviewed by: two anonymous referees

Figure 1 .
Figure1.Mass of droplets per size bin at t = 2500 s.Bins are 1 µm wide.Points depict an averaged result of = 10 4 "one-to-one" simulations with and without linear sampling of collision pairs.Error bars show a 95 % confidence interval.The line depicts a numerical solution of the master equation (see Fig.8in Alfonso and Raga, 2017, data courtesy of L. Alfonso).

Figure 2 .
Figure 2. Relative standard deviation of the mass of the largest droplet in the system.Details of the SDM simulations are given in the caption of Fig. 1.The size of the ensemble of SSA simulations is = 10 3 .The SSA results are taken from Fig. 7 in Alfonso and Raga (2017) (data courtesy of L. Alfonso).

Figure 3 .
Figure 3. Relative standard deviation of the mass of the largest droplet for different cell sizes.Estimated from ensembles of = 10 4 "one-to-one" simulations for each value of N 0 .

Figure 4 .
Figure 4. Mean t 10 % for different cell sizes and different numbers of computational droplets N CD .In SDM simulations, N CD = N SD , and in DNS, N CD = N 0 .The single DNS result is taken from Onishi et al. (2015) (the NoT-HI case therein).Ensemble sizes are ≥ 10 3 for SDM simulations and = 10 2 for DNS.The 95 % confidence intervals are smaller than the plotted points.The rightmost point in each SDM series comes from the "one-to-one" simulations.Other points in SDM series come from the "constant SD" simulations with various values of N SD .The horizontal line is a value obtained by numerically solving the Smoluchowski equation using the flux method fromBott (1998).
equation.It turns out that the Hall kernel gives overly short autoconversion times.The same issue was observed by Onishi et al. (2015) (cf.Fig. 1b therein).

Figure 9 .
Figure 9. Mean concentration of raindrops from the same simulations as in Fig. 8.