2.1 Results for infinite systems in coagulation and percolation theory

The most commonly accepted approach to modeling the collision coalescence process in cloud models with detailed microphysics relies upon the Smoluchowski kinetic equation or kinetic collection equation (KCE), governing the time evolution of the average number of particles. The discrete form of this equation can be written as follows (Pruppacher and Klett, 1997):

where N(i,t) is the average concentration of droplets with mass x_i at time t, and K(i,j) is the coagulation kernel related to the probability of coalescence of two drops of masses x_i and x_j. In Eq. (1), the first term on the right-hand side describes the average rate of production of droplets of mass x_i due to coalescence between pairs of drops, whose masses add up to mass x_i, and the second term describes the average rate of depletion of droplets with mass x_i due to their collision and coalescence with other droplets.

However, the KCE may have a serious limitation in some cases (Lushnikov, 2004) and hence cannot accurately describe the coagulation process. The limitation essentially lies in the fact that the coagulation equation inevitably creates particles with infinite mass. For example, for a multiplicative coagulation kernel ($K(i,j)=C{x}_i{x}_j$), an attempt to calculate the second moment of the droplet mass spectrum:

leads to the result

Thus, after t=T_gel, the second moment may become undefined, and the total mass of the system starts to decrease (see Appendix A for more details). This result applies to infinite (with negligible fluctuations and correlations) coagulating systems in the thermodynamic limit, which is the limit for a large number K of particles where the volume V is taken to grow in proportion with the number of particles. Then, in the limit $K,V\to \mathrm{\infty },K/V\to N<\mathrm{\infty }$. The infinite system interpretation of the sol–gel transition assumes the presence of a gel phase (which is not predicted by the KCE equation) and introduces an additional assumption as to whether or not the gel interacts with the infinite size clusters that are not described by the KCE equation.

The other scenario considers that coagulation takes place in a system with a finite number of monomers in a finite volume. This approach is based on the scheme developed by Markus (1968) and Bayewitz et al. (1974) and was studied by Lushnikov (1978, 2004), Tanaka and Nakazawa (1993, 1994), and Matsoukas (2015) by using analytical tools and more recently by Alfonso (2015) and Alfonso and Raga (2017) numerically. Within this approach there is no mass loss, and the phase transition is manifested in the emergence of a giant particle that contains a finite fraction of the total mass of the system. Solutions in the post-gel regime were obtained analytically by Lushnikov (2004) and Matsoukas (2015) and numerically by Alfonso and Raga (2017).

The sol–gel transition has been observed experimentally. For example, aerogels in aerosol physics (Lushnikov et al., 1990) and in other theoretical models, such as that of percolation (Botet and Płoszajczak, 2005; Kolb and Axelos, 1990), where there is a close analogy between percolation and droplet coagulation. In bond percolation, each lattice corresponds to a monomer, and a proportion p of active bonds is set randomly between sites. Then clusters of size s are defined as an ensemble of s-occupied sites connected by active bonds. For a definite value of p=p_c, a macroscopic cluster appears, corresponding to the sol–gel transition.

Recent results in percolation theory show that the largest cluster follows the Gumbel distribution for subcritical percolation (Bazant, 2000) and, at the critical point, follows the Kolmogorov–Smirnov (K-S) distribution (Botet and Płoszajczak, 2005). The K-S distribution is the distribution of the maximum value of the deviation between the experimental realization of a random process and its theoretical cumulative distribution, and it has following the cumulative distribution:

or the equivalent expression:

Botet and Płoszajczak (2005) also found evidence (from numerical solutions of the KCE equation) that, for multiplicative coalescence (with a collection kernel proportional to the product of the masses), the largest cluster follows the distribution in Eqs. (5a) and (5b) at the time of the phase transition. At this point, a hypothesis is formed in which the results obtained in percolation are extrapolated in order to find the probability distribution of the largest (runaway) droplet at t=T_gel.

2.2 Some theoretical and experimental results for finite systems in coagulation theory and nuclear physics

We will now consider some results obtained for finite systems in coagulation theory (Botet, 2011) and in nuclear physics (Gruyer et al., 2013). Unlike those in infinite systems, fluctuations and correlations in a finite system are not negligible.

We must emphasize that phase transitions cannot take place in a finite system. This is due to the fact that a phase transition is defined as a singularity in the free energy or any thermodynamic property of a system. For finite-sized systems, the free energy is proportional to the logarithm of a finite number of exponentials, which are always positive (Bhattacharjee, 2001). Consequently, those singularities are only possible within infinite systems by taking the thermodynamic limit. Thus, for finite systems, the notion of pseudo-critical region is introduced (which is the finite system equivalent of a sol–gel transition time).

Some interesting simulation and experimental results were obtained for these systems in Botet (2011) for the Smoluchowski model (1) and in Gruyer et al. (2013) for nuclear multi-fragmentation. Botet et al. (2001) found, from stochastic simulations of coagulation process with the product kernel (for a system of N=512 monomers), that the distribution of the largest cluster in the pseudo-critical region can be described as a mixture of the well-known Gaussian and Gumbel distributions:

In Eq. (6), the coefficients θ and (1−θ) are the mixture weights (probabilities associated with each component). The individual distributions $\mathrm{Gumbel}(x,{\mathit{\mu }}_{\mathrm{1}},\mathit{\beta })$ and $\mathrm{Gauss}(x,{\mathit{\mu }}_{\mathrm{2}},\mathit{\sigma })$ are the mixture components.

The Gumbel distribution is one of the asymptotic distributions from extreme value theory (EVT) and has the following form:

where μ is the position parameter and β the scale parameter. The distribution in Eq. (6) has its origin in the fact that, for finite systems in the pseudo-critical zone, the system experiences large fluctuations and the gel distribution is a combination of both distributions, a Gumbel and a Gaussian (Gruyer et al., 2013). A similar result was obtained by Botet (2011) using synthetic data from stochastic simulations, for collision probabilities proportional to the product of the masses.

The fundamental hypothesis of our work is that the gel mass (largest drop) in the initial phase of precipitation formation is distributed as a mixture of two asymptotic distributions: one Gumbel and one Gaussian, following the pseudo-critical clustering scenario described in Gruyer et al. (2013).

3.1 Results for the product kernel ($K(i,j)=C{x}_i{x}_j$)

For synthetic data analysis, the empirical distributions of the largest droplet mass (M_max) were obtained from Monte Carlo simulations, following Botet (2011). The species-accounting formulation (Laurenzi et al., 2002) of the stochastic simulation algorithm (SSA) of Gillespie (1975), which rigorously accounts for fluctuations and correlations in a coalescing system, was used for the stochastic simulation in this work (see Appendix B).

The main difference between the Gillespie's SSA and other Monte Carlo methods based on the simulation particles (SIPs) approach (like the super droplet method developed by Shima et al., 2009) is that the Gillespie's SSA involved the collision of only two physical particles (droplets in our case) per MC cycle, while the approach based on SIP in each MC cycle collides SIP (super droplets, for example), which represents multiple numbers of droplets with the same attributes (radius r or mass in the simplest case) and position. However, Gillespie's SSA works perfectly for our purposes because, due to the finiteness of our systems, our simulations are performed for small volumes with a small number of droplets (in the range 50–300 cm⁻³).

Our methodology of synthetic data analysis consists of generating N-realizations (at each time step) using the algorithm of Gillespie. For each realization, there is a certain distribution of droplets. The largest droplet mass obtained from each distribution at each realization (for a fixed time step) would be the distribution to be fitted to the distribution in Eq. (6). Thus, the sample size would be equal to the number of realizations of the Monte Carlo algorithm.

Simulations were performed for the product kernel ($K(i,j)=C{x}_i{x}_j$), with an initial mono-disperse distribution of 100 droplets of 14 µm in radius (droplet mass $\mathrm{1.15}\times {\mathrm{10}}^{-\mathrm{8}}$ g) in a cloud volume of 1 cm³, with $C=\mathrm{5.49}\times {\mathrm{10}}^{\mathrm{10}}$ cm³ s⁻¹.

The product kernel is proportional to the product of the masses of the colliding droplets. It is widely used because analytical solutions of the KCE or Smoluchowski equation (Eq. 1) have been obtained for this kernel by Golovin (1963), Scott (1968), Drake (1972), and Drake and Wright (1972). The value of the constant C ($C=\mathrm{5.49}\times {\mathrm{10}}^{\mathrm{10}}$ cm³ g⁻² s⁻¹) in the product kernel is the result of the polynomial approximation $K(x,y)=A+B(x+y)+Cxy$ (Long, 1974) of the hydrodynamic collection kernel (Eq. 11).

The empirical distribution of the maxima was obtained for 1000 realizations of the stochastic algorithm. There is no need for a larger number of realizations to get better statistics, since the number of realizations in our Monte Carlo algorithm must be equal to the sample size in the application of the block maxima (BM) approach (see the next section for more details). On the other hand, this number is not much bigger than the number of blocks in the data for which the largest droplet maxima was fitted to fog data.

Figure 1a–d present the largest droplet mass empirical distributions obtained at four different times. Note that Eq. (6) provides a good fit for the distribution of the mass of the largest droplet (M_max) both around and far from the sol–gel transition time (T_gel), which was calculated from Eq. (4) and found to be equal to 1378 s.

Figure 1Panels (a)–(d) (histograms) show the largest droplet mass distributions calculated from Monte Carlo simulations at four different times, for a system with an initial mono-disperse distribution of 100 droplets of 14 µm in radius. The solid line shows the fit using Eq. (6).

Download

Figure 2 presents the time evolution of the coefficient θ, which represents the mixing fraction in Eq. (6) for the time interval [500 s, 2000 s]. Despite the noisy behavior of the coefficient θ (due to the finiteness of the system), there is a decreasing trend with time, showing larger values of θ (∼0.65) for times close to 500 s and values down to 0.2 at the end of the time interval. This figure indicates that, although the largest droplet distribution is adequately described by a mixture of Gaussian and Gumbel distributions, it has a larger Gumbel component (see Eq. 6) during the early stages of the coagulation process. As time progresses, the Gaussian contribution becomes more important (smaller values of θ) in providing a better fit to the largest droplet mass distribution.

Figure 2Time evolution of the coefficient θ in Eq. (6), obtained for a simulation with the product kernel.

Download

These findings are in accordance with Gruyer et al. (2013) and Botet (2011): at an early stage of coagulation development, correlations are negligible, and, consequently, the largest fragments can be considered independent random variables. Therefore, the probability distribution of the largest fragment is given by the limit theorem for extremal variables, which states that the maximum of sample-independent and identically distributed random variables can only converge in distribution in the form of one of three possible distributions: Gumbel, Fréchet or Weibull.

As the coagulation process continues, fluctuations and correlations between droplets increase and the system reaches a critical point (Alfonso and Raga, 2017). Where the largest droplets are no longer independent random variables, the limit theorem for extremal variables no longer applies, and the largest droplet distribution is no longer described by a Gumbel distribution. At later times, away from the pseudo-critical region, the Gaussian contribution is the most important part of the largest droplet mass distribution. This can be explained by the additive nature of the process at this stage (Botet, 2011; Gruyer et al., 2013; Clusel and Bertoin, 2008), and the central limit theorem applies.

In the intermediate zone (which can be defined as the pseudo-critical zone), the distribution is well described by a mixture of Gumbel and Gaussian distributions and the weights associated with each distribution are comparable. It is expected that it can be observed that θ=0.5 at the infinite system critical point, T_gel, found to be 1378 s from Eq. (4). However, due to the finiteness of the system, the critical point corresponds approximately to a value θ=0.35 (see Fig. 2).

We can find whether or not a system is in the pseudo-critical region by defining the following ratio (Botet, 2011; Gruyer et al., 2013):

where w_Gumbel=θ and $w_{\mathrm{Gaussian}}=\mathrm{1}-\mathit{\theta }$ are the relative weights of the Gumbel and Gaussian distributions, respectively (see Eq. 6). By definition, $\mathit{\eta }=+\mathrm{1},-\mathrm{1}$ corresponds to pure Gaussian and Gumbel distributions. If $-\mathrm{1}<\mathit{\eta }<\mathrm{1}$, the system is in the pseudo-critical domain.

Alternatively, Botet (2011) estimates the limits of the pseudo-critical region as the times when the largest droplet mass standard deviation σ(M_max) calculated from Eq. (9) is small.

In Eq. (9), N_r is the number of iterations of the stochastic simulation algorithm of Gillespie (1975), M_max the mass of the largest particle and 〈M_max〉 its ensemble mean over all the realizations.

Even though the second moment of the distribution M₂(τ) diverges (see Eq. 3) for the infinite system, there is no divergence of the second moment for a finite system (with no critical behavior). For that case, the standard deviation for the largest particle mass (σ(M_max)) is expected to reach a maximum in the vicinity of $T_{\mathrm{gel}}={\left[C{M}_{\mathrm{2}}\left(t_{\mathrm{0}}\right)\right]}^{-\mathrm{1}}$. Moreover, computing the time evolution of the normalized standard deviation $\mathit{\sigma }\left(M_{max}\right)/〈M_{max}〉$ instead of σ(M_max) yielded better results in estimating T_gel in Inaba et al. (1999), Alfonso et al. (2008, 2010, 2013), and Alfonso and Raga (2017).

Figure 3a shows the time evolution of $\mathit{\sigma }\left(M_{max}\right)/〈M_{max}〉$, as an example for the system defined at the beginning of this section. Note that the maximum occurs at T=1315 s, close to T_gel=1378 s calculated from Eq. (4), and the time when the maximum of $\mathit{\sigma }\left(M_{max}\right)/〈M_{max}〉$ occurs is a reliable estimate of the sol–gel transition time for the corresponding infinite system.

Figure 3For the finite system, the normalized standard deviation $\mathit{\sigma }\left(M_{max}\right)/〈M_{max}〉$ of the largest droplet mass versus time (a). The initial number of droplets was set equal to N=100 droplets of 14 µm in radius in a volume of 1 cm³. Simulations were performed with the product kernel $K(i,j)=C{x}_i{x}_j$ (with $C=\mathrm{5.49}\times {\mathrm{10}}^{\mathrm{10}}$ cm³ g⁻² s⁻¹), and N_r=1000 realizations of the stochastic algorithm were performed. The maximum value of $\mathit{\sigma }\left(M_{max}\right)/〈M_{max}〉$ is found to be 1315 s (dashed vertical line) and is very close to the sol–gel transition time (continuous vertical line) for the infinite system (1378 s). In panel (b) the small end of the pseudo-critical domain is estimated as the time where $\mathit{\sigma }\left(M_{max}\right)=\mathrm{0.1}{\mathit{\sigma }}_{max}$.

Download

Botet (2011) defines $\mathit{\sigma }=\mathrm{0.1}{\mathit{\sigma }}_{max}$ as the limits of the pseudo-critical interval, which corresponds to $t_{inf}=\mathrm{0.37}{T}_{\mathrm{gel}}$ and $t_{sup}=\mathrm{1.5}{T}_{\mathrm{gel}}$ (see Fig. 3b). While Eq. (8) could be used to determine if a sample collected inside a cloud is in the pseudo-critical region, Eq. (9) implies that the time evolution of σ(M_max) is needed, and therefore a practical application is only viable in the case of synthetic data obtained from stochastic simulations or cloud droplet data collected dynamically at different times or cloud levels.

3.2 Numerical results for turbulent conditions

In our simulations, turbulent effects were considered by implementing the turbulence-induced collision enhancement factor P_Turb(x_i,x_j) that is calculated in Pinsky et al. (2008) for a cumulonimbus with dissipation rate ε=0.1 m² s⁻³ and Reynolds number ${\mathit{\text{Re}}}_{\mathit{\lambda }}=\mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}$ and for cloud droplets with radii ≤21 µm. The turbulent collection kernel has the following form:

where K_g(x_i,x_j) is hydrodynamic kernel, which considers collisions between droplets under pure gravity conditions and has the following form:

The hydrodynamic kernel takes into account the fact that droplets with different masses (x_i and x_j and corresponding radii, r_i and r_j) have different terminal velocities V(x_i), which are functions of their masses. In Eq. (10), E(r_i) and E(r_j) are the collection efficiencies calculated according to Hall (1980).

Monte Carlo simulations were performed with an initial bi-modal distribution (200 droplets of 10 µm in radius and 50 droplets of 12.6 µm) for a cloud volume of 1 cm³.

As we want to perform simulations for small systems (with a small number of particles) for which fluctuations and correlations are relevant, the number of droplets per cubic centimeter used in the simulations are small. They are of the same order of the droplet concentrations for each block obtained from observations, which fluctuate between 0 and 392 cm⁻³, with an average of 146 cm⁻³ (see Fig. 6).

The empirical distribution for the largest droplet mass was generated by extracting the maximum from the droplet distribution at each realization for a fixed time step. Additionally, the ratio $\mathit{\sigma }\left(M_{max}\right)/〈M_{max}〉$ is evaluated from 1000 realizations of the Monte Carlo algorithm (see Fig. 4), which reaches its maximum at around 1815 s and serves as an estimate for the sol–gel transition time for the infinite system. Four empirical probability distributions were fitted to the combined distribution (Eq. 6) for times in the vicinity of T_gel. The results are displayed in Fig. 5a–d. Note also that for this case, the combined distribution (Eq. 6) provides a good fit for the largest droplet mass. Moreover, the coefficient θ decreases in time (check Fig. 5), in concordance with the scenario described in Sect. 3.1.

Figure 4Time evolution of the normalized standard deviation $\mathit{\sigma }\left(M_{max}\right)/〈M_{max}〉$ of the largest droplet mass versus time estimated from the Monte Carlo algorithm. The simulations were performed for the turbulent hydrodynamic kernel with a bi-disperse initial condition (200 droplets of 10 µm in radius and 50 droplets of 12.6 µm) in a volume of 1 cm³.

Download

Figure 5Panels (a)–(d) (histograms) show the simulated M_max distributions in a system with an initial bi-disperse distribution (200 droplets of 10 µm in radius and 50 droplets of 12.6 µm) at four different times. The solid line shows the fit using Eq. (6). The simulations were performed for the turbulent hydrodynamic kernel.

Download